AliPhysics  421aab4 (421aab4)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliJetModelBaseTask.cxx
Go to the documentation of this file.
1 //
2 // Jet modeling task.
3 //
4 // Author: S.Aiola, C.Loizides
5 
6 #include "AliJetModelBaseTask.h"
7 
8 #include <TClonesArray.h>
9 #include <TH1I.h>
10 #include <TLorentzVector.h>
11 #include <TRandom3.h>
12 #include <TList.h>
13 #include <TF2.h>
14 #include <TGrid.h>
15 #include <TFile.h>
16 
17 #include "AliVEvent.h"
18 #include "AliAODCaloCluster.h"
19 #include "AliESDCaloCluster.h"
20 #include "AliVCluster.h"
21 #include "AliEMCALDigit.h"
22 #include "AliEMCALRecPoint.h"
23 #include "AliESDCaloCells.h"
24 #include "AliAODCaloCells.h"
25 #include "AliAODMCParticle.h"
26 #include "AliVCaloCells.h"
27 #include "AliPicoTrack.h"
28 #include "AliEMCALGeometry.h"
29 #include "AliLog.h"
30 #include "AliNamedArrayI.h"
31 
33 
34 //________________________________________________________________________
37  fGeomName(),
38  fTracksName(),
39  fOutTracksName(),
40  fCaloName(),
41  fOutCaloName(),
42  fCellsName(),
43  fOutCellsName(),
44  fMCParticlesName(),
45  fOutMCParticlesName(),
46  fPythiaInfoName(""),
47  fIsMC(kFALSE),
48  fSuffix(),
49  fEtaMin(-1),
50  fEtaMax(1),
51  fPhiMin(0),
52  fPhiMax(TMath::Pi() * 2),
53  fPtMin(0),
54  fPtMax(0),
55  fGenType(0),
56  fCopyArray(kTRUE),
57  fNClusters(0),
58  fNCells(0),
59  fNTracks(0),
60  fMarkMC(99999),
61  fPtSpectrum(0),
62  fPtPhiEvPlDistribution(0),
63  fDensitySpectrum(0),
64  fDifferentialV2(0),
65  fAddV2(kFALSE),
66  fFlowFluctuations(kFALSE),
67  fQAhistos(kFALSE),
68  fPsi(0),
69  fIsInit(0),
70  fGeom(0),
71  fClusters(0),
72  fOutClusters(0),
73  fTracks(0),
74  fOutTracks(0),
75  fCaloCells(0),
76  fOutCaloCells(0),
77  fAddedCells(0),
78  fMCParticles(0),
79  fMCParticlesMap(0),
80  fOutMCParticles(0),
81  fOutMCParticlesMap(0),
82  fMCLabelShift(0),
83  fEsdMode(kFALSE),
84  fOutput(0),
85  fPythiaInfo(0x0),
86  fhpTEmb(0),
87  fhMEmb(0),
88  fhEtaEmb(0),
89  fhPhiEmb(0),
90  fMassFromDistr(kFALSE),
91  fHMassDistrib(0),
92  fHMassPtDistrib(0)
93 {
94  // Default constructor.
95 
96  fVertex[0] = 0;
97  fVertex[1] = 0;
98  fVertex[2] = 0;
99 }
100 
101 //________________________________________________________________________
103  AliAnalysisTaskSE(name),
104  fGeomName("EMCAL_COMPLETE12SMV1"),
105  fTracksName("PicoTracks"),
106  fOutTracksName("PicoTracksEmbedded"),
107  fCaloName("CaloClustersCorr"),
108  fOutCaloName("CaloClustersCorrEmbedded"),
109  fCellsName(""),
110  fOutCellsName(""),
111  fMCParticlesName(""),
112  fOutMCParticlesName(""),
113  fPythiaInfoName(""),
114  fIsMC(kFALSE),
115  fSuffix("Processed"),
116  fEtaMin(-1),
117  fEtaMax(1),
118  fPhiMin(0),
119  fPhiMax(TMath::Pi() * 2),
120  fPtMin(50),
121  fPtMax(60),
122  fGenType(0),
123  fCopyArray(kTRUE),
124  fNClusters(0),
125  fNCells(0),
126  fNTracks(1),
127  fMarkMC(99999),
128  fPtSpectrum(0),
129  fPtPhiEvPlDistribution(0),
130  fDensitySpectrum(0),
131  fDifferentialV2(0),
132  fAddV2(kFALSE),
133  fFlowFluctuations(kFALSE),
134  fQAhistos(drawqa),
135  fPsi(0),
136  fIsInit(0),
137  fGeom(0),
138  fClusters(0),
139  fOutClusters(0),
140  fTracks(0),
141  fOutTracks(0),
142  fCaloCells(0),
143  fOutCaloCells(0),
144  fAddedCells(0),
145  fMCParticles(0),
146  fMCParticlesMap(0),
147  fOutMCParticles(0),
148  fOutMCParticlesMap(0),
149  fMCLabelShift(0),
150  fEsdMode(kFALSE),
151  fOutput(0),
152  fPythiaInfo(0x0),
153  fhpTEmb(0),
154  fhMEmb(0),
155  fhEtaEmb(0),
156  fhPhiEmb(0),
157  fMassFromDistr(kFALSE),
158  fHMassDistrib(0),
159  fHMassPtDistrib(0)
160 {
161  // Standard constructor.
162 
163  if (fQAhistos) {
164  DefineOutput(1, TList::Class());
165  }
166 
167  fVertex[0] = 0;
168  fVertex[1] = 0;
169  fVertex[2] = 0;
170 }
171 
172 //________________________________________________________________________
174 {
175  // Destructor
176 }
177 
178 //________________________________________________________________________
180  Printf("Note: If the code changes this method could give a wrong result");
181  TString futurenamefOutputTrack;
182  if(fCopyArray) futurenamefOutputTrack = Form("%s%s", fTracksName.Data(), fSuffix.Data());
183  else futurenamefOutputTrack = fTracksName;
184  return futurenamefOutputTrack;
185 
186 }
187 //________________________________________________________________________
189 {
190  // Create user output.
191  if (!fQAhistos)
192  return;
193 
194  OpenFile(1);
195  fOutput = new TList();
196  fOutput->SetOwner();
197 
198  fhpTEmb = new TH1F("fhpTEmb","#it{p}_{T} distribution; #it{p}_{T}(GeV/c)", 120, 0., 120);
199  fhpTEmb->Sumw2();
200  fOutput->Add(fhpTEmb);
201 
202  fhMEmb = new TH1F("fhMEmb","Mass distribution; #it{M} (GeV)", 80, 0, 80.);
203  fhMEmb->Sumw2();
204  fOutput->Add(fhMEmb);
205 
206  fhEtaEmb = new TH1F("fhEtaEmb","#eta distribution; #eta", 100, -1, 1);
207  fhEtaEmb->Sumw2();
208  fOutput->Add(fhEtaEmb);
209 
210  fhPhiEmb = new TH1F("fhPhiEmb","#varphi distribution; #varphi", 100, (-1)*TMath::Pi(), TMath::Pi());
211  fhPhiEmb->Sumw2();
212  fOutput->Add(fhPhiEmb);
213 
214  fhEvents = new TH1I("fhEvents", "Number of events", 3, 0, 2);
215  fOutput->Add(fhEvents);
216 
217  PostData(1, fOutput);
218 }
219 
220 //________________________________________________________________________
222 {
223  // Execute per event.
224 
225  if (!fIsInit)
226  fIsInit = ExecOnce();
227 
228  if (!fIsInit)
229  return;
230 
231  fVertex[0] = 0;
232  fVertex[1] = 0;
233  fVertex[2] = 0;
234 
235  const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
236  if (vert)
237  vert->GetXYZ(fVertex);
238 
239  if (fCopyArray) {
240  if (fOutTracks)
241  fOutTracks->Delete();
242  if (fOutClusters)
243  fOutClusters->Delete();
244  if (fOutMCParticles)
245  fOutMCParticles->Delete();
246  }
247 
248  if (fDensitySpectrum) {
249  fNTracks = TMath::Nint(fDensitySpectrum->GetRandom());
250  fNCells = TMath::Nint(fDensitySpectrum->GetRandom());
251  fNClusters = TMath::Nint(fDensitySpectrum->GetRandom());
252  }
253 
254  // Clear map
255  if (fOutMCParticlesMap)
257 
258  AliVCaloCells *tempCaloCells = 0;
259 
260  if (fCaloCells) {
261  fAddedCells = 0;
262  if (!fCopyArray) {
263  tempCaloCells = fCaloCells;
264  fCaloCells = static_cast<AliVCaloCells*>(tempCaloCells->Clone(Form("%s_old",fCaloCells->GetName())));
265  }
266  }
267 
269  fPsi = gRandom->Rndm() * TMath::Pi();
270 
271  Run();
272 
273  if (fCaloCells && !fCopyArray) {
274  delete fCaloCells;
275  fCaloCells = tempCaloCells;
276  }
277 }
278 
279 //________________________________________________________________________
281 {
282  // Init task.
283 
284  delete gRandom;
285  gRandom = new TRandom3(0);
286 
287  fEsdMode = InputEvent()->InheritsFrom("AliESDEvent");
288 
289  if (fPtMax < fPtMin) {
290  AliWarning (Form("PtMax (%f) < PtMin (%f), setting PtMax = PtMin = %f", fPtMax, fPtMin, fPtMin));
291  fPtMax = fPtMin;
292  }
293 
294  if (fEtaMax < fEtaMin) {
295  AliWarning (Form("EtaMax (%f) < EtaMin (%f), setting EtaMax = EtaMin = %f", fEtaMax, fEtaMin, fEtaMin));
296  fEtaMax = fEtaMin;
297  }
298 
299  if (fPhiMax < fPhiMin) {
300  AliWarning (Form("PhiMax (%f) < PhiMin (%f), setting PhiMax = PhiMin = %f", fPhiMax, fPhiMin, fPhiMin));
301  fPhiMax = fPhiMin;
302  }
303 
304  if (!fCellsName.IsNull()) {
305  fCaloCells = dynamic_cast<AliVCaloCells*>(InputEvent()->FindListObject(fCellsName));
306  if (!fCaloCells) {
307  AliWarning(Form("%s: Couldn't retrieve calo cells with name %s!", GetName(), fCellsName.Data()));
308  }
309  else if (!fCaloCells->InheritsFrom("AliVCaloCells")) {
310  AliError(Form("%s: Collection %s does not contain a AliVCaloCells object!", GetName(), fCellsName.Data()));
311  fCaloCells = 0;
312  return kFALSE;
313  }
314 
315  if (!fOutCaloCells) {
317  if (fCopyArray)
319  if (fCopyArray || !fCaloCells) {
320  if (fEsdMode)
321  fOutCaloCells = new AliESDCaloCells(fOutCellsName,fOutCellsName);
322  else
323  fOutCaloCells = new AliAODCaloCells(fOutCellsName,fOutCellsName);
324 
325  if (InputEvent()->FindListObject(fOutCellsName)) {
326  AliFatal(Form("%s: Collection %s is already present in the event!", GetName(), fOutCellsName.Data()));
327  return kFALSE;
328  }
329  else {
330  InputEvent()->AddObject(fOutCaloCells);
331  }
332  }
333  else {
335  }
336  }
337  }
338 
339  if (!fTracksName.IsNull()) {
340  fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName));
341  if (!fTracks) {
342  AliWarning(Form("%s: Couldn't retrieve tracks with name %s!", GetName(), fTracksName.Data()));
343  }
344  else if (!fTracks->GetClass()->GetBaseClass("AliPicoTrack")) {
345  AliError(Form("%s: Collection %s does not contain AliPicoTrack objects!", GetName(), fTracksName.Data()));
346  fTracks = 0;
347  return kFALSE;
348  }
349  if (!fOutTracks) {
351  if (fCopyArray)
353  if (fCopyArray || !fTracks) {
354  fOutTracks = new TClonesArray("AliPicoTrack");
355  //fOutTracks->SetOwner(kTRUE);
356  fOutTracks->SetName(fOutTracksName);
357  if (InputEvent()->FindListObject(fOutTracksName)) {
358  AliFatal(Form("%s: Collection %s is already present in the event!", GetName(), fOutTracksName.Data()));
359  return kFALSE;
360  }
361  else {
362  InputEvent()->AddObject(fOutTracks);
363  }
364  }
365  else {
367  InputEvent()->AddObject(fOutTracks);
368  }
369  }
370  InputEvent()->ls();
371  }
372 
373  if(fAddV2 && (!fDifferentialV2)) {
374  AliWarning(Form("%s: Cannot add v2 without diffential v2!", GetName()));
375  }
376 
377  if (!fCaloName.IsNull()) {
378  fClusters = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCaloName));
379 
380  if (!fClusters) {
381  AliWarning(Form("%s: Couldn't retrieve clusters with name %s!", GetName(), fCaloName.Data()));
382  }
383  else if (!fClusters->GetClass()->GetBaseClass("AliVCluster")) {
384  AliError(Form("%s: Collection %s does not contain AliVCluster objects!", GetName(), fCaloName.Data()));
385  fClusters = 0;
386  return kFALSE;
387  }
388 
389  if (!fOutClusters) {
391  if (fCopyArray)
393  TString className;
394  if (fClusters)
395  className = fClusters->GetClass()->GetName();
396  else if (fEsdMode)
397  className = "AliESDCaloCluster";
398  else
399  className = "AliAODCaloCluster";
400 
401  if (fCopyArray || !fClusters) {
402  fOutClusters = new TClonesArray(className.Data());
403  fOutClusters->SetName(fOutCaloName);
404  if (InputEvent()->FindListObject(fOutCaloName)) {
405  AliFatal(Form("%s: Collection %s is already present in the event!", GetName(), fOutCaloName.Data()));
406  return kFALSE;
407  }
408  else {
409  InputEvent()->AddObject(fOutClusters);
410  }
411  }
412  else {
414  }
415  }
416  }
417 
418  if (!fMCParticlesName.IsNull()) {
419  fMCParticles = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fMCParticlesName));
420  if (!fMCParticles) {
421  AliWarning(Form("%s: Couldn't retrieve MC particles with name %s!", GetName(), fMCParticlesName.Data()));
422  }
423  else {
424  if (!fMCParticles->GetClass()->GetBaseClass("AliAODMCParticle")) {
425  AliError(Form("%s: Collection %s does not contain AliAODMCParticle objects!", GetName(), fMCParticlesName.Data()));
426  fMCParticles = 0;
427  return kFALSE;
428  }
429 
430  fMCParticlesMap = dynamic_cast<AliNamedArrayI*>(InputEvent()->FindListObject(fMCParticlesName + "_Map"));
431 
432  if (!fMCParticlesMap) {
433  AliWarning(Form("%s: Could not retrieve map for MC particles %s! Will assume MC labels consistent with indexes...", GetName(), fMCParticlesName.Data()));
434  fMCParticlesMap = new AliNamedArrayI(fMCParticlesName + "_Map", 99999);
435  for (Int_t i = 0; i < 99999; i++) {
436  fMCParticlesMap->AddAt(i,i);
437  }
438  }
439  }
440 
441  if (!fOutMCParticles) {
443  if (fCopyArray)
445  if (fCopyArray || !fMCParticles) {
446 
447  fOutMCParticles = new TClonesArray("AliAODMCParticle");
449  if (InputEvent()->FindListObject(fOutMCParticlesName)) {
450  AliFatal(Form("%s: Collection %s is already present in the event!", GetName(), fOutMCParticlesName.Data()));
451  return kFALSE;
452  }
453  else {
454  InputEvent()->AddObject(fOutMCParticles);
455  }
456 
458  if (InputEvent()->FindListObject(fOutMCParticlesName + "_Map")) {
459  AliFatal(Form("%s: Map %s_Map is already present in the event!", GetName(), fOutMCParticlesName.Data()));
460  return kFALSE;
461  }
462  else {
463  InputEvent()->AddObject(fOutMCParticlesMap);
464  }
465  }
466  else {
469  }
470  }
471  }
472 
473  if (!fGeom) {
474  if (fGeomName.Length() > 0) {
475  fGeom = AliEMCALGeometry::GetInstance(fGeomName);
476  if (!fGeom) {
477  AliFatal(Form("Could not get geometry with name %s!", fGeomName.Data()));
478  return kFALSE;
479  }
480  } else {
481  fGeom = AliEMCALGeometry::GetInstance();
482  if (!fGeom) {
483  AliFatal("Could not get default geometry!");
484  return kFALSE;
485  }
486  }
487  }
488 
489  return kTRUE;
490 }
491 
492 //________________________________________________________________________
494 {
495  if (fOutCaloCells->GetNumberOfCells() < n) {
496  fOutCaloCells->DeleteContainer();
497  fOutCaloCells->CreateContainer(n);
498  }
499  else {
500  fOutCaloCells->SetNumberOfCells(n);
501  }
502 
503  fAddedCells = 0;
504 
505  return n;
506 }
507 
508 //________________________________________________________________________
510 {
511  // Add a cell to the event.
512 
513  Int_t absId = 0;
514  if (eta < -100 || phi < 0) {
515  GetRandomCell(eta, phi, absId);
516  }
517  else {
518  fGeom->EtaPhiFromIndex(absId, eta, phi);
519  }
520 
521  if (absId == -1) {
522  AliWarning(Form("Unable to embed cell in eta = %f, phi = %f!"
523  " Maybe the eta-phi range is not inside the EMCal acceptance (eta = [%f, %f], phi = [%f, %f])",
524  eta, phi, fEtaMin, fEtaMax, fPhiMin, fPhiMax));
525  return 0;
526  }
527 
528  if (e < 0) {
529  Double_t pt = GetRandomPt();
530  TLorentzVector nPart;
531  nPart.SetPtEtaPhiM(pt, eta, phi, 0);
532  e = nPart.E();
533  }
534 
535  return AddCell(e, absId);
536 }
537 
538 //________________________________________________________________________
540 {
541  // Add a cell to the event.
542 
543  if (label < 0)
544  label = 0;
545 
546  label += fMarkMC + fMCLabelShift;
547 
548  Short_t pos = -1;
549  if (fCaloCells)
550  pos = fCaloCells->GetCellPosition(absId);
551 
552  Double_t efrac = 1;
553  Bool_t increaseOnSuccess = kFALSE;
554 
555  if (pos < 0) {
556  increaseOnSuccess = kTRUE;
557  pos = fAddedCells;
558  }
559  else {
560  Short_t cellNumber = -1;
561  Double_t old_e = 0;
562  Double_t old_time = 0;
563  Int_t old_label = 0;
564  Double_t old_efrac = 0;
565  fOutCaloCells->GetCell(pos, cellNumber, old_e, old_time, old_label, old_efrac);
566 
567  efrac = e / (old_e + e);
568 
569  if (old_label > 0 && e < old_efrac * old_e) {
570  label = old_label;
571  efrac = old_efrac;
572  time = old_time;
573  }
574 
575  e += old_e;
576  }
577 
578  Bool_t r = fOutCaloCells->SetCell(pos, absId, e, time, label, efrac);
579 
580  if (r) {
581  if (increaseOnSuccess)
582  fAddedCells++;
583  return fAddedCells;
584  }
585  else
586  return 0;
587 }
588 
589 //________________________________________________________________________
590 AliVCluster* AliJetModelBaseTask::AddCluster(AliVCluster *oc)
591 {
592  // Add a cluster to the event.
593 
594  const Int_t nClusters = fOutClusters->GetEntriesFast();
595  AliVCluster *dc = static_cast<AliVCluster*>(fOutClusters->New(nClusters));
596  dc->SetType(AliVCluster::kEMCALClusterv1);
597  dc->SetE(oc->E());
598  Float_t pos[3] = {0};
599  oc->GetPosition(pos);
600  dc->SetPosition(pos);
601  dc->SetNCells(oc->GetNCells());
602  dc->SetCellsAbsId(oc->GetCellsAbsId());
603  dc->SetCellsAmplitudeFraction(oc->GetCellsAmplitudeFraction());
604  dc->SetID(oc->GetID());
605  dc->SetDispersion(oc->GetDispersion());
606  dc->SetEmcCpvDistance(-1);
607  dc->SetChi2(-1);
608  dc->SetTOF(oc->GetTOF()); //time-of-flight
609  dc->SetNExMax(oc->GetNExMax()); //number of local maxima
610  dc->SetM02(oc->GetM02());
611  dc->SetM20(oc->GetM20());
612  dc->SetDistanceToBadChannel(oc->GetDistanceToBadChannel());
613 
614  //MC
615  UInt_t nlabels = oc->GetNLabels();
616  Int_t *labels = oc->GetLabels();
617  TArrayI parents;
618  if (nlabels != 0 && labels)
619  parents.Set(nlabels, labels);
620  else {
621  nlabels = 1;
622  parents.Set(1);
623  parents[0] = 0;
624  }
625 
626  if (fMarkMC+fMCLabelShift != 0) {
627  for (UInt_t i = 0; i < nlabels; i++) {
628  parents[i] += fMarkMC+fMCLabelShift;
629  }
630  }
631 
632  AliESDCaloCluster *esdClus = dynamic_cast<AliESDCaloCluster*>(dc);
633  if (esdClus) {
634  esdClus->AddLabels(parents);
635  }
636  else {
637  AliAODCaloCluster *aodClus = dynamic_cast<AliAODCaloCluster*>(dc);
638  if (aodClus)
639  aodClus->SetLabel(parents.GetArray(), nlabels);
640  }
641 
642  return dc;
643 }
644 
645 //________________________________________________________________________
647 {
648  // Add a cluster to the event.
649 
650  Int_t absId = 0;
651  if (eta < -100 || phi < 0) {
652  GetRandomCell(eta, phi, absId);
653  }
654  else {
655  fGeom->EtaPhiFromIndex(absId, eta, phi);
656  }
657 
658  if (absId == -1) {
659  AliWarning(Form("Unable to embed cluster in eta = %f, phi = %f!"
660  " Maybe the eta-phi range is not inside the EMCal acceptance (eta = [%f, %f], phi = [%f, %f])",
661  eta, phi, fEtaMin, fEtaMax, fPhiMin, fPhiMax));
662  return 0;
663  }
664 
665  if (e < 0) {
666  Double_t pt = GetRandomPt();
667  TLorentzVector nPart;
668  nPart.SetPtEtaPhiM(pt, eta, phi, 0);
669  e = nPart.E();
670  }
671 
672  return AddCluster(e, absId, label);
673 }
674 
675 //________________________________________________________________________
677 {
678  // Add a cluster to the event.
679 
680  const Int_t nClusters = fOutClusters->GetEntriesFast();
681 
682  TClonesArray digits("AliEMCALDigit", 1);
683 
684  AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(digits.New(0));
685  digit->SetId(absId);
686  digit->SetIndexInList(0);
687  digit->SetType(AliEMCALDigit::kHG);
688  digit->SetAmplitude(e);
689 
690  AliEMCALRecPoint recPoint("");
691  recPoint.AddDigit(*digit, e, kFALSE);
692  recPoint.EvalGlobalPosition(0, &digits);
693 
694  TVector3 gpos;
695  recPoint.GetGlobalPosition(gpos);
696  Float_t g[3];
697  gpos.GetXYZ(g);
698 
699  AliVCluster *cluster = static_cast<AliVCluster*>(fOutClusters->New(nClusters));
700  cluster->SetType(AliVCluster::kEMCALClusterv1);
701  cluster->SetE(recPoint.GetEnergy());
702  cluster->SetPosition(g);
703  cluster->SetNCells(1);
704  UShort_t shortAbsId = absId;
705  cluster->SetCellsAbsId(&shortAbsId);
706  Double32_t fract = 1;
707  cluster->SetCellsAmplitudeFraction(&fract);
708  cluster->SetID(nClusters);
709  cluster->SetEmcCpvDistance(-1);
710 
711  //MC label
712  if (label < 0)
713  label = 0;
714 
715  label += fMarkMC+fMCLabelShift;
716 
717  if (fEsdMode) {
718  AliESDCaloCluster *esdClus = static_cast<AliESDCaloCluster*>(cluster);
719  TArrayI parents(1, &label);
720  esdClus->AddLabels(parents);
721  }
722  else {
723  AliAODCaloCluster *aodClus = static_cast<AliAODCaloCluster*>(cluster);
724  aodClus->SetLabel(&label, 1);
725  }
726 
727  return cluster;
728 }
729 
730 //________________________________________________________________________
732 {
733  // Add a track to the event.
734  if (pt < 0 && eta < -100 && phi < -100 && mass < -100) {
735  GetRandomMassiveParticle(pt,eta,phi, kFALSE, mass);
736  } else {
737  if (pt < 0 && eta < -100 && phi < -100) {
738  GetRandomParticle(pt,eta,phi);
739  }
740  else {
741  if (pt < -100)
742  pt = GetRandomPt();
743  if (eta < -100)
744  eta = GetRandomEta();
745  if (phi < -100)
746  phi = GetRandomPhi(pt);
747  }
748  }
749  //Printf("Adding LABEL %d", label);
750  if (label >= 0)
751  label += fMarkMC+fMCLabelShift;
752  else if (label < 0)
753  label -= fMarkMC+fMCLabelShift;
754  if(fAddV2) AddV2(phi, pt);
755 
756  const Int_t nTracks = fOutTracks->GetEntriesFast();
757  //Printf("+ %d = %d", fMarkMC, label);
758  AliPicoTrack *track = new ((*fOutTracks)[nTracks]) AliPicoTrack(pt,
759  eta,
760  phi,
761  charge,
762  label,
763  type,
764  etaemc,
765  phiemc,
766  ptemc,
767  ise,
768  mass);
769 
770  return track;
771 }
772 
773 //________________________________________________________________________
774 AliAODMCParticle* AliJetModelBaseTask::AddMCParticle(AliAODMCParticle *part, Int_t origIndex)
775 {
776  const Int_t nPart = fOutMCParticles->GetEntriesFast();
777 
778  AliAODMCParticle *aodpart = new ((*fOutMCParticles)[nPart]) AliAODMCParticle(*part);
779 
780  if (origIndex + fMCLabelShift >= fOutMCParticlesMap->GetSize())
781  fOutMCParticlesMap->Set((origIndex + fMCLabelShift)*2);
782 
783  fOutMCParticlesMap->AddAt(nPart, origIndex + fMCLabelShift);
784  AliDebug(2, Form("Setting bin %d to %d (fMCLabelShift=%d, origIndex=%d)",
785  origIndex + fMCLabelShift, fOutMCParticlesMap->At(origIndex + fMCLabelShift), fMCLabelShift, origIndex));
786 
787  return aodpart;
788 }
789 
790 //_____________________________________________________________________________
792 {
793  // similar to AliFlowTrackSimple::AddV2, except for the flow fluctuations
794  Double_t phi0(phi), v2(0.), f(0.), fp(0.), phiprev(0.);
795  if(fDifferentialV2) v2 = fDifferentialV2->Eval(pt);
796  if(TMath::AreEqualAbs(v2, 0, 1e-5)) return;
797  // introduce flow fluctuations (gaussian)
798  if(fFlowFluctuations) v2 += TMath::Sqrt(2*(v2*.25)*(v2*.25))*TMath::ErfInverse(2*(gRandom->Uniform(0, fFlowFluctuations))-1);
799  for (Int_t i(0); i < 100; i++) {
800  phiprev=phi; //store last value for comparison
801  f = phi-phi0+v2*TMath::Sin(2.*(phi-fPsi));
802  fp = 1.0+2.0*v2*TMath::Cos(2.*(phi-fPsi)); //first derivative
803  phi -= f/fp;
804  if (TMath::AreEqualAbs(phiprev, phi, 1e-10)) break;
805  }
806 }
807 
808 //_____________________________________________________________________________
810 {
811  if (!fCaloCells)
812  return;
813 
814  fAddedCells = 0;
815  fCaloCells->Sort();
816  for (Short_t i = 0; i < fCaloCells->GetNumberOfCells(); i++) {
817  Int_t mclabel = 0;
818  Double_t efrac = 0.;
819  Double_t time = -1;
820  Short_t cellNum = -1;
821  Double_t amp = -1;
822 
823  fCaloCells->GetCell(i, cellNum, amp, time, mclabel, efrac);
824 
825  if (!fIsMC)
826  mclabel = 0;
827 
828  fOutCaloCells->SetCell(i, cellNum, amp, time, mclabel, efrac);
829  fAddedCells++;
830  }
831 
832  AliDebug(2, Form("%d cells from the current event", fAddedCells));
833 }
834 
835 //________________________________________________________________________
837 {
838  // Copy all the clusters in the new collection
839  if (!fClusters)
840  return;
841 
842  const Int_t nClusters = fClusters->GetEntriesFast();
843  Int_t nCopiedClusters = 0;
844 
845  if (fEsdMode) {
846  for (Int_t i = 0; i < nClusters; ++i) {
847  AliESDCaloCluster *esdcluster = static_cast<AliESDCaloCluster*>(fClusters->At(i));
848  if (!esdcluster || !esdcluster->IsEMCAL())
849  continue;
850  AliESDCaloCluster *clus = new ((*fOutClusters)[nCopiedClusters]) AliESDCaloCluster(*esdcluster);
851  if (!fIsMC) {
852  TArrayI *labels = clus->GetLabelsArray();
853  if (labels)
854  labels->Reset();
855  }
856  nCopiedClusters++;
857  }
858  }
859  else {
860  for (Int_t i = 0; i < nClusters; ++i) {
861  AliAODCaloCluster *aodcluster = static_cast<AliAODCaloCluster*>(fClusters->At(i));
862  if (!aodcluster || !aodcluster->IsEMCAL())
863  continue;
864  AliAODCaloCluster *clus = new ((*fOutClusters)[nCopiedClusters]) AliAODCaloCluster(*aodcluster);
865  if (!fIsMC)
866  clus->SetLabel(0,0);
867  nCopiedClusters++;
868  }
869  }
870 }
871 
872 //________________________________________________________________________
874 {
875  // Copy all the tracks in the new collection
876 
877  if (!fTracks)
878  return;
879 
880  const Int_t nTracks = fTracks->GetEntriesFast();
881  Int_t nCopiedTracks = 0;
882  for (Int_t i = 0; i < nTracks; ++i) {
883  AliPicoTrack *picotrack = static_cast<AliPicoTrack*>(fTracks->At(i));
884  if (!picotrack)
885  continue;
886  AliPicoTrack *track = new ((*fOutTracks)[nCopiedTracks]) AliPicoTrack(*picotrack);
887  if (!fIsMC && track->GetLabel() != 0)
888  track->SetLabel(0);
889  nCopiedTracks++;
890  }
891 }
892 
893 //________________________________________________________________________
895 {
896  // Copy all the MC particles in the new collection
897 
898  if (!fMCParticles)
899  return;
900 
901  const Int_t nPart = fMCParticles->GetEntriesFast();
902  Int_t nCopiedPart = 0;
903  for (Int_t i = 0; i < nPart; ++i) {
904  AliAODMCParticle *part = static_cast<AliAODMCParticle*>(fMCParticles->At(i));
905  if (!part)
906  continue;
907  new ((*fOutMCParticles)[nCopiedPart]) AliAODMCParticle(*part);
908 
909  nCopiedPart++;
910  }
911 
913  return;
914 
915  if (fOutMCParticlesMap->GetSize() < fMCParticlesMap->GetSize())
916  fOutMCParticlesMap->Set(fMCParticlesMap->GetSize() * 2);
917 
918  for (Int_t i = 0; i < fMCParticlesMap->GetSize(); i++) {
919  fOutMCParticlesMap->AddAt(fMCParticlesMap->At(i), i);
920  if (fMCParticlesMap->At(i) >= 0)
921  fMCLabelShift = i;
922  }
923 
924  AliDebug(2,Form("MC particles copied. fMCLabelShift=%d",fMCLabelShift));
925 }
926 
927 //________________________________________________________________________
929 {
930  // Get random cell.
931 
932  Int_t repeats = 0;
933  Double_t rndEta = eta;
934  Double_t rndPhi = phi;
935  do {
936  if (eta < -100)
937  rndEta = GetRandomEta(kTRUE);
938  if (phi < 0)
939  rndPhi = GetRandomPhi(kTRUE);
940  fGeom->GetAbsCellIdFromEtaPhi(rndEta, rndPhi, absId);
941  repeats++;
942  } while (absId == -1 && repeats < 100);
943 
944  if (!(absId > -1)) {
945  AliWarning(Form("Could not extract random cluster! Random eta-phi extracted more than 100 times!\n"
946  "eta [%f, %f], phi [%f, %f]\n", fEtaMin, fEtaMax, fPhiMin, fPhiMax));
947  }
948  else {
949  eta = rndEta;
950  phi = rndPhi;
951  }
952 }
953 
954 //________________________________________________________________________
956 {
957  // Get random eta.
958 
961 
962  if (emcal) {
963  const Double_t EmcalMinEta = fGeom->GetArm1EtaMin();
964  const Double_t EmcalMaxEta = fGeom->GetArm1EtaMax();
965 
966  if (etamax > EmcalMaxEta) etamax = EmcalMaxEta;
967  if (etamax < EmcalMinEta) etamax = EmcalMinEta;
968  if (etamin > EmcalMaxEta) etamin = EmcalMaxEta;
969  if (etamin < EmcalMinEta) etamin = EmcalMinEta;
970  }
971 
972  return gRandom->Rndm() * (etamax - etamin) + etamin;
973 }
974 
975 //________________________________________________________________________
977 {
978  // Get random phi.
979 
980  Double_t phimax = fPhiMax;
982 
983  if (emcal) {
984  const Double_t EmcalMinPhi = fGeom->GetArm1PhiMin() * TMath::DegToRad();
985  const Double_t EmcalMaxPhi = fGeom->GetArm1PhiMax() * TMath::DegToRad();
986 
987  if (phimax > EmcalMaxPhi) phimax = EmcalMaxPhi;
988  if (phimax < EmcalMinPhi) phimax = EmcalMinPhi;
989  if (phimin > EmcalMaxPhi) phimin = EmcalMaxPhi;
990  if (phimin < EmcalMinPhi) phimin = EmcalMinPhi;
991  }
992 
993  Double_t result = gRandom->Rndm() * (phimax - phimin) + phimin;
994 
995  return result;
996 }
997 
998 //________________________________________________________________________
1000 {
1001  // Get random pt.
1002 
1003  if (fPtSpectrum)
1004  return fPtSpectrum->GetRandom();
1005  else
1006  return gRandom->Rndm() * (fPtMax - fPtMin) + fPtMin;
1007 }
1008 
1009 //________________________________________________________________________
1010 
1012 
1013  Double_t m = 0;
1014 
1015  if(fMassFromDistr) {
1016  if(fHMassDistrib) m = fHMassDistrib->GetRandom();
1017  else {
1018  AliError("Template distribution for mass of track embedding not found, use 0");
1019  m = 0;
1020  }
1021  }
1022  return m;
1023 }
1024 
1025 //________________________________________________________________________
1026 
1028 
1029  Double_t maxedgex = fHMassPtDistrib->GetXaxis()->GetBinLowEdge(fHMassPtDistrib->GetNbinsX());
1030  Double_t maxedgey = fHMassPtDistrib->GetYaxis()->GetBinLowEdge(fHMassPtDistrib->GetNbinsY());
1031 
1032  if(maxedgex < maxedgey) //condition assuming that the mass axis has smaller range than the pT axis!!
1033  fHMassPtDistrib->GetRandom2(m, pt);
1034  else fHMassPtDistrib->GetRandom2(pt, m);
1035 
1036 }
1037 
1038 //________________________________________________________________________
1040 {
1041  // Get a random particle.
1042 
1043  eta = GetRandomEta(emcal);
1044 
1045  if (fPtPhiEvPlDistribution) {
1046  Double_t phimax = fPhiMax;
1048 
1049  if (emcal) {
1050  const Double_t EmcalMinPhi = fGeom->GetArm1PhiMin() * TMath::DegToRad();
1051  const Double_t EmcalMaxPhi = fGeom->GetArm1PhiMax() * TMath::DegToRad();
1052 
1053  if (phimax > EmcalMaxPhi) phimax = EmcalMaxPhi;
1054  if (phimax < EmcalMinPhi) phimax = EmcalMinPhi;
1055  if (phimin > EmcalMaxPhi) phimin = EmcalMaxPhi;
1056  if (phimin < EmcalMinPhi) phimin = EmcalMinPhi;
1057  }
1058 
1059  if (fPtPhiEvPlDistribution->GetXmin() > phimax || fPtPhiEvPlDistribution->GetXmax() < phimin) {
1060  AliWarning(Form("The hisogram %s does not overlap with the EMCal acceptance limits. It will be ignored.",fPtPhiEvPlDistribution->GetName()));
1061  pt = GetRandomPt();
1062  phi = GetRandomPhi(emcal);
1063  }
1064  else {
1065  do {
1066  fPtPhiEvPlDistribution->GetRandom2(pt,phi);
1067  phi += fPsi;
1068  if (phi > TMath::Pi() * 2) phi -= TMath::Pi() * 2;
1069  } while (phi > phimax || phi < phimin);
1070  }
1071  }
1072  else {
1073  pt = GetRandomPt();
1074  phi = GetRandomPhi(emcal);
1075  }
1076 }
1077 
1078 //________________________________________________________________________
1080 
1083 
1084  phi = GetRandomPhi(emcal);
1085  eta = GetRandomEta(emcal);
1086  GetRandomMvsPt(m, pt);
1087 
1088 }
1089 //________________________________________________________________________
1091 
1092  if(!fHMassPtDistrib) {
1093  GetRandomParticle(pt,eta,phi,emcal);
1094  m = GetRandomM();
1095 
1096  } else GetRandomMvsPtParticle(pt,m,eta,phi,emcal);
1097 
1098 }
1099 
1100 //________________________________________________________________________
1102 {
1103  // Run.
1104 }
1105 //________________________________________________________________________
1107 
1108  if(!fhpTEmb || !fhMEmb || !fhEtaEmb || !fhPhiEmb) {
1109  AliError("Histograms not found, are the QA histograms active?");
1110  }
1111  fhEvents->Fill(0);
1112  // fill histograms
1113  Int_t nentries = fOutTracks->GetEntries();
1114  for(Int_t it = 0; it<fNTracks; it++){
1115  AliPicoTrack *trackEmb = dynamic_cast<AliPicoTrack*>(fOutTracks->At(nentries-it-1));
1116  if(!trackEmb) continue;
1117  if(trackEmb->GetLabel() >= fMarkMC+fMCLabelShift){
1118  fhpTEmb ->Fill(trackEmb->Pt());
1119  fhMEmb ->Fill(trackEmb->M());
1120  fhEtaEmb->Fill(trackEmb->Eta());
1121  fhPhiEmb->Fill(trackEmb->Phi());
1122  }
1123  }
1124 
1125  PostData(1, fOutput);
1126 
1127 }
1128 
1129 //________________________________________________________________________
1130 
1132  if(!hM){
1133  AliError("Null histogram for mass distribution");
1134  return;
1135  }
1136  fMassFromDistr = kTRUE;
1137  fHMassDistrib = hM;
1138  AliInfo("Input mass distribution set");
1139 
1140  return;
1141 }
1142 
1143 //________________________________________________________________________
1144 
1146  if(!hmasspt){
1147  AliError("Null histogram for mass vs pt distribution");
1148  return;
1149  }
1150  fMassFromDistr = kTRUE;
1151  fHMassPtDistrib = hmasspt;
1152  AliInfo("Input mass vs pt distribution set");
1153 
1154  return;
1155 }
1156 
1157 //________________________________________________________________________
1159 
1160  SetDistributionFromFile(filename, histoname, 1);
1161 }
1162 
1163 //________________________________________________________________________
1165  SetDistributionFromFile(filename, histoname, 2);
1166 }
1167 
1168 //________________________________________________________________________
1170  SetDistributionFromFile(filename, histoname, 3);
1171 }
1172 
1173 //________________________________________________________________________
1175 
1176  if(filename.Contains("alien")) {
1177  TGrid::Connect("alien://");
1178  }
1179  TFile *f = TFile::Open(filename);
1180  if(!f){
1181  AliFatal(Form("File %s not found, cannot SetMassDistribution", filename.Data()));
1182  return;
1183  }
1184 
1185  TH1F* h = 0x0;
1186  TH2F* g = 0x0;
1187  if(type < 3){
1188  h = dynamic_cast<TH1F*> (f->Get(histoname));
1189  if(!h) {
1190  AliError("Input file for Mass not found");
1191  f->ls();
1192  }
1193  }
1194  if(type == 3){
1195  g = dynamic_cast<TH2F*> (f->Get(histoname));
1196  if(!g) {
1197  AliError("Input file for Mass not found");
1198  f->ls();
1199  }
1200  }
1201 
1202  if(type == 1) SetPtSpectrum(h);
1203  if(type == 2) SetMassDistribution(h);
1204  if(type == 3) SetMassVsPtDistribution(g);
1205 
1206  //f->Close();
1207  //delete f;
1208 
1209  return;
1210 
1211 }
1212 
1213 //________________________________________________________________________
1214 
1215 void AliJetModelBaseTask::SetMassAndPtDistributionFromFile(TString filenameM, TString filenamepT, TString histonameM, TString histonamepT){
1216  SetMassDistributionFromFile(filenameM, histonameM);
1217  SetpTDistributionFromFile(filenamepT, histonamepT);
1218  return;
1219 }
1220 
1221 
TH1F * fDensitySpectrum
particle density spectrum to extract random density values
Int_t charge
const char * filename
Definition: TestFCM.C:1
double Double_t
Definition: External.C:58
virtual void Run()
intialize task
TClonesArray * fClusters
! cluster collection
Double_t GetRandomPt()
generate a random phi value in the given range
Definition: External.C:236
void SetpTDistributionFromFile(TString filename, TString histoname)
TString fOutCellsName
name of output cells collection
Float_t fPtMin
pt minimum value
Double_t mass
TList * fOutput
! output list for QA histograms
TClonesArray * fOutClusters
! output cluster collection
Int_t fMarkMC
which MC label is to be used (default=100)
TString fOutTracksName
name of output track collection
void SetMassDistribution(TH1F *hM)
void SetMassVsPtDistribution(TH2F *hmasspt)
Bool_t fIsInit
! =true if initialized
TString fMCParticlesName
name of MC particle collection
AliNamedArrayI * fMCParticlesMap
! MC particles mapping
void GetRandomMassiveParticle(Double_t &pt, Double_t &eta, Double_t &phi, Bool_t emcal, Double_t &m)
generate a particle with random eta,phi,pt values
Base class for embedding into an event.
TString GetOutTrackName() const
Bool_t fIsMC
whether the current event is MC or not
TString fCaloName
name of calo cluster collection
TString fCellsName
name of calo cells collection
TRandom * gRandom
AliVCaloCells * fCaloCells
! cells collection
Float_t fEtaMax
eta maximum value
TClonesArray * fOutTracks
! output track collection
TH1I * fhEvents
! store the number of events analysed
const Double_t etamin
TClonesArray * fOutMCParticles
! output MC particles collection
Int_t GetLabel() const
Definition: AliPicoTrack.h:41
Double_t Eta() const
Definition: AliPicoTrack.h:37
Int_t fNTracks
how many tracks are being processed
Float_t fPhiMin
phi minimum value
int Int_t
Definition: External.C:63
Bool_t fEsdMode
! ESD/AOD mode
Bool_t fAddV2
add v2 sampled from a tf1
Double_t GetRandomM()
generate a random pt value in the given range
Definition: External.C:204
void GetRandomParticle(Double_t &pt, Double_t &eta, Double_t &phi, Bool_t emcal=kFALSE)
generate a random m value from a given distribution or take a fixed value
unsigned int UInt_t
Definition: External.C:33
float Float_t
Definition: External.C:68
TString fTracksName
name of track collection
TH1F * fhpTEmb
! embedded tracks pT
Bool_t fFlowFluctuations
introduce gaussian flow fluctuation
TString fOutMCParticlesName
name of output MC particle collection
void GetRandomCell(Double_t &eta, Double_t &phi, Int_t &absId)
TString fOutCaloName
name of output cluster collection
Double_t M() const
Definition: AliPicoTrack.h:36
Int_t fNClusters
how many clusters are being processed
Double_t Phi() const
Definition: AliPicoTrack.h:33
Int_t fMCLabelShift
! MC label shift
void SetLabel(Int_t label)
Definition: AliPicoTrack.h:42
void SetDistributionFromFile(TString filename, TString histoname, Int_t type)
AliAODMCParticle * AddMCParticle(AliAODMCParticle *part, Int_t origIndex)
add a track; if values are -1 generate random parameters
AliVCaloCells * fOutCaloCells
! output cells collection
Bool_t fQAhistos
draw QA histograms
TString fSuffix
suffix to add in the name of new collections
TH1F * fhPhiEmb
! embedded tracks phi
ClassImp(AliJetModelBaseTask) AliJetModelBaseTask
Double_t Pt() const
Definition: AliPicoTrack.h:23
Bool_t fCopyArray
whether or not the array will be copied to a new one before modelling
Int_t SetNumberOfOutCells(Int_t n)
void GetRandomMvsPt(Double_t &m, Double_t &pt)
generate a particle with random eta,phi,pt,mass values
TClonesArray * fMCParticles
! MC particles collection
Int_t AddCell(Double_t e=-1, Double_t eta=-999, Double_t phi=-1)
set the number of cells
short Short_t
Definition: External.C:23
Double_t GetRandomPhi(Bool_t emcal=kFALSE)
generate a random eta value in the given range
TH2F * fHMassPtDistrib
shape of mass vs pt distribution of embedded track
TF2 * fPtPhiEvPlDistribution
pt vs. (phi-psi) distribution to extract random pt/phi values
void FillHistograms()
do jet model action
AliNamedArrayI * fOutMCParticlesMap
! MC particles mapping
Int_t fAddedCells
! number of added cells
Double_t fPsi
! simmetry plane for the elliptic flow
TH1F * fhEtaEmb
! embedded tracks eta
AliVCluster * AddCluster(Double_t e=-1, Double_t eta=-999, Double_t phi=-1, Int_t label=0)
add a cell with given energy, position and times
TF1 * fDifferentialV2
v2 as function of pt
Bool_t fMassFromDistr
draw the particle mass from fHMassDistrib
TH1F * fhMEmb
! embedded tracks M
Float_t fEtaMin
eta minimum value
void Clear(Option_t *option="")
TH1F * fPtSpectrum
pt spectrum to extract random pt values
const Double_t etamax
TString fGeomName
Fill QA histograms.
unsigned short UShort_t
Definition: External.C:28
TH1F * fHMassDistrib
shape of mass distribution of embedded tracks
const char Option_t
Definition: External.C:48
Double_t fVertex[3]
! event vertex
bool Bool_t
Definition: External.C:53
void AddV2(Double_t &phi, Double_t &pt) const
AliPicoTrack * AddTrack(Double_t pt=-999, Double_t eta=-999, Double_t phi=-999, Byte_t type=0, Double_t etaemc=0, Double_t phiemc=0, Double_t ptemc=0, Bool_t ise=kFALSE, Int_t label=0, Short_t charge=1, Double_t mass=0.1396)
add a cluster (copy)
void SetMassDistributionFromFile(TString filename, TString histoname)
Float_t fPhiMax
phi maximum value
void SetMassAndPtDistributionFromFile(TString filenameM, TString filenamepT, TString histonameM, TString histonamepT)
TClonesArray * fTracks
! track collection
void GetRandomMvsPtParticle(Double_t &pt, Double_t &m, Double_t &eta, Double_t &phi, Bool_t emcal=kFALSE)
generate 2 random values for pt and mass from a gived 2D distribution
AliEMCALGeometry * fGeom
! pointer to EMCal geometry
virtual Bool_t ExecOnce()
generate a particle with random eta,phi, and correlated pt,mass values
void SetPtSpectrum(TH1F *f)
Double_t GetRandomEta(Bool_t emcal=kFALSE)
generate a random cell in the calorimeter
Int_t fNCells
how many cells are being processed
const Double_t phimin
TList * OpenFile(const char *fname)
Definition: DrawAnaELoss.C:65
void SetMassVsPtDistributionFromFile(TString filename, TString histoname)
Float_t fPtMax
pt maximum value