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