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