]>
Commit | Line | Data |
---|---|---|
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 --- | |
e957fea8 | 26 | #include "TRandom.h" |
77414c90 | 27 | |
77414c90 | 28 | #include "TLorentzVector.h" |
29 | #include "TList.h" | |
77414c90 | 30 | #include "TParticle.h" |
77414c90 | 31 | #include "AliPHOSGammaJet.h" |
32 | #include "AliPHOSGetter.h" | |
33 | ||
34 | ClassImp(AliPHOSGammaJet) | |
35 | ||
36 | //____________________________________________________________________________ | |
37 | AliPHOSGammaJet::AliPHOSGammaJet() { | |
38 | // ctor | |
39 | } | |
40 | ||
41 | //____________________________________________________________________________ | |
42 | AliPHOSGammaJet::AliPHOSGammaJet(const TString inputfilename) { | |
43 | // ctor | |
44 | AliPHOSGetter::Instance(inputfilename) ; | |
45 | } | |
46 | ||
47 | //____________________________________________________________________________ | |
90cceaf6 | 48 | AliPHOSGammaJet::AliPHOSGammaJet(const AliPHOSGammaJet & gj) : TTask(gj) { |
77414c90 | 49 | } |
50 | ||
51 | //____________________________________________________________________________ | |
52 | AliPHOSGammaJet::~AliPHOSGammaJet() { | |
53 | } | |
54 | ||
55 | //____________________________________________________________________________ | |
e957fea8 | 56 | void AliPHOSGammaJet::Exec(Option_t * opt) |
77414c90 | 57 | { |
58 | // does the job | |
59 | ||
e957fea8 | 60 | if (! opt) |
61 | return ; | |
62 | ||
77414c90 | 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; | |
e957fea8 | 117 | Int_t idg = -1; |
118 | GetGammaJet(particleList,gamma, idg); | |
119 | GetLeadingCharge(particleList,charge, idg); | |
90cceaf6 | 120 | GetLeadingPi0(particleList,pi0); |
e957fea8 | 121 | Info("Pi0Decay", "Gamma: %f %d", gamma.Energy(), idg) ; |
77414c90 | 122 | Info("Pi0Decay", "Charge: %f", charge.Energy()) ; |
123 | Info("Pi0Decay", "Pi0: %f", pi0.Energy()) ; | |
e957fea8 | 124 | // GetLeadingGammaPair(particleList, gammapair, idg, |
77414c90 | 125 | // 0.04,0.2, 1.0,0.13,0.14); |
126 | Info("Pi0Decay", "Pair: %f", gammapair.Energy()) ; | |
127 | }//loop: events | |
128 | } | |
129 | ||
130 | //____________________________________________________________________________ | |
131 | void 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 | //____________________________________________________________________________ | |
170 | void AliPHOSGammaJet::GetGammaJet(TList &particleList, TLorentzVector &gamma, Int_t & id) | |
171 | { | |
e957fea8 | 172 | // Get the lists of jet particles and gamma |
77414c90 | 173 | TParticle *particle = 0x0; |
174 | ||
175 | Int_t iPrimary=-1, id_motherg = -1; | |
176 | ||
177 | TIter next(&particleList); | |
178 | while ( (particle = (TParticle*)next()) ) { | |
179 | iPrimary++; | |
180 | Int_t ksCode = particle->GetStatusCode(); | |
181 | Int_t iMother= particle->GetMother(0); | |
182 | ||
183 | if (ksCode == 21 && iMother == -1) | |
184 | if(particle->GetPdgCode()==22){ | |
185 | id_motherg = iPrimary; | |
186 | //cout<<"idmother "<<id_motherg<<endl; | |
187 | } | |
188 | if(ksCode == 1){ | |
189 | ||
190 | if(id_motherg == particle->GetMother(0)){ | |
191 | particle->Momentum(gamma); | |
192 | id = iPrimary; | |
193 | break; | |
194 | } | |
195 | }// kscode == 1 | |
196 | }// while | |
197 | } | |
198 | ||
199 | //____________________________________________________________________________ | |
200 | void AliPHOSGammaJet::GetLeadingCharge(TList &particleList, TLorentzVector &charge, Int_t & id) | |
201 | { | |
e957fea8 | 202 | // Gets the leading particle from the list of charged particles |
77414c90 | 203 | TParticle *particle = 0x0; |
204 | ||
205 | Int_t iPrimary=-1; | |
e957fea8 | 206 | Double_t ptmax = 0, pti = 0; |
77414c90 | 207 | TIter next(&particleList); |
208 | while ( (particle = (TParticle*)next()) ) { | |
209 | iPrimary++; | |
210 | Int_t ksCode = particle->GetStatusCode(); | |
211 | ||
212 | if((ksCode == 1)&&(id != iPrimary) | |
213 | &&(particle->GetPDG(0)->Charge()!=0)){ | |
e957fea8 | 214 | pti = particle->Pt(); |
215 | if(pti> ptmax){ | |
216 | ptmax = pti; | |
77414c90 | 217 | particle->Momentum(charge); |
e957fea8 | 218 | }//ptmax |
77414c90 | 219 | }// kscode == 1 |
220 | }// while | |
221 | } | |
222 | ||
223 | //____________________________________________________________________________ | |
90cceaf6 | 224 | void AliPHOSGammaJet::GetLeadingPi0(TList &particleList, TLorentzVector &pi0) |
77414c90 | 225 | { |
e957fea8 | 226 | // Gets the leading pi0 from the list of particles |
77414c90 | 227 | TParticle *particle = 0x0; |
228 | ||
229 | Int_t iPrimary=-1; | |
e957fea8 | 230 | Double_t ptmax = 0, pti = 0; |
77414c90 | 231 | TIter next(&particleList); |
232 | while ( (particle = (TParticle*)next()) ) { | |
233 | iPrimary++; | |
234 | Int_t ksCode = particle->GetStatusCode(); | |
235 | ||
236 | if((ksCode == 2)) | |
77414c90 | 237 | { |
e957fea8 | 238 | pti = particle->Pt(); |
239 | if(pti> ptmax){ | |
240 | ptmax = pti; | |
77414c90 | 241 | particle->Momentum(pi0); |
e957fea8 | 242 | }//ptmax |
77414c90 | 243 | }// kscode == 1 |
244 | }// while | |
245 | } | |
246 | ||
247 | //____________________________________________________________________________ | |
90cceaf6 | 248 | // void AliPHOSGammaJet::GetLeadingGammaPair(TList &particleList, TLorentzVector &gammapair, Int_t & id, |
249 | // Double_t & thetacut,Double_t & ratiocut1, Double_t & ratiocut2, | |
250 | // Double_t & invmasscut1,Double_t & invmasscut2) | |
251 | // { | |
77414c90 | 252 | // TParticle *particle = 0x0; |
253 | ||
254 | // Int_t iPrimary=-1; | |
255 | // Double_t El = 0, E12 = 0; | |
256 | // TLorentzVector gamma_i,gamma_j; | |
257 | // TIter next(&particleList); | |
258 | // while ( (particle = (TParticle*)next()) ) { | |
259 | // iPrimary++; | |
260 | // Int_t ksCode = particle->GetStatusCode(); | |
261 | // Int_t ksPdg = particle->GetPdgCode(); | |
262 | // if((ksCode == 1) && (iPrimary != id) && (ksPdg == 22)){ | |
263 | // particle->Momentum(gamma_i); | |
264 | // Int_t jPrimary=-1; | |
265 | // TIter next2(&particleList); | |
266 | // while ( (particle = (TParticle*)next2()) ) { | |
267 | // jPrimary++; | |
268 | // if(jPrimary>iPrimary){ | |
269 | // ksCode = particle->GetStatusCode(); | |
270 | // ksPdg = particle->GetPdgCode(); | |
271 | // if((ksCode == 1) && (iPrimary != id) && (ksPdg == 22)){ | |
272 | // particle->Momentum(gamma_j); | |
273 | // if(gamma_j.Angle(gamma_i.Vect())<thetacut){ | |
274 | // Float_t invmass = (gamma_i+gamma_j).M(); | |
275 | // h_invmass->Fill(Eg,invmass); | |
276 | // if((invmass>invmasscut1) && (invmass<invmasscut2)){ | |
277 | // E12 = (gamma_i+gamma_j).Energy(); | |
278 | // if(E12>El && (E12/Eg>ratiocut1) && (E12/Eg<ratiocut2)){ | |
279 | // //cout<<E12<<" "<<E12/Eg<<endl; | |
280 | // El = E12; | |
281 | // id_i = iPrimary; | |
282 | // id_j = jPrimary; | |
283 | // gammapair = gamma_i+gamma_j; | |
284 | // }//E12>El && (E12/Eg>0.2 && E12/Eg<1.) | |
285 | // }//(invmass>0.125) && (invmass<0.145) | |
286 | // }//gamma_j.Angle(gamma_i.Vect())<0.04 | |
287 | // }//(ksCode == 1) | |
288 | // } | |
289 | // }//while | |
290 | // // cout<<"jPrimary "<<jPrimary<<endl; | |
291 | // }// if kscode 1 | |
292 | // }//while | |
293 | // // cout<<"iPrimary "<<iPrimary<<endl; | |
90cceaf6 | 294 | // } |