Jet and Particle identification tasks moved to different directories
[u/mrichter/AliRoot.git] / PWG4 / PartCorr / AliAnaParticleJetFinderCorrelation.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: AliAnaParticleJetFinderCorrelation.cxx 22232 2007-11-17 16:39:49Z gustavo $ */
16
17 //_________________________________________________________________________
18 // Class for the analysis of particle (direct gamma) -jet (jet found with finder) correlations
19 //*-- Author: Gustavo Conesa (LNF-INFN) 
20 //////////////////////////////////////////////////////////////////////////////
21
22
23 // --- ROOT system ---
24 #include "TH2F.h"
25 #include "Riostream.h"
26
27 //---- AliRoot system ----
28 #include "AliCaloTrackReader.h"
29 #include "AliAODJet.h"
30 #include "AliAnaParticleJetFinderCorrelation.h" 
31 #include "AliLog.h"
32   
33   ClassImp(AliAnaParticleJetFinderCorrelation)
34   
35
36 //____________________________________________________________________________
37   AliAnaParticleJetFinderCorrelation::AliAnaParticleJetFinderCorrelation() : 
38         AliAnaPartCorrBaseClass(),  
39     fDeltaPhiMaxCut(0.), fDeltaPhiMinCut(0.), fRatioMaxCut(0.),  fRatioMinCut(0.), 
40         fConeSize(0), fPtThresholdInCone(0),fUseJetRefTracks(0), fMakeCorrelationInHistoMaker(0),
41         fhDeltaEta(0), fhDeltaPhi(0), fhDeltaPt(0), fhPtRatio(0), fhPt(0),
42         fhFFz(0),fhFFxi(0),fhFFpt(0),fhNTracksInCone(0)
43 {
44         //Default Ctor
45
46         //Initialize parameters
47         InitParameters();
48 }
49
50 //____________________________________________________________________________
51 AliAnaParticleJetFinderCorrelation::AliAnaParticleJetFinderCorrelation(const AliAnaParticleJetFinderCorrelation & pjf) :   
52         AliAnaPartCorrBaseClass(pjf),  
53         fDeltaPhiMaxCut(pjf.fDeltaPhiMaxCut), fDeltaPhiMinCut(pjf.fDeltaPhiMinCut), 
54         fRatioMaxCut(pjf.fRatioMaxCut), fRatioMinCut(pjf.fRatioMinCut), 
55         fConeSize(pjf.fConeSize), fPtThresholdInCone(pjf.fPtThresholdInCone),   
56         fUseJetRefTracks(pjf.fUseJetRefTracks), fMakeCorrelationInHistoMaker(pjf.fMakeCorrelationInHistoMaker),
57         fhDeltaEta(pjf.fhDeltaEta), fhDeltaPhi(pjf.fhDeltaPhi), 
58         fhDeltaPt(pjf.fhDeltaPt), fhPtRatio(pjf.fhPtRatio), fhPt(pjf.fhPt), 
59         fhFFz(pjf.fhFFz),fhFFxi(pjf.fhFFxi),fhFFpt(pjf.fhFFpt),
60         fhNTracksInCone(pjf.fhNTracksInCone)
61 {
62         // cpy ctor
63
64 }
65
66 //_________________________________________________________________________
67 AliAnaParticleJetFinderCorrelation & AliAnaParticleJetFinderCorrelation::operator = (const AliAnaParticleJetFinderCorrelation & pjf)
68 {
69         // assignment operator
70
71         if(this == &pjf)return *this;
72         ((AliAnaPartCorrBaseClass *)this)->operator=(pjf);
73   
74         fDeltaPhiMaxCut    = pjf.fDeltaPhiMaxCut ; 
75         fDeltaPhiMinCut    = pjf.fDeltaPhiMinCut ; 
76         fRatioMaxCut       = pjf.fRatioMaxCut ;
77         fRatioMinCut       = pjf.fRatioMinCut ; 
78         fConeSize          = pjf.fConeSize ; 
79         fPtThresholdInCone = pjf.fPtThresholdInCone ;
80         fUseJetRefTracks   = pjf.fUseJetRefTracks ;
81         fMakeCorrelationInHistoMaker = pjf.fMakeCorrelationInHistoMaker ;  
82           
83         //Histograms
84         fhDeltaEta = pjf.fhDeltaEta;
85         fhDeltaPhi = pjf.fhDeltaPhi;
86         fhDeltaPt  = pjf.fhDeltaPt;
87         fhPtRatio  = pjf.fhPtRatio;
88         fhPt       = pjf.fhPt;
89   
90         fhFFz   = pjf.fhFFz;
91         fhFFxi  = pjf.fhFFxi;
92         fhFFpt  = pjf.fhFFpt;
93         fhNTracksInCone = pjf.fhNTracksInCone;  
94   
95   return *this;
96
97 }
98
99 //____________________________________________________________________________
100 //AliAnaParticleJetFinderCorrelation::~AliAnaParticleJetFinderCorrelation() 
101 //{
102 //   // Remove all pointers except analysis output pointers.
103 // 
104 //}
105
106
107 //________________________________________________________________________
108 TList *  AliAnaParticleJetFinderCorrelation::GetCreateOutputObjects()
109 {  
110   // Create histograms to be saved in output file and 
111   // store them in fOutputContainer
112
113   if(GetDebug()>1) printf("Init histograms \n");
114
115   TList * outputContainer = new TList() ; 
116   outputContainer->SetName("ParticleJetFinderHistos") ; 
117
118   fhDeltaPhi  = new TH2F("DeltaPhi","#phi_{jet} - #phi_{trigger} vs p_{T trigger}",120,0,120,200,-4,4); 
119   fhDeltaPhi->SetYTitle("#Delta #phi");
120   fhDeltaPhi->SetXTitle("p_{T trigger} (GeV/c)");
121   outputContainer->Add(fhDeltaPhi);
122
123   fhDeltaEta  = new TH2F("DeltaEta","#eta_{jet} - #eta_{trigger} vs p_{T trigger}",120,0,120,200,-5,5); 
124   fhDeltaEta->SetYTitle("#Delta #eta");
125   fhDeltaEta->SetXTitle("p_{T trigger} (GeV/c)");
126   outputContainer->Add(fhDeltaEta);
127
128   fhDeltaPt  = new TH2F("DeltaPt","#p_{T trigger} - #p_{T jet} vs p_{T trigger}",120,0,120,200,-100,100); 
129   fhDeltaPt->SetYTitle("#Delta #p_{T}");
130   fhDeltaPt->SetXTitle("p_{T trigger} (GeV/c)"); 
131   outputContainer->Add(fhDeltaPt);
132
133   fhPtRatio  = new TH2F("PtRatio","#p_{T jet} / #p_{T trigger} vs p_{T trigger}",120,0,120,200,0,2.); 
134   fhPtRatio->SetYTitle("ratio");
135   fhPtRatio->SetXTitle("p_{T trigger} (GeV/c)");
136   outputContainer->Add(fhPtRatio);
137   
138   fhPt  = new TH2F("Pt","#p_{T jet} vs p_{T trigger}",120,0,120,120,0,120); 
139   fhPt->SetYTitle("#p_{T jet}(GeV/c)");
140   fhPt->SetXTitle("p_{T trigger} (GeV/c)");
141   outputContainer->Add(fhPt);
142   
143   fhFFz  = new TH2F("FFz","z = p_{T i charged}/p_{T trigger} vs p_{T trigger}", 120,0.,120.,200,0.,2);  
144   fhFFz->SetYTitle("z");
145   fhFFz->SetXTitle("p_{T trigger}");
146   outputContainer->Add(fhFFz) ;
147
148   fhFFxi  = new TH2F("FFxi","#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger}", 120,0.,120.,100,0.,10.); 
149   fhFFxi->SetYTitle("#xi");
150   fhFFxi->SetXTitle("p_{T trigger}");
151   outputContainer->Add(fhFFxi) ;
152   
153   fhFFpt  = new TH2F("FFpt","#xi = p_{T i charged}) vs p_{T trigger}", 120,0.,120.,200,0.,50.); 
154   fhFFpt->SetYTitle("p_{T charged hadron}");
155   fhFFpt->SetXTitle("p_{T trigger}");
156   outputContainer->Add(fhFFpt) ;
157
158   fhNTracksInCone  = new TH2F("NTracksInCone","#xi = p_{T i charged}) vs p_{T trigger}", 120,0.,120.,200,0.,50.); 
159   fhNTracksInCone->SetYTitle("p_{T charged hadron}");
160   fhNTracksInCone->SetXTitle("p_{T trigger}");
161   outputContainer->Add(fhNTracksInCone) ;
162         
163   return outputContainer;
164 }
165
166  //____________________________________________________________________________
167 void AliAnaParticleJetFinderCorrelation::InitParameters()
168 {
169  
170   //Initialize the parameters of the analysis.
171
172   fDeltaPhiMinCut    = 1.5 ;
173   fDeltaPhiMaxCut    = 4.5 ; 
174   fRatioMaxCut       = 1.2 ; 
175   fRatioMinCut       = 0.3 ; 
176   fConeSize          = 0.3 ;
177   fPtThresholdInCone = 0.5 ;
178   fUseJetRefTracks   = kFALSE ;
179   fMakeCorrelationInHistoMaker = kFALSE ;
180   }
181
182 //__________________________________________________________________
183 Int_t  AliAnaParticleJetFinderCorrelation::SelectJet(AliAODParticleCorrelation * particle)
184 {
185         //Returns the index of the jet that is opposite to the particle
186         
187         Int_t njets = (GetReader()->GetAOD())->GetNJets() ;     
188         cout<<">> Loop on Jets"<<endl;
189         AliAODJet * jet = new AliAODJet ;
190         Int_t index = -1;
191         for(Int_t ijet = 0; ijet < njets ; ijet++){
192                 jet = (GetReader()->GetAOD())->GetJet(ijet) ;     
193                 Float_t dphi  = TMath::Abs(particle->Phi()-jet->Phi());
194                 Float_t ratio = jet->Pt()/particle->Pt();
195                 if(GetDebug() > 3)
196                         printf("Jet %d, Ratio pT %2.3f, Delta phi %2.3f\n",ijet,ratio,dphi);      
197                 Float_t dphiprev= 10000;
198                 if((dphi > fDeltaPhiMinCut) && (dphi<fDeltaPhiMaxCut) &&
199                         (ratio > fRatioMinCut) && (ratio < fRatioMaxCut)  &&
200                         (TMath::Abs(dphi-3.14159) < TMath::Abs(dphiprev-3.14159))){//In case there is more than one jet select the most opposite.
201                         dphiprev = dphi;
202                         index = ijet ;  
203                 }//Selection cuts
204         }//AOD Jet loop
205
206         return index ;
207
208 }
209
210 //__________________________________________________________________
211 void  AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() 
212 {  
213         //Particle-Jet Correlation Analysis, fill AODs
214         
215         Int_t ntrig =  GetAODBranch()->GetEntriesFast() ;  
216     if(GetDebug() > -3){
217     printf("Begin jet finder  correlation analysis, fill AODs \n");
218     printf("In particle branch aod entries %d\n", ntrig);
219     printf("In jet      branch aod entries %d\n", (GetReader()->GetAOD())->GetNJets());
220         }
221         
222         //Loop on stored AOD particles, trigger
223         for(Int_t iaod = 0; iaod < ntrig ; iaod++){
224                 AliAODParticleCorrelation* particle =  dynamic_cast<AliAODParticleCorrelation*> (GetAODBranch()->At(iaod));
225                 //cout<<">> Test reference "<<particle->GetRefJet()<<endl;
226                         //Correlate with jets
227                 Int_t ijet = SelectJet(particle);
228                 if(ijet > -1){
229                         if(GetDebug() > 2) printf ("Jet with index %d selected \n",ijet);
230                         AliAODJet *jet = (GetReader()->GetAOD())->GetJet(ijet);
231                         particle->SetRefJet(jet);       
232                 }
233         }                 
234         
235         if(GetDebug() >1)printf("End of jet leading cone analysis, fill AODs \n");
236
237
238 //__________________________________________________________________
239 void  AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() 
240 {
241         //Particle-Jet Correlation Analysis, fill histograms
242         Int_t ntrig   =  GetAODBranch()->GetEntriesFast() ;    
243         if(GetDebug() > 1){
244                 printf("Begin jet finder  correlation analysis, fill histograms \n");
245                 printf("In particle branch aod entries %d\n", ntrig);
246                 printf("In jet      branch aod entries %d\n", (GetReader()->GetAOD())->GetNJets());
247         }
248         
249         //Loop on stored AOD particles, trigger
250         for(Int_t iaod = 0; iaod < ntrig ; iaod++){
251                 AliAODParticleCorrelation* particle =  dynamic_cast<AliAODParticleCorrelation*> (GetAODBranch()->At(iaod));
252                 //cout<<">> Test reference "<<particle->GetRefJet()<<endl;
253                 
254                 //Recover the jet correlated, found previously.
255                 AliAODJet       * jet = particle->GetJet();
256                 //If correlation not made before, do it now.
257                 if(fMakeCorrelationInHistoMaker){
258                         //Correlate with jets
259                         Int_t ijet = SelectJet(particle);
260                         if(ijet > -1){
261                                 if(GetDebug() > 2) printf ("Jet with index %d selected \n",ijet);
262                                 jet = (GetReader()->GetAOD())->GetJet(ijet);
263                                 particle->SetRefJet(jet);       
264                         }
265                 }
266                 
267                 if (!jet) continue ;
268         
269                 //Fill Histograms
270           
271                 Double_t ptTrig = particle->Pt();
272                 Double_t ptJet = jet->Pt();
273                 Double_t phiTrig = particle->Phi();
274                 Double_t phiJet = jet->Phi();
275                 Double_t etaTrig = particle->Eta();
276                 Double_t etaJet = jet->Eta();
277                 //printf("pT trigger %2.3f, pT jet %2.3f, Delta phi %2.3f, Delta eta %2.3f, Delta pT %2.3f, ratio %2.3f \n",
278                 //      ptTrig,ptJet, phiJet-phiTrig, etaJet-etaTrig, ptTrig-ptJet, ptJet/ptTrig);
279                 fhDeltaPt ->Fill(ptTrig, ptTrig-ptJet);
280                 fhDeltaPhi->Fill(ptTrig, phiJet-phiTrig);
281                 fhDeltaEta->Fill(ptTrig, etaJet-etaTrig);
282                 fhPtRatio ->Fill(ptTrig, ptJet/ptTrig);
283                 fhPt      ->Fill(ptTrig, ptJet);
284
285                 //Fragmentation function
286                 Float_t  rad = 0, pt = 0, eta = 0, phi = 0;
287                 Int_t    npartcone = 0;
288                 TVector3 p3;
289                 AliAODTrack* track = new AliAODTrack ;
290         
291                 Int_t ntracks =  0;
292                 if(!fUseJetRefTracks)
293                         ntracks =GetAODCTS()->GetEntriesFast();
294                 else //If you want to use jet tracks from JETAN
295                         ntracks =  (jet->GetRefTracks())->GetEntriesFast();
296         
297                 for(Int_t ipr = 0;ipr < ntracks ; ipr ++ ){
298                         if(!fUseJetRefTracks)
299                                 track = dynamic_cast<AliAODTrack *>(GetAODCTS()->At(ipr)) ; 
300                         else //If you want to use jet tracks from JETAN
301                                 track = dynamic_cast<AliAODTrack *>((jet->GetRefTracks())->At(ipr));
302                                         
303                         p3.SetXYZ(track->Px(),track->Py(),track->Pz());
304                         pt    = p3.Pt();
305                         eta  = p3.Eta();
306                         phi  = p3.Phi() ;
307                         if(phi < 0) phi+=TMath::TwoPi();
308                 
309                         //Check if there is any particle inside cone with pt larger than  fPtThreshold
310                         rad = TMath::Sqrt((eta-etaJet)*(eta-etaJet)+ (phi-phiJet)*(phi-phiJet));
311                         if(rad < fConeSize  && pt > fPtThresholdInCone){        
312                                 //printf("charged in jet cone pt %f, phi %f, eta %f, R %f \n",pt,phi,eta,rad);
313                                 npartcone++;
314                                 fhFFz ->Fill(ptTrig, pt/ptTrig);
315                                 fhFFxi->Fill(ptTrig, TMath::Log(ptTrig/pt));
316                                 fhFFpt->Fill(ptTrig, pt);
317                         }
318                 }//Tracks
319                 fhNTracksInCone->Fill(ptTrig, npartcone);
320         }//AOD trigger particle loop
321         if(GetDebug() >1) printf("End of jet leading cone analysis, fill histograms \n");
322
323
324
325 //__________________________________________________________________
326 void AliAnaParticleJetFinderCorrelation::Print(const Option_t * opt) const
327 {
328
329         //Print some relevant parameters set for the analysis
330         if(! opt)
331                 return;
332   
333         printf("Phi trigger-jet        <     %3.2f\n", fDeltaPhiMaxCut) ; 
334         printf("Phi trigger-jet        >     %3.2f\n", fDeltaPhiMinCut) ;
335         printf("pT Ratio trigger/jet   <     %3.2f\n", fRatioMaxCut) ; 
336         printf("pT Ratio trigger/jet   >     %3.2f\n", fRatioMinCut) ;
337         printf("fConeSize              =     %3.2f\n", fConeSize) ; 
338         printf("fPtThresholdInCone     =     %3.2f\n", fPtThresholdInCone) ;
339         printf("fUseJetRefTracks           =     %d\n",    fUseJetRefTracks) ;
340         printf("fMakeCorrelationInHistoMaker       =     %d\n",    fMakeCorrelationInHistoMaker) ;      
341
342