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