AliPhysics  32e057f (32e057f)
AliCaloTrackMCReader.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 // --- ROOT system ---
17 #include <TParticle.h>
18 #include <TDatabasePDG.h>
19 #include <TRandom.h>
20 #include <TArrayI.h>
21 #include "TParticle.h"
22 
23 //---- ANALYSIS system ----
24 #include "AliCaloTrackMCReader.h"
25 #include "AliGenEventHeader.h"
26 #include "AliMCEvent.h"
27 #include "AliVCluster.h"
28 #include "AliAODTrack.h"
29 #include "AliAODEvent.h"
30 #include "AliFiducialCut.h"
31 #include "AliMCAnalysisUtils.h"
32 #include "AliLog.h"
33 
35 ClassImp(AliCaloTrackMCReader) ;
37 
38 //____________________________________________
40 //____________________________________________
42 AliCaloTrackReader(), fDecayPi0(0),
43 fNeutralParticlesArray(0x0), fChargedParticlesArray(0x0),
44 fStatusArray(0x0), fKeepAllStatus(0),
45 fCheckOverlap(0), fEMCALOverlapAngle(0),fPHOSOverlapAngle(0),
46 fIndex2ndPhoton(0), fOnlyGeneratorParticles(kTRUE),
47 fMomentum(), fPi0Momentum(),
48 fGamDecayMom1(), fGamDecayMom2()
49 {
51 }
52 
53 //___________________________________________
55 //___________________________________________
57 {
59 
62  if(fStatusArray) delete fStatusArray ;
63 }
64 
65 //________________________________________________________
67 //________________________________________________________
69 {
70  TArrayF pv;
71  GetGenEventHeader()->PrimaryVertex(pv);
72  v[0]=pv.At(0);
73  v[1]=pv.At(1);
74  v[2]=pv.At(2);
75 }
76 
77 //_________________________________________
79 //_________________________________________
81 {
82  fDecayPi0 = kFALSE;
83 
85  fChargedParticlesArray->SetAt(11,0);
86  //Int_t pdgarray[]={12,14,16};// skip neutrinos
87  //fNeutralParticlesArray = new TArrayI(3, pdgarray);
89  fNeutralParticlesArray->SetAt(12,0); fNeutralParticlesArray->SetAt(14,1); fNeutralParticlesArray->SetAt(16,2);
90  fStatusArray = new TArrayI(1);
91  fStatusArray->SetAt(1,0);
92 
94  fKeepAllStatus = kTRUE;
95 
96  fCheckOverlap = kFALSE;
97  fEMCALOverlapAngle = 2.5 * TMath::DegToRad();
98  fPHOSOverlapAngle = 0.5 * TMath::DegToRad();
99  fIndex2ndPhoton = -1;
100 
101  fDataType = kMC;
102 
103  //For this reader we own the objects of the arrays
104  fCTSTracks ->SetOwner(kTRUE);
105  fEMCALClusters->SetOwner(kTRUE);
106  fPHOSClusters ->SetOwner(kTRUE);
107 }
108 
109 //____________________________________________________________________________________
111 //____________________________________________________________________________________
113  Int_t & iPrimary, Int_t & index, Int_t & pdg)
114 {
115  if( fIndex2ndPhoton==iPrimary )
116  {
117  fIndex2ndPhoton=-1;
118  return;
119  }
120  else
121  fIndex2ndPhoton=-1;
122 
123  if(pdg!=22) return;
124 
125  AliVParticle *meson = GetMC()->GetTrack(imom);
126  Int_t mepdg = meson->PdgCode();
127  Int_t idaug1 = meson->GetFirstDaughter();
128  if((mepdg == 111 || mepdg == 221 ) && meson->GetNDaughters() == 2)
129  { //Check only decay in 2 photons
130  AliVParticle * d1 = GetMC()->GetTrack(idaug1 );
131  AliVParticle * d2 = GetMC()->GetTrack(idaug1+1);
132  if(d1->PdgCode() == 22 && d2->PdgCode() == 22 )
133  {
134  d1->Momentum(fGamDecayMom1);
135  d2->Momentum(fGamDecayMom2);
136  //printf("angle %2.2f\n",ph1.Angle(ph2.Vect()));
137 
138  if(anglethres > fGamDecayMom1.Angle(fGamDecayMom2.Vect()))
139  {
140  //Keep the meson
141  pdg=mepdg;
142  index=imom;
143  meson->Momentum(fMomentum);
144  //printf("Overlap:: pt %2.2f, phi %2.2f, eta %2.2f\n",mom.Pt(),mom.Phi(),mom.Eta());
145  if(iPrimary == idaug1) iPrimary++; //skip next photon in list
146  }
147  else
148  {
149  //Do not check overlapping for next decay photon from same meson
150  if(iPrimary == idaug1)
151  {
152  fIndex2ndPhoton = idaug1+1;
153  }
154 
155  }
156  }
157  } // Meson Decay with 2 photon daughters
158 }
159 
160 //__________________________________________________________________________________
162 //__________________________________________________________________________________
164 {
165  Char_t ttype = 0;
166  Float_t overAngleLimit = 100;
167 
168  if (fFillPHOS)
169  {
170  if( fMomentum.Pt() < fPHOSPtMin ) return;
171  if( fCheckFidCut && !fFiducialCut->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kPHOS )) return;
172  ttype = AliVCluster::kPHOSNeutral;
173  overAngleLimit = fPHOSOverlapAngle;
174  }
175  else if(fFillEMCAL)
176  {
177  if( fMomentum.Pt() < fEMCALPtMin ) return;
178  if( fCheckFidCut && !fFiducialCut->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kEMCAL)) return;
179  ttype= AliVCluster::kEMCALClusterv1;
180  overAngleLimit = fEMCALOverlapAngle;
181  }
182 
183  Int_t index = iParticle ;
184  if(fCheckOverlap)
185  CheckOverlap(overAngleLimit,motherIndex,index, iParticle, pdg);
186 
187  Int_t labels[] = {index};
188  Double_t x[] = {fMomentum.X(), fMomentum.Y(), fMomentum.Z()};
189 
190  // Create object and write it to file
191  AliAODCaloCluster *calo = new AliAODCaloCluster(index,1,labels,fMomentum.E(), x, NULL, ttype, 0);
192 
193  SetCaloClusterPID(pdg,calo) ;
194 
195  AliDebug(3,Form("PHOS %d?, EMCAL %d? : Selected cluster pdg %d, E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f",
196  fFillPHOS, fFillEMCAL, pdg,fMomentum.E(), fMomentum.Pt(), fMomentum.Phi()*TMath::RadToDeg(),fMomentum.Eta()));
197 
198  // Reference the selected object to the list
199  if(fFillPHOS)fPHOSClusters ->Add(calo);
200  else fEMCALClusters->Add(calo);
201 }
202 
203 //___________________________________________________________________________
205 //___________________________________________________________________________
207  const char * /*currentFileName*/)
208 {
209  fEventNumber = iEntry;
210  //fCurrentFileName = TString(currentFileName);
211 
212  for(Int_t iptCut = 0; iptCut < fTrackMultNPtCut; iptCut++ )
213  {
214  fTrackMult [iptCut] = 0;
215  fTrackSumPt[iptCut] = 0;
216  }
217 
218  // In case of analysis of events with jets, skip those with jet pt > 5 pt hard
219  if(fComparePtHardAndJetPt && GetMC())
220  {
221  if(!ComparePtHardAndJetPt()) return kFALSE ;
222  }
223 
224  // Fill Vertex array
225  FillVertexArray();
226 
227  Int_t iParticle = 0 ;
228  Int_t nparticles = GetMC()->GetNumberOfTracks() ;
229 
230  if(fOnlyGeneratorParticles) nparticles=GetMC()->GetNumberOfPrimaries();
231 
232  for (iParticle = 0 ; iParticle < nparticles ; iParticle++)
233  {
234  AliVParticle * particle = GetMC()->GetTrack(iParticle);
235  Float_t p[3];
236  Float_t x[3];
237 
238  //Keep particles with a given status
239  if(KeepParticleWithStatus(particle->MCStatusCode()) && (particle->Pt() > 0) )
240  {
241  //Skip bizarre particles, they crash when charge is calculated
242  // if(TMath::Abs(pdg) == 3124 || TMath::Abs(pdg) > 10000000) continue ;
243 
244  //charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
245  Int_t charge = particle->Charge();
246  Int_t momIndex = particle->GetMother();
247  Int_t pdg = particle->PdgCode();
248 
249  particle->Momentum(fMomentum);
250 
251  Float_t en = particle->E();
252  Float_t pt = particle->Pt();
253  Float_t eta = particle->Eta();
254  Float_t phi = particle->Phi();
255 
256  //---------- Charged particles ----------------------
257  if(charge != 0)
258  {
259 // if(TMath::Abs(eta)< fTrackMultEtaCut) fTrackMult++;
260  if(TMath::Abs(eta)< fTrackMultEtaCut)
261  {
262  for(Int_t iptCut = 0; iptCut < fTrackMultNPtCut; iptCut++ )
263  {
264  if(pt > fTrackMultPtCut[iptCut])
265  {
266  fTrackMult [iptCut]++;
267  fTrackSumPt[iptCut]+=pt;
268  }
269  }
270  }
271 
272 // if(fFillCTS && (pt > fCTSPtMin))
273  if(fFillCTS && (fCTSPtMin > pt || fCTSPtMax < pt)) continue ;
274  {
275  // Particles in CTS acceptance
276 
277  if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(eta,phi,kCTS)) continue;
278 
279  if(TMath::Abs(pdg) == 11 && (GetMC()->GetTrack(particle->GetMother()))->PdgCode()==22) continue ;
280 
281  AliDebug(2,Form("CTS : Selected tracks E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f",
282  en,pt,phi*TMath::RadToDeg(),eta));
283 
284  x[0] = particle->Xv(); x[1] = particle->Yv(); x[2] = particle->Zv();
285  p[0] = particle->Px(); p[1] = particle->Py(); p[2] = particle->Pz();
286  //Create object and write it to file
287  AliAODTrack *aodTrack = new AliAODTrack(0, iParticle, p, kTRUE, x, kFALSE,NULL, 0, 0,
288  // NULL,
289  0x0,//primary,
290  kFALSE, // No fit performed
291  kFALSE, // No fit performed
293  0);
294  SetTrackChargeAndPID(pdg, aodTrack);
295  fCTSTracks->Add(aodTrack);//reference the selected object to the list
296  }
297 
298  // Keep some charged particles in calorimeters lists
299  if((fFillPHOS || fFillEMCAL) && KeepChargedParticles(pdg))
300  FillCalorimeters(iParticle, momIndex, pdg);
301 
302  }//Charged
303 
304  //-------------Neutral particles ----------------------
305  else if(charge == 0 && (fFillPHOS || fFillEMCAL))
306  {
307  //Skip neutrinos or other neutral particles
308  //if(SkipNeutralParticles(pdg) || particle->IsPrimary()) continue ; // protection added (MG)
309  if(SkipNeutralParticles(pdg)) continue ;
310  //Fill particle/calocluster arrays
311  if(!fDecayPi0)
312  {
313  AliDebug(2,Form("Calo : Selected tracks E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f",
314  en,pt,phi*TMath::RadToDeg(),eta));
315  FillCalorimeters(iParticle, momIndex, pdg);
316  }
317  else
318  {
319  //Sometimes pi0 are stable for the generator, if needed decay it by hand
320  if(pdg == 111 )
321  {
322  if(pt > fPHOSPtMin || pt > fEMCALPtMin)
323  {
324  //Decay
325  //Double_t angle = 0;
326  particle->Momentum(fPi0Momentum);
327 
328  MakePi0Decay();//,angle);
329 
330  // Check if Pi0 is in the acceptance of the calorimeters, if aperture angle is small, keep it
331  // and fill particle/calocluster arrays
332  fMomentum.SetPxPyPzE(fGamDecayMom1.Px(),fGamDecayMom1.Py(), fGamDecayMom1.Pz(),fGamDecayMom1.E());
333  FillCalorimeters(iParticle,iParticle,22);
334 
335  fMomentum.SetPxPyPzE(fGamDecayMom2.Px(),fGamDecayMom2.Py(), fGamDecayMom2.Pz(),fGamDecayMom2.E());
336  FillCalorimeters(iParticle,iParticle,22);
337  }//pt cut
338  }//pi0
339  else
340  {
341  AliDebug(2,Form("Calo : Selected tracks E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f",
342  en,pt,phi*TMath::RadToDeg(),eta));
343  FillCalorimeters(iParticle,momIndex,pdg); //Add the rest
344  }
345  }
346  }//neutral particles
347  } //particle with correct status
348  }//particle loop
349 
350  fIndex2ndPhoton = -1; //In case of overlapping studies, reset for each event
351 
352  return kTRUE;
353 }
354 
355 //_____________________________________________________________________
358 //_____________________________________________________________________
360 {
361  if(!fKeepAllStatus)
362  {
363  for(Int_t i= 0; i < fStatusArray->GetSize(); i++)
364  if(status == fStatusArray->At(i)) return kTRUE ;
365 
366  return kFALSE;
367  }
368  else
369  return kTRUE ;
370 }
371 
372 //________________________________________________________________
375 //________________________________________________________________
377 {
378  for(Int_t i= 0; i < fChargedParticlesArray->GetSize(); i++)
379  {
380  if(TMath::Abs(pdg) == fChargedParticlesArray->At(i)) return kTRUE ;
381  }
382 
383  return kFALSE;
384 }
385 
386 //__________________________________________________________
388 //__________________________________________________________
389 void AliCaloTrackMCReader::Print(const Option_t * opt) const
390 {
391  if(! opt)
392  return;
393 
395 
396  printf("**** Print **** %s %s ****\n", GetName(), GetTitle() ) ;
397 
398  printf("Decay Pi0? : %d\n", fDecayPi0) ;
399  printf("Check Overlap in Calo? : %d\n", fCheckOverlap) ;
400  printf("Keep all status? : %d\n", fKeepAllStatus) ;
401 
402  if(!fKeepAllStatus) printf("Keep particles with status : ");
403  for(Int_t i= 0; i < fStatusArray->GetSize(); i++)
404  printf(" %d ; ", fStatusArray->At(i));
405  printf("\n");
406 
407  printf("Skip neutral particles in calo : ");
408  for(Int_t i= 0; i < fNeutralParticlesArray->GetSize(); i++)
409  printf(" %d ; ", fNeutralParticlesArray->At(i));
410  printf("\n");
411 
412  printf("Keep charged particles in calo : ");
413  for(Int_t i= 0; i < fChargedParticlesArray->GetSize(); i++)
414  printf(" %d ; ", fChargedParticlesArray->At(i));
415  printf("\n");
416 }
417 
418 //_______________________________________
422 //_______________________________________
423 void AliCaloTrackMCReader::MakePi0Decay() //, Double_t &angle)
424 {
425  // cout<<"Boost vector"<<endl;
426  Double_t mPi0 = TDatabasePDG::Instance()->GetParticle(111)->Mass();
427  TVector3 b = fPi0Momentum.BoostVector();
428  //cout<<"Parameters"<<endl;
429  //Double_t mPi0 = fPi0Momentum.M();
430  Double_t phi = TMath::TwoPi() * gRandom->Rndm();
431  Double_t cosThe = 2 * gRandom->Rndm() - 1;
432  Double_t cosPhi = TMath::Cos(phi);
433  Double_t sinPhi = TMath::Sin(phi);
434  Double_t sinThe = TMath::Sqrt(1-cosThe*cosThe);
435  Double_t ePi0 = mPi0/2.;
436  //cout<<"ePi0 "<<ePi0<<endl;
437  //cout<<"Components"<<endl;
438  fGamDecayMom1.SetPx(+ePi0*cosPhi*sinThe);
439  fGamDecayMom1.SetPy(+ePi0*sinPhi*sinThe);
440  fGamDecayMom1.SetPz(+ePi0*cosThe);
441  fGamDecayMom1.SetE(ePi0);
442  //cout<<"fGamDecayMom1: "<<fGamDecayMom1.Px()<<" "<<fGamDecayMom1.Py()<<" "<<fGamDecayMom1.Pz()<<" "<<fGamDecayMom1.E()<<endl;
443  //cout<<"fGamDecayMom1 Mass: "<<fGamDecayMom1.Px()*fGamDecayMom1.Px()+fGamDecayMom1.Py()*fGamDecayMom1.Py()+fGamDecayMom1.Pz()*fGamDecayMom1.Pz()-fGamDecayMom1.E()*fGamDecayMom1.E()<<endl;
444  fGamDecayMom2.SetPx(-ePi0*cosPhi*sinThe);
445  fGamDecayMom2.SetPy(-ePi0*sinPhi*sinThe);
446  fGamDecayMom2.SetPz(-ePi0*cosThe);
447  fGamDecayMom2.SetE(ePi0);
448  //cout<<"fGamDecayMom2: "<<fGamDecayMom2.Px()<<" "<<fGamDecayMom2.Py()<<" "<<fGamDecayMom2.Pz()<<" "<<fGamDecayMom2.E()<<endl;
449  //cout<<"fGamDecayMom2 Mass: "<<fGamDecayMom2.Px()*fGamDecayMom2.Px()+fGamDecayMom2.Py()*fGamDecayMom2.Py()+fGamDecayMom2.Pz()*fGamDecayMom2.Pz()-fGamDecayMom2.E()*fGamDecayMom2.E()<<endl;
450  //cout<<"Boost "<<b.X()<<" "<<b.Y()<<" "<<b.Z()<<endl;
451  fGamDecayMom1.Boost(b);
452  //cout<<"fGamDecayMom1: "<<fGamDecayMom1.Px()<<" "<<fGamDecayMom1.Py()<<" "<<fGamDecayMom1.Pz()<<" "<<fGamDecayMom1.E()<<endl;
453  fGamDecayMom2.Boost(b);
454  //cout<<"fGamDecayMom2: "<<fGamDecayMom2.Px()<<" "<<fGamDecayMom2.Py()<<" "<<fGamDecayMom2.Pz()<<" "<<fGamDecayMom2.E()<<endl;
455  //cout<<"angle"<<endl;
456  //angle = fGamDecayMom1.Angle(fGamDecayMom2.Vect());
457  //cout<<angle<<endl;
458 }
459 
460 //__________________________________________________________________
462 //__________________________________________________________________
464  AliAODEvent* aod,
465  AliMCEvent* mc)
466 {
467  SetMC(mc);
468  SetOutputEvent(aod);
469 }
470 
471 //________________________________________________________________
474 //________________________________________________________________
476 {
477  for(Int_t i= 0; i < fNeutralParticlesArray->GetSize(); i++)
478  {
479  if(TMath::Abs(pdg) == fNeutralParticlesArray->At(i)) return kTRUE ;
480  }
481 
482  return kFALSE;
483 }
484 
485 //_______________________________________________________________________
487 //_______________________________________________________________________
489  AliAODTrack *track) const
490 {
491  Float_t pid[10] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
492 
493  switch (pdgCode)
494  {
495  case 22: // gamma
496  track->SetCharge(0);
497  pid[AliAODTrack::kUnknown] = 1.;
498  track->SetPID(pid);
499  break;
500 
501  case 11: // e-
502  track->SetCharge(-1);
503  pid[AliAODTrack::kElectron] = 1.;
504  track->SetPID(pid);
505  break;
506 
507  case -11: // e+
508  track->SetCharge(+1);
509  pid[AliAODTrack::kElectron] = 1.;
510  track->SetPID(pid);
511  break;
512 
513  case 13: // mu-
514  track->SetCharge(-1);
515  pid[AliAODTrack::kMuon] = 1.;
516  track->SetPID(pid);
517  break;
518 
519  case -13: // mu+
520  track->SetCharge(+1);
521  pid[AliAODTrack::kMuon] = 1.;
522  track->SetPID(pid);
523  break;
524 
525  case 111: // pi0
526  track->SetCharge(0);
527  pid[AliAODTrack::kUnknown] = 1.;
528  track->SetPID(pid);
529  break;
530 
531  case 211: // pi+
532  track->SetCharge(+1);
533  pid[AliAODTrack::kPion] = 1.;
534  track->SetPID(pid);
535  break;
536 
537  case -211: // pi-
538  track->SetCharge(-1);
539  pid[AliAODTrack::kPion] = 1.;
540  track->SetPID(pid);
541  break;
542 
543  case 130: // K0L
544  track->SetCharge(0);
545  pid[AliAODTrack::kUnknown] = 1.;
546  track->SetPID(pid);
547  break;
548 
549  case 321: // K+
550  track->SetCharge(+1);
551  pid[AliAODTrack::kKaon] = 1.;
552  track->SetPID(pid);
553  break;
554 
555  case -321: // K-
556  track->SetCharge(-1);
557  pid[AliAODTrack::kKaon] = 1.;
558  track->SetPID(pid);
559  break;
560 
561  case 2112: // n
562  track->SetCharge(0);
563  pid[AliAODTrack::kUnknown] = 1.;
564  track->SetPID(pid);
565  break;
566 
567  case 2212: // p
568  track->SetCharge(+1);
569  pid[AliAODTrack::kProton] = 1.;
570  track->SetPID(pid);
571  break;
572 
573  case -2212: // anti-p
574  track->SetCharge(-1);
575  pid[AliAODTrack::kProton] = 1.;
576  track->SetPID(pid);
577  break;
578 
579  case 310: // K0S
580  track->SetCharge(0);
581  pid[AliAODTrack::kUnknown] = 1.;
582  track->SetPID(pid);
583  break;
584 
585  case 311: // K0
586  track->SetCharge(0);
587  pid[AliAODTrack::kUnknown] = 1.;
588  track->SetPID(pid);
589  break;
590 
591  case -311: // anti-K0
592  track->SetCharge(0);
593  pid[AliAODTrack::kUnknown] = 1.;
594  track->SetPID(pid);
595  break;
596 
597  case 221: // eta
598  track->SetCharge(0);
599  pid[AliAODTrack::kUnknown] = 1.;
600  track->SetPID(pid);
601  break;
602 
603  case 3122: // lambda
604  track->SetCharge(0);
605  pid[AliAODTrack::kUnknown] = 1.;
606  track->SetPID(pid);
607  break;
608 
609  case 3222: // Sigma+
610  track->SetCharge(+1);
611  pid[AliAODTrack::kUnknown] = 1.;
612  track->SetPID(pid);
613  break;
614 
615  case 3212: // Sigma0
616  track->SetCharge(-1);
617  pid[AliAODTrack::kUnknown] = 1.;
618  track->SetPID(pid);
619  break;
620 
621  case 3112: // Sigma-
622  track->SetCharge(-1);
623  pid[AliAODTrack::kUnknown] = 1.;
624  track->SetPID(pid);
625  break;
626 
627  case 3322: // Xi0
628  track->SetCharge(0);
629  pid[AliAODTrack::kUnknown] = 1.;
630  track->SetPID(pid);
631  break;
632 
633  case 3312: // Xi-
634  track->SetCharge(-1);
635  pid[AliAODTrack::kUnknown] = 1.;
636  track->SetPID(pid);
637  break;
638 
639  case 3334: // Omega-
640  track->SetCharge(-1);
641  pid[AliAODTrack::kUnknown] = 1.;
642  track->SetPID(pid);
643  break;
644 
645  case -2112: // n-bar
646  track->SetCharge(0);
647  pid[AliAODTrack::kUnknown] = 1.;
648  track->SetPID(pid);
649  break;
650 
651  case -3122: // anti-Lambda
652  track->SetCharge(0);
653  pid[AliAODTrack::kUnknown] = 1.;
654  track->SetPID(pid);
655  break;
656 
657  case -3222: // anti-Sigma-
658  track->SetCharge(-1);
659  pid[AliAODTrack::kUnknown] = 1.;
660  track->SetPID(pid);
661  break;
662 
663  case -3212: // anti-Sigma0
664  track->SetCharge(0);
665  pid[AliAODTrack::kUnknown] = 1.;
666  track->SetPID(pid);
667  break;
668 
669  case -3112: // anti-Sigma+
670  track->SetCharge(+1);
671  pid[AliAODTrack::kUnknown] = 1.;
672  track->SetPID(pid);
673  break;
674 
675  case -3322: // anti-Xi0
676  track->SetCharge(0);
677  pid[AliAODTrack::kUnknown] = 1.;
678  track->SetPID(pid);
679  break;
680 
681  case -3312: // anti-Xi+
682  track->SetCharge(+1);
683  break;
684 
685  case -3334: // anti-Omega+
686  track->SetCharge(+1);
687  pid[AliAODTrack::kUnknown] = 1.;
688  track->SetPID(pid);
689  break;
690 
691  case 411: // D+
692  track->SetCharge(+1);
693  pid[AliAODTrack::kUnknown] = 1.;
694  track->SetPID(pid);
695  break;
696 
697  case -411: // D-
698  track->SetCharge(-1);
699  pid[AliAODTrack::kUnknown] = 1.;
700  track->SetPID(pid);
701  break;
702 
703  case 421: // D0
704  track->SetCharge(0);
705  pid[AliAODTrack::kUnknown] = 1.;
706  track->SetPID(pid);
707  break;
708 
709  case -421: // anti-D0
710  track->SetCharge(0);
711  pid[AliAODTrack::kUnknown] = 1.;
712  track->SetPID(pid);
713  break;
714 
715  default : // unknown
716  track->SetCharge(-99);
717  pid[AliAODTrack::kUnknown] = 1.;
718  track->SetPID(pid);
719  }
720 
721  track->SetPID(pid);
722 
723  return;
724 }
725 
726 //____________________________________________________________________
728 //____________________________________________________________________
730  AliVCluster *calo) const
731 {
732  Float_t pid[13] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
733 
734  switch (pdgCode)
735  {
736  case 22: // gamma
737  pid[AliVCluster::kPhoton] = 1.;
738  calo->SetPID(pid);
739  break;
740 
741  case 11: // e-
742  pid[AliVCluster::kElectron] = 1.;
743  calo->SetPID(pid);
744  break;
745 
746  case -11: // e+
747  pid[AliVCluster::kElectron] = 1.;
748  calo->SetPID(pid);
749  break;
750 
751  case 13: // mu-
752  pid[AliVCluster::kCharged] = 1.;
753  calo->SetPID(pid);
754  break;
755 
756  case -13: // mu+
757  pid[AliVCluster::kCharged] = 1.;
758  calo->SetPID(pid);
759  break;
760 
761  case 111: // pi0
762  pid[AliVCluster::kPi0] = 1.;
763  calo->SetPID(pid);
764  break;
765 
766  case 211: // pi+
767  pid[AliVCluster::kCharged] = 1.;
768  calo->SetPID(pid);
769  break;
770 
771  case -211: // pi-
772  pid[AliVCluster::kCharged] = 1.;
773  calo->SetPID(pid);
774  break;
775 
776  case 130: // K0L
777  pid[AliVCluster::kKaon0] = 1.;
778  pid[AliVCluster::kNeutral] = 1;
779  calo->SetPID(pid);
780  break;
781 
782  case 321: // K+
783  pid[AliVCluster::kCharged] = 1.;
784  calo->SetPID(pid);
785  break;
786 
787  case -321: // K-
788  pid[AliVCluster::kCharged] = 1.;
789  calo->SetPID(pid);
790  break;
791 
792  case 2112: // n
793  pid[AliVCluster::kNeutron] = 1.;
794  pid[AliVCluster::kNeutral] = 1.;
795  calo->SetPID(pid);
796  break;
797 
798  case 2212: // p
799  pid[AliVCluster::kCharged] = 1.;
800  calo->SetPID(pid);
801  break;
802 
803  case -2212: // anti-p
804  pid[AliVCluster::kCharged] = 1.;
805  calo->SetPID(pid);
806  break;
807 
808  case 310: // K0S
809  pid[AliVCluster::kKaon0] = 1.;
810  pid[AliVCluster::kNeutral] = 1.;
811  calo->SetPID(pid);
812  break;
813 
814  case 311: // K0
815  pid[AliVCluster::kKaon0] = 1.;
816  pid[AliVCluster::kNeutral] = 1.;
817  calo->SetPID(pid);
818  break;
819 
820  case -311: // anti-K0
821  pid[AliVCluster::kKaon0] = 1.;
822  pid[AliVCluster::kNeutral] = 1.;
823  calo->SetPID(pid);
824  break;
825 
826  case 221: // eta
827  pid[AliVCluster::kNeutral] = 1.;
828  calo->SetPID(pid);
829  break;
830 
831  case 3122: // lambda
832  pid[AliVCluster::kUnknown] = 1.;
833  calo->SetPID(pid);
834  break;
835 
836  case 3222: // Sigma+
837  pid[AliVCluster::kUnknown] = 1.;
838  calo->SetPID(pid);
839  break;
840 
841  case 3212: // Sigma0
842  pid[AliVCluster::kUnknown] = 1.;
843  calo->SetPID(pid);
844  break;
845 
846  case 3112: // Sigma-
847  pid[AliVCluster::kUnknown] = 1.;
848  calo->SetPID(pid);
849  break;
850 
851  case 3322: // Xi0
852  pid[AliVCluster::kUnknown] = 1.;
853  calo->SetPID(pid);
854  break;
855 
856  case 3312: // Xi-
857  pid[AliVCluster::kUnknown] = 1.;
858  calo->SetPID(pid);
859  break;
860 
861  case 3334: // Omega-
862  pid[AliVCluster::kUnknown] = 1.;
863  calo->SetPID(pid);
864  break;
865 
866  case -2112: // n-bar
867  pid[AliVCluster::kNeutron] = 1.;
868  pid[AliVCluster::kNeutral] = 1.;
869  calo->SetPID(pid);
870  break;
871 
872  case -3122: // anti-Lambda
873  pid[AliVCluster::kUnknown] = 1.;
874  calo->SetPID(pid);
875  break;
876 
877  case -3222: // anti-Sigma-
878  pid[AliVCluster::kUnknown] = 1.;
879  calo->SetPID(pid);
880  break;
881 
882  case -3212: // anti-Sigma0
883  pid[AliVCluster::kUnknown] = 1.;
884  calo->SetPID(pid);
885  break;
886 
887  case -3112: // anti-Sigma+
888  pid[AliVCluster::kUnknown] = 1.;
889  calo->SetPID(pid);
890  break;
891 
892  case -3322: // anti-Xi0
893  pid[AliVCluster::kUnknown] = 1.;
894  calo->SetPID(pid);
895  break;
896 
897  case -3312: // anti-Xi+
898  pid[AliVCluster::kUnknown] = 1.;
899  calo->SetPID(pid);
900  break;
901 
902  case -3334: // anti-Omega+
903  pid[AliVCluster::kUnknown] = 1.;
904  calo->SetPID(pid);
905  break;
906 
907  case 411: // D+
908  pid[AliVCluster::kUnknown] = 1.;
909  calo->SetPID(pid);
910  break;
911 
912  case -411: // D-
913  pid[AliVCluster::kUnknown] = 1.;
914  calo->SetPID(pid);
915  break;
916 
917  case 421: // D0
918  pid[AliVCluster::kUnknown] = 1.;
919  calo->SetPID(pid);
920  break;
921 
922  case -421: // anti-D0
923  pid[AliVCluster::kUnknown] = 1.;
924  calo->SetPID(pid);
925  break;
926 
927  default : // unknown
928  pid[AliVCluster::kUnknown] = 1.;
929  calo->SetPID(pid);
930  }
931 
932  return;
933 }
934 
Int_t charge
Int_t pdg
virtual AliMCEvent * GetMC() const
TLorentzVector fGamDecayMom2
! Gamma decay 2 momentum
double Double_t
Definition: External.C:58
virtual ~AliCaloTrackMCReader()
Destructor.
TObjArray * fPHOSClusters
Temporal array with PHOS CaloClusters.
void InitParameters()
Initialize the parameters of the analysis.
void SetCaloClusterPID(Int_t pdgCode, AliVCluster *calo) const
Give a PID weight for CaloClusters equal to 1 depending on the particle type.
virtual void SetMC(AliMCEvent *const mc)
TLorentzVector fGamDecayMom1
! Gamma decay 1 momentum
Bool_t fOnlyGeneratorParticles
Use particles only generated by PYTHIA/HERWIG/... and not by the MC tranport G3/G4/FLUKA ...
void CheckOverlap(Float_t anglethres, Int_t imom, Int_t &iPrimary, Int_t &index, Int_t &pdg)
Check overlap of decay photons.
char Char_t
Definition: External.C:18
void SetInputOutputMCEvent(AliVEvent *esd, AliAODEvent *aod, AliMCEvent *mc)
Connect the input data pointer.
Float_t fEMCALPtMin
pT Threshold on emcal clusters.
Bool_t SkipNeutralParticles(Int_t pdg) const
Bool_t fDecayPi0
If not decayed, decay pi0 by hand.
Bool_t fFillCTS
Use data from CTS.
virtual AliGenEventHeader * GetGenEventHeader() const
TRandom * gRandom
Int_t fTrackMult[10]
Track multiplicity, count for different pT cuts.
AliCaloTrackMCReader()
Default constructor. Initialize parameters.
Float_t fPHOSPtMin
pT Threshold on phos clusters.
TLorentzVector fMomentum
! Momentum
Int_t fEventNumber
Event number.
Float_t fTrackMultPtCut[10]
Track multiplicity and sum pt cuts list.
TLorentzVector fPi0Momentum
! Pi0 momentum
Float_t fCTSPtMax
pT Threshold on charged particles.
Float_t fPHOSOverlapAngle
Aperture angle of photons from decay that is not resolved by PHOS, in radians.
int Int_t
Definition: External.C:63
Int_t fIndex2ndPhoton
Check overlap of first decay photon already done, internal use.
Bool_t FillInputEvent(Int_t iEntry, const char *currentFileName)
Fill the event counter and input lists that are needed, called by AliAnaCaloTrackCorrMaker.
TObjArray * fEMCALClusters
Temporal array with EMCAL CaloClusters.
float Float_t
Definition: External.C:68
virtual void SetOutputEvent(AliAODEvent *aod)
Class for filtering generated MC particles and prepare them as input for the analysis.
virtual Bool_t ComparePtHardAndJetPt()
Bool_t IsInFiducialCut(Float_t eta, Float_t phi, Int_t det) const
Float_t fTrackMultEtaCut
Track multiplicity eta cut.
Base class for event, clusters and tracks filtering and preparation for the analysis.
Bool_t KeepParticleWithStatus(Int_t status) const
void SetTrackChargeAndPID(Int_t pdgCode, AliAODTrack *track) const
Give a PID weight for tracks equal to 1 depending on the particle type.
TArrayI * fNeutralParticlesArray
Do not keep neutral particles of this list in calorimeter.
void Print(const Option_t *opt) const
Print some relevant parameters set for the analysis.
Float_t fEMCALOverlapAngle
Aperture angle of photons from decay that is not resolved by EMCAL, in radians.
virtual void Print(const Option_t *opt) const
Print some relevant parameters set for the analysis.
Bool_t fCheckFidCut
Do analysis for clusters in defined region.
Bool_t fComparePtHardAndJetPt
In MonteCarlo, jet events, reject fake events with wrong jet energy.
AliFiducialCut * fFiducialCut
Acceptance cuts.
TArrayI * fChargedParticlesArray
Keep charged particles of this list in calorimeter.
Bool_t fCheckOverlap
Check of overlapped photons from pi0 enter the calorimeter.
Bool_t fFillEMCAL
Use data from EMCAL.
Bool_t fFillPHOS
Use data from PHOS.
Float_t fCTSPtMin
pT Threshold on charged particles.
void GetVertex(Double_t v[3]) const
Bool_t fKeepAllStatus
Do or do not select particles depending on their status code.
void FillCalorimeters(Int_t &iParticle, Int_t motherIndex, Int_t pdg)
Fill CaloClusters or AliVParticles lists of PHOS or EMCAL.
TArrayI * fStatusArray
Keep particles with status of the list.
virtual void FillVertexArray()
TObjArray * fCTSTracks
Temporal array with tracks.
const char Option_t
Definition: External.C:48
Float_t fTrackSumPt[10]
Track sum pT, count for different pT cuts.
Bool_t KeepChargedParticles(Int_t pdg) const
bool Bool_t
Definition: External.C:53
Int_t fTrackMultNPtCut
Track multiplicty, number of pt cuts.
TString meson
Int_t fDataType
Select MC: Kinematics, Data: ESD/AOD, MCData: Both.