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