]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PWG4/AliGammaMCReader.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Correctly taking into account multiple coulomb scattering effects
[u/mrichter/AliRoot.git] / PWG4 / AliGammaMCReader.cxx
... / ...
CommitLineData
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
51ClassImp(AliGammaMCReader)
52
53//____________________________________________________________________________
54AliGammaMCReader::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//____________________________________________________________________________
66AliGammaMCReader::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//_________________________________________________________________________
74AliGammaMCReader & 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//_________________________________
89AliGammaMCReader::~AliGammaMCReader() {
90 //Dtor
91
92 delete fNeutralParticlesArray ;
93
94}
95
96//___________________________________________________________________________
97void 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//___________________________________________________________________________
153void 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//_______________________________________________________________
178void 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//____________________________________________________________________________
192void 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//___________________________________________________________________________
273void 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//___________________________________________________________________________
323Bool_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//___________________________________________________________________________
337Bool_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//____________________________________________________________________________
352void 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//________________________________________________________________
393void 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//________________________________________________________________
410void 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//________________________________________________________________
429Bool_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}