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