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.1.2.1  2007/07/26 10:32:09  schutz
  21  * new analysis classes in the the new analysis framework
  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
  43 //---- ANALYSIS system ----
  44 #include "AliGammaMCReader.h"
  45 #include "Riostream.h"
  46 #include "AliLog.h"
  47
  48 ClassImp(AliGammaMCReader)
  49
  50 //____________________________________________________________________________
  51 AliGammaMCReader::AliGammaMCReader() :
  52   AliGammaReader(),
  53   fEMCALIPDistance(0.),   fPHOSIPDistance(0.),
  54   fEMCALMinDistance(0.),   fPHOSMinDistance(0.),
  55   fDecayPi0(kFALSE)
  56 {
  57   //Ctor
  58
  59   //Initialize parameters
  60   fDataType = kMC;
  61   InitParameters();
  62
  63
  64 }
  65
  66 //____________________________________________________________________________
  67 AliGammaMCReader::AliGammaMCReader(const AliGammaMCReader & g) :
  68   AliGammaReader(g),
  69   fEMCALIPDistance(g.fEMCALIPDistance),  fPHOSIPDistance(g.fPHOSIPDistance),
  70   fEMCALMinDistance(g.fEMCALMinDistance),  fPHOSMinDistance(g.fPHOSMinDistance),
  71   fDecayPi0(g.fDecayPi0)
  72 {
  73   // cpy ctor
  74 }
  75
  76 //_________________________________________________________________________
  77 AliGammaMCReader & AliGammaMCReader::operator = (const AliGammaMCReader & source)
  78 {
  79   // assignment operator
  80
  81   if(&source == this) return *this;
  82
  83   fEMCALIPDistance = source.fEMCALIPDistance;
  84   fPHOSIPDistance = source.fPHOSIPDistance;
  85   fEMCALMinDistance = source.fEMCALMinDistance;
  86   fPHOSMinDistance = source.fPHOSMinDistance;
  87   fDecayPi0 = source.fDecayPi0;
  88
  89
  90   return *this;
  91
  92 }
  93
  94
  95 //____________________________________________________________________________
  96 void AliGammaMCReader::CreateParticleList(TObject * data, TObject *,
  97                                          TClonesArray * plCh,
  98                                          TClonesArray * plEMCAL,
  99                                           TClonesArray * plPHOS,
 100                                           TClonesArray *plParton)
 101 {
 102   //Create list of particles from EMCAL, PHOS and CTS.
 103   AliStack * stack = (AliStack *) data ;
 104   Int_t indexCh     = plCh->GetEntries() ;
 105   Int_t indexEMCAL = plEMCAL->GetEntries() ;
 106   Int_t indexPHOS = plPHOS->GetEntries() ;
 107   Int_t indexParton = plParton->GetEntries() ;
 108   Int_t iParticle = 0 ;
 109   Double_t charge = 0.;
 110
 111   for (iParticle=0 ; iParticle <  stack->GetNprimary() ; iParticle++) {
 112     TParticle * particle = stack->Particle(iParticle);
 113
 114     //Keep partons
 115     if(particle->GetStatusCode() == 21 && iParticle>=2){//All partons, not nucleus
 116       new((*plParton)[indexParton++])  TParticle(*particle) ;
 117     }
 118
 119     //Keep Stable particles
 120     if((particle->GetStatusCode() == 0) && (particle->Pt() > 0)){
 121
 122       charge = TDatabasePDG::Instance()->GetParticle(particle->GetPdgCode())->Charge();
 123
 124       //---------- Charged particles ----------------------
 125       if((charge != 0) && (particle->Pt() > fChargedPtCut)){
 126         //Particles in CTS acceptance
 127         if(TMath::Abs(particle->Eta())<fCTSEtaCut){
 128           //Fill lists
 129           new((*plCh)[indexCh++])       TParticle(*particle) ;
 130         }
 131       }
 132       //-------------Neutral particles ----------------------
 133       else if((charge == 0) && particle->Pt() > fNeutralPtCut &&
 134               TMath::Abs(particle->GetPdgCode())>16){//Avoid neutrinos
 135
 136         if(particle->GetPdgCode()!=111){
 137           if(IsInPHOS(particle->Phi(),particle->Eta()))
 138             new((*plPHOS)[indexPHOS++])  TParticle(*particle) ;
 139           else if(IsInEMCAL(particle->Phi(),particle->Eta()))
 140             new((*plEMCAL)[indexEMCAL++])  TParticle(*particle) ;
 141         }//no pi0
 142         else{
 143           if(fDecayPi0 == kNoDecay){//keep the pi0 do not decay
 144             if(IsInPHOS(particle->Phi(),particle->Eta()))
 145               new((*plPHOS)[indexPHOS++])  TParticle(*particle) ;
 146             else if(IsInEMCAL(particle->Phi(),particle->Eta()))
 147               new((*plEMCAL)[indexEMCAL++])  TParticle(*particle) ;
 148           }
 149           else if(fDecayPi0 == kDecay)
 150             MakePi0Decay(particle,plEMCAL,indexEMCAL,plPHOS, indexPHOS);
 151           else if(fDecayPi0 == kGeantDecay)
 152             SetGeantDecay(particle, stack,plEMCAL, indexEMCAL, plPHOS, indexPHOS);
 153         }//pi0
 154       }//neutral particle
 155     }//stable particle
 156   }//particle loop
 157 }
 158
 159 //___________________________________________________________________________
 160 Bool_t  AliGammaMCReader::IsInEMCAL(Double_t phi, Double_t eta){
 161   //Check if particle is in EMCAL acceptance
 162   if(phi<0)
 163      phi+=TMath::TwoPi();
 164      if( phi > fPhiEMCALCut[0] && phi < fPhiEMCALCut[1] &&
 165        TMath::Abs(eta)<fEMCALEtaCut) return kTRUE ;
 166   else  return kFALSE;
 167
 168   return kFALSE ;
 169 }
 170
 171 //___________________________________________________________________________
 172 Bool_t  AliGammaMCReader::IsInPHOS(Double_t phi, Double_t eta){
 173   //Check if particle is in EMCAL acceptance
 174   if(phi<0)
 175     phi+=TMath::TwoPi();
 176   if( phi > fPhiPHOSCut[0] && phi < fPhiPHOSCut[1] &&
 177       TMath::Abs(eta)<fPHOSEtaCut) return kTRUE ;
 178   else  return kFALSE;
 179
 180   return kFALSE ;
 181 }
 182
 183 //___________________________________________________________________________
 184 void AliGammaMCReader::SetGeantDecay(TParticle * particle, AliStack * stack,
 185                                      TClonesArray * plEMCAL, Int_t &indexEMCAL,
 186                                      TClonesArray * plPHOS, Int_t &indexPHOS){
 187   //Find decay gamma from pi0 and put them in the list.
 188
 189   Int_t ndaug = particle->GetNDaughters() ;
 190   if(ndaug<=2 && ndaug >0){//At least 1 daugther
 191     TParticle * d1 = stack->Particle(particle->GetDaughter(0));
 192     if(d1->GetPdgCode()==22){
 193       if(IsInEMCAL(d1->Phi(),d1->Eta()))
 194         new((*plEMCAL)[indexEMCAL++])       TParticle(*d1) ;
 195       else if(IsInPHOS(d1->Phi(),d1->Eta()))
 196         new((*plPHOS)[indexPHOS++])       TParticle(*d1) ;
 197     }
 198
 199     if(ndaug>1){//second daugther if present
 200       TParticle * d2 = stack->Particle(particle->GetDaughter(1));
 201       if(IsInEMCAL(d2->Phi(),d2->Eta()))
 202         new((*plEMCAL)[indexEMCAL++])       TParticle(*d2) ;
 203       else if(IsInPHOS(d2->Phi(),d2->Eta()))
 204         new((*plPHOS)[indexPHOS++])       TParticle(*d2) ;
 205     }
 206   }
 207 }
 208
 209 //___________________________________________________________________________
 210 void AliGammaMCReader::MakePi0Decay(TParticle * particle,
 211                                     TClonesArray * plEMCAL, Int_t &indexEMCAL,
 212                                     TClonesArray * plPHOS, Int_t &indexPHOS){
 213
 214   //Decays pi0, see if aperture angle is small and then add the pi0 or the 2 gamma
 215
 216   TLorentzVector pPi0, pGamma1, pGamma2 ;
 217   Double_t angle = 0, cellDistance = 0.;
 218   Bool_t checkPhoton = kTRUE;
 219
 220   pPi0.SetPxPyPzE(particle->Px(),particle->Py(),particle->Pz(),particle->Energy());
 221
 222   //Decay
 223   Pi0Decay(pPi0,pGamma1,pGamma2,angle);
 224
 225   //Check if Pi0 is in the acceptance of the calorimeters, if aperture angle is small, keep it
 226   if(IsInPHOS(particle->Phi(), particle->Eta())){
 227     cellDistance = angle*fPHOSIPDistance;
 228     if (cellDistance < fPHOSMinDistance){
 229       new((*plPHOS)[indexPHOS++])       TParticle(*particle) ;
 230       checkPhoton = kFALSE;
 231     }
 232   }
 233   else if(IsInEMCAL(particle->Phi(), particle->Eta())){
 234     cellDistance = angle*fEMCALIPDistance;
 235     if (cellDistance < fEMCALMinDistance) {
 236       new((*plEMCAL)[indexEMCAL++])       TParticle(*particle) ;
 237       checkPhoton = kFALSE;
 238     }
 239   }
 240   else checkPhoton = kTRUE ;
 241
 242   if (checkPhoton) {
 243     //Gamma Not overlapped
 244     TParticle * photon1 = new TParticle(22,1,0,0,0,0,pGamma1.Px(),pGamma1.Py(),
 245                                         pGamma1.Pz(),pGamma1.E(),0,0,0,0);
 246     if(photon1->Pt() > fNeutralPtCut) {
 247
 248       if(IsInPHOS(photon1->Phi(), photon1->Eta()))
 249         new((*plPHOS)[indexPHOS++])       TParticle(*photon1) ;
 250
 251       else if(IsInEMCAL(photon1->Phi(), photon1->Eta()))
 252         new((*plEMCAL)[indexEMCAL++])       TParticle(*photon1) ;
 253
 254     }// photon 1 of pi0 in acceptance
 255
 256
 257     TParticle * photon2 = new TParticle(22,1,0,0,0,0,pGamma2.Px(),pGamma2.Py(),
 258                                         pGamma2.Pz(),pGamma2.E(),0,0,0,0);
 259
 260     if(photon2->Pt() > fNeutralPtCut) {
 261
 262       if(IsInPHOS(photon2->Phi(), photon2->Eta()))
 263         new((*plPHOS)[indexPHOS++])       TParticle(*photon2) ;
 264
 265       else if(IsInEMCAL(photon2->Phi(), photon2->Eta()))
 266         new((*plEMCAL)[indexEMCAL++])       TParticle(*photon2) ;
 267
 268     }// photon 2 of pi0 in acceptance
 269   }//Not overlapped gamma
 270 }
 271
 272
 273 //_______________________________________________________________
 274 void AliGammaMCReader::InitParameters()
 275 {
 276
 277   //Initialize the parameters of the analysis.
 278
 279   //Fill particle lists when PID is ok
 280   fEMCALMinDistance    = 10. ;
 281   fPHOSMinDistance    = 3.6 ;
 282   fEMCALIPDistance    = 450. ;//cm
 283   fPHOSIPDistance    = 460. ;//cm
 284   fDecayPi0 = kGeantDecay;
 285 }
 286
 287 //____________________________________________________________________________
 288 void AliGammaMCReader::Pi0Decay(TLorentzVector &p0, TLorentzVector &p1,
 289                                 TLorentzVector &p2, Double_t &angle) {
 290   // Perform isotropic decay pi0 -> 2 photons
 291   // p0 is pi0 4-momentum (inut)
 292   // p1 and p2 are photon 4-momenta (output)
 293   //  cout<<"Boost vector"<<endl;
 294   Double_t mPi0 = TDatabasePDG::Instance()->GetParticle(111)->Mass();
 295   TVector3 b = p0.BoostVector();
 296   //cout<<"Parameters"<<endl;
 297   //Double_t mPi0   = p0.M();
 298   Double_t phi    = TMath::TwoPi() * gRandom->Rndm();
 299   Double_t cosThe = 2 * gRandom->Rndm() - 1;
 300   Double_t cosPhi = TMath::Cos(phi);
 301   Double_t sinPhi = TMath::Sin(phi);
 302   Double_t sinThe = TMath::Sqrt(1-cosThe*cosThe);
 303   Double_t ePi0   = mPi0/2.;
 304   //cout<<"ePi0 "<<ePi0<<endl;
 305   //cout<<"Components"<<endl;
 306   p1.SetPx(+ePi0*cosPhi*sinThe);
 307   p1.SetPy(+ePi0*sinPhi*sinThe);
 308   p1.SetPz(+ePi0*cosThe);
 309   p1.SetE(ePi0);
 310   //cout<<"p1: "<<p1.Px()<<" "<<p1.Py()<<" "<<p1.Pz()<<" "<<p1.E()<<endl;
 311   //cout<<"p1 Mass: "<<p1.Px()*p1.Px()+p1.Py()*p1.Py()+p1.Pz()*p1.Pz()-p1.E()*p1.E()<<endl;
 312   p2.SetPx(-ePi0*cosPhi*sinThe);
 313   p2.SetPy(-ePi0*sinPhi*sinThe);
 314   p2.SetPz(-ePi0*cosThe);
 315   p2.SetE(ePi0);
 316   //cout<<"p2: "<<p2.Px()<<" "<<p2.Py()<<" "<<p2.Pz()<<" "<<p2.E()<<endl;
 317   //cout<<"p2 Mass: "<<p2.Px()*p2.Px()+p2.Py()*p2.Py()+p2.Pz()*p2.Pz()-p2.E()*p2.E()<<endl;
 318   //cout<<"Boost "<<b.X()<<" "<<b.Y()<<" "<<b.Z()<<endl;
 319   p1.Boost(b);
 320   //cout<<"p1: "<<p1.Px()<<" "<<p1.Py()<<" "<<p1.Pz()<<" "<<p1.E()<<endl;
 321   p2.Boost(b);
 322   //cout<<"p2: "<<p2.Px()<<" "<<p2.Py()<<" "<<p2.Pz()<<" "<<p2.E()<<endl;
 323   //cout<<"angle"<<endl;
 324   angle = p1.Angle(p2.Vect());
 325   //cout<<angle<<endl;
 326 }
 327
 328 //________________________________________________________________
 329 void AliGammaMCReader::Print(const Option_t * opt) const
 330 {
 331
 332   //Print some relevant parameters set for the analysis
 333   if(! opt)
 334     return;
 335
 336   Info("Print", "%s %s", GetName(), GetTitle() ) ;
 337
 338   printf("IP distance to PHOS         : %f\n", fPHOSIPDistance) ;
 339   printf("IP distance to EMCAL         : %f\n", fEMCALIPDistance) ;
 340   printf("Min gamma decay distance in PHOS         : %f\n", fPHOSMinDistance) ;
 341   printf("Min gamma decay distance in EMCAL         : %f\n", fEMCALMinDistance) ;
 342   printf("Decay Pi0?          : %d\n", fDecayPi0) ;
 343
 344 }