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