1)AliAODParticleCorrelation.h: Data member fIsolation changed from float to bool
[u/mrichter/AliRoot.git] / PWG4 / PartCorrBase / AliCaloTrackMCReader.cxx
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 /* $Id:  $ */
16
17 //_________________________________________________________________________
18 // Class for reading data (Kinematics) in order to do prompt gamma 
19 // or other particle identification and correlations
20 //
21 //*-- Author: Gustavo Conesa (LNF-INFN) 
22 //////////////////////////////////////////////////////////////////////////////
23
24
25 // --- ROOT system ---
26
27 #include <TParticle.h>
28 #include <TRandom.h>
29 #include <TClonesArray.h>
30 //#include "Riostream.h"
31
32 //---- ANALYSIS system ----
33 #include "AliCaloTrackMCReader.h" 
34 #include "AliLog.h"
35 #include "AliGenEventHeader.h"
36 #include "AliStack.h"
37 #include "AliAODCaloCluster.h"
38 #include "AliAODTrack.h"
39
40 ClassImp(AliCaloTrackMCReader)
41
42 //____________________________________________________________________________
43 AliCaloTrackMCReader::AliCaloTrackMCReader() : 
44   AliCaloTrackReader(), fDecayPi0(0), 
45   fNeutralParticlesArray(0x0),    fChargedParticlesArray(0x0), 
46   fStatusArray(0x0), fKeepAllStatus(0), fClonesArrayType(0), 
47   fCheckOverlap(0),  fEMCALOverlapAngle(0),fPHOSOverlapAngle(0),
48   fIndex2ndPhoton(0)
49 {
50   //Ctor
51   
52   //Initialize parameters
53   InitParameters();
54   fDataType = kMC;  
55   
56 }
57
58 //____________________________________________________________________________
59 AliCaloTrackMCReader::AliCaloTrackMCReader(const AliCaloTrackMCReader & g) :   
60   AliCaloTrackReader(g), fDecayPi0(g.fDecayPi0), 
61   fNeutralParticlesArray(g.fNeutralParticlesArray?new TArrayI(*g.fNeutralParticlesArray):0x0),
62   fChargedParticlesArray(g.fChargedParticlesArray?new TArrayI(*g.fChargedParticlesArray):0x0),
63   fStatusArray(g.fStatusArray?new TArrayI(*g.fStatusArray):0x0),
64   fKeepAllStatus(g.fKeepAllStatus), fClonesArrayType(g.fClonesArrayType),
65   fCheckOverlap(g.fCheckOverlap),
66   fEMCALOverlapAngle( g.fEMCALOverlapAngle), fPHOSOverlapAngle(g.fPHOSOverlapAngle),
67   fIndex2ndPhoton(g.fIndex2ndPhoton)
68 {
69   // cpy ctor
70 }
71
72 //_________________________________________________________________________
73 AliCaloTrackMCReader & AliCaloTrackMCReader::operator = (const AliCaloTrackMCReader & source)
74 {
75   // assignment operator
76
77   if(&source == this) return *this;
78
79   fDecayPi0 = source.fDecayPi0; 
80
81   delete fChargedParticlesArray;
82   fChargedParticlesArray = source.fChargedParticlesArray?new TArrayI(*source.fChargedParticlesArray):0x0;
83
84   delete fNeutralParticlesArray;
85   fNeutralParticlesArray = source.fNeutralParticlesArray?new TArrayI(*source.fNeutralParticlesArray):0x0;
86
87   delete fStatusArray;
88   fStatusArray = source.fStatusArray?new TArrayI(*source.fStatusArray):0x0;
89  
90   fKeepAllStatus = source.fKeepAllStatus ;
91   fClonesArrayType = source.fClonesArrayType ;
92
93   return *this;
94
95 }
96
97 //_________________________________
98 AliCaloTrackMCReader::~AliCaloTrackMCReader() {
99   //Dtor
100
101   if(fChargedParticlesArray) delete fChargedParticlesArray ;
102   if(fNeutralParticlesArray) delete fNeutralParticlesArray ;
103   if(fStatusArray) delete fStatusArray ;
104
105 }
106
107 //____________________________________________________________________________
108 void AliCaloTrackMCReader::GetVertex(Double_t  v[3]) const {
109   //Return vertex position
110
111   TArrayF pv;
112   GetGenEventHeader()->PrimaryVertex(pv);
113   v[0]=pv.At(0);
114   v[1]=pv.At(1);
115   v[2]=pv.At(2);
116
117 }
118
119
120 //_______________________________________________________________
121 void AliCaloTrackMCReader::InitParameters()
122 {
123   
124   //Initialize the parameters of the analysis.
125
126   fDecayPi0 = kFALSE;
127
128   fChargedParticlesArray = new TArrayI(1);
129   fChargedParticlesArray->SetAt(11,0);  
130   //Int_t pdgarray[]={12,14,16};// skip neutrinos
131   //fNeutralParticlesArray = new TArrayI(3, pdgarray);
132   fNeutralParticlesArray = new TArrayI(3);
133   fNeutralParticlesArray->SetAt(12,0); fNeutralParticlesArray->SetAt(14,1); fNeutralParticlesArray->SetAt(16,2); 
134   fStatusArray = new TArrayI(1);
135   fStatusArray->SetAt(1,0); 
136  
137   fKeepAllStatus = kTRUE;
138   fClonesArrayType = kAliAOD ;
139
140   fCheckOverlap = kFALSE;
141   fEMCALOverlapAngle = 2.5 * TMath::DegToRad();
142   fPHOSOverlapAngle = 0.5 * TMath::DegToRad();
143   fIndex2ndPhoton = -1;
144 }
145 //____________________________________________________________________________
146 void  AliCaloTrackMCReader::CheckOverlap(const Float_t anglethres, const Int_t imom, Int_t & iPrimary, Int_t & index, TLorentzVector & mom, Int_t & pdg) {
147   //Check overlap of decay photons
148   if( fIndex2ndPhoton==iPrimary ){
149     fIndex2ndPhoton=-1;
150     return;
151   }
152   else fIndex2ndPhoton=-1;
153   
154
155   if(pdg!=22) return;
156   
157   TLorentzVector ph1, ph2;
158   TParticle *meson = GetStack()->Particle(imom);
159   Int_t mepdg = meson->GetPdgCode();
160   Int_t idaug1 = meson->GetFirstDaughter();
161   if((mepdg == 111 || mepdg == 221 ) && meson->GetNDaughters() == 2){ //Check only decay in 2 photons
162     TParticle * d1 = GetStack()->Particle(idaug1);
163     TParticle  *d2 = GetStack()->Particle(idaug1+1);
164     if(d1->GetPdgCode() == 22 && d2->GetPdgCode() == 22 ){
165       d1->Momentum(ph1);
166       d2->Momentum(ph2);
167       //printf("angle %2.2f\n",ph1.Angle(ph2.Vect()));
168       
169       if(anglethres >  ph1.Angle(ph2.Vect())){    
170         //Keep the meson
171         pdg=mepdg;
172         index=imom;
173         meson->Momentum(mom);
174         //printf("Overlap:: pt %2.2f, phi %2.2f, eta %2.2f\n",mom.Pt(),mom.Phi(),mom.Eta());
175         if(iPrimary == idaug1) iPrimary++; //skip next photon in list
176       }
177       else{
178         //Do not check overlapping for next decay photon from same meson
179         if(iPrimary == idaug1) {fIndex2ndPhoton = idaug1+1;
180         }
181
182       }
183     }
184   }//Meson Decay with 2 photon daughters
185 }
186
187 //____________________________________________________________________________
188 void  AliCaloTrackMCReader::FillCalorimeters(Int_t & iParticle, TParticle* particle, TLorentzVector momentum,
189                                              Int_t &indexPHOS, Int_t &indexEMCAL) {
190         //Fill AODCaloClusters or TParticles lists of PHOS or EMCAL
191         //In PHOS
192         if(fFillPHOS && fFidutialCut->IsInFidutialCut(momentum,"PHOS") && momentum.Pt() > fPHOSPtMin){
193                 
194                 if(fClonesArrayType == kTParticle) new((*fAODPHOS)[indexPHOS++])       TParticle(*particle) ;
195                 else{
196                         Int_t index = iParticle ;
197                         Int_t pdg = TMath::Abs(particle->GetPdgCode());
198                         if(fCheckOverlap) 
199                           CheckOverlap(fPHOSOverlapAngle,particle->GetFirstMother(),index, iParticle, momentum, pdg);
200
201                         Char_t ttype= AliAODCluster::kPHOSNeutral;      
202                         Int_t labels[] = {index};
203                         Float_t x[] = {momentum.X(), momentum.Y(), momentum.Z()};
204                         AliAODCaloCluster *calo = new((*fAODPHOS)[indexPHOS++]) 
205                           AliAODCaloCluster(index,1,labels,momentum.E(), x, NULL, ttype, 0);
206                 
207                         SetCaloClusterPID(pdg,calo) ;
208                         if(fDebug > 3 && momentum.Pt() > 0.2)
209                                         printf("Fill MC PHOS :: Selected tracks %s E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
210                                                         particle->GetName(),momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());                        
211                 }
212         }
213         //In EMCAL
214         else if(fFillEMCAL && fFidutialCut->IsInFidutialCut(momentum,"EMCAL") && momentum.Pt() > fEMCALPtMin){
215                 if(fClonesArrayType == kTParticle) new((*fAODEMCAL)[indexEMCAL++])       TParticle(*particle) ;
216                 else{
217                                         Int_t index = iParticle ;
218                         Int_t pdg = TMath::Abs(particle->GetPdgCode());
219                         //Int_t pdgorg=pdg;
220                         if(fCheckOverlap) 
221                           CheckOverlap(fEMCALOverlapAngle,particle->GetFirstMother(),iParticle, index, momentum, pdg);
222
223                         Char_t ttype= AliAODCluster::kEMCALClusterv1;
224                         Int_t labels[] = {index};
225                         Float_t x[] = {momentum.X(), momentum.Y(), momentum.Z()};
226                         AliAODCaloCluster *calo = new((*fAODEMCAL)[indexEMCAL++]) 
227                         AliAODCaloCluster(iParticle,1,labels,momentum.E(), x, NULL, ttype, 0);
228                      
229                         SetCaloClusterPID(pdg,calo) ;
230                         if(fDebug > 3 && momentum.Pt() > 0.2)
231                                         printf("Fill MC EMCAL :: Selected tracks %s E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
232                                                         particle->GetName(),momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());        
233                 }
234         }
235 }
236
237 //____________________________________________________________________________
238 void AliCaloTrackMCReader::FillInputEvent(Int_t iEntry){
239         //Fill the event counter and input lists that are needed, called by the analysis maker.
240         
241         fEventNumber = iEntry;
242         
243         if(fClonesArrayType == kTParticle){
244                 fAODCTS = new TClonesArray("TParticle",0);
245                 fAODEMCAL = new TClonesArray("TParticle",0);
246                 fAODPHOS = new TClonesArray("TParticle",0);
247         }
248         else if(fClonesArrayType == kAliAOD){
249                 fAODCTS = new TClonesArray("AliAODTrack",0);
250                 fAODPHOS = new TClonesArray("AliAODCaloCluster",0);
251                 fAODEMCAL = new TClonesArray("AliAODCaloCluster",0);
252         }
253         else {AliFatal("Wrong clones type");}
254         
255         
256         Int_t indexCh    = 0 ;
257         Int_t indexEMCAL = 0 ;
258         Int_t indexPHOS  = 0 ;
259         
260         Int_t iParticle = 0 ;
261         Double_t charge = 0.;
262
263         for (iParticle = 0 ; iParticle <  GetStack()->GetNtrack() ; iParticle++) {
264                 TParticle * particle = GetStack()->Particle(iParticle);
265                 TLorentzVector momentum;
266                 Float_t p[3];
267                 Float_t x[3];
268                 Int_t pdg = particle->GetPdgCode();                                             
269
270                 //Keep particles with a given status 
271                 if(KeepParticleWithStatus(particle->GetStatusCode()) && (particle->Pt() > 0) ){
272                 
273                     //Skip bizarre particles, they crash when charge is calculated
274                   //    if(TMath::Abs(pdg) == 3124 || TMath::Abs(pdg) > 10000000) continue ;
275                         
276                         charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
277                         particle->Momentum(momentum);
278                         //---------- Charged particles ----------------------
279                         if((charge != 0) && (momentum.Pt() > fCTSPtMin) && (fFidutialCut->IsInFidutialCut(momentum,"CTS"))){
280                                 if(fFillCTS){
281                                         //Particles in CTS acceptance
282                                         if(fDebug > 3 && momentum.Pt() > 0.2)
283                                                 printf("Fill MC CTS :: Selected tracks E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
284                                                                 momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
285                                         
286                                         if(fClonesArrayType == kTParticle) new((*fAODCTS)[indexCh++])       TParticle(*particle) ;
287                                         else{
288                                                 x[0] = particle->Vx(); x[1] = particle->Vy(); x[2] = particle->Vz();
289                                                 p[0] = particle->Px(); p[1] = particle->Py(); p[2] = particle->Pz();
290                                                 AliAODTrack *aodTrack = new((*fAODCTS)[indexCh++]) 
291                                                 AliAODTrack(0, iParticle, p, kTRUE, x, kFALSE,NULL, 0, 0, 
292                                                                         NULL,
293                                                                         0x0,//primary,
294                                                                         kFALSE, // No fit performed
295                                                                         kFALSE, // No fit performed
296                                                                         AliAODTrack::kPrimary, 
297                                                                         0);
298                                                 SetTrackChargeAndPID(pdg, aodTrack);
299                                         }
300                                 }
301                                 //Keep some charged particles in calorimeters lists
302                                 if((fFillPHOS || fFillEMCAL) && KeepChargedParticles(pdg)) FillCalorimeters(iParticle, particle, momentum, indexPHOS, indexEMCAL);
303                                 
304                         }//Charged
305                         
306                         //-------------Neutral particles ----------------------
307                         else if(charge == 0 && (fFillPHOS || fFillEMCAL)){
308                                 //Skip neutrinos or other neutral particles
309                                 //if(SkipNeutralParticles(pdg) || particle->IsPrimary()) continue ; // protection added (MG)
310                                 if(SkipNeutralParticles(pdg)) continue ;
311                                 //Fill particle/calocluster arrays
312                                 if(!fDecayPi0) {
313                                         FillCalorimeters(iParticle, particle, momentum, indexPHOS, indexEMCAL);
314                                 }
315                                 else {
316                                         //Sometimes pi0 are stable for the generator, if needed decay it by hand
317                                         if(pdg == 111 ){
318                                                 if(momentum.Pt() >  fPHOSPtMin || momentum.Pt() >  fEMCALPtMin){
319                                                         TLorentzVector lvGamma1, lvGamma2 ;
320                                                         //Double_t angle = 0;
321                                                         
322                                                         //Decay
323                                                         MakePi0Decay(momentum,lvGamma1,lvGamma2);//,angle);
324                                                         
325                                                         //Check if Pi0 is in the acceptance of the calorimeters, if aperture angle is small, keep it
326                                                         TParticle * pPhoton1 = new TParticle(22,1,iParticle,0,0,0,lvGamma1.Px(),lvGamma1.Py(),
327                                                                                                                                  lvGamma1.Pz(),lvGamma1.E(),0,0,0,0);   
328                                                         TParticle * pPhoton2 = new TParticle(22,1,iParticle,0,0,0,lvGamma2.Px(),lvGamma2.Py(),
329                                                                                                                                  lvGamma2.Pz(),lvGamma2.E(),0,0,0,0);
330                                                         //Fill particle/calocluster arrays
331                                                         FillCalorimeters(iParticle,pPhoton1,lvGamma1, indexPHOS, indexEMCAL);
332                                                         FillCalorimeters(iParticle,pPhoton2,lvGamma2, indexPHOS, indexEMCAL);
333                                                 }//pt cut
334                                         }//pi0
335                                         else FillCalorimeters(iParticle,particle, momentum, indexPHOS, indexEMCAL); //Add the rest
336                                 }
337                         }//neutral particles
338                 } //particle with correct status
339         }//particle loop
340
341         fIndex2ndPhoton = -1; //In case of overlapping studies, reset for each event    
342 }
343
344 //________________________________________________________________
345 Bool_t AliCaloTrackMCReader::KeepParticleWithStatus(Int_t status) const {
346   //Check if status is equal to one of the  list
347   //These particles will be used in analysis.
348   if(!fKeepAllStatus){
349     for(Int_t i= 0; i < fStatusArray->GetSize(); i++)
350       if(status ==  fStatusArray->At(i)) return kTRUE ;
351
352     return kFALSE; 
353     
354   }
355   else
356     return kTRUE ;  
357 }
358
359 //________________________________________________________________
360 Bool_t AliCaloTrackMCReader::KeepChargedParticles(Int_t pdg) const {
361   //Check if pdg is equal to one of the charged particles list
362   //These particles will be added to the calorimeters lists.
363
364   for(Int_t i= 0; i < fChargedParticlesArray->GetSize(); i++)
365     if(TMath::Abs(pdg) ==  fChargedParticlesArray->At(i)) return kTRUE ;
366   
367   return kFALSE; 
368   
369 }
370
371 //________________________________________________________________
372 void AliCaloTrackMCReader::Print(const Option_t * opt) const
373 {
374   //Print some relevant parameters set for the analysis
375   if(! opt)
376     return;
377   
378   Info("**** Print **** ", "%s %s", GetName(), GetTitle() ) ;
379   
380   printf("Decay Pi0?          : %d\n", fDecayPi0) ;
381   printf("TClonesArray type   : %d\n", fClonesArrayType) ;
382   printf("Check Overlap in Calo?    : %d\n", fCheckOverlap) ;
383   printf("Keep all status?    : %d\n", fKeepAllStatus) ;
384   
385   if(!fKeepAllStatus) printf("Keep particles with status : ");
386   for(Int_t i= 0; i < fStatusArray->GetSize(); i++)
387     printf(" %d ; ", fStatusArray->At(i));
388   printf("\n");
389   
390   printf("Skip neutral particles in calo : ");
391   for(Int_t i= 0; i < fNeutralParticlesArray->GetSize(); i++)
392     printf(" %d ; ", fNeutralParticlesArray->At(i));
393   printf("\n");
394   
395   printf("Keep charged particles in calo : ");
396   for(Int_t i= 0; i < fChargedParticlesArray->GetSize(); i++)
397     printf(" %d ; ", fChargedParticlesArray->At(i));
398   printf("\n");
399
400 }
401
402 //____________________________________________________________________________
403 void AliCaloTrackMCReader::MakePi0Decay(TLorentzVector &p0, TLorentzVector &p1, 
404                                 TLorentzVector &p2) const {//, Double_t &angle) {
405   // Perform isotropic decay pi0 -> 2 photons
406   // p0 is pi0 4-momentum (inut)
407   // p1 and p2 are photon 4-momenta (output)
408   //  cout<<"Boost vector"<<endl;
409   Double_t mPi0 = TDatabasePDG::Instance()->GetParticle(111)->Mass();
410   TVector3 b = p0.BoostVector();
411   //cout<<"Parameters"<<endl;
412   //Double_t mPi0   = p0.M();
413   Double_t phi    = TMath::TwoPi() * gRandom->Rndm();
414   Double_t cosThe = 2 * gRandom->Rndm() - 1;
415   Double_t cosPhi = TMath::Cos(phi);
416   Double_t sinPhi = TMath::Sin(phi);
417   Double_t sinThe = TMath::Sqrt(1-cosThe*cosThe);
418   Double_t ePi0   = mPi0/2.;
419   //cout<<"ePi0 "<<ePi0<<endl;
420   //cout<<"Components"<<endl;
421   p1.SetPx(+ePi0*cosPhi*sinThe);
422   p1.SetPy(+ePi0*sinPhi*sinThe);
423   p1.SetPz(+ePi0*cosThe);
424   p1.SetE(ePi0);
425   //cout<<"p1: "<<p1.Px()<<" "<<p1.Py()<<" "<<p1.Pz()<<" "<<p1.E()<<endl;
426   //cout<<"p1 Mass: "<<p1.Px()*p1.Px()+p1.Py()*p1.Py()+p1.Pz()*p1.Pz()-p1.E()*p1.E()<<endl;
427   p2.SetPx(-ePi0*cosPhi*sinThe);
428   p2.SetPy(-ePi0*sinPhi*sinThe);
429   p2.SetPz(-ePi0*cosThe);
430   p2.SetE(ePi0);
431   //cout<<"p2: "<<p2.Px()<<" "<<p2.Py()<<" "<<p2.Pz()<<" "<<p2.E()<<endl;
432   //cout<<"p2 Mass: "<<p2.Px()*p2.Px()+p2.Py()*p2.Py()+p2.Pz()*p2.Pz()-p2.E()*p2.E()<<endl;
433   //cout<<"Boost "<<b.X()<<" "<<b.Y()<<" "<<b.Z()<<endl;
434   p1.Boost(b);
435   //cout<<"p1: "<<p1.Px()<<" "<<p1.Py()<<" "<<p1.Pz()<<" "<<p1.E()<<endl;
436   p2.Boost(b);
437   //cout<<"p2: "<<p2.Px()<<" "<<p2.Py()<<" "<<p2.Pz()<<" "<<p2.E()<<endl;
438   //cout<<"angle"<<endl;
439   //angle = p1.Angle(p2.Vect());
440   //cout<<angle<<endl;
441 }
442
443 //____________________________________________________________________________
444 void AliCaloTrackMCReader::SetInputEvent(TObject* /*esd*/, TObject* aod, TObject* mc) {
445   // Connect the data pointer
446   SetMC((AliMCEvent*) mc);
447   SetAOD((AliAODEvent*) aod);
448 }
449
450
451 //________________________________________________________________
452 Bool_t AliCaloTrackMCReader::SkipNeutralParticles(Int_t pdg) const {
453   //Check if pdg is equal to one of the neutral particles list
454   //These particles will be skipped from analysis.
455
456   for(Int_t i= 0; i < fNeutralParticlesArray->GetSize(); i++)
457     if(TMath::Abs(pdg) ==  fNeutralParticlesArray->At(i)) return kTRUE ;
458   
459   return kFALSE; 
460   
461 }
462
463
464 //____________________________________________________________________
465 void AliCaloTrackMCReader::SetTrackChargeAndPID(const Int_t pdgCode, AliAODTrack *track) const {
466 //Give a PID weight for tracks equal to 1 depending on the particle type
467
468   Float_t pid[10] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
469
470   switch (pdgCode) {
471
472   case 22: // gamma
473     track->SetCharge(0);
474     pid[AliAODTrack::kUnknown] = 1.;
475     track->SetPID(pid);
476     break;
477
478   case 11: // e- 
479     track->SetCharge(-1);
480     pid[AliAODTrack::kElectron] = 1.;
481     track->SetPID(pid);
482     break;
483     
484   case -11: // e+
485     track->SetCharge(+1);
486     pid[AliAODTrack::kElectron] = 1.;
487     track->SetPID(pid);
488     break;
489     
490   case 13: // mu- 
491     track->SetCharge(-1);
492     pid[AliAODTrack::kMuon] = 1.;
493     track->SetPID(pid);
494     break;
495     
496   case -13: // mu+
497     track->SetCharge(+1);
498     pid[AliAODTrack::kMuon] = 1.;
499     track->SetPID(pid);
500     break;
501     
502   case 111: // pi0
503     track->SetCharge(0);
504     pid[AliAODTrack::kUnknown] = 1.;
505     track->SetPID(pid);
506     break;
507     
508   case 211: // pi+
509     track->SetCharge(+1);
510     pid[AliAODTrack::kPion] = 1.;
511     track->SetPID(pid);
512     break;
513     
514   case -211: // pi-
515     track->SetCharge(-1);
516     pid[AliAODTrack::kPion] = 1.;
517     track->SetPID(pid);
518     break;
519     
520   case 130: // K0L
521     track->SetCharge(0);
522     pid[AliAODTrack::kUnknown] = 1.;
523     track->SetPID(pid);
524     break;
525     
526   case 321: // K+
527     track->SetCharge(+1);
528     pid[AliAODTrack::kKaon] = 1.;
529     track->SetPID(pid);
530     break;
531     
532   case -321: // K- 
533     track->SetCharge(-1);
534     pid[AliAODTrack::kKaon] = 1.;
535     track->SetPID(pid);
536     break;
537     
538   case 2112: // n
539     track->SetCharge(0);
540     pid[AliAODTrack::kUnknown] = 1.;
541     track->SetPID(pid);
542     break;
543     
544   case 2212: // p
545     track->SetCharge(+1);
546     pid[AliAODTrack::kProton] = 1.;
547     track->SetPID(pid);
548     break;
549     
550   case -2212: // anti-p
551     track->SetCharge(-1);
552     pid[AliAODTrack::kProton] = 1.;
553     track->SetPID(pid);
554     break;
555
556   case 310: // K0S
557     track->SetCharge(0);
558     pid[AliAODTrack::kUnknown] = 1.;
559     track->SetPID(pid);
560     break;
561     
562   case 311: // K0
563     track->SetCharge(0);
564     pid[AliAODTrack::kUnknown] = 1.;
565     track->SetPID(pid);
566     break;
567     
568   case -311: // anti-K0
569     track->SetCharge(0);
570     pid[AliAODTrack::kUnknown] = 1.;
571     track->SetPID(pid);
572     break;
573     
574   case 221: // eta
575     track->SetCharge(0);
576     pid[AliAODTrack::kUnknown] = 1.;
577     track->SetPID(pid);
578     break;
579
580   case 3122: // lambda
581     track->SetCharge(0);
582     pid[AliAODTrack::kUnknown] = 1.;
583     track->SetPID(pid);
584     break;
585
586   case 3222: // Sigma+
587     track->SetCharge(+1);
588     pid[AliAODTrack::kUnknown] = 1.;
589     track->SetPID(pid);
590     break;
591
592   case 3212: // Sigma0
593     track->SetCharge(-1);
594     pid[AliAODTrack::kUnknown] = 1.;
595     track->SetPID(pid);
596     break;
597
598   case 3112: // Sigma-
599     track->SetCharge(-1);
600     pid[AliAODTrack::kUnknown] = 1.;
601     track->SetPID(pid);
602     break;
603
604   case 3322: // Xi0
605     track->SetCharge(0);
606     pid[AliAODTrack::kUnknown] = 1.;
607     track->SetPID(pid);
608     break;
609
610   case 3312: // Xi-
611     track->SetCharge(-1);
612     pid[AliAODTrack::kUnknown] = 1.;
613     track->SetPID(pid);
614     break;
615
616   case 3334: // Omega-
617     track->SetCharge(-1);
618     pid[AliAODTrack::kUnknown] = 1.;
619     track->SetPID(pid);
620     break;
621
622   case -2112: // n-bar
623     track->SetCharge(0);
624     pid[AliAODTrack::kUnknown] = 1.;
625     track->SetPID(pid);
626     break;
627
628   case -3122: // anti-Lambda
629     track->SetCharge(0);
630     pid[AliAODTrack::kUnknown] = 1.;
631     track->SetPID(pid);
632     break;
633
634   case -3222: // anti-Sigma-
635     track->SetCharge(-1);
636     pid[AliAODTrack::kUnknown] = 1.;
637     track->SetPID(pid);
638     break;
639
640   case -3212: // anti-Sigma0
641     track->SetCharge(0);
642     pid[AliAODTrack::kUnknown] = 1.;
643     track->SetPID(pid);
644     break;
645
646   case -3112: // anti-Sigma+
647     track->SetCharge(+1);
648     pid[AliAODTrack::kUnknown] = 1.;
649     track->SetPID(pid);
650     break;
651
652   case -3322: // anti-Xi0
653     track->SetCharge(0);
654     pid[AliAODTrack::kUnknown] = 1.;
655     track->SetPID(pid);
656     break;
657
658   case -3312: // anti-Xi+
659     track->SetCharge(+1);
660     break;
661
662   case -3334: // anti-Omega+
663     track->SetCharge(+1);
664     pid[AliAODTrack::kUnknown] = 1.;
665     track->SetPID(pid);
666     break;
667
668   case 411: // D+
669     track->SetCharge(+1);
670     pid[AliAODTrack::kUnknown] = 1.;
671     track->SetPID(pid);
672     break;
673
674   case -411: // D- 
675     track->SetCharge(-1);
676     pid[AliAODTrack::kUnknown] = 1.;
677     track->SetPID(pid);
678     break;
679
680   case 421: // D0
681     track->SetCharge(0);
682     pid[AliAODTrack::kUnknown] = 1.;
683     track->SetPID(pid);
684     break;
685
686   case -421: // anti-D0
687     track->SetCharge(0);
688     pid[AliAODTrack::kUnknown] = 1.;
689     track->SetPID(pid);
690     break;
691
692   default : // unknown
693     track->SetCharge(-99);
694     pid[AliAODTrack::kUnknown] = 1.;
695     track->SetPID(pid);
696  }
697
698   track->SetPID(pid);
699
700   return;
701 }
702
703 //____________________________________________________________________
704 void AliCaloTrackMCReader::SetCaloClusterPID(const Int_t pdgCode, AliAODCaloCluster *calo) const {
705 //Give a PID weight for CaloClusters equal to 1 depending on the particle type
706
707   Float_t pid[13] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
708
709   switch (pdgCode) {
710
711   case 22: // gamma
712     pid[AliAODCaloCluster::kPhoton] = 1.;
713     calo->SetPID(pid);
714     break;
715
716   case 11: // e- 
717     pid[AliAODCaloCluster::kElectron] = 1.;
718     calo->SetPID(pid);
719     break;
720     
721   case -11: // e+
722     pid[AliAODCaloCluster::kElectron] = 1.;
723     calo->SetPID(pid);
724     break;
725     
726   case 13: // mu- 
727     pid[AliAODCaloCluster::kCharged] = 1.;
728     calo->SetPID(pid);
729     break;
730     
731   case -13: // mu+
732     pid[AliAODCaloCluster::kCharged] = 1.;
733     calo->SetPID(pid);
734     break;
735     
736   case 111: // pi0
737     pid[AliAODCaloCluster::kPi0] = 1.;
738     calo->SetPID(pid);
739     break;
740     
741   case 211: // pi+
742     pid[AliAODCaloCluster::kCharged] = 1.;
743     calo->SetPID(pid);
744     break;
745     
746   case -211: // pi-
747     pid[AliAODCaloCluster::kCharged] = 1.;
748     calo->SetPID(pid);
749     break;
750     
751   case 130: // K0L
752     pid[AliAODCaloCluster::kKaon0] = 1.;
753     pid[AliAODCaloCluster::kNeutral] = 1;
754     calo->SetPID(pid);
755     break;
756     
757   case 321: // K+
758     pid[AliAODCaloCluster::kCharged] = 1.;
759     calo->SetPID(pid);
760     break;
761     
762   case -321: // K- 
763     pid[AliAODCaloCluster::kCharged] = 1.;
764     calo->SetPID(pid);
765     break;
766     
767   case 2112: // n
768     pid[AliAODCaloCluster::kNeutron] = 1.;
769     pid[AliAODCaloCluster::kNeutral] = 1.;
770     calo->SetPID(pid);
771     break;
772     
773   case 2212: // p
774     pid[AliAODCaloCluster::kCharged] = 1.;
775     calo->SetPID(pid);
776     break;
777     
778   case -2212: // anti-p
779     pid[AliAODCaloCluster::kCharged] = 1.;
780     calo->SetPID(pid);
781     break;
782
783   case 310: // K0S
784     pid[AliAODCaloCluster::kKaon0] = 1.;
785     pid[AliAODCaloCluster::kNeutral] = 1.;
786     calo->SetPID(pid);
787     break;
788     
789   case 311: // K0
790     pid[AliAODCaloCluster::kKaon0] = 1.;
791     pid[AliAODCaloCluster::kNeutral] = 1.;
792     calo->SetPID(pid);
793     break;
794     
795   case -311: // anti-K0
796     pid[AliAODCaloCluster::kKaon0] = 1.;
797     pid[AliAODCaloCluster::kNeutral] = 1.;
798     calo->SetPID(pid);
799     break;
800     
801   case 221: // eta
802     pid[AliAODCaloCluster::kNeutral] = 1.;
803     calo->SetPID(pid);
804     break;
805
806   case 3122: // lambda
807     pid[AliAODCaloCluster::kUnknown] = 1.;
808     calo->SetPID(pid);
809     break;
810
811   case 3222: // Sigma+
812     pid[AliAODCaloCluster::kUnknown] = 1.;
813     calo->SetPID(pid);
814     break;
815
816   case 3212: // Sigma0
817     pid[AliAODCaloCluster::kUnknown] = 1.;
818     calo->SetPID(pid);
819     break;
820
821   case 3112: // Sigma-
822     pid[AliAODCaloCluster::kUnknown] = 1.;
823     calo->SetPID(pid);
824     break;
825
826   case 3322: // Xi0
827     pid[AliAODCaloCluster::kUnknown] = 1.;
828     calo->SetPID(pid);
829     break;
830
831   case 3312: // Xi-
832     pid[AliAODCaloCluster::kUnknown] = 1.;
833     calo->SetPID(pid);
834     break;
835
836   case 3334: // Omega-
837     pid[AliAODCaloCluster::kUnknown] = 1.;
838     calo->SetPID(pid);
839     break;
840
841   case -2112: // n-bar
842     pid[AliAODCaloCluster::kNeutron] = 1.;
843     pid[AliAODCaloCluster::kNeutral] = 1.;
844     calo->SetPID(pid);
845     break;
846
847   case -3122: // anti-Lambda
848     pid[AliAODCaloCluster::kUnknown] = 1.;
849     calo->SetPID(pid);
850     break;
851
852   case -3222: // anti-Sigma-
853     pid[AliAODCaloCluster::kUnknown] = 1.;
854     calo->SetPID(pid);
855     break;
856
857   case -3212: // anti-Sigma0
858     pid[AliAODCaloCluster::kUnknown] = 1.;
859     calo->SetPID(pid);
860     break;
861
862   case -3112: // anti-Sigma+
863     pid[AliAODCaloCluster::kUnknown] = 1.;
864     calo->SetPID(pid);
865     break;
866
867   case -3322: // anti-Xi0
868     pid[AliAODCaloCluster::kUnknown] = 1.;
869     calo->SetPID(pid);
870     break;
871
872   case -3312: // anti-Xi+
873     pid[AliAODCaloCluster::kUnknown] = 1.;
874     calo->SetPID(pid);
875     break;
876
877   case -3334: // anti-Omega+
878     pid[AliAODCaloCluster::kUnknown] = 1.;
879     calo->SetPID(pid);
880     break;
881
882   case 411: // D+
883     pid[AliAODCaloCluster::kUnknown] = 1.;
884     calo->SetPID(pid);
885     break;
886
887   case -411: // D- 
888     pid[AliAODCaloCluster::kUnknown] = 1.;
889     calo->SetPID(pid);
890     break;
891
892   case 421: // D0
893     pid[AliAODCaloCluster::kUnknown] = 1.;
894     calo->SetPID(pid);
895     break;
896
897   case -421: // anti-D0
898     pid[AliAODCaloCluster::kUnknown] = 1.;
899     calo->SetPID(pid);
900     break;
901
902   default : // unknown
903     pid[AliAODCaloCluster::kUnknown] = 1.;
904     calo->SetPID(pid);
905  }
906
907  
908   return;
909 }