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