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