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