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