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