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