]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/AliAnaGammaJetLeadCone.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Changes in AliESDCaloCluster propagated here
[u/mrichter/AliRoot.git] / PWG4 / AliAnaGammaJetLeadCone.cxx
CommitLineData
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
52ClassImp(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//____________________________________________________________________________
106AliAnaGammaJetLeadCone::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//_________________________________________________________________________
156AliAnaGammaJetLeadCone & 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//____________________________________________________________________________
218AliAnaGammaJetLeadCone::~AliAnaGammaJetLeadCone()
219{
3788def4 220 // Remove all pointers except analysis output pointers.
bdcfac30 221
bdcfac30 222}
223
224
225
226//________________________________________________________________________
227TList * 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//__________________________________________________________________
601void 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//__________________________________________________________________
627void 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//____________________________________________________________________________
661Bool_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//____________________________________________________________________________
766void 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//____________________________________________________________________________
798void 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//__________________________________________________________________________-
891Bool_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//____________________________________________________________________________
972void 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//___________________________________________________________________
1130void 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//____________________________________________________________________________
1162void 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}