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