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