8e01cbb289c429518dd01999fcc1b4af16f496d3
[u/mrichter/AliRoot.git] / PWGGA / CaloTrackCorrelations / AliAnaParticleJetLeadingConeCorrelation.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
16 //_________________________________________________________________________
17 // Class that contains the algorithm for the reconstruction of jet, cone around leading particle
18 // The seed is a backward particle (direct photon)
19 // 1) Take the trigger particle stored in AliAODPWG4ParticleCorrelation,
20 // 2) Search for the highest pt leading particle opposite to the photon within a phi, pt window
21 // 3) Take all particles around leading in a cone R with pt larger than threshold and construct the jet
22 //
23 //  Class created from old AliPHOSGammaJet 
24 //  (see AliRoot versions previous Release 4-09)
25 //
26 //*-- Author: Gustavo Conesa (LNF-INFN) 
27 //////////////////////////////////////////////////////////////////////////////
28
29
30 // --- ROOT system ---
31 #include "TH2F.h"
32 #include "TClonesArray.h"
33 #include "TClass.h"
34 //#include "Riostream.h"
35
36 //---- Analysis system ----
37 #include "AliVTrack.h"
38 #include "AliVCluster.h"
39 #include "AliCaloTrackReader.h"
40 #include "AliNeutralMesonSelection.h"
41 #include "AliAnaParticleJetLeadingConeCorrelation.h"  
42 #include "AliCaloPID.h"
43 #include "AliAODPWG4ParticleCorrelation.h"
44 #include "AliFiducialCut.h"
45
46 ClassImp(AliAnaParticleJetLeadingConeCorrelation)
47
48
49 //____________________________________________________________________________
50   AliAnaParticleJetLeadingConeCorrelation::AliAnaParticleJetLeadingConeCorrelation() : 
51     AliAnaCaloTrackCorrBaseClass(), fJetsOnlyInCTS(kFALSE), fPbPb(kFALSE),     
52     fSeveralConeAndPtCuts(0),  fReMakeJet(0),
53     fDeltaPhiMaxCut(0.), fDeltaPhiMinCut(0.), 
54     fLeadingRatioMaxCut(0.),  fLeadingRatioMinCut(0.), 
55     fJetCTSRatioMaxCut(0.), fJetCTSRatioMinCut(0.), 
56     fJetRatioMaxCut(0.),  fJetRatioMinCut(0.), 
57     fJetNCone(0),fJetNPt(0), fJetCone(0), 
58     fJetPtThreshold(0),fJetPtThresPbPb(0),
59     fPtTriggerSelectionCut(0.0), fSelect(0),fSelectIsolated(0),
60     fTrackVector(),fBkgMom(),fJetMom(),fJetConstMom(),
61     fLeadingMom(),fLeadingPi0Mom(),fLeadingPhoMom1(),fLeadingPhoMom2(),fLeadingChargeMom(),
62     //Histograms
63     fOutCont(0x0),
64     //Leading
65     fhChargedLeadingPt(0),fhChargedLeadingPhi(0),fhChargedLeadingEta(0),
66     fhChargedLeadingDeltaPt(0),fhChargedLeadingDeltaPhi(0),fhChargedLeadingDeltaEta(0),
67     fhChargedLeadingRatioPt(0),
68     fhNeutralLeadingPt(0),fhNeutralLeadingPhi(0),fhNeutralLeadingEta(0),
69     fhNeutralLeadingDeltaPt(0),fhNeutralLeadingDeltaPhi(0),fhNeutralLeadingDeltaEta(0),
70     fhNeutralLeadingRatioPt(0),fhChargedLeadingXi(0), fhNeutralLeadingXi(0),
71     fhChargedLeadingDeltaPhiRatioPt30(0), fhNeutralLeadingDeltaPhiRatioPt30(0),
72     fhChargedLeadingDeltaPhiRatioPt50(0), fhNeutralLeadingDeltaPhiRatioPt50(0),
73     //Jet
74     fhJetPt(0),fhJetRatioPt(0),fhJetDeltaPhi(0), fhJetDeltaEta(0),
75     fhJetLeadingRatioPt(0),fhJetLeadingDeltaPhi(0),fhJetLeadingDeltaEta(0),
76     fhJetFFz(0),fhJetFFxi(0),fhJetFFpt(0),fhJetNTracksInCone(0),
77     fhBkgPt(0),fhBkgRatioPt(0),fhBkgDeltaPhi(0), fhBkgDeltaEta(0),
78     fhBkgLeadingRatioPt(0),fhBkgLeadingDeltaPhi(0),fhBkgLeadingDeltaEta(0),
79     fhBkgFFz(0),fhBkgFFxi(0),fhBkgFFpt(0),fhBkgNTracksInCone(0),
80     //Several cones and thres histograms
81     fhJetPts(),fhJetRatioPts(),fhJetDeltaPhis(), fhJetDeltaEtas(),
82     fhJetLeadingRatioPts(),fhJetLeadingDeltaPhis(),fhJetLeadingDeltaEtas(),
83     fhJetFFzs(),fhJetFFxis(),fhJetFFpts(),fhJetNTracksInCones(),
84     fhBkgPts(),fhBkgRatioPts(),fhBkgDeltaPhis(), fhBkgDeltaEtas(),
85     fhBkgLeadingRatioPts(),fhBkgLeadingDeltaPhis(),fhBkgLeadingDeltaEtas(),
86     fhBkgFFzs(),fhBkgFFxis(),fhBkgFFpts(),fhBkgNTracksInCones()
87 {
88   //Default Ctor
89   
90   //Initialize parameters
91
92   for(Int_t i = 0; i<6; i++){
93     fJetXMin1[i]     = 0.0 ;
94     fJetXMin2[i]     = 0.0 ;
95     fJetXMax1[i]     = 0.0 ;
96     fJetXMax2[i]     = 0.0 ;
97     fBkgMean[i]      = 0.0 ;
98     fBkgRMS[i]       = 0.0 ;
99     if( i < 2 ){
100       fJetE1[i]        = 0.0 ;
101       fJetE2[i]        = 0.0 ;
102       fJetSigma1[i]    = 0.0 ;
103       fJetSigma2[i]    = 0.0 ;
104     }
105   }
106   
107   //Several cones and thres histograms
108   for(Int_t i = 0; i<5; i++){
109     fJetCones[i]         = 0.0 ;
110     fJetNameCones[i]     = ""  ;
111     fJetPtThres[i]      = 0.0 ;
112     fJetNamePtThres[i]  = ""  ;
113     for(Int_t j = 0; j<5; j++){
114       fhJetPts[i][j]=0 ;
115       fhJetRatioPts[i][j]=0 ;
116       fhJetDeltaPhis[i][j]=0 ; 
117       fhJetDeltaEtas[i][j]=0 ;
118       fhJetLeadingRatioPts[i][j]=0 ;
119       fhJetLeadingDeltaPhis[i][j]=0 ;
120       fhJetLeadingDeltaEtas[i][j]=0 ;
121       fhJetFFzs[i][j]=0 ;
122       fhJetFFxis[i][j]=0 ;
123       fhJetFFpts[i][j]=0 ;
124       fhJetNTracksInCones[i][j]=0 ;
125       fhBkgPts[i][j]=0 ;
126       fhBkgRatioPts[i][j]=0 ;
127       fhBkgDeltaPhis[i][j]=0 ; 
128       fhBkgDeltaEtas[i][j]=0 ;
129       fhBkgLeadingRatioPts[i][j]=0 ;
130       fhBkgLeadingDeltaPhis[i][j]=0 ;
131       fhBkgLeadingDeltaEtas[i][j]=0 ;
132       fhBkgFFzs[i][j]=0 ;
133       fhBkgFFxis[i][j]=0 ;
134       fhBkgFFpts[i][j]=0 ;
135       fhBkgNTracksInCones[i][j]=0 ;
136     }
137   }
138
139   InitParameters();
140
141 }
142 /*
143 //____________________________________________________________________________
144 AliAnaParticleJetLeadingConeCorrelation::AliAnaParticleJetLeadingConeCorrelation(const AliAnaParticleJetLeadingConeCorrelation & jetlc) :   
145   AliAnaCaloTrackCorrBaseClass(jetlc), fJetsOnlyInCTS(jetlc.fJetsOnlyInCTS), fPbPb(jetlc.fPbPb), 
146   fSeveralConeAndPtCuts(jetlc.fSeveralConeAndPtCuts),  fReMakeJet(jetlc. fReMakeJet),
147   fDeltaPhiMaxCut(jetlc. fDeltaPhiMaxCut), fDeltaPhiMinCut(jetlc.fDeltaPhiMinCut), 
148   fLeadingRatioMaxCut(jetlc.fLeadingRatioMaxCut),  fLeadingRatioMinCut(jetlc.fLeadingRatioMinCut), 
149   fJetCTSRatioMaxCut(jetlc.fJetCTSRatioMaxCut),
150   fJetCTSRatioMinCut(jetlc.fJetCTSRatioMinCut), fJetRatioMaxCut(jetlc.fJetRatioMaxCut),
151   fJetRatioMinCut(jetlc.fJetRatioMinCut),  fJetNCone(jetlc.fJetNCone),
152   fJetNPt(jetlc.fJetNPt), fJetCone(jetlc.fJetCone),
153   fJetPtThreshold(jetlc.fJetPtThreshold),fJetPtThresPbPb(jetlc.fJetPtThresPbPb),
154   fPtTriggerSelectionCut(jetlc.fPtTriggerSelectionCut), fSelect(jetlc.fSelect), 
155   fSelectIsolated(jetlc.fSelectIsolated),    
156   //Histograms
157   fOutCont(jetlc. fOutCont),
158   //Leading
159   fhChargedLeadingPt(jetlc.fhChargedLeadingPt), fhChargedLeadingPhi(jetlc.fhChargedLeadingPhi),
160   fhChargedLeadingEta(jetlc.fhChargedLeadingEta), fhChargedLeadingDeltaPt(jetlc.fhChargedLeadingDeltaPt),
161   fhChargedLeadingDeltaPhi(jetlc.fhChargedLeadingDeltaPhi),fhChargedLeadingDeltaEta(jetlc.fhChargedLeadingDeltaEta),
162   fhChargedLeadingRatioPt(jetlc.fhChargedLeadingRatioPt),
163   fhNeutralLeadingPt(jetlc.fhNeutralLeadingPt),fhNeutralLeadingPhi(jetlc.fhNeutralLeadingPhi),
164   fhNeutralLeadingEta(jetlc.fhNeutralLeadingEta), fhNeutralLeadingDeltaPt(jetlc.fhNeutralLeadingDeltaPt),
165   fhNeutralLeadingDeltaPhi(jetlc.fhNeutralLeadingDeltaPhi),fhNeutralLeadingDeltaEta(jetlc.fhNeutralLeadingDeltaEta),
166   fhNeutralLeadingRatioPt(jetlc.fhNeutralLeadingRatioPt),
167   fhChargedLeadingXi(jetlc.fhChargedLeadingXi), fhNeutralLeadingXi(jetlc.fhNeutralLeadingXi),
168   fhChargedLeadingDeltaPhiRatioPt30(jetlc.fhChargedLeadingDeltaPhiRatioPt30), fhNeutralLeadingDeltaPhiRatioPt30(jetlc.fhNeutralLeadingDeltaPhiRatioPt30),
169   fhChargedLeadingDeltaPhiRatioPt50(jetlc.fhChargedLeadingDeltaPhiRatioPt50), fhNeutralLeadingDeltaPhiRatioPt50(jetlc.fhNeutralLeadingDeltaPhiRatioPt50),
170   //Jet
171   fhJetPt(jetlc.fhJetPt),fhJetRatioPt(jetlc.fhJetRatioPt),fhJetDeltaPhi(jetlc.fhJetDeltaPhi), 
172   fhJetDeltaEta(jetlc.fhJetDeltaEta), fhJetLeadingRatioPt(jetlc.fhJetLeadingRatioPt),
173   fhJetLeadingDeltaPhi(jetlc.fhJetLeadingDeltaPhi),fhJetLeadingDeltaEta(jetlc.fhJetLeadingDeltaEta),
174   fhJetFFz(jetlc.fhJetFFz),fhJetFFxi(jetlc.fhJetFFxi),fhJetFFpt(jetlc.fhJetFFpt),
175   fhJetNTracksInCone(jetlc.fhJetNTracksInCone),
176   fhBkgPt(jetlc.fhBkgPt),fhBkgRatioPt(jetlc.fhBkgRatioPt),fhBkgDeltaPhi(jetlc.fhBkgDeltaPhi), 
177   fhBkgDeltaEta(jetlc.fhBkgDeltaEta), fhBkgLeadingRatioPt(jetlc.fhBkgLeadingRatioPt),
178   fhBkgLeadingDeltaPhi(jetlc.fhBkgLeadingDeltaPhi),fhBkgLeadingDeltaEta(jetlc.fhBkgLeadingDeltaEta),
179   fhBkgFFz(jetlc.fhBkgFFz),fhBkgFFxi(jetlc.fhBkgFFxi),fhBkgFFpt(jetlc.fhBkgFFpt),
180   fhBkgNTracksInCone(jetlc.fhBkgNTracksInCone),
181   //Several cones and thres histograms
182   fhJetPts(),fhJetRatioPts(),fhJetDeltaPhis(), fhJetDeltaEtas(),
183   fhJetLeadingRatioPts(),fhJetLeadingDeltaPhis(),fhJetLeadingDeltaEtas(),
184   fhJetFFzs(),fhJetFFxis(),fhJetFFpts(),fhJetNTracksInCones(),
185   fhBkgPts(),fhBkgRatioPts(),fhBkgDeltaPhis(), fhBkgDeltaEtas(),
186   fhBkgLeadingRatioPts(),fhBkgLeadingDeltaPhis(),fhBkgLeadingDeltaEtas(),
187   fhBkgFFzs(),fhBkgFFxis(),fhBkgFFpts(),fhBkgNTracksInCones()
188 {
189   // cpy ctor
190
191   for(Int_t i = 0; i<6; i++){
192     fJetXMin1[i]       = jetlc.fJetXMin1[i] ;
193     fJetXMin2[i]       = jetlc.fJetXMin2[i] ;
194     fJetXMax1[i]       = jetlc.fJetXMax1[i] ;
195     fJetXMax2[i]       = jetlc.fJetXMax2[i] ;
196     fBkgMean[i]        = jetlc.fBkgMean[i] ;
197     fBkgRMS[i]         = jetlc.fBkgRMS[i] ;
198     if( i < 2 ){
199       fJetE1[i]        = jetlc.fJetE1[i] ;
200       fJetE2[i]        = jetlc.fJetE2[i] ;
201       fJetSigma1[i]    = jetlc.fJetSigma1[i] ;
202       fJetSigma2[i]    = jetlc.fJetSigma2[i] ;
203     }
204   }          
205   
206   //Several cones and thres histograms
207   for(Int_t i = 0; i<5; i++){
208     fJetCones[i]        = jetlc.fJetCones[i] ;
209     fJetNameCones[i]    = jetlc.fJetNameCones[i] ;
210     fJetPtThres[i]      = jetlc.fJetPtThres[i] ;
211     fJetNamePtThres[i]  = jetlc.fJetNamePtThres[i] ;
212     for(Int_t j = 0; j<5; j++){
213       fhJetPts[i][j] = jetlc.fhJetPts[i][j] ;
214       fhJetRatioPts[i][j] = jetlc.fhJetRatioPts[i][j] ;
215       fhJetDeltaPhis[i][j] = jetlc.fhJetDeltaPhis[i][j] ; 
216       fhJetDeltaEtas[i][j] = jetlc.fhJetDeltaEtas[i][j] ;
217       fhJetLeadingRatioPts[i][j] = jetlc.fhJetLeadingRatioPts[i][j] ;
218       fhJetLeadingDeltaPhis[i][j] = jetlc.fhJetLeadingDeltaPhis[i][j] ;
219       fhJetLeadingDeltaEtas[i][j] = jetlc.fhJetLeadingDeltaEtas[i][j] ;
220       fhJetFFzs[i][j] = jetlc.fhJetFFzs[i][j] ;
221       fhJetFFxis[i][j] = jetlc.fhJetFFxis[i][j] ;
222       fhJetFFpts[i][j] = jetlc.fhJetFFpts[i][j] ;
223       fhJetNTracksInCones[i][j] = fhJetNTracksInCones[i][j] ;
224       fhBkgPts[i][j] = jetlc.fhBkgPts[i][j] ;
225       fhBkgRatioPts[i][j] = jetlc.fhBkgRatioPts[i][j] ;
226       fhBkgDeltaPhis[i][j] = jetlc.fhBkgDeltaPhis[i][j] ; 
227       fhBkgDeltaEtas[i][j] = jetlc.fhBkgDeltaEtas[i][j] ;
228       fhBkgLeadingRatioPts[i][j] = jetlc.fhBkgLeadingRatioPts[i][j] ;
229       fhBkgLeadingDeltaPhis[i][j] = jetlc.fhBkgLeadingDeltaPhis[i][j] ;
230       fhBkgLeadingDeltaEtas[i][j] = jetlc.fhBkgLeadingDeltaEtas[i][j] ;
231       fhBkgFFzs[i][j] = jetlc.fhBkgFFzs[i][j] ;
232       fhBkgFFxis[i][j] = jetlc.fhBkgFFxis[i][j] ;
233       fhBkgFFpts[i][j] = jetlc.fhBkgFFpts[i][j] ;
234       fhBkgNTracksInCones[i][j] = jetlc.fhBkgNTracksInCones[i][j] ;
235     }
236   }  
237 }
238
239 //_________________________________________________________________________
240 AliAnaParticleJetLeadingConeCorrelation & AliAnaParticleJetLeadingConeCorrelation::operator = (const AliAnaParticleJetLeadingConeCorrelation & jetlc)
241 {
242   // assignment operator
243   
244   if(this == &jetlc)return *this;
245   ((AliAnaCaloTrackCorrBaseClass *)this)->operator=(jetlc);
246   
247   fSeveralConeAndPtCuts  = jetlc.fSeveralConeAndPtCuts ; 
248   fPbPb                  = jetlc.fPbPb ;
249   fReMakeJet             = jetlc.fReMakeJet ;
250   fJetsOnlyInCTS         = jetlc.fJetsOnlyInCTS;
251
252   fDeltaPhiMaxCut        = jetlc.fDeltaPhiMaxCut ; 
253   fDeltaPhiMinCut        = jetlc.fDeltaPhiMinCut ; 
254   fLeadingRatioMaxCut    = jetlc.fLeadingRatioMaxCut ;
255   fLeadingRatioMinCut    = jetlc.fLeadingRatioMinCut ;
256   
257   fJetCTSRatioMaxCut     = jetlc.fJetCTSRatioMaxCut ;
258   fJetCTSRatioMinCut     = jetlc.fJetCTSRatioMinCut ; 
259   fJetRatioMaxCut        = jetlc.fJetRatioMaxCut ;
260   fJetRatioMinCut        = jetlc.fJetRatioMinCut ; 
261  
262   fJetNCone              = jetlc.fJetNCone ;
263   fJetNPt                = jetlc.fJetNPt ; fJetCone = jetlc.fJetCone ; 
264   fJetPtThreshold        = jetlc.fJetPtThreshold ;
265   fJetPtThresPbPb        = jetlc.fJetPtThresPbPb ;
266   fPtTriggerSelectionCut = jetlc.fPtTriggerSelectionCut ;
267   fSelect                = jetlc.fSelect ;  
268   fSelectIsolated        = jetlc.fSelectIsolated ;
269
270   for(Int_t i = 0; i<6; i++){
271     fJetXMin1[i]       = jetlc.fJetXMin1[i] ;
272     fJetXMin2[i]       = jetlc.fJetXMin2[i] ;
273     fJetXMax1[i]       = jetlc.fJetXMax1[i] ;
274     fJetXMax2[i]       = jetlc.fJetXMax2[i] ;
275     fBkgMean[i]        = jetlc.fBkgMean[i] ;
276     fBkgRMS[i]         = jetlc.fBkgRMS[i] ;
277     if( i < 2 ){
278       fJetE1[i]        = jetlc.fJetE1[i] ;
279       fJetE2[i]        = jetlc.fJetE2[i] ;
280       fJetSigma1[i]    = jetlc.fJetSigma1[i] ;
281       fJetSigma2[i]    = jetlc.fJetSigma2[i] ;
282     }
283   }     
284
285   //Histograms
286   fOutCont = jetlc. fOutCont ;
287   fhChargedLeadingPt = jetlc.fhChargedLeadingPt; fhChargedLeadingPhi = jetlc.fhChargedLeadingPhi;
288   fhChargedLeadingEta = jetlc.fhChargedLeadingEta; fhChargedLeadingDeltaPt = jetlc.fhChargedLeadingDeltaPt;
289   fhChargedLeadingDeltaPhi = jetlc.fhChargedLeadingDeltaPhi;fhChargedLeadingDeltaEta = jetlc.fhChargedLeadingDeltaEta;
290   fhChargedLeadingRatioPt = jetlc.fhChargedLeadingRatioPt;
291   fhNeutralLeadingPt = jetlc.fhNeutralLeadingPt;fhNeutralLeadingPhi = jetlc.fhNeutralLeadingPhi;
292   fhNeutralLeadingEta = jetlc.fhNeutralLeadingEta; fhNeutralLeadingDeltaPt = jetlc.fhNeutralLeadingDeltaPt;
293   fhNeutralLeadingDeltaPhi = jetlc.fhNeutralLeadingDeltaPhi;fhNeutralLeadingDeltaEta = jetlc.fhNeutralLeadingDeltaEta;
294   fhNeutralLeadingRatioPt = jetlc.fhNeutralLeadingRatioPt;
295   fhChargedLeadingXi = jetlc.fhChargedLeadingXi; 
296   fhNeutralLeadingXi = jetlc.fhNeutralLeadingXi;
297
298   fhChargedLeadingDeltaPhiRatioPt30 = jetlc.fhChargedLeadingDeltaPhiRatioPt30; 
299   fhNeutralLeadingDeltaPhiRatioPt30 = jetlc.fhNeutralLeadingDeltaPhiRatioPt30;
300   fhChargedLeadingDeltaPhiRatioPt50 = jetlc.fhChargedLeadingDeltaPhiRatioPt50; 
301   fhNeutralLeadingDeltaPhiRatioPt50 = jetlc.fhNeutralLeadingDeltaPhiRatioPt50;  
302         
303   fhJetPt = jetlc.fhJetPt;fhJetRatioPt = jetlc.fhJetRatioPt;fhJetDeltaPhi = jetlc.fhJetDeltaPhi; 
304   fhJetDeltaEta = jetlc.fhJetDeltaEta; fhJetLeadingRatioPt = jetlc.fhJetLeadingRatioPt;
305   fhJetLeadingDeltaPhi = jetlc.fhJetLeadingDeltaPhi;fhJetLeadingDeltaEta = jetlc.fhJetLeadingDeltaEta;
306   fhJetFFz = jetlc.fhJetFFz;fhJetFFxi = jetlc.fhJetFFxi;fhJetFFpt = jetlc.fhJetFFpt;
307   fhJetNTracksInCone = jetlc.fhJetNTracksInCone;
308   fhBkgPt = jetlc.fhBkgPt;fhBkgRatioPt = jetlc.fhBkgRatioPt;fhBkgDeltaPhi = jetlc.fhBkgDeltaPhi; 
309   fhBkgDeltaEta = jetlc.fhBkgDeltaEta; fhBkgLeadingRatioPt = jetlc.fhBkgLeadingRatioPt;
310   fhBkgLeadingDeltaPhi = jetlc.fhBkgLeadingDeltaPhi;fhBkgLeadingDeltaEta = jetlc.fhBkgLeadingDeltaEta;
311   fhBkgFFz = jetlc.fhBkgFFz;fhBkgFFxi = jetlc.fhBkgFFxi;fhBkgFFpt = jetlc.fhBkgFFpt;
312   fhBkgNTracksInCone = jetlc.fhBkgNTracksInCone;
313
314
315   //Several cones and thres histograms
316   for(Int_t i = 0; i<5; i++){
317     fJetCones[i]        = jetlc.fJetCones[i] ;
318     fJetNameCones[i]    = jetlc.fJetNameCones[i] ;
319     fJetPtThres[i]      = jetlc.fJetPtThres[i] ;
320     fJetNamePtThres[i]  = jetlc.fJetNamePtThres[i] ;
321     
322     for(Int_t j = 0; j<5; j++){
323       fhJetPts[i][j] = jetlc.fhJetPts[i][j] ;
324       fhJetRatioPts[i][j] = jetlc.fhJetRatioPts[i][j] ;
325       fhJetDeltaPhis[i][j] = jetlc.fhJetDeltaPhis[i][j] ; 
326       fhJetDeltaEtas[i][j] = jetlc.fhJetDeltaEtas[i][j] ;
327       fhJetLeadingRatioPts[i][j] = jetlc.fhJetLeadingRatioPts[i][j] ;
328       fhJetLeadingDeltaPhis[i][j] = jetlc.fhJetLeadingDeltaPhis[i][j] ;
329       fhJetLeadingDeltaEtas[i][j] = jetlc.fhJetLeadingDeltaEtas[i][j] ;
330       fhJetFFzs[i][j] = jetlc.fhJetFFzs[i][j] ;
331       fhJetFFxis[i][j] = jetlc.fhJetFFxis[i][j] ;
332       fhJetFFpts[i][j] = jetlc.fhJetFFpts[i][j] ;
333       fhJetNTracksInCones[i][j] = fhJetNTracksInCones[i][j] ;
334       fhBkgPts[i][j] = jetlc.fhBkgPts[i][j] ;
335       fhBkgRatioPts[i][j] = jetlc.fhBkgRatioPts[i][j] ;
336       fhBkgDeltaPhis[i][j] = jetlc.fhBkgDeltaPhis[i][j] ; 
337       fhBkgDeltaEtas[i][j] = jetlc.fhBkgDeltaEtas[i][j] ;
338       fhBkgLeadingRatioPts[i][j] = jetlc.fhBkgLeadingRatioPts[i][j] ;
339       fhBkgLeadingDeltaPhis[i][j] = jetlc.fhBkgLeadingDeltaPhis[i][j] ;
340       fhBkgLeadingDeltaEtas[i][j] = jetlc.fhBkgLeadingDeltaEtas[i][j] ;
341       fhBkgFFzs[i][j] = jetlc.fhBkgFFzs[i][j] ;
342       fhBkgFFxis[i][j] = jetlc.fhBkgFFxis[i][j] ;
343       fhBkgFFpts[i][j] = jetlc.fhBkgFFpts[i][j] ;
344       fhBkgNTracksInCones[i][j] = jetlc.fhBkgNTracksInCones[i][j] ;
345     }
346   }      
347
348   return *this;
349
350 }
351 */
352 //____________________________________________________________________________
353 AliAnaParticleJetLeadingConeCorrelation::~AliAnaParticleJetLeadingConeCorrelation() 
354 {
355    // Remove all pointers except analysis output pointers.
356 }
357
358 //____________________________________________________________________________
359 Double_t AliAnaParticleJetLeadingConeCorrelation::CalculateJetRatioLimit(const Double_t ptg, const Double_t *par, const Double_t *x) const {
360   //Calculate the ratio of the jet and trigger particle limit for the selection
361   //WARNING: need to check what it does
362   //printf("CalculateLimit: x1 %2.3f, x2%2.3f\n",x[0],x[1]);
363   Double_t ePP = par[0] + par[1] * ptg ;
364   Double_t sPP = par[2] + par[3] * ptg ;
365   Double_t f   = x[0]   + x[1]   * ptg ;
366   Double_t ePbPb = ePP + par[4] ;
367   Double_t sPbPb = TMath::Sqrt(sPP*sPP+ par[5]*par[5]) ;
368   Double_t rat = (ePbPb - sPbPb * f) / ptg ;
369   //printf("CalculateLimit: ePP %2.3f, sPP %2.3f, f %2.3f\n", ePP, sPP, f);
370   //printf("CalculateLimit: ePbPb %2.3f, sPbPb %2.3f, rat %2.3f\n", ePbPb, sPbPb, rat);
371   return rat ;
372 }
373
374 //___________________________________________________________________________________________________
375 void AliAnaParticleJetLeadingConeCorrelation::FillJetHistos(AliAODPWG4ParticleCorrelation * particle,
376                                                             const TLorentzVector jet,
377                                                             const TString & type, const TString & lastname)
378 {
379   //Fill jet and background histograms 
380   Double_t ptTrig  = particle->Pt();
381   Double_t ptJet   = jet.Pt();
382   Double_t ptLead  = fLeadingMom.Pt();
383   Double_t phiTrig = particle->Phi();
384   Double_t phiJet  = jet.Phi();
385   if(phiJet < 0) phiJet+=TMath::TwoPi();
386   Double_t phiLead = fLeadingMom.Phi();
387   if(phiLead < 0) phiLead+=TMath::TwoPi();
388   Double_t etaTrig = particle->Eta();
389   Double_t etaJet  = jet.Eta();
390   Double_t etaLead = fLeadingMom.Eta();
391   
392   TH2F *h1 = 0x0;
393   h1 = dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sPt%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
394   if(h1)h1->Fill(ptTrig,ptJet);
395   
396   TH2F *h2 = 0x0;
397   h2 = dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sRatioPt%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
398   if(h2) h2->Fill(ptTrig,ptJet/ptTrig);
399   
400   TH2F *h3 = 0x0;
401   h3 = dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sLeadingRatioPt%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
402   if(h3)h3->Fill(ptTrig,ptLead/ptJet);
403   
404   //   dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sPhi%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())))->
405   //     Fill(ptTrig,phiJet);
406   TH2F *h4 = 0x0;
407   h4 = dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sDeltaPhi%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
408   if(h4) h4->Fill(ptTrig,phiJet-phiTrig);
409   TH2F *h5 = 0x0;
410   h5 = dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sLeadingDeltaPhi%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
411   if(h5) h5->Fill(ptTrig,phiJet-phiLead);
412   
413   //   dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sEta%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())))->
414   //     Fill(ptTrig,etaJet);
415   TH2F *h6 = 0x0;
416   h6 = dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sDeltaEta%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
417   if(h6) h6->Fill(ptTrig,etaJet-etaTrig);
418   TH2F *h7 = 0x0;
419   h7 = dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sLeadingDeltaEta%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
420   if(h7) h7->Fill(ptTrig,etaJet-etaLead);
421   
422   //Construct fragmentation function
423   TObjArray * pl = new TObjArray;
424   
425   if(type == "Jet") pl = particle->GetObjArray(Form("%sTracks",GetAODObjArrayName().Data()));
426   else if(type == "Bkg") particle->GetObjArray(Form("%sTracksBkg",GetAODObjArrayName().Data()));
427   
428   if(!pl) return ;
429   
430   //Different pt cut for jet particles in different collisions systems
431   //Only needed when jet is recalculated from AODs
432   //Float_t ptcut = fJetPtThreshold;
433   //if(fPbPb && !fSeveralConeAndPtCuts && ptTrig > fPtTriggerSelectionCut)  ptcut = fJetPtThresPbPb ;
434   
435   Int_t nTracksInCone = 0;
436   
437   for(Int_t ipr = 0;ipr < pl->GetEntriesFast() ; ipr ++ )
438   {
439     AliVTrack* track = dynamic_cast<AliVTrack *>(pl->At(ipr)) ;
440     if(track)fTrackVector.SetXYZ(track->Px(),track->Py(),track->Pz());
441     else printf("AliAnaParticleJetLeadingConeCorrelation::FillJetHistos() - Track not available\n");
442     
443     //Recheck if particle is in jet cone
444     if(fReMakeJet || fSeveralConeAndPtCuts)
445       if(!IsParticleInJetCone(fTrackVector.Eta(), fTrackVector.Phi(), fLeadingMom.Eta(), fLeadingMom.Phi()) ) continue ;
446     
447     nTracksInCone++; 
448     
449     TH2F *ha =dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sFFz%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
450     if(ha) ha->Fill(ptTrig,fTrackVector.Pt()/ptTrig);
451     TH2F *hb  =dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sFFxi%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
452     if(hb) hb->Fill(ptTrig,TMath::Log(ptTrig/fTrackVector.Pt()));
453     TH2F *hc =dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sFFpt%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
454     if(hc) hc->Fill(ptTrig,fTrackVector.Pt());
455     
456   }//track loop
457   
458   if(nTracksInCone > 0)
459   {
460     TH2F *hd = dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sNTracksInCone%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
461     if(hd)hd->Fill(ptTrig, nTracksInCone);
462   }
463   
464 }
465
466 //________________________________________________________________________
467 TList *  AliAnaParticleJetLeadingConeCorrelation::GetCreateOutputObjects()
468 {  
469   // Create histograms to be saved in output file and 
470   // store them in fOutCont
471   
472   if(GetDebug()>1) printf("AliAnaParticleJetLeadingConeCorrelation::GetCreateOutputObjects() - Init histograms \n");
473   
474   fOutCont = new TList() ; 
475   fOutCont->SetName("ParticleJetLeadingInConeCorrelationHistograms") ; 
476   
477   Int_t nptbins  = GetHistogramRanges()->GetHistoPtBins();
478   Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins();
479   Int_t netabins = GetHistogramRanges()->GetHistoEtaBins();
480   Float_t ptmax  = GetHistogramRanges()->GetHistoPtMax();
481   Float_t phimax = GetHistogramRanges()->GetHistoPhiMax();
482   Float_t etamax = GetHistogramRanges()->GetHistoEtaMax();
483   Float_t ptmin  = GetHistogramRanges()->GetHistoPtMin();
484   Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
485   Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();      
486   
487   fhChargedLeadingPt  = new TH2F("ChargedLeadingPt","p_{T leading charge} vs p_{T trigger}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
488   fhChargedLeadingPt->SetYTitle("p_{T leading charge}");
489   fhChargedLeadingPt->SetXTitle("p_{T trigger} (GeV/c)");
490   
491   fhChargedLeadingPhi  = new TH2F("ChargedLeadingPhi","#phi_{h^{#pm}}  vs p_{T trigger}", nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
492   fhChargedLeadingPhi->SetYTitle("#phi_{h^{#pm}} (rad)");
493   fhChargedLeadingPhi->SetXTitle("p_{T trigger} (GeV/c)");
494   
495   fhChargedLeadingEta  = new TH2F("ChargedLeadingEta","#eta_{h^{#pm}}  vs p_{T trigger}",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
496   fhChargedLeadingEta->SetYTitle("#eta_{h^{#pm}} ");
497   fhChargedLeadingEta->SetXTitle("p_{T trigger} (GeV/c)");
498   
499   fhChargedLeadingDeltaPt  = new TH2F("ChargedLeadingDeltaPt","p_{T trigger} - p_{T h^{#pm}} vs p_{T trigger}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
500   fhChargedLeadingDeltaPt->SetYTitle("#Delta p_{T} (GeV/c)");
501   fhChargedLeadingDeltaPt->SetXTitle("p_{T trigger} (GeV/c)");
502   
503   fhChargedLeadingDeltaPhi  = new TH2F("ChargedLeadingDeltaPhi","#phi_{trigger} - #phi_{h^{#pm}} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,TMath::TwoPi()); 
504   fhChargedLeadingDeltaPhi->SetYTitle("#Delta #phi (rad)");
505   fhChargedLeadingDeltaPhi->SetXTitle("p_{T trigger} (GeV/c)");
506   
507   fhChargedLeadingDeltaEta  = new TH2F("ChargedLeadingDeltaEta","#eta_{trigger} - #eta_{h^{#pm}} vs p_{T trigger}",nptbins,ptmin,ptmax,120,-2,2); 
508   fhChargedLeadingDeltaEta->SetYTitle("#Delta #eta");
509   fhChargedLeadingDeltaEta->SetXTitle("p_{T trigger} (GeV/c)");
510   
511   fhChargedLeadingRatioPt  = new TH2F("ChargedLeadingRatioPt","p_{T leading charge} /p_{T trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,2); 
512   fhChargedLeadingRatioPt->SetYTitle("p_{T lead charge} /p_{T trigger}");
513   fhChargedLeadingRatioPt->SetXTitle("p_{T trigger} (GeV/c)");
514   
515   fhChargedLeadingXi  = new TH2F("ChargedLeadingXi","ln(p_{T trigger} / p_{T leading charge} ) vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,10); 
516   fhChargedLeadingXi->SetYTitle("#xi");
517   fhChargedLeadingXi->SetXTitle("p_{T trigger} (GeV/c)");
518         
519   fOutCont->Add(fhChargedLeadingPt) ;
520   fOutCont->Add(fhChargedLeadingPhi) ;
521   fOutCont->Add(fhChargedLeadingEta) ;
522   fOutCont->Add(fhChargedLeadingDeltaPt) ; 
523   fOutCont->Add(fhChargedLeadingDeltaPhi) ; 
524   fOutCont->Add(fhChargedLeadingDeltaEta) ; 
525   fOutCont->Add(fhChargedLeadingRatioPt) ;
526   fOutCont->Add(fhChargedLeadingXi) ;
527         
528   fhChargedLeadingDeltaPhiRatioPt30  = new TH2F("ChargedLeadingDeltaPhiRatioPt30","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T leading}/p_{T trigger}, charged leading, p_{T trigger} > 30 GeV/c",120,0,TMath::TwoPi(),nptbins,0,1); 
529   fhChargedLeadingDeltaPhiRatioPt30->SetXTitle("#Delta #phi (rad)");
530   fhChargedLeadingDeltaPhiRatioPt30->SetYTitle("p_{T leading} / p_{T trigger}");        
531         
532   fhChargedLeadingDeltaPhiRatioPt50  = new TH2F("ChargedLeadingDeltaPhiRatioPt50","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T leading}/p_{T trigger}, charged leading, p_{T trigger} > 50 GeV/c",120,0,TMath::TwoPi(),nptbins,0,1); 
533   fhChargedLeadingDeltaPhiRatioPt50->SetXTitle("#Delta #phi (rad)");
534   fhChargedLeadingDeltaPhiRatioPt50->SetYTitle("p_{T leading} / p_{T trigger}");          
535         
536   fOutCont->Add(fhChargedLeadingDeltaPhiRatioPt30) ; 
537   fOutCont->Add(fhChargedLeadingDeltaPhiRatioPt50) ; 
538         
539   if(!fJetsOnlyInCTS){
540     
541     fhNeutralLeadingPt  = new TH2F("NeutralLeadingPt","p_{T leading #pi^{0}} vs p_{T trigger}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
542     fhNeutralLeadingPt->SetYTitle("p_{T leading #pi^{0}}");
543     fhNeutralLeadingPt->SetXTitle("p_{T trigger} (GeV/c)");
544     
545     fhNeutralLeadingPhi  = new TH2F("NeutralLeadingPhi","#phi_{#pi^{0}}  vs p_{T trigger}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
546     fhNeutralLeadingPhi->SetYTitle("#phi_{#pi^{0}} (rad)");
547     fhNeutralLeadingPhi->SetXTitle("p_{T trigger} (GeV/c)");
548     
549     fhNeutralLeadingEta  = new TH2F("NeutralLeadingEta","#eta_{#pi^{0}}  vs p_{T trigger}",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
550     fhNeutralLeadingEta->SetYTitle("#eta_{#pi^{0}} ");
551     fhNeutralLeadingEta->SetXTitle("p_{T trigger} (GeV/c)");
552     
553     fhNeutralLeadingDeltaPt  = new TH2F("NeutralLeadingDeltaPt","p_{T trigger} - p_{T #pi^{0}} vs p_{T trigger}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
554     fhNeutralLeadingDeltaPt->SetYTitle("#Delta p_{T} (GeV/c)");
555     fhNeutralLeadingDeltaPt->SetXTitle("p_{T trigger} (GeV/c)");
556     
557     fhNeutralLeadingDeltaPhi  = new TH2F("NeutralLeadingDeltaPhi","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,TMath::TwoPi()); 
558     fhNeutralLeadingDeltaPhi->SetYTitle("#Delta #phi (rad)");
559     fhNeutralLeadingDeltaPhi->SetXTitle("p_{T trigger} (GeV/c)");
560     
561     fhNeutralLeadingDeltaEta  = new TH2F("NeutralLeadingDeltaEta","#eta_{trigger} - #eta_{#pi^{0}} vs p_{T trigger}",nptbins,ptmin,ptmax,120,-2,2); 
562     fhNeutralLeadingDeltaEta->SetYTitle("#Delta #eta");
563     fhNeutralLeadingDeltaEta->SetXTitle("p_{T trigger} (GeV/c)");
564     
565     fhNeutralLeadingRatioPt  = new TH2F("NeutralLeadingRatioPt","p_{T leading #pi^{0}} /p_{T trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,2); 
566     fhNeutralLeadingRatioPt->SetYTitle("p_{T lead #pi^{0}} /p_{T trigger}");
567     fhNeutralLeadingRatioPt->SetXTitle("p_{T trigger} (GeV/c)");
568     
569     fhNeutralLeadingXi  = new TH2F("NeutralLeadingXi","ln(p_{T trigger} / p_{T leading #pi^{0}} ) vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,10); 
570     fhNeutralLeadingXi->SetYTitle("#xi");
571     fhNeutralLeadingXi->SetXTitle("p_{T trigger} (GeV/c)");
572     
573     fOutCont->Add(fhNeutralLeadingPt) ;
574     fOutCont->Add(fhNeutralLeadingPhi) ;
575     fOutCont->Add(fhNeutralLeadingEta) ;
576     fOutCont->Add(fhNeutralLeadingDeltaPt) ; 
577     fOutCont->Add(fhNeutralLeadingDeltaPhi) ; 
578     fOutCont->Add(fhNeutralLeadingDeltaEta) ; 
579     fOutCont->Add(fhNeutralLeadingRatioPt) ;
580     fOutCont->Add(fhNeutralLeadingXi) ;
581           
582         fhNeutralLeadingDeltaPhiRatioPt30  = new TH2F("NeutralLeadingDeltaPhiRatioPt30","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T leading}/p_{T trigger}, neutral leading, p_{T trigger} > 30 GeV/c",120,0,TMath::TwoPi(),nptbins,0,1); 
583         fhNeutralLeadingDeltaPhiRatioPt30->SetXTitle("#Delta #phi (rad)");
584         fhNeutralLeadingDeltaPhiRatioPt30->SetYTitle("p_{T leading} / p_{T trigger}");
585           
586         fhNeutralLeadingDeltaPhiRatioPt50  = new TH2F("NeutralLeadingDeltaPhiRatioPt50","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T leading}/p_{T trigger}, neutral leading, p_{T trigger} > 50 GeV/c",120,0,TMath::TwoPi(),nptbins,0,1); 
587         fhNeutralLeadingDeltaPhiRatioPt50->SetXTitle("#Delta #phi (rad)");
588         fhNeutralLeadingDeltaPhiRatioPt50->SetYTitle("p_{T leading} / p_{T trigger}");
589         fOutCont->Add(fhNeutralLeadingDeltaPhiRatioPt30) ; 
590         fOutCont->Add(fhNeutralLeadingDeltaPhiRatioPt50) ; 
591
592   }
593   
594   if(!fSeveralConeAndPtCuts){// not several cones
595     
596     //Jet Distributions
597     fhJetPt  = new TH2F("JetPt","p_{T  jet} vs p_{T trigger}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
598     fhJetPt->SetYTitle("p_{T  jet}");
599     fhJetPt->SetXTitle("p_{T trigger} (GeV/c)");
600     
601     fhJetRatioPt  = new TH2F("JetRatioPt","p_{T  jet}/p_{T trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,2); 
602     fhJetRatioPt->SetYTitle("p_{T  jet}/p_{T trigger}");
603     fhJetRatioPt->SetXTitle("p_{T trigger} (GeV/c)");
604     
605     fhJetDeltaPhi  = new TH2F("JetDeltaPhi","#phi_{jet} - #phi_{trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,TMath::TwoPi()); 
606     fhJetDeltaPhi->SetYTitle("#Delta #phi (rad)");
607     fhJetDeltaPhi->SetXTitle("p_{T trigger} (GeV/c)");
608     
609     fhJetDeltaEta  = new TH2F("JetDeltaEta","#eta_{jet} - #eta_{trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,120,-2,2); 
610     fhJetDeltaEta->SetYTitle("#Delta #eta");
611     fhJetDeltaEta->SetXTitle("p_{T trigger} (GeV/c)");
612     
613     fhJetLeadingRatioPt  = new TH2F("JetLeadingRatioPt","p_{T  jet} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,2); 
614     fhJetLeadingRatioPt->SetYTitle("p_{T  leading}/p_{T jet}");
615     fhJetLeadingRatioPt->SetXTitle("p_{T trigger} (GeV/c)");
616     
617     fhJetLeadingDeltaPhi  = new TH2F("JetLeadingDeltaPhi","#phi_{jet} - #phi_{leading} vs p_{T trigger}",nptbins,ptmin,ptmax,120,-TMath::Pi(),TMath::Pi()); 
618     fhJetLeadingDeltaPhi->SetYTitle("#Delta #phi (rad)");
619     fhJetLeadingDeltaPhi->SetXTitle("p_{T trigger} (GeV/c)");
620     
621     fhJetLeadingDeltaEta  = new TH2F("JetLeadingDeltaEta","#eta_{jet} - #eta_{leading} vs p_{T trigger}",nptbins,ptmin,ptmax,120,-2,2); 
622     fhJetLeadingDeltaEta->SetYTitle("#Delta #eta");
623     fhJetLeadingDeltaEta->SetXTitle("p_{T trigger} (GeV/c)");
624     
625     fhJetFFz  = new TH2F("JetFFz","z = p_{T i charged}/p_{T trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,200,0.,2); 
626     fhJetFFz->SetYTitle("z");
627     fhJetFFz->SetXTitle("p_{T trigger}");
628     
629     fhJetFFxi  = new TH2F("JetFFxi","#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger}",nptbins,ptmin,ptmax,100,0.,10.); 
630     fhJetFFxi->SetYTitle("#xi");
631     fhJetFFxi->SetXTitle("p_{T trigger}");
632     
633     fhJetFFpt  = new TH2F("JetFFpt","#xi = p_{T i charged}) vs p_{T trigger}",nptbins,ptmin,ptmax,200,0.,50.); 
634     fhJetFFpt->SetYTitle("p_{T charged hadron}");
635     fhJetFFpt->SetXTitle("p_{T trigger}");
636     
637     fhJetNTracksInCone  = new TH2F("JetNTracksInCone","N particles in cone vs p_{T trigger}",nptbins,ptmin,ptmax,5000,0, 5000); 
638     fhJetNTracksInCone->SetYTitle("N tracks in jet cone");
639     fhJetNTracksInCone->SetXTitle("p_{T trigger} (GeV/c)");
640     
641     fOutCont->Add(fhJetPt) ; 
642     fOutCont->Add(fhJetRatioPt) ; 
643     fOutCont->Add(fhJetDeltaPhi) ;
644     fOutCont->Add(fhJetDeltaEta) ;
645     fOutCont->Add(fhJetLeadingRatioPt) ;
646     fOutCont->Add(fhJetLeadingDeltaPhi) ;
647     fOutCont->Add(fhJetLeadingDeltaEta) ;
648     fOutCont->Add(fhJetFFz) ;
649     fOutCont->Add(fhJetFFxi) ;
650     fOutCont->Add(fhJetFFpt) ;
651     fOutCont->Add(fhJetNTracksInCone) ;
652     
653     //Bkg Distributions
654     fhBkgPt  = new TH2F("BkgPt","p_{T  bkg} vs p_{T trigger}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
655     fhBkgPt->SetYTitle("p_{T  bkg}");
656     fhBkgPt->SetXTitle("p_{T trigger} (GeV/c)");
657     
658     fhBkgRatioPt  = new TH2F("BkgRatioPt","p_{T  bkg}/p_{T trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,2); 
659     fhBkgRatioPt->SetYTitle("p_{T  bkg}/p_{T trigger}");
660     fhBkgRatioPt->SetXTitle("p_{T trigger} (GeV/c)");
661     
662     fhBkgDeltaPhi  = new TH2F("BkgDeltaPhi","#phi_{bkg} - #phi_{trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,TMath::TwoPi()); 
663     fhBkgDeltaPhi->SetYTitle("#Delta #phi (rad)");
664     fhBkgDeltaPhi->SetXTitle("p_{T trigger} (GeV/c)");
665     
666     fhBkgDeltaEta  = new TH2F("BkgDeltaEta","#eta_{bkg} - #eta_{trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,120,-2,2); 
667     fhBkgDeltaEta->SetYTitle("#Delta #eta");
668     fhBkgDeltaEta->SetXTitle("p_{T trigger} (GeV/c)");
669     
670     fhBkgLeadingRatioPt  = new TH2F("BkgLeadingRatioPt","p_{T  bkg} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,2); 
671     fhBkgLeadingRatioPt->SetYTitle("p_{T  leading}/p_{T bkg}");
672     fhBkgLeadingRatioPt->SetXTitle("p_{T trigger} (GeV/c)");
673     
674     fhBkgLeadingDeltaPhi  = new TH2F("BkgLeadingDeltaPhi","#phi_{bkg} - #phi_{leading} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,TMath::TwoPi()); 
675     fhBkgLeadingDeltaPhi->SetYTitle("#Delta #phi (rad)");
676     fhBkgLeadingDeltaPhi->SetXTitle("p_{T trigger} (GeV/c)");
677     
678     fhBkgLeadingDeltaEta  = new TH2F("BkgLeadingDeltaEta","#eta_{bkg} - #eta_{leading} vs p_{T trigger}",nptbins,ptmin,ptmax,120,-2,2); 
679     fhBkgLeadingDeltaEta->SetYTitle("#Delta #eta");
680     fhBkgLeadingDeltaEta->SetXTitle("p_{T trigger} (GeV/c)");
681     
682     fhBkgFFz  = new TH2F("BkgFFz","z = p_{T i charged}/p_{T trigger} vs p_{T trigger}", nptbins,ptmin,ptmax,200,0.,2); 
683     fhBkgFFz->SetYTitle("z");
684     fhBkgFFz->SetXTitle("p_{T trigger}");
685     
686     fhBkgFFxi  = new TH2F("BkgFFxi","#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger}", nptbins,ptmin,ptmax,100,0.,10.); 
687     fhBkgFFxi->SetYTitle("#xi");
688     fhBkgFFxi->SetXTitle("p_{T trigger}");
689     
690     fhBkgFFpt  = new TH2F("BkgFFpt","p_{T charged hadron } vs p_{T trigger}", nptbins,ptmin,ptmax,200,0.,50.); 
691     fhBkgFFpt->SetYTitle("p_{T charged} hadron");
692     fhBkgFFpt->SetXTitle("p_{T trigger}");
693     
694     fhBkgNTracksInCone  = new TH2F("BkgNTracksInCone","N particles in cone vs p_{T trigger}",nptbins,ptmin,ptmax,5000,0, 5000); 
695     fhBkgNTracksInCone->SetYTitle("N tracks in bkg cone");
696     fhBkgNTracksInCone->SetXTitle("p_{T trigger} (GeV/c)");
697     
698     fOutCont->Add(fhBkgPt) ; 
699     fOutCont->Add(fhBkgRatioPt) ; 
700     fOutCont->Add(fhBkgDeltaPhi) ;
701     fOutCont->Add(fhBkgDeltaEta) ;
702     fOutCont->Add(fhBkgLeadingRatioPt) ;
703     fOutCont->Add(fhBkgLeadingDeltaPhi) ;
704     fOutCont->Add(fhBkgLeadingDeltaEta) ;
705     fOutCont->Add(fhBkgFFz) ;
706     fOutCont->Add(fhBkgFFxi) ;
707     fOutCont->Add(fhBkgFFpt) ;
708     fOutCont->Add(fhBkgNTracksInCone) ;
709     
710   }//not several cones
711   else{ //If we want to study the jet for different cones and pt
712     for(Int_t icone = 0; icone<fJetNCone; icone++){//icone
713       for(Int_t ipt = 0; ipt<fJetNPt;ipt++){ //ipt
714         
715         TString lastnamehist ="Cone"+ fJetNameCones[icone]+"Pt"+ fJetNamePtThres[ipt];
716         TString lastnametitle =", cone ="+fJetNameCones[icone]+", pt > " +fJetNamePtThres[ipt]+" GeV/c";
717         
718         //Jet Distributions
719         fhJetPts[icone][ipt] = new TH2F(Form("JetPt%s",lastnamehist.Data()),Form("p_{T  jet} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
720         fhJetPts[icone][ipt]->SetYTitle("p_{T  jet}");
721         fhJetPts[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
722         
723         fhJetRatioPts[icone][ipt] = new TH2F(Form("JetRatioPt%s",lastnamehist.Data()),Form("p_{T  jet}/p_{T trigger} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,0,2); 
724         fhJetRatioPts[icone][ipt]->SetYTitle("p_{T  jet}/p_{T trigger}");
725         fhJetRatioPts[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
726         
727         fhJetDeltaPhis[icone][ipt] = new TH2F(Form("JetDeltaPhi%s",lastnamehist.Data()),Form("#phi_{jet} - #phi_{trigger} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,0,TMath::TwoPi()); 
728         fhJetDeltaPhis[icone][ipt]->SetYTitle("#Delta #phi (rad)");
729         fhJetDeltaPhis[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
730         
731         fhJetDeltaEtas[icone][ipt] = new TH2F(Form("JetDeltaEta%s",lastnamehist.Data()),Form("#eta_{jet} - #eta_{trigger} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,-2,2); 
732         fhJetDeltaEtas[icone][ipt]->SetYTitle("#Delta #eta");
733         fhJetDeltaEtas[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
734         
735         fhJetLeadingRatioPts[icone][ipt] = new TH2F(Form("JetLeadingRatioPt%s",lastnamehist.Data()),Form("p_{T  jet} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,0,2); 
736         fhJetLeadingRatioPts[icone][ipt]->SetYTitle("p_{T  leading}/p_{T jet}");
737         fhJetLeadingRatioPts[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
738         
739         fhJetLeadingDeltaPhis[icone][ipt] = new TH2F(Form("JetLeadingDeltaPhi%s",lastnamehist.Data()),Form("#phi_{jet} - #phi_{leading} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,0,TMath::TwoPi()); 
740         fhJetLeadingDeltaPhis[icone][ipt]->SetYTitle("#Delta #phi (rad)");
741         fhJetLeadingDeltaPhis[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
742         
743         fhJetLeadingDeltaEtas[icone][ipt] = new TH2F(Form("JetLeadingDeltaEta%s",lastnamehist.Data()),Form("#eta_{jet} - #eta_{leading} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,-2,2); 
744         fhJetLeadingDeltaEtas[icone][ipt]->SetYTitle("#Delta #eta");
745         fhJetLeadingDeltaEtas[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
746         
747         fhJetFFzs[icone][ipt] = new TH2F(Form("JetFFz%s",lastnamehist.Data()),"z = p_{T i charged}/p_{T trigger} vs p_{T trigger}", 120,0.,120.,200,0.,2); 
748         fhJetFFzs[icone][ipt]->SetYTitle("z");
749         fhJetFFzs[icone][ipt]->SetXTitle("p_{T trigger}");
750         
751         fhJetFFxis[icone][ipt] = new TH2F(Form("JetFFxi%s",lastnamehist.Data()),"#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger}", 120,0.,120.,100,0.,10.); 
752         fhJetFFxis[icone][ipt]->SetYTitle("#xi");
753         fhJetFFxis[icone][ipt]->SetXTitle("p_{T trigger}");
754         
755         fhJetFFpts[icone][ipt] = new TH2F(Form("JetFFpt%s",lastnamehist.Data()),"p_{T charged hadron } in jet vs p_{T trigger}", 120,0.,120.,200,0.,50.); 
756         fhJetFFpts[icone][ipt]->SetYTitle("p_{T charged hadron}");
757         fhJetFFpts[icone][ipt]->SetXTitle("p_{T trigger}");
758         
759         fhJetNTracksInCones[icone][ipt] = new TH2F(Form("JetNTracksInCone%s",lastnamehist.Data()),Form("N particles in cone vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,5000,0, 5000); 
760         fhJetNTracksInCones[icone][ipt]->SetYTitle("N tracks in jet cone");
761         fhJetNTracksInCones[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
762         
763         fOutCont->Add(fhJetPts[icone][ipt]) ; 
764         fOutCont->Add(fhJetRatioPts[icone][ipt]) ; 
765         fOutCont->Add(fhJetDeltaPhis[icone][ipt]) ;
766         fOutCont->Add(fhJetDeltaEtas[icone][ipt]) ;
767         fOutCont->Add(fhJetLeadingRatioPts[icone][ipt]) ;
768         fOutCont->Add(fhJetLeadingDeltaPhis[icone][ipt]) ;
769         fOutCont->Add(fhJetLeadingDeltaEtas[icone][ipt]) ;
770         fOutCont->Add(fhJetFFzs[icone][ipt]) ;
771         fOutCont->Add(fhJetFFxis[icone][ipt]) ;
772         fOutCont->Add(fhJetFFpts[icone][ipt]) ;
773         fOutCont->Add(fhJetNTracksInCones[icone][ipt]) ;
774         
775         //Bkg Distributions
776         fhBkgPts[icone][ipt] = new TH2F(Form("BkgPt%s",lastnamehist.Data()),Form("p_{T  bkg} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
777         fhBkgPts[icone][ipt]->SetYTitle("p_{T  bkg}");
778         fhBkgPts[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
779         
780         fhBkgRatioPts[icone][ipt] = new TH2F(Form("BkgRatioPt%s",lastnamehist.Data()),Form("p_{T  bkg}/p_{T trigger} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,0,2); 
781         fhBkgRatioPts[icone][ipt]->SetYTitle("p_{T  bkg}/p_{T trigger}");
782         fhBkgRatioPts[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
783         
784         fhBkgDeltaPhis[icone][ipt] = new TH2F(Form("BkgDeltaPhi%s",lastnamehist.Data()),Form("#phi_{bkg} - #phi_{trigger} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,0,TMath::TwoPi()); 
785         fhBkgDeltaPhis[icone][ipt]->SetYTitle("#Delta #phi (rad)");
786         fhBkgDeltaPhis[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
787         
788         fhBkgDeltaEtas[icone][ipt] = new TH2F(Form("BkgDeltaEta%s",lastnamehist.Data()),Form("#eta_{bkg} - #eta_{trigger} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,-2,2); 
789         fhBkgDeltaEtas[icone][ipt]->SetYTitle("#Delta #eta");
790         fhBkgDeltaEtas[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
791         
792         fhBkgLeadingRatioPts[icone][ipt] = new TH2F(Form("BkgLeadingRatioPt%s",lastnamehist.Data()),Form("p_{T  bkg} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,0,2); 
793         fhBkgLeadingRatioPts[icone][ipt]->SetYTitle("p_{T  leading}/p_{T bkg}");
794         fhBkgLeadingRatioPts[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
795         
796         fhBkgLeadingDeltaPhis[icone][ipt] = new TH2F(Form("BkgLeadingDeltaPhi%s",lastnamehist.Data()),Form("#phi_{bkg} - #phi_{leading} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,0,TMath::TwoPi()); 
797         fhBkgLeadingDeltaPhis[icone][ipt]->SetYTitle("#Delta #phi (rad)");
798         fhBkgLeadingDeltaPhis[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
799         
800         fhBkgLeadingDeltaEtas[icone][ipt] = new TH2F(Form("BkgLeadingDeltaEta%s",lastnamehist.Data()),Form("#eta_{bkg} - #eta_{leading} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,-2,2); 
801         fhBkgLeadingDeltaEtas[icone][ipt]->SetYTitle("#Delta #eta");
802         fhBkgLeadingDeltaEtas[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
803         
804         fhBkgFFzs[icone][ipt] = new TH2F(Form("BkgFFz%s",lastnamehist.Data()),"z = p_{T i charged}/p_{T trigger} vs p_{T trigger}", 120,0.,120.,200,0.,2); 
805         fhBkgFFzs[icone][ipt]->SetYTitle("z");
806         fhBkgFFzs[icone][ipt]->SetXTitle("p_{T trigger}");
807         
808         fhBkgFFxis[icone][ipt] = new TH2F(Form("BkgFFxi%s",lastnamehist.Data()),"#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger}", 120,0.,120.,100,0.,10.); 
809         fhBkgFFxis[icone][ipt]->SetYTitle("#xi");
810         fhBkgFFxis[icone][ipt]->SetXTitle("p_{T trigger}");
811         
812         fhBkgFFpts[icone][ipt] = new TH2F(Form("BkgFFpt%s",lastnamehist.Data()),"p_{T charged hadron} in jet vs p_{T trigger}", 120,0.,120.,200,0.,50.); 
813         fhBkgFFpts[icone][ipt]->SetYTitle("p_{T charged hadron}");
814         fhBkgFFpts[icone][ipt]->SetXTitle("p_{T trigger}");
815         
816         fhBkgNTracksInCones[icone][ipt] = new TH2F(Form("BkgNTracksInCone%s",lastnamehist.Data()),Form("N particles in cone vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,5000,0, 5000); 
817         fhBkgNTracksInCones[icone][ipt]->SetYTitle("N tracks in bkg cone");
818         fhBkgNTracksInCones[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
819         
820         fOutCont->Add(fhBkgPts[icone][ipt]) ; 
821         fOutCont->Add(fhBkgRatioPts[icone][ipt]) ; 
822         fOutCont->Add(fhBkgDeltaPhis[icone][ipt]) ;
823         fOutCont->Add(fhBkgDeltaEtas[icone][ipt]) ;
824         fOutCont->Add(fhBkgLeadingRatioPts[icone][ipt]) ;
825         fOutCont->Add(fhBkgLeadingDeltaPhis[icone][ipt]) ;
826         fOutCont->Add(fhBkgLeadingDeltaEtas[icone][ipt]) ;
827         fOutCont->Add(fhBkgFFzs[icone][ipt]) ;
828         fOutCont->Add(fhBkgFFxis[icone][ipt]) ;
829         fOutCont->Add(fhBkgFFpts[icone][ipt]) ;
830         fOutCont->Add(fhBkgNTracksInCones[icone][ipt]) ;
831         
832       }//ipt
833     } //icone
834   }//If we want to study any cone or pt threshold
835   
836   //Keep neutral meson selection histograms if requiered
837   //Setting done in AliNeutralMesonSelection
838   if(GetNeutralMesonSelection()){
839     TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
840     if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
841       for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) fOutCont->Add(nmsHistos->At(i)) ;
842         delete nmsHistos;
843   }
844   
845   
846   if(GetDebug() > 2){
847     printf("AliAnaParticleJetLeadingConeCorrelation::GetCreateOutputObjects() - All histograms names : \n");
848     for(Int_t i  = 0 ;  i<  fOutCont->GetEntries(); i++)
849       printf("Histo i %d name %s\n",i,((fOutCont->At(i))->GetName()));
850     //cout<< (fOutCont->At(i))->GetName()<<endl;
851   }
852   
853   return fOutCont;
854
855 }
856
857 //__________________________________________________________________________________________________________
858 Bool_t  AliAnaParticleJetLeadingConeCorrelation::GetLeadingParticle(AliAODPWG4ParticleCorrelation * particle)
859 {
860   //Search Charged or Neutral leading particle, select the highest one and fill AOD
861   
862   
863   GetLeadingCharge(particle) ;
864   if(!fJetsOnlyInCTS) GetLeadingPi0(particle) ;
865   
866   Double_t ptch = fLeadingChargeMom.Pt();
867   Double_t ptpi = fLeadingPi0Mom   .Pt();
868   
869   if (ptch > 0 || ptpi > 0)
870   {
871     if((ptch >= ptpi))
872     {
873       if(GetDebug() > 1)printf("AliAnaParticleJetLeadingConeCorrelation::GetLeadingParticle() - Leading found in CTS \n");
874       
875       fLeadingMom = fLeadingChargeMom;
876       
877       if(GetDebug() > 1) printf("AliAnaParticleJetLeadingConeCorrelation::GetLeadingParticle() - Found Leading: pt %2.3f, phi %2.3f deg, eta %2.3f\n",
878                                 fLeadingMom.Pt(),fLeadingMom.Phi()*TMath::RadToDeg(),fLeadingMom.Eta()) ;
879       
880       //Put leading in AOD
881       particle->SetLeading(fLeadingChargeMom);
882       particle->SetLeadingDetector("CTS");
883       return kTRUE;
884     }
885     else
886     {
887       if(GetDebug() > 1)
888         printf("AliAnaParticleJetLeadingConeCorrelation::GetLeadingParticle() - Leading found in EMCAL \n");
889       
890       fLeadingMom = fLeadingPi0Mom;
891       
892       if(GetDebug() > 1) printf("AliAnaParticleJetLeadingConeCorrelation::GetLeadingParticle() - Found Leading: pt %2.3f, phi %2.3f, eta %2.3f\n",
893                                 fLeadingMom.Pt(),fLeadingMom.Phi()*TMath::RadToDeg(),fLeadingMom.Eta()) ;
894       //Put leading in AOD
895       particle->SetLeading(fLeadingPi0Mom);
896       particle->SetLeadingDetector("EMCAL");
897       return kTRUE;
898     }
899   }
900   
901   if(GetDebug() > 1)printf ("AliAnaParticleJetLeadingConeCorrelation::GetLeadingParticle() - NO LEADING PARTICLE FOUND \n");
902   
903   return kFALSE;
904   
905 }
906
907 //_______________________________________________________________________________________________________
908 void  AliAnaParticleJetLeadingConeCorrelation::GetLeadingCharge(AliAODPWG4ParticleCorrelation * particle)
909 {
910   //Search for the charged particle with highest pt and with
911   //Phi=Phi_trigger-Pi and pT=0.1E_gamma
912   
913   if(!GetCTSTracks()) return;
914   
915   Double_t ptTrig  = particle->Pt();
916   Double_t phiTrig = particle->Phi();
917   Double_t rat     = -100 ;
918   Double_t ptl     = -100 ;
919   Double_t phil    = -100 ;
920   Double_t pt      = -100.;
921   Double_t phi     = -100.;
922   
923   for(Int_t ipr = 0;ipr < GetCTSTracks()->GetEntriesFast() ; ipr ++ )
924   {
925     AliVTrack* track = (AliVTrack *)(GetCTSTracks()->At(ipr)) ;
926     fTrackVector.SetXYZ(track->Px(),track->Py(),track->Pz());
927     pt   = fTrackVector.Pt();
928     phi  = fTrackVector.Phi() ;
929     if(phi < 0) phi+=TMath::TwoPi();
930     rat  = pt/ptTrig ;
931     
932     //printf("AliAnaParticleJetLeadingConeCorrelation::GetLeadingCharge() - Tracks: pt %2.3f eta %2.3f phi %2.3f pt/ptTrig %2.3f \n",
933     //     pt, fTrackVector.Eta(), phi,pt/ptTrig) ;
934     
935     Float_t deltaphi = TMath::Abs(phiTrig-phi);
936     if((deltaphi > fDeltaPhiMinCut)     && (deltaphi < fDeltaPhiMaxCut) &&
937        (rat      > fLeadingRatioMinCut) && (rat      < fLeadingRatioMaxCut)  && (pt  > ptl))
938     {
939       phil = phi ;
940       ptl  = pt ;
941       fLeadingChargeMom.SetVect(fTrackVector);
942     }
943   }// track loop
944   
945   if(GetDebug() > 1 && ptl > 0 )
946     printf("AliAnaParticleJetLeadingConeCorrelation::GetLeadingCharge() - Leading in CTS: pt %2.3f eta %2.3f phi %2.3f pt/ptTrig %2.3f, |phiTrig-phi| %2.3f \n",
947            ptl, fLeadingChargeMom.Eta(), phil,ptl/ptTrig, TMath::Abs(phiTrig-phil)) ;
948   
949 }
950
951 //____________________________________________________________________________________________________
952 void  AliAnaParticleJetLeadingConeCorrelation::GetLeadingPi0(AliAODPWG4ParticleCorrelation * particle)
953 {
954   //Search for the neutral pion with highest pt and with
955   //Phi=Phi_trigger-Pi and pT=0.1E_gamma
956   
957   if(!GetEMCALClusters()) return ;
958   
959   Double_t ptTrig  = particle->Pt();
960   Double_t phiTrig = particle->Phi();
961   Double_t rat     = -100 ;
962   Double_t ptl     = -100 ;
963   Double_t phil    = -100 ;
964   Double_t pt      = -100.;
965   Double_t phi     = -100.;
966   
967   //Get vertex for photon momentum calculation
968   Double_t vertex [] = {0,0,0} ; //vertex
969   if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
970   {
971     GetVertex(vertex);
972   }
973   
974   //Cluster loop, select pairs with good pt, phi and fill AODs or histograms
975   for(Int_t iclus = 0;iclus < GetEMCALClusters()->GetEntriesFast() ; iclus ++ )
976   {
977     AliVCluster * calo = (AliVCluster *)(GetEMCALClusters()->At(iclus)) ;
978     
979     //Cluster selection, not charged, with photon or pi0 id and in fiducial cut
980     Int_t pdgi=0;
981     if(!SelectCluster(calo, vertex,  fLeadingPhoMom1, pdgi))  continue ;
982     
983     if(GetDebug() > 2) printf("AliAnaParticleJetLeadingConeCorrelation::GetLeadingPi0() - Neutral cluster: pt %2.3f, phi %2.3f \n",
984                               fLeadingPhoMom1.Pt(),fLeadingPhoMom1.Phi());
985     
986     //2 gamma overlapped, found with PID
987     if(pdgi == AliCaloPID::kPi0)
988     {
989       if(GetDebug() > 2) printf("AliAnaParticleJetLeadingConeCorrelation::GetLeadingPi0() - Neutral cluster ID as Pi0 \n");
990       
991       pt  = fLeadingPhoMom1.Pt();
992       rat = pt/ptTrig;
993       phi = fLeadingPhoMom1.Phi();
994       if(phi < 0) phi+=TMath::TwoPi();
995       
996       //Selection within angular and energy limits
997       Float_t deltaphi = TMath::Abs(phiTrig-phi);
998       if(pt > ptl  && rat > fLeadingRatioMinCut  && rat < fLeadingRatioMaxCut  &&
999          deltaphi > fDeltaPhiMinCut && deltaphi < fDeltaPhiMaxCut )
1000       {
1001         phil = phi ;
1002         ptl  = pt ;
1003         fLeadingPi0Mom.SetPxPyPzE(fLeadingPhoMom1.Px(),fLeadingPhoMom1.Py(),fLeadingPhoMom1.Pz(),fLeadingPhoMom1.E());
1004       }// cuts
1005     }// pdg = AliCaloPID::kPi0
1006     //Make invariant mass analysis
1007     else if(pdgi == AliCaloPID::kPhoton)
1008     {
1009       //Search the photon companion in case it comes from  a Pi0 decay
1010       //Apply several cuts to select the good pair
1011       for(Int_t jclus = iclus+1; jclus < GetEMCALClusters()->GetEntriesFast() ; jclus++ )
1012       {
1013         AliVCluster * calo2 = (AliVCluster *) (GetEMCALClusters()->At(jclus)) ;
1014         
1015         //Cluster selection, not charged with photon or pi0 id and in fiducial cut
1016         Int_t pdgj=0;
1017         
1018         if     (!SelectCluster(calo2, vertex,  fLeadingPhoMom2, pdgj))  continue ;
1019         
1020         if(pdgj == AliCaloPID::kPhoton )
1021         {
1022           pt  = (fLeadingPhoMom1+fLeadingPhoMom2).Pt();
1023           phi = (fLeadingPhoMom1+fLeadingPhoMom2).Phi();
1024           if(phi < 0) phi+=TMath::TwoPi();
1025           rat = pt/ptTrig;
1026           
1027           //Selection within angular and energy limits
1028           Float_t deltaphi = TMath::Abs(phiTrig-phi);
1029           
1030           if(GetDebug() > 3 ) printf("AliAnaParticleJetLeadingConeCorrelation::GetLeadingPi0() - Neutral Hadron Correlation: gamma pair: pt %2.2f, phi %2.2f, eta %2.2f, |phiTrig-phi| %2.3f, pt/ptTrig %2.3f, M %2.3f\n",
1031                                      pt,phi,(fLeadingPhoMom1+fLeadingPhoMom2).Eta(), deltaphi, rat, (fLeadingPhoMom1+fLeadingPhoMom2).M());
1032           
1033           if(pt > ptl  && rat > fLeadingRatioMinCut  && rat < fLeadingRatioMaxCut  &&
1034              deltaphi > fDeltaPhiMinCut && deltaphi < fDeltaPhiMaxCut ){
1035             //Select good pair (aperture and invariant mass)
1036             if(GetNeutralMesonSelection()->SelectPair(fLeadingPhoMom1, fLeadingPhoMom2,kEMCAL)){
1037               phil = phi ;
1038               ptl  = pt ;
1039               fLeadingPi0Mom=(fLeadingPhoMom1+fLeadingPhoMom2);
1040               
1041               if(GetDebug() > 3 ) printf("AliAnaParticleJetLeadingConeCorrelation::GetLeadingPi0() - Neutral Hadron Correlation: Selected gamma pair: pt %2.2f, phi %2.2f, eta %2.2f, M %2.3f\n",
1042                                          ptl,phil,(fLeadingPhoMom1+fLeadingPhoMom2).Eta(), (fLeadingPhoMom1+fLeadingPhoMom2).M());
1043             }//pi0 selection
1044             
1045             
1046           }//Pair selected as leading
1047         }//if pair of gammas
1048       }//2nd loop
1049     }// if pdg = 22
1050   }// 1st Loop
1051   
1052   if(GetDebug() > 2 && fLeadingPi0Mom.Pt() > 0 )
1053     printf("AliAnaParticleJetLeadingConeCorrelation::GetLeadingPi0() - Leading EMCAL: pt %2.3f eta %2.3f phi %2.3f pt/Eg %2.3f \n",
1054            fLeadingPi0Mom.Pt(), fLeadingPi0Mom.Eta(), phil,  fLeadingPi0Mom.Pt()/ptTrig) ;
1055   
1056 }
1057
1058 //____________________________________________________________________________
1059 void AliAnaParticleJetLeadingConeCorrelation::InitParameters()
1060 {
1061   //Initialize the parameters of the analysis.
1062   SetInputAODName("PWG4Particle");
1063   SetAODObjArrayName("JetLeadingCone");    
1064   AddToHistogramsName("AnaJetLCCorr_");
1065   
1066   fJetsOnlyInCTS      = kFALSE ;
1067   fPbPb               = kFALSE ;
1068   fReMakeJet          = kFALSE ;
1069   
1070   //Leading selection parameters
1071   fDeltaPhiMinCut     = 2.9 ;
1072   fDeltaPhiMaxCut     = 3.4 ; 
1073   fLeadingRatioMinCut = 0.1;
1074   fLeadingRatioMaxCut = 1.5; 
1075   
1076   //Jet selection parameters
1077   //Fixed cut   
1078   fJetRatioMaxCut     = 1.2 ; 
1079   fJetRatioMinCut     = 0.3 ; 
1080   fJetCTSRatioMaxCut  = 1.2 ;
1081   fJetCTSRatioMinCut  = 0.3 ;
1082   fSelect               = 0  ; //0, Accept all jets, 1, selection depends on energy, 2 fixed selection
1083   
1084   fSelectIsolated = kFALSE;
1085   
1086   //Cut depending on gamma energy
1087   fPtTriggerSelectionCut = 10.; //For Low pt jets+BKG, another limits applied
1088   //Reconstructed jet energy dependence parameters 
1089   //e_jet = a1+e_gamma b2. 
1090   //Index 0-> Pt>2 GeV r = 0.3; Index 1-> Pt>0.5 GeV r = 0.3
1091   fJetE1[0] = -5.75; fJetE1[1] = -4.1;
1092   fJetE2[0] = 1.005; fJetE2[1] = 1.05;
1093   
1094   //Reconstructed sigma of jet energy dependence parameters 
1095   //s_jet = a1+e_gamma b2. 
1096   //Index 0-> Pt>2 GeV r = 0.3; Index 1-> Pt>0.5 GeV r = 0.3
1097   fJetSigma1[0] = 2.65;   fJetSigma1[1] = 2.75;
1098   fJetSigma2[0] = 0.0018; fJetSigma2[1] = 0.033;
1099   
1100   //Background mean energy and RMS
1101   //Index 0-> No BKG; Index 1-> BKG > 2 GeV; 
1102   //Index 2-> (low pt jets)BKG > 0.5 GeV;
1103   //Index > 2, same for CTS conf
1104   fBkgMean[0] = 0.; fBkgMean[1] = 8.8 ; fBkgMean[2] = 69.5;
1105   fBkgMean[3] = 0.; fBkgMean[4] = 6.4;  fBkgMean[5] = 48.6;
1106   fBkgRMS[0]  = 0.; fBkgRMS[1]  = 7.5;  fBkgRMS[2]  = 22.0; 
1107   fBkgRMS[3]  = 0.; fBkgRMS[4]  = 5.4;  fBkgRMS[5]  = 13.2; 
1108
1109   //Factor x of min/max = E -+ x * sigma. Obtained after selecting the
1110   //limits for monoenergetic jets.
1111   //Index 0-> No BKG; Index 1-> BKG > 2 GeV; 
1112   //Index 2-> (low pt jets) BKG > 0.5 GeV;
1113   //Index > 2, same for CTS conf
1114
1115   fJetXMin1[0] =-0.69 ; fJetXMin1[1] = 0.39 ; fJetXMin1[2] =-0.88 ; 
1116   fJetXMin1[3] =-2.0  ; fJetXMin1[4] =-0.442 ; fJetXMin1[5] =-1.1  ;
1117   fJetXMin2[0] = 0.066; fJetXMin2[1] = 0.038; fJetXMin2[2] = 0.034; 
1118   fJetXMin2[3] = 0.25 ; fJetXMin2[4] = 0.113; fJetXMin2[5] = 0.077 ;
1119   fJetXMax1[0] =-3.8  ; fJetXMax1[1] =-0.76 ; fJetXMax1[2] =-3.6  ; 
1120   fJetXMax1[3] =-2.7  ; fJetXMax1[4] =-1.21 ; fJetXMax1[5] =-3.7  ;
1121   fJetXMax2[0] =-0.012; fJetXMax2[1] =-0.022; fJetXMax2[2] = 0.016; 
1122   fJetXMax2[3] =-0.024; fJetXMax2[4] =-0.008; fJetXMax2[5] = 0.027;
1123
1124
1125   //Different cones and pt thresholds to construct the jet
1126
1127   fJetCone        = 0.3  ;
1128   fJetPtThreshold = 0.5   ;
1129   fJetPtThresPbPb = 2.   ;
1130   fJetNCone       = 4    ;
1131   fJetNPt         = 4    ;
1132   fJetCones[0]    = 0.2  ; fJetNameCones[0]   = "02" ;
1133   fJetCones[1]    = 0.3  ; fJetNameCones[1]   = "03" ;
1134   fJetCones[2]    = 0.4  ; fJetNameCones[2]   = "04" ;
1135   fJetCones[2]    = 0.5  ; fJetNameCones[2]   = "05" ;
1136
1137   fJetPtThres[0]  = 0.0  ; fJetNamePtThres[0] = "00" ;
1138   fJetPtThres[1]  = 0.5  ; fJetNamePtThres[1] = "05" ;
1139   fJetPtThres[2]  = 1.0  ; fJetNamePtThres[2] = "10" ;
1140   fJetPtThres[3]  = 2.0  ; fJetNamePtThres[3] = "20" ;
1141 }
1142
1143 //__________________________________________________________________________-
1144 Bool_t AliAnaParticleJetLeadingConeCorrelation::IsJetSelected(const Double_t ptTrig, const Double_t ptjet) const {
1145   //Given the pt of the jet and the trigger particle, select the jet or not
1146   //3 options, fSelect=0 accepts all, fSelect=1 selects jets depending on a 
1147   //function energy dependent and fSelect=2 selects on simple fixed cuts
1148   
1149   if(ptjet == 0) return kFALSE;
1150   
1151   Double_t rat = ptTrig / ptjet ;
1152   
1153   //###############################################################
1154   if(fSelect == 0)
1155     return kTRUE; //Accept all jets, no restriction
1156   //###############################################################
1157   else if(fSelect == 1){
1158     //Check if the energy of the reconstructed jet is within an energy window
1159     //WARNING: to be rechecked, don't remember what all the steps mean
1160     Double_t par[6];
1161     Double_t xmax[2];
1162     Double_t xmin[2];
1163     
1164     Int_t iCTS = 0;
1165     if(fJetsOnlyInCTS)
1166       iCTS = 3 ;
1167     
1168     if(!fPbPb){
1169       //Phythia alone, jets with pt_th > 0.2, r = 0.3 
1170       par[0] = fJetE1[0]; par[1] = fJetE2[0]; 
1171       //Energy of the jet peak
1172       //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
1173       par[2] = fJetSigma1[0]; par[3] = fJetSigma2[0];
1174       //Sigma  of the jet peak
1175       //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
1176       par[4] = fBkgMean[0 + iCTS]; par[5] = fBkgRMS[0 + iCTS];
1177       //Parameters reserved for PbPb bkg.
1178       xmax[0] = fJetXMax1[0 + iCTS]; xmax[1] = fJetXMax2[0 + iCTS];
1179       xmin[0] = fJetXMin1[0 + iCTS]; xmin[1] = fJetXMin2[0 + iCTS];
1180       //Factor that multiplies sigma to obtain the best limits, 
1181       //by observation, of mono jet ratios (ptjet/ptTrig)
1182       //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
1183       
1184     }
1185     else{
1186       if(ptTrig > fPtTriggerSelectionCut){
1187         //Phythia +PbPb with  pt_th > 2 GeV/c, r = 0.3 
1188         par[0] = fJetE1[0]; par[1] = fJetE2[0]; 
1189         //Energy of the jet peak, same as in pp
1190         //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
1191         par[2] = fJetSigma1[0]; par[3] = fJetSigma2[0];
1192         //Sigma  of the jet peak, same as in pp
1193         //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
1194         par[4] = fBkgMean[1 + iCTS]; par[5] = fBkgRMS[1 + iCTS];
1195         //Mean value and RMS of PbPb Bkg 
1196         xmax[0] = fJetXMax1[1 + iCTS]; xmax[1] = fJetXMax2[1 + iCTS];
1197         xmin[0] = fJetXMin1[1 + iCTS]; xmin[1] = fJetXMin2[1 + iCTS];
1198         //Factor that multiplies sigma to obtain the best limits, 
1199         //by observation, of mono jet ratios (ptjet/ptTrig) mixed with PbPb Bkg, 
1200         //pt_th > 2 GeV, r = 0.3
1201         //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
1202         
1203       }
1204       else{
1205         //Phythia + PbPb with  pt_th > 0.5 GeV/c, r = 0.3
1206         par[0] = fJetE1[1]; par[1] = fJetE2[1]; 
1207         //Energy of the jet peak, pt_th > 2 GeV/c, r = 0.3 
1208         //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
1209         par[2] = fJetSigma1[1]; par[3] = fJetSigma2[1];
1210         //Sigma  of the jet peak, pt_th > 2 GeV/c, r = 0.3
1211         //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
1212         par[4] = fBkgMean[2 + iCTS]; par[5] = fBkgRMS[2 + iCTS];
1213         //Mean value and RMS of PbPb Bkg in a 0.3 cone, pt > 2 GeV.
1214         xmax[0] = fJetXMax1[2 + iCTS]; xmax[1] = fJetXMax2[2 + iCTS];
1215         xmin[0] = fJetXMin1[2 + iCTS]; xmin[1] = fJetXMin2[2 + iCTS];
1216         //Factor that multiplies sigma to obtain the best limits, 
1217         //by observation, of mono jet ratios (ptjet/ptTrig) mixed with PbPb Bkg, 
1218         //pt_th > 2 GeV, r = 0.3
1219         //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
1220         
1221       }//If low pt jet in bkg
1222     }//if Bkg
1223     
1224     //Calculate minimum and maximum limits of the jet ratio.
1225     Double_t min = CalculateJetRatioLimit(ptTrig, par, xmin);
1226     Double_t max = CalculateJetRatioLimit(ptTrig, par, xmax);
1227     
1228     if(GetDebug() > 3)printf("AliAnaParticleJetLeadingConeCorrelation::IsJetSelected() - Jet selection?  : Limits min %2.3f, max %2.3f,  pt_jet %2.3f,  pt_gamma %2.3f, pt_jet / pt_gamma %2.3f\n",min,max,ptjet,ptTrig,rat);
1229     
1230     if(( min < rat ) && ( max > ptjet/rat))
1231       return kTRUE;
1232     else
1233       return kFALSE;
1234   }//fSelect = 1
1235   //###############################################################
1236   else if(fSelect == 2){
1237     //Simple selection
1238     if(!fJetsOnlyInCTS){
1239       if((rat <  fJetRatioMaxCut) && (rat > fJetRatioMinCut )) return kTRUE;
1240     }
1241     else{
1242       if((rat <  fJetCTSRatioMaxCut) && (rat > fJetCTSRatioMinCut )) return kTRUE;
1243     }
1244   }// fSelect = 2
1245   //###############################################################
1246   else{
1247     printf("AliAnaParticleJetLeadingConeCorrelation::IsJetSelected() - Jet selection option larger than 2, DON'T SELECT JETS\n");
1248     return kFALSE ;
1249   }
1250   
1251   return kFALSE;
1252   
1253 }
1254
1255 //___________________________________________________________________
1256 Bool_t AliAnaParticleJetLeadingConeCorrelation::IsParticleInJetCone(const Double_t eta, Double_t phi, const Double_t etal,  Double_t phil)
1257   const {
1258   //Check if the particle is inside the cone defined by the leading particle
1259   //WARNING: To be rechecked
1260   
1261   if(phi < 0) phi+=TMath::TwoPi();
1262   if(phil < 0) phil+=TMath::TwoPi();  
1263   Double_t  rad = 10000 + fJetCone;
1264   
1265   if(TMath::Abs(phi-phil) <= (TMath::TwoPi() - fJetCone))
1266     rad = TMath::Sqrt(TMath::Power(eta-etal,2)+TMath::Power(phi-phil,2));
1267   else{
1268     if(phi-phil > TMath::TwoPi() - fJetCone)
1269       rad = TMath::Sqrt(TMath::Power(eta-etal,2)+TMath::Power((phi-TMath::TwoPi())-phil,2));
1270     if(phi-phil < -(TMath::TwoPi() - fJetCone))
1271       rad = TMath::Sqrt(TMath::Power(eta-etal,2)+TMath::Power((phi+TMath::TwoPi())-phil,2));
1272   }
1273   
1274   if(rad < fJetCone) return kTRUE ;
1275   else return kFALSE ;
1276   
1277 }
1278
1279 //__________________________________________________________________
1280 void  AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillAOD() 
1281 {  
1282   //Particle-Hadron Correlation Analysis, fill AODs
1283   
1284   if(!GetInputAODBranch()){
1285     printf("AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillAOD() - No input particles in AOD with name branch < %s > \n",
1286            GetInputAODName().Data());
1287     abort();
1288   }     
1289   
1290   if(strcmp(GetInputAODBranch()->GetClass()->GetName(), "AliAODPWG4ParticleCorrelation")){
1291     printf("AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillAOD() - Wrong type of AOD object, change AOD class name in input AOD: It should be <AliAODPWG4ParticleCorrelation> and not <%s> \n",GetInputAODBranch()->GetClass()->GetName());
1292     abort();
1293   }     
1294   
1295   if(GetDebug() > 1){
1296     printf("AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillAOD() - Begin jet leading cone  correlation analysis, fill AODs \n");
1297     printf("AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillAOD() - In particle branch aod entries %d\n", GetInputAODBranch()->GetEntriesFast());
1298     printf("AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillAOD() - In CTS aod entries %d\n", GetCTSTracks()->GetEntriesFast());
1299     printf("AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillAOD() - In EMCAL aod entries %d\n", GetEMCALClusters()->GetEntriesFast());
1300   }
1301   
1302   
1303   //Loop on stored AOD particles, trigger
1304   Int_t naod = GetInputAODBranch()->GetEntriesFast();
1305   for(Int_t iaod = 0; iaod < naod ; iaod++)
1306   {
1307     AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));                  
1308     
1309     //  printf("AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillAOD() - Trigger : pt %3.2f, phi %2.2f, eta %2.2f\n",particle->Pt(), particle->Phi(), particle->Eta());
1310     
1311     //Search leading particles in CTS and EMCAL 
1312     if(GetLeadingParticle(particle))
1313     {
1314       //Construct the jet around the leading, Fill AOD jet particle list, select jet 
1315       //and fill AOD with jet and background     
1316       MakeAODJet(particle);
1317       
1318     }//Leading found
1319   }//AOD trigger particle loop
1320   
1321   if(GetDebug() >1)printf("AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillAOD() - End of jet leading cone analysis, fill AODs \n");
1322   
1323
1324
1325 //_________________________________________________________________________
1326 void  AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillHistograms()
1327 {
1328   
1329   //Particle-Hadron Correlation Analysis, fill histograms
1330   
1331   if(!GetInputAODBranch())
1332   {
1333     printf("AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillHistograms() - No input particles in AOD with name branch < %s > \n",
1334            GetInputAODName().Data());
1335     abort();
1336   }
1337   
1338   if(GetDebug() > 1)
1339   {
1340     printf("AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillHistograms() - Begin jet leading cone  correlation analysis, fill histograms \n");
1341     printf("AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillHistograms() - In particle branch aod entries %d\n", GetInputAODBranch()->GetEntriesFast());
1342     printf("AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillHistograms() - In CTS aod entries %d\n", GetCTSTracks()->GetEntriesFast());
1343     printf("AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillHistograms() - In EMCAL aod entries %d\n", GetEMCALClusters()->GetEntriesFast());
1344   }
1345   
1346   //Loop on stored AOD particles, trigger
1347   Int_t naod = GetInputAODBranch()->GetEntriesFast();
1348   for(Int_t iaod = 0; iaod < naod ; iaod++)
1349   {
1350     AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1351     
1352     if(OnlyIsolated() && !particle->IsIsolated()) continue;
1353     
1354     Double_t pt  = particle->Pt();
1355     Double_t phi = particle->Phi();
1356     Double_t eta = particle->Eta();
1357     
1358     //Get leading particle, fill histograms
1359     fLeadingMom = particle->GetLeading();
1360     TString det = particle->GetLeadingDetector();
1361     
1362     if(det!="" && fLeadingMom.Pt() > 0)
1363     {
1364       Double_t ptL = fLeadingMom.Pt();
1365       Double_t phiL = fLeadingMom.Phi();
1366       if(phiL < 0 ) phiL+=TMath::TwoPi();
1367       Double_t etaL = fLeadingMom.Eta();
1368       
1369       if(GetDebug() > 1)   printf("AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillHistograms() - Trigger with pt %3.2f, phi %2.2f, eta %2.2f\n", pt, phi, eta);
1370       
1371       if(det == "CTS"){
1372         fhChargedLeadingPt->Fill(pt,ptL);
1373         fhChargedLeadingPhi->Fill(pt,phiL);
1374         fhChargedLeadingEta->Fill(pt,etaL);
1375         fhChargedLeadingDeltaPt->Fill(pt,pt-ptL);
1376         fhChargedLeadingDeltaPhi->Fill(pt,TMath::Abs(phi-phiL));
1377         fhChargedLeadingDeltaEta->Fill(pt,eta-etaL);
1378         fhChargedLeadingRatioPt->Fill(pt,ptL/pt);
1379         fhChargedLeadingXi->Fill(pt,TMath::Log(pt/ptL));
1380         if(pt > 30) fhChargedLeadingDeltaPhiRatioPt30->Fill(TMath::Abs(phi-phiL),ptL/pt);
1381         if(pt > 50) fhChargedLeadingDeltaPhiRatioPt50->Fill(TMath::Abs(phi-phiL),ptL/pt);
1382       }
1383       else if(det== "EMCAL"){
1384         fhNeutralLeadingPt->Fill(pt,ptL);
1385         fhNeutralLeadingPhi->Fill(pt,phiL);
1386         fhNeutralLeadingEta->Fill(pt,etaL);
1387         fhNeutralLeadingDeltaPt->Fill(pt,pt-ptL);
1388         fhNeutralLeadingDeltaPhi->Fill(pt,TMath::Abs(phi-phiL));
1389         fhNeutralLeadingDeltaEta->Fill(pt,eta-etaL);
1390         fhNeutralLeadingRatioPt->Fill(pt,ptL/pt);
1391         fhNeutralLeadingXi->Fill(pt,TMath::Log(pt/ptL));
1392         if(pt > 30) fhNeutralLeadingDeltaPhiRatioPt30->Fill(TMath::Abs(phi-phiL),ptL/pt);
1393         if(pt > 50) fhNeutralLeadingDeltaPhiRatioPt50->Fill(TMath::Abs(phi-phiL),ptL/pt);
1394         
1395       }
1396       
1397       //Fill Jet histograms
1398       
1399       if(!fSeveralConeAndPtCuts)
1400       {//just fill histograms
1401         if(!fReMakeJet)
1402         {
1403           fJetMom=particle->GetCorrelatedJet();
1404           fBkgMom=particle->GetCorrelatedBackground();
1405         }
1406         else  MakeJetFromAOD(particle);
1407         
1408         if(fJetMom.Pt() > 0)
1409         {//Jet was found
1410           FillJetHistos(particle, fJetMom,"Jet","");
1411           FillJetHistos(particle, fBkgMom,"Bkg","");
1412         }
1413       }
1414       else if(fSeveralConeAndPtCuts)
1415       {
1416         for(Int_t icone = 0; icone<fJetNCone; icone++)
1417         {
1418           fJetCone=fJetCones[icone];
1419           for(Int_t ipt = 0; ipt<fJetNPt;ipt++)
1420           {
1421             TString lastname ="Cone"+ fJetNameCones[icone]+"Pt"+ fJetNamePtThres[ipt];
1422             fJetPtThreshold=fJetPtThres[ipt];
1423             
1424             MakeJetFromAOD(particle);
1425             
1426             if(fJetMom.Pt() > 0)
1427             {//Jet was found
1428               FillJetHistos(particle, fJetMom,"Jet",lastname);
1429               FillJetHistos(particle, fBkgMom,"Bkg",lastname);
1430             }
1431           }//icone
1432         }//ipt
1433       }//fSeveralConeAndPtCuts
1434     }//Leading
1435   }//AOD trigger particle loop
1436   
1437   if(GetDebug() >1)printf("AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillHistograms() - End of jet leading cone analysis, fill histograms \n");
1438   
1439
1440
1441 //_______________________________________________________________________________________________
1442 void AliAnaParticleJetLeadingConeCorrelation::MakeAODJet(AliAODPWG4ParticleCorrelation *particle)
1443 {
1444   //Fill the jet with the particles around the leading particle with
1445   //R=fJetCone and pt_th = fJetPtThres. Calculate the energy of the jet and
1446   //fill aod with found information
1447   
1448   Double_t ptTrig   = particle->Pt();
1449   Double_t phiTrig  = particle->Phi();
1450   Double_t phil     = fLeadingMom.Phi();
1451   if(phil<0) phil+=TMath::TwoPi();
1452   Double_t etal     = fLeadingMom.Eta();
1453   
1454   //Different pt cut for jet particles in different collisions systems
1455   Float_t ptcut = fJetPtThreshold;
1456   if(fPbPb && !fSeveralConeAndPtCuts && ptTrig > fPtTriggerSelectionCut)  ptcut = fJetPtThresPbPb ;
1457   
1458   //Add charged particles to jet if they are in cone around the leading particle
1459   if(!GetCTSTracks()) {
1460     printf("AliAnaParticleJetLeadingConeCorrelation::MakeAODJet() - Cannot construct jets without tracks, STOP analysis");
1461     return;
1462   }
1463   
1464   //Fill jet with tracks
1465   //Initialize reference arrays that will contain jet and background tracks
1466   TObjArray * reftracks     = new TObjArray;
1467   TObjArray * reftracksbkg  = new TObjArray;
1468   
1469   for(Int_t ipr = 0;ipr < (GetCTSTracks())->GetEntriesFast() ; ipr ++ ){
1470     AliVTrack* track = (AliVTrack *)((GetCTSTracks())->At(ipr)) ;
1471     fTrackVector.SetXYZ(track->Px(),track->Py(),track->Pz());
1472     
1473     //Particles in jet
1474     if(IsParticleInJetCone(fTrackVector.Eta(), fTrackVector.Phi(), etal, phil)){
1475       
1476       reftracks->Add(track);
1477       
1478       if(fTrackVector.Pt() > ptcut )
1479       {
1480         fJetConstMom.SetVect(fTrackVector);
1481         fJetMom+=fJetConstMom;
1482       }
1483     }
1484     
1485     //Background around (phi_gamma-pi, eta_leading)
1486     else if(IsParticleInJetCone(fTrackVector.Eta(),fTrackVector.Phi(),etal, phiTrig)) {
1487       
1488       reftracksbkg->Add(track);
1489       
1490       if(fTrackVector.Pt() > ptcut ){
1491         fJetConstMom.SetVect(fTrackVector);
1492         fBkgMom+=fJetConstMom;
1493       }
1494     }
1495   }//Track loop
1496   
1497   //Add referenced tracks to AOD
1498   if(reftracks->GetEntriesFast() > 0 )
1499   {
1500     reftracks->SetName(Form("%sTracks",GetAODObjArrayName().Data()));
1501     particle->AddObjArray(reftracks);
1502   }
1503   else  if(GetDebug() > 2 ) printf("AliAnaParticleJetLeadingConeCorrelation::MakeAODJet() - No tracks in jet cone\n");
1504   if(reftracksbkg->GetEntriesFast() > 0 )
1505   {
1506     reftracksbkg->SetName(Form("%sTracksBkg",GetAODObjArrayName().Data()));
1507     particle->AddObjArray(reftracksbkg);
1508   }
1509   else  if(GetDebug() > 2 ) printf("AliAnaParticleJetLeadingConeCorrelation::MakeAODJet() - No background tracks in jet cone\n");
1510   
1511   //Add neutral particles to jet
1512   //Initialize reference arrays that will contain jet and background tracks
1513   TObjArray * refclusters     = new TObjArray;
1514   TObjArray * refclustersbkg  = new TObjArray;
1515   if(!fJetsOnlyInCTS && GetEMCALClusters()){
1516     
1517     //Get vertex for photon momentum calculation
1518     Double_t vertex[]  = {0,0,0} ; //vertex
1519     if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC)
1520     {
1521       GetReader()->GetVertex(vertex);
1522       //if(GetReader()->GetSecondInputAODTree()) GetReader()->GetSecondInputAODVertex(vertex2);
1523     }
1524     
1525     for(Int_t iclus = 0;iclus < (GetEMCALClusters())->GetEntriesFast() ; iclus ++ )
1526     {
1527       AliVCluster * calo = (AliVCluster *) (GetEMCALClusters()->At(iclus)) ;
1528       
1529       //Cluster selection, not charged
1530       if(IsTrackMatched(calo,GetReader()->GetInputEvent())) continue ;
1531       
1532       //Get Momentum vector,
1533       calo->GetMomentum(fJetConstMom,vertex) ;//Assume that come from vertex in straight line
1534       
1535       //Particles in jet
1536       if(IsParticleInJetCone(fJetConstMom.Eta(),fJetConstMom.Phi(), etal, phil)){
1537         
1538         refclusters->Add(calo);
1539         
1540         if(fJetConstMom.Pt() > ptcut )  fJetMom+=fJetConstMom;
1541       }
1542       //Background around (phi_gamma-pi, eta_leading)
1543       else if(IsParticleInJetCone(fJetConstMom.Eta(),fJetConstMom.Phi(),etal, phiTrig)){
1544         
1545         
1546         refclustersbkg->Add(calo);
1547         
1548         if(fJetConstMom.Pt() > ptcut )  fBkgMom+=fJetConstMom;
1549       }
1550     }//cluster loop
1551   }//jets with neutral particles
1552   
1553   //Add referenced clusters to AOD
1554   if(refclusters->GetEntriesFast() > 0 ){
1555     refclusters->SetName(Form("%sClusters",GetAODObjArrayName().Data()));
1556     particle->AddObjArray(refclusters);
1557   }
1558   else  if(GetDebug() > 2 ) printf("AliAnaParticleJetLeadingConeCorrelation::MakeAODJet() - No clusters in jet cone\n");
1559   if(refclustersbkg->GetEntriesFast() > 0 )
1560   {
1561     refclustersbkg->SetName(Form("%sClustersBkg",GetAODObjArrayName().Data()));
1562     particle->AddObjArray(refclustersbkg);
1563   }
1564   else if(GetDebug() > 2 ) printf("AliAnaParticleJetLeadingConeCorrelation::MakeAODJet() - No background clusters in jet cone\n");
1565   
1566   //If there is any jet found, select after some criteria and
1567   //and fill AOD with corresponding TLorentzVector kinematics
1568   if(IsJetSelected(particle->Pt(), fJetMom.Pt()))
1569   {
1570     particle->SetCorrelatedJet(fJetMom);
1571     particle->SetCorrelatedBackground(fBkgMom);
1572     if(GetDebug()>1) printf("AliAnaParticleJetLeadingConeCorrelation::MakeAODJet() - Found jet: Trigger  pt %2.3f, Jet pt %2.3f, Bkg pt %2.3f\n",ptTrig,fJetMom.Pt(),fBkgMom.Pt());
1573   }  
1574   
1575 }
1576
1577 //______________________________________________________________________________________________________
1578 void AliAnaParticleJetLeadingConeCorrelation::MakeJetFromAOD(AliAODPWG4ParticleCorrelation *particle)
1579 {
1580   //Fill the jet with the particles around the leading particle with
1581   //R=fJetCone and pt_th = fJetPtThres. Calculate the energy of the jet and
1582   //fill aod tlorentzvectors with jet and bakcground found
1583   
1584   Double_t ptTrig  = particle->Pt();
1585   Double_t phiTrig = particle->Phi();
1586   Double_t etal    = fLeadingMom.Eta();
1587   Double_t phil    = fLeadingMom.Phi();
1588   if(phil < 0) phil+=TMath::TwoPi();
1589   
1590   TObjArray * refclusters    = particle->GetObjArray(Form("Clusters%s"   ,GetAODObjArrayName().Data()));
1591   TObjArray * reftracks      = particle->GetObjArray(Form("Tracks%s"     ,GetAODObjArrayName().Data()));
1592   TObjArray * refclustersbkg = particle->GetObjArray(Form("ClustersBkg%s",GetAODObjArrayName().Data()));
1593   TObjArray * reftracksbkg   = particle->GetObjArray(Form("TracksBkg%s"  ,GetAODObjArrayName().Data()));
1594   
1595   //Different pt cut for jet particles in different collisions systems
1596   Float_t ptcut = fJetPtThreshold;
1597   if(fPbPb && !fSeveralConeAndPtCuts && ptTrig > fPtTriggerSelectionCut)  ptcut = fJetPtThresPbPb ;
1598   
1599   //Fill jet with tracks
1600   //Particles in jet
1601   if(reftracks)
1602   {
1603     for(Int_t ipr = 0;ipr < reftracks->GetEntriesFast() ; ipr ++ )
1604     {
1605       AliVTrack* track = (AliVTrack *) reftracks->At(ipr) ;
1606       fTrackVector.SetXYZ(track->Px(),track->Py(),track->Pz());
1607       Float_t phi = fTrackVector.Phi();
1608       if(phi < 0) phi+=TMath::TwoPi();
1609       if(fTrackVector.Pt() > ptcut && IsParticleInJetCone(fTrackVector.Eta(), phi, etal, phil) )
1610       {
1611         fJetConstMom.SetVect(fTrackVector);
1612         fJetMom+=fJetConstMom;
1613       }
1614     }//jet Track loop
1615   }
1616   //Particles in background
1617   if(reftracksbkg){
1618     for(Int_t ipr = 0;ipr < reftracksbkg->GetEntriesFast() ; ipr ++ )
1619     {
1620       AliVTrack* track = (AliVTrack *) reftracksbkg->At(ipr) ;
1621       fTrackVector.SetXYZ(track->Px(),track->Py(),track->Pz());
1622       if(fTrackVector.Pt() > ptcut && IsParticleInJetCone(fTrackVector.Eta(),fTrackVector.Phi(),etal, phiTrig) )
1623       {
1624         fJetConstMom.SetVect(fTrackVector);
1625         fBkgMom+=fJetConstMom;
1626       }
1627     }//background Track loop
1628   }
1629   
1630   //Add neutral particles to jet
1631   if(!fJetsOnlyInCTS && refclusters)
1632   {
1633     //Get vertex for photon momentum calculation
1634     Double_t vertex[]  = {0,0,0} ; //vertex
1635     if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC)
1636     {
1637       GetReader()->GetVertex(vertex);
1638     }
1639     
1640     //Loop on jet particles
1641     if(refclusters){
1642       for(Int_t iclus = 0;iclus < refclusters->GetEntriesFast() ; iclus ++ )
1643       {
1644         AliVCluster * calo = (AliVCluster *) refclusters->At(iclus) ;
1645         
1646         calo->GetMomentum(fJetConstMom,vertex) ;//Assume that come from vertex in straight line
1647         
1648         if(fJetConstMom.Pt() > ptcut && IsParticleInJetCone(fJetConstMom.Eta(),fJetConstMom.Phi(), etal, phil)) fJetMom+=fJetConstMom;
1649       }//jet cluster loop
1650     }
1651     
1652     //Loop on background particles
1653     if(refclustersbkg)
1654     {
1655       for(Int_t iclus = 0;iclus < refclustersbkg->GetEntriesFast() ; iclus ++ )
1656       {
1657         AliVCluster * calo = (AliVCluster *) refclustersbkg->At(iclus) ;
1658         
1659         calo->GetMomentum(fJetConstMom,vertex) ;//Assume that come from vertex in straight line
1660         
1661         if( fJetConstMom.Pt() > ptcut && IsParticleInJetCone(fJetConstMom.Eta(),fJetConstMom.Phi(),etal, phiTrig)) fBkgMom+=fJetConstMom;
1662       }//background cluster loop
1663     }
1664   }//clusters in jet
1665   
1666   //If there is any jet found, leave jet and bkg as they are,
1667   //if not set them to 0.
1668   if(!IsJetSelected(particle->Pt(), fJetMom.Pt()))
1669   {
1670     fJetMom.SetPxPyPzE(0.,0.,0.,0.);
1671     fBkgMom.SetPxPyPzE(0.,0.,0.,0.);
1672   }
1673   else
1674     if(GetDebug()>1) printf("AliAnaParticleJetLeadingConeCorrelation::MakeJetFromAOD() - Found jet: Trigger  pt %2.3f, Jet pt %2.3f, Bkg pt %2.3f\n",
1675                             ptTrig,fJetMom.Pt(),fBkgMom.Pt());
1676   
1677 }
1678
1679 //__________________________________________________________________
1680 void AliAnaParticleJetLeadingConeCorrelation::Print(const Option_t * opt) const
1681 {
1682   
1683   //Print some relevant parameters set for the analysis
1684   if(! opt)
1685     return;
1686   
1687   printf("**** Print  %s %s ****\n", GetName(), GetTitle() ) ;
1688   AliAnaCaloTrackCorrBaseClass::Print(" ");
1689   
1690   if(fJetsOnlyInCTS)printf("Jets reconstructed in CTS \n");
1691   else printf("Jets reconstructed in CTS+EMCAL \n");
1692   
1693   if(fPbPb) printf("PbPb events, pT cut in jet cone energy reconstruction %2.1f \n", fJetPtThreshold);
1694   else printf("pp events, pT cut in jet cone energy reconstruction %2.1f \n", fJetPtThresPbPb);
1695   
1696   printf("If pT of trigger < %2.3f, select jets as in pp? \n", fPtTriggerSelectionCut);
1697   
1698   printf("Phi gamma-Leading        <     %3.2f\n", fDeltaPhiMaxCut) ; 
1699   printf("Phi gamma-Leading        >     %3.2f\n", fDeltaPhiMinCut) ;
1700   printf("pT Leading / pT Trigger  <     %3.2f\n", fLeadingRatioMaxCut) ; 
1701   printf("pT Leading / pT Trigger  >     %3.2f\n", fLeadingRatioMinCut) ;
1702   
1703   if(fSelect == 2){
1704     printf("pT Jet / pT Gamma                     <    %3.2f\n", fJetRatioMaxCut) ; 
1705     printf("pT Jet / pT Gamma                     >    %3.2f\n", fJetRatioMinCut) ;
1706     printf("pT Jet (Only CTS)/ pT Trigger   <    %3.2f\n", fJetCTSRatioMaxCut) ; 
1707     printf("pT Jet (Only CTS)/ pT Trigger   >    %3.2f\n", fJetCTSRatioMinCut) ;
1708   }
1709   else if(fSelect == 0)
1710     printf("Accept all reconstructed jets \n") ;
1711   else   if(fSelect == 1)
1712     printf("Accept jets depending on trigger energy \n") ;
1713   else 
1714     printf("Wrong jet selection option:   %d \n", fSelect) ;
1715   
1716   printf("Isolated Trigger?  %d\n", fSelectIsolated) ;
1717
1718