Introducing new class : AliPHOSGammaJet
[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
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
37 ClassImp(AliPHOSGammaJet)
38
39 //____________________________________________________________________________
40 AliPHOSGammaJet::AliPHOSGammaJet() {
41   // ctor
42 }
43
44 //____________________________________________________________________________
45 AliPHOSGammaJet::AliPHOSGammaJet(const TString inputfilename) {
46   // ctor 
47   AliPHOSGetter::Instance(inputfilename) ; 
48 }
49
50 //____________________________________________________________________________
51 AliPHOSGammaJet::AliPHOSGammaJet(const AliPHOSGammaJet * gj) {
52 }
53
54 //____________________________________________________________________________
55 AliPHOSGammaJet::~AliPHOSGammaJet() {
56 }
57
58 //____________________________________________________________________________
59 void 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 //____________________________________________________________________________
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   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 //____________________________________________________________________________
199 void 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 //____________________________________________________________________________
222 void 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 //____________________________________________________________________________
246 void 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 }