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