Introducing new class : AliPHOSGammaJet
[u/mrichter/AliRoot.git] / PHOS / AliPHOSGammaJet.cxx
CommitLineData
77414c90 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//_________________________________________________________________________
18// Class for the analysis of gamma-jet correlations
19//
20//
21//*-- Author: Gustavo Conesa & Yves Schutz (IFIC, CERN)
22//////////////////////////////////////////////////////////////////////////////
23
24
25// --- ROOT system ---
26
27#include "TString.h"
28#include "TFile.h"
29#include "TLorentzVector.h"
30#include "TList.h"
31#include "TTree.h"
32#include "TParticle.h"
33#include "AliRun.h"
34#include "AliPHOSGammaJet.h"
35#include "AliPHOSGetter.h"
36
37ClassImp(AliPHOSGammaJet)
38
39//____________________________________________________________________________
40AliPHOSGammaJet::AliPHOSGammaJet() {
41 // ctor
42}
43
44//____________________________________________________________________________
45AliPHOSGammaJet::AliPHOSGammaJet(const TString inputfilename) {
46 // ctor
47 AliPHOSGetter::Instance(inputfilename) ;
48}
49
50//____________________________________________________________________________
51AliPHOSGammaJet::AliPHOSGammaJet(const AliPHOSGammaJet * gj) {
52}
53
54//____________________________________________________________________________
55AliPHOSGammaJet::~AliPHOSGammaJet() {
56}
57
58//____________________________________________________________________________
59void AliPHOSGammaJet::Exec(Option_t *option)
60{
61 // does the job
62
63 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
64 Int_t nEvent1 = gime->MaxEvent() ;
65 Int_t iEvent = 0 ;
66 for ( iEvent = 0 ; iEvent < nEvent1 ; iEvent++) {
67 if (iEvent <= 100 || iEvent%10 == 0)
68 Info("Exec", "Event %d", iEvent) ;
69 TParticle * particle = 0 ;
70 //-----------------Fill list with particles--------------------
71 TLorentzVector pPi0, pGamma1, pGamma2 ;
72 Double_t angle = 0 ;
73 TList particleList ;
74 Int_t n = -1;
75 gime->Event(iEvent, "X") ;
76 Int_t nparticles = gime->NPrimaries() ;
77 Int_t iParticle=0 ;
78 for (iParticle=0 ; iParticle < nparticles ; iParticle++) {
79 //Keep original partons
80 particle = gime->Primary(iParticle) ;
81 if((particle->GetStatusCode()== 21)){
82 particleList.Add(particle);
83 n++;
84 }
85 //Keep Stable particles within eta range
86 Float_t etacut = 0. ;
87 if((particle->GetStatusCode() == 1)&&
88 (TMath::Abs(particle->Eta())<etacut)){
89 // Keep particles different from Pi0
90 if(particle->GetPdgCode() != 111){
91 particleList.Add(particle);
92 n++;
93 }
94 //Decay Pi0 and keep it with different status name
95 //Keep decay photons
96 if(particle->GetPdgCode() == 111) {
97 //cout<<"Pi0 "<<endl;
98 n += 3 ;
99 particle->Momentum(pPi0);
100 Pi0Decay(particle->GetMass(),pPi0,pGamma1,pGamma2,angle);
101 TParticle * photon1 = new TParticle(22,1,0,0,0,0,pGamma1.Px(),pGamma1.Py(),
102 pGamma1.Pz(),pGamma1.E(),0,0,0,0);
103 TParticle * photon2 = new TParticle(22,1,0,0,0,0,pGamma2.Px(),pGamma2.Py(),
104 pGamma2.Pz(),pGamma2.E(),0,0,0,0);
105 photon1->SetWeight(1);
106 photon2->SetWeight(2);
107 particle->SetStatusCode(2);
108 particleList.Add(particle);
109 particleList.Add(photon1);
110 particleList.Add(photon2);
111 //photon1->Print();
112 //photon2->Print();
113 }//if pi0
114 }//final particle etacut
115 }//for (iParticle<nParticle)
116 TLorentzVector gamma,charge,pi0, gammapair;
117 Int_t id_g = -1;
118 GetGammaJet(particleList,gamma, id_g);
119 GetLeadingCharge(particleList,charge, id_g);
120 GetLeadingPi0(particleList,pi0, id_g);
121 Info("Pi0Decay", "Gamma: %f %d", gamma.Energy(), id_g) ;
122 Info("Pi0Decay", "Charge: %f", charge.Energy()) ;
123 Info("Pi0Decay", "Pi0: %f", pi0.Energy()) ;
124 // GetLeadingGammaPair(particleList, gammapair, id_g,
125 // 0.04,0.2, 1.0,0.13,0.14);
126 Info("Pi0Decay", "Pair: %f", gammapair.Energy()) ;
127 }//loop: events
128}
129
130//____________________________________________________________________________
131void AliPHOSGammaJet::Pi0Decay(Double_t mPi0, TLorentzVector &p0,
132 TLorentzVector &p1, TLorentzVector &p2, Double_t &angle) {
133 // Perform isotropic decay pi0 -> 2 photons
134 // p0 is pi0 4-momentum (input)
135 // p1 and p2 are photon 4-momenta (output)
136 // cout<<"Boost vector"<<endl;
137 TVector3 b = p0.BoostVector();
138 //cout<<"Parameters"<<endl;
139 //Double_t mPi0 = p0.M();
140 Double_t phi = TMath::TwoPi() * gRandom->Rndm();
141 Double_t cosThe = 2 * gRandom->Rndm() - 1;
142 Double_t cosPhi = TMath::Cos(phi);
143 Double_t sinPhi = TMath::Sin(phi);
144 Double_t sinThe = TMath::Sqrt(1-cosThe*cosThe);
145 Double_t ePi0 = mPi0/2.;
146 //cout<<"ePi0 "<<ePi0<<endl;
147 //cout<<"Components"<<endl;
148 p1.SetPx(+ePi0*cosPhi*sinThe);
149 p1.SetPy(+ePi0*sinPhi*sinThe);
150 p1.SetPz(+ePi0*cosThe);
151 p1.SetE(ePi0);
152 //cout<<"p1: "<<p1.Px()<<" "<<p1.Py()<<" "<<p1.Pz()<<" "<<p1.E()<<endl;
153 //cout<<"p1 Mass: "<<p1.Px()*p1.Px()+p1.Py()*p1.Py()+p1.Pz()*p1.Pz()-p1.E()*p1.E()<<endl;
154 p2.SetPx(-ePi0*cosPhi*sinThe);
155 p2.SetPy(-ePi0*sinPhi*sinThe);
156 p2.SetPz(-ePi0*cosThe);
157 p2.SetE(ePi0);
158 //cout<<"p2: "<<p2.Px()<<" "<<p2.Py()<<" "<<p2.Pz()<<" "<<p2.E()<<endl;
159 //cout<<"p2 Mass: "<<p2.Px()*p2.Px()+p2.Py()*p2.Py()+p2.Pz()*p2.Pz()-p2.E()*p2.E()<<endl;
160 //cout<<"Boost "<<b.X()<<" "<<b.Y()<<" "<<b.Z()<<endl;
161 p1.Boost(b);
162 //cout<<"p1: "<<p1.Px()<<" "<<p1.Py()<<" "<<p1.Pz()<<" "<<p1.E()<<endl;
163 p2.Boost(b);
164 //cout<<"p2: "<<p2.Px()<<" "<<p2.Py()<<" "<<p2.Pz()<<" "<<p2.E()<<endl;
165 //cout<<"angle"<<endl;
166 angle = p1.Angle(p2.Vect());
167 //cout<<angle<<endl;
168}
169//____________________________________________________________________________
170void AliPHOSGammaJet::GetGammaJet(TList &particleList, TLorentzVector &gamma, Int_t & id)
171{
172 TParticle *particle = 0x0;
173
174 Int_t iPrimary=-1, id_motherg = -1;
175
176 TIter next(&particleList);
177 while ( (particle = (TParticle*)next()) ) {
178 iPrimary++;
179 Int_t ksCode = particle->GetStatusCode();
180 Int_t iMother= particle->GetMother(0);
181
182 if (ksCode == 21 && iMother == -1)
183 if(particle->GetPdgCode()==22){
184 id_motherg = iPrimary;
185 //cout<<"idmother "<<id_motherg<<endl;
186 }
187 if(ksCode == 1){
188
189 if(id_motherg == particle->GetMother(0)){
190 particle->Momentum(gamma);
191 id = iPrimary;
192 break;
193 }
194 }// kscode == 1
195 }// while
196}
197
198//____________________________________________________________________________
199void AliPHOSGammaJet::GetLeadingCharge(TList &particleList, TLorentzVector &charge, Int_t & id)
200{
201 TParticle *particle = 0x0;
202
203 Int_t iPrimary=-1;
204 Double_t pt_max = 0, pt_i = 0;
205 TIter next(&particleList);
206 while ( (particle = (TParticle*)next()) ) {
207 iPrimary++;
208 Int_t ksCode = particle->GetStatusCode();
209
210 if((ksCode == 1)&&(id != iPrimary)
211 &&(particle->GetPDG(0)->Charge()!=0)){
212 pt_i = particle->Pt();
213 if(pt_i> pt_max){
214 pt_max = pt_i;
215 particle->Momentum(charge);
216 }//pt_max
217 }// kscode == 1
218 }// while
219}
220
221//____________________________________________________________________________
222void AliPHOSGammaJet::GetLeadingPi0(TList &particleList, TLorentzVector &pi0, Int_t & id)
223{
224 TParticle *particle = 0x0;
225
226 Int_t iPrimary=-1;
227 Double_t pt_max = 0, pt_i = 0;
228 TIter next(&particleList);
229 while ( (particle = (TParticle*)next()) ) {
230 iPrimary++;
231 Int_t ksCode = particle->GetStatusCode();
232
233 if((ksCode == 2))
234 //&&(id != iPrimary)&&(particle->GetPdgCode() == 111))
235 {
236 pt_i = particle->Pt();
237 if(pt_i> pt_max){
238 pt_max = pt_i;
239 particle->Momentum(pi0);
240 }//pt_max
241 }// kscode == 1
242 }// while
243}
244
245//____________________________________________________________________________
246void AliPHOSGammaJet::GetLeadingGammaPair(TList &particleList, TLorentzVector &gammapair, Int_t & id,
247 Double_t & thetacut,Double_t & ratiocut1, Double_t & ratiocut2,
248 Double_t & invmasscut1,Double_t & invmasscut2)
249{
250// TParticle *particle = 0x0;
251
252// Int_t iPrimary=-1;
253// Double_t El = 0, E12 = 0;
254// TLorentzVector gamma_i,gamma_j;
255// TIter next(&particleList);
256// while ( (particle = (TParticle*)next()) ) {
257// iPrimary++;
258// Int_t ksCode = particle->GetStatusCode();
259// Int_t ksPdg = particle->GetPdgCode();
260// if((ksCode == 1) && (iPrimary != id) && (ksPdg == 22)){
261// particle->Momentum(gamma_i);
262// Int_t jPrimary=-1;
263// TIter next2(&particleList);
264// while ( (particle = (TParticle*)next2()) ) {
265// jPrimary++;
266// if(jPrimary>iPrimary){
267// ksCode = particle->GetStatusCode();
268// ksPdg = particle->GetPdgCode();
269// if((ksCode == 1) && (iPrimary != id) && (ksPdg == 22)){
270// particle->Momentum(gamma_j);
271// if(gamma_j.Angle(gamma_i.Vect())<thetacut){
272// Float_t invmass = (gamma_i+gamma_j).M();
273// h_invmass->Fill(Eg,invmass);
274// if((invmass>invmasscut1) && (invmass<invmasscut2)){
275// E12 = (gamma_i+gamma_j).Energy();
276// if(E12>El && (E12/Eg>ratiocut1) && (E12/Eg<ratiocut2)){
277// //cout<<E12<<" "<<E12/Eg<<endl;
278// El = E12;
279// id_i = iPrimary;
280// id_j = jPrimary;
281// gammapair = gamma_i+gamma_j;
282// }//E12>El && (E12/Eg>0.2 && E12/Eg<1.)
283// }//(invmass>0.125) && (invmass<0.145)
284// }//gamma_j.Angle(gamma_i.Vect())<0.04
285// }//(ksCode == 1)
286// }
287// }//while
288// // cout<<"jPrimary "<<jPrimary<<endl;
289// }// if kscode 1
290// }//while
291// // cout<<"iPrimary "<<iPrimary<<endl;
292}