update chamber position
[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 ---
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
34ClassImp(AliPHOSGammaJet)
35
36//____________________________________________________________________________
37AliPHOSGammaJet::AliPHOSGammaJet() {
38 // ctor
39}
40
41//____________________________________________________________________________
42AliPHOSGammaJet::AliPHOSGammaJet(const TString inputfilename) {
43 // ctor
44 AliPHOSGetter::Instance(inputfilename) ;
45}
46
47//____________________________________________________________________________
90cceaf6 48AliPHOSGammaJet::AliPHOSGammaJet(const AliPHOSGammaJet & gj) : TTask(gj) {
77414c90 49}
50
51//____________________________________________________________________________
52AliPHOSGammaJet::~AliPHOSGammaJet() {
53}
54
55//____________________________________________________________________________
e957fea8 56void 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//____________________________________________________________________________
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{
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//____________________________________________________________________________
200void 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 224void 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// }