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