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