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