Correct and clean the vertex retrieval in case of SE or ME analysis
[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 #include "TClonesArray.h"
34 #include "TClass.h"
35 //#include "Riostream.h"
36
37 //---- Analysis system ----
38 #include "AliAODTrack.h"
39 #include "AliVCluster.h"
40 #include "AliCaloTrackReader.h"
41 #include "AliNeutralMesonSelection.h"
42 #include "AliAnaParticleJetLeadingConeCorrelation.h"  
43 #include "AliCaloPID.h"
44 #include "AliAODPWG4ParticleCorrelation.h"
45 #include "AliFiducialCut.h"
46 #include "AliAODTrack.h"
47
48 ClassImp(AliAnaParticleJetLeadingConeCorrelation)
49
50
51 //____________________________________________________________________________
52   AliAnaParticleJetLeadingConeCorrelation::AliAnaParticleJetLeadingConeCorrelation() : 
53     AliAnaPartCorrBaseClass(), fJetsOnlyInCTS(kFALSE), fPbPb(kFALSE),     
54     fSeveralConeAndPtCuts(0),  fReMakeJet(0),
55     fDeltaPhiMaxCut(0.), fDeltaPhiMinCut(0.), 
56     fLeadingRatioMaxCut(0.),  fLeadingRatioMinCut(0.), 
57     fJetCTSRatioMaxCut(0.), fJetCTSRatioMinCut(0.), 
58     fJetRatioMaxCut(0.),  fJetRatioMinCut(0.), 
59     fJetNCone(0),fJetNPt(0), fJetCone(0), 
60     fJetPtThreshold(0),fJetPtThresPbPb(0),
61     fPtTriggerSelectionCut(0.0), fSelect(0),
62     fSelectIsolated(0), 
63     //Histograms
64     fOutCont(0x0),
65     //Leading
66     fhChargedLeadingPt(0),fhChargedLeadingPhi(0),fhChargedLeadingEta(0),
67     fhChargedLeadingDeltaPt(0),fhChargedLeadingDeltaPhi(0),fhChargedLeadingDeltaEta(0),
68     fhChargedLeadingRatioPt(0),
69     fhNeutralLeadingPt(0),fhNeutralLeadingPhi(0),fhNeutralLeadingEta(0),
70     fhNeutralLeadingDeltaPt(0),fhNeutralLeadingDeltaPhi(0),fhNeutralLeadingDeltaEta(0),
71     fhNeutralLeadingRatioPt(0),fhChargedLeadingXi(0), fhNeutralLeadingXi(0),
72     fhChargedLeadingDeltaPhiRatioPt30(0), fhNeutralLeadingDeltaPhiRatioPt30(0),
73     fhChargedLeadingDeltaPhiRatioPt50(0), fhNeutralLeadingDeltaPhiRatioPt50(0),
74     //Jet
75     fhJetPt(0),fhJetRatioPt(0),fhJetDeltaPhi(0), fhJetDeltaEta(0),
76     fhJetLeadingRatioPt(0),fhJetLeadingDeltaPhi(0),fhJetLeadingDeltaEta(0),
77     fhJetFFz(0),fhJetFFxi(0),fhJetFFpt(0),fhJetNTracksInCone(0),
78     fhBkgPt(0),fhBkgRatioPt(0),fhBkgDeltaPhi(0), fhBkgDeltaEta(0),
79     fhBkgLeadingRatioPt(0),fhBkgLeadingDeltaPhi(0),fhBkgLeadingDeltaEta(0),
80     fhBkgFFz(0),fhBkgFFxi(0),fhBkgFFpt(0),fhBkgNTracksInCone(0),
81     //Several cones and thres histograms
82     fhJetPts(),fhJetRatioPts(),fhJetDeltaPhis(), fhJetDeltaEtas(),
83     fhJetLeadingRatioPts(),fhJetLeadingDeltaPhis(),fhJetLeadingDeltaEtas(),
84     fhJetFFzs(),fhJetFFxis(),fhJetFFpts(),fhJetNTracksInCones(),
85     fhBkgPts(),fhBkgRatioPts(),fhBkgDeltaPhis(), fhBkgDeltaEtas(),
86     fhBkgLeadingRatioPts(),fhBkgLeadingDeltaPhis(),fhBkgLeadingDeltaEtas(),
87     fhBkgFFzs(),fhBkgFFxis(),fhBkgFFpts(),fhBkgNTracksInCones()
88 {
89   //Default Ctor
90   
91   //Initialize parameters
92
93   for(Int_t i = 0; i<6; i++){
94     fJetXMin1[i]     = 0.0 ;
95     fJetXMin2[i]     = 0.0 ;
96     fJetXMax1[i]     = 0.0 ;
97     fJetXMax2[i]     = 0.0 ;
98     fBkgMean[i]      = 0.0 ;
99     fBkgRMS[i]       = 0.0 ;
100     if( i < 2 ){
101       fJetE1[i]        = 0.0 ;
102       fJetE2[i]        = 0.0 ;
103       fJetSigma1[i]    = 0.0 ;
104       fJetSigma2[i]    = 0.0 ;
105     }
106   }
107   
108   //Several cones and thres histograms
109   for(Int_t i = 0; i<5; i++){
110     fJetCones[i]         = 0.0 ;
111     fJetNameCones[i]     = ""  ;
112     fJetPtThres[i]      = 0.0 ;
113     fJetNamePtThres[i]  = ""  ;
114     for(Int_t j = 0; j<5; j++){
115       fhJetPts[i][j]=0 ;
116       fhJetRatioPts[i][j]=0 ;
117       fhJetDeltaPhis[i][j]=0 ; 
118       fhJetDeltaEtas[i][j]=0 ;
119       fhJetLeadingRatioPts[i][j]=0 ;
120       fhJetLeadingDeltaPhis[i][j]=0 ;
121       fhJetLeadingDeltaEtas[i][j]=0 ;
122       fhJetFFzs[i][j]=0 ;
123       fhJetFFxis[i][j]=0 ;
124       fhJetFFpts[i][j]=0 ;
125       fhJetNTracksInCones[i][j]=0 ;
126       fhBkgPts[i][j]=0 ;
127       fhBkgRatioPts[i][j]=0 ;
128       fhBkgDeltaPhis[i][j]=0 ; 
129       fhBkgDeltaEtas[i][j]=0 ;
130       fhBkgLeadingRatioPts[i][j]=0 ;
131       fhBkgLeadingDeltaPhis[i][j]=0 ;
132       fhBkgLeadingDeltaEtas[i][j]=0 ;
133       fhBkgFFzs[i][j]=0 ;
134       fhBkgFFxis[i][j]=0 ;
135       fhBkgFFpts[i][j]=0 ;
136       fhBkgNTracksInCones[i][j]=0 ;
137     }
138   }
139
140   InitParameters();
141
142 }
143 /*
144 //____________________________________________________________________________
145 AliAnaParticleJetLeadingConeCorrelation::AliAnaParticleJetLeadingConeCorrelation(const AliAnaParticleJetLeadingConeCorrelation & jetlc) :   
146   AliAnaPartCorrBaseClass(jetlc), fJetsOnlyInCTS(jetlc.fJetsOnlyInCTS), fPbPb(jetlc.fPbPb), 
147   fSeveralConeAndPtCuts(jetlc.fSeveralConeAndPtCuts),  fReMakeJet(jetlc. fReMakeJet),
148   fDeltaPhiMaxCut(jetlc. fDeltaPhiMaxCut), fDeltaPhiMinCut(jetlc.fDeltaPhiMinCut), 
149   fLeadingRatioMaxCut(jetlc.fLeadingRatioMaxCut),  fLeadingRatioMinCut(jetlc.fLeadingRatioMinCut), 
150   fJetCTSRatioMaxCut(jetlc.fJetCTSRatioMaxCut),
151   fJetCTSRatioMinCut(jetlc.fJetCTSRatioMinCut), fJetRatioMaxCut(jetlc.fJetRatioMaxCut),
152   fJetRatioMinCut(jetlc.fJetRatioMinCut),  fJetNCone(jetlc.fJetNCone),
153   fJetNPt(jetlc.fJetNPt), fJetCone(jetlc.fJetCone),
154   fJetPtThreshold(jetlc.fJetPtThreshold),fJetPtThresPbPb(jetlc.fJetPtThresPbPb),
155   fPtTriggerSelectionCut(jetlc.fPtTriggerSelectionCut), fSelect(jetlc.fSelect), 
156   fSelectIsolated(jetlc.fSelectIsolated),    
157   //Histograms
158   fOutCont(jetlc. fOutCont),
159   //Leading
160   fhChargedLeadingPt(jetlc.fhChargedLeadingPt), fhChargedLeadingPhi(jetlc.fhChargedLeadingPhi),
161   fhChargedLeadingEta(jetlc.fhChargedLeadingEta), fhChargedLeadingDeltaPt(jetlc.fhChargedLeadingDeltaPt),
162   fhChargedLeadingDeltaPhi(jetlc.fhChargedLeadingDeltaPhi),fhChargedLeadingDeltaEta(jetlc.fhChargedLeadingDeltaEta),
163   fhChargedLeadingRatioPt(jetlc.fhChargedLeadingRatioPt),
164   fhNeutralLeadingPt(jetlc.fhNeutralLeadingPt),fhNeutralLeadingPhi(jetlc.fhNeutralLeadingPhi),
165   fhNeutralLeadingEta(jetlc.fhNeutralLeadingEta), fhNeutralLeadingDeltaPt(jetlc.fhNeutralLeadingDeltaPt),
166   fhNeutralLeadingDeltaPhi(jetlc.fhNeutralLeadingDeltaPhi),fhNeutralLeadingDeltaEta(jetlc.fhNeutralLeadingDeltaEta),
167   fhNeutralLeadingRatioPt(jetlc.fhNeutralLeadingRatioPt),
168   fhChargedLeadingXi(jetlc.fhChargedLeadingXi), fhNeutralLeadingXi(jetlc.fhNeutralLeadingXi),
169   fhChargedLeadingDeltaPhiRatioPt30(jetlc.fhChargedLeadingDeltaPhiRatioPt30), fhNeutralLeadingDeltaPhiRatioPt30(jetlc.fhNeutralLeadingDeltaPhiRatioPt30),
170   fhChargedLeadingDeltaPhiRatioPt50(jetlc.fhChargedLeadingDeltaPhiRatioPt50), fhNeutralLeadingDeltaPhiRatioPt50(jetlc.fhNeutralLeadingDeltaPhiRatioPt50),
171   //Jet
172   fhJetPt(jetlc.fhJetPt),fhJetRatioPt(jetlc.fhJetRatioPt),fhJetDeltaPhi(jetlc.fhJetDeltaPhi), 
173   fhJetDeltaEta(jetlc.fhJetDeltaEta), fhJetLeadingRatioPt(jetlc.fhJetLeadingRatioPt),
174   fhJetLeadingDeltaPhi(jetlc.fhJetLeadingDeltaPhi),fhJetLeadingDeltaEta(jetlc.fhJetLeadingDeltaEta),
175   fhJetFFz(jetlc.fhJetFFz),fhJetFFxi(jetlc.fhJetFFxi),fhJetFFpt(jetlc.fhJetFFpt),
176   fhJetNTracksInCone(jetlc.fhJetNTracksInCone),
177   fhBkgPt(jetlc.fhBkgPt),fhBkgRatioPt(jetlc.fhBkgRatioPt),fhBkgDeltaPhi(jetlc.fhBkgDeltaPhi), 
178   fhBkgDeltaEta(jetlc.fhBkgDeltaEta), fhBkgLeadingRatioPt(jetlc.fhBkgLeadingRatioPt),
179   fhBkgLeadingDeltaPhi(jetlc.fhBkgLeadingDeltaPhi),fhBkgLeadingDeltaEta(jetlc.fhBkgLeadingDeltaEta),
180   fhBkgFFz(jetlc.fhBkgFFz),fhBkgFFxi(jetlc.fhBkgFFxi),fhBkgFFpt(jetlc.fhBkgFFpt),
181   fhBkgNTracksInCone(jetlc.fhBkgNTracksInCone),
182   //Several cones and thres histograms
183   fhJetPts(),fhJetRatioPts(),fhJetDeltaPhis(), fhJetDeltaEtas(),
184   fhJetLeadingRatioPts(),fhJetLeadingDeltaPhis(),fhJetLeadingDeltaEtas(),
185   fhJetFFzs(),fhJetFFxis(),fhJetFFpts(),fhJetNTracksInCones(),
186   fhBkgPts(),fhBkgRatioPts(),fhBkgDeltaPhis(), fhBkgDeltaEtas(),
187   fhBkgLeadingRatioPts(),fhBkgLeadingDeltaPhis(),fhBkgLeadingDeltaEtas(),
188   fhBkgFFzs(),fhBkgFFxis(),fhBkgFFpts(),fhBkgNTracksInCones()
189 {
190   // cpy ctor
191
192   for(Int_t i = 0; i<6; i++){
193     fJetXMin1[i]       = jetlc.fJetXMin1[i] ;
194     fJetXMin2[i]       = jetlc.fJetXMin2[i] ;
195     fJetXMax1[i]       = jetlc.fJetXMax1[i] ;
196     fJetXMax2[i]       = jetlc.fJetXMax2[i] ;
197     fBkgMean[i]        = jetlc.fBkgMean[i] ;
198     fBkgRMS[i]         = jetlc.fBkgRMS[i] ;
199     if( i < 2 ){
200       fJetE1[i]        = jetlc.fJetE1[i] ;
201       fJetE2[i]        = jetlc.fJetE2[i] ;
202       fJetSigma1[i]    = jetlc.fJetSigma1[i] ;
203       fJetSigma2[i]    = jetlc.fJetSigma2[i] ;
204     }
205   }          
206   
207   //Several cones and thres histograms
208   for(Int_t i = 0; i<5; i++){
209     fJetCones[i]        = jetlc.fJetCones[i] ;
210     fJetNameCones[i]    = jetlc.fJetNameCones[i] ;
211     fJetPtThres[i]      = jetlc.fJetPtThres[i] ;
212     fJetNamePtThres[i]  = jetlc.fJetNamePtThres[i] ;
213     for(Int_t j = 0; j<5; j++){
214       fhJetPts[i][j] = jetlc.fhJetPts[i][j] ;
215       fhJetRatioPts[i][j] = jetlc.fhJetRatioPts[i][j] ;
216       fhJetDeltaPhis[i][j] = jetlc.fhJetDeltaPhis[i][j] ; 
217       fhJetDeltaEtas[i][j] = jetlc.fhJetDeltaEtas[i][j] ;
218       fhJetLeadingRatioPts[i][j] = jetlc.fhJetLeadingRatioPts[i][j] ;
219       fhJetLeadingDeltaPhis[i][j] = jetlc.fhJetLeadingDeltaPhis[i][j] ;
220       fhJetLeadingDeltaEtas[i][j] = jetlc.fhJetLeadingDeltaEtas[i][j] ;
221       fhJetFFzs[i][j] = jetlc.fhJetFFzs[i][j] ;
222       fhJetFFxis[i][j] = jetlc.fhJetFFxis[i][j] ;
223       fhJetFFpts[i][j] = jetlc.fhJetFFpts[i][j] ;
224       fhJetNTracksInCones[i][j] = fhJetNTracksInCones[i][j] ;
225       fhBkgPts[i][j] = jetlc.fhBkgPts[i][j] ;
226       fhBkgRatioPts[i][j] = jetlc.fhBkgRatioPts[i][j] ;
227       fhBkgDeltaPhis[i][j] = jetlc.fhBkgDeltaPhis[i][j] ; 
228       fhBkgDeltaEtas[i][j] = jetlc.fhBkgDeltaEtas[i][j] ;
229       fhBkgLeadingRatioPts[i][j] = jetlc.fhBkgLeadingRatioPts[i][j] ;
230       fhBkgLeadingDeltaPhis[i][j] = jetlc.fhBkgLeadingDeltaPhis[i][j] ;
231       fhBkgLeadingDeltaEtas[i][j] = jetlc.fhBkgLeadingDeltaEtas[i][j] ;
232       fhBkgFFzs[i][j] = jetlc.fhBkgFFzs[i][j] ;
233       fhBkgFFxis[i][j] = jetlc.fhBkgFFxis[i][j] ;
234       fhBkgFFpts[i][j] = jetlc.fhBkgFFpts[i][j] ;
235       fhBkgNTracksInCones[i][j] = jetlc.fhBkgNTracksInCones[i][j] ;
236     }
237   }  
238 }
239
240 //_________________________________________________________________________
241 AliAnaParticleJetLeadingConeCorrelation & AliAnaParticleJetLeadingConeCorrelation::operator = (const AliAnaParticleJetLeadingConeCorrelation & jetlc)
242 {
243   // assignment operator
244   
245   if(this == &jetlc)return *this;
246   ((AliAnaPartCorrBaseClass *)this)->operator=(jetlc);
247   
248   fSeveralConeAndPtCuts  = jetlc.fSeveralConeAndPtCuts ; 
249   fPbPb                  = jetlc.fPbPb ;
250   fReMakeJet             = jetlc.fReMakeJet ;
251   fJetsOnlyInCTS         = jetlc.fJetsOnlyInCTS;
252
253   fDeltaPhiMaxCut        = jetlc.fDeltaPhiMaxCut ; 
254   fDeltaPhiMinCut        = jetlc.fDeltaPhiMinCut ; 
255   fLeadingRatioMaxCut    = jetlc.fLeadingRatioMaxCut ;
256   fLeadingRatioMinCut    = jetlc.fLeadingRatioMinCut ;
257   
258   fJetCTSRatioMaxCut     = jetlc.fJetCTSRatioMaxCut ;
259   fJetCTSRatioMinCut     = jetlc.fJetCTSRatioMinCut ; 
260   fJetRatioMaxCut        = jetlc.fJetRatioMaxCut ;
261   fJetRatioMinCut        = jetlc.fJetRatioMinCut ; 
262  
263   fJetNCone              = jetlc.fJetNCone ;
264   fJetNPt                = jetlc.fJetNPt ; fJetCone = jetlc.fJetCone ; 
265   fJetPtThreshold        = jetlc.fJetPtThreshold ;
266   fJetPtThresPbPb        = jetlc.fJetPtThresPbPb ;
267   fPtTriggerSelectionCut = jetlc.fPtTriggerSelectionCut ;
268   fSelect                = jetlc.fSelect ;  
269   fSelectIsolated        = jetlc.fSelectIsolated ;
270
271   for(Int_t i = 0; i<6; i++){
272     fJetXMin1[i]       = jetlc.fJetXMin1[i] ;
273     fJetXMin2[i]       = jetlc.fJetXMin2[i] ;
274     fJetXMax1[i]       = jetlc.fJetXMax1[i] ;
275     fJetXMax2[i]       = jetlc.fJetXMax2[i] ;
276     fBkgMean[i]        = jetlc.fBkgMean[i] ;
277     fBkgRMS[i]         = jetlc.fBkgRMS[i] ;
278     if( i < 2 ){
279       fJetE1[i]        = jetlc.fJetE1[i] ;
280       fJetE2[i]        = jetlc.fJetE2[i] ;
281       fJetSigma1[i]    = jetlc.fJetSigma1[i] ;
282       fJetSigma2[i]    = jetlc.fJetSigma2[i] ;
283     }
284   }     
285
286   //Histograms
287   fOutCont = jetlc. fOutCont ;
288   fhChargedLeadingPt = jetlc.fhChargedLeadingPt; fhChargedLeadingPhi = jetlc.fhChargedLeadingPhi;
289   fhChargedLeadingEta = jetlc.fhChargedLeadingEta; fhChargedLeadingDeltaPt = jetlc.fhChargedLeadingDeltaPt;
290   fhChargedLeadingDeltaPhi = jetlc.fhChargedLeadingDeltaPhi;fhChargedLeadingDeltaEta = jetlc.fhChargedLeadingDeltaEta;
291   fhChargedLeadingRatioPt = jetlc.fhChargedLeadingRatioPt;
292   fhNeutralLeadingPt = jetlc.fhNeutralLeadingPt;fhNeutralLeadingPhi = jetlc.fhNeutralLeadingPhi;
293   fhNeutralLeadingEta = jetlc.fhNeutralLeadingEta; fhNeutralLeadingDeltaPt = jetlc.fhNeutralLeadingDeltaPt;
294   fhNeutralLeadingDeltaPhi = jetlc.fhNeutralLeadingDeltaPhi;fhNeutralLeadingDeltaEta = jetlc.fhNeutralLeadingDeltaEta;
295   fhNeutralLeadingRatioPt = jetlc.fhNeutralLeadingRatioPt;
296   fhChargedLeadingXi = jetlc.fhChargedLeadingXi; 
297   fhNeutralLeadingXi = jetlc.fhNeutralLeadingXi;
298
299   fhChargedLeadingDeltaPhiRatioPt30 = jetlc.fhChargedLeadingDeltaPhiRatioPt30; 
300   fhNeutralLeadingDeltaPhiRatioPt30 = jetlc.fhNeutralLeadingDeltaPhiRatioPt30;
301   fhChargedLeadingDeltaPhiRatioPt50 = jetlc.fhChargedLeadingDeltaPhiRatioPt50; 
302   fhNeutralLeadingDeltaPhiRatioPt50 = jetlc.fhNeutralLeadingDeltaPhiRatioPt50;  
303         
304   fhJetPt = jetlc.fhJetPt;fhJetRatioPt = jetlc.fhJetRatioPt;fhJetDeltaPhi = jetlc.fhJetDeltaPhi; 
305   fhJetDeltaEta = jetlc.fhJetDeltaEta; fhJetLeadingRatioPt = jetlc.fhJetLeadingRatioPt;
306   fhJetLeadingDeltaPhi = jetlc.fhJetLeadingDeltaPhi;fhJetLeadingDeltaEta = jetlc.fhJetLeadingDeltaEta;
307   fhJetFFz = jetlc.fhJetFFz;fhJetFFxi = jetlc.fhJetFFxi;fhJetFFpt = jetlc.fhJetFFpt;
308   fhJetNTracksInCone = jetlc.fhJetNTracksInCone;
309   fhBkgPt = jetlc.fhBkgPt;fhBkgRatioPt = jetlc.fhBkgRatioPt;fhBkgDeltaPhi = jetlc.fhBkgDeltaPhi; 
310   fhBkgDeltaEta = jetlc.fhBkgDeltaEta; fhBkgLeadingRatioPt = jetlc.fhBkgLeadingRatioPt;
311   fhBkgLeadingDeltaPhi = jetlc.fhBkgLeadingDeltaPhi;fhBkgLeadingDeltaEta = jetlc.fhBkgLeadingDeltaEta;
312   fhBkgFFz = jetlc.fhBkgFFz;fhBkgFFxi = jetlc.fhBkgFFxi;fhBkgFFpt = jetlc.fhBkgFFpt;
313   fhBkgNTracksInCone = jetlc.fhBkgNTracksInCone;
314
315
316   //Several cones and thres histograms
317   for(Int_t i = 0; i<5; i++){
318     fJetCones[i]        = jetlc.fJetCones[i] ;
319     fJetNameCones[i]    = jetlc.fJetNameCones[i] ;
320     fJetPtThres[i]      = jetlc.fJetPtThres[i] ;
321     fJetNamePtThres[i]  = jetlc.fJetNamePtThres[i] ;
322     
323     for(Int_t j = 0; j<5; j++){
324       fhJetPts[i][j] = jetlc.fhJetPts[i][j] ;
325       fhJetRatioPts[i][j] = jetlc.fhJetRatioPts[i][j] ;
326       fhJetDeltaPhis[i][j] = jetlc.fhJetDeltaPhis[i][j] ; 
327       fhJetDeltaEtas[i][j] = jetlc.fhJetDeltaEtas[i][j] ;
328       fhJetLeadingRatioPts[i][j] = jetlc.fhJetLeadingRatioPts[i][j] ;
329       fhJetLeadingDeltaPhis[i][j] = jetlc.fhJetLeadingDeltaPhis[i][j] ;
330       fhJetLeadingDeltaEtas[i][j] = jetlc.fhJetLeadingDeltaEtas[i][j] ;
331       fhJetFFzs[i][j] = jetlc.fhJetFFzs[i][j] ;
332       fhJetFFxis[i][j] = jetlc.fhJetFFxis[i][j] ;
333       fhJetFFpts[i][j] = jetlc.fhJetFFpts[i][j] ;
334       fhJetNTracksInCones[i][j] = fhJetNTracksInCones[i][j] ;
335       fhBkgPts[i][j] = jetlc.fhBkgPts[i][j] ;
336       fhBkgRatioPts[i][j] = jetlc.fhBkgRatioPts[i][j] ;
337       fhBkgDeltaPhis[i][j] = jetlc.fhBkgDeltaPhis[i][j] ; 
338       fhBkgDeltaEtas[i][j] = jetlc.fhBkgDeltaEtas[i][j] ;
339       fhBkgLeadingRatioPts[i][j] = jetlc.fhBkgLeadingRatioPts[i][j] ;
340       fhBkgLeadingDeltaPhis[i][j] = jetlc.fhBkgLeadingDeltaPhis[i][j] ;
341       fhBkgLeadingDeltaEtas[i][j] = jetlc.fhBkgLeadingDeltaEtas[i][j] ;
342       fhBkgFFzs[i][j] = jetlc.fhBkgFFzs[i][j] ;
343       fhBkgFFxis[i][j] = jetlc.fhBkgFFxis[i][j] ;
344       fhBkgFFpts[i][j] = jetlc.fhBkgFFpts[i][j] ;
345       fhBkgNTracksInCones[i][j] = jetlc.fhBkgNTracksInCones[i][j] ;
346     }
347   }      
348
349   return *this;
350
351 }
352 */
353 //____________________________________________________________________________
354 AliAnaParticleJetLeadingConeCorrelation::~AliAnaParticleJetLeadingConeCorrelation() 
355 {
356    // Remove all pointers except analysis output pointers.
357 }
358
359 //____________________________________________________________________________
360 Double_t AliAnaParticleJetLeadingConeCorrelation::CalculateJetRatioLimit(const Double_t ptg, const Double_t *par, const Double_t *x) const {
361   //Calculate the ratio of the jet and trigger particle limit for the selection
362   //WARNING: need to check what it does
363   //printf("CalculateLimit: x1 %2.3f, x2%2.3f\n",x[0],x[1]);
364   Double_t ePP = par[0] + par[1] * ptg ;
365   Double_t sPP = par[2] + par[3] * ptg ;
366   Double_t f   = x[0]   + x[1]   * ptg ;
367   Double_t ePbPb = ePP + par[4] ;
368   Double_t sPbPb = TMath::Sqrt(sPP*sPP+ par[5]*par[5]) ;
369   Double_t rat = (ePbPb - sPbPb * f) / ptg ;
370   //printf("CalculateLimit: ePP %2.3f, sPP %2.3f, f %2.3f\n", ePP, sPP, f);
371   //printf("CalculateLimit: ePbPb %2.3f, sPbPb %2.3f, rat %2.3f\n", ePbPb, sPbPb, rat);
372   return rat ;
373 }
374
375 //____________________________________________________________________________
376 void AliAnaParticleJetLeadingConeCorrelation::FillJetHistos(AliAODPWG4ParticleCorrelation * particle, const TLorentzVector  leading, const TLorentzVector jet, const TString type, const TString lastname)
377 {
378   //Fill jet and background histograms 
379   Double_t ptTrig = particle->Pt();
380   Double_t ptJet = jet.Pt();
381   Double_t ptLead = leading.Pt();
382   Double_t phiTrig = particle->Phi();
383   Double_t phiJet = jet.Phi();
384   if(phiJet < 0) phiJet+=TMath::TwoPi();
385   Double_t phiLead = leading.Phi();
386   if(phiLead < 0) phiLead+=TMath::TwoPi();
387   Double_t etaTrig = particle->Eta();
388   Double_t etaJet = jet.Eta();
389   Double_t etaLead = leading.Eta();
390   
391   TH2F *h1 = 0x0;
392   h1 = dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sPt%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
393   if(h1)h1->Fill(ptTrig,ptJet);
394   
395   TH2F *h2 = 0x0;
396   h2 = dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sRatioPt%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
397   if(h2) h2->Fill(ptTrig,ptJet/ptTrig);
398   
399   TH2F *h3 = 0x0;
400   h3 = dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sLeadingRatioPt%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
401   if(h3)h3->Fill(ptTrig,ptLead/ptJet);
402   
403   //   dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sPhi%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())))->
404   //     Fill(ptTrig,phiJet);
405   TH2F *h4 = 0x0;
406   h4 = dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sDeltaPhi%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
407   if(h4) h4->Fill(ptTrig,phiJet-phiTrig);
408   TH2F *h5 = 0x0;
409   h5 = dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sLeadingDeltaPhi%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
410   if(h5) h5->Fill(ptTrig,phiJet-phiLead);
411   
412   //   dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sEta%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())))->
413   //     Fill(ptTrig,etaJet);
414   TH2F *h6 = 0x0;
415   h6 = dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sDeltaEta%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
416   if(h6) h6->Fill(ptTrig,etaJet-etaTrig);
417   TH2F *h7 = 0x0;
418   h7 = dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sLeadingDeltaEta%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
419   if(h7) h7->Fill(ptTrig,etaJet-etaLead);
420   
421   //Construct fragmentation function
422   TObjArray * pl = new TObjArray;
423   
424   if(type == "Jet") pl = particle->GetObjArray(Form("%sTracks",GetAODObjArrayName().Data()));
425   else if(type == "Bkg") particle->GetObjArray(Form("%sTracksBkg",GetAODObjArrayName().Data()));
426   
427   if(!pl) return ;
428   
429   //Different pt cut for jet particles in different collisions systems
430   //Only needed when jet is recalculated from AODs
431   Float_t ptcut = fJetPtThreshold;
432   if(fPbPb && !fSeveralConeAndPtCuts && ptTrig > fPtTriggerSelectionCut)  ptcut = fJetPtThresPbPb ;
433   
434   TVector3 p3;
435   Int_t nTracksInCone = 0; 
436   
437   for(Int_t ipr = 0;ipr < pl->GetEntriesFast() ; ipr ++ ){
438     AliAODTrack* track = dynamic_cast<AliAODTrack *>(pl->At(ipr)) ;
439     p3.SetXYZ(track->Px(),track->Py(),track->Pz());
440     
441     //Recheck if particle is in jet cone
442     if(fReMakeJet || fSeveralConeAndPtCuts)
443       if(!IsParticleInJetCone(p3.Eta(), p3.Phi(), leading.Eta(), leading.Phi()) ) continue ; 
444     
445     nTracksInCone++; 
446     
447     TH2F *ha = 0x0;
448     ha =dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sFFz%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
449     if(ha) ha->Fill(ptTrig,p3.Pt()/ptTrig);
450     TH2F *hb = 0x0;
451     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/p3.Pt()));
453     TH2F *hc = 0x0;
454     hc =dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sFFpt%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
455     if(hc) hc->Fill(ptTrig,p3.Pt());
456     
457   }//track loop
458   
459   if(nTracksInCone > 0) dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sNTracksInCone%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())))
460     ->Fill(ptTrig, nTracksInCone);
461   
462 }
463
464 //________________________________________________________________________
465 TList *  AliAnaParticleJetLeadingConeCorrelation::GetCreateOutputObjects()
466 {  
467   // Create histograms to be saved in output file and 
468   // store them in fOutCont
469   
470   if(GetDebug()>1) printf("AliAnaParticleJetLeadingConeCorrelation::GetCreateOutputObjects() - Init histograms \n");
471   
472   fOutCont = new TList() ; 
473   fOutCont->SetName("ParticleJetLeadingInConeCorrelationHistograms") ; 
474   
475   Int_t nptbins  = GetHistoPtBins();
476   Int_t nphibins = GetHistoPhiBins();
477   Int_t netabins = GetHistoEtaBins();
478   Float_t ptmax  = GetHistoPtMax();
479   Float_t phimax = GetHistoPhiMax();
480   Float_t etamax = GetHistoEtaMax();
481   Float_t ptmin  = GetHistoPtMin();
482   Float_t phimin = GetHistoPhiMin();
483   Float_t etamin = GetHistoEtaMin();    
484   
485   fhChargedLeadingPt  = new TH2F("ChargedLeadingPt","p_{T leading charge} vs p_{T trigger}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
486   fhChargedLeadingPt->SetYTitle("p_{T leading charge}");
487   fhChargedLeadingPt->SetXTitle("p_{T trigger} (GeV/c)");
488   
489   fhChargedLeadingPhi  = new TH2F("ChargedLeadingPhi","#phi_{h^{#pm}}  vs p_{T trigger}", nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
490   fhChargedLeadingPhi->SetYTitle("#phi_{h^{#pm}} (rad)");
491   fhChargedLeadingPhi->SetXTitle("p_{T trigger} (GeV/c)");
492   
493   fhChargedLeadingEta  = new TH2F("ChargedLeadingEta","#eta_{h^{#pm}}  vs p_{T trigger}",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
494   fhChargedLeadingEta->SetYTitle("#eta_{h^{#pm}} ");
495   fhChargedLeadingEta->SetXTitle("p_{T trigger} (GeV/c)");
496   
497   fhChargedLeadingDeltaPt  = new TH2F("ChargedLeadingDeltaPt","p_{T trigger} - p_{T h^{#pm}} vs p_{T trigger}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
498   fhChargedLeadingDeltaPt->SetYTitle("#Delta p_{T} (GeV/c)");
499   fhChargedLeadingDeltaPt->SetXTitle("p_{T trigger} (GeV/c)");
500   
501   fhChargedLeadingDeltaPhi  = new TH2F("ChargedLeadingDeltaPhi","#phi_{trigger} - #phi_{h^{#pm}} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,TMath::TwoPi()); 
502   fhChargedLeadingDeltaPhi->SetYTitle("#Delta #phi (rad)");
503   fhChargedLeadingDeltaPhi->SetXTitle("p_{T trigger} (GeV/c)");
504   
505   fhChargedLeadingDeltaEta  = new TH2F("ChargedLeadingDeltaEta","#eta_{trigger} - #eta_{h^{#pm}} vs p_{T trigger}",nptbins,ptmin,ptmax,120,-2,2); 
506   fhChargedLeadingDeltaEta->SetYTitle("#Delta #eta");
507   fhChargedLeadingDeltaEta->SetXTitle("p_{T trigger} (GeV/c)");
508   
509   fhChargedLeadingRatioPt  = new TH2F("ChargedLeadingRatioPt","p_{T leading charge} /p_{T trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,2); 
510   fhChargedLeadingRatioPt->SetYTitle("p_{T lead charge} /p_{T trigger}");
511   fhChargedLeadingRatioPt->SetXTitle("p_{T trigger} (GeV/c)");
512   
513   fhChargedLeadingXi  = new TH2F("ChargedLeadingXi","ln(p_{T trigger} / p_{T leading charge} ) vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,10); 
514   fhChargedLeadingXi->SetYTitle("#xi");
515   fhChargedLeadingXi->SetXTitle("p_{T trigger} (GeV/c)");
516         
517   fOutCont->Add(fhChargedLeadingPt) ;
518   fOutCont->Add(fhChargedLeadingPhi) ;
519   fOutCont->Add(fhChargedLeadingEta) ;
520   fOutCont->Add(fhChargedLeadingDeltaPt) ; 
521   fOutCont->Add(fhChargedLeadingDeltaPhi) ; 
522   fOutCont->Add(fhChargedLeadingDeltaEta) ; 
523   fOutCont->Add(fhChargedLeadingRatioPt) ;
524   fOutCont->Add(fhChargedLeadingXi) ;
525         
526   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); 
527   fhChargedLeadingDeltaPhiRatioPt30->SetXTitle("#Delta #phi (rad)");
528   fhChargedLeadingDeltaPhiRatioPt30->SetYTitle("p_{T leading} / p_{T trigger}");        
529         
530   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); 
531   fhChargedLeadingDeltaPhiRatioPt50->SetXTitle("#Delta #phi (rad)");
532   fhChargedLeadingDeltaPhiRatioPt50->SetYTitle("p_{T leading} / p_{T trigger}");          
533         
534   fOutCont->Add(fhChargedLeadingDeltaPhiRatioPt30) ; 
535   fOutCont->Add(fhChargedLeadingDeltaPhiRatioPt50) ; 
536         
537   if(!fJetsOnlyInCTS){
538     
539     fhNeutralLeadingPt  = new TH2F("NeutralLeadingPt","p_{T leading #pi^{0}} vs p_{T trigger}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
540     fhNeutralLeadingPt->SetYTitle("p_{T leading #pi^{0}}");
541     fhNeutralLeadingPt->SetXTitle("p_{T trigger} (GeV/c)");
542     
543     fhNeutralLeadingPhi  = new TH2F("NeutralLeadingPhi","#phi_{#pi^{0}}  vs p_{T trigger}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
544     fhNeutralLeadingPhi->SetYTitle("#phi_{#pi^{0}} (rad)");
545     fhNeutralLeadingPhi->SetXTitle("p_{T trigger} (GeV/c)");
546     
547     fhNeutralLeadingEta  = new TH2F("NeutralLeadingEta","#eta_{#pi^{0}}  vs p_{T trigger}",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
548     fhNeutralLeadingEta->SetYTitle("#eta_{#pi^{0}} ");
549     fhNeutralLeadingEta->SetXTitle("p_{T trigger} (GeV/c)");
550     
551     fhNeutralLeadingDeltaPt  = new TH2F("NeutralLeadingDeltaPt","p_{T trigger} - p_{T #pi^{0}} vs p_{T trigger}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
552     fhNeutralLeadingDeltaPt->SetYTitle("#Delta p_{T} (GeV/c)");
553     fhNeutralLeadingDeltaPt->SetXTitle("p_{T trigger} (GeV/c)");
554     
555     fhNeutralLeadingDeltaPhi  = new TH2F("NeutralLeadingDeltaPhi","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,TMath::TwoPi()); 
556     fhNeutralLeadingDeltaPhi->SetYTitle("#Delta #phi (rad)");
557     fhNeutralLeadingDeltaPhi->SetXTitle("p_{T trigger} (GeV/c)");
558     
559     fhNeutralLeadingDeltaEta  = new TH2F("NeutralLeadingDeltaEta","#eta_{trigger} - #eta_{#pi^{0}} vs p_{T trigger}",nptbins,ptmin,ptmax,120,-2,2); 
560     fhNeutralLeadingDeltaEta->SetYTitle("#Delta #eta");
561     fhNeutralLeadingDeltaEta->SetXTitle("p_{T trigger} (GeV/c)");
562     
563     fhNeutralLeadingRatioPt  = new TH2F("NeutralLeadingRatioPt","p_{T leading #pi^{0}} /p_{T trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,2); 
564     fhNeutralLeadingRatioPt->SetYTitle("p_{T lead #pi^{0}} /p_{T trigger}");
565     fhNeutralLeadingRatioPt->SetXTitle("p_{T trigger} (GeV/c)");
566     
567     fhNeutralLeadingXi  = new TH2F("NeutralLeadingXi","ln(p_{T trigger} / p_{T leading #pi^{0}} ) vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,10); 
568     fhNeutralLeadingXi->SetYTitle("#xi");
569     fhNeutralLeadingXi->SetXTitle("p_{T trigger} (GeV/c)");
570     
571     fOutCont->Add(fhNeutralLeadingPt) ;
572     fOutCont->Add(fhNeutralLeadingPhi) ;
573     fOutCont->Add(fhNeutralLeadingEta) ;
574     fOutCont->Add(fhNeutralLeadingDeltaPt) ; 
575     fOutCont->Add(fhNeutralLeadingDeltaPhi) ; 
576     fOutCont->Add(fhNeutralLeadingDeltaEta) ; 
577     fOutCont->Add(fhNeutralLeadingRatioPt) ;
578     fOutCont->Add(fhNeutralLeadingXi) ;
579           
580         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); 
581         fhNeutralLeadingDeltaPhiRatioPt30->SetXTitle("#Delta #phi (rad)");
582         fhNeutralLeadingDeltaPhiRatioPt30->SetYTitle("p_{T leading} / p_{T trigger}");
583           
584         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); 
585         fhNeutralLeadingDeltaPhiRatioPt50->SetXTitle("#Delta #phi (rad)");
586         fhNeutralLeadingDeltaPhiRatioPt50->SetYTitle("p_{T leading} / p_{T trigger}");
587         fOutCont->Add(fhNeutralLeadingDeltaPhiRatioPt30) ; 
588         fOutCont->Add(fhNeutralLeadingDeltaPhiRatioPt50) ; 
589
590   }
591   
592   if(!fSeveralConeAndPtCuts){// not several cones
593     
594     //Jet Distributions
595     fhJetPt  = new TH2F("JetPt","p_{T  jet} vs p_{T trigger}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
596     fhJetPt->SetYTitle("p_{T  jet}");
597     fhJetPt->SetXTitle("p_{T trigger} (GeV/c)");
598     
599     fhJetRatioPt  = new TH2F("JetRatioPt","p_{T  jet}/p_{T trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,2); 
600     fhJetRatioPt->SetYTitle("p_{T  jet}/p_{T trigger}");
601     fhJetRatioPt->SetXTitle("p_{T trigger} (GeV/c)");
602     
603     fhJetDeltaPhi  = new TH2F("JetDeltaPhi","#phi_{jet} - #phi_{trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,TMath::TwoPi()); 
604     fhJetDeltaPhi->SetYTitle("#Delta #phi (rad)");
605     fhJetDeltaPhi->SetXTitle("p_{T trigger} (GeV/c)");
606     
607     fhJetDeltaEta  = new TH2F("JetDeltaEta","#eta_{jet} - #eta_{trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,120,-2,2); 
608     fhJetDeltaEta->SetYTitle("#Delta #eta");
609     fhJetDeltaEta->SetXTitle("p_{T trigger} (GeV/c)");
610     
611     fhJetLeadingRatioPt  = new TH2F("JetLeadingRatioPt","p_{T  jet} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,2); 
612     fhJetLeadingRatioPt->SetYTitle("p_{T  leading}/p_{T jet}");
613     fhJetLeadingRatioPt->SetXTitle("p_{T trigger} (GeV/c)");
614     
615     fhJetLeadingDeltaPhi  = new TH2F("JetLeadingDeltaPhi","#phi_{jet} - #phi_{leading} vs p_{T trigger}",nptbins,ptmin,ptmax,120,-TMath::Pi(),TMath::Pi()); 
616     fhJetLeadingDeltaPhi->SetYTitle("#Delta #phi (rad)");
617     fhJetLeadingDeltaPhi->SetXTitle("p_{T trigger} (GeV/c)");
618     
619     fhJetLeadingDeltaEta  = new TH2F("JetLeadingDeltaEta","#eta_{jet} - #eta_{leading} vs p_{T trigger}",nptbins,ptmin,ptmax,120,-2,2); 
620     fhJetLeadingDeltaEta->SetYTitle("#Delta #eta");
621     fhJetLeadingDeltaEta->SetXTitle("p_{T trigger} (GeV/c)");
622     
623     fhJetFFz  = new TH2F("JetFFz","z = p_{T i charged}/p_{T trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,200,0.,2); 
624     fhJetFFz->SetYTitle("z");
625     fhJetFFz->SetXTitle("p_{T trigger}");
626     
627     fhJetFFxi  = new TH2F("JetFFxi","#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger}",nptbins,ptmin,ptmax,100,0.,10.); 
628     fhJetFFxi->SetYTitle("#xi");
629     fhJetFFxi->SetXTitle("p_{T trigger}");
630     
631     fhJetFFpt  = new TH2F("JetFFpt","#xi = p_{T i charged}) vs p_{T trigger}",nptbins,ptmin,ptmax,200,0.,50.); 
632     fhJetFFpt->SetYTitle("p_{T charged hadron}");
633     fhJetFFpt->SetXTitle("p_{T trigger}");
634     
635     fhJetNTracksInCone  = new TH2F("JetNTracksInCone","N particles in cone vs p_{T trigger}",nptbins,ptmin,ptmax,5000,0, 5000); 
636     fhJetNTracksInCone->SetYTitle("N tracks in jet cone");
637     fhJetNTracksInCone->SetXTitle("p_{T trigger} (GeV/c)");
638     
639     fOutCont->Add(fhJetPt) ; 
640     fOutCont->Add(fhJetRatioPt) ; 
641     fOutCont->Add(fhJetDeltaPhi) ;
642     fOutCont->Add(fhJetDeltaEta) ;
643     fOutCont->Add(fhJetLeadingRatioPt) ;
644     fOutCont->Add(fhJetLeadingDeltaPhi) ;
645     fOutCont->Add(fhJetLeadingDeltaEta) ;
646     fOutCont->Add(fhJetFFz) ;
647     fOutCont->Add(fhJetFFxi) ;
648     fOutCont->Add(fhJetFFpt) ;
649     fOutCont->Add(fhJetNTracksInCone) ;
650     
651     //Bkg Distributions
652     fhBkgPt  = new TH2F("BkgPt","p_{T  bkg} vs p_{T trigger}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
653     fhBkgPt->SetYTitle("p_{T  bkg}");
654     fhBkgPt->SetXTitle("p_{T trigger} (GeV/c)");
655     
656     fhBkgRatioPt  = new TH2F("BkgRatioPt","p_{T  bkg}/p_{T trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,2); 
657     fhBkgRatioPt->SetYTitle("p_{T  bkg}/p_{T trigger}");
658     fhBkgRatioPt->SetXTitle("p_{T trigger} (GeV/c)");
659     
660     fhBkgDeltaPhi  = new TH2F("BkgDeltaPhi","#phi_{bkg} - #phi_{trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,TMath::TwoPi()); 
661     fhBkgDeltaPhi->SetYTitle("#Delta #phi (rad)");
662     fhBkgDeltaPhi->SetXTitle("p_{T trigger} (GeV/c)");
663     
664     fhBkgDeltaEta  = new TH2F("BkgDeltaEta","#eta_{bkg} - #eta_{trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,120,-2,2); 
665     fhBkgDeltaEta->SetYTitle("#Delta #eta");
666     fhBkgDeltaEta->SetXTitle("p_{T trigger} (GeV/c)");
667     
668     fhBkgLeadingRatioPt  = new TH2F("BkgLeadingRatioPt","p_{T  bkg} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,2); 
669     fhBkgLeadingRatioPt->SetYTitle("p_{T  leading}/p_{T bkg}");
670     fhBkgLeadingRatioPt->SetXTitle("p_{T trigger} (GeV/c)");
671     
672     fhBkgLeadingDeltaPhi  = new TH2F("BkgLeadingDeltaPhi","#phi_{bkg} - #phi_{leading} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,TMath::TwoPi()); 
673     fhBkgLeadingDeltaPhi->SetYTitle("#Delta #phi (rad)");
674     fhBkgLeadingDeltaPhi->SetXTitle("p_{T trigger} (GeV/c)");
675     
676     fhBkgLeadingDeltaEta  = new TH2F("BkgLeadingDeltaEta","#eta_{bkg} - #eta_{leading} vs p_{T trigger}",nptbins,ptmin,ptmax,120,-2,2); 
677     fhBkgLeadingDeltaEta->SetYTitle("#Delta #eta");
678     fhBkgLeadingDeltaEta->SetXTitle("p_{T trigger} (GeV/c)");
679     
680     fhBkgFFz  = new TH2F("BkgFFz","z = p_{T i charged}/p_{T trigger} vs p_{T trigger}", nptbins,ptmin,ptmax,200,0.,2); 
681     fhBkgFFz->SetYTitle("z");
682     fhBkgFFz->SetXTitle("p_{T trigger}");
683     
684     fhBkgFFxi  = new TH2F("BkgFFxi","#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger}", nptbins,ptmin,ptmax,100,0.,10.); 
685     fhBkgFFxi->SetYTitle("#xi");
686     fhBkgFFxi->SetXTitle("p_{T trigger}");
687     
688     fhBkgFFpt  = new TH2F("BkgFFpt","p_{T charged hadron } vs p_{T trigger}", nptbins,ptmin,ptmax,200,0.,50.); 
689     fhBkgFFpt->SetYTitle("p_{T charged} hadron");
690     fhBkgFFpt->SetXTitle("p_{T trigger}");
691     
692     fhBkgNTracksInCone  = new TH2F("BkgNTracksInCone","N particles in cone vs p_{T trigger}",nptbins,ptmin,ptmax,5000,0, 5000); 
693     fhBkgNTracksInCone->SetYTitle("N tracks in bkg cone");
694     fhBkgNTracksInCone->SetXTitle("p_{T trigger} (GeV/c)");
695     
696     fOutCont->Add(fhBkgPt) ; 
697     fOutCont->Add(fhBkgRatioPt) ; 
698     fOutCont->Add(fhBkgDeltaPhi) ;
699     fOutCont->Add(fhBkgDeltaEta) ;
700     fOutCont->Add(fhBkgLeadingRatioPt) ;
701     fOutCont->Add(fhBkgLeadingDeltaPhi) ;
702     fOutCont->Add(fhBkgLeadingDeltaEta) ;
703     fOutCont->Add(fhBkgFFz) ;
704     fOutCont->Add(fhBkgFFxi) ;
705     fOutCont->Add(fhBkgFFpt) ;
706     fOutCont->Add(fhBkgNTracksInCone) ;
707     
708   }//not several cones
709   else{ //If we want to study the jet for different cones and pt
710     for(Int_t icone = 0; icone<fJetNCone; icone++){//icone
711       for(Int_t ipt = 0; ipt<fJetNPt;ipt++){ //ipt
712         
713         TString lastnamehist ="Cone"+ fJetNameCones[icone]+"Pt"+ fJetNamePtThres[ipt];
714         TString lastnametitle =", cone ="+fJetNameCones[icone]+", pt > " +fJetNamePtThres[ipt]+" GeV/c";
715         
716         //Jet Distributions
717         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); 
718         fhJetPts[icone][ipt]->SetYTitle("p_{T  jet}");
719         fhJetPts[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
720         
721         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); 
722         fhJetRatioPts[icone][ipt]->SetYTitle("p_{T  jet}/p_{T trigger}");
723         fhJetRatioPts[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
724         
725         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()); 
726         fhJetDeltaPhis[icone][ipt]->SetYTitle("#Delta #phi (rad)");
727         fhJetDeltaPhis[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
728         
729         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); 
730         fhJetDeltaEtas[icone][ipt]->SetYTitle("#Delta #eta");
731         fhJetDeltaEtas[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
732         
733         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); 
734         fhJetLeadingRatioPts[icone][ipt]->SetYTitle("p_{T  leading}/p_{T jet}");
735         fhJetLeadingRatioPts[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
736         
737         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()); 
738         fhJetLeadingDeltaPhis[icone][ipt]->SetYTitle("#Delta #phi (rad)");
739         fhJetLeadingDeltaPhis[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
740         
741         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); 
742         fhJetLeadingDeltaEtas[icone][ipt]->SetYTitle("#Delta #eta");
743         fhJetLeadingDeltaEtas[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
744         
745         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); 
746         fhJetFFzs[icone][ipt]->SetYTitle("z");
747         fhJetFFzs[icone][ipt]->SetXTitle("p_{T trigger}");
748         
749         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.); 
750         fhJetFFxis[icone][ipt]->SetYTitle("#xi");
751         fhJetFFxis[icone][ipt]->SetXTitle("p_{T trigger}");
752         
753         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.); 
754         fhJetFFpts[icone][ipt]->SetYTitle("p_{T charged hadron}");
755         fhJetFFpts[icone][ipt]->SetXTitle("p_{T trigger}");
756         
757         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); 
758         fhJetNTracksInCones[icone][ipt]->SetYTitle("N tracks in jet cone");
759         fhJetNTracksInCones[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
760         
761         fOutCont->Add(fhJetPts[icone][ipt]) ; 
762         fOutCont->Add(fhJetRatioPts[icone][ipt]) ; 
763         fOutCont->Add(fhJetDeltaPhis[icone][ipt]) ;
764         fOutCont->Add(fhJetDeltaEtas[icone][ipt]) ;
765         fOutCont->Add(fhJetLeadingRatioPts[icone][ipt]) ;
766         fOutCont->Add(fhJetLeadingDeltaPhis[icone][ipt]) ;
767         fOutCont->Add(fhJetLeadingDeltaEtas[icone][ipt]) ;
768         fOutCont->Add(fhJetFFzs[icone][ipt]) ;
769         fOutCont->Add(fhJetFFxis[icone][ipt]) ;
770         fOutCont->Add(fhJetFFpts[icone][ipt]) ;
771         fOutCont->Add(fhJetNTracksInCones[icone][ipt]) ;
772         
773         //Bkg Distributions
774         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); 
775         fhBkgPts[icone][ipt]->SetYTitle("p_{T  bkg}");
776         fhBkgPts[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
777         
778         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); 
779         fhBkgRatioPts[icone][ipt]->SetYTitle("p_{T  bkg}/p_{T trigger}");
780         fhBkgRatioPts[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
781         
782         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()); 
783         fhBkgDeltaPhis[icone][ipt]->SetYTitle("#Delta #phi (rad)");
784         fhBkgDeltaPhis[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
785         
786         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); 
787         fhBkgDeltaEtas[icone][ipt]->SetYTitle("#Delta #eta");
788         fhBkgDeltaEtas[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
789         
790         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); 
791         fhBkgLeadingRatioPts[icone][ipt]->SetYTitle("p_{T  leading}/p_{T bkg}");
792         fhBkgLeadingRatioPts[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
793         
794         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()); 
795         fhBkgLeadingDeltaPhis[icone][ipt]->SetYTitle("#Delta #phi (rad)");
796         fhBkgLeadingDeltaPhis[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
797         
798         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); 
799         fhBkgLeadingDeltaEtas[icone][ipt]->SetYTitle("#Delta #eta");
800         fhBkgLeadingDeltaEtas[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
801         
802         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); 
803         fhBkgFFzs[icone][ipt]->SetYTitle("z");
804         fhBkgFFzs[icone][ipt]->SetXTitle("p_{T trigger}");
805         
806         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.); 
807         fhBkgFFxis[icone][ipt]->SetYTitle("#xi");
808         fhBkgFFxis[icone][ipt]->SetXTitle("p_{T trigger}");
809         
810         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.); 
811         fhBkgFFpts[icone][ipt]->SetYTitle("p_{T charged hadron}");
812         fhBkgFFpts[icone][ipt]->SetXTitle("p_{T trigger}");
813         
814         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); 
815         fhBkgNTracksInCones[icone][ipt]->SetYTitle("N tracks in bkg cone");
816         fhBkgNTracksInCones[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
817         
818         fOutCont->Add(fhBkgPts[icone][ipt]) ; 
819         fOutCont->Add(fhBkgRatioPts[icone][ipt]) ; 
820         fOutCont->Add(fhBkgDeltaPhis[icone][ipt]) ;
821         fOutCont->Add(fhBkgDeltaEtas[icone][ipt]) ;
822         fOutCont->Add(fhBkgLeadingRatioPts[icone][ipt]) ;
823         fOutCont->Add(fhBkgLeadingDeltaPhis[icone][ipt]) ;
824         fOutCont->Add(fhBkgLeadingDeltaEtas[icone][ipt]) ;
825         fOutCont->Add(fhBkgFFzs[icone][ipt]) ;
826         fOutCont->Add(fhBkgFFxis[icone][ipt]) ;
827         fOutCont->Add(fhBkgFFpts[icone][ipt]) ;
828         fOutCont->Add(fhBkgNTracksInCones[icone][ipt]) ;
829         
830       }//ipt
831     } //icone
832   }//If we want to study any cone or pt threshold
833   
834   //Keep neutral meson selection histograms if requiered
835   //Setting done in AliNeutralMesonSelection
836   if(GetNeutralMesonSelection()){
837     TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
838     if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
839       for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) fOutCont->Add(nmsHistos->At(i)) ;
840         delete nmsHistos;
841   }
842   
843   
844   if(GetDebug() > 2){
845     printf("AliAnaParticleJetLeadingConeCorrelation::GetCreateOutputObjects() - All histograms names : \n");
846     for(Int_t i  = 0 ;  i<  fOutCont->GetEntries(); i++)
847       printf("Histo i %d name %s\n",i,((fOutCont->At(i))->GetName()));
848     //cout<< (fOutCont->At(i))->GetName()<<endl;
849   }
850   
851   return fOutCont;
852
853 }
854
855 //____________________________________________________________________________
856 Bool_t  AliAnaParticleJetLeadingConeCorrelation::GetLeadingParticle(AliAODPWG4ParticleCorrelation *particle, TLorentzVector & pLeading) 
857 {
858   //Search Charged or Neutral leading particle, select the highest one and fill AOD
859   
860   TLorentzVector pLeadingCh(0,0,0,0) ;
861   TLorentzVector pLeadingPi0(0,0,0,0) ;
862   
863   GetLeadingCharge(particle, pLeadingCh) ;
864   if(!fJetsOnlyInCTS) GetLeadingPi0(particle, pLeadingPi0) ;
865   
866   Double_t ptch = pLeadingCh.Pt(); 
867   Double_t ptpi = pLeadingPi0.Pt(); 
868   if (ptch > 0 || ptpi > 0){
869     if((ptch >= ptpi)){
870       if(GetDebug() > 1)printf("AliAnaParticleJetLeadingConeCorrelation::GetLeadingParticle() - Leading found in CTS \n");
871       pLeading = pLeadingCh;
872       if(GetDebug() > 1) printf("AliAnaParticleJetLeadingConeCorrelation::GetLeadingParticle() - Found Leading: pt %2.3f, phi %2.3f deg, eta %2.3f\n", 
873                                 pLeading.Pt(),pLeading.Phi()*TMath::RadToDeg(),pLeading.Eta()) ;
874       //Put leading in AOD
875       particle->SetLeading(pLeadingCh);
876       particle->SetLeadingDetector("CTS");
877       return kTRUE;
878     }
879     else{
880       if(GetDebug() > 1)
881         printf("AliAnaParticleJetLeadingConeCorrelation::GetLeadingParticle() - Leading found in EMCAL \n");
882       pLeading = pLeadingPi0;
883       if(GetDebug() > 1) printf("AliAnaParticleJetLeadingConeCorrelation::GetLeadingParticle() - Found Leading: pt %2.3f, phi %2.3f, eta %2.3f\n", 
884                                 pLeading.Pt(),pLeading.Phi()*TMath::RadToDeg(),pLeading.Eta()) ;
885       //Put leading in AOD
886       particle->SetLeading(pLeadingPi0);
887       particle->SetLeadingDetector("EMCAL");
888       return kTRUE;   
889     }
890   }  
891   
892   if(GetDebug() > 1)printf ("AliAnaParticleJetLeadingConeCorrelation::GetLeadingParticle() - NO LEADING PARTICLE FOUND \n");
893   
894   return kFALSE; 
895   
896 }
897
898 //____________________________________________________________________________
899 void  AliAnaParticleJetLeadingConeCorrelation::GetLeadingCharge(AliAODPWG4ParticleCorrelation* const particle, TLorentzVector & pLeading) const
900 {  
901   //Search for the charged particle with highest pt and with 
902   //Phi=Phi_trigger-Pi and pT=0.1E_gamma 
903   
904   if(GetAODCTS()){
905     Double_t ptTrig  = particle->Pt();
906     Double_t phiTrig = particle->Phi();
907     Double_t rat     = -100 ;
908     Double_t ptl     = -100 ;
909     Double_t phil    = -100 ;
910     Double_t pt      = -100.;
911     Double_t phi     = -100.;
912     TVector3 p3;
913     
914     for(Int_t ipr = 0;ipr < GetAODCTS()->GetEntriesFast() ; ipr ++ ){
915       AliAODTrack* track = (AliAODTrack *)(GetAODCTS()->At(ipr)) ;
916       p3.SetXYZ(track->Px(),track->Py(),track->Pz());
917       pt   = p3.Pt();
918       phi  = p3.Phi() ;
919       if(phi < 0) phi+=TMath::TwoPi();
920       rat  = pt/ptTrig ;
921       //printf("AliAnaParticleJetLeadingConeCorrelation::GetLeadingCharge() - Tracks: pt %2.3f eta %2.3f phi %2.3f pt/ptTrig %2.3f \n", 
922       //           pt, p3.Eta(), phi,pt/ptTrig) ; 
923       Float_t deltaphi = TMath::Abs(phiTrig-phi);
924       if((deltaphi > fDeltaPhiMinCut) && (deltaphi < fDeltaPhiMaxCut) &&
925          (rat > fLeadingRatioMinCut) && (rat < fLeadingRatioMaxCut)  && (pt  > ptl)) {
926         phil = phi ;
927         ptl  = pt ;
928         pLeading.SetVect(p3);
929       }
930     }// track loop
931     
932     if(GetDebug() > 1 && ptl > 0 ) printf("AliAnaParticleJetLeadingConeCorrelation::GetLeadingCharge() - Leading in CTS: pt %2.3f eta %2.3f phi %2.3f pt/ptTrig %2.3f, |phiTrig-phi| %2.3f \n", 
933                                           ptl, pLeading.Eta(), phil,ptl/ptTrig, TMath::Abs(phiTrig-phil)) ;
934     
935   }//CTS list exist
936 }
937
938 //____________________________________________________________________________
939 void  AliAnaParticleJetLeadingConeCorrelation::GetLeadingPi0(AliAODPWG4ParticleCorrelation* const particle, TLorentzVector & pLeading)
940 {  
941   //Search for the neutral pion with highest pt and with 
942   //Phi=Phi_trigger-Pi and pT=0.1E_gamma
943   
944   if(GetAODEMCAL()){
945     Double_t ptTrig  = particle->Pt();
946     Double_t phiTrig = particle->Phi();
947     Double_t rat     = -100 ;
948     Double_t ptl     = -100 ;
949     Double_t phil    = -100 ;
950     Double_t pt      = -100.;
951     Double_t phi     = -100.;
952     
953     TLorentzVector gammai;
954     TLorentzVector gammaj;
955           
956     //Get vertex for photon momentum calculation
957     Double_t vertex [] = {0,0,0} ; //vertex 
958     //Double_t vertex2[] = {0,0,0} ; //vertex of second input AOD 
959     if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) 
960     {
961       GetVertex(vertex);
962       //if(GetReader()->GetSecondInputAODTree()) GetReader()->GetSecondInputAODVertex(vertex2);
963     }
964           
965     //Cluster loop, select pairs with good pt, phi and fill AODs or histograms
966     for(Int_t iclus = 0;iclus < GetAODEMCAL()->GetEntriesFast() ; iclus ++ ){
967       AliVCluster * calo = (AliVCluster *)(GetAODEMCAL()->At(iclus)) ;
968       
969       //Input from second AOD?
970       //Int_t inputi = 0;
971       //          if     (particle->GetDetector() == "EMCAL" && GetReader()->GetAODEMCALNormalInputEntries() <= iclus) inputi = 1 ;
972       //          else if(particle->GetDetector() == "PHOS"  && GetReader()->GetAODPHOSNormalInputEntries()  <= iclus) inputi = 1;
973       
974       //Cluster selection, not charged, with photon or pi0 id and in fiducial cut
975       Int_t pdgi=0;
976       //if     (inputi == 0 && !SelectCluster(calo, vertex,  gammai, pdgi))  continue ;
977       //else if(inputi == 1 && !SelectCluster(calo, vertex2, gammai, pdgi))  continue ; 
978       if(!SelectCluster(calo, vertex,  gammai, pdgi))  continue ;
979       
980       if(GetDebug() > 2) printf("AliAnaParticleJetLeadingConeCorrelation::GetLeadingPi0() - Neutral cluster: pt %2.3f, phi %2.3f \n", 
981                                 gammai.Pt(),gammai.Phi());
982       
983       //2 gamma overlapped, found with PID
984       if(pdgi == AliCaloPID::kPi0){ 
985         
986         if(GetDebug() > 2) printf("AliAnaParticleJetLeadingConeCorrelation::GetLeadingPi0() - Neutral cluster ID as Pi0 \n");
987         
988         pt  = gammai.Pt();
989         rat = pt/ptTrig;
990         phi = gammai.Phi();
991         if(phi < 0) phi+=TMath::TwoPi();
992         
993         //Selection within angular and energy limits
994         Float_t deltaphi = TMath::Abs(phiTrig-phi);
995         if(pt > ptl  && rat > fLeadingRatioMinCut  && rat < fLeadingRatioMaxCut  && 
996            deltaphi > fDeltaPhiMinCut && deltaphi < fDeltaPhiMaxCut )
997         {
998           phil = phi ;
999           ptl  = pt ;
1000           pLeading.SetPxPyPzE(gammai.Px(),gammai.Py(),gammai.Pz(),gammai.E());
1001         }// cuts
1002       }// pdg = AliCaloPID::kPi0
1003       //Make invariant mass analysis
1004       else if(pdgi == AliCaloPID::kPhoton){     
1005         //Search the photon companion in case it comes from  a Pi0 decay
1006         //Apply several cuts to select the good pair
1007         for(Int_t jclus = iclus+1; jclus < GetAODEMCAL()->GetEntriesFast() ; jclus ++ ){
1008           AliVCluster * calo2 = (AliVCluster *) (GetAODEMCAL()->At(jclus)) ;
1009           
1010           //Input from second AOD?
1011           //Int_t inputj = 0;
1012           //      if     (particle->GetDetector() == "EMCAL" && GetReader()->GetAODEMCALNormalInputEntries() <= jclus) inputj = 1;
1013           //      else if(particle->GetDetector() == "PHOS"  && GetReader()->GetAODPHOSNormalInputEntries()  <= jclus) inputj = 1;
1014           
1015           //Cluster selection, not charged with photon or pi0 id and in fiducial cut
1016           Int_t pdgj=0;
1017           //if     (inputj == 0 && !SelectCluster(calo2, vertex,  gammaj, pdgj))  continue ;
1018           //else if(inputj == 1 && !SelectCluster(calo2, vertex2, gammaj, pdgj))  continue ;
1019           if     (!SelectCluster(calo2, vertex,  gammaj, pdgj))  continue ;
1020
1021           if(pdgj == AliCaloPID::kPhoton ){
1022             
1023             pt  = (gammai+gammaj).Pt();
1024             phi = (gammai+gammaj).Phi();
1025             if(phi < 0) phi+=TMath::TwoPi();  
1026             rat = pt/ptTrig;
1027             
1028             //Selection within angular and energy limits
1029             Float_t deltaphi = TMath::Abs(phiTrig-phi);  
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,(gammai+gammaj).Eta(), deltaphi, rat, (gammai+gammaj).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(gammai, gammaj)){
1037                 phil = phi ;
1038                 ptl  = pt ;
1039                 pLeading=(gammai+gammaj);       
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,(gammai+gammaj).Eta(), (gammai+gammaj).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 && pLeading.Pt() > 0 ) printf("AliAnaParticleJetLeadingConeCorrelation::GetLeadingPi0() - Leading EMCAL: pt %2.3f eta %2.3f phi %2.3f pt/Eg %2.3f \n",  
1053                                                     pLeading.Pt(), pLeading.Eta(), phil,  pLeading.Pt()/ptTrig) ;
1054     
1055   }//EMCAL list exists
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", GetAODCTS()->GetEntriesFast());
1299     printf("AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillAOD() - In EMCAL aod entries %d\n", GetAODEMCAL()->GetEntriesFast());
1300   }
1301   
1302   TLorentzVector pLeading(0,0,0,0); //It will contain the kinematics of the found leading particle
1303   
1304   //Loop on stored AOD particles, trigger
1305   Int_t naod = GetInputAODBranch()->GetEntriesFast();
1306   for(Int_t iaod = 0; iaod < naod ; iaod++){
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, pLeading)){
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, pLeading);
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     printf("AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillHistograms() - No input particles in AOD with name branch < %s > \n",
1333            GetInputAODName().Data());   
1334     abort();
1335   }
1336   if(GetDebug() > 1){
1337     printf("AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillHistograms() - Begin jet leading cone  correlation analysis, fill histograms \n");
1338     printf("AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillHistograms() - In particle branch aod entries %d\n", GetInputAODBranch()->GetEntriesFast());
1339     printf("AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillHistograms() - In CTS aod entries %d\n", GetAODCTS()->GetEntriesFast());
1340     printf("AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillHistograms() - In EMCAL aod entries %d\n", GetAODEMCAL()->GetEntriesFast());
1341   }
1342   
1343   TLorentzVector pLeading(0,0,0,0) ;
1344   
1345   //Loop on stored AOD particles, trigger
1346   Int_t naod = GetInputAODBranch()->GetEntriesFast();
1347   for(Int_t iaod = 0; iaod < naod ; iaod++){
1348     AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1349     
1350     if(OnlyIsolated() && !particle->IsIsolated()) continue;
1351     
1352     Double_t pt = particle->Pt();
1353     Double_t phi = particle->Phi();
1354     Double_t eta = particle->Eta();
1355     
1356     //Get leading particle, fill histograms
1357     pLeading = particle->GetLeading();
1358     TString det = particle->GetLeadingDetector();      
1359     
1360     if(det!="" && pLeading.Pt() > 0){
1361       Double_t ptL = pLeading.Pt(); 
1362       Double_t phiL = pLeading.Phi(); 
1363       if(phiL < 0 ) phiL+=TMath::TwoPi();
1364       Double_t etaL = pLeading.Eta(); 
1365       
1366       if(GetDebug() > 1)   printf("AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillHistograms() - Trigger with pt %3.2f, phi %2.2f, eta %2.2f\n", pt, phi, eta);       
1367
1368       if(det == "CTS"){
1369         fhChargedLeadingPt->Fill(pt,ptL);
1370         fhChargedLeadingPhi->Fill(pt,phiL);
1371         fhChargedLeadingEta->Fill(pt,etaL);
1372         fhChargedLeadingDeltaPt->Fill(pt,pt-ptL);
1373         fhChargedLeadingDeltaPhi->Fill(pt,TMath::Abs(phi-phiL));
1374         fhChargedLeadingDeltaEta->Fill(pt,eta-etaL);
1375         fhChargedLeadingRatioPt->Fill(pt,ptL/pt);
1376         fhChargedLeadingXi->Fill(pt,TMath::Log(pt/ptL));
1377         if(pt > 30) fhChargedLeadingDeltaPhiRatioPt30->Fill(TMath::Abs(phi-phiL),ptL/pt);
1378         if(pt > 50) fhChargedLeadingDeltaPhiRatioPt50->Fill(TMath::Abs(phi-phiL),ptL/pt);
1379       }
1380       else if(det== "EMCAL"){
1381         fhNeutralLeadingPt->Fill(pt,ptL);
1382         fhNeutralLeadingPhi->Fill(pt,phiL);
1383         fhNeutralLeadingEta->Fill(pt,etaL);
1384         fhNeutralLeadingDeltaPt->Fill(pt,pt-ptL);
1385         fhNeutralLeadingDeltaPhi->Fill(pt,TMath::Abs(phi-phiL));
1386         fhNeutralLeadingDeltaEta->Fill(pt,eta-etaL);
1387         fhNeutralLeadingRatioPt->Fill(pt,ptL/pt);
1388         fhNeutralLeadingXi->Fill(pt,TMath::Log(pt/ptL));
1389         if(pt > 30) fhNeutralLeadingDeltaPhiRatioPt30->Fill(TMath::Abs(phi-phiL),ptL/pt);
1390         if(pt > 50) fhNeutralLeadingDeltaPhiRatioPt50->Fill(TMath::Abs(phi-phiL),ptL/pt);
1391
1392       }
1393       
1394       //Fill Jet histograms
1395       TLorentzVector bkg(0,0,0,0);
1396       TLorentzVector jet(0,0,0,0);
1397       if(!fSeveralConeAndPtCuts){//just fill histograms      
1398         if(!fReMakeJet){
1399           jet=particle->GetCorrelatedJet();
1400           bkg=particle->GetCorrelatedBackground();
1401         }
1402         else  MakeJetFromAOD(particle, pLeading, jet,bkg);
1403         
1404         if(jet.Pt() > 0){//Jet was found
1405           FillJetHistos(particle, pLeading, jet,"Jet","");
1406           FillJetHistos(particle, pLeading, bkg,"Bkg","");
1407         }
1408       }
1409       else if(fSeveralConeAndPtCuts){
1410         for(Int_t icone = 0; icone<fJetNCone; icone++) {
1411           fJetCone=fJetCones[icone];      
1412           for(Int_t ipt = 0; ipt<fJetNPt;ipt++) {  
1413             TString lastname ="Cone"+ fJetNameCones[icone]+"Pt"+ fJetNamePtThres[ipt];
1414             fJetPtThreshold=fJetPtThres[ipt];
1415             MakeJetFromAOD(particle, pLeading, jet,bkg);
1416             if(jet.Pt() > 0) {//Jet was found
1417               FillJetHistos(particle, pLeading, jet,"Jet",lastname);
1418               FillJetHistos(particle, pLeading, bkg,"Bkg",lastname);
1419             }
1420           }//icone
1421         }//ipt
1422       }//fSeveralConeAndPtCuts
1423     }//Leading
1424   }//AOD trigger particle loop
1425   
1426   if(GetDebug() >1)printf("AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillHistograms() - End of jet leading cone analysis, fill histograms \n");
1427   
1428
1429
1430 //____________________________________________________________________________
1431 void AliAnaParticleJetLeadingConeCorrelation::MakeAODJet(AliAODPWG4ParticleCorrelation *particle, const TLorentzVector  pLeading)
1432   const {
1433   //Fill the jet with the particles around the leading particle with 
1434   //R=fJetCone and pt_th = fJetPtThres. Calculate the energy of the jet and 
1435   //fill aod with found information 
1436   
1437   TLorentzVector bkg(0,0,0,0);
1438   TLorentzVector jet(0,0,0,0);
1439   TLorentzVector lv (0,0,0,0); //Temporal container for jet particles kinematics 
1440   
1441   Double_t ptTrig   = particle->Pt();
1442   Double_t phiTrig  = particle->Phi();
1443   Double_t phil     = pLeading.Phi();
1444   if(phil<0) phil+=TMath::TwoPi();
1445   Double_t etal     = pLeading.Eta();
1446
1447   //Different pt cut for jet particles in different collisions systems
1448   Float_t ptcut = fJetPtThreshold;
1449   if(fPbPb && !fSeveralConeAndPtCuts && ptTrig > fPtTriggerSelectionCut)  ptcut = fJetPtThresPbPb ;
1450   
1451   //Add charged particles to jet if they are in cone around the leading particle
1452   if(!GetAODCTS()) {
1453     printf("AliAnaParticleJetLeadingConeCorrelation::MakeAODJet() - Cannot construct jets without tracks, STOP analysis");
1454     return;
1455   }
1456   
1457   //Fill jet with tracks
1458   TVector3 p3;
1459   //Initialize reference arrays that will contain jet and background tracks
1460   TObjArray * reftracks  = new TObjArray;
1461   TObjArray * reftracksbkg  = new TObjArray;
1462   
1463   for(Int_t ipr = 0;ipr < (GetAODCTS())->GetEntriesFast() ; ipr ++ ){
1464     AliAODTrack* track = (AliAODTrack *)((GetAODCTS())->At(ipr)) ;
1465     p3.SetXYZ(track->Px(),track->Py(),track->Pz());
1466     
1467     //Particles in jet 
1468     if(IsParticleInJetCone(p3.Eta(), p3.Phi(), etal, phil)){
1469       
1470       reftracks->Add(track); 
1471       
1472       if(p3.Pt() > ptcut ){
1473         lv.SetVect(p3);
1474         jet+=lv;
1475       }
1476     } 
1477
1478     //Background around (phi_gamma-pi, eta_leading)
1479     else if(IsParticleInJetCone(p3.Eta(),p3.Phi(),etal, phiTrig)) { 
1480       
1481       reftracksbkg->Add(track); 
1482
1483       if(p3.Pt() > ptcut ){
1484         lv.SetVect(p3);
1485         bkg+=lv;
1486       }
1487     }
1488   }//Track loop
1489   
1490   //Add referenced tracks to AOD
1491   if(reftracks->GetEntriesFast() > 0 ){
1492     reftracks->SetName(Form("%sTracks",GetAODObjArrayName().Data()));
1493     particle->AddObjArray(reftracks);
1494   }
1495   else  if(GetDebug() > 2 ) printf("AliAnaParticleJetLeadingConeCorrelation::MakeAODJet() - No tracks in jet cone\n");
1496   if(reftracksbkg->GetEntriesFast() > 0 ){
1497     reftracksbkg->SetName(Form("%sTracksBkg",GetAODObjArrayName().Data()));
1498     particle->AddObjArray(reftracksbkg);
1499   }
1500   else  if(GetDebug() > 2 ) printf("AliAnaParticleJetLeadingConeCorrelation::MakeAODJet() - No background tracks in jet cone\n");
1501   
1502   //Add neutral particles to jet
1503   //Initialize reference arrays that will contain jet and background tracks
1504   TObjArray * refclusters  = new TObjArray;
1505   TObjArray * refclustersbkg  = new TObjArray;
1506   if(!fJetsOnlyInCTS && GetAODEMCAL()){
1507     
1508         //Get vertex for photon momentum calculation
1509         Double_t vertex[]  = {0,0,0} ; //vertex 
1510         //Double_t vertex2[] = {0,0,0} ; //vertex of second input aod
1511         if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) 
1512         {
1513                 GetReader()->GetVertex(vertex);
1514                 //if(GetReader()->GetSecondInputAODTree()) GetReader()->GetSecondInputAODVertex(vertex2);
1515         }
1516           
1517     for(Int_t iclus = 0;iclus < (GetAODEMCAL())->GetEntriesFast() ; iclus ++ ){
1518       AliVCluster * calo = (AliVCluster *) (GetAODEMCAL()->At(iclus)) ;
1519       
1520       //Cluster selection, not charged
1521       if(IsTrackMatched(calo)) continue ;
1522       
1523           //Input from second AOD?
1524           Int_t input = 0;
1525 //        if     (particle->GetDetector() == "EMCAL" && GetReader()->GetAODEMCALNormalInputEntries() <= iclus) input = 1 ;
1526 //        else if(particle->GetDetector() == "PHOS"  && GetReader()->GetAODPHOSNormalInputEntries()  <= iclus) input = 1;
1527 //              
1528           //Get Momentum vector, 
1529           if     (input == 0) calo->GetMomentum(lv,vertex) ;//Assume that come from vertex in straight line
1530           //else if(input == 1) calo->GetMomentum(lv,vertex2);//Assume that come from vertex in straight line  
1531                 
1532       //Particles in jet 
1533       if(IsParticleInJetCone(lv.Eta(),lv.Phi(), etal, phil)){
1534
1535         refclusters->Add(calo); 
1536         
1537         if(lv.Pt() > ptcut )  jet+=lv;
1538       }
1539       //Background around (phi_gamma-pi, eta_leading)
1540       else if(IsParticleInJetCone(lv.Eta(),lv.Phi(),etal, phiTrig)){
1541         
1542         
1543         refclustersbkg->Add(calo); 
1544         
1545         if(lv.Pt() > ptcut )  bkg+=lv;
1546       }
1547     }//cluster loop
1548   }//jets with neutral particles
1549   
1550   //Add referenced clusters to AOD
1551   if(refclusters->GetEntriesFast() > 0 ){
1552     refclusters->SetName(Form("%sClusters",GetAODObjArrayName().Data()));
1553     particle->AddObjArray(refclusters);
1554   }
1555   else  if(GetDebug() > 2 ) printf("AliAnaParticleJetLeadingConeCorrelation::MakeAODJet() - No clusters in jet cone\n");
1556   if(refclustersbkg->GetEntriesFast() > 0 ){
1557     refclustersbkg->SetName(Form("%sClustersBkg",GetAODObjArrayName().Data()));
1558     particle->AddObjArray(refclustersbkg);
1559   }
1560   else if(GetDebug() > 2 ) printf("AliAnaParticleJetLeadingConeCorrelation::MakeAODJet() - No background clusters in jet cone\n");
1561   
1562   //If there is any jet found, select after some criteria and 
1563   //and fill AOD with corresponding TLorentzVector kinematics
1564   if(IsJetSelected(particle->Pt(), jet.Pt())) {
1565     particle->SetCorrelatedJet(jet);
1566     particle->SetCorrelatedBackground(bkg);
1567     if(GetDebug()>1) printf("AliAnaParticleJetLeadingConeCorrelation::MakeAODJet() - Found jet: Trigger  pt %2.3f, Jet pt %2.3f, Bkg pt %2.3f\n",ptTrig,jet.Pt(),bkg.Pt());
1568   }  
1569   
1570 }
1571
1572 //____________________________________________________________________________
1573 void AliAnaParticleJetLeadingConeCorrelation::MakeJetFromAOD(AliAODPWG4ParticleCorrelation *particle, const TLorentzVector  pLeading, TLorentzVector & jet, TLorentzVector & bkg)
1574   const {
1575   //Fill the jet with the particles around the leading particle with 
1576   //R=fJetCone and pt_th = fJetPtThres. Calculate the energy of the jet and 
1577   //fill aod tlorentzvectors with jet and bakcground found
1578   
1579   TLorentzVector lv (0,0,0,0); //Temporal container for jet particles kinematics 
1580   
1581   Double_t ptTrig  = particle->Pt();
1582   Double_t phiTrig = particle->Phi();
1583   Double_t phil = pLeading.Phi();
1584   if(phil < 0) phil+=TMath::TwoPi();
1585   Double_t etal = pLeading.Eta();
1586   
1587   TObjArray * refclusters    = particle->GetObjArray(Form("Clusters%s"   ,GetAODObjArrayName().Data()));
1588   TObjArray * reftracks      = particle->GetObjArray(Form("Tracks%s"     ,GetAODObjArrayName().Data()));
1589   TObjArray * refclustersbkg = particle->GetObjArray(Form("ClustersBkg%s",GetAODObjArrayName().Data()));
1590   TObjArray * reftracksbkg   = particle->GetObjArray(Form("TracksBkg%s"  ,GetAODObjArrayName().Data()));
1591   
1592   //Different pt cut for jet particles in different collisions systems
1593   Float_t ptcut = fJetPtThreshold;
1594   if(fPbPb && !fSeveralConeAndPtCuts && ptTrig > fPtTriggerSelectionCut)  ptcut = fJetPtThresPbPb ;
1595   
1596   //Fill jet with tracks
1597   //Particles in jet   
1598   TVector3 p3;
1599   if(reftracks){
1600     for(Int_t ipr = 0;ipr < reftracks->GetEntriesFast() ; ipr ++ ){
1601       AliAODTrack* track = (AliAODTrack *) reftracks->At(ipr) ;
1602       p3.SetXYZ(track->Px(),track->Py(),track->Pz());
1603       Float_t phi = p3.Phi();
1604       if(phi < 0) phi+=TMath::TwoPi();
1605       if(p3.Pt() > ptcut && IsParticleInJetCone(p3.Eta(), phi, etal, phil) ){
1606         lv.SetVect(p3);
1607         jet+=lv;
1608       }
1609     }//jet Track loop
1610   }
1611   //Particles in background   
1612   if(reftracksbkg){
1613     for(Int_t ipr = 0;ipr < reftracksbkg->GetEntriesFast() ; ipr ++ ){
1614       AliAODTrack* track = (AliAODTrack *) reftracksbkg->At(ipr) ;
1615       p3.SetXYZ(track->Px(),track->Py(),track->Pz());
1616       if(p3.Pt() > ptcut && IsParticleInJetCone(p3.Eta(),p3.Phi(),etal, phiTrig) ) {  
1617         lv.SetVect(p3);
1618         bkg+=lv;   
1619       }
1620     }//background Track loop
1621   }
1622   
1623   //Add neutral particles to jet
1624   if(!fJetsOnlyInCTS && refclusters){
1625     
1626         //Get vertex for photon momentum calculation
1627         Double_t vertex[]  = {0,0,0} ; //vertex 
1628         //Double_t vertex2[] = {0,0,0} ; //vertex of second input aod
1629         if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) 
1630         {
1631                 GetReader()->GetVertex(vertex);
1632                 //if(GetReader()->GetSecondInputAODTree()) GetReader()->GetSecondInputAODVertex(vertex2);
1633         }
1634           
1635     //Loop on jet particles
1636     if(refclusters){
1637       for(Int_t iclus = 0;iclus < refclusters->GetEntriesFast() ; iclus ++ ){
1638                   AliVCluster * calo = (AliVCluster *) refclusters->At(iclus) ;
1639                   
1640                   //Input from second AOD?
1641                   Int_t input = 0;
1642         //        if     (particle->GetDetector() == "EMCAL" && GetReader()->GetAODEMCALNormalInputEntries() <= iclus) input = 1 ;
1643 //                else if(particle->GetDetector() == "PHOS"  && GetReader()->GetAODPHOSNormalInputEntries()  <= iclus) input = 1;
1644                   
1645                   //Get Momentum vector, 
1646                   if     (input == 0) calo->GetMomentum(lv,vertex) ;//Assume that come from vertex in straight line
1647                   //else if(input == 1) calo->GetMomentum(lv,vertex2);//Assume that come from vertex in straight line  
1648                   
1649                   if(lv.Pt() > ptcut && IsParticleInJetCone(lv.Eta(),lv.Phi(), etal, phil)) jet+=lv;   
1650       }//jet cluster loop
1651     }
1652     
1653     //Loop on background particles
1654     if(refclustersbkg){
1655       for(Int_t iclus = 0;iclus < refclustersbkg->GetEntriesFast() ; iclus ++ ){
1656                   AliVCluster * calo = (AliVCluster *) refclustersbkg->At(iclus) ;
1657
1658                   //Input from second AOD?
1659                   Int_t input = 0;
1660 //                if     (particle->GetDetector() == "EMCAL" && GetReader()->GetAODEMCALNormalInputEntries() <= iclus) input = 1 ;
1661 //                else if(particle->GetDetector() == "PHOS"  && GetReader()->GetAODPHOSNormalInputEntries()  <= iclus) input = 1;
1662                   
1663                   //Get Momentum vector, 
1664                   if     (input == 0) calo->GetMomentum(lv,vertex) ;//Assume that come from vertex in straight line
1665                   //else if(input == 1) calo->GetMomentum(lv,vertex2);//Assume that come from vertex in straight line  
1666                   
1667                   if( lv.Pt() > ptcut && IsParticleInJetCone(lv.Eta(),lv.Phi(),etal, phiTrig)) bkg+=lv;
1668       }//background cluster loop 
1669     }
1670   }//clusters in jet
1671   
1672   //If there is any jet found, leave jet and bkg as they are, 
1673   //if not set them to 0.
1674   if(!IsJetSelected(particle->Pt(), jet.Pt())) {
1675     jet.SetPxPyPzE(0.,0.,0.,0.);
1676     bkg.SetPxPyPzE(0.,0.,0.,0.);
1677   }
1678   else
1679     if(GetDebug()>1) printf("AliAnaParticleJetLeadingConeCorrelation::MakeJetFromAOD() - Found jet: Trigger  pt %2.3f, Jet pt %2.3f, Bkg pt %2.3f\n",ptTrig,jet.Pt(),bkg.Pt());
1680   
1681 }
1682
1683 //____________________________________________________________________________
1684 //Bool_t  AliAnaParticleJetLeadingConeCorrelation::SelectCluster(AliVCluster * calo, Double_t *vertex, TLorentzVector & mom, Int_t & pdg) {
1685 //  //Select cluster depending on its pid and acceptance selections
1686 //  
1687 //  //Skip matched clusters with tracks
1688 //  if(IsTrackMatched(calo)) return kFALSE;
1689 //  
1690 //  //Check PID
1691 //  calo->GetMomentum(mom,vertex);//Assume that come from vertex in straight line
1692 //  pdg = AliCaloPID::kPhoton;   
1693 //  if(IsCaloPIDOn()){
1694 //    //Get most probable PID, 2 options check PID weights (in MC this option is mandatory)
1695 //    //or redo PID, recommended option for EMCal.
1696 //    if(!IsCaloPIDRecalculationOn() || GetReader()->GetDataType() == AliCaloTrackReader::kMC )
1697 //      pdg = GetCaloPID()->GetPdg("EMCAL",calo->GetPID(),mom.E());//PID with weights
1698 //    else
1699 //      pdg = GetCaloPID()->GetPdg("EMCAL",mom,calo);//PID recalculated
1700 //    
1701 //  //  if(GetDebug() > 3) printf("AliAnaParticleJetLeadingConeCorrelation::SelectCluster() - PDG of identified particle %d\n",pdg);
1702 //    //If it does not pass pid, skip
1703 //    if(pdg != AliCaloPID::kPhoton && pdg != AliCaloPID::kPi0) 
1704 //      return kFALSE ;
1705 //  }//CaloPID
1706 //  
1707 //   //Check acceptance selection
1708 //  if(IsFiducialCutOn()){
1709 //    Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,"EMCAL") ;
1710 //    if(! in ) return kFALSE ;
1711 //  }
1712 //  
1713 //  //if(GetDebug() > 3) printf("AliAnaParticleJetLeadingConeCorrelation::SelectCluster() - Cluster selection cuts passed: pT %3.2f, pdg %d\n",mom.Pt(), pdg);
1714 //  
1715 //  return kTRUE; 
1716 //  
1717 //}
1718
1719 //__________________________________________________________________
1720 void AliAnaParticleJetLeadingConeCorrelation::Print(const Option_t * opt) const
1721 {
1722   
1723   //Print some relevant parameters set for the analysis
1724   if(! opt)
1725     return;
1726   
1727   printf("**** Print  %s %s ****\n", GetName(), GetTitle() ) ;
1728   AliAnaPartCorrBaseClass::Print(" ");
1729   
1730   if(fJetsOnlyInCTS)printf("Jets reconstructed in CTS \n");
1731   else printf("Jets reconstructed in CTS+EMCAL \n");
1732   
1733   if(fPbPb) printf("PbPb events, pT cut in jet cone energy reconstruction %2.1f \n", fJetPtThreshold);
1734   else printf("pp events, pT cut in jet cone energy reconstruction %2.1f \n", fJetPtThresPbPb);
1735   
1736   printf("If pT of trigger < %2.3f, select jets as in pp? \n", fPtTriggerSelectionCut);
1737   
1738   printf("Phi gamma-Leading        <     %3.2f\n", fDeltaPhiMaxCut) ; 
1739   printf("Phi gamma-Leading        >     %3.2f\n", fDeltaPhiMinCut) ;
1740   printf("pT Leading / pT Trigger  <     %3.2f\n", fLeadingRatioMaxCut) ; 
1741   printf("pT Leading / pT Trigger  >     %3.2f\n", fLeadingRatioMinCut) ;
1742   
1743   if(fSelect == 2){
1744     printf("pT Jet / pT Gamma                     <    %3.2f\n", fJetRatioMaxCut) ; 
1745     printf("pT Jet / pT Gamma                     >    %3.2f\n", fJetRatioMinCut) ;
1746     printf("pT Jet (Only CTS)/ pT Trigger   <    %3.2f\n", fJetCTSRatioMaxCut) ; 
1747     printf("pT Jet (Only CTS)/ pT Trigger   >    %3.2f\n", fJetCTSRatioMinCut) ;
1748   }
1749   else if(fSelect == 0)
1750     printf("Accept all reconstructed jets \n") ;
1751   else   if(fSelect == 1)
1752     printf("Accept jets depending on trigger energy \n") ;
1753   else 
1754     printf("Wrong jet selection option:   %d \n", fSelect) ;
1755   
1756   printf("Isolated Trigger?  %d\n", fSelectIsolated) ;
1757
1758