1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 /* History of cvs commits:
20 * Revision 1.1.2.1 2007/07/26 10:32:09 schutz
21 * new analysis classes in the the new analysis framework
26 //_________________________________________________________________________
27 // Class for the analysis of gamma-jet correlations, jet reconstructed in cone around leading
29 // Class created from old AliPHOSGammaJet
30 // (see AliRoot versions previous Release 4-09)
32 //*-- Author: Gustavo Conesa (LNF-INFN)
33 //////////////////////////////////////////////////////////////////////////////
36 // --- ROOT system ---
37 #include "Riostream.h"
39 //---- AliRoot system ----
41 #include "AliNeutralMesonSelection.h"
42 #include "AliAnaGammaJetLeadCone.h"
44 ClassImp(AliAnaGammaJetLeadCone)
47 //____________________________________________________________________________
48 AliAnaGammaJetLeadCone::AliAnaGammaJetLeadCone() :
49 AliAnaGammaCorrelation(), fPbPb(kFALSE),
50 fSeveralConeAndPtCuts(0),
52 fJetCTSRatioMaxCut(0.), fJetCTSRatioMinCut(0.), fJetRatioMaxCut(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)
68 //Initialize parameters
70 SetCorrelationType(kJetLeadCone);
72 for(Int_t i = 0; i<10; i++){
74 fJetNameCones[i] = "" ;
75 fJetPtThres[i] = 0.0 ;
76 fJetNamePtThres[i] = "" ;
89 fPhiEMCALCut[i] = 0.0 ;
97 //____________________________________________________________________________
98 AliAnaGammaJetLeadCone::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)
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] ;
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] ;
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] ;
147 //_________________________________________________________________________
148 AliAnaGammaJetLeadCone & AliAnaGammaJetLeadCone::operator = (const AliAnaGammaJetLeadCone & source)
150 // assignment operator
152 if(this == &source)return *this;
153 ((AliAnaGammaCorrelation *)this)->operator=(source);
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 ;
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 ;
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 ;
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] ;
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] ;
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] ;
209 //____________________________________________________________________________
210 AliAnaGammaJetLeadCone::~AliAnaGammaJetLeadCone()
213 delete fhChargedRatio ;
214 delete fhNeutralRatio ;
216 delete fhPhiCharged ;
217 delete fhPhiNeutral ;
218 delete fhEtaCharged ;
219 delete fhEtaNeutral ;
220 delete fhDeltaPhiGammaCharged ;
221 delete fhDeltaPhiGammaNeutral ;
222 delete fhDeltaEtaGammaCharged ;
223 delete fhDeltaEtaGammaNeutral ;
225 delete fhAnglePairLeading ;
226 delete fhInvMassPairLeading ;
234 delete fhJetFragment ;
235 delete fhBkgFragment ;
239 delete [] fhJetRatios;
241 delete [] fhBkgRatios;
244 delete [] fhNLeadings;
248 delete [] fhJetFragments;
249 delete [] fhBkgFragments;
250 delete [] fhJetPtDists;
251 delete [] fhBkgPtDists;
256 //________________________________________________________________________
257 TList * AliAnaGammaJetLeadCone::GetCreateOutputObjects()
259 // Create histograms to be saved in output file and
260 // store them in fOutputContainer
262 AliDebug(1,"Init jet in leading cone histograms");
264 TList * outputContainer = new TList() ;
265 outputContainer->SetName("GammaJetCorrelationHistos") ;
267 fhChargedRatio = new TH2F
268 ("ChargedRatio","p_{T leading charge} /p_{T #gamma} vs p_{T #gamma}",
270 fhChargedRatio->SetYTitle("p_{T lead charge} /p_{T #gamma}");
271 fhChargedRatio->SetXTitle("p_{T #gamma} (GeV/c)");
273 fhPhiCharged = new TH2F
274 ("PhiCharged","#phi_{#pi^{#pm}} vs p_{T #gamma}",
276 fhPhiCharged->SetYTitle("#phi_{#pi^{#pm}} (rad)");
277 fhPhiCharged->SetXTitle("p_{T #gamma} (GeV/c)");
279 fhEtaCharged = new TH2F
280 ("EtaCharged","#eta_{#pi^{#pm}} vs p_{T #gamma}",
282 fhEtaCharged->SetYTitle("#eta_{#pi^{#pm}} (rad)");
283 fhEtaCharged->SetXTitle("p_{T #gamma} (GeV/c)");
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)");
291 fhDeltaEtaGammaCharged = new TH2F
292 ("DeltaEtaGammaCharged","#eta_{#gamma} - #eta_{#pi^{#pm}} vs p_{T #gamma}",
294 fhDeltaEtaGammaCharged->SetYTitle("#Delta #eta");
295 fhDeltaEtaGammaCharged->SetXTitle("p_{T #gamma} (GeV/c)");
297 outputContainer->Add(fhPhiCharged) ;
298 outputContainer->Add(fhEtaCharged) ;
299 outputContainer->Add(fhChargedRatio) ;
300 outputContainer->Add(fhDeltaPhiGammaCharged) ;
301 outputContainer->Add(fhDeltaEtaGammaCharged) ;
303 if(!AreJetsOnlyInCTS()){
305 fhNeutralRatio = new TH2F
306 ("NeutralRatio","p_{T leading #pi^{0}} /p_{T #gamma} vs p_{T #gamma}",
308 fhNeutralRatio->SetYTitle("p_{T lead #pi^{0}} /p_{T #gamma}");
309 fhNeutralRatio->SetXTitle("p_{T #gamma} (GeV/c)");
312 fhAnglePairLeading = new TH2F
314 "Angle between all #gamma pair finally selected vs p_{T #pi^{0}}",
316 fhAnglePairLeading->SetYTitle("Angle (rad)");
317 fhAnglePairLeading->SetXTitle("E_{ #pi^{0}} (GeV/c)");
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)");
326 fhPhiNeutral = new TH2F
327 ("PhiNeutral","#phi_{#pi^{0}} vs p_{T #gamma}",
329 fhPhiNeutral->SetYTitle("#phi_{#pi^{0}} (rad)");
330 fhPhiNeutral->SetXTitle("p_{T #gamma} (GeV/c)");
332 fhEtaNeutral = new TH2F
333 ("EtaNeutral","#eta_{#pi^{0}} vs p_{T #gamma}",
335 fhEtaNeutral->SetYTitle("#eta_{#pi^{0}} (rad)");
336 fhEtaNeutral->SetXTitle("p_{T #gamma} (GeV/c)");
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)");
344 fhDeltaEtaGammaNeutral = new TH2F
345 ("DeltaEtaGammaNeutral","#eta_{#gamma} - #eta_{#pi^{#pm}} vs p_{T #gamma}",
347 fhDeltaEtaGammaNeutral->SetYTitle("#Delta #eta");
348 fhDeltaEtaGammaNeutral->SetXTitle("p_{T #gamma} (GeV/c)");
350 outputContainer->Add(fhPhiNeutral) ;
351 outputContainer->Add(fhEtaNeutral) ;
352 outputContainer->Add(fhNeutralRatio) ;
353 outputContainer->Add(fhDeltaPhiGammaNeutral) ;
354 outputContainer->Add(fhDeltaEtaGammaNeutral) ;
356 outputContainer->Add(fhInvMassPairLeading) ;
357 outputContainer->Add(fhAnglePairLeading) ;
360 if(!fSeveralConeAndPtCuts){// not several cones
363 fhNBkg = new TH1F("NBkg","bkg multiplicity",9000,0,9000);
364 fhNBkg->SetYTitle("counts");
365 fhNBkg->SetXTitle("N");
366 outputContainer->Add(fhNBkg) ;
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) ;
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) ;
379 //Ratios and Pt dist of reconstructed (not selected) jets
381 fhJetRatio = new TH2F
382 ("JetRatio","p_{T jet lead}/p_{T #gamma} vs p_{T #gamma}",
384 fhJetRatio->SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
385 fhJetRatio->SetXTitle("p_{T #gamma} (GeV/c)");
386 outputContainer->Add(fhJetRatio) ;
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) ;
396 fhBkgRatio = new TH2F
397 ("BkgRatio","p_{T bkg lead}/p_{T #gamma} vs p_{T #gamma}",
399 fhBkgRatio->SetYTitle("p_{T bkg lead charge}/p_{T #gamma}");
400 fhBkgRatio->SetXTitle("p_{T #gamma} (GeV/c)");
401 outputContainer->Add(fhBkgRatio) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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
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",
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]) ;
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]) ;
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",
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]) ;
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]) ;
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]) ;
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]) ;
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]) ;
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]) ;
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]) ;
526 //Jet particle distribution
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]) ;
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]) ;
544 }//If we want to study any cone or pt threshold
546 SetOutputContainer(outputContainer);
548 return outputContainer;
551 //____________________________________________________________________________
552 void AliAnaGammaJetLeadCone::InitParameters()
555 //Initialize the parameters of the analysis.
556 SetJetsOnlyInCTS(kFALSE) ;
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) ;
564 //Jet selection parameters
566 fJetRatioMaxCut = 1.2 ;
567 fJetRatioMinCut = 0.3 ;
568 fJetCTSRatioMaxCut = 1.2 ;
569 fJetCTSRatioMinCut = 0.3 ;
572 //Cut depending on gamma energy
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;
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;
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;
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
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;
612 //Different cones and pt thresholds to construct the jet
615 fJetPtThreshold = 0.5 ;
616 fJetPtThresPbPb = 2. ;
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" ;
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" ;
630 //__________________________________________________________________
631 void AliAnaGammaJetLeadCone::Print(const Option_t * opt) const
634 //Print some relevant parameters set for the analysis
638 Info("Print", "%s %s", GetName(), GetTitle() ) ;
639 printf("Correlation analysis = %d\n",kJetLeadCone) ;
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()) ;
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) ;
656 //__________________________________________________________________
657 void AliAnaGammaJetLeadCone::MakeGammaCorrelation(TParticle * pGamma, TClonesArray * plCTS, TClonesArray * plNe)
659 //Gamma Jet Correlation Analysis
660 AliDebug(2, "Make correlation");
662 TParticle *pLeading = new TParticle; //It will contain the kinematics of the found leading particle
664 //Search leading particles in CTS and EMCAL
665 if(GetLeadingParticle(plCTS, plNe, pGamma, pLeading)){
667 AliDebug(1,Form("Leading: pt %f, phi %f, eta %f", pLeading->Pt(),pLeading->Phi(),pLeading->Eta())) ;
670 if(!fSeveralConeAndPtCuts)
671 MakeJet(plCTS, plNe, pGamma, pLeading, "");
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);
681 }//fSeveralConeAndPtCuts
684 AliDebug(2, "End of analysis, delete pointers");
690 //____________________________________________________________________________
691 Bool_t AliAnaGammaJetLeadCone::GetLeadingParticle(TClonesArray * plCTS, TClonesArray * plNe,
692 TParticle * pGamma, TParticle * pLeading)
694 //Search Charged or Neutral leading particle, select the highest one.
696 TParticle * pLeadingCh = new TParticle();
697 TParticle * pLeadingPi0 = new TParticle();
699 Double_t ptg = pGamma->Pt();
700 Double_t phig = pGamma->Phi();
701 Double_t etag = pGamma->Eta();
703 if(!AreJetsOnlyInCTS())
705 AliDebug(3, "GetLeadingPi0");
706 GetLeadingPi0(plNe, pGamma, pLeadingPi0) ;
707 AliDebug(3, "GetLeadingCharge");
708 GetLeadingCharge(plCTS, pGamma, pLeadingCh) ;
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();
717 //Is leading cone inside EMCAL acceptance?
719 Bool_t insidech = kFALSE ;
720 if((phich - fJetCone) > fPhiEMCALCut[0] &&
721 (phich + fJetCone) < fPhiEMCALCut[1] &&
722 (etach-fJetCone) < fEtaEMCALCut )
725 Bool_t insidepi = kFALSE ;
726 if((phipi - fJetCone) > fPhiEMCALCut[0] &&
727 (phipi + fJetCone) < fPhiEMCALCut[1] &&
728 (etapi-fJetCone) < fEtaEMCALCut )
731 AliDebug(2,Form("Leading: charged pt %f, pi0 pt %f",ptch,ptpi)) ;
733 if (ptch > 0 || ptpi > 0)
735 if(insidech && (ptch > ptpi))
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);
749 else if((ptpi > ptch) && insidepi)
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);
763 AliDebug(3,"NO LEADING PARTICLE FOUND");}
767 AliDebug(3,"NO LEADING PARTICLE FOUND");
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();
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)) ;
789 AliDebug(3,"NO LEADING PARTICLE FOUND");
795 //____________________________________________________________________________
796 void AliAnaGammaJetLeadCone::GetLeadingCharge(TClonesArray * pl, TParticle * pGamma, TParticle * pLeading) const
798 //Search for the charged particle with highest with
799 //Phi=Phi_gamma-Pi and pT=0.1E_gamma
801 Double_t phi = -100.;
803 for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
805 TParticle * particle = dynamic_cast<TParticle *>(pl->At(ipr)) ;
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();
812 //Selection within angular and energy limits
813 if(((phig-phil)> GetDeltaPhiMinCut()) && ((phig-phil)<GetDeltaPhiMaxCut()) &&
814 (rat > GetRatioMinCut()) && (rat < GetRatioMaxCut()) && (ptl > pt)) {
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)) ;
822 AliDebug(1,Form("Leading in CTS: pt %f eta %f phi %f pt/Eg %f \n", pt, pLeading->Eta(), phi,pt/pGamma->Pt())) ;
827 //____________________________________________________________________________
828 void AliAnaGammaJetLeadCone::GetLeadingPi0(TClonesArray * pl, TParticle * pGamma, TParticle * pLeading)
831 //Search for the neutral pion with highest with
832 //Phi=Phi_gamma-Pi and pT=0.1E_gamma
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();
842 TParticle * particlei = 0;
843 TParticle * particlej = 0;
844 TLorentzVector gammai;
845 TLorentzVector gammaj;
848 Double_t angle = 0., e = 0., invmass = 0.;
852 while ( (particlei = (TParticle*)next()) ) {
854 ksPdg = particlei->GetPdgCode();
856 ptl = particlei->Pt();
857 //2 gamma overlapped, found with PID
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() )
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())) ;
872 //Make invariant mass analysis
873 else if(ksPdg == 22){// gamma i
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);
880 while ( (particlej = (TParticle*)next2()) ) {
882 if(jPrimary>iPrimary){
883 ksPdg = particlej->GetPdgCode();
884 particlej->Momentum(gammaj);
885 if(ksPdg == 22 && GetNeutralMesonSelection()->SelectPair(pGamma, gammai, gammaj))
887 //Info("GetLeadingPi0","Egammai %f, Egammaj %f",
888 //gammai.Pt(),gammaj.Pt());
889 ptl = (gammai+gammaj).Pt();
890 phil = (gammai+gammaj).Phi();
892 phil+=TMath::TwoPi();
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())) ;
904 }//if pdg = 22, pair selected
910 //Final pi0 found, highest pair energy, fill histograms
912 fhInvMassPairLeading->Fill(e,invmass);
913 fhAnglePairLeading->Fill(e,angle);
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())) ;
920 //__________________________________________________________________________-
921 Bool_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
929 if(AreJetsOnlyInCTS())
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
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
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
985 }//If low pt jet in bkg
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);
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));
994 if(( min < ptj/ptg ) && ( max > ptj/ptg))
1001 //____________________________________________________________________________
1002 void AliAnaGammaJetLeadCone::MakeJet(TClonesArray * plCTS, TClonesArray * plNe,
1003 TParticle * pGamma, TParticle* pLeading,TString lastname)
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
1009 TClonesArray * jetList = new TClonesArray("TParticle",1000);
1010 TClonesArray * bkgList = new TClonesArray("TParticle",1000);
1012 TLorentzVector jet (0,0,0,0);
1013 TLorentzVector bkg(0,0,0,0);
1014 TLorentzVector lv (0,0,0,0);
1016 Double_t ptjet = 0.0;
1017 Double_t ptbkg = 0.0;
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();
1029 Float_t ptcut = fJetPtThreshold;
1030 if(fPbPb && !fSeveralConeAndPtCuts && ptg > fPtJetSelectionCut) ptcut = fJetPtThresPbPb ;
1032 //Add charged particles to jet
1033 TIter nextch(plCTS) ;
1034 TParticle * particle = 0 ;
1035 while ( (particle = dynamic_cast<TParticle*>(nextch())) ) {
1041 SetJet(particle, b0, fJetCone, etal, phil) ;
1044 new((*jetList)[n0++]) TParticle(*particle) ;
1045 particle->Momentum(lv);
1046 if(particle->Pt() > ptcut ){
1048 ptjet+=particle->Pt();
1052 //Background around (phi_gamma-pi, eta_leading)
1053 SetJet(particle, b1, fJetCone,etal, phig) ;
1056 new((*bkgList)[n1++]) TParticle(*particle) ;
1057 particle->Momentum(lv);
1058 if(particle->Pt() > ptcut ){
1060 ptbkg+=particle->Pt();
1065 //Add neutral particles to jet
1066 TIter nextne(plNe) ;
1068 while ( (particle = dynamic_cast<TParticle*>(nextne())) ) {
1074 SetJet(particle, b0, fJetCone, etal, phil) ;
1077 new((*jetList)[n0++]) TParticle(*particle) ;
1078 particle->Momentum(lv);
1079 if(particle->Pt() > ptcut ){
1081 ptjet+=particle->Pt();
1085 //Background around (phi_gamma-pi, eta_leading)
1086 SetJet(particle, b1, fJetCone,etal, phig) ;
1089 new((*bkgList)[n1++]) TParticle(*particle) ;
1090 particle->Momentum(lv);
1091 if(particle->Pt() > ptcut ){
1093 ptbkg+=particle->Pt();
1103 AliDebug(2,Form("Gamma pt %f, Jet pt %f, Bkg pt %f",ptg,ptjet,ptbkg));
1107 Double_t ratjet = ptjet/ptg ;
1108 Double_t ratbkg = ptbkg/ptg ;
1111 (GetOutputContainer()->FindObject("JetRatio"+lastname))
1114 (GetOutputContainer()->FindObject("JetPt"+lastname))
1118 (GetOutputContainer()->FindObject("BkgRatio"+lastname))
1122 (GetOutputContainer()->FindObject("BkgPt"+lastname))
1126 Bool_t kSelect = kFALSE;
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;
1133 else if(fSelect == 2){
1135 if(!AreJetsOnlyInCTS()){
1136 if((ratjet < fJetRatioMaxCut) && (ratjet > fJetRatioMinCut )) kSelect = kTRUE;
1139 if((ratjet < fJetCTSRatioMaxCut) && (ratjet > fJetCTSRatioMinCut )) kSelect = kTRUE;
1143 AliError("Jet selection option larger than 2, DON'T SELECT JETS");
1147 AliDebug(1,Form("Jet Selected: pt %f ", ptjet)) ;
1149 FillJetHistos(jetList, ptg, ptl,"Jet",lastname);
1150 FillJetHistos(bkgList, ptg, ptl, "Bkg",lastname);
1159 //___________________________________________________________________
1160 void AliAnaGammaJetLeadCone::SetJet(TParticle * part, Bool_t & b, Float_t cone,
1161 Double_t eta, Double_t phi)
1164 //Check if the particle is inside the cone defined by the leading particle
1167 if(phi > TMath::TwoPi())
1168 phi-=TMath::TwoPi();
1170 phi+=TMath::TwoPi();
1172 Double_t rad = 10000 + cone;
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));
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));
1191 //____________________________________________________________________________
1192 void AliAnaGammaJetLeadCone::FillJetHistos(TClonesArray * pl, Double_t ptg, Double_t ptl, TString type, TString lastname)
1194 //Fill histograms wth jet fragmentation
1195 //and number of selected jets and leading particles
1196 //and the background multiplicity
1197 TParticle * particle = 0 ;
1202 while ( (particle = dynamic_cast<TParticle*>(next())) ) {
1204 Double_t pt = particle->Pt();
1206 charge = TDatabasePDG::Instance()
1207 ->GetParticle(particle->GetPdgCode())->Charge();
1208 if(charge != 0){//Only jet Charged particles
1210 (GetOutputContainer()->FindObject(type+"Fragment"+lastname))
1213 (GetOutputContainer()->FindObject(type+"PtDist"+lastname))
1221 (GetOutputContainer()->FindObject("NBkg"+lastname))
1225 (GetOutputContainer()->FindObject("NJet"+lastname))->
1228 (GetOutputContainer()->FindObject("NLeading"+lastname))