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