PWG4/AliGammaMCReader.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 /* History of cvs commits:
  18  *
  19  * $Log$
  20  * Revision 1.4  2007/10/29 13:48:42  gustavo
  21  * Corrected coding violations
  22  *
  23  *
  24  */
  25
  26 //_________________________________________________________________________
  27 // Class for reading data (Kinematics) in order to do prompt gamma correlations
  28 //  Class created from old AliPHOSGammaJet
  29 //  (see AliRoot versions previous Release 4-09)
  30 //
  31 //*-- Author: Gustavo Conesa (LNF-INFN)
  32 //////////////////////////////////////////////////////////////////////////////
  33
  34
  35 // --- ROOT system ---
  36
  37 #include <TFile.h>
  38 #include <TParticle.h>
  39 #include <TH2.h>
  40 #include <TChain.h>
  41 #include <TRandom.h>
  42 #include <TArrayI.h>
  43 #include <TClonesArray.h>
  44
  45 //---- ANALYSIS system ----
  46 #include "AliGammaMCReader.h"
  47 #include "Riostream.h"
  48 #include "AliLog.h"
  49 #include "AliStack.h"
  50
  51 ClassImp(AliGammaMCReader)
  52
  53 //____________________________________________________________________________
  54 AliGammaMCReader::AliGammaMCReader() :
  55   AliGammaReader(), fDecayPi0(0), fCheckOverlapping(kFALSE),  fNeutralParticlesArray(0x0)
  56 {
  57   //Ctor
  58
  59   //Initialize parameters
  60   InitParameters();
  61   fDataType = kMC;
  62
  63 }
  64
  65 //____________________________________________________________________________
  66 AliGammaMCReader::AliGammaMCReader(const AliGammaMCReader & g) :
  67   AliGammaReader(g), fDecayPi0(g.fDecayPi0), fCheckOverlapping(g.fCheckOverlapping),
  68   fNeutralParticlesArray(g.fNeutralParticlesArray?new TArrayI(*g.fNeutralParticlesArray):0x0)
  69 {
  70   // cpy ctor
  71 }
  72
  73 //_________________________________________________________________________
  74 AliGammaMCReader & AliGammaMCReader::operator = (const AliGammaMCReader & source)
  75 {
  76   // assignment operator
  77
  78   if(&source == this) return *this;
  79
  80   fDecayPi0 = source.fDecayPi0;
  81   fCheckOverlapping = source.fCheckOverlapping;
  82   delete fNeutralParticlesArray;
  83   fNeutralParticlesArray = source.fNeutralParticlesArray?new TArrayI(*source.fNeutralParticlesArray):0x0;
  84   return *this;
  85
  86 }
  87
  88 //_________________________________
  89 AliGammaMCReader::~AliGammaMCReader() {
  90   //Dtor
  91
  92   delete fNeutralParticlesArray ;
  93
  94 }
  95
  96 //___________________________________________________________________________
  97 void AliGammaMCReader::CaseDecayGamma(Int_t index, TParticle * particle, AliStack * sta,
  98                                   TClonesArray * plEMCAL, Int_t &indexEMCAL,
  99                                   TClonesArray * plPHOS, Int_t &indexPHOS){
 100   //In case pi0 are decayed by pythia, check if mother is pi0 and in such case look if
 101   //there is overlapping. Send particle=pi0 if there is overlapping.
 102
 103   TParticle * pmother =sta->Particle(particle->GetFirstMother());
 104   if(pmother->GetPdgCode() == 111 && pmother->GetNDaughters() == 2) {//Do not consider 3 particle decay case
 105     Int_t idaug0 = pmother->GetDaughter(0);
 106     Int_t idaug1 = pmother->GetDaughter(1);
 107     TParticle * pdaug0 = sta -> Particle(idaug0);
 108     TParticle * pdaug1 = sta -> Particle(idaug1);
 109
 110     if((index ==  idaug0 &&  pdaug0->Pt() > fNeutralPtCut) ||
 111        (index ==  idaug1 &&  pdaug0->Pt() <= fNeutralPtCut)){//Check decay when first daughter arrives, do nothing with second.
 112
 113       FillListWithDecayGammaOrPi0(pmother, pdaug0, pdaug1, plEMCAL, indexEMCAL, plPHOS, indexPHOS);
 114
 115     }
 116   }//mother is a pi0 with 2 daughters
 117   else{
 118
 119     if(IsInPHOS(particle->Phi(),particle->Eta())) {
 120       TParticle * pmother =sta->Particle(particle->GetFirstMother());
 121       SetPhotonStatus(particle,pmother);
 122       new((*plPHOS)[indexPHOS++])  TParticle(*particle) ;
 123     }
 124     else if(IsInEMCAL(particle->Phi(),particle->Eta())){
 125       TParticle * pmother =sta->Particle(particle->GetFirstMother());
 126       SetPhotonStatus(particle,pmother);
 127       new((*plEMCAL)[indexEMCAL++])  TParticle(*particle) ;
 128     }
 129   }
 130 }
 131
 132 //___________________________________________________________________________
 133  void AliGammaMCReader::CaseGeantDecay(TParticle * particle, AliStack * stack,
 134                                    TClonesArray * plEMCAL, Int_t &indexEMCAL,
 135                                    TClonesArray * plPHOS, Int_t &indexPHOS){
 136    //Find decay gamma from pi0, decayed by GEANT and put them in the list.
 137
 138    Int_t ndaug = particle->GetNDaughters() ;
 139    TParticle * d1 = new TParticle();
 140    TParticle * d2 = new TParticle();
 141
 142    if(ndaug > 0){//At least 1 daugther
 143      d1 = stack->Particle(particle->GetDaughter(0));
 144      if (ndaug > 1 )
 145        d2 = stack->Particle(particle->GetDaughter(1));
 146
 147      FillListWithDecayGammaOrPi0(particle, d1, d2, plEMCAL, indexEMCAL, plPHOS, indexPHOS);
 148
 149    }// ndaugh > 0
 150  }
 151
 152 //___________________________________________________________________________
 153 void AliGammaMCReader::CasePi0Decay(TParticle * pPi0,
 154                                     TClonesArray * plEMCAL, Int_t &indexEMCAL,
 155                                     TClonesArray * plPHOS, Int_t &indexPHOS){
 156
 157   //Decays pi0, see if aperture angle is small and then add the pi0 or the 2 gamma
 158
 159   TLorentzVector lvPi0, lvGamma1, lvGamma2 ;
 160   //Double_t angle = 0;
 161
 162   lvPi0.SetPxPyPzE(pPi0->Px(),pPi0->Py(),pPi0->Pz(),pPi0->Energy());
 163
 164   //Decay
 165   MakePi0Decay(lvPi0,lvGamma1,lvGamma2);//,angle);
 166
 167   //Check if Pi0 is in the acceptance of the calorimeters, if aperture angle is small, keep it
 168   TParticle * pPhoton1 = new TParticle(22,1,0,0,0,0,lvGamma1.Px(),lvGamma1.Py(),
 169                                       lvGamma1.Pz(),lvGamma1.E(),0,0,0,0);
 170   TParticle * pPhoton2 = new TParticle(22,1,0,0,0,0,lvGamma2.Px(),lvGamma2.Py(),
 171                                       lvGamma2.Pz(),lvGamma2.E(),0,0,0,0);
 172
 173   FillListWithDecayGammaOrPi0(pPi0, pPhoton1, pPhoton2, plEMCAL, indexEMCAL, plPHOS, indexPHOS);
 174
 175 }
 176
 177 //_______________________________________________________________
 178 void AliGammaMCReader::InitParameters()
 179 {
 180
 181   //Initialize the parameters of the analysis.
 182
 183   fDecayPi0 = kGeantDecay;
 184   fCheckOverlapping = kTRUE ;
 185
 186   Int_t pdgarray[]={12,14,16};// skip neutrinos
 187   fNeutralParticlesArray = new TArrayI(3, pdgarray);
 188
 189 }
 190
 191 //____________________________________________________________________________
 192 void AliGammaMCReader::CreateParticleList(TObject * data, TObject *,
 193                                           TClonesArray * plCh,
 194                                           TClonesArray * plEMCAL,
 195                                           TClonesArray * plPHOS,
 196                                           TClonesArray *plParton,TClonesArray *,TClonesArray *)
 197 {
 198   //Create list of particles from EMCAL, PHOS and CTS.
 199   AliStack * stack = (AliStack *) data ;
 200   Int_t indexCh     = plCh->GetEntries() ;
 201   Int_t indexEMCAL = plEMCAL->GetEntries() ;
 202   Int_t indexPHOS = plPHOS->GetEntries() ;
 203   Int_t indexParton = plParton->GetEntries() ;
 204   Int_t iParticle = 0 ;
 205   Double_t charge = 0.;
 206
 207   for (iParticle=0 ; iParticle <  stack->GetNprimary() ; iParticle++) {
 208     TParticle * particle = stack->Particle(iParticle);
 209
 210     //Keep partons
 211     if(particle->GetStatusCode() == 21 && iParticle>=2){//All partons, not nucleus
 212       new((*plParton)[indexParton++])  TParticle(*particle) ;
 213     }
 214
 215     //Keep Stable particles
 216     if((particle->GetStatusCode() == 1) && (particle->Pt() > 0)){
 217
 218       charge = TDatabasePDG::Instance()->GetParticle(particle->GetPdgCode())->Charge();
 219
 220       //---------- Charged particles ----------------------
 221       if( fSwitchOnCTS && (charge != 0) && (particle->Pt() > fChargedPtCut)){
 222         //Particles in CTS acceptance
 223         if(TMath::Abs(particle->Eta())<fCTSEtaCut){
 224           //Fill lists
 225           new((*plCh)[indexCh++])       TParticle(*particle) ;
 226         }
 227       }
 228
 229       //-------------Neutral particles ----------------------
 230       else if((charge == 0) && particle->Pt() > fNeutralPtCut ){ //&&
 231         //TMath::Abs(particle->GetPdgCode())>16){//Avoid neutrinos
 232
 233         if(particle->GetPdgCode()!=111){
 234           //In case that we force PYTHIA to decay pi0, and we want to check the overlapping of
 235           // the decay gamma.
 236           if(particle->GetPdgCode() == 22 && fDecayPi0==kDecayGamma){
 237             CaseDecayGamma(iParticle,particle, stack,plEMCAL, indexEMCAL, plPHOS, indexPHOS); //If pythia decays pi0
 238           }
 239           else{
 240             //Take out some particles like neutrinos
 241             if(!SkipNeutralParticles(particle->GetPdgCode()) && !particle->IsPrimary()){ // protection added (MG)
 242               TParticle * pmother =stack->Particle(particle->GetFirstMother());
 243               if(IsInPHOS(particle->Phi(),particle->Eta())){
 244                 if(particle->GetPdgCode()==22)SetPhotonStatus(particle,pmother);
 245                 new((*plPHOS)[indexPHOS++])  TParticle(*particle) ;
 246               }
 247               else if(IsInEMCAL(particle->Phi(),particle->Eta())){
 248                 if(particle->GetPdgCode()==22) SetPhotonStatus(particle,pmother);
 249                 new((*plEMCAL)[indexEMCAL++])  TParticle(*particle) ;
 250               }
 251             }//skip neutral particles
 252           }//Case kDecayGamma
 253         }//no pi0
 254         else{
 255           if(fDecayPi0 == kNoDecay){//keep the pi0 do not decay
 256             if(IsInPHOS(particle->Phi(),particle->Eta()))
 257               new((*plPHOS)[indexPHOS++])  TParticle(*particle) ;
 258             else if(IsInEMCAL(particle->Phi(),particle->Eta()))
 259               new((*plEMCAL)[indexEMCAL++])  TParticle(*particle) ;
 260           }//nodecay
 261           else if(fDecayPi0 == kDecay)
 262             CasePi0Decay(particle,plEMCAL,indexEMCAL,plPHOS, indexPHOS);
 263           else if(fDecayPi0 == kGeantDecay)
 264             CaseGeantDecay(particle, stack,plEMCAL, indexEMCAL, plPHOS, indexPHOS);
 265         }//pi0
 266       }//neutral particle
 267     }//stable particle
 268   }//particle loop
 269 }
 270
 271
 272 //___________________________________________________________________________
 273 void AliGammaMCReader::FillListWithDecayGammaOrPi0(TParticle * pPi0, TParticle * pdaug0, TParticle * pdaug1,
 274                                    TClonesArray * plEMCAL, Int_t &indexEMCAL,
 275                                    TClonesArray * plPHOS, Int_t &indexPHOS){
 276
 277   //Check if decay gamma overlapp in calorimeter, in such case keep the pi0, if not keep both photons.
 278   Bool_t  overlap = kFALSE ;
 279   TLorentzVector lv1 , lv2 ;
 280   pdaug0->Momentum(lv1);
 281   pdaug1->Momentum(lv2);
 282   Double_t angle = lv1.Angle(lv2.Vect());
 283
 284   if(IsInEMCAL(pPi0->Phi(), pPi0->Eta()))
 285     {
 286       if (angle < fEMCALMinAngle){
 287         new((*plEMCAL)[indexEMCAL++])       TParticle(*pPi0) ;
 288         overlap = kTRUE;
 289       }
 290     }
 291
 292   else if(IsInPHOS(pPi0->Phi(), pPi0->Eta())){
 293     if (angle < fPHOSMinAngle){
 294       new((*plPHOS)[indexPHOS++])       TParticle(*pPi0) ;
 295       overlap = kTRUE;
 296     }
 297   }//fCheckOverlapping
 298
 299   //Fill with gammas if not overlapp
 300   if(!overlap){
 301     if(pdaug0->GetPdgCode() == 22 || TMath::Abs(pdaug0->GetPdgCode() ) == 11 ){
 302       pdaug0->SetStatusCode(kPi0DecayPhoton);
 303       if(IsInEMCAL(pdaug0->Phi(),pdaug0->Eta()) &&  pdaug0->Pt() > fNeutralPtCut)
 304         new((*plEMCAL)[indexEMCAL++])       TParticle(*pdaug0) ;
 305       else if(IsInPHOS(pdaug0->Phi(),pdaug0->Eta()) &&  pdaug0->Pt() > fNeutralPtCut)
 306         new((*plPHOS)[indexPHOS++])       TParticle(*pdaug0) ;
 307     }
 308
 309     if(pdaug1->GetPdgCode() == 22 || TMath::Abs(pdaug1->GetPdgCode() ) == 11 ){
 310       pdaug1->SetStatusCode(kPi0DecayPhoton);
 311       if(IsInEMCAL(pdaug1->Phi(),pdaug1->Eta()) &&  pdaug1->Pt() > fNeutralPtCut)
 312         new((*plEMCAL)[indexEMCAL++])       TParticle(*pdaug1) ;
 313       else if(IsInPHOS(pdaug1->Phi(),pdaug1->Eta()) &&  pdaug1->Pt() > fNeutralPtCut)
 314         new((*plPHOS)[indexPHOS++])       TParticle(*pdaug1) ;
 315     }
 316   }// overlap?
 317
 318  }
 319
 320
 321
 322 //___________________________________________________________________________
 323 Bool_t  AliGammaMCReader::IsInEMCAL(Double_t phi, Double_t eta) const {
 324   //Check if particle is in EMCAL acceptance
 325
 326   if(fSwitchOnEMCAL){
 327     if(phi<0)
 328       phi+=TMath::TwoPi();
 329     if( phi > fPhiEMCALCut[0] && phi < fPhiEMCALCut[1] &&
 330         TMath::Abs(eta)<fEMCALEtaCut) return kTRUE ;
 331     else  return kFALSE;
 332   }
 333   return kFALSE ;
 334 }
 335
 336 //___________________________________________________________________________
 337 Bool_t  AliGammaMCReader::IsInPHOS(Double_t phi, Double_t eta) const {
 338   //Check if particle is in PHOS acceptance
 339
 340   if(fSwitchOnPHOS){
 341     if(phi<0)
 342       phi+=TMath::TwoPi();
 343     if( phi > fPhiPHOSCut[0] && phi < fPhiPHOSCut[1] &&
 344         TMath::Abs(eta)<fPHOSEtaCut) return kTRUE ;
 345     else  return kFALSE;
 346   }
 347
 348   return kFALSE ;
 349 }
 350
 351 //____________________________________________________________________________
 352 void AliGammaMCReader::MakePi0Decay(TLorentzVector &p0, TLorentzVector &p1,
 353                                 TLorentzVector &p2){//, Double_t &angle) {
 354   // Perform isotropic decay pi0 -> 2 photons
 355   // p0 is pi0 4-momentum (inut)
 356   // p1 and p2 are photon 4-momenta (output)
 357   //  cout<<"Boost vector"<<endl;
 358   Double_t mPi0 = TDatabasePDG::Instance()->GetParticle(111)->Mass();
 359   TVector3 b = p0.BoostVector();
 360   //cout<<"Parameters"<<endl;
 361   //Double_t mPi0   = p0.M();
 362   Double_t phi    = TMath::TwoPi() * gRandom->Rndm();
 363   Double_t cosThe = 2 * gRandom->Rndm() - 1;
 364   Double_t cosPhi = TMath::Cos(phi);
 365   Double_t sinPhi = TMath::Sin(phi);
 366   Double_t sinThe = TMath::Sqrt(1-cosThe*cosThe);
 367   Double_t ePi0   = mPi0/2.;
 368   //cout<<"ePi0 "<<ePi0<<endl;
 369   //cout<<"Components"<<endl;
 370   p1.SetPx(+ePi0*cosPhi*sinThe);
 371   p1.SetPy(+ePi0*sinPhi*sinThe);
 372   p1.SetPz(+ePi0*cosThe);
 373   p1.SetE(ePi0);
 374   //cout<<"p1: "<<p1.Px()<<" "<<p1.Py()<<" "<<p1.Pz()<<" "<<p1.E()<<endl;
 375   //cout<<"p1 Mass: "<<p1.Px()*p1.Px()+p1.Py()*p1.Py()+p1.Pz()*p1.Pz()-p1.E()*p1.E()<<endl;
 376   p2.SetPx(-ePi0*cosPhi*sinThe);
 377   p2.SetPy(-ePi0*sinPhi*sinThe);
 378   p2.SetPz(-ePi0*cosThe);
 379   p2.SetE(ePi0);
 380   //cout<<"p2: "<<p2.Px()<<" "<<p2.Py()<<" "<<p2.Pz()<<" "<<p2.E()<<endl;
 381   //cout<<"p2 Mass: "<<p2.Px()*p2.Px()+p2.Py()*p2.Py()+p2.Pz()*p2.Pz()-p2.E()*p2.E()<<endl;
 382   //cout<<"Boost "<<b.X()<<" "<<b.Y()<<" "<<b.Z()<<endl;
 383   p1.Boost(b);
 384   //cout<<"p1: "<<p1.Px()<<" "<<p1.Py()<<" "<<p1.Pz()<<" "<<p1.E()<<endl;
 385   p2.Boost(b);
 386   //cout<<"p2: "<<p2.Px()<<" "<<p2.Py()<<" "<<p2.Pz()<<" "<<p2.E()<<endl;
 387   //cout<<"angle"<<endl;
 388   //angle = p1.Angle(p2.Vect());
 389   //cout<<angle<<endl;
 390 }
 391
 392 //________________________________________________________________
 393 void AliGammaMCReader::Print(const Option_t * opt) const
 394 {
 395
 396   //Print some relevant parameters set for the analysis
 397   if(! opt)
 398     return;
 399
 400   Info("Print", "%s %s", GetName(), GetTitle() ) ;
 401
 402   printf("Decay Pi0?          : %d\n", fDecayPi0) ;
 403   printf("Check Overlapping?          : %d\n", fCheckOverlapping) ;
 404
 405 }
 406
 407
 408
 409 //________________________________________________________________
 410 void AliGammaMCReader::SetPhotonStatus(TParticle* pphoton, TParticle* pmother)
 411 {
 412
 413   //Check the origin of the photon and tag it  as decay from pi0, from eta, from other mesons, or prompt or fragmentation.
 414
 415   if(pmother->GetStatusCode() != 21 ){
 416     if(pmother->GetStatusCode() ==11){
 417       if(pmother->GetPdgCode()==111) pphoton->SetStatusCode(kPi0DecayPhoton);//decay gamma from pi0
 418       if(pmother->GetPdgCode()==221) pphoton->SetStatusCode(kEtaDecayPhoton);//decay gamma from eta
 419       else pphoton->SetStatusCode(kOtherDecayPhoton);// decay gamma from other pphotons
 420     }
 421     else         pphoton->SetStatusCode(kUnknown);
 422   }
 423   else if(pmother->GetPdgCode()==22)    pphoton->SetStatusCode(kPromptPhoton);//Prompt photon
 424   else pphoton->SetStatusCode(kFragmentPhoton); //Fragmentation photon
 425
 426 }
 427
 428 //________________________________________________________________
 429 Bool_t AliGammaMCReader::SkipNeutralParticles(Int_t pdg) const {
 430   //Check if pdg is equal to one of the neutral particles list
 431   //These particles will be skipped from analysis.
 432
 433   for(Int_t i= 0; i < fNeutralParticlesArray->GetSize(); i++)
 434     if(TMath::Abs(pdg) ==  fNeutralParticlesArray->At(i)) return kTRUE ;
 435
 436   return kFALSE;
 437
 438 }