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