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