coding convention
[u/mrichter/AliRoot.git] / PHOS / AliPHOSGammaJet.cxx
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 #include "TRandom.h"
27
28 #include "TLorentzVector.h"
29 #include "TList.h"
30 #include "TParticle.h"
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 //____________________________________________________________________________
48 AliPHOSGammaJet::AliPHOSGammaJet(const AliPHOSGammaJet & gj) : TTask(gj) {
49 }
50
51 //____________________________________________________________________________
52 AliPHOSGammaJet::~AliPHOSGammaJet() {
53 }
54
55 //____________________________________________________________________________
56 void AliPHOSGammaJet::Exec(Option_t * opt) 
57 {
58   // does the job
59   
60   if (! opt) 
61     return ; 
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 idg = -1;
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()) ;
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 {
172   // Get the lists of jet particles and gamma
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 {
202   // Gets the leading particle from the list of charged particles
203   TParticle *particle = 0x0;
204   
205   Int_t iPrimary=-1;
206   Double_t ptmax = 0, pti = 0;
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)){
214       pti = particle->Pt(); 
215       if(pti> ptmax){
216         ptmax = pti;
217         particle->Momentum(charge);
218       }//ptmax   
219     }// kscode == 1
220   }// while
221 }
222
223 //____________________________________________________________________________
224 void AliPHOSGammaJet::GetLeadingPi0(TList &particleList, TLorentzVector &pi0) 
225 {
226   // Gets the leading pi0 from the list of particles
227   TParticle *particle = 0x0;
228   
229   Int_t iPrimary=-1;
230   Double_t ptmax = 0, pti = 0;
231   TIter next(&particleList);
232   while ( (particle = (TParticle*)next()) ) {
233     iPrimary++;  
234     Int_t ksCode = particle->GetStatusCode();
235     
236     if((ksCode == 2))
237       {
238         pti = particle->Pt(); 
239         if(pti> ptmax){
240           ptmax = pti;
241           particle->Momentum(pi0);
242         }//ptmax   
243       }// kscode == 1
244   }// while
245 }
246
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) 
251 //  {
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;
294 //  }