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