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