]>
Commit | Line | Data |
---|---|---|
bdcfac30 | 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 | /* History of cvs commits: | |
18 | * | |
19 | * $Log$ | |
3788def4 | 20 | * Revision 1.3 2007/10/29 13:48:42 gustavo |
21 | * Corrected coding violations | |
22 | * | |
3bb2c538 | 23 | * Revision 1.2 2007/08/17 12:40:04 schutz |
24 | * New analysis classes by Gustavo Conesa | |
25 | * | |
bdcfac30 | 26 | * Revision 1.1.2.1 2007/07/26 10:32:09 schutz |
27 | * new analysis classes in the the new analysis framework | |
28 | * | |
29 | * | |
30 | */ | |
31 | ||
32 | //_________________________________________________________________________ | |
3bb2c538 | 33 | // Class for the analysis of gamma-jet correlations: |
34 | // 1)Take the prompt photon found with AliAnaGammaDirect, | |
35 | // 2) Search for the highest pt leading particle opposite to the photon within a phi, pt window | |
36 | // 3) Take all particles around leading in a cone R with pt larger than threshold and construct the jet | |
bdcfac30 | 37 | // |
38 | // Class created from old AliPHOSGammaJet | |
39 | // (see AliRoot versions previous Release 4-09) | |
40 | // | |
41 | //*-- Author: Gustavo Conesa (LNF-INFN) | |
42 | ////////////////////////////////////////////////////////////////////////////// | |
43 | ||
44 | ||
45 | // --- ROOT system --- | |
bdcfac30 | 46 | |
47 | //---- AliRoot system ---- | |
bdcfac30 | 48 | #include "AliNeutralMesonSelection.h" |
49 | #include "AliAnaGammaJetLeadCone.h" | |
3bb2c538 | 50 | #include "AliLog.h" |
bdcfac30 | 51 | |
52 | ClassImp(AliAnaGammaJetLeadCone) | |
53 | ||
54 | ||
55 | //____________________________________________________________________________ | |
56 | AliAnaGammaJetLeadCone::AliAnaGammaJetLeadCone() : | |
57 | AliAnaGammaCorrelation(), fPbPb(kFALSE), | |
58 | fSeveralConeAndPtCuts(0), | |
59 | fEtaEMCALCut(0.), | |
60 | fJetCTSRatioMaxCut(0.), fJetCTSRatioMinCut(0.), fJetRatioMaxCut(0.), | |
61 | fJetRatioMinCut(0.), | |
62 | fJetNCone(0),fJetNPt(0), fJetCone(0), | |
63 | fJetPtThreshold(0),fJetPtThresPbPb(0), | |
64 | fPtJetSelectionCut(0.0), fSelect(0), | |
65 | fhPhiCharged(0), fhPhiNeutral(0), fhEtaCharged(0), fhEtaNeutral(0), | |
66 | fhDeltaPhiGammaCharged(0), fhDeltaPhiGammaNeutral(0), | |
67 | fhDeltaEtaGammaCharged(0), fhDeltaEtaGammaNeutral(0), | |
68 | fhAnglePairLeading(), fhInvMassPairLeading(), | |
69 | fhChargedRatio(0), fhNeutralRatio (0), | |
70 | fhNBkg (0), fhNLeading(0), fhNJet(0), fhJetRatio(0), fhJetPt (0), | |
71 | fhBkgRatio (0), fhBkgPt(0), fhJetFragment(0), fhBkgFragment(0), | |
72 | fhJetPtDist(0), fhBkgPtDist(0) | |
73 | { | |
74 | //Default Ctor | |
75 | ||
76 | //Initialize parameters | |
77 | ||
78 | SetCorrelationType(kJetLeadCone); | |
79 | ||
80 | for(Int_t i = 0; i<10; i++){ | |
81 | fJetCones[i] = 0.0 ; | |
82 | fJetNameCones[i] = "" ; | |
83 | fJetPtThres[i] = 0.0 ; | |
84 | fJetNamePtThres[i] = "" ; | |
85 | if( i < 6 ){ | |
86 | fJetXMin1[i] = 0.0 ; | |
87 | fJetXMin2[i] = 0.0 ; | |
88 | fJetXMax1[i] = 0.0 ; | |
89 | fJetXMax2[i] = 0.0 ; | |
90 | fBkgMean[i] = 0.0 ; | |
91 | fBkgRMS[i] = 0.0 ; | |
92 | if( i < 2 ){ | |
93 | fJetE1[i] = 0.0 ; | |
94 | fJetE2[i] = 0.0 ; | |
95 | fJetSigma1[i] = 0.0 ; | |
96 | fJetSigma2[i] = 0.0 ; | |
97 | fPhiEMCALCut[i] = 0.0 ; | |
98 | } | |
99 | } | |
100 | } | |
101 | ||
102 | InitParameters(); | |
103 | } | |
104 | ||
105 | //____________________________________________________________________________ | |
106 | AliAnaGammaJetLeadCone::AliAnaGammaJetLeadCone(const AliAnaGammaJetLeadCone & g) : | |
107 | AliAnaGammaCorrelation(g), fPbPb(g.fPbPb), | |
108 | fSeveralConeAndPtCuts(g.fSeveralConeAndPtCuts), | |
109 | fEtaEMCALCut(g.fEtaEMCALCut), | |
110 | fJetCTSRatioMaxCut(g.fJetCTSRatioMaxCut), | |
111 | fJetCTSRatioMinCut(g.fJetCTSRatioMinCut), fJetRatioMaxCut(g.fJetRatioMaxCut), | |
112 | fJetRatioMinCut(g.fJetRatioMinCut), fJetNCone(g.fJetNCone), | |
113 | fJetNPt(g.fJetNPt), fJetCone(g.fJetCone), | |
114 | fJetPtThreshold(g.fJetPtThreshold),fJetPtThresPbPb(g.fJetPtThresPbPb), | |
115 | fPtJetSelectionCut(g.fPtJetSelectionCut), fSelect(g.fSelect), | |
116 | fhPhiCharged(g.fhPhiCharged), fhPhiNeutral(g.fhPhiNeutral), | |
117 | fhEtaCharged(g.fhEtaCharged), fhEtaNeutral(g.fhEtaNeutral), | |
118 | fhDeltaPhiGammaCharged(g.fhDeltaPhiGammaCharged), | |
119 | fhDeltaPhiGammaNeutral(g.fhDeltaPhiGammaNeutral), | |
120 | fhDeltaEtaGammaCharged(g.fhDeltaEtaGammaCharged), | |
121 | fhDeltaEtaGammaNeutral(g.fhDeltaEtaGammaNeutral), | |
122 | fhAnglePairLeading(g.fhAnglePairLeading), fhInvMassPairLeading(g.fhInvMassPairLeading), | |
123 | fhChargedRatio(g.fhChargedRatio), fhNeutralRatio(g.fhNeutralRatio), | |
124 | fhNBkg(g. fhNBkg), fhNLeading(g. fhNLeading), fhNJet(g.fhNJet), fhJetRatio(g.fhJetRatio), fhJetPt(g.fhJetPt), | |
125 | fhBkgRatio (g.fhBkgRatio), fhBkgPt(g.fhBkgPt), fhJetFragment(g.fhJetFragment), fhBkgFragment(g.fhBkgFragment), | |
126 | fhJetPtDist(g.fhJetPtDist), fhBkgPtDist(g.fhBkgPtDist) | |
127 | { | |
128 | // cpy ctor | |
129 | ||
130 | for(Int_t i = 0; i<10; i++){ | |
131 | fJetCones[i] = g.fJetCones[i] ; | |
132 | fJetNameCones[i] = g.fJetNameCones[i] ; | |
133 | fJetPtThres[i] = g.fJetPtThres[i] ; | |
134 | fJetNamePtThres[i] = g.fJetNamePtThres[i] ; | |
135 | if( i < 6 ){ | |
136 | fJetXMin1[i] = g.fJetXMin1[i] ; | |
137 | fJetXMin2[i] = g.fJetXMin2[i] ; | |
138 | fJetXMax1[i] = g.fJetXMax1[i] ; | |
139 | fJetXMax2[i] = g.fJetXMax2[i] ; | |
140 | fBkgMean[i] = g.fBkgMean[i] ; | |
141 | fBkgRMS[i] = g.fBkgRMS[i] ; | |
142 | if( i < 2 ){ | |
143 | fJetE1[i] = g.fJetE1[i] ; | |
144 | fJetE2[i] = g.fJetE2[i] ; | |
145 | fJetSigma1[i] = g.fJetSigma1[i] ; | |
146 | fJetSigma2[i] = g.fJetSigma2[i] ; | |
147 | fPhiEMCALCut[i] = g.fPhiEMCALCut[i] ; | |
148 | } | |
149 | } | |
150 | } | |
151 | ||
152 | ||
153 | } | |
154 | ||
155 | //_________________________________________________________________________ | |
156 | AliAnaGammaJetLeadCone & AliAnaGammaJetLeadCone::operator = (const AliAnaGammaJetLeadCone & source) | |
157 | { | |
158 | // assignment operator | |
159 | ||
160 | if(this == &source)return *this; | |
161 | ((AliAnaGammaCorrelation *)this)->operator=(source); | |
162 | ||
163 | fSeveralConeAndPtCuts = source.fSeveralConeAndPtCuts ; | |
164 | fPbPb = source.fPbPb ; | |
165 | fEtaEMCALCut = source.fEtaEMCALCut ; | |
166 | fJetCTSRatioMaxCut = source.fJetCTSRatioMaxCut ; | |
167 | fJetCTSRatioMinCut = source.fJetCTSRatioMinCut ; fJetRatioMaxCut = source.fJetRatioMaxCut ; | |
168 | fJetRatioMinCut = source.fJetRatioMinCut ; fJetNCone = source.fJetNCone ; | |
169 | fJetNPt = source.fJetNPt ; fJetCone = source.fJetCone ; | |
170 | fJetPtThreshold = source.fJetPtThreshold ; | |
171 | fJetPtThresPbPb = source.fJetPtThresPbPb ; | |
172 | fPtJetSelectionCut = source.fPtJetSelectionCut ; | |
173 | fSelect = source.fSelect ; fhChargedRatio = source.fhChargedRatio ; fhNeutralRatio = source.fhNeutralRatio ; | |
174 | ||
175 | fhPhiCharged = source.fhPhiCharged ; fhPhiNeutral = source.fhPhiNeutral ; | |
176 | fhEtaCharged = source.fhEtaCharged ; fhEtaNeutral = source.fhEtaNeutral ; | |
177 | fhDeltaPhiGammaCharged = source.fhDeltaPhiGammaCharged ; | |
178 | fhDeltaPhiGammaNeutral = source.fhDeltaPhiGammaNeutral ; | |
179 | fhDeltaEtaGammaCharged = source.fhDeltaEtaGammaCharged ; | |
180 | fhDeltaEtaGammaNeutral = source.fhDeltaEtaGammaNeutral ; | |
181 | ||
182 | fhAnglePairLeading = source.fhAnglePairLeading ; | |
183 | fhInvMassPairLeading = source.fhInvMassPairLeading ; | |
184 | fhNBkg = source. fhNBkg ; fhNLeading = source. fhNLeading ; | |
185 | fhNJet = source.fhNJet ; fhJetRatio = source.fhJetRatio ; fhJetPt = source.fhJetPt ; | |
186 | fhBkgRatio = source.fhBkgRatio ; fhBkgPt = source.fhBkgPt ; | |
187 | fhJetFragment = source.fhJetFragment ; fhBkgFragment = source.fhBkgFragment ; | |
188 | fhJetPtDist = source.fhJetPtDist ; fhBkgPtDist = source.fhBkgPtDist ; | |
189 | ||
190 | ||
191 | for(Int_t i = 0; i<10; i++){ | |
192 | fJetCones[i] = source.fJetCones[i] ; | |
193 | fJetNameCones[i] = source.fJetNameCones[i] ; | |
194 | fJetPtThres[i] = source.fJetPtThres[i] ; | |
195 | fJetNamePtThres[i] = source.fJetNamePtThres[i] ; | |
196 | if( i < 6 ){ | |
197 | fJetXMin1[i] = source.fJetXMin1[i] ; | |
198 | fJetXMin2[i] = source.fJetXMin2[i] ; | |
199 | fJetXMax1[i] = source.fJetXMax1[i] ; | |
200 | fJetXMax2[i] = source.fJetXMax2[i] ; | |
201 | fBkgMean[i] = source.fBkgMean[i] ; | |
202 | fBkgRMS[i] = source.fBkgRMS[i] ; | |
203 | if( i < 2 ){ | |
204 | fJetE1[i] = source.fJetE1[i] ; | |
205 | fJetE2[i] = source.fJetE2[i] ; | |
206 | fJetSigma1[i] = source.fJetSigma1[i] ; | |
207 | fJetSigma2[i] = source.fJetSigma2[i] ; | |
208 | fPhiEMCALCut[i] = source.fPhiEMCALCut[i] ; | |
209 | } | |
210 | } | |
211 | } | |
212 | ||
213 | return *this; | |
214 | ||
215 | } | |
216 | ||
217 | //____________________________________________________________________________ | |
218 | AliAnaGammaJetLeadCone::~AliAnaGammaJetLeadCone() | |
219 | { | |
3788def4 | 220 | // Remove all pointers except analysis output pointers. |
bdcfac30 | 221 | |
bdcfac30 | 222 | } |
223 | ||
224 | ||
225 | ||
226 | //________________________________________________________________________ | |
227 | TList * AliAnaGammaJetLeadCone::GetCreateOutputObjects() | |
228 | { | |
229 | // Create histograms to be saved in output file and | |
230 | // store them in fOutputContainer | |
231 | ||
232 | AliDebug(1,"Init jet in leading cone histograms"); | |
233 | ||
234 | TList * outputContainer = new TList() ; | |
235 | outputContainer->SetName("GammaJetCorrelationHistos") ; | |
236 | ||
237 | fhChargedRatio = new TH2F | |
238 | ("ChargedRatio","p_{T leading charge} /p_{T #gamma} vs p_{T #gamma}", | |
239 | 120,0,120,120,0,1); | |
240 | fhChargedRatio->SetYTitle("p_{T lead charge} /p_{T #gamma}"); | |
241 | fhChargedRatio->SetXTitle("p_{T #gamma} (GeV/c)"); | |
242 | ||
243 | fhPhiCharged = new TH2F | |
244 | ("PhiCharged","#phi_{#pi^{#pm}} vs p_{T #gamma}", | |
245 | 120,0,120,120,0,7); | |
246 | fhPhiCharged->SetYTitle("#phi_{#pi^{#pm}} (rad)"); | |
247 | fhPhiCharged->SetXTitle("p_{T #gamma} (GeV/c)"); | |
248 | ||
249 | fhEtaCharged = new TH2F | |
250 | ("EtaCharged","#eta_{#pi^{#pm}} vs p_{T #gamma}", | |
251 | 120,0,120,120,-1,1); | |
252 | fhEtaCharged->SetYTitle("#eta_{#pi^{#pm}} (rad)"); | |
253 | fhEtaCharged->SetXTitle("p_{T #gamma} (GeV/c)"); | |
254 | ||
255 | fhDeltaPhiGammaCharged = new TH2F | |
256 | ("DeltaPhiGammaCharged","#phi_{#gamma} - #phi_{charged #pi} vs p_{T #gamma}", | |
257 | 200,0,120,200,0,6.4); | |
258 | fhDeltaPhiGammaCharged->SetYTitle("#Delta #phi"); | |
259 | fhDeltaPhiGammaCharged->SetXTitle("p_{T #gamma} (GeV/c)"); | |
260 | ||
261 | fhDeltaEtaGammaCharged = new TH2F | |
262 | ("DeltaEtaGammaCharged","#eta_{#gamma} - #eta_{#pi^{#pm}} vs p_{T #gamma}", | |
263 | 200,0,120,200,-2,2); | |
264 | fhDeltaEtaGammaCharged->SetYTitle("#Delta #eta"); | |
265 | fhDeltaEtaGammaCharged->SetXTitle("p_{T #gamma} (GeV/c)"); | |
266 | ||
267 | outputContainer->Add(fhPhiCharged) ; | |
268 | outputContainer->Add(fhEtaCharged) ; | |
269 | outputContainer->Add(fhChargedRatio) ; | |
270 | outputContainer->Add(fhDeltaPhiGammaCharged) ; | |
271 | outputContainer->Add(fhDeltaEtaGammaCharged) ; | |
272 | ||
273 | if(!AreJetsOnlyInCTS()){ | |
274 | ||
275 | fhNeutralRatio = new TH2F | |
276 | ("NeutralRatio","p_{T leading #pi^{0}} /p_{T #gamma} vs p_{T #gamma}", | |
277 | 120,0,120,120,0,1); | |
278 | fhNeutralRatio->SetYTitle("p_{T lead #pi^{0}} /p_{T #gamma}"); | |
279 | fhNeutralRatio->SetXTitle("p_{T #gamma} (GeV/c)"); | |
280 | ||
281 | // | |
282 | fhAnglePairLeading = new TH2F | |
283 | ("AnglePairLeading", | |
284 | "Angle between all #gamma pair finally selected vs p_{T #pi^{0}}", | |
285 | 200,0,50,200,0,0.2); | |
286 | fhAnglePairLeading->SetYTitle("Angle (rad)"); | |
287 | fhAnglePairLeading->SetXTitle("E_{ #pi^{0}} (GeV/c)"); | |
288 | ||
289 | fhInvMassPairLeading = new TH2F | |
290 | ("InvMassPairLeading", | |
291 | "Invariant Mass of #gamma pair selected vs p_{T #gamma}", | |
292 | 120,0,120,360,0,0.5); | |
293 | fhInvMassPairLeading->SetYTitle("Invariant Mass (GeV/c^{2})"); | |
294 | fhInvMassPairLeading->SetXTitle("p_{T #gamma} (GeV/c)"); | |
295 | ||
296 | fhPhiNeutral = new TH2F | |
297 | ("PhiNeutral","#phi_{#pi^{0}} vs p_{T #gamma}", | |
298 | 120,0,120,120,0,7); | |
299 | fhPhiNeutral->SetYTitle("#phi_{#pi^{0}} (rad)"); | |
300 | fhPhiNeutral->SetXTitle("p_{T #gamma} (GeV/c)"); | |
301 | ||
302 | fhEtaNeutral = new TH2F | |
303 | ("EtaNeutral","#eta_{#pi^{0}} vs p_{T #gamma}", | |
304 | 120,0,120,120,-1,1); | |
305 | fhEtaNeutral->SetYTitle("#eta_{#pi^{0}} (rad)"); | |
306 | fhEtaNeutral->SetXTitle("p_{T #gamma} (GeV/c)"); | |
307 | ||
308 | fhDeltaPhiGammaNeutral = new TH2F | |
309 | ("DeltaPhiGammaNeutral","#phi_{#gamma} - #phi_{#pi^{0}} vs p_{T #gamma}", | |
310 | 200,0,120,200,0,6.4); | |
311 | fhDeltaPhiGammaNeutral->SetYTitle("#Delta #phi"); | |
312 | fhDeltaPhiGammaNeutral->SetXTitle("p_{T #gamma} (GeV/c)"); | |
313 | ||
314 | fhDeltaEtaGammaNeutral = new TH2F | |
315 | ("DeltaEtaGammaNeutral","#eta_{#gamma} - #eta_{#pi^{#pm}} vs p_{T #gamma}", | |
316 | 200,0,120,200,-2,2); | |
317 | fhDeltaEtaGammaNeutral->SetYTitle("#Delta #eta"); | |
318 | fhDeltaEtaGammaNeutral->SetXTitle("p_{T #gamma} (GeV/c)"); | |
319 | ||
320 | outputContainer->Add(fhPhiNeutral) ; | |
321 | outputContainer->Add(fhEtaNeutral) ; | |
322 | outputContainer->Add(fhNeutralRatio) ; | |
323 | outputContainer->Add(fhDeltaPhiGammaNeutral) ; | |
324 | outputContainer->Add(fhDeltaEtaGammaNeutral) ; | |
325 | ||
326 | outputContainer->Add(fhInvMassPairLeading) ; | |
327 | outputContainer->Add(fhAnglePairLeading) ; | |
328 | } | |
329 | ||
330 | if(!fSeveralConeAndPtCuts){// not several cones | |
331 | ||
332 | //Count | |
333 | fhNBkg = new TH1F("NBkg","bkg multiplicity",9000,0,9000); | |
334 | fhNBkg->SetYTitle("counts"); | |
335 | fhNBkg->SetXTitle("N"); | |
336 | outputContainer->Add(fhNBkg) ; | |
337 | ||
338 | fhNLeading = new TH2F | |
339 | ("NLeading","Accepted Jet Leading", 240,0,120,240,0,120); | |
340 | fhNLeading->SetYTitle("p_{T charge} (GeV/c)"); | |
341 | fhNLeading->SetXTitle("p_{T #gamma}(GeV/c)"); | |
342 | outputContainer->Add(fhNLeading) ; | |
343 | ||
344 | fhNJet = new TH1F("NJet","Accepted jets",240,0,120); | |
345 | fhNJet->SetYTitle("N"); | |
346 | fhNJet->SetXTitle("p_{T #gamma}(GeV/c)"); | |
347 | outputContainer->Add(fhNJet) ; | |
348 | ||
349 | //Ratios and Pt dist of reconstructed (not selected) jets | |
350 | //Jet | |
351 | fhJetRatio = new TH2F | |
352 | ("JetRatio","p_{T jet lead}/p_{T #gamma} vs p_{T #gamma}", | |
353 | 240,0,120,200,0,10); | |
354 | fhJetRatio->SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}"); | |
355 | fhJetRatio->SetXTitle("p_{T #gamma} (GeV/c)"); | |
356 | outputContainer->Add(fhJetRatio) ; | |
357 | ||
358 | fhJetPt = new TH2F | |
359 | ("JetPt", "p_{T jet lead} vs p_{T #gamma}",240,0,120,400,0,200); | |
360 | fhJetPt->SetYTitle("p_{T jet}"); | |
361 | fhJetPt->SetXTitle("p_{T #gamma} (GeV/c)"); | |
362 | outputContainer->Add(fhJetPt) ; | |
363 | ||
364 | //Bkg | |
365 | ||
366 | fhBkgRatio = new TH2F | |
367 | ("BkgRatio","p_{T bkg lead}/p_{T #gamma} vs p_{T #gamma}", | |
368 | 240,0,120,200,0,10); | |
369 | fhBkgRatio->SetYTitle("p_{T bkg lead charge}/p_{T #gamma}"); | |
370 | fhBkgRatio->SetXTitle("p_{T #gamma} (GeV/c)"); | |
371 | outputContainer->Add(fhBkgRatio) ; | |
372 | ||
373 | fhBkgPt = new TH2F | |
374 | ("BkgPt","p_{T jet lead} vs p_{T #gamma}",240,0,120,400,0,200); | |
375 | fhBkgPt->SetYTitle("p_{T jet lead charge}/p_{T #gamma}"); | |
376 | fhBkgPt->SetXTitle("p_{T #gamma} (GeV/c)"); | |
377 | outputContainer->Add(fhBkgPt) ; | |
378 | ||
379 | //Jet Distributions | |
380 | ||
381 | fhJetFragment = | |
382 | new TH2F("JetFragment","x = p_{T i charged}/p_{T #gamma}", | |
383 | 240,0.,120.,1000,0.,1.2); | |
384 | fhJetFragment->SetYTitle("x_{T}"); | |
385 | fhJetFragment->SetXTitle("p_{T #gamma}"); | |
386 | outputContainer->Add(fhJetFragment) ; | |
387 | ||
388 | fhBkgFragment = new TH2F | |
389 | ("BkgFragment","x = p_{T i charged}/p_{T #gamma}", | |
390 | 240,0.,120.,1000,0.,1.2); | |
391 | fhBkgFragment->SetYTitle("x_{T}"); | |
392 | fhBkgFragment->SetXTitle("p_{T #gamma}"); | |
393 | outputContainer->Add(fhBkgFragment) ; | |
394 | ||
395 | fhJetPtDist = | |
396 | new TH2F("JetPtDist","x = p_{T i charged}",240,0.,120.,400,0.,200.); | |
397 | fhJetPtDist->SetXTitle("p_{T #gamma} (GeV/c)"); | |
398 | outputContainer->Add(fhJetPtDist) ; | |
399 | ||
400 | fhBkgPtDist = new TH2F | |
401 | ("BkgPtDist","x = p_{T i charged}",240,0.,120.,400,0.,200.); | |
402 | fhBkgPtDist->SetXTitle("p_{T #gamma} (GeV/c)"); | |
403 | outputContainer->Add(fhBkgPtDist) ; | |
404 | ||
405 | }//not several cones | |
406 | else{ //If we want to study the jet for different cones and pt | |
407 | for(Int_t icone = 0; icone<fJetNCone; icone++){//icone | |
408 | for(Int_t ipt = 0; ipt<fJetNPt;ipt++){ //ipt | |
409 | ||
410 | //Jet | |
411 | ||
412 | fhJetRatios[icone][ipt] = new TH2F | |
413 | ("JetRatioCone"+fJetNameCones[icone]+"Pt"+fJetNamePtThres[ipt], | |
414 | "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone =" | |
415 | +fJetNameCones[icone]+", pt>" +fJetNamePtThres[ipt]+" GeV/c", | |
416 | 240,0,120,200,0,10); | |
417 | fhJetRatios[icone][ipt]-> | |
418 | SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}"); | |
419 | fhJetRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)"); | |
420 | outputContainer->Add(fhJetRatios[icone][ipt]) ; | |
421 | ||
422 | ||
423 | fhJetPts[icone][ipt] = new TH2F | |
424 | ("JetPtCone"+fJetNameCones[icone]+"Pt"+fJetNamePtThres[ipt], | |
425 | "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone =" | |
426 | +fJetNameCones[icone]+", pt>" +fJetNamePtThres[ipt]+" GeV/c", | |
427 | 240,0,120,400,0,200); | |
428 | fhJetPts[icone][ipt]-> | |
429 | SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}"); | |
430 | fhJetPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)"); | |
431 | outputContainer->Add(fhJetPts[icone][ipt]) ; | |
432 | ||
433 | //Bkg | |
434 | fhBkgRatios[icone][ipt] = new TH2F | |
435 | ("BkgRatioCone"+fJetNameCones[icone]+"Pt"+fJetNamePtThres[ipt], | |
436 | "p_{T bkg lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone =" | |
437 | +fJetNameCones[icone]+", pt>" +fJetNamePtThres[ipt]+" GeV/c", | |
438 | 240,0,120,200,0,10); | |
439 | fhBkgRatios[icone][ipt]-> | |
440 | SetYTitle("p_{T bkg lead #pi^{0}}/p_{T #gamma}"); | |
441 | fhBkgRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)"); | |
442 | outputContainer->Add(fhBkgRatios[icone][ipt]) ; | |
443 | ||
444 | fhBkgPts[icone][ipt] = new TH2F | |
445 | ("BkgPtCone"+fJetNameCones[icone]+"Pt"+fJetNamePtThres[ipt], | |
446 | "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone =" | |
447 | +fJetNameCones[icone]+", pt>" +fJetNamePtThres[ipt]+" GeV/c", | |
448 | 240,0,120,400,0,200); | |
449 | fhBkgPts[icone][ipt]-> | |
450 | SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}"); | |
451 | fhBkgPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)"); | |
452 | outputContainer->Add(fhBkgPts[icone][ipt]) ; | |
453 | ||
454 | //Counts | |
455 | fhNBkgs[icone][ipt] = new TH1F | |
456 | ("NBkgCone"+fJetNameCones[icone]+"Pt"+fJetNamePtThres[ipt], | |
457 | "bkg multiplicity cone ="+fJetNameCones[icone]+", pt>" | |
458 | +fJetNamePtThres[ipt]+" GeV/c",9000,0,9000); | |
459 | fhNBkgs[icone][ipt]->SetYTitle("counts"); | |
460 | fhNBkgs[icone][ipt]->SetXTitle("N"); | |
461 | outputContainer->Add(fhNBkgs[icone][ipt]) ; | |
462 | ||
463 | fhNLeadings[icone][ipt] = new TH2F | |
464 | ("NLeadingCone"+fJetNameCones[icone]+"Pt"+fJetNamePtThres[ipt], | |
465 | "p_{T #gamma} vs p_{T #pi^{0}} cone ="+fJetNameCones[icone]+", pt>" | |
466 | +fJetNamePtThres[ipt]+" GeV/c",120,0,120,120,0,120); | |
467 | fhNLeadings[icone][ipt]->SetYTitle("p_{T #pi^{0}}(GeV/c)"); | |
468 | fhNLeadings[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)"); | |
469 | outputContainer->Add(fhNLeadings[icone][ipt]) ; | |
470 | ||
471 | fhNJets[icone][ipt] = new TH1F | |
472 | ("NJetCone"+fJetNameCones[icone]+"Pt"+fJetNamePtThres[ipt], | |
473 | "Number of neutral jets, cone ="+fJetNameCones[icone]+", pt>" | |
474 | +fJetNamePtThres[ipt]+" GeV/c",120,0,120); | |
475 | fhNJets[icone][ipt]->SetYTitle("N"); | |
476 | fhNJets[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)"); | |
477 | outputContainer->Add(fhNJets[icone][ipt]) ; | |
478 | ||
479 | //Fragmentation Function | |
480 | fhJetFragments[icone][ipt] = new TH2F | |
481 | ("JetFragmentCone"+fJetNameCones[icone]+"Pt"+fJetNamePtThres[ipt], | |
482 | "x_{T} = p_{T i}/p_{T #gamma}, cone ="+fJetNameCones[icone]+", pt>" | |
483 | +fJetNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2); | |
484 | fhJetFragments[icone][ipt]->SetYTitle("x_{T}"); | |
485 | fhJetFragments[icone][ipt]->SetXTitle("p_{T #gamma}"); | |
486 | outputContainer->Add(fhJetFragments[icone][ipt]) ; | |
487 | ||
488 | fhBkgFragments[icone][ipt] = new TH2F | |
489 | ("BkgFragmentCone"+fJetNameCones[icone]+"Pt"+fJetNamePtThres[ipt], | |
490 | "x_{T} = p_{T i}/p_{T #gamma}, cone ="+fJetNameCones[icone]+", pt>" | |
491 | +fJetNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2); | |
492 | fhBkgFragments[icone][ipt]->SetYTitle("x_{T}"); | |
493 | fhBkgFragments[icone][ipt]->SetXTitle("p_{T #gamma}"); | |
494 | outputContainer->Add(fhBkgFragments[icone][ipt]) ; | |
495 | ||
496 | //Jet particle distribution | |
497 | ||
498 | fhJetPtDists[icone][ipt] = new TH2F | |
499 | ("JetPtDistCone"+fJetNameCones[icone]+"Pt"+fJetNamePtThres[ipt], | |
500 | "p_{T i}, cone ="+fJetNameCones[icone]+", pt>" +fJetNamePtThres[ipt]+ | |
501 | " GeV/c",120,0.,120.,120,0.,120.); | |
502 | fhJetPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)"); | |
503 | outputContainer->Add(fhJetPtDists[icone][ipt]) ; | |
504 | ||
505 | fhBkgPtDists[icone][ipt] = new TH2F | |
506 | ("BkgPtDistCone"+fJetNameCones[icone]+"Pt"+fJetNamePtThres[ipt], | |
507 | "p_{T i}, cone ="+fJetNameCones[icone]+", pt>" +fJetNamePtThres[ipt]+ | |
508 | " GeV/c",120,0.,120.,120,0.,120.); | |
509 | fhBkgPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)"); | |
510 | outputContainer->Add(fhBkgPtDists[icone][ipt]) ; | |
511 | ||
512 | }//ipt | |
513 | } //icone | |
514 | }//If we want to study any cone or pt threshold | |
515 | ||
516 | SetOutputContainer(outputContainer); | |
517 | ||
518 | return outputContainer; | |
519 | } | |
520 | ||
521 | //____________________________________________________________________________ | |
522 | void AliAnaGammaJetLeadCone::InitParameters() | |
523 | { | |
524 | ||
525 | //Initialize the parameters of the analysis. | |
526 | SetJetsOnlyInCTS(kFALSE) ; | |
527 | fPbPb = kFALSE ; | |
528 | fEtaEMCALCut = 0.7 ; | |
529 | fPhiEMCALCut[0] = 80. *TMath::Pi()/180.; | |
530 | fPhiEMCALCut[1] = 190.*TMath::Pi()/180.; | |
531 | SetDeltaPhiCutRange(2.9,3.4) ; | |
532 | SetRatioCutRange(0.1,1.0) ; | |
533 | ||
534 | //Jet selection parameters | |
535 | //Fixed cut | |
536 | fJetRatioMaxCut = 1.2 ; | |
537 | fJetRatioMinCut = 0.3 ; | |
538 | fJetCTSRatioMaxCut = 1.2 ; | |
539 | fJetCTSRatioMinCut = 0.3 ; | |
540 | fSelect = 0 ; | |
541 | ||
542 | //Cut depending on gamma energy | |
543 | ||
544 | fPtJetSelectionCut = 20.; //For Low pt jets+BKG, another limits applied | |
545 | //Reconstructed jet energy dependence parameters | |
546 | //e_jet = a1+e_gamma b2. | |
547 | //Index 0-> Pt>2 GeV r = 0.3; Index 1-> Pt>0.5 GeV r = 0.3 | |
548 | fJetE1[0] = -5.75; fJetE1[1] = -4.1; | |
549 | fJetE2[0] = 1.005; fJetE2[1] = 1.05; | |
550 | ||
551 | //Reconstructed sigma of jet energy dependence parameters | |
552 | //s_jet = a1+e_gamma b2. | |
553 | //Index 0-> Pt>2 GeV r = 0.3; Index 1-> Pt>0.5 GeV r = 0.3 | |
554 | fJetSigma1[0] = 2.65; fJetSigma1[1] = 2.75; | |
555 | fJetSigma2[0] = 0.0018; fJetSigma2[1] = 0.033; | |
556 | ||
557 | //Background mean energy and RMS | |
558 | //Index 0-> No BKG; Index 1-> BKG > 2 GeV; | |
559 | //Index 2-> (low pt jets)BKG > 0.5 GeV; | |
560 | //Index > 2, same for CTS conf | |
561 | fBkgMean[0] = 0.; fBkgMean[1] = 8.8 ; fBkgMean[2] = 69.5; | |
562 | fBkgMean[3] = 0.; fBkgMean[4] = 6.4; fBkgMean[5] = 48.6; | |
563 | fBkgRMS[0] = 0.; fBkgRMS[1] = 7.5; fBkgRMS[2] = 22.0; | |
564 | fBkgRMS[3] = 0.; fBkgRMS[4] = 5.4; fBkgRMS[5] = 13.2; | |
565 | ||
566 | //Factor x of min/max = E -+ x * sigma. Obtained after selecting the | |
567 | //limits for monoenergetic jets. | |
568 | //Index 0-> No BKG; Index 1-> BKG > 2 GeV; | |
569 | //Index 2-> (low pt jets) BKG > 0.5 GeV; | |
570 | //Index > 2, same for CTS conf | |
571 | ||
572 | fJetXMin1[0] =-0.69 ; fJetXMin1[1] = 0.39 ; fJetXMin1[2] =-0.88 ; | |
573 | fJetXMin1[3] =-2.0 ; fJetXMin1[4] =-0.442 ; fJetXMin1[5] =-1.1 ; | |
574 | fJetXMin2[0] = 0.066; fJetXMin2[1] = 0.038; fJetXMin2[2] = 0.034; | |
575 | fJetXMin2[3] = 0.25 ; fJetXMin2[4] = 0.113; fJetXMin2[5] = 0.077 ; | |
576 | fJetXMax1[0] =-3.8 ; fJetXMax1[1] =-0.76 ; fJetXMax1[2] =-3.6 ; | |
577 | fJetXMax1[3] =-2.7 ; fJetXMax1[4] =-1.21 ; fJetXMax1[5] =-3.7 ; | |
578 | fJetXMax2[0] =-0.012; fJetXMax2[1] =-0.022; fJetXMax2[2] = 0.016; | |
579 | fJetXMax2[3] =-0.024; fJetXMax2[4] =-0.008; fJetXMax2[5] = 0.027; | |
580 | ||
581 | ||
582 | //Different cones and pt thresholds to construct the jet | |
583 | ||
584 | fJetCone = 0.3 ; | |
585 | fJetPtThreshold = 0.5 ; | |
586 | fJetPtThresPbPb = 2. ; | |
587 | fJetNCone = 4 ; | |
588 | fJetNPt = 4 ; | |
589 | fJetCones[0] = 0.2 ; fJetNameCones[0] = "02" ; | |
590 | fJetCones[1] = 0.3 ; fJetNameCones[1] = "03" ; | |
591 | fJetCones[2] = 0.4 ; fJetNameCones[2] = "04" ; | |
592 | fJetCones[2] = 0.5 ; fJetNameCones[2] = "05" ; | |
593 | ||
594 | fJetPtThres[0] = 0.0 ; fJetNamePtThres[0] = "00" ; | |
595 | fJetPtThres[1] = 0.5 ; fJetNamePtThres[1] = "05" ; | |
596 | fJetPtThres[2] = 1.0 ; fJetNamePtThres[2] = "10" ; | |
597 | fJetPtThres[3] = 2.0 ; fJetNamePtThres[3] = "20" ; | |
598 | } | |
599 | ||
600 | //__________________________________________________________________ | |
601 | void AliAnaGammaJetLeadCone::Print(const Option_t * opt) const | |
602 | { | |
603 | ||
604 | //Print some relevant parameters set for the analysis | |
605 | if(! opt) | |
606 | return; | |
607 | ||
608 | Info("Print", "%s %s", GetName(), GetTitle() ) ; | |
609 | printf("Correlation analysis = %d\n",kJetLeadCone) ; | |
610 | ||
611 | ||
612 | printf("Phi gamma-Leading < %f\n", GetDeltaPhiMaxCut()) ; | |
613 | printf("Phi gamma-Leading > %f\n", GetDeltaPhiMinCut()) ; | |
614 | printf("pT Leading / pT Gamma < %f\n", GetRatioMaxCut()) ; | |
615 | printf("pT Leading / pT Gamma > %f\n", GetRatioMinCut()) ; | |
616 | ||
617 | if(fSelect == 2){ | |
618 | printf("pT Jet / pT Gamma < %f\n", fJetRatioMaxCut) ; | |
619 | printf("pT Jet / pT Gamma > %f\n", fJetRatioMinCut) ; | |
620 | printf("pT Jet (Only CTS)/ pT Gamma < %f\n", fJetCTSRatioMaxCut) ; | |
621 | printf("pT Jet (Only CTS)/ pT Gamma > %f\n", fJetCTSRatioMinCut) ; | |
622 | } | |
623 | ||
624 | } | |
625 | ||
626 | //__________________________________________________________________ | |
627 | void AliAnaGammaJetLeadCone::MakeGammaCorrelation(TParticle * pGamma, TClonesArray * plCTS, TClonesArray * plNe) | |
628 | { | |
629 | //Gamma Jet Correlation Analysis | |
630 | AliDebug(2, "Make correlation"); | |
631 | ||
632 | TParticle *pLeading = new TParticle; //It will contain the kinematics of the found leading particle | |
633 | ||
634 | //Search leading particles in CTS and EMCAL | |
635 | if(GetLeadingParticle(plCTS, plNe, pGamma, pLeading)){ | |
636 | ||
637 | AliDebug(1,Form("Leading: pt %f, phi %f, eta %f", pLeading->Pt(),pLeading->Phi(),pLeading->Eta())) ; | |
638 | ||
639 | //Search Jet | |
640 | if(!fSeveralConeAndPtCuts) | |
641 | MakeJet(plCTS, plNe, pGamma, pLeading, ""); | |
642 | else{ | |
643 | for(Int_t icone = 0; icone<fJetNCone; icone++) { | |
644 | for(Int_t ipt = 0; ipt<fJetNPt;ipt++) { | |
645 | TString lastname ="Cone"+ fJetNameCones[icone]+"Pt"+ fJetNamePtThres[ipt]; | |
646 | fJetCone=fJetCones[icone]; | |
647 | fJetPtThreshold=fJetPtThres[ipt]; | |
648 | MakeJet(plCTS, plNe, pGamma, pLeading, lastname); | |
649 | }//icone | |
650 | }//ipt | |
651 | }//fSeveralConeAndPtCuts | |
652 | }//Leading | |
653 | ||
654 | AliDebug(2, "End of analysis, delete pointers"); | |
655 | ||
656 | pLeading->Delete(); | |
657 | ||
658 | } | |
659 | ||
660 | //____________________________________________________________________________ | |
661 | Bool_t AliAnaGammaJetLeadCone::GetLeadingParticle(TClonesArray * plCTS, TClonesArray * plNe, | |
662 | TParticle * pGamma, TParticle * pLeading) | |
663 | { | |
664 | //Search Charged or Neutral leading particle, select the highest one. | |
665 | ||
666 | TParticle * pLeadingCh = new TParticle(); | |
667 | TParticle * pLeadingPi0 = new TParticle(); | |
668 | ||
669 | Double_t ptg = pGamma->Pt(); | |
670 | Double_t phig = pGamma->Phi(); | |
671 | Double_t etag = pGamma->Eta(); | |
672 | ||
673 | if(!AreJetsOnlyInCTS()) | |
674 | { | |
675 | AliDebug(3, "GetLeadingPi0"); | |
676 | GetLeadingPi0(plNe, pGamma, pLeadingPi0) ; | |
677 | AliDebug(3, "GetLeadingCharge"); | |
678 | GetLeadingCharge(plCTS, pGamma, pLeadingCh) ; | |
679 | ||
680 | Double_t ptch = pLeadingCh->Pt(); | |
681 | Double_t phich = pLeadingCh->Phi(); | |
682 | Double_t etach = pLeadingCh->Eta(); | |
683 | Double_t ptpi = pLeadingPi0->Pt(); | |
684 | Double_t phipi = pLeadingPi0->Phi(); | |
685 | Double_t etapi = pLeadingPi0->Eta(); | |
686 | ||
687 | //Is leading cone inside EMCAL acceptance? | |
688 | ||
689 | Bool_t insidech = kFALSE ; | |
690 | if((phich - fJetCone) > fPhiEMCALCut[0] && | |
691 | (phich + fJetCone) < fPhiEMCALCut[1] && | |
692 | (etach-fJetCone) < fEtaEMCALCut ) | |
693 | insidech = kTRUE ; | |
694 | ||
695 | Bool_t insidepi = kFALSE ; | |
696 | if((phipi - fJetCone) > fPhiEMCALCut[0] && | |
697 | (phipi + fJetCone) < fPhiEMCALCut[1] && | |
698 | (etapi-fJetCone) < fEtaEMCALCut ) | |
699 | insidepi = kTRUE ; | |
700 | ||
701 | AliDebug(2,Form("Leading: charged pt %f, pi0 pt %f",ptch,ptpi)) ; | |
702 | ||
703 | if (ptch > 0 || ptpi > 0) | |
704 | { | |
705 | if(insidech && (ptch > ptpi)) | |
706 | { | |
707 | ||
708 | AliDebug(1,"Leading found in CTS"); | |
709 | pLeading->SetMomentum(pLeadingCh->Px(),pLeadingCh->Py(),pLeadingCh->Pz(),pLeadingCh->Energy()); | |
710 | AliDebug(3,Form("Final leading found in CTS, pt %f, phi %f, eta %f",ptch,phich,etach)) ; | |
711 | fhChargedRatio->Fill(ptg,ptch/pGamma->Pt()); | |
712 | fhEtaCharged->Fill(ptg,etach); | |
713 | fhPhiCharged->Fill(ptg,phich); | |
714 | fhDeltaPhiGammaCharged->Fill(ptg,pGamma->Phi()-phich); | |
715 | fhDeltaEtaGammaCharged->Fill(ptg,pGamma->Eta()-etach); | |
716 | return 1; | |
717 | } | |
718 | ||
719 | else if((ptpi > ptch) && insidepi) | |
720 | { | |
721 | AliDebug(1,"Leading found in EMCAL"); | |
722 | pLeading->SetMomentum(pLeadingPi0->Px(),pLeadingPi0->Py(),pLeadingPi0->Pz(),pLeadingPi0->Energy()); | |
723 | AliDebug(3,Form("Final leading found in EMCAL, pt %f, phi %f, eta %f",ptpi,phipi,etapi)) ; | |
724 | fhNeutralRatio ->Fill(ptg,ptpi/ptg); | |
725 | fhEtaNeutral->Fill(ptg,etapi); | |
726 | fhPhiNeutral->Fill(ptg,phipi); | |
727 | fhDeltaPhiGammaNeutral->Fill(ptg,phig-phipi); | |
728 | fhDeltaEtaGammaNeutral->Fill(ptg,etag-etapi); | |
729 | return 1; | |
730 | } | |
731 | ||
732 | else{ | |
733 | AliDebug(3,"NO LEADING PARTICLE FOUND");} | |
734 | return 0; | |
735 | } | |
736 | else{ | |
737 | AliDebug(3,"NO LEADING PARTICLE FOUND"); | |
738 | return 0; | |
739 | } | |
740 | } | |
741 | else | |
742 | { | |
743 | //No calorimeter present for Leading particle detection | |
744 | GetLeadingCharge(plCTS, pGamma, pLeading) ; | |
745 | Double_t ptch = pLeading->Pt(); | |
746 | Double_t phich = pLeading->Phi(); | |
747 | Double_t etach = pLeading->Eta(); | |
748 | if(ptch > 0){ | |
749 | fhChargedRatio->Fill(ptg,ptch/ptg); | |
750 | fhEtaCharged->Fill(ptg,etach); | |
751 | fhPhiCharged->Fill(ptg,phich); | |
752 | fhDeltaPhiGammaCharged->Fill(ptg,phig-phich); | |
753 | fhDeltaEtaGammaCharged->Fill(ptg,etag-etach); | |
754 | AliDebug(3,Form("Leading found : pt %f, phi %f, eta %f",ptch,phich,etach)) ; | |
755 | return 1; | |
756 | } | |
757 | else | |
758 | { | |
759 | AliDebug(3,"NO LEADING PARTICLE FOUND"); | |
760 | return 0; | |
761 | } | |
762 | } | |
763 | } | |
764 | ||
765 | //____________________________________________________________________________ | |
766 | void AliAnaGammaJetLeadCone::GetLeadingCharge(TClonesArray * pl, TParticle * pGamma, TParticle * pLeading) const | |
767 | { | |
768 | //Search for the charged particle with highest with | |
769 | //Phi=Phi_gamma-Pi and pT=0.1E_gamma | |
770 | Double_t pt = -100.; | |
771 | Double_t phi = -100.; | |
772 | ||
773 | for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){ | |
774 | ||
775 | TParticle * particle = dynamic_cast<TParticle *>(pl->At(ipr)) ; | |
776 | ||
777 | Double_t ptl = particle->Pt(); | |
778 | Double_t rat = ptl/pGamma->Pt() ; | |
779 | Double_t phil = particle->Phi() ; | |
780 | Double_t phig = pGamma->Phi(); | |
781 | ||
782 | //Selection within angular and energy limits | |
783 | if(((phig-phil)> GetDeltaPhiMinCut()) && ((phig-phil)<GetDeltaPhiMaxCut()) && | |
784 | (rat > GetRatioMinCut()) && (rat < GetRatioMaxCut()) && (ptl > pt)) { | |
785 | phi = phil ; | |
786 | pt = ptl ; | |
787 | pLeading->SetMomentum(particle->Px(),particle->Py(),particle->Pz(),particle->Energy()); | |
788 | AliDebug(4,Form("Charge in CTS: pt %f eta %f phi %f pt/Eg %f \n", pt, particle->Eta(), phi,rat)) ; | |
789 | } | |
790 | } | |
791 | ||
792 | AliDebug(1,Form("Leading in CTS: pt %f eta %f phi %f pt/Eg %f \n", pt, pLeading->Eta(), phi,pt/pGamma->Pt())) ; | |
793 | ||
794 | } | |
795 | ||
796 | ||
797 | //____________________________________________________________________________ | |
798 | void AliAnaGammaJetLeadCone::GetLeadingPi0(TClonesArray * pl, TParticle * pGamma, TParticle * pLeading) | |
799 | { | |
800 | ||
801 | //Search for the neutral pion with highest with | |
802 | //Phi=Phi_gamma-Pi and pT=0.1E_gamma | |
803 | Double_t pt = -100.; | |
804 | Double_t phi = -100.; | |
805 | Double_t ptl = -100.; | |
806 | Double_t rat = -100.; | |
807 | Double_t phil = -100. ; | |
808 | Double_t ptg = pGamma->Pt(); | |
809 | Double_t phig = pGamma->Phi(); | |
810 | ||
811 | TIter next(pl); | |
812 | TParticle * particlei = 0; | |
813 | TParticle * particlej = 0; | |
814 | TLorentzVector gammai; | |
815 | TLorentzVector gammaj; | |
816 | ||
817 | Int_t iPrimary = -1; | |
818 | Double_t angle = 0., e = 0., invmass = 0.; | |
819 | Int_t ksPdg = 0; | |
820 | Int_t jPrimary=-1; | |
821 | ||
822 | while ( (particlei = (TParticle*)next()) ) { | |
823 | iPrimary++; | |
824 | ksPdg = particlei->GetPdgCode(); | |
825 | ||
826 | ptl = particlei->Pt(); | |
827 | //2 gamma overlapped, found with PID | |
828 | if(ksPdg == 111){ | |
829 | rat = ptl/ptg ; | |
830 | phil = particlei->Phi() ; | |
831 | //Selection within angular and energy limits | |
832 | if(ptl > pt && rat > GetRatioMinCut() && rat < GetRatioMaxCut() && | |
833 | (phig-phil) > GetDeltaPhiMinCut() && (phig-phil) < GetDeltaPhiMaxCut() ) | |
834 | { | |
835 | phi = phil ; | |
836 | pt = ptl ; | |
837 | pLeading->SetMomentum(particlei->Px(),particlei->Py(),particlei->Pz(),particlei->Energy()); | |
838 | AliDebug(4,Form("Pi0 candidate: pt %f eta %f phi %f pt/Eg %f \n", pLeading->Pt(), pLeading->Eta(), pLeading->Phi(), pLeading->Pt()/pGamma->Pt())) ; | |
839 | }// cuts | |
840 | }// pdg = 111 | |
841 | ||
842 | //Make invariant mass analysis | |
843 | else if(ksPdg == 22){// gamma i | |
844 | ||
845 | //Search the photon companion in case it comes from a Pi0 decay | |
846 | //Apply several cuts to select the good pair | |
847 | particlei->Momentum(gammai); | |
848 | jPrimary=-1; | |
849 | TIter next2(pl); | |
850 | while ( (particlej = (TParticle*)next2()) ) { | |
851 | jPrimary++; | |
852 | if(jPrimary>iPrimary){ | |
853 | ksPdg = particlej->GetPdgCode(); | |
854 | particlej->Momentum(gammaj); | |
855 | if(ksPdg == 22 && GetNeutralMesonSelection()->SelectPair(pGamma, gammai, gammaj)) | |
856 | { | |
857 | //Info("GetLeadingPi0","Egammai %f, Egammaj %f", | |
858 | //gammai.Pt(),gammaj.Pt()); | |
859 | ptl = (gammai+gammaj).Pt(); | |
860 | phil = (gammai+gammaj).Phi(); | |
861 | if(phil < 0) | |
862 | phil+=TMath::TwoPi(); | |
863 | rat = ptl/ptg ; | |
864 | ||
865 | if(ptl > pt ){ | |
866 | pt = ptl; | |
867 | phi = phil ; | |
868 | e = (gammai+gammaj).E(); | |
869 | angle = gammaj.Angle(gammai.Vect()); | |
870 | invmass = (gammai+gammaj).M(); | |
871 | pLeading->SetMomentum( (gammai+gammaj).Px(), (gammai+gammaj).Py(), (gammai+gammaj).Pz(), (gammai+gammaj).E()); | |
872 | AliDebug(4,Form("Pi0 candidate: pt %f eta %f phi %f pt/Eg %f \n", pLeading->Pt(), pLeading->Eta(), pLeading->Phi(), pLeading->Pt()/pGamma->Pt())) ; | |
873 | } | |
874 | }//if pdg = 22, pair selected | |
875 | }//iprimary<jprimary | |
876 | }//while | |
877 | }// if pdg = 22 | |
878 | }//while | |
879 | ||
880 | //Final pi0 found, highest pair energy, fill histograms | |
881 | if(e > 0.){ | |
882 | fhInvMassPairLeading->Fill(e,invmass); | |
883 | fhAnglePairLeading->Fill(e,angle); | |
884 | } | |
885 | ||
886 | AliDebug(3,Form("Leading EMCAL: pt %f eta %f phi %f pt/Eg %f \n", pLeading->Pt(), pLeading->Eta(), pLeading->Phi(), pLeading->Pt()/pGamma->Pt())) ; | |
887 | } | |
888 | ||
889 | ||
890 | //__________________________________________________________________________- | |
891 | Bool_t AliAnaGammaJetLeadCone::IsJetSelected(const Double_t ptg, const Double_t ptj){ | |
892 | //Check if the energy of the reconstructed jet is within an energy window | |
893 | ||
894 | Double_t par[6]; | |
895 | Double_t xmax[2]; | |
896 | Double_t xmin[2]; | |
897 | ||
898 | Int_t iCTS = 0; | |
899 | if(AreJetsOnlyInCTS()) | |
900 | iCTS = 3 ; | |
901 | ||
902 | if(!fPbPb){ | |
903 | //Phythia alone, jets with pt_th > 0.2, r = 0.3 | |
904 | par[0] = fJetE1[0]; par[1] = fJetE2[0]; | |
905 | //Energy of the jet peak | |
906 | //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit | |
907 | par[2] = fJetSigma1[0]; par[3] = fJetSigma2[0]; | |
908 | //Sigma of the jet peak | |
909 | //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit | |
910 | par[4] = fBkgMean[0 + iCTS]; par[5] = fBkgRMS[0 + iCTS]; | |
911 | //Parameters reserved for PbPb bkg. | |
912 | xmax[0] = fJetXMax1[0 + iCTS]; xmax[1] = fJetXMax2[0 + iCTS]; | |
913 | xmin[0] = fJetXMin1[0 + iCTS]; xmin[1] = fJetXMin2[0 + iCTS]; | |
914 | //Factor that multiplies sigma to obtain the best limits, | |
915 | //by observation, of mono jet ratios (ptjet/ptg) | |
916 | //X_jet = fJetX1[0]+fJetX2[0]*e_gamma | |
917 | ||
918 | } | |
919 | else{ | |
920 | if(ptg > fPtJetSelectionCut){ | |
921 | //Phythia +PbPb with pt_th > 2 GeV/c, r = 0.3 | |
922 | par[0] = fJetE1[0]; par[1] = fJetE2[0]; | |
923 | //Energy of the jet peak, same as in pp | |
924 | //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit | |
925 | par[2] = fJetSigma1[0]; par[3] = fJetSigma2[0]; | |
926 | //Sigma of the jet peak, same as in pp | |
927 | //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit | |
928 | par[4] = fBkgMean[1 + iCTS]; par[5] = fBkgRMS[1 + iCTS]; | |
929 | //Mean value and RMS of PbPb Bkg | |
930 | xmax[0] = fJetXMax1[1 + iCTS]; xmax[1] = fJetXMax2[1 + iCTS]; | |
931 | xmin[0] = fJetXMin1[1 + iCTS]; xmin[1] = fJetXMin2[1 + iCTS]; | |
932 | //Factor that multiplies sigma to obtain the best limits, | |
933 | //by observation, of mono jet ratios (ptjet/ptg) mixed with PbPb Bkg, | |
934 | //pt_th > 2 GeV, r = 0.3 | |
935 | //X_jet = fJetX1[0]+fJetX2[0]*e_gamma | |
936 | ||
937 | } | |
938 | else{ | |
939 | //Phythia + PbPb with pt_th > 0.5 GeV/c, r = 0.3 | |
940 | par[0] = fJetE1[1]; par[1] = fJetE2[1]; | |
941 | //Energy of the jet peak, pt_th > 2 GeV/c, r = 0.3 | |
942 | //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit | |
943 | par[2] = fJetSigma1[1]; par[3] = fJetSigma2[1]; | |
944 | //Sigma of the jet peak, pt_th > 2 GeV/c, r = 0.3 | |
945 | //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit | |
946 | par[4] = fBkgMean[2 + iCTS]; par[5] = fBkgRMS[2 + iCTS]; | |
947 | //Mean value and RMS of PbPb Bkg in a 0.3 cone, pt > 2 GeV. | |
948 | xmax[0] = fJetXMax1[2 + iCTS]; xmax[1] = fJetXMax2[2 + iCTS]; | |
949 | xmin[0] = fJetXMin1[2 + iCTS]; xmin[1] = fJetXMin2[2 + iCTS]; | |
950 | //Factor that multiplies sigma to obtain the best limits, | |
951 | //by observation, of mono jet ratios (ptjet/ptg) mixed with PbPb Bkg, | |
952 | //pt_th > 2 GeV, r = 0.3 | |
953 | //X_jet = fJetX1[0]+fJetX2[0]*e_gamma | |
954 | ||
955 | }//If low pt jet in bkg | |
956 | }//if Bkg | |
957 | ||
958 | //Calculate minimum and maximum limits of the jet ratio. | |
959 | Double_t min = CalculateJetRatioLimit(ptg, par, xmin); | |
960 | Double_t max = CalculateJetRatioLimit(ptg, par, xmax); | |
961 | ||
962 | AliDebug(3,Form("Jet selection? : Limits min %f, max %f, pt_jet %f, pt_gamma %f, pt_jet / pt_gamma %f",min,max,ptj,ptg,ptj/ptg)); | |
963 | ||
964 | if(( min < ptj/ptg ) && ( max > ptj/ptg)) | |
965 | return kTRUE; | |
966 | else | |
967 | return kFALSE; | |
968 | ||
969 | } | |
970 | ||
971 | //____________________________________________________________________________ | |
972 | void AliAnaGammaJetLeadCone::MakeJet(TClonesArray * plCTS, TClonesArray * plNe, | |
973 | TParticle * pGamma, TParticle* pLeading,TString lastname) | |
974 | { | |
975 | //Fill the jet with the particles around the leading particle with | |
976 | //R=fJetCone and pt_th = fJetPtThres. Calculate the energy of the jet and | |
977 | //check if we select it. Fill jet histograms | |
978 | ||
979 | TClonesArray * jetList = new TClonesArray("TParticle",1000); | |
980 | TClonesArray * bkgList = new TClonesArray("TParticle",1000); | |
981 | ||
982 | TLorentzVector jet (0,0,0,0); | |
983 | TLorentzVector bkg(0,0,0,0); | |
984 | TLorentzVector lv (0,0,0,0); | |
985 | ||
986 | Double_t ptjet = 0.0; | |
987 | Double_t ptbkg = 0.0; | |
988 | Int_t n0 = 0; | |
989 | Int_t n1 = 0; | |
990 | Bool_t b1 = kFALSE; | |
991 | Bool_t b0 = kFALSE; | |
992 | ||
993 | Double_t ptg = pGamma->Pt(); | |
994 | Double_t phig = pGamma->Phi(); | |
995 | Double_t ptl = pLeading->Pt(); | |
996 | Double_t phil = pLeading->Phi(); | |
997 | Double_t etal = pLeading->Eta(); | |
998 | ||
999 | Float_t ptcut = fJetPtThreshold; | |
1000 | if(fPbPb && !fSeveralConeAndPtCuts && ptg > fPtJetSelectionCut) ptcut = fJetPtThresPbPb ; | |
1001 | ||
1002 | //Add charged particles to jet | |
1003 | TIter nextch(plCTS) ; | |
1004 | TParticle * particle = 0 ; | |
1005 | while ( (particle = dynamic_cast<TParticle*>(nextch())) ) { | |
1006 | ||
1007 | b0 = kFALSE; | |
1008 | b1 = kFALSE; | |
1009 | ||
1010 | //Particles in jet | |
1011 | SetJet(particle, b0, fJetCone, etal, phil) ; | |
1012 | ||
1013 | if(b0){ | |
1014 | new((*jetList)[n0++]) TParticle(*particle) ; | |
1015 | particle->Momentum(lv); | |
1016 | if(particle->Pt() > ptcut ){ | |
1017 | jet+=lv; | |
1018 | ptjet+=particle->Pt(); | |
1019 | } | |
1020 | } | |
1021 | ||
1022 | //Background around (phi_gamma-pi, eta_leading) | |
1023 | SetJet(particle, b1, fJetCone,etal, phig) ; | |
1024 | ||
1025 | if(b1) { | |
1026 | new((*bkgList)[n1++]) TParticle(*particle) ; | |
1027 | particle->Momentum(lv); | |
1028 | if(particle->Pt() > ptcut ){ | |
1029 | bkg+=lv; | |
1030 | ptbkg+=particle->Pt(); | |
1031 | } | |
1032 | } | |
1033 | } | |
1034 | ||
1035 | //Add neutral particles to jet | |
1036 | TIter nextne(plNe) ; | |
1037 | particle = 0 ; | |
1038 | while ( (particle = dynamic_cast<TParticle*>(nextne())) ) { | |
1039 | ||
1040 | b0 = kFALSE; | |
1041 | b1 = kFALSE; | |
1042 | ||
1043 | //Particles in jet | |
1044 | SetJet(particle, b0, fJetCone, etal, phil) ; | |
1045 | ||
1046 | if(b0){ | |
1047 | new((*jetList)[n0++]) TParticle(*particle) ; | |
1048 | particle->Momentum(lv); | |
1049 | if(particle->Pt() > ptcut ){ | |
1050 | jet+=lv; | |
1051 | ptjet+=particle->Pt(); | |
1052 | } | |
1053 | } | |
1054 | ||
1055 | //Background around (phi_gamma-pi, eta_leading) | |
1056 | SetJet(particle, b1, fJetCone,etal, phig) ; | |
1057 | ||
1058 | if(b1) { | |
1059 | new((*bkgList)[n1++]) TParticle(*particle) ; | |
1060 | particle->Momentum(lv); | |
1061 | if(particle->Pt() > ptcut ){ | |
1062 | bkg+=lv; | |
1063 | ptbkg+=particle->Pt(); | |
1064 | } | |
1065 | } | |
1066 | } | |
1067 | ||
1068 | ptjet = jet.Pt(); | |
1069 | ptbkg = bkg.Pt(); | |
1070 | ||
1071 | if(ptjet > 0.) { | |
1072 | ||
1073 | AliDebug(2,Form("Gamma pt %f, Jet pt %f, Bkg pt %f",ptg,ptjet,ptbkg)); | |
1074 | ||
1075 | //Fill histograms | |
1076 | ||
1077 | Double_t ratjet = ptjet/ptg ; | |
1078 | Double_t ratbkg = ptbkg/ptg ; | |
1079 | ||
1080 | dynamic_cast<TH2F*> | |
1081 | (GetOutputContainer()->FindObject("JetRatio"+lastname)) | |
1082 | ->Fill(ptg,ratjet); | |
1083 | dynamic_cast<TH2F*> | |
1084 | (GetOutputContainer()->FindObject("JetPt"+lastname)) | |
1085 | ->Fill(ptg,ptjet); | |
1086 | ||
1087 | dynamic_cast<TH2F*> | |
1088 | (GetOutputContainer()->FindObject("BkgRatio"+lastname)) | |
1089 | ->Fill(ptg,ratbkg); | |
1090 | ||
1091 | dynamic_cast<TH2F*> | |
1092 | (GetOutputContainer()->FindObject("BkgPt"+lastname)) | |
1093 | ->Fill(ptg,ptbkg); | |
1094 | ||
1095 | //Jet selection | |
1096 | Bool_t kSelect = kFALSE; | |
1097 | if(fSelect == 0) | |
1098 | kSelect = kTRUE; //Accept all jets, no restriction | |
1099 | else if(fSelect == 1){ | |
1100 | //Selection with parametrized cuts | |
1101 | if(IsJetSelected(ptg,ptjet)) kSelect = kTRUE; | |
1102 | } | |
1103 | else if(fSelect == 2){ | |
1104 | //Simple selection | |
1105 | if(!AreJetsOnlyInCTS()){ | |
1106 | if((ratjet < fJetRatioMaxCut) && (ratjet > fJetRatioMinCut )) kSelect = kTRUE; | |
1107 | } | |
1108 | else{ | |
1109 | if((ratjet < fJetCTSRatioMaxCut) && (ratjet > fJetCTSRatioMinCut )) kSelect = kTRUE; | |
1110 | } | |
1111 | } | |
1112 | else | |
1113 | AliError("Jet selection option larger than 2, DON'T SELECT JETS"); | |
1114 | ||
1115 | ||
1116 | if(kSelect){ | |
1117 | AliDebug(1,Form("Jet Selected: pt %f ", ptjet)) ; | |
1118 | ||
1119 | FillJetHistos(jetList, ptg, ptl,"Jet",lastname); | |
1120 | FillJetHistos(bkgList, ptg, ptl, "Bkg",lastname); | |
1121 | } | |
1122 | } //ptjet > 0 | |
1123 | ||
1124 | jetList ->Delete(); | |
1125 | bkgList ->Delete(); | |
1126 | ||
1127 | } | |
1128 | ||
1129 | //___________________________________________________________________ | |
1130 | void AliAnaGammaJetLeadCone::SetJet(TParticle * part, Bool_t & b, Float_t cone, | |
1131 | Double_t eta, Double_t phi) | |
1132 | { | |
1133 | ||
1134 | //Check if the particle is inside the cone defined by the leading particle | |
1135 | b = kFALSE; | |
1136 | ||
1137 | if(phi > TMath::TwoPi()) | |
1138 | phi-=TMath::TwoPi(); | |
1139 | if(phi < 0.) | |
1140 | phi+=TMath::TwoPi(); | |
1141 | ||
1142 | Double_t rad = 10000 + cone; | |
1143 | ||
1144 | if(TMath::Abs(part->Phi()-phi) <= (TMath::TwoPi() - cone)) | |
1145 | rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+ | |
1146 | TMath::Power(part->Phi()-phi,2)); | |
1147 | else{ | |
1148 | if(part->Phi()-phi > TMath::TwoPi() - cone) | |
1149 | rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+ | |
1150 | TMath::Power((part->Phi()-TMath::TwoPi())-phi,2)); | |
1151 | if(part->Phi()-phi < -(TMath::TwoPi() - cone)) | |
1152 | rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+ | |
1153 | TMath::Power((part->Phi()+TMath::TwoPi())-phi,2)); | |
1154 | } | |
1155 | ||
1156 | if(rad < cone ) | |
1157 | b = kTRUE; | |
1158 | ||
1159 | } | |
1160 | ||
1161 | //____________________________________________________________________________ | |
1162 | void AliAnaGammaJetLeadCone::FillJetHistos(TClonesArray * pl, Double_t ptg, Double_t ptl, TString type, TString lastname) | |
1163 | { | |
1164 | //Fill histograms wth jet fragmentation | |
1165 | //and number of selected jets and leading particles | |
1166 | //and the background multiplicity | |
1167 | TParticle * particle = 0 ; | |
1168 | Int_t ipr = 0; | |
1169 | Float_t charge = 0; | |
1170 | ||
1171 | TIter next(pl) ; | |
1172 | while ( (particle = dynamic_cast<TParticle*>(next())) ) { | |
1173 | ipr++ ; | |
1174 | Double_t pt = particle->Pt(); | |
1175 | ||
1176 | charge = TDatabasePDG::Instance() | |
1177 | ->GetParticle(particle->GetPdgCode())->Charge(); | |
1178 | if(charge != 0){//Only jet Charged particles | |
1179 | dynamic_cast<TH2F*> | |
1180 | (GetOutputContainer()->FindObject(type+"Fragment"+lastname)) | |
1181 | ->Fill(ptg,pt/ptg); | |
1182 | dynamic_cast<TH2F*> | |
1183 | (GetOutputContainer()->FindObject(type+"PtDist"+lastname)) | |
1184 | ->Fill(ptg,pt); | |
1185 | }//charged | |
1186 | ||
1187 | }//while | |
1188 | ||
1189 | if(type == "Bkg") | |
1190 | dynamic_cast<TH1F*> | |
1191 | (GetOutputContainer()->FindObject("NBkg"+lastname)) | |
1192 | ->Fill(ipr); | |
1193 | else{ | |
1194 | dynamic_cast<TH1F*> | |
1195 | (GetOutputContainer()->FindObject("NJet"+lastname))-> | |
1196 | Fill(ptg); | |
1197 | dynamic_cast<TH2F*> | |
1198 | (GetOutputContainer()->FindObject("NLeading"+lastname)) | |
1199 | ->Fill(ptg,ptl); | |
1200 | } | |
1201 | ||
1202 | } |