AliPhysics  aaf9c62 (aaf9c62)
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  // Main header
220 
221  if ( fGenEventHeader )
222  {
223  AliInfo(Form("Selected event header class <%s>, name <%s>; cocktail %p",
224  fGenEventHeader->ClassName(),
225  fGenEventHeader->GetName(),
226  GetMC()->GetCocktailList()));
227  }
228 
229  //----------------------------------------------------------------
230  // Reject the event if the event header name is not
231  // the one requested among the possible generators.
232  // Needed in case of cocktail MC generation with multiple options.
233  //----------------------------------------------------------------
235  {
236  if ( !fGenEventHeader ) return kFALSE;
237 
238  AliDebug(1,"Pass Event header selection");
239 
240  fhNEventsAfterCut->Fill(15.5);
241  }
242 
243  // Pythia header
244  TString pyGenName = "";
245  TString pyProcessName = "";
246  Int_t pyProcess = 0;
247  Int_t pyFirstGenPart = 0;
248  Int_t pythiaVersion = 0;
249 
252  pyGenName,pyProcessName,pyProcess,pyFirstGenPart,pythiaVersion);
254  {
255  AliInfo(Form("Pythia v%d name <%s>, process %d <%s>, first generated particle %d",
256  pythiaVersion, pyGenName.Data(), pyProcess, pyProcessName.Data(), pyFirstGenPart));
257  }
258 
259  //---------------------------------------------------------------------------
260  // In case of analysis of events with jets, skip those with jet pt > 5 pt hard
261  // To be used on for MC data in pT hard bins
262  //---------------------------------------------------------------------------
263 
265  {
266  if ( !ComparePtHardAndJetPt(pyProcess, pyProcessName) ) return kFALSE ;
267 
268  AliDebug(1,"Pass Pt Hard - Jet rejection");
269  }
270 
272  {
273  if ( !ComparePtHardAndClusterPt(pyProcess, pyProcessName) ) return kFALSE ;
274 
275  AliDebug(1,"Pass Pt Hard - Cluster rejection");
276  }
277 
278  // Fill Vertex array
279  FillVertexArray();
280 
281  Int_t iParticle = 0 ;
282  Int_t nparticles = GetMC()->GetNumberOfTracks() ;
283 
284  if(fOnlyGeneratorParticles) nparticles=GetMC()->GetNumberOfPrimaries();
285 
286  for (iParticle = 0 ; iParticle < nparticles ; iParticle++)
287  {
288  AliVParticle * particle = GetMC()->GetTrack(iParticle);
289  Float_t p[3];
290  Float_t x[3];
291 
292  //Keep particles with a given status
293  if(KeepParticleWithStatus(particle->MCStatusCode()) && (particle->Pt() > 0) )
294  {
295  //Skip bizarre particles, they crash when charge is calculated
296  // if(TMath::Abs(pdg) == 3124 || TMath::Abs(pdg) > 10000000) continue ;
297 
298  //charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
299  Int_t charge = particle->Charge();
300  Int_t momIndex = particle->GetMother();
301  Int_t pdg = particle->PdgCode();
302 
303  particle->Momentum(fMomentum);
304 
305  Float_t en = particle->E();
306  Float_t pt = particle->Pt();
307  Float_t eta = particle->Eta();
308  Float_t phi = particle->Phi();
309 
310  //---------- Charged particles ----------------------
311  if(charge != 0)
312  {
313 // if(TMath::Abs(eta)< fTrackMultEtaCut) fTrackMult++;
314  if(TMath::Abs(eta)< fTrackMultEtaCut)
315  {
316  for(Int_t iptCut = 0; iptCut < fTrackMultNPtCut; iptCut++ )
317  {
318  if(pt > fTrackMultPtCut[iptCut])
319  {
320  fTrackMult [iptCut]++;
321  fTrackSumPt[iptCut]+=pt;
322  }
323  }
324  }
325 
326 // if(fFillCTS && (pt > fCTSPtMin))
327  if(fFillCTS && (fCTSPtMin > pt || fCTSPtMax < pt)) continue ;
328  {
329  // Particles in CTS acceptance
330 
331  if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(eta,phi,kCTS)) continue;
332 
333  if(TMath::Abs(pdg) == 11 && (GetMC()->GetTrack(particle->GetMother()))->PdgCode()==22) continue ;
334 
335  AliDebug(2,Form("CTS : Selected tracks E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f",
336  en,pt,phi*TMath::RadToDeg(),eta));
337 
338  x[0] = particle->Xv(); x[1] = particle->Yv(); x[2] = particle->Zv();
339  p[0] = particle->Px(); p[1] = particle->Py(); p[2] = particle->Pz();
340  //Create object and write it to file
341  AliAODTrack *aodTrack = new AliAODTrack(0, iParticle, p, kTRUE, x, kFALSE,NULL, 0, 0,
342  // NULL,
343  0x0,//primary,
344  kFALSE, // No fit performed
345  kFALSE, // No fit performed
347  0);
348  SetTrackChargeAndPID(pdg, aodTrack);
349  fCTSTracks->Add(aodTrack);//reference the selected object to the list
350  }
351 
352  // Keep some charged particles in calorimeters lists
353  if((fFillPHOS || fFillEMCAL) && KeepChargedParticles(pdg))
354  FillCalorimeters(iParticle, momIndex, pdg);
355 
356  }//Charged
357 
358  //-------------Neutral particles ----------------------
359  else if(charge == 0 && (fFillPHOS || fFillEMCAL))
360  {
361  //Skip neutrinos or other neutral particles
362  //if(SkipNeutralParticles(pdg) || particle->IsPrimary()) continue ; // protection added (MG)
363  if(SkipNeutralParticles(pdg)) continue ;
364  //Fill particle/calocluster arrays
365  if(!fDecayPi0)
366  {
367  AliDebug(2,Form("Calo : Selected tracks E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f",
368  en,pt,phi*TMath::RadToDeg(),eta));
369  FillCalorimeters(iParticle, momIndex, pdg);
370  }
371  else
372  {
373  //Sometimes pi0 are stable for the generator, if needed decay it by hand
374  if(pdg == 111 )
375  {
376  if(pt > fPHOSPtMin || pt > fEMCALPtMin)
377  {
378  //Decay
379  //Double_t angle = 0;
380  particle->Momentum(fPi0Momentum);
381 
382  MakePi0Decay();//,angle);
383 
384  // Check if Pi0 is in the acceptance of the calorimeters, if aperture angle is small, keep it
385  // and fill particle/calocluster arrays
386  fMomentum.SetPxPyPzE(fGamDecayMom1.Px(),fGamDecayMom1.Py(), fGamDecayMom1.Pz(),fGamDecayMom1.E());
387  FillCalorimeters(iParticle,iParticle,22);
388 
389  fMomentum.SetPxPyPzE(fGamDecayMom2.Px(),fGamDecayMom2.Py(), fGamDecayMom2.Pz(),fGamDecayMom2.E());
390  FillCalorimeters(iParticle,iParticle,22);
391  }//pt cut
392  }//pi0
393  else
394  {
395  AliDebug(2,Form("Calo : Selected tracks E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f",
396  en,pt,phi*TMath::RadToDeg(),eta));
397  FillCalorimeters(iParticle,momIndex,pdg); //Add the rest
398  }
399  }
400  }//neutral particles
401  } //particle with correct status
402  }//particle loop
403 
404  fIndex2ndPhoton = -1; //In case of overlapping studies, reset for each event
405 
406  return kTRUE;
407 }
408 
409 //_____________________________________________________________________
412 //_____________________________________________________________________
414 {
415  if(!fKeepAllStatus)
416  {
417  for(Int_t i= 0; i < fStatusArray->GetSize(); i++)
418  if(status == fStatusArray->At(i)) return kTRUE ;
419 
420  return kFALSE;
421  }
422  else
423  return kTRUE ;
424 }
425 
426 //________________________________________________________________
429 //________________________________________________________________
431 {
432  for(Int_t i= 0; i < fChargedParticlesArray->GetSize(); i++)
433  {
434  if(TMath::Abs(pdg) == fChargedParticlesArray->At(i)) return kTRUE ;
435  }
436 
437  return kFALSE;
438 }
439 
440 //__________________________________________________________
442 //__________________________________________________________
444 {
445  if(! opt)
446  return;
447 
449 
450  printf("**** Print **** %s %s ****\n", GetName(), GetTitle() ) ;
451 
452  printf("Decay Pi0? : %d\n", fDecayPi0) ;
453  printf("Check Overlap in Calo? : %d\n", fCheckOverlap) ;
454  printf("Keep all status? : %d\n", fKeepAllStatus) ;
455 
456  if(!fKeepAllStatus) printf("Keep particles with status : ");
457  for(Int_t i= 0; i < fStatusArray->GetSize(); i++)
458  printf(" %d ; ", fStatusArray->At(i));
459  printf("\n");
460 
461  printf("Skip neutral particles in calo : ");
462  for(Int_t i= 0; i < fNeutralParticlesArray->GetSize(); i++)
463  printf(" %d ; ", fNeutralParticlesArray->At(i));
464  printf("\n");
465 
466  printf("Keep charged particles in calo : ");
467  for(Int_t i= 0; i < fChargedParticlesArray->GetSize(); i++)
468  printf(" %d ; ", fChargedParticlesArray->At(i));
469  printf("\n");
470 }
471 
472 //_______________________________________
476 //_______________________________________
477 void AliCaloTrackMCReader::MakePi0Decay() //, Double_t &angle)
478 {
479  // cout<<"Boost vector"<<endl;
480  Double_t mPi0 = TDatabasePDG::Instance()->GetParticle(111)->Mass();
481  TVector3 b = fPi0Momentum.BoostVector();
482  //cout<<"Parameters"<<endl;
483  //Double_t mPi0 = fPi0Momentum.M();
484  Double_t phi = TMath::TwoPi() * gRandom->Rndm();
485  Double_t cosThe = 2 * gRandom->Rndm() - 1;
486  Double_t cosPhi = TMath::Cos(phi);
487  Double_t sinPhi = TMath::Sin(phi);
488  Double_t sinThe = TMath::Sqrt(1-cosThe*cosThe);
489  Double_t ePi0 = mPi0/2.;
490  //cout<<"ePi0 "<<ePi0<<endl;
491  //cout<<"Components"<<endl;
492  fGamDecayMom1.SetPx(+ePi0*cosPhi*sinThe);
493  fGamDecayMom1.SetPy(+ePi0*sinPhi*sinThe);
494  fGamDecayMom1.SetPz(+ePi0*cosThe);
495  fGamDecayMom1.SetE(ePi0);
496  //cout<<"fGamDecayMom1: "<<fGamDecayMom1.Px()<<" "<<fGamDecayMom1.Py()<<" "<<fGamDecayMom1.Pz()<<" "<<fGamDecayMom1.E()<<endl;
497  //cout<<"fGamDecayMom1 Mass: "<<fGamDecayMom1.Px()*fGamDecayMom1.Px()+fGamDecayMom1.Py()*fGamDecayMom1.Py()+fGamDecayMom1.Pz()*fGamDecayMom1.Pz()-fGamDecayMom1.E()*fGamDecayMom1.E()<<endl;
498  fGamDecayMom2.SetPx(-ePi0*cosPhi*sinThe);
499  fGamDecayMom2.SetPy(-ePi0*sinPhi*sinThe);
500  fGamDecayMom2.SetPz(-ePi0*cosThe);
501  fGamDecayMom2.SetE(ePi0);
502  //cout<<"fGamDecayMom2: "<<fGamDecayMom2.Px()<<" "<<fGamDecayMom2.Py()<<" "<<fGamDecayMom2.Pz()<<" "<<fGamDecayMom2.E()<<endl;
503  //cout<<"fGamDecayMom2 Mass: "<<fGamDecayMom2.Px()*fGamDecayMom2.Px()+fGamDecayMom2.Py()*fGamDecayMom2.Py()+fGamDecayMom2.Pz()*fGamDecayMom2.Pz()-fGamDecayMom2.E()*fGamDecayMom2.E()<<endl;
504  //cout<<"Boost "<<b.X()<<" "<<b.Y()<<" "<<b.Z()<<endl;
505  fGamDecayMom1.Boost(b);
506  //cout<<"fGamDecayMom1: "<<fGamDecayMom1.Px()<<" "<<fGamDecayMom1.Py()<<" "<<fGamDecayMom1.Pz()<<" "<<fGamDecayMom1.E()<<endl;
507  fGamDecayMom2.Boost(b);
508  //cout<<"fGamDecayMom2: "<<fGamDecayMom2.Px()<<" "<<fGamDecayMom2.Py()<<" "<<fGamDecayMom2.Pz()<<" "<<fGamDecayMom2.E()<<endl;
509  //cout<<"angle"<<endl;
510  //angle = fGamDecayMom1.Angle(fGamDecayMom2.Vect());
511  //cout<<angle<<endl;
512 }
513 
514 //__________________________________________________________________
516 //__________________________________________________________________
518  AliAODEvent* aod,
519  AliMCEvent* mc)
520 {
521  SetMC(mc);
522  SetOutputEvent(aod);
523 }
524 
525 //________________________________________________________________
528 //________________________________________________________________
530 {
531  for(Int_t i= 0; i < fNeutralParticlesArray->GetSize(); i++)
532  {
533  if(TMath::Abs(pdg) == fNeutralParticlesArray->At(i)) return kTRUE ;
534  }
535 
536  return kFALSE;
537 }
538 
539 //_______________________________________________________________________
541 //_______________________________________________________________________
543  AliAODTrack *track) const
544 {
545  Float_t pid[10] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
546 
547  switch (pdgCode)
548  {
549  case 22: // gamma
550  track->SetCharge(0);
551  pid[AliAODTrack::kUnknown] = 1.;
552  track->SetPID(pid);
553  break;
554 
555  case 11: // e-
556  track->SetCharge(-1);
557  pid[AliAODTrack::kElectron] = 1.;
558  track->SetPID(pid);
559  break;
560 
561  case -11: // e+
562  track->SetCharge(+1);
563  pid[AliAODTrack::kElectron] = 1.;
564  track->SetPID(pid);
565  break;
566 
567  case 13: // mu-
568  track->SetCharge(-1);
569  pid[AliAODTrack::kMuon] = 1.;
570  track->SetPID(pid);
571  break;
572 
573  case -13: // mu+
574  track->SetCharge(+1);
575  pid[AliAODTrack::kMuon] = 1.;
576  track->SetPID(pid);
577  break;
578 
579  case 111: // pi0
580  track->SetCharge(0);
581  pid[AliAODTrack::kUnknown] = 1.;
582  track->SetPID(pid);
583  break;
584 
585  case 211: // pi+
586  track->SetCharge(+1);
587  pid[AliAODTrack::kPion] = 1.;
588  track->SetPID(pid);
589  break;
590 
591  case -211: // pi-
592  track->SetCharge(-1);
593  pid[AliAODTrack::kPion] = 1.;
594  track->SetPID(pid);
595  break;
596 
597  case 130: // K0L
598  track->SetCharge(0);
599  pid[AliAODTrack::kUnknown] = 1.;
600  track->SetPID(pid);
601  break;
602 
603  case 321: // K+
604  track->SetCharge(+1);
605  pid[AliAODTrack::kKaon] = 1.;
606  track->SetPID(pid);
607  break;
608 
609  case -321: // K-
610  track->SetCharge(-1);
611  pid[AliAODTrack::kKaon] = 1.;
612  track->SetPID(pid);
613  break;
614 
615  case 2112: // n
616  track->SetCharge(0);
617  pid[AliAODTrack::kUnknown] = 1.;
618  track->SetPID(pid);
619  break;
620 
621  case 2212: // p
622  track->SetCharge(+1);
623  pid[AliAODTrack::kProton] = 1.;
624  track->SetPID(pid);
625  break;
626 
627  case -2212: // anti-p
628  track->SetCharge(-1);
629  pid[AliAODTrack::kProton] = 1.;
630  track->SetPID(pid);
631  break;
632 
633  case 310: // K0S
634  track->SetCharge(0);
635  pid[AliAODTrack::kUnknown] = 1.;
636  track->SetPID(pid);
637  break;
638 
639  case 311: // K0
640  track->SetCharge(0);
641  pid[AliAODTrack::kUnknown] = 1.;
642  track->SetPID(pid);
643  break;
644 
645  case -311: // anti-K0
646  track->SetCharge(0);
647  pid[AliAODTrack::kUnknown] = 1.;
648  track->SetPID(pid);
649  break;
650 
651  case 221: // eta
652  track->SetCharge(0);
653  pid[AliAODTrack::kUnknown] = 1.;
654  track->SetPID(pid);
655  break;
656 
657  case 3122: // lambda
658  track->SetCharge(0);
659  pid[AliAODTrack::kUnknown] = 1.;
660  track->SetPID(pid);
661  break;
662 
663  case 3222: // Sigma+
664  track->SetCharge(+1);
665  pid[AliAODTrack::kUnknown] = 1.;
666  track->SetPID(pid);
667  break;
668 
669  case 3212: // Sigma0
670  track->SetCharge(-1);
671  pid[AliAODTrack::kUnknown] = 1.;
672  track->SetPID(pid);
673  break;
674 
675  case 3112: // Sigma-
676  track->SetCharge(-1);
677  pid[AliAODTrack::kUnknown] = 1.;
678  track->SetPID(pid);
679  break;
680 
681  case 3322: // Xi0
682  track->SetCharge(0);
683  pid[AliAODTrack::kUnknown] = 1.;
684  track->SetPID(pid);
685  break;
686 
687  case 3312: // Xi-
688  track->SetCharge(-1);
689  pid[AliAODTrack::kUnknown] = 1.;
690  track->SetPID(pid);
691  break;
692 
693  case 3334: // Omega-
694  track->SetCharge(-1);
695  pid[AliAODTrack::kUnknown] = 1.;
696  track->SetPID(pid);
697  break;
698 
699  case -2112: // n-bar
700  track->SetCharge(0);
701  pid[AliAODTrack::kUnknown] = 1.;
702  track->SetPID(pid);
703  break;
704 
705  case -3122: // anti-Lambda
706  track->SetCharge(0);
707  pid[AliAODTrack::kUnknown] = 1.;
708  track->SetPID(pid);
709  break;
710 
711  case -3222: // anti-Sigma-
712  track->SetCharge(-1);
713  pid[AliAODTrack::kUnknown] = 1.;
714  track->SetPID(pid);
715  break;
716 
717  case -3212: // anti-Sigma0
718  track->SetCharge(0);
719  pid[AliAODTrack::kUnknown] = 1.;
720  track->SetPID(pid);
721  break;
722 
723  case -3112: // anti-Sigma+
724  track->SetCharge(+1);
725  pid[AliAODTrack::kUnknown] = 1.;
726  track->SetPID(pid);
727  break;
728 
729  case -3322: // anti-Xi0
730  track->SetCharge(0);
731  pid[AliAODTrack::kUnknown] = 1.;
732  track->SetPID(pid);
733  break;
734 
735  case -3312: // anti-Xi+
736  track->SetCharge(+1);
737  break;
738 
739  case -3334: // anti-Omega+
740  track->SetCharge(+1);
741  pid[AliAODTrack::kUnknown] = 1.;
742  track->SetPID(pid);
743  break;
744 
745  case 411: // D+
746  track->SetCharge(+1);
747  pid[AliAODTrack::kUnknown] = 1.;
748  track->SetPID(pid);
749  break;
750 
751  case -411: // D-
752  track->SetCharge(-1);
753  pid[AliAODTrack::kUnknown] = 1.;
754  track->SetPID(pid);
755  break;
756 
757  case 421: // D0
758  track->SetCharge(0);
759  pid[AliAODTrack::kUnknown] = 1.;
760  track->SetPID(pid);
761  break;
762 
763  case -421: // anti-D0
764  track->SetCharge(0);
765  pid[AliAODTrack::kUnknown] = 1.;
766  track->SetPID(pid);
767  break;
768 
769  default : // unknown
770  track->SetCharge(-99);
771  pid[AliAODTrack::kUnknown] = 1.;
772  track->SetPID(pid);
773  }
774 
775  track->SetPID(pid);
776 
777  return;
778 }
779 
780 //____________________________________________________________________
782 //____________________________________________________________________
784  AliVCluster *calo) const
785 {
786  Float_t pid[13] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
787 
788  switch (pdgCode)
789  {
790  case 22: // gamma
791  pid[AliVCluster::kPhoton] = 1.;
792  calo->SetPID(pid);
793  break;
794 
795  case 11: // e-
796  pid[AliVCluster::kElectron] = 1.;
797  calo->SetPID(pid);
798  break;
799 
800  case -11: // e+
801  pid[AliVCluster::kElectron] = 1.;
802  calo->SetPID(pid);
803  break;
804 
805  case 13: // mu-
806  pid[AliVCluster::kCharged] = 1.;
807  calo->SetPID(pid);
808  break;
809 
810  case -13: // mu+
811  pid[AliVCluster::kCharged] = 1.;
812  calo->SetPID(pid);
813  break;
814 
815  case 111: // pi0
816  pid[AliVCluster::kPi0] = 1.;
817  calo->SetPID(pid);
818  break;
819 
820  case 211: // pi+
821  pid[AliVCluster::kCharged] = 1.;
822  calo->SetPID(pid);
823  break;
824 
825  case -211: // pi-
826  pid[AliVCluster::kCharged] = 1.;
827  calo->SetPID(pid);
828  break;
829 
830  case 130: // K0L
831  pid[AliVCluster::kKaon0] = 1.;
832  pid[AliVCluster::kNeutral] = 1;
833  calo->SetPID(pid);
834  break;
835 
836  case 321: // K+
837  pid[AliVCluster::kCharged] = 1.;
838  calo->SetPID(pid);
839  break;
840 
841  case -321: // K-
842  pid[AliVCluster::kCharged] = 1.;
843  calo->SetPID(pid);
844  break;
845 
846  case 2112: // n
847  pid[AliVCluster::kNeutron] = 1.;
848  pid[AliVCluster::kNeutral] = 1.;
849  calo->SetPID(pid);
850  break;
851 
852  case 2212: // p
853  pid[AliVCluster::kCharged] = 1.;
854  calo->SetPID(pid);
855  break;
856 
857  case -2212: // anti-p
858  pid[AliVCluster::kCharged] = 1.;
859  calo->SetPID(pid);
860  break;
861 
862  case 310: // K0S
863  pid[AliVCluster::kKaon0] = 1.;
864  pid[AliVCluster::kNeutral] = 1.;
865  calo->SetPID(pid);
866  break;
867 
868  case 311: // K0
869  pid[AliVCluster::kKaon0] = 1.;
870  pid[AliVCluster::kNeutral] = 1.;
871  calo->SetPID(pid);
872  break;
873 
874  case -311: // anti-K0
875  pid[AliVCluster::kKaon0] = 1.;
876  pid[AliVCluster::kNeutral] = 1.;
877  calo->SetPID(pid);
878  break;
879 
880  case 221: // eta
881  pid[AliVCluster::kNeutral] = 1.;
882  calo->SetPID(pid);
883  break;
884 
885  case 3122: // lambda
886  pid[AliVCluster::kUnknown] = 1.;
887  calo->SetPID(pid);
888  break;
889 
890  case 3222: // Sigma+
891  pid[AliVCluster::kUnknown] = 1.;
892  calo->SetPID(pid);
893  break;
894 
895  case 3212: // Sigma0
896  pid[AliVCluster::kUnknown] = 1.;
897  calo->SetPID(pid);
898  break;
899 
900  case 3112: // Sigma-
901  pid[AliVCluster::kUnknown] = 1.;
902  calo->SetPID(pid);
903  break;
904 
905  case 3322: // Xi0
906  pid[AliVCluster::kUnknown] = 1.;
907  calo->SetPID(pid);
908  break;
909 
910  case 3312: // Xi-
911  pid[AliVCluster::kUnknown] = 1.;
912  calo->SetPID(pid);
913  break;
914 
915  case 3334: // Omega-
916  pid[AliVCluster::kUnknown] = 1.;
917  calo->SetPID(pid);
918  break;
919 
920  case -2112: // n-bar
921  pid[AliVCluster::kNeutron] = 1.;
922  pid[AliVCluster::kNeutral] = 1.;
923  calo->SetPID(pid);
924  break;
925 
926  case -3122: // anti-Lambda
927  pid[AliVCluster::kUnknown] = 1.;
928  calo->SetPID(pid);
929  break;
930 
931  case -3222: // anti-Sigma-
932  pid[AliVCluster::kUnknown] = 1.;
933  calo->SetPID(pid);
934  break;
935 
936  case -3212: // anti-Sigma0
937  pid[AliVCluster::kUnknown] = 1.;
938  calo->SetPID(pid);
939  break;
940 
941  case -3112: // anti-Sigma+
942  pid[AliVCluster::kUnknown] = 1.;
943  calo->SetPID(pid);
944  break;
945 
946  case -3322: // anti-Xi0
947  pid[AliVCluster::kUnknown] = 1.;
948  calo->SetPID(pid);
949  break;
950 
951  case -3312: // anti-Xi+
952  pid[AliVCluster::kUnknown] = 1.;
953  calo->SetPID(pid);
954  break;
955 
956  case -3334: // anti-Omega+
957  pid[AliVCluster::kUnknown] = 1.;
958  calo->SetPID(pid);
959  break;
960 
961  case 411: // D+
962  pid[AliVCluster::kUnknown] = 1.;
963  calo->SetPID(pid);
964  break;
965 
966  case -411: // D-
967  pid[AliVCluster::kUnknown] = 1.;
968  calo->SetPID(pid);
969  break;
970 
971  case 421: // D0
972  pid[AliVCluster::kUnknown] = 1.;
973  calo->SetPID(pid);
974  break;
975 
976  case -421: // anti-D0
977  pid[AliVCluster::kUnknown] = 1.;
978  calo->SetPID(pid);
979  break;
980 
981  default : // unknown
982  pid[AliVCluster::kUnknown] = 1.;
983  calo->SetPID(pid);
984  }
985 
986  return;
987 }
988 
Int_t charge
Int_t pdg
virtual AliMCEvent * GetMC() const
Bool_t fComparePtHardAndClusterPt
In MonteCarlo, jet events, reject events with too large cluster energy.
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.
AliGenEventHeader * GetGenEventHeader() const
static AliGenPythiaEventHeader * GetPythiaEventHeader(AliMCEvent *mcevent, TString selecHeaderName, TString &genName, TString &processName, Int_t &process, Int_t &firstParticle, Int_t &pythiaVersion)
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.
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.
AliGenPythiaEventHeader * fGenPythiaEventHeader
! Event header casted to pythia
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.
AliGenEventHeader * fGenEventHeader
! Event header
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.
Bool_t IsInFiducialCut(Float_t eta, Float_t phi, Int_t det) const
virtual Bool_t ComparePtHardAndClusterPt(Int_t process, TString processName)
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 Bool_t ComparePtHardAndJetPt(Int_t process, TString processName)
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.
TH1I * fhNEventsAfterCut
! Each bin represents number of events resulting after a given selection cut: vertex, trigger, ...
TString fMCGenerEventHeaderToAccept
Accept events that contain at least this event header name.
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.