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