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