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