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