AliPhysics  6cf2591 (6cf2591)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnaParticleJetFinderCorrelation.cxx
Go to the documentation of this file.
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 
16 // --- ROOT system ---
17 #include "TH2F.h"
18 #include "TClonesArray.h"
19 #include "TClass.h"
20 //#include "Riostream.h"
21 
22 // ---- AliRoot system ----
23 #include "AliCaloTrackReader.h"
24 #include "AliAODJet.h"
26 #include "AliAODPWG4ParticleCorrelation.h"
27 #include "AliVTrack.h"
28 #include "AliAODCaloCluster.h"
29 #include "AliAODEvent.h"
30 
31 // ---- Jets ----
32 #include "AliAODJetEventBackground.h"
33 #include "TRandom2.h"
34 //MC access
35 #include "AliStack.h"
36 #include "AliMCAnalysisUtils.h"
37 #include "AliAODMCParticle.h"
38 
39 // --- Detectors ---
40 #include "AliEMCALGeometry.h"
41 
45 
46 //________________________________________________________________________
48 //________________________________________________________________________
51  fDeltaPhiMaxCut(0.), fDeltaPhiMinCut(0.), fRatioMaxCut(0.), fRatioMinCut(0.),
52  fConeSize(0.), fPtThresholdInCone(0.),fUseJetRefTracks(kTRUE),
53  fMakeCorrelationInHistoMaker(kFALSE), fSelectIsolated(kTRUE),
54  fJetConeSize(0.4),fJetMinPt(5),fJetMinPtBkgSub(-100.),fJetAreaFraction(0.6),
55 //fNonStandardJetFromReader(kTRUE),
56  fJetBranchName("jets"),
57  fBackgroundJetFromReader(kTRUE),
58  fBkgJetBranchName("jets"),
59  fGammaConeSize(0.3),fUseBackgroundSubtractionGamma(kFALSE),fSaveGJTree(kTRUE),
60  fMostEnergetic(kFALSE),fMostOpposite(kTRUE), fUseHistogramJetBkg(kTRUE),
61  fUseHistogramTracks(kTRUE),fUseHistogramJetTracks(kTRUE),fMCStudies(kFALSE),fGenerator(0),
62  fMomentum(),
63  fhDeltaEta(0), /*fhDeltaPhi(0),*/fhDeltaPhiCorrect(0),fhDeltaPhi0PiCorrect(0), fhDeltaPt(0), fhPtRatio(0), fhPt(0),
64  fhFFz(0),fhFFxi(0),fhFFpt(0),fhNTracksInCone(0),
65  fhJetFFz(0),fhJetFFxi(0),fhJetFFpt(0),fhJetFFzCor(0),fhJetFFxiCor(0),
66  fhGamPtPerTrig(0),fhPtGamPtJet(0),
67  fhBkgFFz(),fhBkgFFxi(),fhBkgFFpt(),fhBkgNTracksInCone(),fhBkgSumPtInCone(),fhBkgSumPtnTracksInCone(),
68  fhNjetsNgammas(0),fhCuts(0),
69  fhDeltaEtaBefore(0),fhDeltaPhiBefore(0),fhDeltaPtBefore(0),fhPtRatioBefore(0),
70  fhPtBefore(0),fhDeltaPhi0PiCorrectBefore(0),
71  fhJetPtBefore(0),fhJetPtBeforeCut(0),fhJetPt(0),fhJetPtMostEne(0),fhJetPhi(0),fhJetEta(0),fhJetEtaVsPt(0),
72  fhJetPhiVsEta(0),fhJetEtaVsNpartInJet(0),fhJetEtaVsNpartInJetBkg(0),fhJetChBkgEnergyVsPt(0),fhJetChAreaVsPt(0),/*fhJetNjet(0),*/
73  fhTrackPhiVsEta(0),fhTrackAveTrackPt(0),fhJetNjetOverPtCut(),
74 /*fhJetChBkgEnergyVsPtEtaGt05(0),fhJetChBkgEnergyVsPtEtaLe05(0),fhJetChAreaVsPtEtaGt05(0),fhJetChAreaVsPtEtaLe05(0),*/
75  fhJetChBkgEnergyVsArea(0),fhJetRhoVsPt(0),fhJetRhoVsCentrality(0),//fhJetBkgRho(0),
76  fhJetNparticlesInJet(0),fhJetDeltaEtaDeltaPhi(0),fhJetDeltaEtaDeltaPhiAllTracks(0),
77  fhJetAveTrackPt(0),fhJetNtracksInJetAboveThr(),fhJetRatioNTrkAboveToNTrk(),fhJetNtrackRatioMostEne(),
78  fhJetNtrackRatioJet5GeV(),fhJetNtrackRatioLead5GeV(),
79  fhBkgJetBackground(),fhBkgJetSigma(),fhBkgJetArea(),fhPhotonPtMostEne(0),
80  fhPhotonAverageEnergy(0),fhPhotonRatioAveEneToMostEne(0),fhPhotonAverageEnergyMinus1(0),fhPhotonRatioAveEneMinus1ToMostEne(0),
81  fhPhotonNgammaMoreAverageToNgamma(0),fhPhotonNgammaMoreAverageMinus1ToNgamma(0),fhPhotonNgammaOverPtCut(),
82  fhPhotonBkgRhoVsNtracks(0),fhPhotonBkgRhoVsNclusters(0),fhPhotonBkgRhoVsCentrality(0),
83  fhPhotonBkgRhoVsNcells(0),fhPhotonPt(0),fhPhotonPtCorrected(0),fhPhotonPtCorrectedZoom(0),fhPhotonPtDiff(0),
84  fhPhotonPtDiffVsCentrality(0),fhPhotonPtDiffVsNcells(0),fhPhotonPtDiffVsNtracks(0),fhPhotonPtDiffVsNclusters(0),
85  fhPhotonSumPtInCone(0),fhPhotonSumPtCorrectInCone(0),fhPhotonSumPtChargedInCone(0),
86  fhSelectedJetPhiVsEta(0),fhSelectedJetChBkgEnergyVsPtJet(0),fhSelectedJetChAreaVsPtJet(0),fhSelectedJetNjet(0),fhSelectedNtracks(0),
87  fhSelectedTrackPhiVsEta(0),fhCuts2(0),
88  fhSelectedPhotonNLMVsPt(0),fhSelectedPhotonLambda0VsPt(0), fhRandomPhiEta(),
89  fhMCPhotonCuts(0),fhMCPhotonPt(0),fhMCPhotonEtaPhi(0),fhMCJetOrigin(0),
90  fhMCJetNPartVsPt(0),fhMCJetChNPartVsPt(0),fhMCJetNPart150VsPt(0),fhMCJetChNPart150VsPt(0),fhMCJetChNPart150ConeVsPt(0),
91  fhMCJetRatioChFull(0),fhMCJetRatioCh150Ch(0),
92  fhMCJetEtaPhi(0),fhMCJetChEtaPhi(0),fhMCJet150EtaPhi(0),fhMCJetCh150EtaPhi(0),fhMCJetCh150ConeEtaPhi(0),
93 fTreeGJ (0),
94 fGamPt (0),
95 fGamLambda0 (0),
96 fGamNLM (0),
97 fGamSumPtCh (0),
98 fGamTime (0),
99 fGamNcells (0),
100 fGamEta (0),
101 fGamPhi (0),
102 fGamSumPtNeu(0),
103 fGamNtracks (0),
104 fGamNclusters(0),
105 fGamAvEne (0),
106 fJetPhi (0),
107 fJetEta (0),
108 fJetPt (0),
109 fJetBkgChEne(0),
110 fJetArea (0),
111 fJetNtracks (0),
112 fJetNtracks1(0),
113 fJetNtracks2(0),
114 fJetRho(0),
115 fEventNumber(0),
116 fNtracks (0),
117 fZvertex (0),
118 fCentrality (0),
119 fIso(0),
120 fGamRho(0),
121 fMCGamPt (0),
122 fMCGamEta (0),
123 fMCGamPhi (0),
124 fMCPartonType (0),
125 fMCJetPt (0),
126 fMCJetChPt (0),
127 fMCJet150Pt (0),
128 fMCJetCh150Pt (0),
129 fMCJetNPart (0),
130 fMCJetChNPart (0),
131 fMCJet150NPart (0),
132 fMCJetCh150NPart(0),
133 fMCJetEta (0),
134 fMCJetPhi (0),
135 fMCJetChEta (0),
136 fMCJetChPhi (0),
137 fMCJet150Eta (0),
138 fMCJet150Phi (0),
139 fMCJetCh150Eta (0),
140 fMCJetCh150Phi (0),
141 fMCJetCh150ConePt(0),
142 fMCJetCh150ConeNPart(0),
143 fMCJetCh150ConeEta(0),
144 fMCJetCh150ConePhi(0)
145 {
146  InitParameters();
147  for(Int_t i=0;i<10;i++)
148  {
149  fhJetNjetOverPtCut [i] = 0;
151  }
152 
153  fGenerator = new TRandom2();
154  fGenerator->SetSeed(0);
155 }
156 
157 //___________________________________________________________________
159 //___________________________________________________________________
161 {
162  delete fGenerator;
163 }
164 
165 //___________________________________________________________________
168 //___________________________________________________________________
170 {
171  TList * outputContainer = new TList() ;
172  outputContainer->SetName("ParticleJetFinderHistos") ;
173 
175  // Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins();
176  // Int_t netabins = GetHistogramRanges()->GetHistoEtaBins();
178  // Float_t phimax = GetHistogramRanges()->GetHistoPhiMax();
179  // Float_t etamax = GetHistogramRanges()->GetHistoEtaMax();
181  // Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
182 // Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
183 
184 // fhDeltaPhi = new TH2F("DeltaPhi","#phi_{trigger} - #phi_{jet} vs p_{T trigger}",nptbins,ptmin,ptmax,100,-6,4);
185 // fhDeltaPhi->SetYTitle("#Delta #phi");
186 // fhDeltaPhi->SetXTitle("p_{T trigger} (GeV/c)");
187 // outputContainer->Add(fhDeltaPhi);
188 
189  fhDeltaPhiCorrect = new TH2F("DeltaPhiCorrect","#phi_{trigger} - #phi_{jet} vs p_{T trigger}",nptbins,ptmin,ptmax,100,0,6.5);
190  fhDeltaPhiCorrect->SetYTitle("#Delta #phi");
191  fhDeltaPhiCorrect->SetXTitle("p_{T trigger} (GeV/c)");
192  outputContainer->Add(fhDeltaPhiCorrect);
193 
194  fhDeltaPhi0PiCorrect = new TH2F("DeltaPhi0PiCorrect","#phi_{trigger} - #phi_{jet} (0,#pi) vs p_{T trigger}",nptbins,ptmin,ptmax,100,0,3.5);
195  fhDeltaPhi0PiCorrect->SetYTitle("#Delta #phi");
196  fhDeltaPhi0PiCorrect->SetXTitle("p_{T trigger} (GeV/c)");
197  outputContainer->Add(fhDeltaPhi0PiCorrect);
198 
199 
200  fhDeltaEta = new TH2F("DeltaEta","#eta_{trigger} - #eta_{jet} vs p_{T trigger}",nptbins,ptmin,ptmax,100,-2,2);
201  fhDeltaEta->SetYTitle("#Delta #eta");
202  fhDeltaEta->SetXTitle("p_{T trigger} (GeV/c)");
203  outputContainer->Add(fhDeltaEta);
204 
205  fhDeltaPt = new TH2F("DeltaPt","p_{T trigger} - p_{T jet} vs p_{T trigger}",nptbins,ptmin,ptmax,150,-50,100);
206  fhDeltaPt->SetYTitle("#Delta p_{T}");
207  fhDeltaPt->SetXTitle("p_{T trigger} (GeV/c)");
208  outputContainer->Add(fhDeltaPt);
209 
210  fhPtRatio = new TH2F("PtRatio","p_{T jet} / p_{T trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,200,0,2.);
211  fhPtRatio->SetYTitle("ratio");
212  fhPtRatio->SetXTitle("p_{T trigger} (GeV/c)");
213  outputContainer->Add(fhPtRatio);
214 
215  fhPt = new TH2F("Pt","p_{T jet} vs p_{T trigger}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
216  fhPt->SetYTitle("p_{T jet}(GeV/c)");
217  fhPt->SetXTitle("p_{T trigger} (GeV/c)");
218  outputContainer->Add(fhPt);
219 
220  fhFFz = new TH2F("FFz","z = p_{T i charged}/p_{T trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,200,0.,2);
221  fhFFz->SetYTitle("z");
222  fhFFz->SetXTitle("p_{T trigger}");
223  outputContainer->Add(fhFFz) ;
224 
225  fhFFxi = new TH2F("FFxi","#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger}", nptbins,ptmin,ptmax,100,0.,10.);
226  fhFFxi->SetYTitle("#xi");
227  fhFFxi->SetXTitle("p_{T trigger}");
228  outputContainer->Add(fhFFxi) ;
229 
230  fhFFpt = new TH2F("FFpt","p_{T i charged} vs p_{T trigger}", nptbins,ptmin,ptmax,100,0.,50.);
231  fhFFpt->SetYTitle("p_{T charged hadron}");
232  fhFFpt->SetXTitle("p_{T trigger}");
233  outputContainer->Add(fhFFpt) ;
234 
235  fhNTracksInCone = new TH2F("NTracksInCone","Number of tracks in cone vs p_{T trigger}", nptbins,ptmin,ptmax,100,0.,150.);
236  fhNTracksInCone->SetYTitle("p_{T charged hadron}");
237  fhNTracksInCone->SetXTitle("p_{T trigger}");
238  outputContainer->Add(fhNTracksInCone) ;
239 
240  //FF according to jet axis
241  fhJetFFz = new TH2F("JetFFz","z = p_{T i charged}/p_{T jet} vs p_{T jet}",nptbins,ptmin,ptmax,200,0.,2);
242  fhJetFFz->SetYTitle("z");
243  fhJetFFz->SetXTitle("p_{T jet}");
244  outputContainer->Add(fhJetFFz) ;
245 
246  fhJetFFxi = new TH2F("JetFFxi","#xi = ln(p_{T jet}/p_{T i charged}) vs p_{T jet}", nptbins,ptmin,ptmax,100,0.,10.);
247  fhJetFFxi->SetYTitle("#xi");
248  fhJetFFxi->SetXTitle("p_{T jet}");
249  outputContainer->Add(fhJetFFxi) ;
250 
251  fhJetFFpt = new TH2F("JetFFpt","p_{T i charged} vs p_{T jet}", nptbins,ptmin,ptmax,100,0.,50.);
252  fhJetFFpt->SetYTitle("p_{T charged hadron}");
253  fhJetFFpt->SetXTitle("p_{T jet}");
254  outputContainer->Add(fhJetFFpt) ;
255 
256  fhJetFFzCor = new TH2F("JetFFzCor","z = -cos(#alpha(jet,trig))*p_{T i charged}/p_{T jet} vs p_{T jet}",nptbins,ptmin,ptmax,200,0.,2);
257  fhJetFFzCor->SetYTitle("z");
258  fhJetFFzCor->SetXTitle("p_{T jet}");
259  outputContainer->Add(fhJetFFzCor) ;
260 
261  fhJetFFxiCor = new TH2F("JetFFxiCor","#xi = ln(p_{T jet}/(-cos(#alpha(jet,trig))*p_{T i charged})) vs p_{T jet}", nptbins,ptmin,ptmax,100,0.,10.);
262  fhJetFFxiCor->SetYTitle("#xi");
263  fhJetFFxiCor->SetXTitle("p_{T jet}");
264  outputContainer->Add(fhJetFFxiCor) ;
265 
266  fhGamPtPerTrig = new TH1F("GamPtPerTrig","GamPtPerTrig", nptbins,ptmin,ptmax);
267  fhGamPtPerTrig->SetYTitle("Counts");
268  fhGamPtPerTrig->SetXTitle("p_{T, #gamma}");
269  outputContainer->Add(fhGamPtPerTrig) ;
270 
271  fhPtGamPtJet = new TH2F("PtGamPtJet","p_{T #gamma} vs p_{T jet}", nptbins,ptmin,ptmax,150,-50.,100.);
272  fhPtGamPtJet->SetXTitle("p_{T #gamma}");
273  fhPtGamPtJet->SetYTitle("p_{T jet}");
274  outputContainer->Add(fhPtGamPtJet) ;
275 
276 
277  //background FF
278  fhBkgFFz[0] = new TH2F("BkgFFzRC", "z = p_{T i charged}/p_{T trigger} vs p_{T trigger} Bkg RC" ,nptbins,ptmin,ptmax,200,0.,2);
279  fhBkgFFz[1] = new TH2F("BkgFFzPCG", "z = p_{T i charged}/p_{T trigger} vs p_{T trigger} Bkg PCG" ,nptbins,ptmin,ptmax,200,0.,2);
280  fhBkgFFz[2] = new TH2F("BkgFFzPCJ", "z = p_{T i charged}/p_{T trigger} vs p_{T trigger} Bkg PCJ" ,nptbins,ptmin,ptmax,200,0.,2);
281  fhBkgFFz[3] = new TH2F("BkgFFzMP", "z = p_{T i charged}/p_{T trigger} vs p_{T trigger} Bkg MP" ,nptbins,ptmin,ptmax,200,0.,2);
282  fhBkgFFz[4] = new TH2F("BkgFFzTest","z = p_{T i charged}/p_{T trigger} vs p_{T trigger} Bkg Test",nptbins,ptmin,ptmax,200,0.,2);
283  for(Int_t i=0;i<5;i++){
284  fhBkgFFz[i]->SetYTitle("z");
285  fhBkgFFz[i]->SetXTitle("p_{T trigger}");
286  outputContainer->Add(fhBkgFFz[i]) ;
287  }
288 
289  fhBkgFFxi[0] = new TH2F("BkgFFxiRC", "#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger} Bkg RC", nptbins,ptmin,ptmax,100,0.,10.);
290  fhBkgFFxi[1] = new TH2F("BkgFFxiPCG", "#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger} Bkg PCG", nptbins,ptmin,ptmax,100,0.,10.);
291  fhBkgFFxi[2] = new TH2F("BkgFFxiPCJ", "#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger} Bkg PCJ", nptbins,ptmin,ptmax,100,0.,10.);
292  fhBkgFFxi[3] = new TH2F("BkgFFxiMP", "#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger} Bkg MP", nptbins,ptmin,ptmax,100,0.,10.);
293  fhBkgFFxi[4] = new TH2F("BkgFFxiTest","#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger} Bkg Test",nptbins,ptmin,ptmax,100,0.,10.);
294  for(Int_t i=0;i<5;i++){
295  fhBkgFFxi[i]->SetYTitle("#xi");
296  fhBkgFFxi[i]->SetXTitle("p_{T trigger}");
297  outputContainer->Add(fhBkgFFxi[i]) ;
298  }
299 
300  fhBkgFFpt[0] = new TH2F("BkgFFptRC", "p_{T i charged} vs p_{T trigger} Bkg RC", nptbins,ptmin,ptmax,100,0.,50.);
301  fhBkgFFpt[1] = new TH2F("BkgFFptPCG", "p_{T i charged} vs p_{T trigger} Bkg PCG", nptbins,ptmin,ptmax,100,0.,50.);
302  fhBkgFFpt[2] = new TH2F("BkgFFptPCJ", "p_{T i charged} vs p_{T trigger} Bkg PCJ", nptbins,ptmin,ptmax,100,0.,50.);
303  fhBkgFFpt[3] = new TH2F("BkgFFptMP", "p_{T i charged} vs p_{T trigger} Bkg MP", nptbins,ptmin,ptmax,100,0.,50.);
304  fhBkgFFpt[4] = new TH2F("BkgFFptTest","p_{T i charged} vs p_{T trigger} Bkg Test", nptbins,ptmin,ptmax,100,0.,50.);
305  for(Int_t i=0;i<5;i++){
306  fhBkgFFpt[i]->SetYTitle("p_{T charged hadron}");
307  fhBkgFFpt[i]->SetXTitle("p_{T trigger}");
308  outputContainer->Add(fhBkgFFpt[i]) ;
309  }
310 
311  fhBkgNTracksInCone[0] = new TH2F("BkgNTracksInConeRC", "Number of tracks in cone vs p_{T trigger} Bkg RC", nptbins,ptmin,ptmax,100,0.,150.);
312  fhBkgNTracksInCone[1] = new TH2F("BkgNTracksInConePCG", "Number of tracks in cone vs p_{T trigger} Bkg PCG", nptbins,ptmin,ptmax,100,0.,150.);
313  fhBkgNTracksInCone[2] = new TH2F("BkgNTracksInConePCJ", "Number of tracks in cone vs p_{T trigger} Bkg PCJ", nptbins,ptmin,ptmax,100,0.,150.);
314  fhBkgNTracksInCone[3] = new TH2F("BkgNTracksInConeMP", "Number of tracks in cone vs p_{T trigger} Bkg MP", nptbins,ptmin,ptmax,100,0.,150.);
315  fhBkgNTracksInCone[4] = new TH2F("BkgNTracksInConeTest","Number of tracks in cone vs p_{T trigger} Bkg Test", nptbins,ptmin,ptmax,100,0.,150.);
316  for(Int_t i=0;i<5;i++){
317  fhBkgNTracksInCone[i]->SetYTitle("Number of tracks");
318  fhBkgNTracksInCone[i]->SetXTitle("p_{T trigger}");
319  outputContainer->Add(fhBkgNTracksInCone[i]) ;
320  }
321 
322  fhBkgSumPtInCone[0] = new TH2F("BkgSumPtInConeRC", "Sum P_{T} in cone vs p_{T trigger} Bkg RC", nptbins,ptmin,ptmax,100,0.,100.);
323  fhBkgSumPtInCone[1] = new TH2F("BkgSumPtInConePCG", "Sum P_{T} in cone vs p_{T trigger} Bkg PCG", nptbins,ptmin,ptmax,100,0.,100.);
324  fhBkgSumPtInCone[2] = new TH2F("BkgSumPtInConePCJ", "Sum P_{T} in cone vs p_{T trigger} Bkg PCJ", nptbins,ptmin,ptmax,100,0.,100.);
325  fhBkgSumPtInCone[3] = new TH2F("BkgSumPtInConeMP", "Sum P_{T} in cone vs p_{T trigger} Bkg MP", nptbins,ptmin,ptmax,100,0.,100.);
326  fhBkgSumPtInCone[4] = new TH2F("BkgSumPtInConeTest","Sum P_{T} in cone vs p_{T trigger} Bkg Test", nptbins,ptmin,ptmax,100,0.,100.);
327  for(Int_t i=0;i<5;i++){
328  fhBkgSumPtInCone[i]->SetYTitle("Sum P_{T}");
329  fhBkgSumPtInCone[i]->SetXTitle("p_{T trigger}");
330  outputContainer->Add(fhBkgSumPtInCone[i]) ;
331  }
332 
333  fhBkgSumPtnTracksInCone[0] = new TH2F("BkgSumPtnTracksInConeRC", "Sum p_{T} / Number of tracks in cone vs p_{T trigger} Bkg RC", nptbins,ptmin,ptmax,100,0.,20.);
334  fhBkgSumPtnTracksInCone[1] = new TH2F("BkgSumPtnTracksInConePCG", "Sum p_{T} / Number of tracks in cone vs p_{T trigger} Bkg PCG", nptbins,ptmin,ptmax,100,0.,20.);
335  fhBkgSumPtnTracksInCone[2] = new TH2F("BkgSumPtnTracksInConePCJ", "Sum p_{T} / Number of tracks in cone vs p_{T trigger} Bkg PCJ", nptbins,ptmin,ptmax,100,0.,20.);
336  fhBkgSumPtnTracksInCone[3] = new TH2F("BkgSumPtnTracksInConeMP", "Sum p_{T} / Number of tracks in cone vs p_{T trigger} Bkg MP", nptbins,ptmin,ptmax,100,0.,20.);
337  fhBkgSumPtnTracksInCone[4] = new TH2F("BkgSumPtnTracksInConeTest","Sum p_{T} / Number of tracks in cone vs p_{T trigger} Bkg Test", nptbins,ptmin,ptmax,100,0.,20.);
338  for(Int_t i=0;i<5;i++){
339  fhBkgSumPtnTracksInCone[i]->SetYTitle("Sum p_{T}/Number of tracks");
340  fhBkgSumPtnTracksInCone[i]->SetXTitle("p_{T trigger}");
341  outputContainer->Add(fhBkgSumPtnTracksInCone[i]) ;
342  }
343 
344 
345  //temporary histograms
346  fhNjetsNgammas = new TH2F("NjetsNgammas"," Number of jets vs number of gammas in event",20,0.,100.,10,0.,80.);
347  fhNjetsNgammas->SetYTitle("Number of gammas");
348  fhNjetsNgammas->SetXTitle("Number of jets");
349  outputContainer->Add(fhNjetsNgammas) ;
350 
351  fhCuts = new TH1F("Cuts"," Cuts",10,0.,10.);
352  fhCuts->SetYTitle("Counts");
353  fhCuts->SetXTitle("Cut number");
354  outputContainer->Add(fhCuts) ;
355 
356  fhDeltaPhiBefore = new TH2F("DeltaPhiBefore","#phi_{trigger} - #phi_{jet} vs p_{T trigger}",nptbins,ptmin,ptmax,100,0,6.5);
357  fhDeltaPhiBefore->SetYTitle("#Delta #phi");
358  fhDeltaPhiBefore->SetXTitle("p_{T trigger} (GeV/c)");
359  outputContainer->Add(fhDeltaPhiBefore);
360 
361  fhDeltaEtaBefore = new TH2F("DeltaEtaBefore","#eta_{trigger} - #eta_{jet} vs p_{T trigger}",nptbins,ptmin,ptmax,100,-2,2);
362  fhDeltaEtaBefore->SetYTitle("#Delta #eta");
363  fhDeltaEtaBefore->SetXTitle("p_{T trigger} (GeV/c)");
364  outputContainer->Add(fhDeltaEtaBefore);
365 
366  fhDeltaPtBefore = new TH2F("DeltaPtBefore","p_{T trigger} - p_{T jet} vs p_{T trigger}",nptbins,ptmin,ptmax,100,-50,50);
367  fhDeltaPtBefore->SetYTitle("#Delta p_{T}");
368  fhDeltaPtBefore->SetXTitle("p_{T trigger} (GeV/c)");
369  outputContainer->Add(fhDeltaPtBefore);
370 
371  fhPtRatioBefore = new TH2F("PtRatioBefore","p_{T jet} / p_{T trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,200,0,2.);
372  fhPtRatioBefore->SetYTitle("ratio");
373  fhPtRatioBefore->SetXTitle("p_{T trigger} (GeV/c)");
374  outputContainer->Add(fhPtRatioBefore);
375 
376  fhPtBefore = new TH2F("PtBefore","p_{T jet} vs p_{T trigger}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
377  fhPtBefore->SetYTitle("p_{T jet}(GeV/c)");
378  fhPtBefore->SetXTitle("p_{T trigger} (GeV/c)");
379  outputContainer->Add(fhPtBefore);
380 
381  fhDeltaPhi0PiCorrectBefore = new TH2F("DeltaPhi0PiCorrectBefore","#phi_{trigger} - #phi_{jet} (0,#pi) vs p_{T trigger}",nptbins,ptmin,ptmax,100,0,3.5);
382  fhDeltaPhi0PiCorrectBefore->SetYTitle("#Delta #phi");
383  fhDeltaPhi0PiCorrectBefore->SetXTitle("p_{T trigger} (GeV/c)");
384  outputContainer->Add(fhDeltaPhi0PiCorrectBefore);
385 
386  //temporary jet histograms
387  fhJetPtBefore = new TH1F("JetPtBefore","JetPtBefore",150,-50,100);
388  fhJetPtBefore->SetYTitle("Counts");
389  fhJetPtBefore->SetXTitle("p_{T jet}(GeV/c)");
390  outputContainer->Add(fhJetPtBefore) ;
391 
392  fhJetPtBeforeCut = new TH1F("JetPtBeforeCut","JetPtBeforeCut",150,-50,100);
393  fhJetPtBeforeCut->SetYTitle("Counts");
394  fhJetPtBeforeCut->SetXTitle("p_{T jet}(GeV/c)");
395  outputContainer->Add(fhJetPtBeforeCut) ;
396 
397  fhJetPt = new TH1F("JetPt","JetPt",150,-50,100);
398  fhJetPt->SetYTitle("Counts");
399  fhJetPt->SetXTitle("p_{T jet}(GeV/c)");
400  outputContainer->Add(fhJetPt) ;
401 
402  fhJetPtMostEne = new TH1F("JetPtMostEne","JetPtMostEne",150,0,150);
403  fhJetPtMostEne->SetYTitle("Counts");
404  fhJetPtMostEne->SetXTitle("p_{T jet}(GeV/c)");
405  outputContainer->Add(fhJetPtMostEne) ;
406 
407  fhJetPhi = new TH1F("JetPhi","JetPhi",130,0,6.5);
408  fhJetPhi->SetYTitle("Counts");
409  fhJetPhi->SetXTitle("#phi_{jet}");
410  outputContainer->Add(fhJetPhi) ;
411 
412  fhJetEta = new TH1F("JetEta","JetEta",100,-1,1);
413  fhJetEta->SetYTitle("Counts");
414  fhJetEta->SetXTitle("#eta_{jet}");
415  outputContainer->Add(fhJetEta) ;
416 
417  fhJetEtaVsPt = new TH2F("JetEtaVsPt","JetEtaVsPt",100,0,100,50,-1,1);
418  fhJetEtaVsPt->SetYTitle("#eta_{jet}");
419  fhJetEtaVsPt->SetXTitle("p_{T,jet}(GeV/c)");
420  outputContainer->Add(fhJetEtaVsPt) ;
421 
422  fhJetPhiVsEta = new TH2F("JetPhiVsEta","JetPhiVsEta",65,0,6.5,50,-1,1);
423  fhJetPhiVsEta->SetYTitle("#eta_{jet}");
424  fhJetPhiVsEta->SetXTitle("#phi_{jet}");
425  outputContainer->Add(fhJetPhiVsEta) ;
426 
427  fhJetEtaVsNpartInJet= new TH2F("JetEtaVsNpartInJet","JetEtaVsNpartInJet",50,-1,1,100,0.,200.);
428  fhJetEtaVsNpartInJet->SetYTitle("N_{tracks-in-jet}");
429  fhJetEtaVsNpartInJet->SetXTitle("#eta_{jet}");
430  outputContainer->Add(fhJetEtaVsNpartInJet) ;
431 
432  fhJetEtaVsNpartInJetBkg= new TH2F("JetEtaVsNpartInJetBkg","JetEtaVsNpartInJetBkg",50,-1,1,100,0.,200.);
433  fhJetEtaVsNpartInJetBkg->SetYTitle("N_{tracks-in-jet}");
434  fhJetEtaVsNpartInJetBkg->SetXTitle("#eta_{jet}");
435  outputContainer->Add(fhJetEtaVsNpartInJetBkg) ;
436 
437  fhJetChBkgEnergyVsPt = new TH2F("JetBkgChEnergyVsPt","JetBkgChEnergyVsPt",100,0,100,90,0,90);
438  fhJetChBkgEnergyVsPt->SetYTitle("Jet Bkg Energy (GeV)");
439  fhJetChBkgEnergyVsPt->SetXTitle("p_{T,jet} (GeV/c)");
440  outputContainer->Add(fhJetChBkgEnergyVsPt);
441 
442  fhJetChAreaVsPt = new TH2F("JetChAreaVsPt","JetChAreaVsPt",100,0,100,50,0,1);
443  fhJetChAreaVsPt->SetYTitle("Jet Area");
444  fhJetChAreaVsPt->SetXTitle("p_{T,jet} (GeV/c)");
445  outputContainer->Add(fhJetChAreaVsPt);
446 
447  if(IsHistogramTracks()){
448  fhTrackPhiVsEta = new TH2F("TrackPhiVsEta","TrackPhiVsEta",65,0,6.5,50,-1,1);
449  fhTrackPhiVsEta->SetYTitle("#eta_{track}");
450  fhTrackPhiVsEta->SetXTitle("#phi_{track}");
451  outputContainer->Add(fhTrackPhiVsEta) ;
452 
453  fhTrackAveTrackPt = new TH1F("TrackAveTrackPt","TrackAveTrackPt",45,0,1.5);
454  fhTrackAveTrackPt->SetYTitle("Counts");
455  fhTrackAveTrackPt->SetXTitle("Average p_{T,track} (GeV/c)");
456  outputContainer->Add(fhTrackAveTrackPt);
457 
458  }//end of IsHistogramTracks()
459 
460  for(Int_t i=0;i<10;i++){
461  fhJetNjetOverPtCut[i] = new TH1F(Form("JetNjetOverPtCut%d", i),Form("JetNjetOverPtCut%d", i),100,0,100);
462  fhJetNjetOverPtCut[i]->SetYTitle("Counts");
463  fhJetNjetOverPtCut[i]->SetXTitle("N_{jets} over threshold");
464  outputContainer->Add(fhJetNjetOverPtCut[i]);
465  }
466 
467  fhJetChBkgEnergyVsArea = new TH2F("JetBkgChEnergyVsArea","JetBkgChEnergyVsArea",100,0,100,70,0,0.7);
468  fhJetChBkgEnergyVsArea->SetXTitle("Jet Bkg Energy (GeV)");
469  fhJetChBkgEnergyVsArea->SetYTitle("Area");
470  outputContainer->Add(fhJetChBkgEnergyVsArea);
471 
472  fhJetRhoVsPt = new TH2F("JetRhoVsPt","JetRhoVsPt",100,0,100,100,0,150);
473  fhJetRhoVsPt->SetYTitle("Rho");
474  fhJetRhoVsPt->SetXTitle("p_{T,jet} (GeV/c)");
475  outputContainer->Add(fhJetRhoVsPt);
476 
477  if(IsHistogramJetBkg()){
478  fhJetRhoVsCentrality = new TH2F("JetRhoVsCentrality","JetRhoVsCentrality",100,0,100,100,0,200);
479  fhJetRhoVsCentrality->SetYTitle("Rho");
480  fhJetRhoVsCentrality->SetXTitle("Centrality");
481  outputContainer->Add(fhJetRhoVsCentrality);
482  }
483 
484  fhJetNparticlesInJet = new TH1F("JetNparticlesInJet","JetNparticlesInJet",100,0,200);
485  fhJetNparticlesInJet->SetXTitle("N^{particles}");
486  fhJetNparticlesInJet->SetYTitle("N^{jets}");
487  outputContainer->Add(fhJetNparticlesInJet);
488 
489  fhJetDeltaEtaDeltaPhi = new TH2F("JetDeltaEtaDeltaPhi","#Delta #eta^{jet-track} vs. #Delta #phi^{jet-track} for jet tracks",100,-0.8,0.8,100,-0.8,0.8);
490  fhJetDeltaEtaDeltaPhi->SetXTitle("#Delta #eta^{jet-track}");
491  fhJetDeltaEtaDeltaPhi->SetYTitle("#Delta #phi^{jet-track}");
492  outputContainer->Add(fhJetDeltaEtaDeltaPhi );
493 
494 
495  fhJetDeltaEtaDeltaPhiAllTracks = new TH2F("JetDeltaEtaDeltaPhiAllTracks","#Delta #eta^{jet-track} vs. #Delta #phi^{jet-track} for all tracks",100,-3.2,3.2,100,-3.2,3.2);
496  fhJetDeltaEtaDeltaPhiAllTracks->SetXTitle("#Delta #eta^{jet-track}");
497  fhJetDeltaEtaDeltaPhiAllTracks->SetYTitle("#Delta #phi^{jet-track}");
498  outputContainer->Add(fhJetDeltaEtaDeltaPhiAllTracks);
499 
500 
501  if(IsHistogramJetTracks()){
502  fhJetAveTrackPt = new TH1F("JetAveTrackPt","JetAveTrackPt",45,0.,1.5);
503  fhJetAveTrackPt->SetXTitle("Average p_{T,track} (GeV/c)");
504  fhJetAveTrackPt->SetYTitle("Counts");
505  outputContainer->Add(fhJetAveTrackPt);
506 
507  for(Int_t i=0;i<6;i++){
508  if(i==0) fhJetNtracksInJetAboveThr[i] = new TH2F(Form("JetNtracksInJetAboveThr%d", i),Form("JetNtracksInJetAboveThr%d", i),100,0,100,100,0,200);
509  else fhJetNtracksInJetAboveThr[i] = new TH2F(Form("JetNtracksInJetAboveThr%d", i),Form("JetNtracksInJetAboveThr%d", i),100,0,100,100,0,100);
510  fhJetNtracksInJetAboveThr[i]->SetXTitle("p_{T,jet} (GeV/c)");
511  fhJetNtracksInJetAboveThr[i]->SetYTitle("N_{tracks} over threshold");
512  outputContainer->Add(fhJetNtracksInJetAboveThr[i]);
513  }
514 
515  for(Int_t i=0;i<5;i++){
516  fhJetRatioNTrkAboveToNTrk[i] = new TH2F(Form("JetRatioNTrkAboveToNTrk%d", i),Form("JetRatioNTrkAboveToNTrk%d", i),100,0,100,40,0,1);
517  fhJetRatioNTrkAboveToNTrk[i]->SetXTitle("p_{T,jet} (GeV/c)");
518  fhJetRatioNTrkAboveToNTrk[i]->SetYTitle("Ratio N_{tracks} over threshold to N_{tracks}");
519  outputContainer->Add(fhJetRatioNTrkAboveToNTrk[i]);
520 
521  fhJetNtrackRatioMostEne[i] = new TH2F(Form("JetNtrackRatioMostEne%d", i),Form("JetNtrackRatioMostEne%d", i),100,0,100,40,0,1);
522  fhJetNtrackRatioMostEne[i]->SetXTitle("p_{T,jet} (GeV/c)");
523  fhJetNtrackRatioMostEne[i]->SetYTitle("Ratio N_{tracks} over threshold to N_{tracks}");
524  outputContainer->Add(fhJetNtrackRatioMostEne[i]);
525 
526  fhJetNtrackRatioJet5GeV[i] = new TH2F(Form("JetNtrackRatioJet5GeV%d", i),Form("JetNtrackRatioJet5GeV%d", i),100,0,100,40,0,1);
527  fhJetNtrackRatioJet5GeV[i]->SetXTitle("p_{T,jet} (GeV/c)");
528  fhJetNtrackRatioJet5GeV[i]->SetYTitle("Ratio N_{tracks} over threshold to N_{tracks}");
529  outputContainer->Add(fhJetNtrackRatioJet5GeV[i]);
530 
531  fhJetNtrackRatioLead5GeV[i] = new TH2F(Form("JetNtrackRatioLead5GeV%d", i),Form("JetNtrackRatioLead5GeV%d", i),100,0,100,40,0,1);
532  fhJetNtrackRatioLead5GeV[i]->SetXTitle("p_{T,jet} (GeV/c)");
533  fhJetNtrackRatioLead5GeV[i]->SetYTitle("Ratio N_{tracks} over threshold to N_{tracks}");
534  outputContainer->Add(fhJetNtrackRatioLead5GeV[i]);
535  }
536  }//end of if IsHistogramJetTracks
537 
538  //temporary background jets histograms
539  if(IsHistogramJetBkg()){
540  for(Int_t i=0;i<4;i++){
541  fhBkgJetBackground[i] = new TH1F(Form("BkgJetBackground%d", i),Form("BkgJetBackground%d", i),100,0,200);
542  fhBkgJetBackground[i]->SetXTitle("<#rho> (GeV/c)");
543  fhBkgJetBackground[i]->SetYTitle("Counts");
544  outputContainer->Add(fhBkgJetBackground[i]);
545 
546  fhBkgJetSigma[i] = new TH1F(Form("BkgJetSigma%d", i),Form("BkgJetSigma%d", i),100,0,50);
547  fhBkgJetSigma[i]->SetXTitle("#sigma (GeV/c)");
548  fhBkgJetSigma[i]->SetYTitle("Counts");
549  outputContainer->Add(fhBkgJetSigma[i]);
550 
551  fhBkgJetArea[i] = new TH1F(Form("BkgJetArea%d", i),Form("BkgJetArea%d", i),100,0,1);
552  fhBkgJetArea[i]->SetXTitle("<A>");
553  fhBkgJetArea[i]->SetYTitle("Counts");
554  outputContainer->Add(fhBkgJetArea[i]);
555  }
556  }
557 
558  //temporary photon histograms
559  fhPhotonPtMostEne = new TH1F("PhotonPtMostEne","PhotonPtMostEne",100,0,100);
560  fhPhotonPtMostEne->SetYTitle("Counts");
561  fhPhotonPtMostEne->SetXTitle("p_{T,#gamma} (GeV/c)");
562  outputContainer->Add(fhPhotonPtMostEne);
563 
564 // fhPhotonIndexMostEne = new TH1F("PhotonIndexMostEne","PhotonIndexMostEne",100,0,100);
565 // fhPhotonIndexMostEne->SetYTitle("Counts");
566 // fhPhotonIndexMostEne->SetXTitle("Index");
567 // outputContainer->Add(fhPhotonIndexMostEne);
568 
569  fhPhotonAverageEnergy = new TH1F("PhotonAverageEnergy","PhotonAverageEnergy",100,0,10);
570  fhPhotonAverageEnergy->SetYTitle("Counts");
571  fhPhotonAverageEnergy->SetXTitle("p_{T,#gamma} (GeV/c)");
572  outputContainer->Add(fhPhotonAverageEnergy);
573 
574  fhPhotonRatioAveEneToMostEne = new TH1F("PhotonRatioAveEneToMostEne","PhotonRatioAveEneToMostEne",100,0,1);
575  fhPhotonRatioAveEneToMostEne->SetYTitle("Counts");
576  fhPhotonRatioAveEneToMostEne->SetXTitle("Ratio");
577  outputContainer->Add(fhPhotonRatioAveEneToMostEne);
578 
579  fhPhotonAverageEnergyMinus1 = new TH1F("PhotonAverageEnergyMinus1","PhotonAverageEnergyMinus1",100,0,10);
580  fhPhotonAverageEnergyMinus1->SetYTitle("Counts");
581  fhPhotonAverageEnergyMinus1->SetXTitle("p_{T,#gamma} (GeV/c)");
582  outputContainer->Add(fhPhotonAverageEnergyMinus1);
583 
584  fhPhotonRatioAveEneMinus1ToMostEne = new TH1F("PhotonRatioAveEneMinus1ToMostEne","PhotonRatioAveEneMinus1ToMostEne",100,0,1);
585  fhPhotonRatioAveEneMinus1ToMostEne->SetYTitle("Counts");
586  fhPhotonRatioAveEneMinus1ToMostEne->SetXTitle("Ratio");
587  outputContainer->Add(fhPhotonRatioAveEneMinus1ToMostEne);
588 
589  fhPhotonNgammaMoreAverageToNgamma = new TH1F("PhotonNgammaMoreAverageToNgamma","PhotonNgammaMoreAverageToNgamma",100,0,1);
590  fhPhotonNgammaMoreAverageToNgamma->SetYTitle("Counts");
591  fhPhotonNgammaMoreAverageToNgamma->SetXTitle("Ratio");
592  outputContainer->Add(fhPhotonNgammaMoreAverageToNgamma);
593 
594  fhPhotonNgammaMoreAverageMinus1ToNgamma = new TH1F("PhotonNgammaMoreAverageMinus1ToNgamma","PhotonNgammaMoreAverageMinus1ToNgamma",100,0,1);
595  fhPhotonNgammaMoreAverageMinus1ToNgamma->SetYTitle("Counts");
596  fhPhotonNgammaMoreAverageMinus1ToNgamma->SetXTitle("Ratio");
597  outputContainer->Add(fhPhotonNgammaMoreAverageMinus1ToNgamma);
598 
599  for(Int_t i=0;i<10;i++){
600  fhPhotonNgammaOverPtCut[i] = new TH1F(Form("PhotonNgammaOverPtCut%d",i),Form("PhotonNgammaOverPtCut%d",i),100,0,100);
601  fhPhotonNgammaOverPtCut[i]->SetYTitle("Counts");
602  fhPhotonNgammaOverPtCut[i]->SetXTitle("N_{#gamma} over threshold");
603  outputContainer->Add(fhPhotonNgammaOverPtCut[i]);
604  }
605 
606  fhPhotonBkgRhoVsNtracks = new TH2F("PhotonBkgRhoVsNtracks","PhotonBkgRhoVsNtracks",200,0,2500,75,0,1.5);
607  //fhPhotonBkgRhoVsNtracks->SetXTitle("Counts");
608  fhPhotonBkgRhoVsNtracks->SetXTitle("Ntracks");
609  fhPhotonBkgRhoVsNtracks->SetYTitle("Rho");
610  outputContainer->Add(fhPhotonBkgRhoVsNtracks);
611 
612  fhPhotonBkgRhoVsNclusters = new TH2F("PhotonBkgRhoVsNclusters","PhotonBkgRhoVsNclusters",50,0,100,75,0,1.5);
613  fhPhotonBkgRhoVsNclusters->SetXTitle("Nclusters");
614  fhPhotonBkgRhoVsNclusters->SetYTitle("Rho");
615  outputContainer->Add(fhPhotonBkgRhoVsNclusters);
616 
617  fhPhotonBkgRhoVsCentrality = new TH2F("PhotonBkgRhoVsCentrality","PhotonBkgRhoVsCentrality",100,0,100,75,0,1.5);
618  fhPhotonBkgRhoVsCentrality->SetXTitle("Centrality");
619  fhPhotonBkgRhoVsCentrality->SetYTitle("Rho");
620  outputContainer->Add(fhPhotonBkgRhoVsCentrality);
621 
622  fhPhotonBkgRhoVsNcells = new TH2F("PhotonBkgRhoVsNcells","PhotonBkgRhoVsNcells",100,0,200,75,0,1.5);
623  fhPhotonBkgRhoVsNcells->SetXTitle("N_{cells}");
624  fhPhotonBkgRhoVsNcells->SetYTitle("Rho");
625  outputContainer->Add(fhPhotonBkgRhoVsNcells);
626 
627  fhPhotonPt = new TH1F("PhotonPt","PhotonPt",220,-10,100);
628  fhPhotonPt->SetXTitle("p_{T,#gamma} (GeV/c)");
629  fhPhotonPt->SetYTitle("Counts");
630  outputContainer->Add(fhPhotonPt);
631 
632  fhPhotonPtCorrected = new TH1F("PhotonPtCorrected","PhotonPtCorrected",220,-10,100);
633  fhPhotonPtCorrected->SetXTitle("p_{T,#gamma} (GeV/c)");
634  fhPhotonPtCorrected->SetYTitle("Counts");
635  outputContainer->Add(fhPhotonPtCorrected);
636 
637  fhPhotonPtDiff = new TH1F("PhotonPtDiff","PhotonPtDiff",50,0,10);
638  fhPhotonPtDiff->SetXTitle("p_{T,#gamma} (GeV/c)");
639  fhPhotonPtDiff->SetYTitle("Counts");
640  outputContainer->Add(fhPhotonPtDiff);
641 
642  fhPhotonPtDiffVsNtracks = new TH2F("PhotonPtDiffVsNtracks","PhotonPtDiffVsNtracks",200,0,2500,50,0,5);
643  fhPhotonPtDiffVsNtracks->SetXTitle("N_{tracks}");
644  fhPhotonPtDiffVsNtracks->SetYTitle("<#rho^{#gamma}>*N_{cells}");
645  outputContainer->Add(fhPhotonPtDiffVsNtracks);
646 
647  fhPhotonPtDiffVsNclusters = new TH2F("PhotonPtDiffVsNclusters","PhotonPtDiffVsNclusters",50,0,100,50,0,5);
648  fhPhotonPtDiffVsNclusters->SetXTitle("N_{clusters}");
649  fhPhotonPtDiffVsNclusters->SetYTitle("<#rho^{#gamma}>*N_{cells}");
650  outputContainer->Add(fhPhotonPtDiffVsNclusters);
651 
652  fhPhotonPtDiffVsCentrality = new TH2F("PhotonPtDiffVsCentrality","PhotonPtDiffVsCentrality",100,0,100,50,0,5);
653  fhPhotonPtDiffVsCentrality->SetXTitle("Centrality");
654  fhPhotonPtDiffVsCentrality->SetYTitle("<#rho^{#gamma}>*N_{cells}");
655  outputContainer->Add(fhPhotonPtDiffVsCentrality);
656 
657  fhPhotonPtDiffVsNcells = new TH2F("PhotonPtDiffVsNcells","PhotonPtDiffVsNcells",100,0,200,50,0,5);
658  fhPhotonPtDiffVsNcells->SetXTitle("N_{cells}");
659  fhPhotonPtDiffVsNcells->SetYTitle("<#rho^{#gamma}>*N_{cells}");
660  outputContainer->Add(fhPhotonPtDiffVsNcells);
661 
662 
663  fhPhotonPtCorrectedZoom = new TH1F("PhotonPtCorrectedZoom","PhotonPtCorrectedZoom",200,-5,5);
664  fhPhotonPtCorrectedZoom->SetXTitle("p_{T,#gamma} (GeV/c)");
665  fhPhotonPtCorrectedZoom->SetYTitle("Counts");
666  outputContainer->Add(fhPhotonPtCorrectedZoom);
667 
668  fhPhotonSumPtInCone = new TH1F("PhotonSumPtInCone","PhotonSumPtInCone",100,0,100);
669  fhPhotonSumPtInCone->SetXTitle("#Sigma p_{T,#gamma} (GeV/c)");
670  fhPhotonSumPtInCone->SetYTitle("Counts");
671  outputContainer->Add(fhPhotonSumPtInCone);
672 
673  fhPhotonSumPtCorrectInCone = new TH1F("PhotonSumPtCorrectInCone","PhotonSumPtCorrectInCone",100,-20,80);
674  fhPhotonSumPtCorrectInCone->SetXTitle("#Sigma p_{T,#gamma} (GeV/c)");
675  fhPhotonSumPtCorrectInCone->SetYTitle("Counts");
676  outputContainer->Add(fhPhotonSumPtCorrectInCone);
677 
678  fhPhotonSumPtChargedInCone = new TH1F("PhotonSumPtChargedInCone","PhotonSumPtChargedInCone",100,0,100);
679  fhPhotonSumPtChargedInCone->SetXTitle("#Sigma p_{T,#gamma}^{ch} (GeV/c)");
680  fhPhotonSumPtChargedInCone->SetYTitle("Counts");
681  outputContainer->Add(fhPhotonSumPtChargedInCone);
682 
683 
684  //temporary jet histograms after selection
685  fhSelectedJetPhiVsEta = new TH2F("SelectedJetSelectedPhiVsEta","SelectedJetPhiVsEta",65,0,6.5,50,-1,1);
686  fhSelectedJetPhiVsEta->SetYTitle("#eta_{jet}");
687  fhSelectedJetPhiVsEta->SetXTitle("#phi_{jet}");
688  outputContainer->Add(fhSelectedJetPhiVsEta) ;
689 
690  fhSelectedJetChBkgEnergyVsPtJet = new TH2F("SelectedJetBkgChEnergyVsPtJet","SelectedJetBkgChEnergyVsPtJet",100,0,100,90,0,90);
691  fhSelectedJetChBkgEnergyVsPtJet->SetYTitle("Jet Bkg Energy (GeV)");
692  fhSelectedJetChBkgEnergyVsPtJet->SetXTitle("p_{T,jet} (GeV/c)");
693  outputContainer->Add(fhSelectedJetChBkgEnergyVsPtJet);
694 
695  fhSelectedJetChAreaVsPtJet = new TH2F("SelectedJetChAreaVsPtJet","SelectedJetChAreaVsPtJet",100,0,100,50,0,1);
696  fhSelectedJetChAreaVsPtJet->SetYTitle("Jet Area");
697  fhSelectedJetChAreaVsPtJet->SetXTitle("p_{T,jet} (GeV/c)");
698  outputContainer->Add(fhSelectedJetChAreaVsPtJet);
699 
700  fhSelectedJetNjet = new TH1F("SelectedJetNjet","SelectedJetNjet",100,0,100);
701  fhSelectedJetNjet->SetYTitle("Counts");
702  fhSelectedJetNjet->SetXTitle("N_{jets} per event");
703  outputContainer->Add(fhSelectedJetNjet);
704 
705  fhSelectedNtracks = new TH1F("SelectedNtracks","SelectedNtracks",100,0,2000);
706  fhSelectedNtracks->SetYTitle("Counts");
707  fhSelectedNtracks->SetXTitle("N_{tracks} per event");
708  outputContainer->Add(fhSelectedNtracks);
709 
710  fhSelectedTrackPhiVsEta = new TH2F("SelectedTrackPhiVsEta","SelectedTrackPhiVsEta",65,0,6.5,50,-1,1);
711  fhSelectedTrackPhiVsEta->SetYTitle("#eta_{track}");
712  fhSelectedTrackPhiVsEta->SetXTitle("#phi_{track}");
713  outputContainer->Add(fhSelectedTrackPhiVsEta) ;
714 
715  fhCuts2 = new TH1F("Cuts2","Cuts2",10,0.,10.);
716  fhCuts2->SetYTitle("Counts");
717  fhCuts2->SetXTitle("Cut number");
718  outputContainer->Add(fhCuts2);
719 
720  fhSelectedPhotonNLMVsPt = new TH2F("SelectedPhotonNLMVsPt","SelectedPhotonNLMVsPt",100,0,100,10,0,10);
721  fhSelectedPhotonNLMVsPt->SetYTitle("NLM");
722  fhSelectedPhotonNLMVsPt->SetXTitle("p_{T,#gamma} (GeV/c)");
723  outputContainer->Add(fhSelectedPhotonNLMVsPt);
724 
725  fhSelectedPhotonLambda0VsPt = new TH2F("SelectedPhotonLambda0VsPt","SelectedPhotonLambda0VsPt",100,0,100,50,0,5);
726  fhSelectedPhotonLambda0VsPt->SetYTitle("#lambda_{0}");
727  fhSelectedPhotonLambda0VsPt->SetXTitle("p_{T,#gamma} (GeV/c)");
728  outputContainer->Add(fhSelectedPhotonLambda0VsPt);
729 
730  //random
731  fhRandomPhiEta[0] = new TH2F("RandomPhiEtaRC", "RandomPhiEtaRC", 100,0,6.5,100,-1.,1.);
732  fhRandomPhiEta[1] = new TH2F("RandomPhiEtaPCG", "RandomPhiEtaPerpConePhoton",100,0,6.5,100,-1.,1.);
733  fhRandomPhiEta[2] = new TH2F("RandomPhiEtaPCJ", "RandomPhiEtaPerpConeJet", 100,0,6.5,100,-1.,1.);
734  fhRandomPhiEta[3] = new TH2F("RandomPhiEtaMP", "RandomPhiEtaMidPoint", 100,0,6.5,100,-1.,1.);
735  fhRandomPhiEta[4] = new TH2F("RandomPhiEtaTest","RandomPhiEtaTest", 100,0,6.5,100,-1.,1.);
736  for(Int_t i=0;i<5;i++){
737  fhRandomPhiEta[i]->SetXTitle("#phi");
738  fhRandomPhiEta[i]->SetYTitle("#eta");
739  outputContainer->Add(fhRandomPhiEta[i]);
740  }
741 
742  //MC generated photons and jets information
743  if(fMCStudies) {
744  fhMCPhotonCuts = new TH1F("MCPhotonCuts","MCPhotonCuts",10,0.,10.);
745  fhMCPhotonCuts->SetYTitle("Counts");
746  fhMCPhotonCuts->SetXTitle("Cut number");
747  outputContainer->Add(fhMCPhotonCuts);
748 
749  fhMCPhotonPt = new TH1F("MCPhotonPt","MCPhotonPt",100,0.,100.);
750  fhMCPhotonPt->SetYTitle("Counts");
751  fhMCPhotonPt->SetXTitle("p_{T,#gamma}^{gen} (GeV/c)");
752  outputContainer->Add(fhMCPhotonPt);
753 
754  fhMCPhotonEtaPhi = new TH2F("MCPhotonEtaPhi","MCPhotonEtaPhi",65,0,6.5,50,-1,1);
755  fhMCPhotonEtaPhi->SetYTitle("#eta_{#gamma}");
756  fhMCPhotonEtaPhi->SetXTitle("#phi_{#gamma}");
757  outputContainer->Add(fhMCPhotonEtaPhi) ;
758 
759  fhMCJetOrigin = new TH1F("MCJetOrigin","MCJetOrigin",35,-10.,25.);
760  fhMCJetOrigin->SetYTitle("Counts");
761  fhMCJetOrigin->SetXTitle("PDG number");
762  outputContainer->Add(fhMCJetOrigin);
763 
764  fhMCJetNPartVsPt = new TH2F("MCJetNPartVsPt","MCJetNPartVsPt",100,0.,100.,100,0.,100.);
765  fhMCJetNPartVsPt->SetYTitle("N_{particles}");
766  fhMCJetNPartVsPt->SetXTitle("p_{T,full-jet}^{gen} (GeV/c)");
767  outputContainer->Add(fhMCJetNPartVsPt) ;
768 
769  fhMCJetChNPartVsPt = new TH2F("MCJetChNPartVsPt","MCJetChNPartVsPt",100,0.,100.,100,0.,100.);
770  fhMCJetChNPartVsPt->SetYTitle("N_{particles}");
771  fhMCJetChNPartVsPt->SetXTitle("p_{T,charged-jet}^{gen} (GeV/c)");
772  outputContainer->Add(fhMCJetChNPartVsPt) ;
773 
774  fhMCJetNPart150VsPt = new TH2F("MCJetNPart150VsPt","MCJetNPart150VsPt",100,0.,100.,100,0.,100.);
775  fhMCJetNPart150VsPt->SetYTitle("N_{particles (p_{T}>150 MeV/c)}");
776  fhMCJetNPart150VsPt->SetXTitle("p_{T,full-jet}^{gen} (GeV/c)");
777  outputContainer->Add(fhMCJetNPart150VsPt) ;
778 
779  fhMCJetChNPart150VsPt = new TH2F("MCJetChNPart150VsPt","MCJetChNPart150VsPt",100,0.,100.,100,0.,100.);
780  fhMCJetChNPart150VsPt->SetYTitle("N_{particles (p_{T}>150 MeV/c)}");
781  fhMCJetChNPart150VsPt->SetXTitle("p_{T,charged-jet}^{gen} (GeV/c)");
782  outputContainer->Add(fhMCJetChNPart150VsPt) ;
783 
784  fhMCJetChNPart150ConeVsPt = new TH2F("MCJetChNPart150ConeVsPt","MCJetChNPart150ConeVsPt",100,0.,100.,50,0.,50.);
785  fhMCJetChNPart150ConeVsPt->SetYTitle("N_{particles (p_{T}>150 MeV/c)}");
786  fhMCJetChNPart150ConeVsPt->SetXTitle("p_{T,charged-jet,R=0.4}^{gen} (GeV/c)");
787  outputContainer->Add(fhMCJetChNPart150ConeVsPt) ;
788 
789  fhMCJetRatioChFull = new TH1F("MCJetRatioChFull","MCJetRatioChFull",100,0.,1.);
790  fhMCJetRatioChFull->SetYTitle("Counts");
791  fhMCJetRatioChFull->SetXTitle("p_{T,charged-jet}^{gen}/p_{T,full-jet}^{gen}");
792  outputContainer->Add(fhMCJetRatioChFull);
793 
794  fhMCJetRatioCh150Ch = new TH1F("MCJetRatioCh150Ch","MCJetRatioCh150Ch",100,0.,1.);
795  fhMCJetRatioCh150Ch->SetYTitle("Counts");
796  fhMCJetRatioCh150Ch->SetXTitle("p_{T,charged-jet,pT>150MeV/c}^{gen}/p_{T,charged-jet}^{gen}");
797  outputContainer->Add(fhMCJetRatioCh150Ch);
798 
799  fhMCJetEtaPhi = new TH2F("MCJetEtaPhi","MCJetEtaPhi",65,0,6.5,50,-1,1);
800  fhMCJetEtaPhi->SetYTitle("#eta_{jet}");
801  fhMCJetEtaPhi->SetXTitle("#phi_{jet}");
802  outputContainer->Add(fhMCJetEtaPhi) ;
803 
804  fhMCJetChEtaPhi = new TH2F("MCJetChEtaPhi","MCJetChEtaPhi",65,0,6.5,50,-1,1);
805  fhMCJetChEtaPhi->SetYTitle("#eta_{jet}");
806  fhMCJetChEtaPhi->SetXTitle("#phi_{jet}");
807  outputContainer->Add(fhMCJetChEtaPhi) ;
808 
809  fhMCJet150EtaPhi = new TH2F("MCJet150EtaPhi","MCJet150EtaPhi",65,0,6.5,50,-1,1);
810  fhMCJet150EtaPhi->SetYTitle("#eta_{jet}");
811  fhMCJet150EtaPhi->SetXTitle("#phi_{jet}");
812  outputContainer->Add(fhMCJet150EtaPhi) ;
813 
814  fhMCJetCh150EtaPhi = new TH2F("MCJetCh150EtaPhi","MCJetCh150EtaPhi",65,0,6.5,50,-1,1);
815  fhMCJetCh150EtaPhi->SetYTitle("#eta_{jet}");
816  fhMCJetCh150EtaPhi->SetXTitle("#phi_{jet}");
817  outputContainer->Add(fhMCJetCh150EtaPhi) ;
818 
819  fhMCJetCh150ConeEtaPhi = new TH2F("MCJetCh150ConeEtaPhi","MCJetCh150ConeEtaPhi",65,0,6.5,50,-1,1);
820  fhMCJetCh150ConeEtaPhi->SetYTitle("#eta_{jet}");
821  fhMCJetCh150ConeEtaPhi->SetXTitle("#phi_{jet}");
822  outputContainer->Add(fhMCJetCh150ConeEtaPhi) ;
823  }
824 
825  //tree with data
826  if(fSaveGJTree){
827  fTreeGJ=new TTree("fTreeGJ","fTreeGJ");
828  fTreeGJ->Branch("fGamPt" ,&fGamPt ,"fGamPt/D");
829  fTreeGJ->Branch("fGamLambda0" ,&fGamLambda0 ,"fGamLambda0/D");
830  fTreeGJ->Branch("fGamNLM" ,&fGamNLM ,"fGamNLM/I");
831  fTreeGJ->Branch("fGamSumPtCh" ,&fGamSumPtCh ,"fGamSumPtCh/D");
832  fTreeGJ->Branch("fGamNtracks" ,&fGamNtracks ,"fGamNtracks/I");
833  fTreeGJ->Branch("fGamTime" ,&fGamTime ,"fGamTime/D");
834  fTreeGJ->Branch("fGamNcells" ,&fGamNcells ,"fGamNcells/I");
835  fTreeGJ->Branch("fGamEta" ,&fGamEta ,"fGamEta/D");
836  fTreeGJ->Branch("fGamPhi" ,&fGamPhi ,"fGamPhi/D");
837  fTreeGJ->Branch("fGamSumPtNeu",&fGamSumPtNeu ,"fGamSumPtNeu/D");
838  fTreeGJ->Branch("fGamNclusters" ,&fGamNclusters ,"fGamNclusters/I");
839  fTreeGJ->Branch("fGamAvEne" ,&fGamAvEne ,"fGamAvEne/D");
840  fTreeGJ->Branch("fJetPhi" ,&fJetPhi ,"fJetPhi/D");
841  fTreeGJ->Branch("fJetEta" ,&fJetEta ,"fJetEta/D");
842  fTreeGJ->Branch("fJetPt" ,&fJetPt ,"fJetPt/D");
843  fTreeGJ->Branch("fJetBkgChEne",&fJetBkgChEne ,"fJetBkgChEne/D");
844  fTreeGJ->Branch("fJetArea" ,&fJetArea ,"fJetArea/D");
845  fTreeGJ->Branch("fJetNtracks" ,&fJetNtracks ,"fJetNtracks/I");
846  fTreeGJ->Branch("fJetNtracks1" ,&fJetNtracks1 ,"fJetNtracks1/I");
847  fTreeGJ->Branch("fJetNtracks2" ,&fJetNtracks2 ,"fJetNtracks2/I");
848  fTreeGJ->Branch("fJetRho" ,&fJetRho ,"fJetRho/D");
849  fTreeGJ->Branch("fEventNumber",&fEventNumber ,"fEventNumber/I");
850  fTreeGJ->Branch("fNtracks" ,&fNtracks ,"fNtracks/I");
851  fTreeGJ->Branch("fZvertex" ,&fZvertex ,"fZvertex/D");
852  fTreeGJ->Branch("fCentrality" ,&fCentrality ,"fCentrality/D");
853  fTreeGJ->Branch("fIso" ,&fIso ,"fIso/O");
854  fTreeGJ->Branch("fGamRho" ,&fGamRho ,"fGamRho/D");
855 
856  fTreeGJ->Branch("fMCGamPt" ,&fMCGamPt ,"fMCGamPt/D" );
857  fTreeGJ->Branch("fMCGamEta" ,&fMCGamEta ,"fMCGamEta/D" );
858  fTreeGJ->Branch("fMCGamPhi" ,&fMCGamPhi ,"fMCGamPhi/D" );
859  fTreeGJ->Branch("fMCPartonType" ,&fMCPartonType ,"fMCPartonType/I" );
860  fTreeGJ->Branch("fMCJetPt" ,&fMCJetPt ,"fMCJetPt/D" );
861  fTreeGJ->Branch("fMCJetChPt" ,&fMCJetChPt ,"fMCJetChPt/D" );
862  fTreeGJ->Branch("fMCJet150Pt" ,&fMCJet150Pt ,"fMCJet150Pt/D" );
863  fTreeGJ->Branch("fMCJetCh150Pt" ,&fMCJetCh150Pt ,"fMCJetCh150Pt/D" );
864  fTreeGJ->Branch("fMCJetNPart" ,&fMCJetNPart ,"fMCJetNPart/I" );
865  fTreeGJ->Branch("fMCJetChNPart" ,&fMCJetChNPart ,"fMCJetChNPart/I" );
866  fTreeGJ->Branch("fMCJet150NPart" ,&fMCJet150NPart ,"fMCJet150NPart/I" );
867  fTreeGJ->Branch("fMCJetCh150NPart" ,&fMCJetCh150NPart ,"fMCJetCh150NPart/I");
868  fTreeGJ->Branch("fMCJetEta" ,&fMCJetEta ,"fMCJetEta/D" );
869  fTreeGJ->Branch("fMCJetPhi" ,&fMCJetPhi ,"fMCJetPhi/D" );
870  fTreeGJ->Branch("fMCJetChEta" ,&fMCJetChEta ,"fMCJetChEta/D" );
871  fTreeGJ->Branch("fMCJetChPhi" ,&fMCJetChPhi ,"fMCJetChPhi/D" );
872  fTreeGJ->Branch("fMCJet150Eta" ,&fMCJet150Eta ,"fMCJet150Eta/D" );
873  fTreeGJ->Branch("fMCJet150Phi" ,&fMCJet150Phi ,"fMCJet150Phi/D" );
874  fTreeGJ->Branch("fMCJetCh150Eta" ,&fMCJetCh150Eta ,"fMCJetCh150Eta/D" );
875  fTreeGJ->Branch("fMCJetCh150Phi" ,&fMCJetCh150Phi ,"fMCJetCh150Phi/D" );
876 
877  fTreeGJ->Branch("fMCJetCh150ConePt" ,&fMCJetCh150ConePt ,"fMCJetCh150ConePt/D" );
878  fTreeGJ->Branch("fMCJetCh150ConeNPart" ,&fMCJetCh150ConeNPart ,"fMCJetCh150ConeNPart/I");
879  fTreeGJ->Branch("fMCJetCh150ConeEta" ,&fMCJetCh150ConeEta ,"fMCJetCh150ConeEta/D" );
880  fTreeGJ->Branch("fMCJetCh150ConePhi" ,&fMCJetCh150ConePhi ,"fMCJetCh150ConePhi/D" );
881 
882  outputContainer->Add(fTreeGJ);
883  }
884 
885  return outputContainer;
886 }
887 
888 //_______________________________________________________
890 //_______________________________________________________
892 {
893  SetInputAODName("PWG4Particle");
894  AddToHistogramsName("AnaJetFinderCorr_");
895 
896  fDeltaPhiMinCut = 2.6 ;
897  fDeltaPhiMaxCut = 3.7 ;
898  fRatioMaxCut = 1.2 ;
899  fRatioMinCut = 0.3 ;
900  fConeSize = 0.4 ;
901  fPtThresholdInCone = 0.5 ;
902  fUseJetRefTracks = kFALSE ;
904  fSelectIsolated = kTRUE;
905  fJetConeSize = 0.4 ;
906  fJetMinPt = 15. ; //GeV/c
907  fJetMinPtBkgSub = -100. ;//GeV/c
908  fJetAreaFraction = 0.6 ;
909  fJetBranchName = "jets";
910  fBkgJetBranchName = "jets";
911  fGammaConeSize = 0.4;
913  fSaveGJTree = kTRUE;
914  fMostEnergetic = kFALSE;
915  fMostOpposite = kTRUE;
916  fUseHistogramJetBkg = kTRUE;
917  fUseHistogramTracks = kTRUE;
918  fUseHistogramJetTracks = kTRUE;
919  fMCStudies = kFALSE;
920 }
921 
922 //__________________________________________________________________
925 //__________________________________________________________________
926 Int_t AliAnaParticleJetFinderCorrelation::SelectJet(AliAODPWG4Particle * particle, TClonesArray *aodRecJets)
927 {
928  Double_t particlePt=particle->Pt();
930  particlePt-=(fGamRho*particle->GetNCells());
931 
932 // Int_t clusterIDtmp = particle->GetCaloLabel(0) ;
933 // Int_t nCells=0;
934 // AliVCluster *cluster=0;
935 // if(!(clusterIDtmp<0) ){
936 // Int_t iclustmp = -1;
937 // TObjArray* clusters = GetEMCALClusters();
938 // cluster = FindCluster(clusters,clusterIDtmp,iclustmp);
939 // nCells = cluster->GetNCells();
940 // }
941 // particlePt-=(fGamRho*nCells);
942  }
943  if(particlePt<=0) {
944  //printf("Particle with negative or 0 pt\n");
945  return -2;
946  }
947 
948  Int_t njets = aodRecJets->GetEntriesFast();
949  AliAODJet * jet = 0 ;
950  Int_t index = -1;
951  Double_t dphiprev= 10000;
952  Double_t particlePhi=particle->Phi();
953  Double_t deltaPhi=-10000.;// in the range (0; 2*pi)
954  Double_t jetPt=0.;
955 
956  Bool_t photonOnlyOnce=kTRUE;
957 
958  for(Int_t ijet = 0; ijet < njets ; ijet++)
959  {
960  jet = dynamic_cast<AliAODJet*>(aodRecJets->At(ijet));
961 
962  if(!jet)
963  {
964  AliInfo("Jet not in container");
965  continue;
966  }
967 
968  fhCuts2->Fill(2., GetEventWeight());
969 
970  jetPt=jet->Pt();
971  if(jetPt<fJetMinPt) continue;
972 
973  fhCuts2->Fill(3., GetEventWeight());
974  //put jet eta requirement here |eta_jet|<0.9-jet_cone_size
975  if(TMath::Abs(jet->Eta()) > (0.9 - fJetConeSize) ) continue;
976 
977  fhCuts2->Fill(4., GetEventWeight());
978 
979  if(jet->EffectiveAreaCharged()<fJetAreaFraction*TMath::Pi()*fJetConeSize*fJetConeSize) continue;
980 
981  fhCuts2->Fill(5., GetEventWeight());
982 
984  {
985  jetPt-= (fJetRho * jet->EffectiveAreaCharged() );
986  }
987 
988  if(jetPt<fJetMinPtBkgSub) continue;
989 
990  fhCuts2->Fill(6., GetEventWeight());
991 
992  //printf("jet found\n");
993  Double_t deltaPhi0pi = TMath::Abs(particle->Phi()-jet->Phi());
994  Double_t ratio = jetPt/particlePt;
995 
996  deltaPhi = particlePhi - jet->Phi() ;
997  if ( deltaPhi0pi > TMath::Pi() ) deltaPhi0pi = 2. * TMath::Pi() - deltaPhi0pi ;
998  if(deltaPhi<0) deltaPhi +=(TMath::Pi()*2.);
999 
1000  // new histogram for Leticia x-check
1001  // isolated photon + jet(s)
1002  if(OnlyIsolated() && !particle->IsIsolated() &&
1003  (deltaPhi > fDeltaPhiMinCut) && (deltaPhi < fDeltaPhiMaxCut) )
1004  {
1005  //fill 1D photon + 2D photon+jets
1006  if(photonOnlyOnce)
1007  {
1008  fhGamPtPerTrig->Fill(particlePt, GetEventWeight());
1009  photonOnlyOnce=kFALSE;
1010  }
1011 
1012  fhPtGamPtJet->Fill(particlePt, jetPt, GetEventWeight());
1013  }
1014 
1015 
1016  fhDeltaPtBefore ->Fill(particlePt, particlePt - jetPt , GetEventWeight());
1017  fhDeltaPhiBefore->Fill(particlePt, deltaPhi , GetEventWeight());
1018  fhDeltaEtaBefore->Fill(particlePt, particle->Eta() - jet->Eta(), GetEventWeight());
1019  fhPtRatioBefore ->Fill(particlePt, jetPt/particlePt , GetEventWeight());
1020  fhPtBefore ->Fill(particlePt, jetPt , GetEventWeight());
1021 
1022  fhDeltaPhi0PiCorrectBefore->Fill(particlePt, deltaPhi0pi , GetEventWeight());//good
1023 
1024  AliDebug(5,Form("Jet %d, Ratio pT %2.3f, Delta phi %2.3f",ijet,ratio,deltaPhi));
1025 
1026  if((deltaPhi > fDeltaPhiMinCut) && (deltaPhi < fDeltaPhiMaxCut) &&
1027  (ratio > fRatioMinCut) && (ratio < fRatioMaxCut) &&
1028  (TMath::Abs(deltaPhi-TMath::Pi()) < TMath::Abs(dphiprev-TMath::Pi()) )){//In case there is more than one jet select the most opposite.
1029  dphiprev = deltaPhi;
1030  index = ijet ;
1031  }//Selection cuts
1032  }//AOD Jet loop
1033 
1034  return index ;
1035 }
1036 
1037 //__________________________________________________________________
1039 //__________________________________________________________________
1041 {
1042  // Get the event, check if there are AODs available, if not, skip this analysis
1043  AliAODEvent * event = NULL;
1044 
1045  // cout<<"GetReader()->GetOutputEvent() "<<GetReader()->GetOutputEvent()<<endl;
1046  // cout<<"GetReader()->GetDataType() "<<GetReader()->GetDataType() <<endl;
1047  // cout<<"AliCaloTrackReader::kAOD "<<AliCaloTrackReader::kAOD<<endl;
1048  // cout<<"GetReader()->GetInputEvent() "<<GetReader()->GetInputEvent()<<endl;
1049 
1050  if(GetReader()->GetOutputEvent())
1051  {
1052  event = dynamic_cast<AliAODEvent*>(GetReader()->GetOutputEvent());
1053  }
1055  {
1056  event = dynamic_cast<AliAODEvent*>(GetReader()->GetInputEvent());
1057  }
1058  else
1059  {
1060  AliDebug(1,"There are no jets available for this analysis");
1061  return;
1062  }
1063 
1064  if(!GetInputAODBranch() || !event)
1065  {
1066  AliFatal(Form("No input particles in AOD with name branch < %s > \n",
1067  GetInputAODName().Data()));
1068  return; // Trick coverity
1069  }
1070 
1071  if(strcmp(GetInputAODBranch()->GetClass()->GetName(), "AliAODPWG4ParticleCorrelation"))
1072  {
1073  AliFatal(Form("Wrong type of AOD object, change AOD class name in input AOD: It should be <AliAODPWG4ParticleCorrelation> and not <%s>",
1074  GetInputAODBranch()->GetClass()->GetName()));
1075  return; // Trick coverity
1076  }
1077 
1078  //
1079  // non-standard jets
1080  //
1081  Int_t nJets=-1;
1082  TClonesArray *aodRecJets = 0;
1083  //if(IsNonStandardJetFromReader()){//jet branch from reader
1084  AliDebug(3,Form("GetNonStandardJets function (from reader) is called"));
1085  aodRecJets = GetNonStandardJets();
1086  AliDebug(3,Form("aodRecJets %p",aodRecJets));
1087  if(aodRecJets==0x0)
1088  {
1089  if(GetDebug() > 3) event->Print();
1090  AliFatal("List of jets is null");
1091  return;
1092  }
1093  nJets=aodRecJets->GetEntries();
1094  AliDebug(3,Form("nJets %d",nJets));
1095  //}
1096  //end of new part
1097 
1098  if(nJets==0) {
1099  //printf("Why number of jets = 0? Check what type of collision it is. If PbPb -problem.\n");
1100  return;
1101  }
1102 
1103  //
1104  // Background jets
1105  //
1106  AliAODJetEventBackground* aodBkgJets = 0;
1107  if(IsBackgroundJetFromReader()){//jet branch from reader
1108  AliDebug(3,"GetBackgroundJets function is called");
1109  aodBkgJets = GetBackgroundJets();
1110  AliDebug(3,Form("aodBkgJets %p",aodBkgJets));
1111  if(aodBkgJets==0x0){
1112  if(GetDebug() > 3) event->Print();
1113  AliFatal("No jet background found\n");
1114  return; // Trick coverity
1115  }
1116  if(GetDebug() > 3) aodBkgJets->Print("c");
1117  }
1118 
1119  Double_t rhoEvent=0.;
1120  if(aodBkgJets && IsBackgroundJetFromReader() )
1121  {
1122  rhoEvent = aodBkgJets->GetBackground(2);//hardcoded
1123  }
1124 
1125  fJetRho = rhoEvent;
1126 
1127 
1128  Int_t ntrig = GetInputAODBranch()->GetEntriesFast() ;
1129 
1130  AliDebug(3,"Begin jet finder correlation analysis, fill AODs");
1131  AliDebug(3,Form("In particle branch aod entries %d\n", ntrig));
1132  AliDebug(3,Form("In standard jet branch aod entries %d\n", event->GetNJets()));
1133  AliDebug(3,Form("In non standard jet branch aod entries %d\n", nJets));
1134 
1135  //if(nJets==0) return;//to speed up
1136  // cout<<"ntrig po return "<<ntrig<<endl;
1137 
1138  //
1139  //calculate average cell energy without most energetic photon
1140  //
1141  Double_t medianPhotonRho=0.;
1142  //TObjArray* clusters = GetEMCALClusters();
1143  //Int_t clusterIDtmp;
1144  //Int_t iclustmp = -1;
1145  //AliVCluster *cluster=0;
1146 
1148  //
1149  // Find most energetic photon without background subtraction
1150  //
1151  Double_t maxPt=0.;
1152  Int_t maxIndex=-1;
1153  AliAODPWG4ParticleCorrelation* particlecorr =0;
1154  for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1155  particlecorr = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1156  if(particlecorr->Pt() > maxPt) {
1157  maxPt = particlecorr->Pt();
1158  maxIndex = iaod;
1159  }
1160  }
1161 
1162  //
1163  // calculate background energy per cell
1164  //
1165  Int_t numberOfcells=0;
1166  if(ntrig>1){
1167  Double_t *photonRhoArr=new Double_t[ntrig-1];
1168  Int_t photonRhoArrayIndex=0;
1169  for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1170  particlecorr = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1171  if(iaod==maxIndex) continue;
1172 // clusterIDtmp = particlecorr->GetCaloLabel(0) ;
1173 // if(clusterIDtmp < 0) continue;
1174 // cluster = FindCluster(clusters,clusterIDtmp,iclustmp);
1175 // photonRhoArr[photonRhoArrayIndex]=particlecorr->Pt()/ cluster->GetNCells();
1176 // numberOfcells+=cluster->GetNCells();
1177  photonRhoArr[photonRhoArrayIndex]=particlecorr->Pt()/ particlecorr->GetNCells();
1178  numberOfcells+=particlecorr->GetNCells();
1179 
1180  photonRhoArrayIndex++;
1181  }
1182  if(photonRhoArrayIndex>0) medianPhotonRho=TMath::Median(photonRhoArrayIndex,photonRhoArr);
1183  delete [] photonRhoArr;
1184  }
1185  }//end of if background calculation for gamma
1186 
1187  fGamRho = medianPhotonRho;
1188 
1189  //
1190  //take most energetic photon and most energetic jet and combine
1191  //
1192  if(fMostEnergetic){
1193  //
1194  // most energetic photon with background subtraction
1195  //
1196  Double_t mostEnePhotonPt=0.;
1197  Int_t indexMostEnePhoton=-1;
1198  AliAODPWG4ParticleCorrelation* particle =0;
1199  Double_t ptCorrect=0.;
1200 // Int_t nCells=0;
1201  for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1202  particle = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1203 // clusterIDtmp = particle->GetCaloLabel(0) ;
1204 // if(!(clusterIDtmp<0)){
1205 // cluster = FindCluster(clusters,clusterIDtmp,iclustmp);
1206 // nCells = cluster->GetNCells();
1207 // }
1208 // ptCorrect = particle->Pt() - medianPhotonRho * nCells;
1209  ptCorrect = particle->Pt() - medianPhotonRho * particle->GetNCells();
1210 
1211  if( ptCorrect > mostEnePhotonPt ){
1212  mostEnePhotonPt = ptCorrect;
1213  indexMostEnePhoton = iaod ;
1214  }
1215  }
1216  // printf ("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - Photon with index %d selected \n",indexMostEnePhoton);
1217  //
1218  // most energetic jet with background subtraction
1219  //
1220  Double_t mostEneJetPt=0.;
1221  Int_t indexMostEneJet=-1;
1222  AliAODJet * jet = 0 ;
1223  //Double_t ptCorrect=0.;
1224  for(Int_t ijet = 0; ijet < nJets ; ijet++){
1225  jet = dynamic_cast<AliAODJet*>(aodRecJets->At(ijet));
1226  if(!jet) continue;
1227  if(jet->Pt()<fJetMinPt) continue;
1228  if(TMath::Abs(jet->Eta()) > (0.9 - fJetConeSize) ) continue;
1229  if(jet->EffectiveAreaCharged()<fJetAreaFraction*TMath::Pi()*fJetConeSize*fJetConeSize) continue;
1230  ptCorrect = jet->Pt() - rhoEvent * jet->EffectiveAreaCharged();
1231  if(ptCorrect > mostEneJetPt){
1232  mostEneJetPt = ptCorrect;
1233  indexMostEneJet = ijet;
1234  }
1235  }
1236  // printf ("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - Jet with index %d selected \n",indexMostEneJet);
1237 
1238  //
1239  // assign jet to photon
1240  //
1241  if(indexMostEneJet>=0 && indexMostEnePhoton>=0){
1242  particle = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(indexMostEnePhoton));
1243  jet = dynamic_cast<AliAODJet*>(aodRecJets-> At(indexMostEneJet));
1244  if(jet)particle->SetRefJet(jet);
1245  }
1246  }//end of take most energetic photon and most ene. jet after bkg subtraction
1247 
1248  if(fMostOpposite){
1249  //Bool_t isJetFound=kFALSE;//new
1250  //Loop on stored AOD particles, trigger
1251  for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1252  AliAODPWG4ParticleCorrelation* particle = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1253 
1254  //Correlate with jets
1255  Int_t ijet = SelectJet(particle,aodRecJets);//input for jets is TClonesArray
1256  if(ijet > -1){
1257  //isJetFound=kTRUE;
1258  AliDebug(2,Form("Jet with index %d selected",ijet));
1259  AliAODJet *jet = dynamic_cast<AliAODJet*>(aodRecJets-> At(ijet));
1260  if(jet)particle->SetRefJet(jet);
1261  //printf("Most opposite found\n");
1262  }
1263  } // input aod loop
1264  // if(GetReader()->WriteDeltaAODToFile() && isJetFound) WriteJetsToOutputBranch(aodRecJets);
1265  }//end of take most opposite photon and jet after bkg subtraction
1266 
1267  AliDebug(1," End fill AODs \n");
1268 }
1269 
1270 //__________________________________________________________________
1272 //__________________________________________________________________
1274 {
1275  AliDebug(3,"I use MakeAnalysisFillHistograms");
1276  AliDebug(3,Form("ntracks before iso %d\n",GetCTSTracks()->GetEntriesFast()));
1277 
1278  // Get the event, check if there are AODs available, if not, skip this analysis
1279  AliAODEvent * event = NULL;
1280 
1281  //printf("GetReader()->GetOutputEvent() %d\n",GetReader()->GetOutputEvent() );
1282  //printf("GetReader()->GetDataType() %d\n",GetReader()->GetDataType() );
1283  //printf("AliCaloTrackReader::kAOD %d\n",AliCaloTrackReader::kAOD );
1284  //printf("GetReader()->GetInputEvent() %d\n",GetReader()->GetInputEvent() );
1285 
1286  if(GetReader()->GetOutputEvent())
1287  {
1288  event = dynamic_cast<AliAODEvent*>(GetReader()->GetOutputEvent());
1289  }
1291  {
1292  event = dynamic_cast<AliAODEvent*>(GetReader()->GetInputEvent());
1293  }
1294  else
1295  {
1296  AliDebug(3,"There are no jets available for this analysis");
1297  return;
1298  }
1299 
1300  if(!GetInputAODBranch() || !event)
1301  {
1302  AliFatal(Form("No input particles in AOD with name branch < %s >",
1303  GetInputAODName().Data()));
1304  return; // Trick coverity
1305  }
1306 
1307  Int_t nJets=-1;
1308  TClonesArray *aodRecJets = 0;
1309  //if(IsNonStandardJetFromReader()){//branch read in reader from reader
1310  AliDebug(3,"GetNonStandardJets function (from reader) is called");
1311  aodRecJets = GetNonStandardJets();
1312  if(aodRecJets==0x0)
1313  {
1314  if(GetDebug() > 3) event->Print();
1315  AliFatal("Jets container not found\n");
1316  return; // trick coverity
1317  }
1318  nJets=aodRecJets->GetEntries();
1319  //}
1320  if(nJets==0)
1321  {
1322  // printf("Why number of jets = 0? Check what type of collision it is. If PbPb -problem.\n");
1324  aodRecJets = GetNonStandardJets();
1325  if(aodRecJets) nJets=aodRecJets->GetEntries();
1326  // printf("nJets = %d\n",nJets);
1327  return;
1328  }
1329 
1330  //
1331  // Background jets
1332  //
1333  AliAODJetEventBackground* aodBkgJets = 0;
1334  if(IsBackgroundJetFromReader()){//jet branch from reader
1335  AliDebug(3,"GetBackgroundJets function is called");
1336  aodBkgJets = GetBackgroundJets();
1337  AliDebug(3,Form("aodBkgJets %p",aodBkgJets));
1338  if(aodBkgJets==0x0)
1339  {
1340  if(GetDebug() > 3) event->Print();
1341  AliFatal("No jet background container found");
1342  return; // trick coverity
1343  }
1344  if(GetDebug() > 3) aodBkgJets->Print("c");
1345  }
1346 
1347  //
1348  // only background jets informations
1349  //
1350  // Float_t pTback = 0;
1351  Double_t rhoEvent=0.;
1352  if(aodBkgJets) {
1353  if(IsBackgroundJetFromReader() ) rhoEvent = aodBkgJets->GetBackground(2);
1354  if(IsHistogramJetBkg())
1355  {
1357 
1358  for(Int_t i=0;i<4;i++)
1359  {
1360  fhBkgJetBackground[i]->Fill(aodBkgJets->GetBackground(i), GetEventWeight());
1361  fhBkgJetSigma [i]->Fill(aodBkgJets->GetSigma(i) , GetEventWeight());
1362  fhBkgJetArea [i]->Fill(aodBkgJets->GetMeanarea(i) , GetEventWeight());
1363  }
1364  }//end of if fill HistogramJetBkg
1365  }//end if aodBkgJets exists
1366 
1367  //
1368  // only track information
1369  //
1370  Int_t nCTSTracks = GetCTSTracks()->GetEntriesFast();
1371  AliAODTrack *aodtrack;
1372  Int_t itrack = 0;
1373  if(IsHistogramTracks())
1374  {
1375  Double_t sumTrackPt=0;
1376  for(itrack = 0; itrack < nCTSTracks ; itrack++)
1377  {
1378  aodtrack = dynamic_cast <AliAODTrack*>(GetCTSTracks()->At(itrack));
1379 
1380  if(!aodtrack) continue;
1381 
1382  fhTrackPhiVsEta->Fill(aodtrack->Phi(), aodtrack->Eta(), GetEventWeight());
1383 
1384  sumTrackPt+=aodtrack->Pt();
1385  }
1386 
1387  if(nCTSTracks)
1388  fhTrackAveTrackPt->Fill(sumTrackPt/nCTSTracks, GetEventWeight());
1389  }//end if IsHistogramTracks
1390 
1391  //
1392  // only jet informations
1393  //
1394  AliAODJet * jettmp = 0 ;
1395  Double_t jetPttmp = 0.;
1396  Double_t var1tmp = 0.;
1397  Double_t var2tmp = 0.;
1398  // fhJetNjet->Fill(nJets, GetEventWeight());
1399  Double_t ptMostEne=0;
1400  // Int_t indexMostEne=-1;
1401  Int_t nJetsOverThreshold[10]={nJets,0,0,0,0,0,0,0,0,0};
1402  Int_t iCounter=0;
1403  Double_t sumJetTrackPt=0.;
1404  Int_t sumNJetTrack=0;
1405  Int_t nTracksInJet=0;
1406  Int_t itrk=0;
1407  for(Int_t ijet = 0; ijet < nJets ; ijet++)
1408  {
1409  jettmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ijet));
1410 
1411  if(!jettmp) continue;
1412 
1413  fhJetPtBefore->Fill(jettmp->Pt(), GetEventWeight());
1414 
1415  jetPttmp = jettmp->Pt() - rhoEvent * jettmp->EffectiveAreaCharged();//<<---changed here
1416 
1417  //calculate number of tracks above 1,2,3,4,5 GeV in jet
1418  AliVTrack* jettrack = 0x0 ;
1419 
1420  Int_t nTrackThrGeV[5]={0,0,0,0,0};
1421 
1422  nTracksInJet=(jettmp->GetRefTracks())->GetEntriesFast();
1423 
1424  fhJetNparticlesInJet->Fill(nTracksInJet, GetEventWeight());
1425 
1426  if(nTracksInJet==0) continue;
1427 
1428  sumNJetTrack+=nTracksInJet;
1429 
1430  for(itrack=0;itrack<nTracksInJet;itrack++)
1431  {
1432  jettrack=(AliVTrack *) ((jettmp->GetRefTracks())->At(itrack));
1433 
1434  if(!jettrack) continue;
1435 
1436  fhJetDeltaEtaDeltaPhi->Fill(jettmp->Eta()-jettrack->Eta(), jettmp->Phi()-jettrack->Phi(), GetEventWeight());
1437 
1438  sumJetTrackPt+=jettrack->Pt();
1439 
1440  if(IsHistogramJetTracks())
1441  {
1442  if(jettrack->Pt()>1.) nTrackThrGeV[0]++;
1443  if(jettrack->Pt()>2.) nTrackThrGeV[1]++;
1444  if(jettrack->Pt()>3.) nTrackThrGeV[2]++;
1445  if(jettrack->Pt()>4.) nTrackThrGeV[3]++;
1446  if(jettrack->Pt()>5.) nTrackThrGeV[4]++;
1447  }//end of if IsHistogramJetTracks
1448  }//end of jet track loop
1449  //printf("jetPt %f ntrks %d ntrks>1,2,3,4,5GeV in jet %d, %d, %d, %d, %d\n",jetPttmp,nTracksInJet,nTrackThrGeV[0],nTrackThrGeV[1],nTrackThrGeV[2],nTrackThrGeV[3],nTrackThrGeV[4]);
1450 
1451  for(itrack = 0; itrack < nCTSTracks ; itrack++)
1452  {
1453  aodtrack = dynamic_cast <AliAODTrack*>(GetCTSTracks()->At(itrack));
1454 
1455  if(aodtrack) fhJetDeltaEtaDeltaPhiAllTracks->Fill(jettmp->Eta()-aodtrack->Eta(),
1456  jettmp->Phi()-aodtrack->Phi(), GetEventWeight());
1457  }
1458 
1459 
1460  if(IsHistogramJetTracks()){
1461  fhJetNtracksInJetAboveThr[0]->Fill(jetPttmp, nTracksInJet, GetEventWeight());//all jets
1462 
1463  for(itrk=0;itrk<5;itrk++) {
1464  fhJetNtracksInJetAboveThr[itrk+1]->Fill(jetPttmp, nTrackThrGeV[itrk], GetEventWeight());//all jets
1465  fhJetRatioNTrkAboveToNTrk[itrk] ->Fill(jetPttmp, (Double_t)nTrackThrGeV[itrk]/(Double_t)nTracksInJet, GetEventWeight());//all jets
1466  }
1467  if(ijet==0){//most ene jet
1468  for(itrk=0;itrk<5;itrk++)
1469  fhJetNtrackRatioMostEne[itrk]->Fill(jetPttmp, (Double_t)nTrackThrGeV[itrk]/(Double_t)nTracksInJet, GetEventWeight());
1470  }
1471  if(jetPttmp>5){//jet with pt>5GeV/c
1472  for(itrk=0;itrk<5;itrk++)
1473  fhJetNtrackRatioJet5GeV[itrk]->Fill(jetPttmp, (Double_t)nTrackThrGeV[itrk]/(Double_t)nTracksInJet, GetEventWeight());
1474  }
1475  if(nTrackThrGeV[4]>0){//jet with leading particle pt>5GeV/c
1476  for(itrk=0;itrk<5;itrk++)
1477  fhJetNtrackRatioLead5GeV[itrk]->Fill(jetPttmp, (Double_t)nTrackThrGeV[itrk]/(Double_t)nTracksInJet, GetEventWeight());
1478  }
1479  }//end of if IsHistogramJetTracks
1480 
1481 
1482  Double_t effectiveChargedBgEnergy=(IsBackgroundJetFromReader()?rhoEvent * jettmp->EffectiveAreaCharged():jettmp->ChargedBgEnergy());
1483 
1484 
1485  fhJetChBkgEnergyVsArea->Fill(effectiveChargedBgEnergy, jettmp->EffectiveAreaCharged(), GetEventWeight());
1486  //if(jettmp->EffectiveAreaCharged()>0)
1487  fhJetRhoVsPt->Fill(jetPttmp, jettmp->ChargedBgEnergy()*jettmp->EffectiveAreaCharged(), GetEventWeight());
1488 
1489  if(jetPttmp>ptMostEne)
1490  {
1491  ptMostEne = jetPttmp;
1492  //indexMostEne=ijet;
1493  }
1494  if(jettmp->Pt()>=fJetMinPt)
1495  fhJetPtBeforeCut->Fill(jetPttmp, GetEventWeight());
1496 
1497  fhJetPt->Fill(jetPttmp, GetEventWeight());
1498  fhJetChBkgEnergyVsPt->Fill(jetPttmp,effectiveChargedBgEnergy, GetEventWeight());
1499  fhJetChAreaVsPt->Fill(jetPttmp,jettmp->EffectiveAreaCharged(), GetEventWeight());
1500 
1501  AliDebug(5,Form("ChargedBgEnergy %f EffectiveAreaCharged %f\n", jettmp->ChargedBgEnergy(),jettmp->EffectiveAreaCharged()));
1502 
1503  for(iCounter=1;iCounter<10;iCounter++)
1504  {
1505  if(jetPttmp>iCounter) nJetsOverThreshold[iCounter]++;
1506  }
1507 
1508  var1tmp = jettmp->Phi();
1509  var2tmp = jettmp->Eta();
1510 
1511  fhJetPhi->Fill(var1tmp, GetEventWeight());
1512  fhJetEta->Fill(var2tmp, GetEventWeight());
1513 
1514  fhJetPhiVsEta->Fill(var1tmp, var2tmp, GetEventWeight());
1515  fhJetEtaVsPt->Fill(jetPttmp, var2tmp, GetEventWeight());
1516 
1517  fhJetEtaVsNpartInJet->Fill(var2tmp, nTracksInJet, GetEventWeight());
1518 
1519  if(jetPttmp>0)
1520  fhJetEtaVsNpartInJetBkg->Fill(var2tmp, nTracksInJet, GetEventWeight());
1521 
1522  }//end of jet loop
1523 
1524  if(IsHistogramJetTracks())
1525  {
1526  if(sumNJetTrack>0)
1527  {
1528  //printf("average track pt %f\n",sumJetTrackPt/sumNJetTrack);
1529  fhJetAveTrackPt->Fill(sumJetTrackPt/sumNJetTrack, GetEventWeight());
1530  }
1531  }//end of if IsHistogramJetTracks
1532 
1533 
1534  fhJetPtMostEne->Fill(ptMostEne, GetEventWeight());
1535  for(iCounter=0;iCounter<10;iCounter++)
1536  {
1537  fhJetNjetOverPtCut[iCounter]->Fill(nJetsOverThreshold[iCounter], GetEventWeight());
1538  nJetsOverThreshold[iCounter]=0;
1539  }
1540 
1541  //end of only jet part
1542 
1543  //
1544  // Photons
1545  //
1546  Int_t ntrig = GetInputAODBranch()->GetEntriesFast() ;
1547 
1548  AliDebug(1,"Begin jet finder correlation analysis, fill histograms");
1549  AliDebug(1,Form("In particle branch aod entries %d\n", ntrig));
1550  AliDebug(1,Form("In jet output branch aod entries %d\n", event->GetNJets()));
1551 
1552  fhNjetsNgammas->Fill(nJets, ntrig, GetEventWeight());
1553 
1554  //if(nJets==0) return;//to speed up
1555  //printf("ntrig %d, njets %d\n",ntrig,nJets);
1556 
1557  //
1558  //Fill only temporary photon histograms
1559  //
1560  Double_t maxPt=0.;
1561  Int_t maxIndex=-1;
1562  Double_t tmpPt=0.;
1563  Double_t sumPt=0.;
1564  nJetsOverThreshold[0]=ntrig;
1565 
1566  for(Int_t iaod = 0; iaod < ntrig ; iaod++)
1567  {
1568  AliAODPWG4ParticleCorrelation* particlecorr = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1569  tmpPt = particlecorr->Pt();
1570  if(tmpPt>maxPt)
1571  {
1572  maxPt = tmpPt;
1573  maxIndex = iaod;
1574  }
1575 
1576  for(iCounter=1;iCounter<10;iCounter++)
1577  {
1578  if(tmpPt>iCounter) nJetsOverThreshold[iCounter]++;
1579  }
1580  sumPt+=tmpPt;
1581  }
1582 
1583  for(iCounter=0;iCounter<10;iCounter++)
1584  {
1585  fhPhotonNgammaOverPtCut[iCounter]->Fill(nJetsOverThreshold[iCounter], GetEventWeight());
1586  // nJetsOverThreshold[iCounter]=0;
1587  }
1588 
1589  fGamAvEne=0;
1590  //TObjArray* clusters = GetEMCALClusters();
1591  //printf("calculate median bkg energy for photons ");
1592  Double_t medianPhotonRho=0.;
1593  //Int_t clusterID;
1594  //Int_t iclustmp = -1;
1595  Int_t numberOfcells=0;
1596  //AliVCluster *cluster = 0;
1597  if(ntrig>1){
1598  Double_t *photonRhoArr=new Double_t[ntrig-1];
1599  fhPhotonPtMostEne->Fill(maxPt, GetEventWeight());
1600  // fhPhotonIndexMostEne->Fill(indexMaxPt, GetEventWeight());
1601  fhPhotonAverageEnergy ->Fill(sumPt/ntrig, GetEventWeight());
1602  fhPhotonRatioAveEneToMostEne->Fill(sumPt/(ntrig*maxPt), GetEventWeight());
1603  fhPhotonAverageEnergyMinus1 ->Fill((sumPt-maxPt)/(ntrig-1), GetEventWeight());
1604  fGamAvEne=(sumPt-maxPt)/(ntrig-1);
1605  fhPhotonRatioAveEneMinus1ToMostEne->Fill((sumPt-maxPt)/((ntrig-1)*maxPt), GetEventWeight());
1606 
1607  Int_t counterGamma=0;
1608  Int_t counterGammaMinus1=0;
1609 
1610  Int_t photonRhoArrayIndex=0;
1611  //TObjArray* clusterstmp = GetEMCALClusters();
1612  for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1613  AliAODPWG4ParticleCorrelation* particlecorr = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1614  if( particlecorr->Pt() > sumPt/ntrig ) counterGamma++;
1615  if( particlecorr->Pt() > (sumPt-maxPt)/(ntrig-1) ) counterGammaMinus1++;
1616 
1617  if(iaod==maxIndex) continue;
1618 // clusterID = particlecorr->GetCaloLabel(0) ;
1619 // if(clusterID < 0) continue;
1620 // cluster = FindCluster(clusters,clusterID,iclustmp);
1621 // photonRhoArr[photonRhoArrayIndex]=particlecorr->Pt()/ cluster->GetNCells();
1622 // numberOfcells+=cluster->GetNCells();
1623  photonRhoArr[photonRhoArrayIndex]=particlecorr->Pt()/ particlecorr->GetNCells();
1624  numberOfcells+=particlecorr->GetNCells();
1625  photonRhoArrayIndex++;
1626  }
1627  if(photonRhoArrayIndex>0) medianPhotonRho=TMath::Median(photonRhoArrayIndex,photonRhoArr);
1628  delete [] photonRhoArr;
1629  fhPhotonNgammaMoreAverageToNgamma ->Fill((Double_t)counterGamma / (Double_t)ntrig, GetEventWeight());
1630  fhPhotonNgammaMoreAverageMinus1ToNgamma->Fill((Double_t)counterGammaMinus1 / (Double_t)ntrig, GetEventWeight());
1631  }
1632 
1633  //printf("median = %f\n",medianPhotonRho);
1634 
1635  fhPhotonBkgRhoVsNtracks ->Fill(GetCTSTracks()->GetEntriesFast(), medianPhotonRho, GetEventWeight());
1636  fhPhotonBkgRhoVsNclusters ->Fill(ntrig, medianPhotonRho, GetEventWeight());
1637  fhPhotonBkgRhoVsCentrality->Fill(GetEventCentrality(), medianPhotonRho, GetEventWeight());
1638  fhPhotonBkgRhoVsNcells ->Fill(numberOfcells, medianPhotonRho, GetEventWeight());
1639 
1640  //AliVCluster *cluster2 = 0;
1641  Double_t photon2Corrected=0;
1642  Double_t sumPtTmp=0.;
1643  Double_t sumPtCorrectTmp=0.;
1644  AliVTrack* trackTmp = 0x0 ;
1645  TVector3 p3Tmp;
1646 
1647  for(Int_t iaod = 0; iaod < ntrig ; iaod++)
1648  {
1649  AliAODPWG4ParticleCorrelation* particlecorr = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1650 // clusterID = particlecorr->GetCaloLabel(0) ;
1651 // if(clusterID < 0) continue;
1652 // cluster = FindCluster(clusters,clusterID,iclustmp);
1653 // Int_t ncells = cluster->GetNCells();
1654  Int_t ncells = particlecorr->GetNCells();
1655  fhPhotonPt ->Fill(particlecorr->Pt(), GetEventWeight());
1656  fhPhotonPtCorrected ->Fill(particlecorr->Pt() - ncells * medianPhotonRho, GetEventWeight());
1657  fhPhotonPtDiff ->Fill(ncells * medianPhotonRho, GetEventWeight());
1658  fhPhotonPtDiffVsCentrality->Fill(GetEventCentrality(),ncells * medianPhotonRho, GetEventWeight());
1659  fhPhotonPtDiffVsNcells ->Fill(numberOfcells,ncells * medianPhotonRho, GetEventWeight());
1660  fhPhotonPtDiffVsNtracks ->Fill(GetCTSTracks()->GetEntriesFast(),ncells * medianPhotonRho, GetEventWeight());
1661  fhPhotonPtDiffVsNclusters ->Fill(ntrig,ncells * medianPhotonRho, GetEventWeight());
1662  fhPhotonPtCorrectedZoom ->Fill(particlecorr->Pt() - ncells * medianPhotonRho, GetEventWeight());
1663 
1664  //test: sum_pt in the cone 0.3 for each photon
1665  //should be: random fake gamma from MB
1666  //is: each gamma for EMCEGA
1667  sumPtTmp=0.;
1668  sumPtCorrectTmp=0.;
1669 
1670  for(Int_t iaod2 = 0; iaod2 < ntrig ; iaod2++)
1671  {
1672  if(iaod==iaod2) continue;
1673 
1674  AliAODPWG4ParticleCorrelation* particlecorr2 = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod2));
1675 // clusterID = particlecorr2->GetCaloLabel(0) ;
1676 // if(clusterID < 0) continue;
1677 // cluster2 = FindCluster(clusters,clusterID,iclustmp);
1678 // photon2Corrected = particlecorr2->Pt() - cluster2->GetNCells() * medianPhotonRho;
1679  photon2Corrected = particlecorr2->Pt() - particlecorr2->GetNCells() * medianPhotonRho;
1680 
1681  //if(Pt()<0.5) continue; //<<hardcoded here //FIXME
1682  if( TMath::Sqrt((particlecorr->Eta()-particlecorr2->Eta())*(particlecorr->Eta()-particlecorr2->Eta()) +
1683  (particlecorr->Phi()-particlecorr2->Phi())*(particlecorr->Phi()-particlecorr2->Phi()) )<fGammaConeSize ){//if(/*cone is correct*/){
1684  sumPtTmp+= particlecorr2->Pt();
1685  sumPtCorrectTmp+=photon2Corrected;
1686  }
1687  }
1688 
1689  fhPhotonSumPtInCone->Fill(sumPtTmp, GetEventWeight());
1690  fhPhotonSumPtCorrectInCone->Fill(sumPtCorrectTmp, GetEventWeight());
1691 
1692  //test: sum_pt in the cone 0.3 for each track
1693  //should be: random fake gamma from MB
1694  //is: each gamma for EMCEGA
1695  sumPtTmp=0.;
1696  for(Int_t ipr = 0;ipr < GetCTSTracks()->GetEntriesFast() ; ipr ++){
1697  trackTmp = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
1698  p3Tmp.SetXYZ(trackTmp->Px(),trackTmp->Py(),trackTmp->Pz());
1699  if( TMath::Sqrt((particlecorr->Eta()-p3Tmp.Eta())*(particlecorr->Eta()-p3Tmp.Eta()) +
1700  (particlecorr->Phi()-p3Tmp.Phi())*(particlecorr->Phi()-p3Tmp.Phi()) )<fGammaConeSize ){
1701  sumPtTmp+=p3Tmp.Pt();
1702  }
1703  }//end of loop over tracks
1704 
1705  fhPhotonSumPtChargedInCone->Fill(sumPtTmp, GetEventWeight());
1706  }
1707 
1708  //End of Fill temporary photon histograms
1709 
1710  //
1711  // Apply background subtraction for photons
1712  //
1713  fGamRho = medianPhotonRho;
1714  if(!IsBackgroundSubtractionGamma()) medianPhotonRho=0;
1715 
1716 
1717  //
1718  //Get vertex for cluster momentum calculation <<----new here
1719  //
1720  Double_t vertex[] = {0,0,0} ; //vertex ;
1721  if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
1722  GetReader()->GetVertex(vertex);
1723  fZvertex = vertex[2];
1724 
1725  //
1726  //Loop on stored AOD particles, trigger
1727  //
1728  for(Int_t iaod = 0; iaod < ntrig ; iaod++)
1729  {
1730  AliAODPWG4ParticleCorrelation* particlecorr = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1731  fhCuts ->Fill(0., GetEventWeight());
1732  fhCuts2->Fill(0., (Double_t)nJets * GetEventWeight());
1733  AliDebug(1,Form("OnlyIsolated %d !particlecorr->IsIsolated() %d \n",OnlyIsolated(), !particlecorr->IsIsolated()));
1734 
1735  if(OnlyIsolated() && !particlecorr->IsIsolated()) continue;
1736 
1737  fhCuts ->Fill(1., GetEventWeight());
1738  fhCuts2->Fill(1., nJets * GetEventWeight());
1739 
1740  if(nJets>0)
1741  {
1742  fhCuts->Fill(2., GetEventWeight());
1743  }
1744 
1745  //Recover the jet correlated, found previously.
1746  AliAODJet * jet = particlecorr->GetJet();
1747  //If correlation not made before, do it now.
1749  //Correlate with jets
1750  Int_t ijet = SelectJet(particlecorr,aodRecJets);//input for jets is TClonesArray
1751  if(ijet > -1)
1752  {
1753  AliDebug(1,Form("Jet with index %d selected \n",ijet));
1754  //jet = event->GetJet(ijet);
1755  jet = dynamic_cast<AliAODJet*>(aodRecJets-> At(ijet));
1756 
1757  if(jet) particlecorr->SetRefJet(jet);
1758 
1759  }
1760  }
1761 
1762  if (!jet) continue ;
1763 
1764  fhCuts ->Fill(3., GetEventWeight());
1765  fhCuts2->Fill(7., GetEventWeight());
1766 
1767  //check MC genereted information
1768  if(fMCStudies) FindMCgenInfo();
1769 
1770  //
1771  //Fill Correlation Histograms
1772  //
1773 // clusterID = particlecorr->GetCaloLabel(0) ;
1774 // if(!(clusterID<0)){
1775 // cluster = FindCluster(clusters,clusterID,iclustmp);
1776 // //fill tree variables
1777 // fGamNcells = cluster->GetNCells();
1778 // }
1779 
1780  fGamNcells = particlecorr->GetNCells();
1781 
1782  Double_t ptTrig = particlecorr->Pt() - medianPhotonRho * fGamNcells;//<<---changed here
1783  Double_t ptJet = jet->Pt() - rhoEvent * jet->EffectiveAreaCharged();//<<---changed here
1784  Double_t phiTrig = particlecorr->Phi();
1785  Double_t phiJet = jet->Phi();
1786  Double_t etaTrig = particlecorr->Eta();
1787  Double_t etaJet = jet->Eta();
1788  Double_t deltaPhi=phiTrig-phiJet;
1789  if(deltaPhi<0)deltaPhi+=(TMath::Pi()*2.);
1790  //printf("pT trigger %2.3f, pT jet %2.3f, Delta phi %2.3f, Delta eta %2.3f, Delta pT %2.3f, ratio %2.3f \n",
1791  // ptTrig,ptJet, phiJet-phiTrig, etaJet-etaTrig, ptTrig-ptJet, ptJet/ptTrig);
1792  fhDeltaPt ->Fill(ptTrig, ptTrig-ptJet, GetEventWeight());
1793  // fhDeltaPhi->Fill(ptTrig, phiTrig-phiJet);//need to be shifted by 2pi
1794 
1795  fhDeltaPhiCorrect->Fill(ptTrig, deltaPhi, GetEventWeight());//correct
1796 
1797  Double_t deltaPhiCorrect = TMath::Abs( particlecorr->Phi() - jet->Phi() );
1798  if ( deltaPhiCorrect > TMath::Pi() ) deltaPhiCorrect = 2. * TMath::Pi() - deltaPhiCorrect ;
1799  fhDeltaPhi0PiCorrect->Fill(ptTrig, deltaPhiCorrect, GetEventWeight());
1800 
1801  fhDeltaEta->Fill(ptTrig, etaTrig-etaJet, GetEventWeight());
1802  fhPtRatio ->Fill(ptTrig, ptJet/ptTrig , GetEventWeight());
1803  fhPt ->Fill(ptTrig, ptJet , GetEventWeight());
1804 
1805  fhSelectedJetPhiVsEta->Fill(phiJet, etaJet, GetEventWeight());
1806 
1807  fhSelectedJetChBkgEnergyVsPtJet->Fill(ptJet, (IsBackgroundJetFromReader()?rhoEvent * jet->EffectiveAreaCharged():jet->ChargedBgEnergy()), GetEventWeight());
1808 
1809  fhSelectedJetChAreaVsPtJet->Fill(ptJet, jet->EffectiveAreaCharged(), GetEventWeight());
1810  fhSelectedJetNjet ->Fill(nJets, GetEventWeight());
1811  fhSelectedNtracks ->Fill(GetCTSTracks()->GetEntriesFast(), GetEventWeight());//to be checked
1812  fhSelectedPhotonNLMVsPt ->Fill(ptTrig, particlecorr->GetNLM(), GetEventWeight());
1813 
1814 // if(clusterID < 0 ){
1815 // fhSelectedPhotonLambda0VsPt->Fill(ptTrig, -1, GetEventWeight());
1816 // //fill tree variables
1817 // fGamLambda0 = -1;
1818 // fGamTime = -1;
1819 // fGamNcells = 0;
1820 // fGamSumPtNeu=0;
1821 // }
1822 // else
1823 // {
1824  //Int_t iclus = -1;
1825  TObjArray* clusters = GetEMCALClusters();
1826  //cluster = FindCluster(clusters,clusterID,iclustmp);
1827  Double_t lambda0=particlecorr->GetM02();
1828  fhSelectedPhotonLambda0VsPt->Fill(ptTrig, lambda0, GetEventWeight());
1829  //fill tree variables
1830  fGamLambda0 = lambda0;
1831  fGamTime = particlecorr->GetTime();
1832  //fGamNcells = cluster->GetNCells();
1833 
1834  fGamSumPtNeu=0;
1835  fGamNclusters=0;
1836  //TVector3 p3Tmp;
1837  //Double_t scalarProduct=0;
1838  //Double_t vectorLength=particlecorr->P();
1839  for(Int_t icalo=0; icalo <clusters->GetEntriesFast(); icalo++){
1840  AliVCluster* calo = (AliVCluster *) clusters->At(icalo);
1841  //if(clusterID==calo->GetID()) continue;//the same cluster as trigger
1842  calo->GetMomentum(fMomentum,vertex) ;//Assume that come from vertex in straight line
1843  //printf("min pt %f\n",GetMinPt());
1844  if(fMomentum.Pt()<GetMinPt()) continue; //<<hardcoded here //FIXME 0.5 check if correct
1845  p3Tmp.SetXYZ(fMomentum.Px(),fMomentum.Py(),fMomentum.Pz());
1846  //calculate sum pt in the cone
1847  if( TMath::Sqrt((particlecorr->Eta()-p3Tmp.Eta())*(particlecorr->Eta()-p3Tmp.Eta()) +
1848  (particlecorr->Phi()-p3Tmp.Phi())*(particlecorr->Phi()-p3Tmp.Phi()) )<fGammaConeSize ){
1849  //scalarProduct = particlecorr->Px()*fMomentum.Px() + particlecorr->Py()*fMomentum.Py() + particlecorr->Pz()*fMomentum.Pz();
1850  //scalarProduct/=fMomentum.P();
1851  //scalarProduct/=vectorLength;
1852  //if(scalarProduct>TMath::Cos(0.3)) {//FIXME photon radius
1853  fGamSumPtNeu+=fMomentum.Pt();
1854  fGamNclusters++;
1855  }
1856  }
1857 // }
1858 
1859  //sum pt of charged tracks in the gamma isolation cone
1860  //starts here
1861  fGamSumPtCh=0;
1862  fGamNtracks=0;
1863  for(itrack = 0; itrack < nCTSTracks ; itrack++){
1864  aodtrack = dynamic_cast <AliAODTrack*>(GetCTSTracks()->At(itrack));
1865  if(!aodtrack) continue;
1866  fhSelectedTrackPhiVsEta->Fill(aodtrack->Phi(), aodtrack->Eta(), GetEventWeight());//fill histogram here
1867  // if(aodtrack->Pt()<0.15) continue;//hardcoded
1868  if(aodtrack->Pt()<fPtThresholdInCone) continue;
1869  if(!aodtrack->IsHybridGlobalConstrainedGlobal()) continue;
1870  if(TMath::Sqrt((particlecorr->Phi() - aodtrack->Phi())*(particlecorr->Phi() - aodtrack->Phi()) +
1871  (particlecorr->Eta() - aodtrack->Eta())*(particlecorr->Eta() - aodtrack->Eta()) ) <fGammaConeSize ) {
1872  fGamSumPtCh+=aodtrack->Pt();
1873  fGamNtracks++;
1874  }
1875  }
1876  //ends here
1877 
1878  // for(Int_t itrack = 0; itrack < nCTSTracks ; itrack++){
1879  // aodtrack = dynamic_cast <AliAODTrack*>(GetCTSTracks()->At(itrack));
1880  // fhSelectedTrackPhiVsEta->Fill(aodtrack->Phi(), aodtrack->Eta(), GetEventWeight());
1881  // }
1882 
1883  //
1884  // Background Fragmentation function
1885  //
1886  TVector3 gammaVector,jetVector;
1887  gammaVector.SetXYZ(particlecorr->Px(),particlecorr->Py(),particlecorr->Pz());
1888  jetVector.SetXYZ(jet->Px(),jet->Py(),jet->Pz());
1889  CalculateBkg(gammaVector,jetVector,vertex,1);//jet perp
1890  CalculateBkg(gammaVector,jetVector,vertex,2);//RC
1891  CalculateBkg(gammaVector,jetVector,vertex,3);//mid point
1892  CalculateBkg(gammaVector,jetVector,vertex,4);//gamma perp
1893  //CalculateBkg(gammaVector,jetVector,vertex,5);/test
1894  Double_t angleJetGam = gammaVector.Angle(jetVector);
1895  //printf("angleJetGam %f\n",angleJetGam*180/TMath::Pi());
1896 
1897  //
1898  // Fragmentation function
1899  //
1900  Float_t rad = 0, pt = 0, eta = 0, phi = 0;
1901  Int_t npartcone = 0;
1902  TVector3 p3;
1903 
1904  Int_t ntracks = 0;
1905 
1906  AliDebug(1,Form("fUseJetRefTracks %d" ,fUseJetRefTracks ));
1907  AliDebug(1,Form("jet->GetRefTracks() %p",jet->GetRefTracks()));
1908  AliDebug(1,Form("GetCTSTracks() %p" ,GetCTSTracks() ));
1909 
1910  if(!fUseJetRefTracks)
1911  ntracks =GetCTSTracks()->GetEntriesFast();
1912  else //If you want to use jet tracks from JETAN
1913  ntracks = (jet->GetRefTracks())->GetEntriesFast();
1914 
1915  AliDebug(3,Form("ntracks %d\n",ntracks));
1916  AliVTrack* track = 0x0 ;
1917  for(Int_t ipr = 0;ipr < ntracks ; ipr ++ ){
1918  if(!fUseJetRefTracks)
1919  track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
1920  else //If you want to use jet tracks from JETAN
1921  track = (AliVTrack *) ((jet->GetRefTracks())->At(ipr));
1922 
1923  p3.SetXYZ(track->Px(),track->Py(),track->Pz());
1924  pt = p3.Pt();
1925  eta = p3.Eta();
1926  phi = p3.Phi() ;
1927  if(phi < 0) phi+=TMath::TwoPi();
1928 
1929  //Check if there is any particle inside cone with pt larger than fPtThreshold
1930  rad = TMath::Sqrt((eta-etaJet)*(eta-etaJet)+ (phi-phiJet)*(phi-phiJet));
1931  if(rad < fConeSize && pt > fPtThresholdInCone)
1932  {
1933  //printf("charged in jet cone pt %f, phi %f, eta %f, R %f \n",pt,phi,eta,rad);
1934  npartcone++;
1935  fhFFz ->Fill(ptTrig, pt/ptTrig, GetEventWeight());
1936  fhFFxi->Fill(ptTrig, TMath::Log(ptTrig/pt), GetEventWeight());
1937  fhFFpt->Fill(ptTrig, pt, GetEventWeight());
1938 
1939  //according to jet axis
1940  fhJetFFz ->Fill(ptJet, pt/ptJet, GetEventWeight());
1941  fhJetFFxi->Fill(ptJet, TMath::Log(ptJet/pt), GetEventWeight());
1942  fhJetFFpt->Fill(ptJet, pt, GetEventWeight());
1943 
1944 
1945  if(TMath::Cos(angleJetGam)<0 && ptJet!=0 && pt!=0 )
1946  {
1947  fhJetFFzCor ->Fill(ptJet, -pt*TMath::Cos(angleJetGam)/ptJet, GetEventWeight());
1948  fhJetFFxiCor->Fill(ptJet, TMath::Log(ptJet/(-pt*TMath::Cos(angleJetGam))), GetEventWeight());
1949  }
1950  }
1951  }//Tracks
1952  fhNTracksInCone->Fill(ptTrig, npartcone, GetEventWeight());
1953  //fill tree here for each photon-jet (isolation required usually)
1954 
1955  fGamPt = ptTrig;
1956  //fGamLambda0 = ;//filled earlier
1957  fGamNLM = particlecorr->GetNLM();
1958  //fGamSumPtCh = ;//filled earlier
1959  //fGamTime = particlecorr->GetTOF();//filled earlier
1960  //fGamNcells = particlecorr->GetNCells();//filled earlier
1961  fGamEta = etaTrig;
1962  fGamPhi = phiTrig;
1963  //fGamSumPtNeu = ;//filled earlier
1964  //fGamNtracks = ;//filled earlier
1965  //fGamNclusters= ;//filled earlier
1966  //fGamAvEne = ;//filled earlier
1967  fJetPhi = phiJet;
1968  fJetEta = etaJet;
1969  fJetPt = ptJet;
1970  fJetBkgChEne = (IsBackgroundJetFromReader()?rhoEvent * jet->EffectiveAreaCharged():jet->ChargedBgEnergy());
1971  fJetArea = jet->EffectiveAreaCharged();
1972  fJetNtracks = (jet->GetRefTracks())->GetEntriesFast();
1973  fEventNumber = 0;
1974  fNtracks = GetCTSTracks()->GetEntriesFast();
1976  fIso = particlecorr->IsIsolated();
1977 
1978  Int_t nTrk1GeV=0;
1979  Int_t nTrk2GeV=0;
1980  for(itrack=0;itrack < fJetNtracks;itrack++){
1981  track = (AliVTrack *) ((jet->GetRefTracks())->At(itrack));
1982  if(track->Pt()>1.) nTrk1GeV++;
1983  if(track->Pt()>2.) nTrk2GeV++;
1984  }
1985 
1986  fJetNtracks1 = nTrk1GeV;
1987  fJetNtracks2 = nTrk2GeV;
1988 
1989  if(fSaveGJTree) fTreeGJ->Fill();
1990  }//AOD trigger particle loop
1991  AliDebug(1,"End fill histograms");
1992 }
1993 
1994 //__________________________________________________________________
1996 //__________________________________________________________________
1998 {
1999  if(! opt)
2000  return;
2001 
2002  printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
2004 
2005  printf("Phi trigger-jet < %3.2f\n", fDeltaPhiMaxCut) ;
2006  printf("Phi trigger-jet > %3.2f\n", fDeltaPhiMinCut) ;
2007  printf("pT Ratio trigger/jet < %3.2f\n", fRatioMaxCut) ;
2008  printf("pT Ratio trigger/jet > %3.2f\n", fRatioMinCut) ;
2009  printf("fConeSize = %3.2f\n", fConeSize) ;
2010  printf("fPtThresholdInCone = %3.2f\n", fPtThresholdInCone) ;
2011  printf("fUseJetRefTracks = %d\n", fUseJetRefTracks) ;
2012  printf("fMakeCorrelationInHistoMaker = %d\n", fMakeCorrelationInHistoMaker) ;
2013  printf("Isolated Trigger? %d\n", fSelectIsolated) ;
2014  printf("Reconstructed jet cone size = %3.2f\n", fJetConeSize) ;
2015  printf("Reconstructed jet minimum pt before background subtraction = %3.2f\n", fJetMinPt) ;
2016  printf("Reconstructed jet minimum pt after background subtraction = %3.2f\n", fJetMinPtBkgSub) ;
2017  printf("Reconstructed jet minimum area fraction = %3.2f\n", fJetAreaFraction) ;
2018 
2019  //if(!fNonStandardJetFromReader){
2020  printf("fJetBranchName = %s\n", fJetBranchName.Data()) ;
2021  //}
2023  printf("fBkgJetBranchName = %s\n", fBkgJetBranchName.Data()) ;
2024  }
2025 
2026  printf("Isolation cone size = %3.2f\n", fGammaConeSize) ;
2027  printf("fUseBackgroundSubtractionGamma = %d\n",fUseBackgroundSubtractionGamma);
2028  printf("fSaveGJTree = %d\n",fSaveGJTree);
2029  printf("fMostEnergetic = %d\n",fMostEnergetic);
2030  printf("fMostOpposite = %d\n",fMostOpposite);
2031 
2032  printf("fUseHistogramJetBkg = %d\n",fUseHistogramJetBkg);
2033  printf("fUseHistogramTracks = %d\n",fUseHistogramTracks);
2034  printf("fUseHistogramJetTracks = %d\n",fUseHistogramJetTracks);
2035  printf("fMCStudies = %d\n",fMCStudies);
2036 }
2037 
2038 //__________________________________________________________________
2044 //__________________________________________________________________
2045 void AliAnaParticleJetFinderCorrelation::CalculateBkg(TVector3 gamma, TVector3 jet,Double_t vertex[3],Int_t type=2)
2046 {
2047  //
2048  // implementation of 2 works, 1 and 4 works
2049  //
2050  Double_t gammaPt = gamma.Pt();
2051  Double_t gammaEta = gamma.Eta();
2052  Double_t gammaPhi = gamma.Phi();
2053  Double_t jetEta = jet.Eta();
2054  Double_t jetPhi = jet.Phi();
2055 
2056  //refference direction of background
2057  Double_t refEta=0.,refPhi=0.;
2058  Double_t rad = 0,rad2 = 0.;
2059  if(type==1){//perpendicular to jet axis
2060  //printf("vertex: %f %f %f \n",vertex[0],vertex[1],vertex[2]);
2061 
2062  Double_t xVar;
2063  Double_t yVar;
2064  Double_t newX=0.;
2065  Double_t newY=0.;
2066  Double_t newZ=0.;
2067  //Z axis vector: [jet.Px(),jet.Py(),jet.Pz()]
2068  Double_t jx=jet.Px();
2069  Double_t jy=jet.Py();
2070  Double_t jz=jet.Pz();
2071  //if(jz==0)printf("problem\n");
2072  //X axis
2073  Double_t Xx=1.0-vertex[0];
2074  Double_t Xy=-1.0*vertex[1];
2075  Double_t Xz=jx/jz*(vertex[0]-1.)+vertex[1]*jy/jz;
2076  //Y axis
2077  Double_t Yx=jy*Xz-jz*Xy;
2078  Double_t Yy=jz*Xx-jx*Xz;
2079  Double_t Yz=jx*Xy-jy*Xx;
2080  //Determinant
2081  Double_t det = Xx*Yy*jz + Xy*Yz*jx + Xz*Yx*jy - Xx*Yz*jy - Xy*Yx*jz - Xz*Yy*jx;
2082  if(det==0)AliWarning("problem det==0\n");
2083  Double_t detX = 0.;
2084  Double_t detY = 0.;
2085  Double_t detZ = 0.;
2086 
2087 // Double_t tmpScalar = jx*Xx+jy*Xy+jz*Xz;
2088 // printf("scalar jet.P o X: %f\n",tmpScalar);
2089 // tmpScalar = jet.Px()*Yx+jet.Py()*Yy+jet.Pz()*Yz;
2090 // printf("scalar jet.P o Y: %f\n",tmpScalar);
2091 // tmpScalar = Xx*Yx+Xy*Yy+Xz*Yz;
2092 // printf("scalar X o Y: %f\n",tmpScalar);
2093 
2094  TVector3 perp;
2095  //randomise
2096  do
2097  {
2098  refPhi=fGenerator->Rndm()*TMath::Pi()*2.;
2099  //refPhi=fGenerator->Uniform(-1,1)*TMath::Pi();
2100  xVar=TMath::Cos(refPhi);
2101  yVar=TMath::Sin(refPhi);
2102  //yVar=TMath::Cos(-TMath::Pi()/2.+refPhi);
2103  //zVar=0 in new surface frame
2104  detX = xVar*Yy*jz + Xz*yVar*jy - xVar*Yz*jy - Xy*yVar*jz;
2105  detY = Xx*yVar*jz + xVar*Yz*jx - xVar*Yx*jz - Xz*yVar*jx;
2106  detZ = Xy*yVar*jx + xVar*Yx*jy - Xx*yVar*jy - xVar*Yy*jx;
2107 
2108  newX=detX/det;
2109  newY=detY/det;
2110  newZ=detZ/det;
2111 
2112  perp.SetXYZ(newX,newY,newZ);
2113  refEta = perp.Eta();
2114  refPhi = perp.Phi();//output in <-pi, pi> range
2115  if(refPhi<0)refPhi+=2*TMath::Pi();
2116  rad = TMath::Sqrt((gammaEta-refEta)*(gammaEta-refEta) + (gammaPhi-refPhi)*(gammaPhi-refPhi));
2117  rad2 = TMath::Sqrt((jetEta-refEta)*(jetEta-refEta) + (jetPhi-refPhi)*(jetPhi-refPhi));
2118  //printf("refEta,refPhi,rad,rad2: %f %f %f %f\n",refEta,refPhi,rad,rad2);
2119  } while (rad<fJetConeSize+fGammaConeSize || rad2<2.*fJetConeSize || TMath::Abs(refEta)>0.9-fJetConeSize);
2120  fhRandomPhiEta[2]->Fill(refPhi, refEta, GetEventWeight());
2121 
2122  }
2123  else if(type==2)
2124  {//random cone
2125  //randomise
2126  do
2127  {
2128  refPhi=fGenerator->Rndm()*TMath::Pi()*2.;
2129  refEta=fGenerator->Uniform(-(0.9-fJetConeSize),0.9-fJetConeSize);
2130  rad = TMath::Sqrt((gammaEta-refEta)*(gammaEta-refEta) + (gammaPhi-refPhi)*(gammaPhi-refPhi));
2131  rad2 = TMath::Sqrt((jetEta-refEta)*(jetEta-refEta) + (jetPhi-refPhi)*(jetPhi-refPhi));
2132  //check if reference is not within jet cone or gamma cone +0.4
2133  //example: FF fConeSize=0.4, fJetConeSize=0.4, fIsoGammaConeSize=0.3
2134  } while (rad<fJetConeSize+fGammaConeSize || rad2<2.*fJetConeSize);
2135  //photon:0.7=0.4+0.3; jets:0.8=0.4 +0.4 //rad<fConeSize+fJetConeSize rad2<fConeSize+0.3
2136  fhRandomPhiEta[0]->Fill(refPhi, refEta, GetEventWeight());
2137  }
2138  else if(type==4){//perpendicular to photon axis
2139  Double_t xVar;
2140  Double_t yVar;
2141  Double_t newX=0.;
2142  Double_t newY=0.;
2143  Double_t newZ=0.;
2144  //Z axis vector: [jet.Px(),jet.Py(),jet.Pz()]
2145  Double_t jx=gamma.Px();
2146  Double_t jy=gamma.Py();
2147  Double_t jz=gamma.Pz();
2148  //if(jz==0)printf("problem\n");
2149  //X axis
2150  Double_t Xx=1.0-vertex[0];
2151  Double_t Xy=-1.0*vertex[1];
2152  Double_t Xz=jx/jz*(vertex[0]-1.)+vertex[1]*jy/jz;
2153  //Y axis
2154  Double_t Yx=jy*Xz-jz*Xy;
2155  Double_t Yy=jz*Xx-jx*Xz;
2156  Double_t Yz=jx*Xy-jy*Xx;
2157  //Determinant
2158  Double_t det = Xx*Yy*jz + Xy*Yz*jx + Xz*Yx*jy - Xx*Yz*jy - Xy*Yx*jz - Xz*Yy*jx;
2159  if(det==0)AliWarning("problem det==0");
2160  Double_t detX = 0.;
2161  Double_t detY = 0.;
2162  Double_t detZ = 0.;
2163 
2164 // Double_t tmpScalar = jx*Xx+jy*Xy+jz*Xz;
2165 // printf("scalar jet.P o X: %f\n",tmpScalar);
2166 // tmpScalar = jet.Px()*Yx+jet.Py()*Yy+jet.Pz()*Yz;
2167 // printf("scalar jet.P o Y: %f\n",tmpScalar);
2168 // tmpScalar = Xx*Yx+Xy*Yy+Xz*Yz;
2169 // printf("scalar X o Y: %f\n",tmpScalar);
2170 
2171  TVector3 perp;
2172  //randomise
2173  do{
2174  refPhi=fGenerator->Rndm()*TMath::Pi()*2.;
2175  //refPhi=fGenerator->Uniform(-1,1)*TMath::Pi();
2176  xVar=TMath::Cos(refPhi);
2177  yVar=TMath::Sin(refPhi);
2178  //yVar=TMath::Cos(-TMath::Pi()/2.+refPhi);
2179  //zVar=0 in new surface frame
2180  detX = xVar*Yy*jz + Xz*yVar*jy - xVar*Yz*jy - Xy*yVar*jz;
2181  detY = Xx*yVar*jz + xVar*Yz*jx - xVar*Yx*jz - Xz*yVar*jx;
2182  detZ = Xy*yVar*jx + xVar*Yx*jy - Xx*yVar*jy - xVar*Yy*jx;
2183 
2184  newX=detX/det;
2185  newY=detY/det;
2186  newZ=detZ/det;
2187 
2188  perp.SetXYZ(newX,newY,newZ);
2189  refEta = perp.Eta();
2190  refPhi = perp.Phi();//output in <-pi, pi> range
2191  if(refPhi<0)refPhi+=2*TMath::Pi();
2192  rad = TMath::Sqrt((gammaEta-refEta)*(gammaEta-refEta) + (gammaPhi-refPhi)*(gammaPhi-refPhi));
2193  rad2 = TMath::Sqrt((jetEta-refEta)*(jetEta-refEta) + (jetPhi-refPhi)*(jetPhi-refPhi));
2194  //printf("refEta,refPhi,rad,rad2: %f %f %f %f\n",refEta,refPhi,rad,rad2);
2195  } while (rad<fJetConeSize+fGammaConeSize || rad2<2.*fJetConeSize || TMath::Abs(refEta)>0.9-fJetConeSize);
2196  fhRandomPhiEta[1]->Fill(refPhi, refEta, GetEventWeight());
2197 
2198  }
2199  else if(type==3){//mid point
2200 
2201  Double_t jx=jet.Px();
2202  Double_t jy=jet.Py();
2203  Double_t jz=jet.Pz();
2204  // if(jz==0)printf("problem\n");
2205  Double_t gx=gamma.Px();
2206  Double_t gy=gamma.Py();
2207  Double_t gz=gamma.Pz();
2208 
2209  Double_t cosAlpha=(jx*gx+jy*gy+jz*gz)/(jet.Mag()*gamma.Mag());
2210  Double_t cosinus=TMath::Sqrt((cosAlpha+1.)/2.);
2211  //perpendicular axis
2212  Double_t Zx=gy*jz-gz*jy;
2213  Double_t Zy=gz*jx-gx*jz;
2214  Double_t Zz=gx*jy-gy*jx;
2215 
2216  //Determinant
2217  Double_t det = Zx*gy*jz + Zy*gz*jx + Zz*gx*jy - Zz*gy*jx - Zy*gx*jz - Zx*gz*jy;
2218 
2219  Double_t newX=0.;
2220  Double_t newY=0.;
2221  Double_t newZ=0.;
2222  if(det!=0) {
2223  Double_t detX = -Zy*gz*cosinus +Zz*cosinus*jy + Zz*gy*cosinus - Zy*cosinus*jz;
2224  Double_t detY = Zx*cosinus*jz - Zz*gx*cosinus - Zz*cosinus*jx + Zx*gz*cosinus;
2225  Double_t detZ = -Zx*gy*cosinus + Zy*cosinus*jx + Zy*gx*cosinus - Zx*cosinus*jy;
2226 
2227  newX=detX/det;
2228  newY=detY/det;
2229  newZ=detZ/det;
2230  }
2231 
2232  TVector3 perp;
2233  perp.SetXYZ(newX,newY,newZ);
2234  refEta = perp.Eta();
2235  refPhi = perp.Phi();//output in <-pi, pi> range
2236  if(refPhi<0)refPhi+=2*TMath::Pi();
2237  rad = TMath::Sqrt((gammaEta-refEta)*(gammaEta-refEta) + (gammaPhi-refPhi)*(gammaPhi-refPhi));
2238  rad2 = TMath::Sqrt((jetEta-refEta)*(jetEta-refEta) + (jetPhi-refPhi)*(jetPhi-refPhi));
2239  //printf("refEta,refPhi,rad,rad2: %f %f %f %f\n",refEta,refPhi,rad,rad2);
2240 
2241  if (rad<fJetConeSize+fGammaConeSize || rad2<2.*fJetConeSize || TMath::Abs(refEta)>0.9-fJetConeSize)
2242  fhRandomPhiEta[3]->Fill(refPhi, refEta, GetEventWeight());
2243  }
2244  else if(type==5){//tmp
2245  //printf("vertex: %f %f %f \n",vertex[0],vertex[1],vertex[2]);
2246 
2247  Double_t xVar;
2248  Double_t newX=0.;
2249  Double_t newY=0.;
2250  Double_t newZ=0.;
2251  //Z axis vector: [jet.Px(),jet.Py(),jet.Pz()]
2252  Double_t jx=jet.Px();
2253  Double_t jy=jet.Py();
2254  Double_t jz=jet.Pz();
2255  // if(jz==0)printf("problem\n");
2256  //X axis
2257  Double_t Xx=1.0-vertex[0];
2258  Double_t Xy=-1.0*vertex[1];
2259  Double_t Xz=jx/jz*(vertex[0]-1.)+vertex[1]*jy/jz;
2260  //Y axis
2261  Double_t Yx=jy*Xz-jz*Xy;
2262  Double_t Yy=jz*Xx-jx*Xz;
2263  Double_t Yz=jx*Xy-jy*Xx;
2264 
2265  // X and Y length
2266  Double_t Xlength=TMath::Sqrt(Xx*Xx+Xy*Xy+Xz*Xz);
2267  Double_t Ylength=TMath::Sqrt(Yx*Yx+Yy*Yy+Yz*Yz);
2268  Double_t ratio=Ylength/Xlength;
2269 
2270  TVector3 perp;
2271  //randomise
2272  do{
2273  refPhi=fGenerator->Rndm()*TMath::Pi()*2.;
2274  xVar=TMath::Tan(refPhi)/ratio;
2275  newX=xVar*Yx+Xx;
2276  newY=xVar*Yy+Xy;
2277  newZ=xVar*Yz+Xz;
2278 
2279  perp.SetXYZ(newX,newY,newZ);
2280  refEta = perp.Eta();
2281  refPhi = perp.Phi();//output in <-pi, pi> range
2282  if(refPhi<0)refPhi+=2*TMath::Pi();
2283  rad = TMath::Sqrt((gammaEta-refEta)*(gammaEta-refEta) + (gammaPhi-refPhi)*(gammaPhi-refPhi));
2284  rad2 = TMath::Sqrt((jetEta-refEta)*(jetEta-refEta) + (jetPhi-refPhi)*(jetPhi-refPhi));
2285  //printf("refEta,refPhi,rad,rad2: %f %f %f %f\n",refEta,refPhi,rad,rad2);
2286  } while (rad<fJetConeSize+fGammaConeSize || rad2<2.*fJetConeSize || TMath::Abs(refEta)>0.9-fJetConeSize);
2287  fhRandomPhiEta[4]->Fill(refPhi, refEta, GetEventWeight());
2288  }
2289 
2290  // calculate FF in background
2291  Int_t ntracks = 0;
2292  ntracks =GetCTSTracks()->GetEntriesFast();
2293  AliVTrack* track = 0x0 ;
2294  TVector3 p3;
2295 
2296  Double_t pt = 0, eta = 0, phi = 0;
2297  Int_t npartcone = 0;
2298  Double_t sumPt=0.;
2299  for(Int_t ipr = 0;ipr < ntracks ; ipr ++ ){
2300  track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
2301  p3.SetXYZ(track->Px(),track->Py(),track->Pz());
2302  pt = p3.Pt();
2303  if(pt<fPtThresholdInCone) {//0.150
2304  //printf("problem: track pt < %f MeV/c \n",fPtThresholdInCone);
2305  continue;
2306  }
2307  eta = p3.Eta() ;
2308  phi = p3.Phi() ;
2309  if(phi < 0) phi+=TMath::TwoPi();
2310  //Check if there is any particle inside cone with pt larger than fPtThreshold
2311  rad = TMath::Sqrt((eta-refEta)*(eta-refEta) + (phi-refPhi)*(phi-refPhi));
2312  if(rad < fConeSize && pt > fPtThresholdInCone){
2313  //printf("charged in jet cone pt %f, phi %f, eta %f, R %f \n",pt,phi,eta,rad);
2314  npartcone++;
2315  sumPt+=pt;
2316  if(type==1){//perp jet
2317  fhBkgFFz[1] ->Fill(gammaPt, pt/gammaPt, GetEventWeight());
2318  fhBkgFFxi[1]->Fill(gammaPt, TMath::Log(gammaPt/pt), GetEventWeight());
2319  fhBkgFFpt[1]->Fill(gammaPt, pt, GetEventWeight());
2320  }
2321  else if(type==2){//RC
2322  fhBkgFFz[0] ->Fill(gammaPt, pt/gammaPt, GetEventWeight());
2323  fhBkgFFxi[0]->Fill(gammaPt, TMath::Log(gammaPt/pt), GetEventWeight());
2324  fhBkgFFpt[0]->Fill(gammaPt, pt, GetEventWeight());
2325  }
2326  else if(type==3){//mid point
2327  fhBkgFFz[3] ->Fill(gammaPt, pt/gammaPt, GetEventWeight());
2328  fhBkgFFxi[3]->Fill(gammaPt, TMath::Log(gammaPt/pt), GetEventWeight());
2329  fhBkgFFpt[3]->Fill(gammaPt, pt, GetEventWeight());
2330  }
2331  else if(type==4){//perp jet
2332  fhBkgFFz[2] ->Fill(gammaPt, pt/gammaPt, GetEventWeight());
2333  fhBkgFFxi[2]->Fill(gammaPt, TMath::Log(gammaPt/pt), GetEventWeight());
2334  fhBkgFFpt[2]->Fill(gammaPt, pt, GetEventWeight());
2335  }
2336  else if(type==5){//test
2337  fhBkgFFz[4] ->Fill(gammaPt, pt/gammaPt, GetEventWeight());
2338  fhBkgFFxi[4]->Fill(gammaPt, TMath::Log(gammaPt/pt), GetEventWeight());
2339  fhBkgFFpt[4]->Fill(gammaPt, pt, GetEventWeight());
2340  }
2341  }
2342  }//end of loop over tracks
2343 
2344  Double_t sumOverTracks=0.;
2345  if(npartcone!=0) sumOverTracks = sumPt/npartcone;
2346  if(type==1)
2347  {
2348  fhBkgNTracksInCone [1]->Fill(gammaPt, npartcone , GetEventWeight());
2349  fhBkgSumPtInCone [1]->Fill(gammaPt, sumPt , GetEventWeight());
2350  fhBkgSumPtnTracksInCone[1]->Fill(gammaPt, sumOverTracks, GetEventWeight());
2351  }
2352  else if(type==2)
2353  {
2354  fhBkgNTracksInCone [0]->Fill(gammaPt, npartcone , GetEventWeight());
2355  fhBkgSumPtInCone [0]->Fill(gammaPt, sumPt , GetEventWeight());
2356  fhBkgSumPtnTracksInCone[0]->Fill(gammaPt, sumOverTracks, GetEventWeight());
2357  }
2358  else if(type==3)
2359  {
2360  fhBkgNTracksInCone [3]->Fill(gammaPt, npartcone , GetEventWeight());
2361  fhBkgSumPtInCone [3]->Fill(gammaPt, sumPt , GetEventWeight());
2362  fhBkgSumPtnTracksInCone[3]->Fill(gammaPt, sumOverTracks, GetEventWeight());
2363  }
2364  else if(type==4)
2365  {
2366  fhBkgNTracksInCone [2]->Fill(gammaPt, npartcone , GetEventWeight());
2367  fhBkgSumPtInCone [2]->Fill(gammaPt,sumPt , GetEventWeight());
2368  fhBkgSumPtnTracksInCone[2]->Fill(gammaPt,sumOverTracks, GetEventWeight());
2369  }
2370  else if(type==5)
2371  {
2372  fhBkgNTracksInCone [4]->Fill(gammaPt, npartcone , GetEventWeight());
2373  fhBkgSumPtInCone [4]->Fill(gammaPt, sumPt , GetEventWeight());
2374  fhBkgSumPtnTracksInCone[4]->Fill(gammaPt, sumOverTracks, GetEventWeight());
2375  }
2376 }
2377 
2378 //__________________________________________________________________
2380 //__________________________________________________________________
2382 {
2383  // frequently used variables
2384  Int_t pdg = 0 ;
2385  Int_t mother = -1 ;
2386  Int_t absID = 0 ;
2387 
2388  //Double_t photonY = -100 ;
2389  //Double_t photonE = -1 ;
2390  Double_t photonPt = -1 ;
2391  Double_t photonPhi = 100 ;
2392  Double_t photonEta = -1 ;
2393  Bool_t inacceptance = kFALSE;
2394  AliAODMCParticle * primTmp = NULL;
2395 
2396  // jet counters
2397  Int_t nParticlesInJet=0;
2398  Int_t nChargedParticlesInJet=0;
2399  Int_t nParticlesInJet150=0;
2400  Int_t nChargedParticlesInJet150=0;
2401  Int_t nChargedParticlesInJet150Cone=0;
2402 
2403  Double_t eneParticlesInJet=0.;
2404  Double_t eneChargedParticlesInJet=0.;
2405  Double_t eneParticlesInJet150=0.;
2406  Double_t eneChargedParticlesInJet150=0.;
2407  Double_t eneChargedParticlesInJet150Cone=0.;
2408 
2409  Double_t pxParticlesInJet=0.;
2410  Double_t pxChargedParticlesInJet=0.;
2411  Double_t pxParticlesInJet150=0.;
2412  Double_t pxChargedParticlesInJet150=0.;
2413  Double_t pxChargedParticlesInJet150Cone=0.;
2414 
2415  Double_t pyParticlesInJet=0.;
2416  Double_t pyChargedParticlesInJet=0.;
2417  Double_t pyParticlesInJet150=0.;
2418  Double_t pyChargedParticlesInJet150=0.;
2419  Double_t pyChargedParticlesInJet150Cone=0.;
2420 
2421  Double_t etaParticlesInJet=0.;
2422  Double_t etaChargedParticlesInJet=0.;
2423  Double_t etaParticlesInJet150=0.;
2424  Double_t etaChargedParticlesInJet150=0.;
2425  Double_t etaChargedParticlesInJet150Cone=0.;
2426 
2427  Double_t phiParticlesInJet=0.;
2428  Double_t phiChargedParticlesInJet=0.;
2429  Double_t phiParticlesInJet150=0.;
2430  Double_t phiChargedParticlesInJet150=0.;
2431  Double_t phiChargedParticlesInJet150Cone=0.;
2432 
2433  Double_t ptParticlesInJet=0.;
2434  Double_t ptChargedParticlesInJet=0.;
2435  Double_t ptParticlesInJet150=0.;
2436  Double_t ptChargedParticlesInJet150=0.;
2437  Double_t ptChargedParticlesInJet150Cone=0.;
2438 
2439  Double_t coneJet=0.;
2440  Double_t coneChargedJet=0.;
2441  Double_t coneJet150=0.;
2442  Double_t coneChargedJet150=0.;
2443 
2444  std::vector<Int_t> jetParticleIndex;
2445 
2446  if(GetReader()->ReadStack()) {//ESD
2447  AliDebug(3,"I use stack");
2448  }//end Stack
2449  else if(GetReader()->ReadAODMCParticles()) {//AOD
2450  TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
2451  if(mcparticles){
2452  //jet origin
2453  //index =6 and 7 is hard scattering (jet-quark or photon)
2454  primTmp = (AliAODMCParticle *) mcparticles->At(6);
2455  pdg=primTmp->GetPdgCode();
2456  AliDebug(3,Form("id 6 pdg %d, pt %f ",pdg,primTmp->Pt() ));
2457  if(TMath::Abs(pdg)<=6 ||pdg==21) {
2458  fhMCJetOrigin->Fill(pdg, GetEventWeight());
2460  }
2461  primTmp = (AliAODMCParticle *) mcparticles->At(7);
2462  pdg=primTmp->GetPdgCode();
2463  AliDebug(3,Form("id 7 pdg %d, pt %f",pdg,primTmp->Pt() ));
2464  if(TMath::Abs(pdg)<=6 ||pdg==21) {
2465  fhMCJetOrigin->Fill(pdg, GetEventWeight());
2467  }
2468  //end of jet origin
2469 
2470  Int_t nprim = mcparticles->GetEntriesFast();
2471  for(Int_t i=0; i < nprim; i++) {
2472 
2473  if ( !GetReader()->AcceptParticleMCLabel( i ) ) continue ;
2474 
2475  AliAODMCParticle * prim = (AliAODMCParticle *) mcparticles->At(i);
2476  pdg = prim->GetPdgCode();
2477  mother=prim->GetMother();
2478  //photon=22, gluon=21, quarks=(1,...,6), antiquarks=(-1,...,-6)
2479  if(pdg == 22){//photon
2480  fhMCPhotonCuts->Fill(0., GetEventWeight());
2481 
2482  if(prim->GetStatus()!=1) continue;
2483 
2484  fhMCPhotonCuts->Fill(1., GetEventWeight());
2485 
2486  AliDebug(5,Form("id %d, prim %d, physPrim %d, status %d\n",i,prim->IsPrimary(),prim->IsPhysicalPrimary(),prim->GetStatus()));
2487  while(mother>7)
2488  {
2489  primTmp = (AliAODMCParticle *) mcparticles->At(mother);
2490  mother=primTmp->GetMother();
2491  }
2492 
2493  if(mother<6)continue;
2494 
2495  fhMCPhotonCuts->Fill(2., GetEventWeight());
2496 
2497  primTmp = (AliAODMCParticle *) mcparticles->At(mother);
2498 
2499  if(primTmp->GetPdgCode()!=22)continue;
2500 
2501  fhMCPhotonCuts->Fill(3., GetEventWeight());
2502 
2503  //Get photon kinematics
2504  photonPt = prim->Pt() ;
2505  photonPhi = prim->Phi() ;
2506  if(photonPhi < 0) photonPhi+=TMath::TwoPi();
2507  photonEta = prim->Eta() ;
2508 
2509  fhMCPhotonPt ->Fill(photonPt, GetEventWeight());
2510  fhMCPhotonEtaPhi->Fill(photonPhi, photonEta, GetEventWeight());
2511 
2512  //Check if photons hit the Calorimeter
2513  fMomentum.SetPxPyPzE(prim->Px(),prim->Py(),prim->Pz(),prim->E());
2514  inacceptance = kFALSE;
2516  {
2517  fhMCPhotonCuts->Fill(4, GetEventWeight());
2518 
2519  //check if in EMCAL
2520  if(GetEMCALGeometry())
2521  {
2522  GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(prim->Eta(),prim->Phi(),absID);
2523  if(absID >= 0) inacceptance = kTRUE;
2524  AliDebug(3,Form("In EMCAL Real acceptance? %d",inacceptance));
2525  }
2526  else{
2527  if(GetFiducialCut()->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kEMCAL)) inacceptance = kTRUE ;
2528  AliDebug(1,Form("In EMCAL fiducial cut acceptance? %d",inacceptance));
2529  }
2530  }
2531  else
2532  {//no EMCAL nor EMCALGeoMatrixSet
2533  AliWarning("not EMCALGeoMatrix set");
2534  }//end of check if EMCAL
2535 
2536  if(inacceptance)fhMCPhotonCuts->Fill(5, GetEventWeight());
2537 
2538  AliDebug(5,Form("Photon Energy %f, Pt %f",prim->E(),prim->Pt()));
2539  fMCGamPt=photonPt;
2540  fMCGamEta=photonEta;
2541  fMCGamPhi=photonPhi;
2542  }//end of check if photon
2543  else
2544  {//not photon
2545  if(prim->GetStatus()!=1) continue;
2546 
2547  AliDebug(5,Form("id %d, prim %d, physPrim %d, status %d, pdg %d, E %f",
2548  i,prim->IsPrimary(),prim->IsPhysicalPrimary(),prim->GetStatus(),prim->GetPdgCode(),prim->E()));
2549 
2550  while(mother>7)
2551  {
2552  primTmp = (AliAODMCParticle *) mcparticles->At(mother);
2553  mother=primTmp->GetMother();
2554  AliDebug(5,Form("next mother %d",mother));
2555  }
2556  if(mother<6)continue;//soft part
2557 
2558  primTmp = (AliAODMCParticle *) mcparticles->At(mother);
2559  pdg=primTmp->GetPdgCode();
2560  if( !(TMath::Abs(pdg)<=6 || pdg==21) ) continue;//origin not hard q or g
2561 
2562  //jetParticleIndex.Add(&i);
2563  jetParticleIndex.push_back(i);
2564 
2565  nParticlesInJet++;
2566  eneParticlesInJet+=prim->E();
2567  pxParticlesInJet+=prim->Px();
2568  pyParticlesInJet+=prim->Py();
2569  etaParticlesInJet+=(prim->E()*prim->Eta());
2570  photonPhi = prim->Phi() ;
2571  if(photonPhi < 0) photonPhi+=TMath::TwoPi();
2572  phiParticlesInJet+=(prim->E()*photonPhi);
2573 
2574  if(prim->Charge()!=0) {
2575  nChargedParticlesInJet++;
2576  eneChargedParticlesInJet+=prim->E();
2577  pxChargedParticlesInJet+=prim->Px();
2578  pyChargedParticlesInJet+=prim->Py();
2579  etaChargedParticlesInJet+=(prim->E()*prim->Eta());
2580  phiChargedParticlesInJet+=(prim->E()*photonPhi);
2581  }
2582  if(prim->Pt()>0.150 && TMath::Abs(prim->Eta())<0.9 ) {//in TPC acceptance :)
2583  nParticlesInJet150++;
2584  eneParticlesInJet150+=prim->E();
2585  pxParticlesInJet150+=prim->Px();
2586  pyParticlesInJet150+=prim->Py();
2587  etaParticlesInJet150+=(prim->E()*prim->Eta());
2588  phiParticlesInJet150+=(prim->E()*photonPhi);
2589  }
2590  if(prim->Charge()!=0 && prim->Pt()>0.150 && TMath::Abs(prim->Eta())<0.9 ) {//in TPC acceptance :)
2591  nChargedParticlesInJet150++;
2592  eneChargedParticlesInJet150+=prim->E();
2593  pxChargedParticlesInJet150+=prim->Px();
2594  pyChargedParticlesInJet150+=prim->Py();
2595  etaChargedParticlesInJet150+=(prim->E()*prim->Eta());
2596  phiChargedParticlesInJet150+=(prim->E()*photonPhi);
2597  }
2598 
2599  }//end of check pdg
2600  }//end of loop over primaries
2601 
2602  if(eneParticlesInJet != 0.) {
2603  etaParticlesInJet/=eneParticlesInJet ;
2604  phiParticlesInJet/=eneParticlesInJet ;
2605  }
2606  if(eneChargedParticlesInJet != 0) {
2607  etaChargedParticlesInJet/=eneChargedParticlesInJet;
2608  phiChargedParticlesInJet/=eneChargedParticlesInJet;
2609  }
2610  if(eneParticlesInJet150 != 0) {
2611  etaParticlesInJet150/=eneParticlesInJet150;
2612  phiParticlesInJet150/=eneParticlesInJet150;
2613  }
2614  if(eneChargedParticlesInJet150 != 0) {
2615  etaChargedParticlesInJet150/=eneChargedParticlesInJet150;
2616  phiChargedParticlesInJet150/=eneChargedParticlesInJet150;
2617  }
2618 
2619  ptParticlesInJet=TMath::Sqrt(pxParticlesInJet*pxParticlesInJet+pyParticlesInJet*pyParticlesInJet);
2620  ptChargedParticlesInJet=TMath::Sqrt(pxChargedParticlesInJet*pxChargedParticlesInJet+pyChargedParticlesInJet*pyChargedParticlesInJet);
2621  ptParticlesInJet150=TMath::Sqrt(pxParticlesInJet150*pxParticlesInJet150+pyParticlesInJet150*pyParticlesInJet150);
2622  ptChargedParticlesInJet150=TMath::Sqrt(pxChargedParticlesInJet150*pxChargedParticlesInJet150+pyChargedParticlesInJet150*pyChargedParticlesInJet150);
2623 
2624  Double_t distance=0.;
2625  Double_t eta=0.;
2626  Double_t phi=0.;
2627  Double_t mostPtCharged=0.;
2628  Int_t mostmostPtChargedId=-1;
2629  std::vector<Int_t>::iterator it;
2630  for( it=jetParticleIndex.begin(); it!=jetParticleIndex.end(); ++it ){
2631  AliAODMCParticle * prim = (AliAODMCParticle *) mcparticles->At(*it);
2632  eta = prim->Eta();
2633  phi = prim->Phi();
2634  if(phi < 0) phi+=TMath::TwoPi();
2635  //full jet
2636  distance=TMath::Sqrt((eta-etaParticlesInJet)*(eta-etaParticlesInJet)+(phi-phiParticlesInJet)*(phi-phiParticlesInJet));
2637  if(distance>coneJet) coneJet=distance;
2638  //charged jet
2639  distance=TMath::Sqrt((eta-etaChargedParticlesInJet)*(eta-etaChargedParticlesInJet)+(phi-phiChargedParticlesInJet)*(phi-phiChargedParticlesInJet));
2640  if(distance>coneChargedJet) coneChargedJet=distance;
2641  //
2642  distance=TMath::Sqrt((eta-etaParticlesInJet150)*(eta-etaParticlesInJet150)+(phi-phiParticlesInJet150)*(phi-phiParticlesInJet150));
2643  if(distance>coneJet150 && TMath::Abs(eta)<0.9 ) coneJet150=distance;
2644  //
2645  distance=TMath::Sqrt((eta-etaChargedParticlesInJet150)*(eta-etaChargedParticlesInJet150)+(phi-phiChargedParticlesInJet150)*(phi-phiChargedParticlesInJet150));
2646  if(distance>coneChargedJet150 && TMath::Abs(eta)<0.9) coneChargedJet150=distance;
2647 
2648  if(prim->Charge()!=0 && prim->Pt()>0.150 && TMath::Abs(eta)<0.9) {
2649  if(prim->Pt()>mostPtCharged) {
2650  mostPtCharged=prim->Pt();
2651  mostmostPtChargedId=(*it);
2652  }
2653  }
2654 
2655  if(distance<=0.4){
2656  if(prim->Charge()!=0 && prim->Pt()>0.150 && TMath::Abs(eta)<0.9) {
2657  nChargedParticlesInJet150Cone++;
2658  eneChargedParticlesInJet150Cone+=prim->E();
2659  pxChargedParticlesInJet150Cone+=prim->Px();
2660  pyChargedParticlesInJet150Cone+=prim->Py();
2661  etaChargedParticlesInJet150Cone+=(prim->E()*eta);
2662  phiChargedParticlesInJet150Cone+=(prim->E()*phi);
2663  }
2664 
2665  }
2666 
2667  }//end of loop over jet particle indexes
2668  if(eneChargedParticlesInJet150Cone != 0) {
2669  etaChargedParticlesInJet150Cone/=eneChargedParticlesInJet150Cone;
2670  phiChargedParticlesInJet150Cone/=eneChargedParticlesInJet150Cone;
2671  }
2672  ptChargedParticlesInJet150Cone=TMath::Sqrt(pxChargedParticlesInJet150Cone*pxChargedParticlesInJet150Cone+pyChargedParticlesInJet150Cone*pyChargedParticlesInJet150Cone);
2673  if(nChargedParticlesInJet150>0 && nChargedParticlesInJet150Cone<1){//no particles in cone? take the most energetic one
2674  nChargedParticlesInJet150Cone=1;
2675  etaChargedParticlesInJet150Cone=((AliAODMCParticle *)mcparticles->At(mostmostPtChargedId))->Eta();
2676  phiChargedParticlesInJet150Cone=((AliAODMCParticle *)mcparticles->At(mostmostPtChargedId))->Phi();
2677  ptChargedParticlesInJet150Cone=((AliAODMCParticle *)mcparticles->At(mostmostPtChargedId))->Pt();
2678  }
2679 
2680 
2681  }//mc array exists and data is MC
2682  }// read AOD MC
2683 
2684  jetParticleIndex.clear();
2685 
2686  //printouts
2687 
2688  AliDebug(3,Form("cone full %f, charged %f, full150 %f, charged150 %f",coneJet,coneChargedJet,coneJet150,coneChargedJet150));
2689  AliDebug(3,Form("Npart %d, NchPart %d, Npart(pt>150M) %d, NchPart(pt>150M) %d, NchPart(pt>150M)Cone %d\n",nParticlesInJet,nChargedParticlesInJet,nParticlesInJet150,nChargedParticlesInJet150,nChargedParticlesInJet150Cone));
2690  AliDebug(3,Form("Etot %f, Ech %f, E(pt>150M) %f, Ech(pt>150M) %f\n",eneParticlesInJet,eneChargedParticlesInJet,eneParticlesInJet150,eneChargedParticlesInJet150));
2691  AliDebug(3,Form("pt %f, ptch %f, pt(pt>150M) %f,ptch(pt>150M) %f,ptch(pt>150M)Cone %f\n",ptParticlesInJet,ptChargedParticlesInJet,ptParticlesInJet150,ptChargedParticlesInJet150,ptChargedParticlesInJet150Cone));
2692  AliDebug(3,Form("eta/phi tot %f/%f, ch %f/%f, tot150 %f/%f, ch150 %f/%f, ch150cone %f/%f\n",etaParticlesInJet,phiParticlesInJet,etaChargedParticlesInJet,phiChargedParticlesInJet,etaParticlesInJet150,phiParticlesInJet150,etaChargedParticlesInJet150,phiChargedParticlesInJet150,etaChargedParticlesInJet150Cone,phiChargedParticlesInJet150Cone));
2693 
2694  //fill histograms
2695  if(ptParticlesInJet) fhMCJetRatioChFull->Fill(ptChargedParticlesInJet/ptParticlesInJet, GetEventWeight());
2696  if(ptChargedParticlesInJet) fhMCJetRatioCh150Ch->Fill(ptChargedParticlesInJet150/ptChargedParticlesInJet, GetEventWeight());
2697 
2698  fhMCJetNPartVsPt ->Fill(ptParticlesInJet,nParticlesInJet, GetEventWeight());
2699  fhMCJetChNPartVsPt ->Fill(ptChargedParticlesInJet,nChargedParticlesInJet, GetEventWeight());
2700  fhMCJetNPart150VsPt ->Fill(ptParticlesInJet150,nParticlesInJet150, GetEventWeight());
2701  fhMCJetChNPart150VsPt->Fill(ptChargedParticlesInJet150,nChargedParticlesInJet150, GetEventWeight());
2702  fhMCJetChNPart150ConeVsPt->Fill(ptChargedParticlesInJet150Cone,nChargedParticlesInJet150Cone, GetEventWeight());
2703 
2704  fhMCJetEtaPhi->Fill(phiParticlesInJet,etaParticlesInJet, GetEventWeight());
2705  fhMCJetChEtaPhi->Fill(phiChargedParticlesInJet,etaChargedParticlesInJet, GetEventWeight());
2706  fhMCJet150EtaPhi->Fill(phiParticlesInJet150,etaParticlesInJet150, GetEventWeight());
2707  fhMCJetCh150EtaPhi->Fill(phiChargedParticlesInJet150,etaChargedParticlesInJet150, GetEventWeight());
2708  fhMCJetCh150ConeEtaPhi->Fill(phiChargedParticlesInJet150Cone,etaChargedParticlesInJet150Cone, GetEventWeight());
2709 
2710  //fill tree
2711  fMCJetPt = ptParticlesInJet;
2712  fMCJetChPt = ptChargedParticlesInJet;
2713  fMCJet150Pt = ptParticlesInJet150;
2714  fMCJetCh150Pt = ptChargedParticlesInJet150;
2715  fMCJetNPart = nParticlesInJet;
2716  fMCJetChNPart = nChargedParticlesInJet;
2717  fMCJet150NPart = nParticlesInJet150;
2718  fMCJetCh150NPart = nChargedParticlesInJet150;
2719  fMCJetEta = etaParticlesInJet ;
2720  fMCJetPhi = phiParticlesInJet ;
2721  fMCJetChEta = etaChargedParticlesInJet ;
2722  fMCJetChPhi = phiChargedParticlesInJet ;
2723  fMCJet150Eta = etaParticlesInJet150 ;
2724  fMCJet150Phi = phiParticlesInJet150 ;
2725  fMCJetCh150Eta = etaChargedParticlesInJet150;
2726  fMCJetCh150Phi = phiChargedParticlesInJet150;
2727 
2728  fMCJetCh150ConePt = ptChargedParticlesInJet150Cone;
2729  fMCJetCh150ConeNPart = nChargedParticlesInJet150Cone;
2730  fMCJetCh150ConeEta = etaChargedParticlesInJet150Cone;
2731  fMCJetCh150ConePhi = phiChargedParticlesInJet150Cone;
2732 
2733 } // end of method FindMCgenInfo
2734 
2735 
Float_t GetHistoPtMax() const
Int_t fMCJetNPart
MC gen number of full jet particles.
TH2F * fhJetNtrackRatioMostEne[5]
! the same for most energetic jet
Int_t pdg
Int_t fGamNtracks
number of tracks in iso cone
Bool_t fSaveGJTree
flag to save gamma-jet tree
TH2F * fhNjetsNgammas
! Number of jets vs number of photons in the event
Float_t GetHistoPtMin() const
virtual TObjArray * GetCTSTracks() const
Bool_t fMakeCorrelationInHistoMaker
Make particle-jet correlation in histogram maker.
double Double_t
Definition: External.C:58
virtual void FillInputNonStandardJets()
Bool_t fSelectIsolated
Select only trigger particles isolated.
TH2F * fhMCJetChNPart150VsPt
! generated N parts (pt>150 MeV/c) vs pt charged jet
Double_t fMCJet150Eta
MC gen full jet eta (pt>150MeV/c)
virtual void AddToHistogramsName(TString add)
Int_t fMCJet150NPart
MC gen number of full jet particles (pt>150MeV/c)
Definition: External.C:236
TH2F * fhPhotonPtDiffVsNtracks
! correction vs Ntracks
Bool_t fMostEnergetic
flag to choose gamma-jet pairs most energetic
Double_t fJetConeSize
Reconstructed jet cone size.
Double_t fGamRho
background energy for photons per cell in EMCal
void InitParameters()
Initialize the parameters of the analysis.
TH2F * fhMCJetEtaPhi
! generated jet eta vs phi for full jet
TH2F * fhDeltaPhiBefore
! Difference of jet phi and trigger particle phi as function of trigger particle pT ...
TH1F * fhJetPtMostEne
! Pt of the most energetic jet
TH1F * fhBkgJetSigma[4]
! sigma of jet in backgroud branch
TH2F * fhSelectedPhotonNLMVsPt
! nlm vs pt for selected photons
virtual AliAODJetEventBackground * GetBackgroundJets() const
Double_t fMCJetCh150Eta
MC gen charged jet eta (pt>150MeV/c)
TH2F * fhRandomPhiEta[5]
! eta and phi from random generator
TH2F * fhSelectedTrackPhiVsEta
! Phi vs eta of all chosen tracks in selected events
virtual AliVEvent * GetInputEvent() const
virtual void SetInputAODName(TString name)
TH1F * fhPhotonAverageEnergyMinus1
! average energy of photon w/o most ene photon
TH2F * fhBkgFFz[5]
! Background fragmentation function, z=ptjet/pttrig
TH2F * fhJetChAreaVsPt
! area of each charged jet vs jet pt
Int_t fJetNtracks2
number of jet tracks with pt>2 GeV/c
TH2F * fhPtRatio
! Ratio of jet pT and trigger particle pT as function of trigger particle pT
virtual TClonesArray * GetNonStandardJets() const
void FindMCgenInfo()
Find information about photon and (quark or gluon) on generated level.
Double_t fMCJet150Phi
MC gen full jet phi (pt>150MeV/c)
TH2F * fhSelectedPhotonLambda0VsPt
! lambda0 vs pt for selected photons
TH2F * fhDeltaPhi0PiCorrect
! Difference of jet phi and trigger particle phi as function of trigger particle pT ...
TH1F * fhJetPt
! Pt of all jets after bkg correction
Int_t fMCJetCh150NPart
MC gen number of charged jet particles (pt>150MeV/c)
TH1F * fhPhotonSumPtInCone
! sum pt in cone before correction
TH2F * fhFFz
! Accepted reconstructed jet fragmentation function, z=pt^particle,jet/pttrig
TH2F * fhJetChBkgEnergyVsPt
! background energy of each charged jet vs jet pt
TH1F * fhPhotonSumPtChargedInCone
! sum pt of charged tracks in the cone before correction
Int_t fMCJetCh150ConeNPart
MC gen number of charged jet particles (pt>150MeV/c),R=0.4.
Double_t fGamSumPtNeu
energy in isolation cone neutral
TH2F * fhMCJetChEtaPhi
! generated jet eta vs phi for charged jet
TH1F * fhPhotonPt
! pt of gamma before bkg correction
TH1F * fhJetNjetOverPtCut[10]
! number of reconstructed jets in event over pT threshold
Int_t fMCPartonType
MC gen parton type origin of jet.
virtual void GetVertex(Double_t v[3]) const
TH2F * fhBkgSumPtInCone[5]
! Background sum pt in cone
TH2F * fhPtRatioBefore
! Ratio of jet pT and trigger particle pT as function of trigger particle pT
TH2F * fhMCJetNPartVsPt
! generated N parts vs pt full jet
TH2F * fhMCJetNPart150VsPt
! generated N parts (pt>150 MeV/c) vs pt full jet
TH2F * fhJetNtrackRatioJet5GeV[5]
! the same for pt jet above 5 GeV
TH2F * fhJetRatioNTrkAboveToNTrk[5]
! ratio tracks in jet with pt above 1,2,3,4,5GeV to ntracks
TH2F * fhJetNtrackRatioLead5GeV[5]
! the same for jet with leading particle pt>5GeV
TH1F * fhJetPtBeforeCut
! Pt of all jets after bkg correction, raw jet pt>fJetMinPt
TString fJetBranchName
name of jet branch not set in reader part //new
TH2F * fhDeltaPt
! Difference of jet pT and trigger particle pT as function of trigger particle pT ...
TRandom2 * fGenerator
! pointer to random generator object
Base class for CaloTrackCorr analysis algorithms.
Double_t fMCJet150Pt
MC gen full jet (pt^particles>150MeV/c) pt.
virtual AliFiducialCut * GetFiducialCut()
AliAnaParticleJetFinderCorrelation()
Default constructor. Initialize parameters.
TH2F * fhDeltaPhi0PiCorrectBefore
! Difference of jet phi and trigger particle phi (0,pi) as function of trigger particle pT ...
TH1F * fhMCJetRatioChFull
! generated ratio pt charged/full jet
int Int_t
Definition: External.C:63
virtual TClonesArray * GetInputAODBranch() const
Bool_t fUseHistogramJetTracks
flag to save jet tracks features
void MakeAnalysisFillHistograms()
Particle-Jet Correlation Analysis, fill histograms.
virtual TClonesArray * GetAODMCParticles() const
virtual AliHistogramRanges * GetHistogramRanges()
TH2F * fhPhotonPtDiffVsCentrality
! correction vs centrality
TH2F * fhJetDeltaEtaDeltaPhi
! delta eta vs delta phi for (jet-track) <-0.8,0.8>
float Float_t
Definition: External.C:68
TH1F * fhGamPtPerTrig
! per trigger normalisation
TH2F * fhJetEtaVsNpartInJetBkg
! Eta vs number of particles in jet for background subtracted jets
void CalculateBkg(TVector3 gamma, TVector3 jet, Double_t *vector, Int_t type)
const Double_t ptmax
virtual AliEMCALGeometry * GetEMCALGeometry() const
TH2F * fhBkgSumPtnTracksInCone[5]
! Background sum pt over ntracks in cone
TH2F * fhSelectedJetPhiVsEta
! phi vs eta of selected jet
TH2F * fhDeltaEta
! Difference of jet eta and trigger particle eta as function of trigger particle pT ...
virtual AliAODEvent * GetOutputEvent() const
TH1F * fhMCJetRatioCh150Ch
! generated ratio pt charged(pt>150MeV/c)/charged jet
Double_t fMCJetCh150Pt
MC gen charged jet (pt^particles>150MeV/c) pt.
Double_t fRatioMaxCut
Jet/particle Ratio cut maximum.
TH1F * fhCuts
! Number of events after cuts
Bool_t fUseHistogramJetBkg
flag to save bkg jet histograms
virtual AliCalorimeterUtils * GetCaloUtils() const
TH1F * fhPhotonNgammaMoreAverageMinus1ToNgamma
! number of gammas with ene. more than average ene (w/o most ene gamma) divided by no...
TH1F * fhPhotonPtCorrectedZoom
! pt of gamma after background correction in +-5 GeV/c
TH2F * fhJetFFz
! Accepted reconstructed jet fragmentation function, z=pt^particle,jet/ptjet
const Double_t ptmin
TH2F * fhFFpt
! Jet particle pt distribution in cone
TH2F * fhJetChBkgEnergyVsArea
! area of each charged jet vs jet background
TH1F * fhSelectedJetNjet
! number of jets in selected event
TH1F * fhBkgJetArea[4]
! area of jet in bkg branch
TH2F * fhMCJet150EtaPhi
! generated jet eta vs phi full jet (pt>150 MeV/c)
virtual Double_t GetEventWeight() const
Double_t fJetMinPt
Minumum jet pt, default 5GeV/c.
Bool_t fMostOpposite
flag to choose gamma-jet pairs most opposite
Bool_t fBackgroundJetFromReader
use background jet from reader //new
TH1F * fhSelectedNtracks
! number of tracks in selected event
TH2F * fhDeltaPtBefore
! Difference of jet pT and trigger particle pT as function of trigger particle pT ...
Double_t fGamAvEne
average energy of photons (without most ene)
TH2F * fhPhotonBkgRhoVsNclusters
! average energy in one cell vs n clusters
Double_t fMCJetCh150ConePhi
MC gen charged jet phi (pt>150MeV/c),R=0.4.
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
TH2F * fhMCJetCh150EtaPhi
! generated jet eta vs phi charged jet (pt>150 MeV/c)
Double_t fConeSize
Jet cone size to calculate fragmentation function.
TH2F * fhFFxi
! Accepted reconstructed jet fragmentation function, xsi = ln(pttrig/pt^particle,jet) ...
TH2F * fhJetNtracksInJetAboveThr[6]
! number of tracks in jet with pt above 0,1,2,3,4,5GeV
TH2F * fhSelectedJetChBkgEnergyVsPtJet
! background energy of selected charged jet vs jet pt
TH1F * fhTrackAveTrackPt
! average track pt in event
Double_t fMCJetCh150Phi
MC gen charged jet phi (pt>150MeV/c)
Int_t fGamNclusters
number of clusters in iso cone
TH2F * fhJetFFzCor
! Accepted reconstructed jet fragmentation function, z=pt^particle,jet*-cos(jet,trig)/ptjet ...
TH1F * fhPhotonPtDiff
! bkg correction = n_cells * median_rho
void Print(const Option_t *opt) const
Print some relevant parameters set for the analysis.
TH2F * fhJetDeltaEtaDeltaPhiAllTracks
! delta eta vs delta phi for (jet-track) <-pi,pi>
TH1F * fhBkgJetBackground[4]
! background from jet bkg branch
TH2F * fhNTracksInCone
! jet multiplicity in cone
TH2F * fhPt
! jet pT vs trigger particle pT
Int_t GetHistoPtBins() const
TH2F * fhPhotonBkgRhoVsNtracks
! average energy in one cell vs n tracks
TH2F * fhPhotonPtDiffVsNclusters
! correction vs Nclustres
TH2F * fhJetFFxi
! Accepted reconstructed jet fragmentation function, xsi = ln(ptjet/pt^particle,jet) ...
TH2F * fhMCPhotonEtaPhi
! generated direct photon eta vs phi
TH1F * fhMCPhotonPt
! generated direct photon pt
TH2F * fhSelectedJetChAreaVsPtJet
! area of selected charged jet vs jet pt
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
TH2F * fhPtGamPtJet
! gamma jet correlation filling
Bool_t fUseJetRefTracks
Use track references from JETAN not the AOD tracks to calculate fragmentation function.
Int_t fJetNtracks1
number of jet tracks with pt>1 GeV/c
Bool_t fUseHistogramTracks
flag to save CTS tracks features
Double_t fDeltaPhiMaxCut
Minimum Delta Phi Gamma-Leading.
TH2F * fhJetFFxiCor
! Accepted reconstructed jet fragmentation function, xsi = ln(ptjet/pt^particle*-cos(jet,trig),jet)
TH1F * fhPhotonAverageEnergy
! average energy of photon
TH2F * fhPtBefore
! jet pT vs trigger particle pT
TH2F * fhDeltaPhiCorrect
! Difference of jet phi and trigger particle phi as function of trigger particle pT ...
virtual void Print(const Option_t *) const
Print some relevant parameters set for the analysis.
TH1F * fhJetAveTrackPt
! average track from jets pt in event
TH1F * fhJetNparticlesInJet
! number of particles in jets
TH2F * fhBkgFFxi[5]
! Background fragmentation function, xsi = ln(pttrig/ptjet)
TString fBkgJetBranchName
name of background jet branch not set in reader part //new
TH2F * fhBkgFFpt[5]
! Background particle pt distribution in cone
Double_t fJetMinPtBkgSub
Minumum jet pt after bkg subtraction, default -100 GeV/c.
TH2F * fhBkgNTracksInCone[5]
! Background multiplicity in cone
const char Option_t
Definition: External.C:48
Double_t fMCJetCh150ConeEta
MC gen charged jet eta (pt>150MeV/c),R=0.4.
Double_t fGamSumPtCh
energy in isolation cone charged
virtual Int_t GetDataType() const
TH1F * fhPhotonSumPtCorrectInCone
! sum pt in cone afrer correction
Double_t fRatioMinCut
Jet/particle Ratio cut minimum.
Bool_t fUseBackgroundSubtractionGamma
flag to use backgrouind subtraction for photons or not
Double_t fPtThresholdInCone
Jet pT threshold in jet cone.
TH1F * fhPhotonNgammaOverPtCut[10]
! number of photons in event over pT threshold
virtual AliCaloTrackReader * GetReader() const
bool Bool_t
Definition: External.C:53
Bool_t IsEMCALGeoMatrixSet() const
TH2F * fhJetRhoVsCentrality
! jet energy density vs centrality
TH2F * fhPhotonBkgRhoVsCentrality
! average energy in one cell vs centrality
Int_t SelectJet(AliAODPWG4Particle *particle, TClonesArray *aodRecJets)
TH1F * fhPhotonPtCorrected
! pt of gamma after background correction
TH1F * fhPhotonRatioAveEneMinus1ToMostEne
! ratio average energy of photon w/o most ene photon to most energetic photon
TH2F * fhJetFFpt
! Jet particle pt distribution in jet cone
virtual TObjArray * GetEMCALClusters() const
Double_t fDeltaPhiMinCut
Maximum Delta Phi Gamma-Leading.
Int_t fMCJetChNPart
MC gen number of charged jet particles.
TH2F * fhTrackPhiVsEta
! Phi vs eta of all chosen tracks in all events
TH2F * fhJetRhoVsPt
! jet energy density vs jet pt
TH2F * fhDeltaEtaBefore
! Difference of jet eta and trigger particle eta as function of trigger particle pT ...
TH1F * fhPhotonNgammaMoreAverageToNgamma
! number of gammas with ene. more than average ene divided by no. of gammas
Int_t nptbins
TH2F * fhMCJetChNPartVsPt
! generated N parts vs pt charged jet
TH1F * fhPhotonRatioAveEneToMostEne
! ratio average energy to most energetic photon
TH2F * fhJetEtaVsNpartInJet
! Eta vs number of particles in jet for all jets
Double_t fJetAreaFraction
Jet area fraction X in X*pi*R^2, default 0.6.
TH2F * fhPhotonBkgRhoVsNcells
! average energy in one cell vs n cells
Double_t fMCJetCh150ConePt
MC gen charged jet (pt^particles>150MeV/c),R=0.4 pt.
TH2F * fhMCJetCh150ConeEtaPhi
! generated jet eta vs phi charged jet (pt>150 MeV/c) R=0.4
void MakeAnalysisFillAOD()
Particle-Jet Correlation Analysis, fill AODs.
TH2F * fhMCJetChNPart150ConeVsPt
! generated N parts (pt>150 MeV/c) vs pt charged jet R=0.4