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 //_________________________________________________________________________
18 // Class for the analysis of gamma-jet correlations
21 //*-- Author: Gustavo Conesa & Yves Schutz (IFIC, CERN)
22 //////////////////////////////////////////////////////////////////////////////
25 // --- ROOT system ---
28 #include "TLorentzVector.h"
30 #include "TParticle.h"
31 #include "AliPHOSGammaJet.h"
32 #include "AliPHOSGetter.h"
34 ClassImp(AliPHOSGammaJet)
36 //____________________________________________________________________________
37 AliPHOSGammaJet::AliPHOSGammaJet() {
41 //____________________________________________________________________________
42 AliPHOSGammaJet::AliPHOSGammaJet(const TString inputfilename) {
44 AliPHOSGetter::Instance(inputfilename) ;
47 //____________________________________________________________________________
48 AliPHOSGammaJet::AliPHOSGammaJet(const AliPHOSGammaJet & gj) : TTask(gj) {
51 //____________________________________________________________________________
52 AliPHOSGammaJet::~AliPHOSGammaJet() {
55 //____________________________________________________________________________
56 void AliPHOSGammaJet::Exec(Option_t * opt)
63 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
64 Int_t nEvent1 = gime->MaxEvent() ;
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 ;
75 gime->Event(iEvent, "X") ;
76 Int_t nparticles = gime->NPrimaries() ;
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);
85 //Keep Stable particles within eta range
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);
94 //Decay Pi0 and keep it with different status name
96 if(particle->GetPdgCode() == 111) {
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);
114 }//final particle etacut
115 }//for (iParticle<nParticle)
116 TLorentzVector gamma,charge,pi0, gammapair;
118 GetGammaJet(particleList,gamma, idg);
119 GetLeadingCharge(particleList,charge, idg);
120 GetLeadingPi0(particleList,pi0);
121 Info("Pi0Decay", "Gamma: %f %d", gamma.Energy(), idg) ;
122 Info("Pi0Decay", "Charge: %f", charge.Energy()) ;
123 Info("Pi0Decay", "Pi0: %f", pi0.Energy()) ;
124 // GetLeadingGammaPair(particleList, gammapair, idg,
125 // 0.04,0.2, 1.0,0.13,0.14);
126 Info("Pi0Decay", "Pair: %f", gammapair.Energy()) ;
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);
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);
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;
162 //cout<<"p1: "<<p1.Px()<<" "<<p1.Py()<<" "<<p1.Pz()<<" "<<p1.E()<<endl;
164 //cout<<"p2: "<<p2.Px()<<" "<<p2.Py()<<" "<<p2.Pz()<<" "<<p2.E()<<endl;
165 //cout<<"angle"<<endl;
166 angle = p1.Angle(p2.Vect());
169 //____________________________________________________________________________
170 void AliPHOSGammaJet::GetGammaJet(TList &particleList, TLorentzVector &gamma, Int_t & id)
172 // Get the lists of jet particles and gamma
173 TParticle *particle = 0x0;
175 Int_t iPrimary=-1, id_motherg = -1;
177 TIter next(&particleList);
178 while ( (particle = (TParticle*)next()) ) {
180 Int_t ksCode = particle->GetStatusCode();
181 Int_t iMother= particle->GetMother(0);
183 if (ksCode == 21 && iMother == -1)
184 if(particle->GetPdgCode()==22){
185 id_motherg = iPrimary;
186 //cout<<"idmother "<<id_motherg<<endl;
190 if(id_motherg == particle->GetMother(0)){
191 particle->Momentum(gamma);
199 //____________________________________________________________________________
200 void AliPHOSGammaJet::GetLeadingCharge(TList &particleList, TLorentzVector &charge, Int_t & id)
202 // Gets the leading particle from the list of charged particles
203 TParticle *particle = 0x0;
206 Double_t ptmax = 0, pti = 0;
207 TIter next(&particleList);
208 while ( (particle = (TParticle*)next()) ) {
210 Int_t ksCode = particle->GetStatusCode();
212 if((ksCode == 1)&&(id != iPrimary)
213 &&(particle->GetPDG(0)->Charge()!=0)){
214 pti = particle->Pt();
217 particle->Momentum(charge);
223 //____________________________________________________________________________
224 void AliPHOSGammaJet::GetLeadingPi0(TList &particleList, TLorentzVector &pi0)
226 // Gets the leading pi0 from the list of particles
227 TParticle *particle = 0x0;
230 Double_t ptmax = 0, pti = 0;
231 TIter next(&particleList);
232 while ( (particle = (TParticle*)next()) ) {
234 Int_t ksCode = particle->GetStatusCode();
238 pti = particle->Pt();
241 particle->Momentum(pi0);
247 //____________________________________________________________________________
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)
252 // TParticle *particle = 0x0;
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()) ) {
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()) ) {
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;
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
290 // // cout<<"jPrimary "<<jPrimary<<endl;
293 // // cout<<"iPrimary "<<iPrimary<<endl;