1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 /* History of cvs commits:
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)
23 * Revision 1.4 2007/10/29 13:48:42 gustavo
24 * Corrected coding violations
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)
34 //*-- Author: Gustavo Conesa (LNF-INFN)
35 //////////////////////////////////////////////////////////////////////////////
38 // --- ROOT system ---
41 #include <TParticle.h>
46 #include <TClonesArray.h>
48 //---- ANALYSIS system ----
49 #include "AliGammaMCReader.h"
50 #include "Riostream.h"
54 ClassImp(AliGammaMCReader)
56 //____________________________________________________________________________
57 AliGammaMCReader::AliGammaMCReader() :
58 AliGammaReader(), fDecayPi0(0), fCheckOverlapping(kFALSE), fNeutralParticlesArray(0x0)
62 //Initialize parameters
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)
76 //_________________________________________________________________________
77 AliGammaMCReader & AliGammaMCReader::operator = (const AliGammaMCReader & source)
79 // assignment operator
81 if(&source == this) return *this;
83 fDecayPi0 = source.fDecayPi0;
84 fCheckOverlapping = source.fCheckOverlapping;
85 delete fNeutralParticlesArray;
86 fNeutralParticlesArray = source.fNeutralParticlesArray?new TArrayI(*source.fNeutralParticlesArray):0x0;
91 //_________________________________
92 AliGammaMCReader::~AliGammaMCReader() {
95 delete fNeutralParticlesArray ;
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.
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);
113 if((index == idaug0 && pdaug0->Pt() > fNeutralPtCut) ||
114 (index == idaug1 && pdaug0->Pt() <= fNeutralPtCut)){//Check decay when first daughter arrives, do nothing with second.
116 FillListWithDecayGammaOrPi0(pmother, pdaug0, pdaug1, plEMCAL, indexEMCAL, plPHOS, indexPHOS);
119 }//mother is a pi0 with 2 daughters
122 if(IsInPHOS(particle->Phi(),particle->Eta())) {
123 TParticle * pmother =sta->Particle(particle->GetFirstMother());
124 SetPhotonStatus(particle,pmother);
125 new((*plPHOS)[indexPHOS++]) TParticle(*particle) ;
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) ;
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.
141 Int_t ndaug = particle->GetNDaughters() ;
142 TParticle * d1 = new TParticle();
143 TParticle * d2 = new TParticle();
145 if(ndaug > 0){//At least 1 daugther
146 d1 = stack->Particle(particle->GetDaughter(0));
148 d2 = stack->Particle(particle->GetDaughter(1));
150 FillListWithDecayGammaOrPi0(particle, d1, d2, plEMCAL, indexEMCAL, plPHOS, indexPHOS);
155 //___________________________________________________________________________
156 void AliGammaMCReader::CasePi0Decay(TParticle * pPi0,
157 TClonesArray * plEMCAL, Int_t &indexEMCAL,
158 TClonesArray * plPHOS, Int_t &indexPHOS){
160 //Decays pi0, see if aperture angle is small and then add the pi0 or the 2 gamma
162 TLorentzVector lvPi0, lvGamma1, lvGamma2 ;
163 //Double_t angle = 0;
165 lvPi0.SetPxPyPzE(pPi0->Px(),pPi0->Py(),pPi0->Pz(),pPi0->Energy());
168 MakePi0Decay(lvPi0,lvGamma1,lvGamma2);//,angle);
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);
176 FillListWithDecayGammaOrPi0(pPi0, pPhoton1, pPhoton2, plEMCAL, indexEMCAL, plPHOS, indexPHOS);
180 //_______________________________________________________________
181 void AliGammaMCReader::InitParameters()
184 //Initialize the parameters of the analysis.
186 fDecayPi0 = kGeantDecay;
187 fCheckOverlapping = kTRUE ;
189 Int_t pdgarray[]={12,14,16};// skip neutrinos
190 fNeutralParticlesArray = new TArrayI(3, pdgarray);
194 //____________________________________________________________________________
195 void AliGammaMCReader::CreateParticleList(TObject * data, TObject *,
197 TClonesArray * plEMCAL,
198 TClonesArray * plPHOS,
199 TClonesArray *plParton,TClonesArray *,TClonesArray *)
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.;
210 for (iParticle=0 ; iParticle < stack->GetNprimary() ; iParticle++) {
211 TParticle * particle = stack->Particle(iParticle);
214 if(particle->GetStatusCode() == 21 && iParticle>=2){//All partons, not nucleus
215 new((*plParton)[indexParton++]) TParticle(*particle) ;
218 //Keep Stable particles
219 if((particle->GetStatusCode() == 1) && (particle->Pt() > 0)){
221 charge = TDatabasePDG::Instance()->GetParticle(particle->GetPdgCode())->Charge();
223 //---------- Charged particles ----------------------
224 if( fSwitchOnCTS && (charge != 0) && (particle->Pt() > fChargedPtCut)){
225 //Particles in CTS acceptance
226 if(TMath::Abs(particle->Eta())<fCTSEtaCut){
228 new((*plCh)[indexCh++]) TParticle(*particle) ;
232 //-------------Neutral particles ----------------------
233 else if((charge == 0) && particle->Pt() > fNeutralPtCut ){ //&&
234 //TMath::Abs(particle->GetPdgCode())>16){//Avoid neutrinos
236 if(particle->GetPdgCode()!=111){
237 //In case that we force PYTHIA to decay pi0, and we want to check the overlapping of
239 if(particle->GetPdgCode() == 22 && fDecayPi0==kDecayGamma){
240 CaseDecayGamma(iParticle,particle, stack,plEMCAL, indexEMCAL, plPHOS, indexPHOS); //If pythia decays pi0
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) ;
250 else if(IsInEMCAL(particle->Phi(),particle->Eta())){
251 if(particle->GetPdgCode()==22) SetPhotonStatus(particle,pmother);
252 new((*plEMCAL)[indexEMCAL++]) TParticle(*particle) ;
254 }//skip neutral particles
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) ;
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);
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) ;
281 //___________________________________________________________________________
282 void AliGammaMCReader::FillListWithDecayGammaOrPi0(TParticle * pPi0, TParticle * pdaug0, TParticle * pdaug1,
283 TClonesArray * plEMCAL, Int_t &indexEMCAL,
284 TClonesArray * plPHOS, Int_t &indexPHOS){
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);
292 if(fCheckOverlapping){
293 Double_t angle = lv1.Angle(lv2.Vect());
295 if(IsInEMCAL(pPi0->Phi(), pPi0->Eta()))
297 if (angle < fEMCALMinAngle){
298 pPi0->SetStatusCode(1);
299 new((*plEMCAL)[indexEMCAL++]) TParticle(*pPi0) ;
304 else if(IsInPHOS(pPi0->Phi(), pPi0->Eta())){
305 if (angle < fPHOSMinAngle){
306 pPi0->SetStatusCode(1);
307 new((*plPHOS)[indexPHOS++]) TParticle(*pPi0) ;
314 //Fill with gammas if not overlapp
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) ;
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) ;
337 //___________________________________________________________________________
338 Bool_t AliGammaMCReader::IsInEMCAL(Double_t phi, Double_t eta) const {
339 //Check if particle is in EMCAL acceptance
344 if( phi > fPhiEMCALCut[0] && phi < fPhiEMCALCut[1] &&
345 TMath::Abs(eta)<fEMCALEtaCut) return kTRUE ;
351 //___________________________________________________________________________
352 Bool_t AliGammaMCReader::IsInPHOS(Double_t phi, Double_t eta) const {
353 //Check if particle is in PHOS acceptance
358 if( phi > fPhiPHOSCut[0] && phi < fPhiPHOSCut[1] &&
359 TMath::Abs(eta)<fPHOSEtaCut) return kTRUE ;
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);
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);
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;
399 //cout<<"p1: "<<p1.Px()<<" "<<p1.Py()<<" "<<p1.Pz()<<" "<<p1.E()<<endl;
401 //cout<<"p2: "<<p2.Px()<<" "<<p2.Py()<<" "<<p2.Pz()<<" "<<p2.E()<<endl;
402 //cout<<"angle"<<endl;
403 //angle = p1.Angle(p2.Vect());
407 //________________________________________________________________
408 void AliGammaMCReader::Print(const Option_t * opt) const
411 //Print some relevant parameters set for the analysis
415 Info("Print", "%s %s", GetName(), GetTitle() ) ;
417 printf("Decay Pi0? : %d\n", fDecayPi0) ;
418 printf("Check Overlapping? : %d\n", fCheckOverlapping) ;
424 //________________________________________________________________
425 void AliGammaMCReader::SetPhotonStatus(TParticle* pphoton, TParticle* pmother)
428 //Check the origin of the photon and tag it as decay from pi0, from eta, from other mesons, or prompt or fragmentation.
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
436 else pphoton->SetStatusCode(kUnknown);
438 else if(pmother->GetPdgCode()==22) pphoton->SetStatusCode(kPromptPhoton);//Prompt photon
439 else pphoton->SetStatusCode(kFragmentPhoton); //Fragmentation photon
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.
448 for(Int_t i= 0; i < fNeutralParticlesArray->GetSize(); i++)
449 if(TMath::Abs(pdg) == fNeutralParticlesArray->At(i)) return kTRUE ;