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