]>
Commit | Line | Data |
---|---|---|
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 | ||
61 | ClassImp(AliAnaGammaJet) | |
62 | ||
63 | //____________________________________________________________________________ | |
64 | AliAnaGammaJet::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 | //____________________________________________________________________________ | |
127 | AliAnaGammaJet::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 | //_________________________________________________________________________ |
183 | AliAnaGammaJet & 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 | //____________________________________________________________________________ |
240 | AliAnaGammaJet::~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 | //____________________________________________________________________________ | |
293 | Double_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 | //______________________________________________________________________________ |
311 | void 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 | //________________________________________________________________________ |
319 | void 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 | //____________________________________________________________________________ |
667 | void 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 | //____________________________________________________________________________ |
764 | void 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 | //____________________________________________________________________________ |
807 | void 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 | 839 | void 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 | //____________________________________________________________________________ | |
950 | Bool_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 | //____________________________________________________________________________ | |
1051 | void 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 | //__________________________________________________________________________- | |
1145 | Bool_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 | //__________________________________________________________________________- | |
1165 | Bool_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 | //____________________________________________________________________________ | |
1246 | void 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 | //____________________________________________________________________________ | |
1374 | void 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 | //___________________________________________________________________ | |
1402 | void 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 | ||
1434 | void 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 | } |