]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGGA/CaloTrackCorrelations/AliAnaParticleJetFinderCorrelation.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
comment out unused data member, fix some coverity messages
[u/mrichter/AliRoot.git] / PWGGA / CaloTrackCorrelations / AliAnaParticleJetFinderCorrelation.cxx
CommitLineData
1c5acb87 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 **************************************************************************/
1c5acb87 15
16//_________________________________________________________________________
17// Class for the analysis of particle (direct gamma) -jet (jet found with finder) correlations
9717eff6 18//*-- Author: Gustavo Conesa (LNF-INFN)
19//*-- Modified: Adam Matyja (INP PAN Cracow)
1c5acb87 20//////////////////////////////////////////////////////////////////////////////
21
22
23// --- ROOT system ---
24#include "TH2F.h"
477d6cee 25#include "TClonesArray.h"
9415d854 26#include "TClass.h"
1c5acb87 27//#include "Riostream.h"
28
29//---- AliRoot system ----
30#include "AliCaloTrackReader.h"
31#include "AliAODJet.h"
32#include "AliAnaParticleJetFinderCorrelation.h"
1c5acb87 33#include "AliAODPWG4ParticleCorrelation.h"
88f9563f 34#include "AliVTrack.h"
477d6cee 35#include "AliAODCaloCluster.h"
36#include "AliAODEvent.h"
1c5acb87 37
290c1e1b 38//jets
39#include "AliAODJetEventBackground.h"
40#include "TRandom2.h"
41
477d6cee 42ClassImp(AliAnaParticleJetFinderCorrelation)
1c5acb87 43
44
745913ae 45//________________________________________________________________________
290c1e1b 46AliAnaParticleJetFinderCorrelation::AliAnaParticleJetFinderCorrelation() :
47AliAnaCaloTrackCorrBaseClass(),
48 fDeltaPhiMaxCut(0.), fDeltaPhiMinCut(0.), fRatioMaxCut(0.), fRatioMinCut(0.),
49 fConeSize(0), fPtThresholdInCone(0),fUseJetRefTracks(0),
50 fMakeCorrelationInHistoMaker(0), fSelectIsolated(0),
51 fJetConeSize(0.4),fJetMinPt(5),fJetAreaFraction(0.6),
9717eff6 52 //fNonStandardJetFromReader(kTRUE),
53 fJetBranchName("jets"),
290c1e1b 54 fBackgroundJetFromReader(kTRUE),fBkgJetBranchName("jets"),
55 fGammaConeSize(0.3),fUseBackgroundSubtractionGamma(0),fSaveGJTree(0),
56 fMostEnergetic(kFALSE),fMostOpposite(kTRUE), fUseHistogramJetBkg(kTRUE),
57 fUseHistogramTracks(kTRUE),fUseHistogramJetTracks(kTRUE),fGenerator(0),
58 fhDeltaEta(0), /*fhDeltaPhi(0),*/fhDeltaPhiCorrect(0),fhDeltaPhi0PiCorrect(0), fhDeltaPt(0), fhPtRatio(0), fhPt(0),
59 fhFFz(0),fhFFxi(0),fhFFpt(0),fhNTracksInCone(0),
60 fhJetFFz(0),fhJetFFxi(0),fhJetFFpt(0),fhJetFFzCor(0),fhJetFFxiCor(0),
61 fhBkgFFz(),fhBkgFFxi(),fhBkgFFpt(),fhBkgNTracksInCone(),fhBkgSumPtInCone(),fhBkgSumPtnTracksInCone(),//<<---new
62 fhNjetsNgammas(0),fhCuts(0),
63 fhDeltaEtaBefore(0),fhDeltaPhiBefore(0),fhDeltaPtBefore(0),fhPtRatioBefore(0),
64 fhPtBefore(0),fhDeltaPhi0PiCorrectBefore(0),
65 fhJetPtBefore(0),fhJetPt(0),fhJetPtMostEne(0),fhJetPhi(0),fhJetEta(0),fhJetEtaVsPt(0),
66 fhJetPhiVsEta(0),fhJetEtaVsNpartInJet(0),fhJetEtaVsNpartInJetBkg(0),fhJetChBkgEnergyVsPt(0),fhJetChAreaVsPt(0),/*fhJetNjet(0),*/
67 fhTrackPhiVsEta(0),fhTrackAveTrackPt(0),fhJetNjetOverPtCut(),
68/*fhJetChBkgEnergyVsPtEtaGt05(0),fhJetChBkgEnergyVsPtEtaLe05(0),fhJetChAreaVsPtEtaGt05(0),fhJetChAreaVsPtEtaLe05(0),*/
69 fhJetChBkgEnergyVsArea(0),fhJetRhoVsPt(0),fhJetRhoVsCentrality(0),//fhJetBkgRho(0),
70 fhJetNparticlesInJet(0),fhJetDeltaEtaDeltaPhi(0),fhJetDeltaEtaDeltaPhiAllTracks(0),
71 fhJetAveTrackPt(0),fhJetNtracksInJetAboveThr(),fhJetRatioNTrkAboveToNTrk(),fhJetNtrackRatioMostEne(),
72 fhJetNtrackRatioJet5GeV(),fhJetNtrackRatioLead5GeV(),
73 fhBkgJetBackground(),fhBkgJetSigma(),fhBkgJetArea(),fhPhotonPtMostEne(0),
74 fhPhotonAverageEnergy(0),fhPhotonRatioAveEneToMostEne(0),fhPhotonAverageEnergyMinus1(0),fhPhotonRatioAveEneMinus1ToMostEne(0),
75 fhPhotonNgammaMoreAverageToNgamma(0),fhPhotonNgammaMoreAverageMinus1ToNgamma(0),fhPhotonNgammaOverPtCut(),
76 fhPhotonBkgRhoVsNtracks(0),fhPhotonBkgRhoVsNclusters(0),fhPhotonBkgRhoVsCentrality(0),
77 fhPhotonBkgRhoVsNcells(0),fhPhotonPt(0),fhPhotonPtCorrected(0),fhPhotonPtCorrectedZoom(0),fhPhotonPtDiff(0),
78 fhPhotonPtDiffVsCentrality(0),fhPhotonPtDiffVsNcells(0),fhPhotonPtDiffVsNtracks(0),fhPhotonPtDiffVsNclusters(0),
79 fhPhotonSumPtInCone(0),fhPhotonSumPtCorrectInCone(0),fhPhotonSumPtChargedInCone(0),
80 fhSelectedJetPhiVsEta(0),fhSelectedJetChBkgEnergyVsPtJet(0),fhSelectedJetChAreaVsPtJet(0),fhSelectedJetNjet(0),fhSelectedNtracks(0),
81 fhSelectedTrackPhiVsEta(0),fhCuts2(0),
82 fhSelectedPhotonNLMVsPt(0),fhSelectedPhotonLambda0VsPt(0), fhRandomPhiEta(),
83fTreeGJ (0),
84fGamPt (0),
85fGamLambda0 (0),
86fGamNLM (0),
87fGamSumPtCh (0),
88fGamTime (0),
89fGamNcells (0),
90fGamEta (0),
91fGamPhi (0),
92fGamSumPtNeu(0),
93fGamNtracks (0),
94fGamNclusters(0),
95fGamAvEne (0),
96fJetPhi (0),
97fJetEta (0),
98fJetPt (0),
99fJetBkgChEne(0),
100fJetArea (0),
101fJetNtracks (0),
102fJetNtracks1(0),
103fJetNtracks2(0),
104fJetRho(0),
105fEventNumber(0),
106fNtracks (0),
107fZvertex (0),
108fCentrality (0),
109fIso(0),
110fGamRho(0)
111
1c5acb87 112{
477d6cee 113 //Default Ctor
290c1e1b 114 // printf("constructor\n");
477d6cee 115
116 //Initialize parameters
117 InitParameters();
290c1e1b 118 for(Int_t i=0;i<10;i++){
119 fhJetNjetOverPtCut[i] = 0;
120 fhPhotonNgammaOverPtCut[i] = 0;
121 }
122 fGenerator = new TRandom2();
123 fGenerator->SetSeed(0);
1c5acb87 124}
1c5acb87 125
290c1e1b 126//___________________________________________________________________
127AliAnaParticleJetFinderCorrelation::~AliAnaParticleJetFinderCorrelation(){
128 delete fGenerator;
129}
130
131
745913ae 132//___________________________________________________________________
1c5acb87 133TList * AliAnaParticleJetFinderCorrelation::GetCreateOutputObjects()
134{
477d6cee 135 // Create histograms to be saved in output file and
136 // store them in fOutputContainer
290c1e1b 137 // printf("GetCreateOutputObjects\n");
138
477d6cee 139 TList * outputContainer = new TList() ;
140 outputContainer->SetName("ParticleJetFinderHistos") ;
141
745913ae 142 Int_t nptbins = GetHistogramRanges()->GetHistoPtBins();
143 // Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins();
144 // Int_t netabins = GetHistogramRanges()->GetHistoEtaBins();
145 Float_t ptmax = GetHistogramRanges()->GetHistoPtMax();
146 // Float_t phimax = GetHistogramRanges()->GetHistoPhiMax();
147 // Float_t etamax = GetHistogramRanges()->GetHistoEtaMax();
148 Float_t ptmin = GetHistogramRanges()->GetHistoPtMin();
149 // Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
150// Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
477d6cee 151
290c1e1b 152// fhDeltaPhi = new TH2F("DeltaPhi","#phi_{trigger} - #phi_{jet} vs p_{T trigger}",nptbins,ptmin,ptmax,100,-6,4);
153// fhDeltaPhi->SetYTitle("#Delta #phi");
154// fhDeltaPhi->SetXTitle("p_{T trigger} (GeV/c)");
155// outputContainer->Add(fhDeltaPhi);
156
157 fhDeltaPhiCorrect = new TH2F("DeltaPhiCorrect","#phi_{trigger} - #phi_{jet} vs p_{T trigger}",nptbins,ptmin,ptmax,100,0,6.5);
158 fhDeltaPhiCorrect->SetYTitle("#Delta #phi");
159 fhDeltaPhiCorrect->SetXTitle("p_{T trigger} (GeV/c)");
160 outputContainer->Add(fhDeltaPhiCorrect);
161
162 fhDeltaPhi0PiCorrect = new TH2F("DeltaPhi0PiCorrect","#phi_{trigger} - #phi_{jet} (0,#pi) vs p_{T trigger}",nptbins,ptmin,ptmax,100,0,3.5);
163 fhDeltaPhi0PiCorrect->SetYTitle("#Delta #phi");
164 fhDeltaPhi0PiCorrect->SetXTitle("p_{T trigger} (GeV/c)");
165 outputContainer->Add(fhDeltaPhi0PiCorrect);
166
167
168 fhDeltaEta = new TH2F("DeltaEta","#eta_{trigger} - #eta_{jet} vs p_{T trigger}",nptbins,ptmin,ptmax,100,-2,2);
477d6cee 169 fhDeltaEta->SetYTitle("#Delta #eta");
170 fhDeltaEta->SetXTitle("p_{T trigger} (GeV/c)");
171 outputContainer->Add(fhDeltaEta);
172
290c1e1b 173 fhDeltaPt = new TH2F("DeltaPt","p_{T trigger} - p_{T jet} vs p_{T trigger}",nptbins,ptmin,ptmax,150,-50,100);
174 fhDeltaPt->SetYTitle("#Delta p_{T}");
477d6cee 175 fhDeltaPt->SetXTitle("p_{T trigger} (GeV/c)");
176 outputContainer->Add(fhDeltaPt);
177
290c1e1b 178 fhPtRatio = new TH2F("PtRatio","p_{T jet} / p_{T trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,200,0,2.);
477d6cee 179 fhPtRatio->SetYTitle("ratio");
180 fhPtRatio->SetXTitle("p_{T trigger} (GeV/c)");
181 outputContainer->Add(fhPtRatio);
182
9415d854 183 fhPt = new TH2F("Pt","p_{T jet} vs p_{T trigger}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
184 fhPt->SetYTitle("p_{T jet}(GeV/c)");
477d6cee 185 fhPt->SetXTitle("p_{T trigger} (GeV/c)");
186 outputContainer->Add(fhPt);
187
188 fhFFz = new TH2F("FFz","z = p_{T i charged}/p_{T trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,200,0.,2);
189 fhFFz->SetYTitle("z");
190 fhFFz->SetXTitle("p_{T trigger}");
191 outputContainer->Add(fhFFz) ;
1c5acb87 192
477d6cee 193 fhFFxi = new TH2F("FFxi","#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger}", nptbins,ptmin,ptmax,100,0.,10.);
194 fhFFxi->SetYTitle("#xi");
195 fhFFxi->SetXTitle("p_{T trigger}");
196 outputContainer->Add(fhFFxi) ;
197
290c1e1b 198 fhFFpt = new TH2F("FFpt","p_{T i charged} vs p_{T trigger}", nptbins,ptmin,ptmax,100,0.,50.);
477d6cee 199 fhFFpt->SetYTitle("p_{T charged hadron}");
200 fhFFpt->SetXTitle("p_{T trigger}");
201 outputContainer->Add(fhFFpt) ;
202
290c1e1b 203 fhNTracksInCone = new TH2F("NTracksInCone","Number of tracks in cone vs p_{T trigger}", nptbins,ptmin,ptmax,100,0.,150.);
477d6cee 204 fhNTracksInCone->SetYTitle("p_{T charged hadron}");
205 fhNTracksInCone->SetXTitle("p_{T trigger}");
206 outputContainer->Add(fhNTracksInCone) ;
207
290c1e1b 208 //FF according to jet axis
209 fhJetFFz = new TH2F("JetFFz","z = p_{T i charged}/p_{T jet} vs p_{T jet}",nptbins,ptmin,ptmax,200,0.,2);
210 fhJetFFz->SetYTitle("z");
211 fhJetFFz->SetXTitle("p_{T jet}");
212 outputContainer->Add(fhJetFFz) ;
213
214 fhJetFFxi = new TH2F("JetFFxi","#xi = ln(p_{T jet}/p_{T i charged}) vs p_{T jet}", nptbins,ptmin,ptmax,100,0.,10.);
215 fhJetFFxi->SetYTitle("#xi");
216 fhJetFFxi->SetXTitle("p_{T jet}");
217 outputContainer->Add(fhJetFFxi) ;
218
219 fhJetFFpt = new TH2F("JetFFpt","p_{T i charged} vs p_{T jet}", nptbins,ptmin,ptmax,100,0.,50.);
220 fhJetFFpt->SetYTitle("p_{T charged hadron}");
221 fhJetFFpt->SetXTitle("p_{T jet}");
222 outputContainer->Add(fhJetFFpt) ;
223
224 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);
225 fhJetFFzCor->SetYTitle("z");
226 fhJetFFzCor->SetXTitle("p_{T jet}");
227 outputContainer->Add(fhJetFFzCor) ;
228
229 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.);
230 fhJetFFxiCor->SetYTitle("#xi");
231 fhJetFFxiCor->SetXTitle("p_{T jet}");
232 outputContainer->Add(fhJetFFxiCor) ;
233
234
235 //background FF
236 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);
237 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);
238 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);
239 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);
240 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);
241 for(Int_t i=0;i<5;i++){
242 fhBkgFFz[i]->SetYTitle("z");
243 fhBkgFFz[i]->SetXTitle("p_{T trigger}");
244 outputContainer->Add(fhBkgFFz[i]) ;
245 }
246
247 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.);
248 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.);
249 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.);
250 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.);
251 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.);
252 for(Int_t i=0;i<5;i++){
253 fhBkgFFxi[i]->SetYTitle("#xi");
254 fhBkgFFxi[i]->SetXTitle("p_{T trigger}");
255 outputContainer->Add(fhBkgFFxi[i]) ;
256 }
257
258 fhBkgFFpt[0] = new TH2F("BkgFFptRC", "p_{T i charged} vs p_{T trigger} Bkg RC", nptbins,ptmin,ptmax,100,0.,50.);
259 fhBkgFFpt[1] = new TH2F("BkgFFptPCG", "p_{T i charged} vs p_{T trigger} Bkg PCG", nptbins,ptmin,ptmax,100,0.,50.);
260 fhBkgFFpt[2] = new TH2F("BkgFFptPCJ", "p_{T i charged} vs p_{T trigger} Bkg PCJ", nptbins,ptmin,ptmax,100,0.,50.);
261 fhBkgFFpt[3] = new TH2F("BkgFFptMP", "p_{T i charged} vs p_{T trigger} Bkg MP", nptbins,ptmin,ptmax,100,0.,50.);
262 fhBkgFFpt[4] = new TH2F("BkgFFptTest","p_{T i charged} vs p_{T trigger} Bkg Test", nptbins,ptmin,ptmax,100,0.,50.);
263 for(Int_t i=0;i<5;i++){
264 fhBkgFFpt[i]->SetYTitle("p_{T charged hadron}");
265 fhBkgFFpt[i]->SetXTitle("p_{T trigger}");
266 outputContainer->Add(fhBkgFFpt[i]) ;
267 }
268
269 fhBkgNTracksInCone[0] = new TH2F("BkgNTracksInConeRC", "Number of tracks in cone vs p_{T trigger} Bkg RC", nptbins,ptmin,ptmax,100,0.,150.);
270 fhBkgNTracksInCone[1] = new TH2F("BkgNTracksInConePCG", "Number of tracks in cone vs p_{T trigger} Bkg PCG", nptbins,ptmin,ptmax,100,0.,150.);
271 fhBkgNTracksInCone[2] = new TH2F("BkgNTracksInConePCJ", "Number of tracks in cone vs p_{T trigger} Bkg PCJ", nptbins,ptmin,ptmax,100,0.,150.);
272 fhBkgNTracksInCone[3] = new TH2F("BkgNTracksInConeMP", "Number of tracks in cone vs p_{T trigger} Bkg MP", nptbins,ptmin,ptmax,100,0.,150.);
273 fhBkgNTracksInCone[4] = new TH2F("BkgNTracksInConeTest","Number of tracks in cone vs p_{T trigger} Bkg Test", nptbins,ptmin,ptmax,100,0.,150.);
274 for(Int_t i=0;i<5;i++){
275 fhBkgNTracksInCone[i]->SetYTitle("Number of tracks");
276 fhBkgNTracksInCone[i]->SetXTitle("p_{T trigger}");
277 outputContainer->Add(fhBkgNTracksInCone[i]) ;
278 }
279
280 fhBkgSumPtInCone[0] = new TH2F("BkgSumPtInConeRC", "Sum P_{T} in cone vs p_{T trigger} Bkg RC", nptbins,ptmin,ptmax,100,0.,100.);
281 fhBkgSumPtInCone[1] = new TH2F("BkgSumPtInConePCG", "Sum P_{T} in cone vs p_{T trigger} Bkg PCG", nptbins,ptmin,ptmax,100,0.,100.);
282 fhBkgSumPtInCone[2] = new TH2F("BkgSumPtInConePCJ", "Sum P_{T} in cone vs p_{T trigger} Bkg PCJ", nptbins,ptmin,ptmax,100,0.,100.);
283 fhBkgSumPtInCone[3] = new TH2F("BkgSumPtInConeMP", "Sum P_{T} in cone vs p_{T trigger} Bkg MP", nptbins,ptmin,ptmax,100,0.,100.);
284 fhBkgSumPtInCone[4] = new TH2F("BkgSumPtInConeTest","Sum P_{T} in cone vs p_{T trigger} Bkg Test", nptbins,ptmin,ptmax,100,0.,100.);
285 for(Int_t i=0;i<5;i++){
286 fhBkgSumPtInCone[i]->SetYTitle("Sum P_{T}");
287 fhBkgSumPtInCone[i]->SetXTitle("p_{T trigger}");
288 outputContainer->Add(fhBkgSumPtInCone[i]) ;
289 }
290
291 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.);
292 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.);
293 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.);
294 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.);
295 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.);
296 for(Int_t i=0;i<5;i++){
297 fhBkgSumPtnTracksInCone[i]->SetYTitle("Sum p_{T}/Number of tracks");
298 fhBkgSumPtnTracksInCone[i]->SetXTitle("p_{T trigger}");
299 outputContainer->Add(fhBkgSumPtnTracksInCone[i]) ;
300 }
301
302
303 //temporary histograms
304 fhNjetsNgammas = new TH2F("NjetsNgammas"," Number of jets vs number of gammas in event",20,0.,100.,10,0.,80.);
305 fhNjetsNgammas->SetYTitle("Number of gammas");
306 fhNjetsNgammas->SetXTitle("Number of jets");
307 outputContainer->Add(fhNjetsNgammas) ;
308
309 fhCuts = new TH1F("Cuts"," Cuts",10,0.,10.);
310 fhCuts->SetYTitle("Counts");
311 fhCuts->SetXTitle("Cut number");
312 outputContainer->Add(fhCuts) ;
313
314 fhDeltaPhiBefore = new TH2F("DeltaPhiBefore","#phi_{trigger} - #phi_{jet} vs p_{T trigger}",nptbins,ptmin,ptmax,100,0,6.5);
315 fhDeltaPhiBefore->SetYTitle("#Delta #phi");
316 fhDeltaPhiBefore->SetXTitle("p_{T trigger} (GeV/c)");
317 outputContainer->Add(fhDeltaPhiBefore);
318
319 fhDeltaEtaBefore = new TH2F("DeltaEtaBefore","#eta_{trigger} - #eta_{jet} vs p_{T trigger}",nptbins,ptmin,ptmax,100,-2,2);
320 fhDeltaEtaBefore->SetYTitle("#Delta #eta");
321 fhDeltaEtaBefore->SetXTitle("p_{T trigger} (GeV/c)");
322 outputContainer->Add(fhDeltaEtaBefore);
323
324 fhDeltaPtBefore = new TH2F("DeltaPtBefore","p_{T trigger} - p_{T jet} vs p_{T trigger}",nptbins,ptmin,ptmax,100,-50,50);
325 fhDeltaPtBefore->SetYTitle("#Delta p_{T}");
326 fhDeltaPtBefore->SetXTitle("p_{T trigger} (GeV/c)");
327 outputContainer->Add(fhDeltaPtBefore);
328
329 fhPtRatioBefore = new TH2F("PtRatioBefore","p_{T jet} / p_{T trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,200,0,2.);
330 fhPtRatioBefore->SetYTitle("ratio");
331 fhPtRatioBefore->SetXTitle("p_{T trigger} (GeV/c)");
332 outputContainer->Add(fhPtRatioBefore);
333
334 fhPtBefore = new TH2F("PtBefore","p_{T jet} vs p_{T trigger}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
335 fhPtBefore->SetYTitle("p_{T jet}(GeV/c)");
336 fhPtBefore->SetXTitle("p_{T trigger} (GeV/c)");
337 outputContainer->Add(fhPtBefore);
338
339 fhDeltaPhi0PiCorrectBefore = new TH2F("DeltaPhi0PiCorrectBefore","#phi_{trigger} - #phi_{jet} (0,#pi) vs p_{T trigger}",nptbins,ptmin,ptmax,100,0,3.5);
340 fhDeltaPhi0PiCorrectBefore->SetYTitle("#Delta #phi");
341 fhDeltaPhi0PiCorrectBefore->SetXTitle("p_{T trigger} (GeV/c)");
342 outputContainer->Add(fhDeltaPhi0PiCorrectBefore);
343
344 //temporary jet histograms
345 fhJetPtBefore = new TH1F("JetPtBefore","JetPtBefore",150,-50,100);
346 fhJetPtBefore->SetYTitle("Counts");
347 fhJetPtBefore->SetXTitle("p_{T jet}(GeV/c)");
348 outputContainer->Add(fhJetPtBefore) ;
349
350 fhJetPt = new TH1F("JetPt","JetPt",150,-50,100);
351 fhJetPt->SetYTitle("Counts");
352 fhJetPt->SetXTitle("p_{T jet}(GeV/c)");
353 outputContainer->Add(fhJetPt) ;
354
355 fhJetPtMostEne = new TH1F("JetPtMostEne","JetPtMostEne",150,0,150);
356 fhJetPtMostEne->SetYTitle("Counts");
357 fhJetPtMostEne->SetXTitle("p_{T jet}(GeV/c)");
358 outputContainer->Add(fhJetPtMostEne) ;
359
360 fhJetPhi = new TH1F("JetPhi","JetPhi",130,0,6.5);
361 fhJetPhi->SetYTitle("Counts");
362 fhJetPhi->SetXTitle("#phi_{jet}");
363 outputContainer->Add(fhJetPhi) ;
364
365 fhJetEta = new TH1F("JetEta","JetEta",100,-1,1);
366 fhJetEta->SetYTitle("Counts");
367 fhJetEta->SetXTitle("#eta_{jet}");
368 outputContainer->Add(fhJetEta) ;
369
370 fhJetEtaVsPt = new TH2F("JetEtaVsPt","JetEtaVsPt",100,0,100,50,-1,1);
371 fhJetEtaVsPt->SetYTitle("#eta_{jet}");
372 fhJetEtaVsPt->SetXTitle("p_{T,jet}(GeV/c)");
373 outputContainer->Add(fhJetEtaVsPt) ;
374
375 fhJetPhiVsEta = new TH2F("JetPhiVsEta","JetPhiVsEta",65,0,6.5,50,-1,1);
376 fhJetPhiVsEta->SetYTitle("#eta_{jet}");
377 fhJetPhiVsEta->SetXTitle("#phi_{jet}");
378 outputContainer->Add(fhJetPhiVsEta) ;
379
380 fhJetEtaVsNpartInJet= new TH2F("JetEtaVsNpartInJet","JetEtaVsNpartInJet",50,-1,1,100,0.,200.);
381 fhJetEtaVsNpartInJet->SetYTitle("N_{tracks-in-jet}");
382 fhJetEtaVsNpartInJet->SetXTitle("#eta_{jet}");
383 outputContainer->Add(fhJetEtaVsNpartInJet) ;
384
385 fhJetEtaVsNpartInJetBkg= new TH2F("JetEtaVsNpartInJetBkg","JetEtaVsNpartInJetBkg",50,-1,1,100,0.,200.);
386 fhJetEtaVsNpartInJetBkg->SetYTitle("N_{tracks-in-jet}");
387 fhJetEtaVsNpartInJetBkg->SetXTitle("#eta_{jet}");
388 outputContainer->Add(fhJetEtaVsNpartInJetBkg) ;
389
390 fhJetChBkgEnergyVsPt = new TH2F("JetBkgChEnergyVsPt","JetBkgChEnergyVsPt",100,0,100,90,0,90);
391 fhJetChBkgEnergyVsPt->SetYTitle("Jet Bkg Energy (GeV)");
392 fhJetChBkgEnergyVsPt->SetXTitle("p_{T,jet} (GeV/c)");
393 outputContainer->Add(fhJetChBkgEnergyVsPt);
394
395 fhJetChAreaVsPt = new TH2F("JetChAreaVsPt","JetChAreaVsPt",100,0,100,50,0,1);
396 fhJetChAreaVsPt->SetYTitle("Jet Area");
397 fhJetChAreaVsPt->SetXTitle("p_{T,jet} (GeV/c)");
398 outputContainer->Add(fhJetChAreaVsPt);
399
400 if(IsHistogramTracks()){
401 fhTrackPhiVsEta = new TH2F("TrackPhiVsEta","TrackPhiVsEta",65,0,6.5,50,-1,1);
402 fhTrackPhiVsEta->SetYTitle("#eta_{track}");
403 fhTrackPhiVsEta->SetXTitle("#phi_{track}");
404 outputContainer->Add(fhTrackPhiVsEta) ;
405
406 fhTrackAveTrackPt = new TH1F("TrackAveTrackPt","TrackAveTrackPt",45,0,1.5);
407 fhTrackAveTrackPt->SetYTitle("Counts");
408 fhTrackAveTrackPt->SetXTitle("Average p_{T,track} (GeV/c)");
409 outputContainer->Add(fhTrackAveTrackPt);
410
411 }//end of IsHistogramTracks()
412
413 for(Int_t i=0;i<10;i++){
414 fhJetNjetOverPtCut[i] = new TH1F(Form("JetNjetOverPtCut%d", i),Form("JetNjetOverPtCut%d", i),100,0,100);
415 fhJetNjetOverPtCut[i]->SetYTitle("Counts");
416 fhJetNjetOverPtCut[i]->SetXTitle("N_{jets} over threshold");
417 outputContainer->Add(fhJetNjetOverPtCut[i]);
418 }
419
420 fhJetChBkgEnergyVsArea = new TH2F("JetBkgChEnergyVsArea","JetBkgChEnergyVsArea",100,0,100,70,0,0.7);
421 fhJetChBkgEnergyVsArea->SetXTitle("Jet Bkg Energy (GeV)");
422 fhJetChBkgEnergyVsArea->SetYTitle("Area");
423 outputContainer->Add(fhJetChBkgEnergyVsArea);
424
425 fhJetRhoVsPt = new TH2F("JetRhoVsPt","JetRhoVsPt",100,0,100,100,0,150);
426 fhJetRhoVsPt->SetYTitle("Rho");
427 fhJetRhoVsPt->SetXTitle("p_{T,jet} (GeV/c)");
428 outputContainer->Add(fhJetRhoVsPt);
429
430 if(IsHistogramJetBkg()){
431 fhJetRhoVsCentrality = new TH2F("JetRhoVsCentrality","JetRhoVsCentrality",100,0,100,100,0,200);
432 fhJetRhoVsCentrality->SetYTitle("Rho");
433 fhJetRhoVsCentrality->SetXTitle("Centrality");
434 outputContainer->Add(fhJetRhoVsCentrality);
435 }
436
437 fhJetNparticlesInJet = new TH1F("JetNparticlesInJet","JetNparticlesInJet",100,0,200);
438 fhJetNparticlesInJet->SetXTitle("N^{particles}");
439 fhJetNparticlesInJet->SetYTitle("N^{jets}");
440 outputContainer->Add(fhJetNparticlesInJet);
441
442 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);
443 fhJetDeltaEtaDeltaPhi->SetXTitle("#Delta #eta^{jet-track}");
444 fhJetDeltaEtaDeltaPhi->SetYTitle("#Delta #phi^{jet-track}");
445 outputContainer->Add(fhJetDeltaEtaDeltaPhi );
446
447
448 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);
449 fhJetDeltaEtaDeltaPhiAllTracks->SetXTitle("#Delta #eta^{jet-track}");
450 fhJetDeltaEtaDeltaPhiAllTracks->SetYTitle("#Delta #phi^{jet-track}");
451 outputContainer->Add(fhJetDeltaEtaDeltaPhiAllTracks);
452
453
454 if(IsHistogramJetTracks()){
455 fhJetAveTrackPt = new TH1F("JetAveTrackPt","JetAveTrackPt",45,0.,1.5);
456 fhJetAveTrackPt->SetXTitle("Average p_{T,track} (GeV/c)");
457 fhJetAveTrackPt->SetYTitle("Counts");
458 outputContainer->Add(fhJetAveTrackPt);
459
460 for(Int_t i=0;i<6;i++){
461 if(i==0) fhJetNtracksInJetAboveThr[i] = new TH2F(Form("JetNtracksInJetAboveThr%d", i),Form("JetNtracksInJetAboveThr%d", i),100,0,100,100,0,200);
462 else fhJetNtracksInJetAboveThr[i] = new TH2F(Form("JetNtracksInJetAboveThr%d", i),Form("JetNtracksInJetAboveThr%d", i),100,0,100,100,0,100);
463 fhJetNtracksInJetAboveThr[i]->SetXTitle("p_{T,jet} (GeV/c)");
464 fhJetNtracksInJetAboveThr[i]->SetYTitle("N_{tracks} over threshold");
465 outputContainer->Add(fhJetNtracksInJetAboveThr[i]);
466 }
467
468 for(Int_t i=0;i<5;i++){
469 fhJetRatioNTrkAboveToNTrk[i] = new TH2F(Form("JetRatioNTrkAboveToNTrk%d", i),Form("JetRatioNTrkAboveToNTrk%d", i),100,0,100,40,0,1);
470 fhJetRatioNTrkAboveToNTrk[i]->SetXTitle("p_{T,jet} (GeV/c)");
471 fhJetRatioNTrkAboveToNTrk[i]->SetYTitle("Ratio N_{tracks} over threshold to N_{tracks}");
472 outputContainer->Add(fhJetRatioNTrkAboveToNTrk[i]);
473
474 fhJetNtrackRatioMostEne[i] = new TH2F(Form("JetNtrackRatioMostEne%d", i),Form("JetNtrackRatioMostEne%d", i),100,0,100,40,0,1);
475 fhJetNtrackRatioMostEne[i]->SetXTitle("p_{T,jet} (GeV/c)");
476 fhJetNtrackRatioMostEne[i]->SetYTitle("Ratio N_{tracks} over threshold to N_{tracks}");
477 outputContainer->Add(fhJetNtrackRatioMostEne[i]);
478
479 fhJetNtrackRatioJet5GeV[i] = new TH2F(Form("JetNtrackRatioJet5GeV%d", i),Form("JetNtrackRatioJet5GeV%d", i),100,0,100,40,0,1);
480 fhJetNtrackRatioJet5GeV[i]->SetXTitle("p_{T,jet} (GeV/c)");
481 fhJetNtrackRatioJet5GeV[i]->SetYTitle("Ratio N_{tracks} over threshold to N_{tracks}");
482 outputContainer->Add(fhJetNtrackRatioJet5GeV[i]);
483
484 fhJetNtrackRatioLead5GeV[i] = new TH2F(Form("JetNtrackRatioLead5GeV%d", i),Form("JetNtrackRatioLead5GeV%d", i),100,0,100,40,0,1);
485 fhJetNtrackRatioLead5GeV[i]->SetXTitle("p_{T,jet} (GeV/c)");
486 fhJetNtrackRatioLead5GeV[i]->SetYTitle("Ratio N_{tracks} over threshold to N_{tracks}");
487 outputContainer->Add(fhJetNtrackRatioLead5GeV[i]);
488 }
489 }//end of if IsHistogramJetTracks
490
491 //temporary background jets histograms
492 if(IsHistogramJetBkg()){
493 for(Int_t i=0;i<4;i++){
494 fhBkgJetBackground[i] = new TH1F(Form("BkgJetBackground%d", i),Form("BkgJetBackground%d", i),100,0,200);
495 fhBkgJetBackground[i]->SetXTitle("<#rho> (GeV/c)");
496 fhBkgJetBackground[i]->SetYTitle("Counts");
497 outputContainer->Add(fhBkgJetBackground[i]);
498
499 fhBkgJetSigma[i] = new TH1F(Form("BkgJetSigma%d", i),Form("BkgJetSigma%d", i),100,0,50);
500 fhBkgJetSigma[i]->SetXTitle("#sigma (GeV/c)");
501 fhBkgJetSigma[i]->SetYTitle("Counts");
502 outputContainer->Add(fhBkgJetSigma[i]);
503
504 fhBkgJetArea[i] = new TH1F(Form("BkgJetArea%d", i),Form("BkgJetArea%d", i),100,0,1);
505 fhBkgJetArea[i]->SetXTitle("<A>");
506 fhBkgJetArea[i]->SetYTitle("Counts");
507 outputContainer->Add(fhBkgJetArea[i]);
508 }
509 }
510
511 //temporary photon histograms
512 fhPhotonPtMostEne = new TH1F("PhotonPtMostEne","PhotonPtMostEne",100,0,100);
513 fhPhotonPtMostEne->SetYTitle("Counts");
514 fhPhotonPtMostEne->SetXTitle("p_{T,#gamma} (GeV/c)");
515 outputContainer->Add(fhPhotonPtMostEne);
516
517// fhPhotonIndexMostEne = new TH1F("PhotonIndexMostEne","PhotonIndexMostEne",100,0,100);
518// fhPhotonIndexMostEne->SetYTitle("Counts");
519// fhPhotonIndexMostEne->SetXTitle("Index");
520// outputContainer->Add(fhPhotonIndexMostEne);
521
522 fhPhotonAverageEnergy = new TH1F("PhotonAverageEnergy","PhotonAverageEnergy",100,0,10);
523 fhPhotonAverageEnergy->SetYTitle("Counts");
524 fhPhotonAverageEnergy->SetXTitle("p_{T,#gamma} (GeV/c)");
525 outputContainer->Add(fhPhotonAverageEnergy);
526
527 fhPhotonRatioAveEneToMostEne = new TH1F("PhotonRatioAveEneToMostEne","PhotonRatioAveEneToMostEne",100,0,1);
528 fhPhotonRatioAveEneToMostEne->SetYTitle("Counts");
529 fhPhotonRatioAveEneToMostEne->SetXTitle("Ratio");
530 outputContainer->Add(fhPhotonRatioAveEneToMostEne);
531
532 fhPhotonAverageEnergyMinus1 = new TH1F("PhotonAverageEnergyMinus1","PhotonAverageEnergyMinus1",100,0,10);
533 fhPhotonAverageEnergyMinus1->SetYTitle("Counts");
534 fhPhotonAverageEnergyMinus1->SetXTitle("p_{T,#gamma} (GeV/c)");
535 outputContainer->Add(fhPhotonAverageEnergyMinus1);
536
537 fhPhotonRatioAveEneMinus1ToMostEne = new TH1F("PhotonRatioAveEneMinus1ToMostEne","PhotonRatioAveEneMinus1ToMostEne",100,0,1);
538 fhPhotonRatioAveEneMinus1ToMostEne->SetYTitle("Counts");
539 fhPhotonRatioAveEneMinus1ToMostEne->SetXTitle("Ratio");
540 outputContainer->Add(fhPhotonRatioAveEneMinus1ToMostEne);
541
542 fhPhotonNgammaMoreAverageToNgamma = new TH1F("PhotonNgammaMoreAverageToNgamma","PhotonNgammaMoreAverageToNgamma",100,0,1);
543 fhPhotonNgammaMoreAverageToNgamma->SetYTitle("Counts");
544 fhPhotonNgammaMoreAverageToNgamma->SetXTitle("Ratio");
545 outputContainer->Add(fhPhotonNgammaMoreAverageToNgamma);
546
547 fhPhotonNgammaMoreAverageMinus1ToNgamma = new TH1F("PhotonNgammaMoreAverageMinus1ToNgamma","PhotonNgammaMoreAverageMinus1ToNgamma",100,0,1);
548 fhPhotonNgammaMoreAverageMinus1ToNgamma->SetYTitle("Counts");
549 fhPhotonNgammaMoreAverageMinus1ToNgamma->SetXTitle("Ratio");
550 outputContainer->Add(fhPhotonNgammaMoreAverageMinus1ToNgamma);
551
552 for(Int_t i=0;i<10;i++){
553 fhPhotonNgammaOverPtCut[i] = new TH1F(Form("PhotonNgammaOverPtCut%d",i),Form("PhotonNgammaOverPtCut%d",i),100,0,100);
554 fhPhotonNgammaOverPtCut[i]->SetYTitle("Counts");
555 fhPhotonNgammaOverPtCut[i]->SetXTitle("N_{#gamma} over threshold");
556 outputContainer->Add(fhPhotonNgammaOverPtCut[i]);
557 }
558
559 fhPhotonBkgRhoVsNtracks = new TH2F("PhotonBkgRhoVsNtracks","PhotonBkgRhoVsNtracks",200,0,2500,75,0,1.5);
560 //fhPhotonBkgRhoVsNtracks->SetXTitle("Counts");
561 fhPhotonBkgRhoVsNtracks->SetXTitle("Ntracks");
562 fhPhotonBkgRhoVsNtracks->SetYTitle("Rho");
563 outputContainer->Add(fhPhotonBkgRhoVsNtracks);
564
565 fhPhotonBkgRhoVsNclusters = new TH2F("PhotonBkgRhoVsNclusters","PhotonBkgRhoVsNclusters",50,0,100,75,0,1.5);
566 fhPhotonBkgRhoVsNclusters->SetXTitle("Nclusters");
567 fhPhotonBkgRhoVsNclusters->SetYTitle("Rho");
568 outputContainer->Add(fhPhotonBkgRhoVsNclusters);
569
570 fhPhotonBkgRhoVsCentrality = new TH2F("PhotonBkgRhoVsCentrality","PhotonBkgRhoVsCentrality",100,0,100,75,0,1.5);
571 fhPhotonBkgRhoVsCentrality->SetXTitle("Centrality");
572 fhPhotonBkgRhoVsCentrality->SetYTitle("Rho");
573 outputContainer->Add(fhPhotonBkgRhoVsCentrality);
574
575 fhPhotonBkgRhoVsNcells = new TH2F("PhotonBkgRhoVsNcells","PhotonBkgRhoVsNcells",100,0,200,75,0,1.5);
576 fhPhotonBkgRhoVsNcells->SetXTitle("N_{cells}");
577 fhPhotonBkgRhoVsNcells->SetYTitle("Rho");
578 outputContainer->Add(fhPhotonBkgRhoVsNcells);
579
580 fhPhotonPt = new TH1F("PhotonPt","PhotonPt",220,-10,100);
581 fhPhotonPt->SetXTitle("p_{T,#gamma} (GeV/c)");
582 fhPhotonPt->SetYTitle("Counts");
583 outputContainer->Add(fhPhotonPt);
584
585 fhPhotonPtCorrected = new TH1F("PhotonPtCorrected","PhotonPtCorrected",220,-10,100);
586 fhPhotonPtCorrected->SetXTitle("p_{T,#gamma} (GeV/c)");
587 fhPhotonPtCorrected->SetYTitle("Counts");
588 outputContainer->Add(fhPhotonPtCorrected);
589
590 fhPhotonPtDiff = new TH1F("PhotonPtDiff","PhotonPtDiff",50,0,10);
591 fhPhotonPtDiff->SetXTitle("p_{T,#gamma} (GeV/c)");
592 fhPhotonPtDiff->SetYTitle("Counts");
593 outputContainer->Add(fhPhotonPtDiff);
594
595 fhPhotonPtDiffVsNtracks = new TH2F("PhotonPtDiffVsNtracks","PhotonPtDiffVsNtracks",200,0,2500,50,0,5);
596 fhPhotonPtDiffVsNtracks->SetXTitle("N_{tracks}");
597 fhPhotonPtDiffVsNtracks->SetYTitle("<#rho^{#gamma}>*N_{cells}");
598 outputContainer->Add(fhPhotonPtDiffVsNtracks);
599
600 fhPhotonPtDiffVsNclusters = new TH2F("PhotonPtDiffVsNclusters","PhotonPtDiffVsNclusters",50,0,100,50,0,5);
601 fhPhotonPtDiffVsNclusters->SetXTitle("N_{clusters}");
602 fhPhotonPtDiffVsNclusters->SetYTitle("<#rho^{#gamma}>*N_{cells}");
603 outputContainer->Add(fhPhotonPtDiffVsNclusters);
604
605 fhPhotonPtDiffVsCentrality = new TH2F("PhotonPtDiffVsCentrality","PhotonPtDiffVsCentrality",100,0,100,50,0,5);
606 fhPhotonPtDiffVsCentrality->SetXTitle("Centrality");
607 fhPhotonPtDiffVsCentrality->SetYTitle("<#rho^{#gamma}>*N_{cells}");
608 outputContainer->Add(fhPhotonPtDiffVsCentrality);
609
610 fhPhotonPtDiffVsNcells = new TH2F("PhotonPtDiffVsNcells","PhotonPtDiffVsNcells",100,0,200,50,0,5);
611 fhPhotonPtDiffVsNcells->SetXTitle("N_{cells}");
612 fhPhotonPtDiffVsNcells->SetYTitle("<#rho^{#gamma}>*N_{cells}");
613 outputContainer->Add(fhPhotonPtDiffVsNcells);
614
615
616 fhPhotonPtCorrectedZoom = new TH1F("PhotonPtCorrectedZoom","PhotonPtCorrectedZoom",200,-5,5);
617 fhPhotonPtCorrectedZoom->SetXTitle("p_{T,#gamma} (GeV/c)");
618 fhPhotonPtCorrectedZoom->SetYTitle("Counts");
619 outputContainer->Add(fhPhotonPtCorrectedZoom);
620
621 fhPhotonSumPtInCone = new TH1F("PhotonSumPtInCone","PhotonSumPtInCone",100,0,100);
622 fhPhotonSumPtInCone->SetXTitle("#Sigma p_{T,#gamma} (GeV/c)");
623 fhPhotonSumPtInCone->SetYTitle("Counts");
624 outputContainer->Add(fhPhotonSumPtInCone);
625
626 fhPhotonSumPtCorrectInCone = new TH1F("PhotonSumPtCorrectInCone","PhotonSumPtCorrectInCone",100,-20,80);
627 fhPhotonSumPtCorrectInCone->SetXTitle("#Sigma p_{T,#gamma} (GeV/c)");
628 fhPhotonSumPtCorrectInCone->SetYTitle("Counts");
629 outputContainer->Add(fhPhotonSumPtCorrectInCone);
630
631 fhPhotonSumPtChargedInCone = new TH1F("PhotonSumPtChargedInCone","PhotonSumPtChargedInCone",100,0,100);
632 fhPhotonSumPtChargedInCone->SetXTitle("#Sigma p_{T,#gamma}^{ch} (GeV/c)");
633 fhPhotonSumPtChargedInCone->SetYTitle("Counts");
634 outputContainer->Add(fhPhotonSumPtChargedInCone);
635
636
637 //temporary jet histograms after selection
638 fhSelectedJetPhiVsEta = new TH2F("SelectedJetSelectedPhiVsEta","SelectedJetPhiVsEta",65,0,6.5,50,-1,1);
639 fhSelectedJetPhiVsEta->SetYTitle("#eta_{jet}");
640 fhSelectedJetPhiVsEta->SetXTitle("#phi_{jet}");
641 outputContainer->Add(fhSelectedJetPhiVsEta) ;
642
643 fhSelectedJetChBkgEnergyVsPtJet = new TH2F("SelectedJetBkgChEnergyVsPtJet","SelectedJetBkgChEnergyVsPtJet",100,0,100,90,0,90);
644 fhSelectedJetChBkgEnergyVsPtJet->SetYTitle("Jet Bkg Energy (GeV)");
645 fhSelectedJetChBkgEnergyVsPtJet->SetXTitle("p_{T,jet} (GeV/c)");
646 outputContainer->Add(fhSelectedJetChBkgEnergyVsPtJet);
647
648 fhSelectedJetChAreaVsPtJet = new TH2F("SelectedJetChAreaVsPtJet","SelectedJetChAreaVsPtJet",100,0,100,50,0,1);
649 fhSelectedJetChAreaVsPtJet->SetYTitle("Jet Area");
650 fhSelectedJetChAreaVsPtJet->SetXTitle("p_{T,jet} (GeV/c)");
651 outputContainer->Add(fhSelectedJetChAreaVsPtJet);
652
653 fhSelectedJetNjet = new TH1F("SelectedJetNjet","SelectedJetNjet",100,0,100);
654 fhSelectedJetNjet->SetYTitle("Counts");
655 fhSelectedJetNjet->SetXTitle("N_{jets} per event");
656 outputContainer->Add(fhSelectedJetNjet);
657
658 fhSelectedNtracks = new TH1F("SelectedNtracks","SelectedNtracks",100,0,2000);
659 fhSelectedNtracks->SetYTitle("Counts");
660 fhSelectedNtracks->SetXTitle("N_{tracks} per event");
661 outputContainer->Add(fhSelectedNtracks);
662
663 fhSelectedTrackPhiVsEta = new TH2F("SelectedTrackPhiVsEta","SelectedTrackPhiVsEta",65,0,6.5,50,-1,1);
664 fhSelectedTrackPhiVsEta->SetYTitle("#eta_{track}");
665 fhSelectedTrackPhiVsEta->SetXTitle("#phi_{track}");
666 outputContainer->Add(fhSelectedTrackPhiVsEta) ;
667
668 fhCuts2 = new TH1F("Cuts2","Cuts2",10,0.,10.);
669 fhCuts2->SetYTitle("Counts");
670 fhCuts2->SetXTitle("Cut number");
671 outputContainer->Add(fhCuts2);
672
673 fhSelectedPhotonNLMVsPt = new TH2F("SelectedPhotonNLMVsPt","SelectedPhotonNLMVsPt",100,0,100,10,0,10);
674 fhSelectedPhotonNLMVsPt->SetYTitle("NLM");
675 fhSelectedPhotonNLMVsPt->SetXTitle("p_{T,#gamma} (GeV/c)");
676 outputContainer->Add(fhSelectedPhotonNLMVsPt);
677
678 fhSelectedPhotonLambda0VsPt = new TH2F("SelectedPhotonLambda0VsPt","SelectedPhotonLambda0VsPt",100,0,100,50,0,5);
679 fhSelectedPhotonLambda0VsPt->SetYTitle("#lambda_{0}");
680 fhSelectedPhotonLambda0VsPt->SetXTitle("p_{T,#gamma} (GeV/c)");
681 outputContainer->Add(fhSelectedPhotonLambda0VsPt);
682
683 //random
684 fhRandomPhiEta[0] = new TH2F("RandomPhiEtaRC", "RandomPhiEtaRC", 100,0,6.5,100,-1.,1.);
685 fhRandomPhiEta[1] = new TH2F("RandomPhiEtaPCG", "RandomPhiEtaPerpConePhoton",100,0,6.5,100,-1.,1.);
686 fhRandomPhiEta[2] = new TH2F("RandomPhiEtaPCJ", "RandomPhiEtaPerpConeJet", 100,0,6.5,100,-1.,1.);
687 fhRandomPhiEta[3] = new TH2F("RandomPhiEtaMP", "RandomPhiEtaMidPoint", 100,0,6.5,100,-1.,1.);
688 fhRandomPhiEta[4] = new TH2F("RandomPhiEtaTest","RandomPhiEtaTest", 100,0,6.5,100,-1.,1.);
689 for(Int_t i=0;i<5;i++){
690 fhRandomPhiEta[i]->SetXTitle("#phi");
691 fhRandomPhiEta[i]->SetYTitle("#eta");
692 outputContainer->Add(fhRandomPhiEta[i]);
693 }
694
695
696 //tree with data
697 if(fSaveGJTree){
698 fTreeGJ=new TTree("fTreeGJ","fTreeGJ");
699 fTreeGJ->Branch("fGamPt" ,&fGamPt ,"fGamPt/D");//! fGamPt
700 fTreeGJ->Branch("fGamLambda0" ,&fGamLambda0 ,"fGamLambda0/D");
701 fTreeGJ->Branch("fGamNLM" ,&fGamNLM ,"fGamNLM/I");
702 fTreeGJ->Branch("fGamSumPtCh" ,&fGamSumPtCh ,"fGamSumPtCh/D");
703 fTreeGJ->Branch("fGamNtracks" ,&fGamNtracks ,"fGamNtracks/I");
704 fTreeGJ->Branch("fGamTime" ,&fGamTime ,"fGamTime/D");
705 fTreeGJ->Branch("fGamNcells" ,&fGamNcells ,"fGamNcells/I");
706 fTreeGJ->Branch("fGamEta" ,&fGamEta ,"fGamEta/D");
707 fTreeGJ->Branch("fGamPhi" ,&fGamPhi ,"fGamPhi/D");
708 fTreeGJ->Branch("fGamSumPtNeu",&fGamSumPtNeu ,"fGamSumPtNeu/D");
709 fTreeGJ->Branch("fGamNclusters" ,&fGamNclusters ,"fGamNclusters/I");
710 fTreeGJ->Branch("fGamAvEne" ,&fGamAvEne ,"fGamAvEne/D");
711 fTreeGJ->Branch("fJetPhi" ,&fJetPhi ,"fJetPhi/D");
712 fTreeGJ->Branch("fJetEta" ,&fJetEta ,"fJetEta/D");
713 fTreeGJ->Branch("fJetPt" ,&fJetPt ,"fJetPt/D");
714 fTreeGJ->Branch("fJetBkgChEne",&fJetBkgChEne ,"fJetBkgChEne/D");
715 fTreeGJ->Branch("fJetArea" ,&fJetArea ,"fJetArea/D");
716 fTreeGJ->Branch("fJetNtracks" ,&fJetNtracks ,"fJetNtracks/I");
717 fTreeGJ->Branch("fJetNtracks1" ,&fJetNtracks1 ,"fJetNtracks1/I");
718 fTreeGJ->Branch("fJetNtracks2" ,&fJetNtracks2 ,"fJetNtracks2/I");
719 fTreeGJ->Branch("fJetRho" ,&fJetRho ,"fJetRho/D");
720 fTreeGJ->Branch("fEventNumber",&fEventNumber ,"fEventNumber/I");
721 fTreeGJ->Branch("fNtracks" ,&fNtracks ,"fNtracks/I");
722 fTreeGJ->Branch("fZvertex" ,&fZvertex ,"fZvertex/D");
723 fTreeGJ->Branch("fCentrality" ,&fCentrality ,"fCentrality/D");
724 fTreeGJ->Branch("fIso" ,&fIso ,"fIso/O");
725 fTreeGJ->Branch("fGamRho" ,&fGamRho ,"fGamRho/D");
726
727 outputContainer->Add(fTreeGJ);
728 }
729
477d6cee 730 return outputContainer;
a3aebfff 731
1c5acb87 732}
733
c5693f62 734//_______________________________________________________
1c5acb87 735void AliAnaParticleJetFinderCorrelation::InitParameters()
736{
290c1e1b 737 // printf("InitParameters\n");
1c5acb87 738 //Initialize the parameters of the analysis.
a3aebfff 739 SetInputAODName("PWG4Particle");
740 AddToHistogramsName("AnaJetFinderCorr_");
741
1c5acb87 742 fDeltaPhiMinCut = 1.5 ;
743 fDeltaPhiMaxCut = 4.5 ;
744 fRatioMaxCut = 1.2 ;
745 fRatioMinCut = 0.3 ;
746 fConeSize = 0.3 ;
747 fPtThresholdInCone = 0.5 ;
748 fUseJetRefTracks = kFALSE ;
749 fMakeCorrelationInHistoMaker = kFALSE ;
477d6cee 750 fSelectIsolated = kFALSE;
290c1e1b 751 fJetConeSize = 0.4 ;
752 fJetMinPt = 5 ; //GeV/c
753 fJetAreaFraction = 0.6 ;
754 fJetBranchName = "jets";
755 fBkgJetBranchName = "jets";
756 fGammaConeSize = 0.3;
757 fUseBackgroundSubtractionGamma = kFALSE;
758 fSaveGJTree = kTRUE;
759 fMostEnergetic = kFALSE;
760 fMostOpposite = kTRUE;
761 fUseHistogramJetBkg = kTRUE;
762 fUseHistogramTracks = kTRUE;
763 fUseHistogramJetTracks = kTRUE;
764
c5693f62 765}
1c5acb87 766
290c1e1b 767//__________________________________________________________________
768Int_t AliAnaParticleJetFinderCorrelation::SelectJet(AliAODPWG4Particle * particle, TClonesArray *aodRecJets)
1c5acb87 769{
290c1e1b 770 //Input for jets is TClonesArray *aodRecJets
477d6cee 771 //Returns the index of the jet that is opposite to the particle
290c1e1b 772 // printf(" SelectJet ");
477d6cee 773
290c1e1b 774 Double_t particlePt=particle->Pt();
775 if(fUseBackgroundSubtractionGamma) {
776 Int_t clusterIDtmp = particle->GetCaloLabel(0) ;
777 Int_t nCells=0;
778 AliVCluster *cluster=0;
779 if(!(clusterIDtmp<0) ){
780 Int_t iclustmp = -1;
781 TObjArray* clusters = GetEMCALClusters();
782 cluster = FindCluster(clusters,clusterIDtmp,iclustmp);
783 nCells = cluster->GetNCells();
784 }
785 particlePt-=(fGamRho*nCells);
786 }
787 if(particlePt<=0) {
788 //printf("Particle with negative or 0 pt\n");
789 return -2;
790 }
477d6cee 791
290c1e1b 792 Int_t njets = aodRecJets->GetEntriesFast();
76c58de9 793 AliAODJet * jet = 0 ;
477d6cee 794 Int_t index = -1;
290c1e1b 795 Double_t dphiprev= 10000;
796 Double_t particlePhi=particle->Phi();
797 Double_t deltaPhi=-10000.;// in the range (0; 2*pi)
798 Double_t jetPt=0.;
799
477d6cee 800 for(Int_t ijet = 0; ijet < njets ; ijet++){
290c1e1b 801 jet = dynamic_cast<AliAODJet*>(aodRecJets->At(ijet));
3098a292 802 if(!jet)
803 {
804 AliInfo("Jet not in container");
805 continue;
806 }
290c1e1b 807 fhCuts2->Fill(2.,1.);
808 jetPt=jet->Pt();
809 if(fBackgroundJetFromReader ){
810 jetPt-= (fJetRho * jet->EffectiveAreaCharged() );
811 }
812 if(jetPt<0.) continue;
813 //put jet eta requirement here |eta_jet|<0.9-jet_cone_size
814 fhCuts2->Fill(3.,1.);
815 if(TMath::Abs(jet->Eta()) > (0.9 - fJetConeSize) ) continue;
816 fhCuts2->Fill(4.,1.);
817 // if(jet->Pt()<5) continue;
818 if(jetPt<fJetMinPt) continue;
819 fhCuts2->Fill(5.,1.);
820 if(jet->EffectiveAreaCharged()<fJetAreaFraction*TMath::Pi()*fJetConeSize*fJetConeSize) continue;
821 fhCuts2->Fill(6.,1.);
822 //printf("jet found\n");
823 Double_t deltaPhi0pi = TMath::Abs(particle->Phi()-jet->Phi());
824 Double_t ratio = jetPt/particlePt;
825
826 deltaPhi = particlePhi - jet->Phi() ;
827 if ( deltaPhi0pi > TMath::Pi() ) deltaPhi0pi = 2. * TMath::Pi() - deltaPhi0pi ;
828 if(deltaPhi<0) deltaPhi +=(TMath::Pi()*2.);
829
830 fhDeltaPtBefore ->Fill(particlePt, particlePt - jetPt);
831 fhDeltaPhiBefore->Fill(particlePt, deltaPhi);
832 fhDeltaEtaBefore->Fill(particlePt, particle->Eta() - jet->Eta());
833 fhPtRatioBefore ->Fill(particlePt, jetPt/particlePt);
834 fhPtBefore ->Fill(particlePt, jetPt);
835
836 fhDeltaPhi0PiCorrectBefore->Fill(particlePt, deltaPhi0pi);//good
837
838 if(GetDebug() > 5)
839 printf("AliAnaParticleJetFinderCorrelation::SelectJet() - Jet %d, Ratio pT %2.3f, Delta phi %2.3f\n",ijet,ratio,deltaPhi);
840
841 if((deltaPhi > fDeltaPhiMinCut) && (deltaPhi < fDeltaPhiMaxCut) &&
477d6cee 842 (ratio > fRatioMinCut) && (ratio < fRatioMaxCut) &&
290c1e1b 843 (TMath::Abs(deltaPhi-TMath::Pi()) < TMath::Abs(dphiprev-TMath::Pi()) )){//In case there is more than one jet select the most opposite.
844 dphiprev = deltaPhi;
845 index = ijet ;
846 }//Selection cuts
847 }//AOD Jet loop
848
849 return index ;
850
1c5acb87 851}
852
853//__________________________________________________________________
290c1e1b 854void AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD()
855{
477d6cee 856 //Particle-Jet Correlation Analysis, fill AODs
290c1e1b 857 // printf("I use MakeAnalysisFillAOD\n");
f37fa8d2 858 //Get the event, check if there are AODs available, if not, skip this analysis
859 AliAODEvent * event = NULL;
290c1e1b 860
861 // cout<<"GetReader()->GetOutputEvent() "<<GetReader()->GetOutputEvent()<<endl;
862 // cout<<"GetReader()->GetDataType() "<<GetReader()->GetDataType() <<endl;
863 // cout<<"AliCaloTrackReader::kAOD "<<AliCaloTrackReader::kAOD<<endl;
864 // cout<<"GetReader()->GetInputEvent() "<<GetReader()->GetInputEvent()<<endl;
865
866 if(GetReader()->GetOutputEvent())
f37fa8d2 867 {
290c1e1b 868 event = dynamic_cast<AliAODEvent*>(GetReader()->GetOutputEvent());
f37fa8d2 869 }
290c1e1b 870 else if(GetReader()->GetDataType() == AliCaloTrackReader::kAOD)
871 {
872 event = dynamic_cast<AliAODEvent*>(GetReader()->GetInputEvent());
f37fa8d2 873 }
290c1e1b 874 else
f37fa8d2 875 {
7672b8a0 876 if(GetDebug() > 3) AliInfo("There are no jets available for this analysis");
f37fa8d2 877 return;
878 }
290c1e1b 879
7672b8a0 880 if(!GetInputAODBranch() || !event)
881 {
882 AliFatal(Form("No input particles in AOD with name branch < %s > \n",
883 GetInputAODName().Data()));
884 return; // Trick coverity
477d6cee 885 }
a3aebfff 886
7672b8a0 887 if(strcmp(GetInputAODBranch()->GetClass()->GetName(), "AliAODPWG4ParticleCorrelation"))
888 {
889 AliFatal(Form("Wrong type of AOD object, change AOD class name in input AOD: It should be <AliAODPWG4ParticleCorrelation> and not <%s>",
890 GetInputAODBranch()->GetClass()->GetName()));
891 return; // Trick coverity
9415d854 892 }
290c1e1b 893
894 //
895 // non-standard jets
896 //
897 Int_t nJets=-1;
898 TClonesArray *aodRecJets = 0;
9717eff6 899 //if(IsNonStandardJetFromReader()){//jet branch from reader
900 if(GetDebug() > 3) AliInfo(Form("GetNonStandardJets function (from reader) is called"));
901 aodRecJets = GetNonStandardJets();
902 if(GetDebug() > 3) AliInfo(Form("aodRecJets %p",aodRecJets));
903 if(aodRecJets==0x0)
7672b8a0 904 {
290c1e1b 905 if(GetDebug() > 3) event->Print();
7672b8a0 906 AliFatal("List of jets is null");
907 return;
290c1e1b 908 }
9717eff6 909 nJets=aodRecJets->GetEntries();
910 if(GetDebug() > 3) printf("nJets %d\n",nJets);
911 //}
290c1e1b 912 //end of new part
913
914 if(nJets==0) {
915 //printf("Why number of jets = 0? Check what type of collision it is. If PbPb -problem.\n");
916 return;
917 }
918
919 //
920 //Background jets
921 //
922 AliAODJetEventBackground* aodBkgJets = 0;
923 if(IsBackgroundJetFromReader()){//jet branch from reader
924 if(GetDebug() > 3) printf("GetBackgroundJets function is called\n");
925 aodBkgJets = GetBackgroundJets();
926 if(GetDebug() > 3) printf("aodBkgJets %p\n",aodBkgJets);
927 if(aodBkgJets==0x0){
928 if(GetDebug() > 3) event->Print();
929 abort();
930 }
931 if(GetDebug() > 3) aodBkgJets->Print("c");
932 }
933
934 Double_t rhoEvent=0.;
935 if(aodBkgJets && IsBackgroundJetFromReader() ) {
936 rhoEvent = aodBkgJets->GetBackground(2);//hardcoded
937 }
938 fJetRho = rhoEvent;
939
940
941 Int_t ntrig = GetInputAODBranch()->GetEntriesFast() ;
477d6cee 942 if(GetDebug() > 3){
a3aebfff 943 printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - Begin jet finder correlation analysis, fill AODs \n");
944 printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - In particle branch aod entries %d\n", ntrig);
290c1e1b 945 printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - In standard jet branch aod entries %d\n", event->GetNJets());
946 printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - In non standard jet branch aod entries %d\n", nJets);
477d6cee 947 }
f37fa8d2 948
290c1e1b 949 //if(nJets==0) return;//to speed up
950 // cout<<"ntrig po return "<<ntrig<<endl;
951
952 //
953 //calculate average cell energy without most energetic photon
954 //
955 Double_t medianPhotonRho=0.;
956 TObjArray* clusters = GetEMCALClusters();
957 Int_t clusterIDtmp;
958 Int_t iclustmp = -1;
959 AliVCluster *cluster=0;
960
961 if(IsBackgroundSubtractionGamma()){
962 //
963 // Find most energetic photon without background subtraction
964 //
965 Double_t maxPt=0.;
966 Int_t maxIndex=-1;
967 AliAODPWG4ParticleCorrelation* particlecorr =0;
968 for(Int_t iaod = 0; iaod < ntrig ; iaod++){
969 particlecorr = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
970 if(particlecorr->Pt() > maxPt) {
971 maxPt = particlecorr->Pt();
972 maxIndex = iaod;
973 }
477d6cee 974 }
290c1e1b 975
976 //
977 // calculate background energy per cell
978 //
979 Int_t numberOfcells=0;
980 if(ntrig>1){
981 Double_t *photonRhoArr=new Double_t[ntrig-1];
982 Int_t photonRhoArrayIndex=0;
983 for(Int_t iaod = 0; iaod < ntrig ; iaod++){
984 particlecorr = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
985 if(iaod==maxIndex) continue;
986 clusterIDtmp = particlecorr->GetCaloLabel(0) ;
987 if(clusterIDtmp < 0) continue;
988 cluster = FindCluster(clusters,clusterIDtmp,iclustmp);
989 photonRhoArr[photonRhoArrayIndex]=particlecorr->Pt()/ cluster->GetNCells();
990 numberOfcells+=cluster->GetNCells();
991 photonRhoArrayIndex++;
992 }
993 if(photonRhoArrayIndex>0) medianPhotonRho=TMath::Median(photonRhoArrayIndex,photonRhoArr);
e410f14e 994 delete [] photonRhoArr;
290c1e1b 995 }
996 }//end of if background calculation for gamma
997 fGamRho = medianPhotonRho;
998
999
1000 //
1001 //take most energetic photon and most energetic jet and combine
1002 //
1003 if(fMostEnergetic){
1004 //
1005 // most energetic photon with background subtraction
1006 //
1007 Double_t mostEnePhotonPt=0.;
1008 Int_t indexMostEnePhoton=-1;
1009 AliAODPWG4ParticleCorrelation* particle =0;
1010 Double_t ptCorrect=0.;
1011 Int_t nCells=0;
1012 for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1013 particle = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1014 clusterIDtmp = particle->GetCaloLabel(0) ;
1015 if(!(clusterIDtmp<0)){
1016 cluster = FindCluster(clusters,clusterIDtmp,iclustmp);
1017 nCells = cluster->GetNCells();
1018 }
1019 ptCorrect = particle->Pt() - medianPhotonRho * nCells;
1020 if( ptCorrect > mostEnePhotonPt ){
1021 mostEnePhotonPt = ptCorrect;
1022 indexMostEnePhoton = iaod ;
1023 }
1024 }
1025 // printf ("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - Photon with index %d selected \n",indexMostEnePhoton);
1026 //
1027 // most energetic jet with background subtraction
1028 //
1029 Double_t mostEneJetPt=0.;
1030 Int_t indexMostEneJet=-1;
1031 AliAODJet * jet = 0 ;
1032 //Double_t ptCorrect=0.;
1033 for(Int_t ijet = 0; ijet < nJets ; ijet++){
1034 jet = dynamic_cast<AliAODJet*>(aodRecJets->At(ijet));
1cc70ab5 1035 if(!jet) continue;
290c1e1b 1036 if(TMath::Abs(jet->Eta()) > (0.9 - fJetConeSize) ) continue;
1037 if(jet->EffectiveAreaCharged()<fJetAreaFraction*TMath::Pi()*fJetConeSize*fJetConeSize) continue;
1038 if(jet->Pt()<fJetMinPt) continue;
1039 ptCorrect = jet->Pt() - rhoEvent * jet->EffectiveAreaCharged();
1040 if(ptCorrect > mostEneJetPt){
1041 mostEneJetPt = ptCorrect;
1042 indexMostEneJet = ijet;
1043 }
1044 }
1045 // printf ("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - Jet with index %d selected \n",indexMostEneJet);
1046
1047 //
1048 // assign jet to photon
1049 //
1050 if(indexMostEneJet>=0 && indexMostEnePhoton>=0){
1051 particle = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(indexMostEnePhoton));
1052 jet = dynamic_cast<AliAODJet*>(aodRecJets-> At(indexMostEneJet));
8c8f7d7a 1053 if(jet)particle->SetRefJet(jet);
290c1e1b 1054 }
1055 }//end of take most energetic photon and most ene. jet after bkg subtraction
1056
1057 if(fMostOpposite){
1058 //Bool_t isJetFound=kFALSE;//new
1059 //Loop on stored AOD particles, trigger
1060 for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1061 AliAODPWG4ParticleCorrelation* particle = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1062
1063 //Correlate with jets
1064 Int_t ijet = SelectJet(particle,aodRecJets);//input for jets is TClonesArray
1065 if(ijet > -1){
1066 //isJetFound=kTRUE;
1067 if(GetDebug() > 2) printf ("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - Jet with index %d selected \n",ijet);
1068 AliAODJet *jet = dynamic_cast<AliAODJet*>(aodRecJets-> At(ijet));
8c8f7d7a 1069 if(jet)particle->SetRefJet(jet);
290c1e1b 1070 //printf("Most opposite found\n");
1071 }
1072 } // input aod loop
1073 // if(GetReader()->WriteDeltaAODToFile() && isJetFound) WriteJetsToOutputBranch(aodRecJets);
1074 }//end of take most opposite photon and jet after bkg subtraction
1075
477d6cee 1076
a3aebfff 1077 if(GetDebug() > 1 ) printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - End fill AODs \n");
1c5acb87 1078}
1079
1080//__________________________________________________________________
290c1e1b 1081void AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms()
1c5acb87 1082{
477d6cee 1083 //Particle-Jet Correlation Analysis, fill histograms
290c1e1b 1084 if(GetDebug() > 3 ) {
1085 printf("I use MakeAnalysisFillHistograms\n");
1086 printf("ntracks before iso %d\n",GetCTSTracks()->GetEntriesFast() );
1087 }
477d6cee 1088
f37fa8d2 1089 //Get the event, check if there are AODs available, if not, skip this analysis
1090 AliAODEvent * event = NULL;
290c1e1b 1091
1092 //printf("GetReader()->GetOutputEvent() %d\n",GetReader()->GetOutputEvent() );
1093 //printf("GetReader()->GetDataType() %d\n",GetReader()->GetDataType() );
1094 //printf("AliCaloTrackReader::kAOD %d\n",AliCaloTrackReader::kAOD );
1095 //printf("GetReader()->GetInputEvent() %d\n",GetReader()->GetInputEvent() );
1096
1097 if(GetReader()->GetOutputEvent())
f37fa8d2 1098 {
290c1e1b 1099 event = dynamic_cast<AliAODEvent*>(GetReader()->GetOutputEvent());
f37fa8d2 1100 }
290c1e1b 1101 else if(GetReader()->GetDataType() == AliCaloTrackReader::kAOD)
1102 {
1103 event = dynamic_cast<AliAODEvent*>(GetReader()->GetInputEvent());
f37fa8d2 1104 }
1105 else {
1106 if(GetDebug() > 3) printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - There are no jets available for this analysis\n");
1107 return;
1108 }
1109
f8006433 1110 if(!GetInputAODBranch() || !event){
a3aebfff 1111 printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - No input particles in AOD with name branch < %s > \n",GetInputAODName().Data());
477d6cee 1112 abort();
1113 }
f37fa8d2 1114
290c1e1b 1115 Int_t nJets=-1;
1116 TClonesArray *aodRecJets = 0;
9717eff6 1117 //if(IsNonStandardJetFromReader()){//branch read in reader from reader
1118 if (GetDebug() > 3) printf("GetNonStandardJets function (from reader) is called\n");
1119 aodRecJets = GetNonStandardJets();
1120 if(aodRecJets==0x0){
1121 if(GetDebug() > 3) event->Print();
1122 AliFatal("Jets container not found");
1123 return; // trick coverity
290c1e1b 1124 }
9717eff6 1125 nJets=aodRecJets->GetEntries();
1126 //}
290c1e1b 1127 if(nJets==0) {
1128 // printf("Why number of jets = 0? Check what type of collision it is. If PbPb -problem.\n");
1129 GetReader()->FillInputNonStandardJets();
1130 aodRecJets = GetNonStandardJets();
3098a292 1131 if(aodRecJets) nJets=aodRecJets->GetEntries();
290c1e1b 1132 // printf("nJets = %d\n",nJets);
1133 return;
1134 }
1135
1136 //
1137 //Background jets
1138 //
1139 AliAODJetEventBackground* aodBkgJets = 0;
1140 if(IsBackgroundJetFromReader()){//jet branch from reader
1141 if(GetDebug() > 3) printf("GetBackgroundJets function is called\n");
1142 aodBkgJets = GetBackgroundJets();
1143 if(GetDebug() > 3) printf("aodBkgJets %p\n",aodBkgJets);
1144 if(aodBkgJets==0x0){
1145 if(GetDebug() > 3) event->Print();
1146 abort();
1147 }
1148 if(GetDebug() > 3) aodBkgJets->Print("c");
1149 }
1150
1151
1152 //
1153 // only background jets informations
1154 //
1155 // Float_t pTback = 0;
1156 Double_t rhoEvent=0.;
1157 if(aodBkgJets) {
1158 if(IsBackgroundJetFromReader() ) rhoEvent = aodBkgJets->GetBackground(2);
1159 if(IsHistogramJetBkg()) {
1160 fhJetRhoVsCentrality->Fill(GetEventCentrality(),rhoEvent);
1161 for(Int_t i=0;i<4;i++){
1162 fhBkgJetBackground[i]->Fill(aodBkgJets->GetBackground(i));
1163 fhBkgJetSigma[i]->Fill(aodBkgJets->GetSigma(i));
1164 fhBkgJetArea[i]->Fill(aodBkgJets->GetMeanarea(i));
1165 }
1166 }//end of if fill HistogramJetBkg
1167 }//end if aodBkgJets exists
1168
1169 //
1170 //only track information
1171 //
1172 Int_t nCTSTracks = GetCTSTracks()->GetEntriesFast();
1173 AliAODTrack *aodtrack;
1174 Int_t itrack = 0;
1175 if(IsHistogramTracks()) {
1176 Double_t sumTrackPt=0;
1177 for(itrack = 0; itrack < nCTSTracks ; itrack++){
1178 aodtrack = dynamic_cast <AliAODTrack*>(GetCTSTracks()->At(itrack));
1cc70ab5 1179 if(!aodtrack) continue;
290c1e1b 1180 fhTrackPhiVsEta->Fill(aodtrack->Phi(),aodtrack->Eta());
1181 sumTrackPt+=aodtrack->Pt();
1182 }
1183 if(nCTSTracks)
1184 fhTrackAveTrackPt->Fill(sumTrackPt/nCTSTracks);
1185 }//end if IsHistogramTracks
1186
1187 //
1188 //only jet informations
1189 //
1190 AliAODJet * jettmp = 0 ;
1191 Double_t jetPttmp = 0.;
1192 Double_t var1tmp = 0.;
1193 Double_t var2tmp = 0.;
1194 // fhJetNjet->Fill(nJets);
1195 Double_t ptMostEne=0;
1196 // Int_t indexMostEne=-1;
1197 Int_t nJetsOverThreshold[10]={nJets,0,0,0,0,0,0,0,0,0};
1198 Int_t iCounter=0;
1199 Double_t sumJetTrackPt=0.;
1200 Int_t sumNJetTrack=0;
1201 Int_t nTracksInJet=0;
1202 Int_t itrk=0;
1203 for(Int_t ijet = 0; ijet < nJets ; ijet++){
1204 jettmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ijet));
1cc70ab5 1205 if(!jettmp) continue;
290c1e1b 1206 fhJetPtBefore->Fill(jettmp->Pt());
1207 jetPttmp = jettmp->Pt() - rhoEvent * jettmp->EffectiveAreaCharged();//<<---changed here
1208
1209 //calculate number of tracks above 1,2,3,4,5 GeV in jet
1210 AliVTrack* jettrack = 0x0 ;
1211
1212 Int_t nTrackThrGeV[5]={0,0,0,0,0};
1213 nTracksInJet=(jettmp->GetRefTracks())->GetEntriesFast();
1214 fhJetNparticlesInJet->Fill(nTracksInJet);
1215 if(nTracksInJet==0) continue;
1216 sumNJetTrack+=nTracksInJet;
1217 for(itrack=0;itrack<nTracksInJet;itrack++){
1218 jettrack=(AliVTrack *) ((jettmp->GetRefTracks())->At(itrack));
1219 if(!jettrack) continue;
1220
1221 fhJetDeltaEtaDeltaPhi->Fill(jettmp->Eta()-jettrack->Eta(),jettmp->Phi()-jettrack->Phi());
1222 sumJetTrackPt+=jettrack->Pt();
1223 if(IsHistogramJetTracks()){
1224 if(jettrack->Pt()>1.) nTrackThrGeV[0]++;
1225 if(jettrack->Pt()>2.) nTrackThrGeV[1]++;
1226 if(jettrack->Pt()>3.) nTrackThrGeV[2]++;
1227 if(jettrack->Pt()>4.) nTrackThrGeV[3]++;
1228 if(jettrack->Pt()>5.) nTrackThrGeV[4]++;
1229 }//end of if IsHistogramJetTracks
1230 }//end of jet track loop
1231 //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]);
1232
1233 for(itrack = 0; itrack < nCTSTracks ; itrack++){
1234 aodtrack = dynamic_cast <AliAODTrack*>(GetCTSTracks()->At(itrack));
3098a292 1235 if(aodtrack) fhJetDeltaEtaDeltaPhiAllTracks->Fill(jettmp->Eta()-aodtrack->Eta(),jettmp->Phi()-aodtrack->Phi());
290c1e1b 1236 }
1237
1238
1239 if(IsHistogramJetTracks()){
1240 fhJetNtracksInJetAboveThr[0]->Fill(jetPttmp,nTracksInJet);//all jets
1241
1242 for(itrk=0;itrk<5;itrk++) {
1243 fhJetNtracksInJetAboveThr[itrk+1]->Fill(jetPttmp,nTrackThrGeV[itrk]);//all jets
1244 fhJetRatioNTrkAboveToNTrk[itrk]->Fill(jetPttmp,(Double_t)nTrackThrGeV[itrk]/(Double_t)nTracksInJet);//all jets
1245 }
1246 if(ijet==0){//most ene jet
1247 for(itrk=0;itrk<5;itrk++)
1248 fhJetNtrackRatioMostEne[itrk]->Fill(jetPttmp,(Double_t)nTrackThrGeV[itrk]/(Double_t)nTracksInJet);
1249 }
1250 if(jetPttmp>5){//jet with pt>5GeV/c
1251 for(itrk=0;itrk<5;itrk++)
1252 fhJetNtrackRatioJet5GeV[itrk]->Fill(jetPttmp,(Double_t)nTrackThrGeV[itrk]/(Double_t)nTracksInJet);
1253 }
1254 if(nTrackThrGeV[4]>0){//jet with leading particle pt>5GeV/c
1255 for(itrk=0;itrk<5;itrk++)
1256 fhJetNtrackRatioLead5GeV[itrk]->Fill(jetPttmp,(Double_t)nTrackThrGeV[itrk]/(Double_t)nTracksInJet);
1257 }
1258 }//end of if IsHistogramJetTracks
1259
1260
1261 Double_t effectiveChargedBgEnergy=(IsBackgroundJetFromReader()?rhoEvent * jettmp->EffectiveAreaCharged():jettmp->ChargedBgEnergy());
1262
1263
1264 fhJetChBkgEnergyVsArea->Fill(effectiveChargedBgEnergy,jettmp->EffectiveAreaCharged());
1265 //if(jettmp->EffectiveAreaCharged()>0)
1266 fhJetRhoVsPt->Fill(jetPttmp,jettmp->ChargedBgEnergy()*jettmp->EffectiveAreaCharged());
1267
1268 if(jetPttmp>ptMostEne) {
1269 ptMostEne = jetPttmp;
1270 //indexMostEne=ijet;
1271 }
1272 fhJetPt->Fill(jetPttmp);
1273 fhJetChBkgEnergyVsPt->Fill(jetPttmp,effectiveChargedBgEnergy);
1274 fhJetChAreaVsPt->Fill(jetPttmp,jettmp->EffectiveAreaCharged());
1275 if(GetDebug()>5) printf("ChargedBgEnergy %f EffectiveAreaCharged %f\n", jettmp->ChargedBgEnergy(),jettmp->EffectiveAreaCharged());
1276 for(iCounter=1;iCounter<10;iCounter++){
1277 if(jetPttmp>iCounter) nJetsOverThreshold[iCounter]++;
1278 }
1279
1280 var1tmp = jettmp->Phi();
1281 var2tmp = jettmp->Eta();
1282 fhJetPhi->Fill(var1tmp);
1283 fhJetEta->Fill(var2tmp);
1284 fhJetPhiVsEta->Fill(var1tmp,var2tmp);
1285 fhJetEtaVsPt->Fill(jetPttmp,var2tmp);
1286 fhJetEtaVsNpartInJet->Fill(var2tmp,nTracksInJet);
1287 if(jetPttmp>0)
1288 fhJetEtaVsNpartInJetBkg->Fill(var2tmp,nTracksInJet);
1289
1290 }//end of jet loop
1291 if(IsHistogramJetTracks()){
1292 if(sumNJetTrack>0){
1293 //printf("average track pt %f\n",sumJetTrackPt/sumNJetTrack);
1294 fhJetAveTrackPt->Fill(sumJetTrackPt/sumNJetTrack);
1295 }
1296 }//end of if IsHistogramJetTracks
1297
1298
1299 fhJetPtMostEne->Fill(ptMostEne);
1300 for(iCounter=0;iCounter<10;iCounter++){
1301 fhJetNjetOverPtCut[iCounter]->Fill(nJetsOverThreshold[iCounter]);
1302 nJetsOverThreshold[iCounter]=0;
1303 }
1304
1305 //end of only jet part
1306
1307 //
1308 // Photons
1309 //
1310 Int_t ntrig = GetInputAODBranch()->GetEntriesFast() ;
477d6cee 1311 if(GetDebug() > 1){
a3aebfff 1312 printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - Begin jet finder correlation analysis, fill histograms \n");
1313 printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - In particle branch aod entries %d\n", ntrig);
f37fa8d2 1314 printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - In jet output branch aod entries %d\n", event->GetNJets());
477d6cee 1315 }
290c1e1b 1316 fhNjetsNgammas->Fill(nJets,ntrig);
1317
1318 //if(nJets==0) return;//to speed up
1319 //printf("ntrig %d, njets %d\n",ntrig,nJets);
1320
1321 //
1322 //Fill only temporary photon histograms
1323 //
1324 Double_t maxPt=0.;
1325 Int_t maxIndex=-1;
1326 Double_t tmpPt=0.;
1327 Double_t sumPt=0.;
1328 nJetsOverThreshold[0]=ntrig;
1329 for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1330 AliAODPWG4ParticleCorrelation* particlecorr = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1331 tmpPt = particlecorr->Pt();
1332 if(tmpPt>maxPt) {
1333 maxPt = tmpPt;
1334 maxIndex = iaod;
1335 }
1336 for(iCounter=1;iCounter<10;iCounter++){
1337 if(tmpPt>iCounter) nJetsOverThreshold[iCounter]++;
1338 }
1339 sumPt+=tmpPt;
1340 }
1341 for(iCounter=0;iCounter<10;iCounter++){
1342 fhPhotonNgammaOverPtCut[iCounter]->Fill(nJetsOverThreshold[iCounter]);
1343 // nJetsOverThreshold[iCounter]=0;
1344 }
1345
1346 fGamAvEne=0;
1347 TObjArray* clusters = GetEMCALClusters();
1348 //printf("calculate median bkg energy for photons ");
1349 Double_t medianPhotonRho=0.;
1350 Int_t clusterID;
1351 Int_t iclustmp = -1;
1352 Int_t numberOfcells=0;
1353 AliVCluster *cluster = 0;
1354 if(ntrig>1){
1355 Double_t *photonRhoArr=new Double_t[ntrig-1];
1356 fhPhotonPtMostEne->Fill(maxPt);
1357 // fhPhotonIndexMostEne->Fill(indexMaxPt);
1358 fhPhotonAverageEnergy->Fill(sumPt/ntrig);
1359 fhPhotonRatioAveEneToMostEne->Fill(sumPt/(ntrig*maxPt));
1360 fhPhotonAverageEnergyMinus1->Fill((sumPt-maxPt)/(ntrig-1));
1361 fGamAvEne=(sumPt-maxPt)/(ntrig-1);
1362 fhPhotonRatioAveEneMinus1ToMostEne->Fill((sumPt-maxPt)/((ntrig-1)*maxPt));
1363
1364 Int_t counterGamma=0;
1365 Int_t counterGammaMinus1=0;
1366
1367 Int_t photonRhoArrayIndex=0;
1368 //TObjArray* clusterstmp = GetEMCALClusters();
1369 for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1370 AliAODPWG4ParticleCorrelation* particlecorr = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1371 if( particlecorr->Pt() > sumPt/ntrig ) counterGamma++;
1372 if( particlecorr->Pt() > (sumPt-maxPt)/(ntrig-1) ) counterGammaMinus1++;
1373
1374 if(iaod==maxIndex) continue;
1375 clusterID = particlecorr->GetCaloLabel(0) ;
1376 if(clusterID < 0) continue;
1377 cluster = FindCluster(clusters,clusterID,iclustmp);
1378 photonRhoArr[photonRhoArrayIndex]=particlecorr->Pt()/ cluster->GetNCells();
1379 numberOfcells+=cluster->GetNCells();
1380 photonRhoArrayIndex++;
1381 }
1382 if(photonRhoArrayIndex>0) medianPhotonRho=TMath::Median(photonRhoArrayIndex,photonRhoArr);
e410f14e 1383 delete [] photonRhoArr;
290c1e1b 1384 fhPhotonNgammaMoreAverageToNgamma->Fill((Double_t)counterGamma / (Double_t)ntrig);
1385 fhPhotonNgammaMoreAverageMinus1ToNgamma->Fill((Double_t)counterGammaMinus1 / (Double_t)ntrig);
1386 }
1387 //printf("median = %f\n",medianPhotonRho);
1388 fhPhotonBkgRhoVsNtracks->Fill(GetCTSTracks()->GetEntriesFast(),medianPhotonRho);
1389 fhPhotonBkgRhoVsNclusters->Fill(ntrig,medianPhotonRho);
1390 fhPhotonBkgRhoVsCentrality->Fill(GetEventCentrality(),medianPhotonRho);
1391 fhPhotonBkgRhoVsNcells->Fill(numberOfcells,medianPhotonRho);
1392
1393
1394 AliVCluster *cluster2 = 0;
1395 Double_t photon2Corrected=0;
1396 Double_t sumPtTmp=0.;
1397 Double_t sumPtCorrectTmp=0.;
1398 AliVTrack* trackTmp = 0x0 ;
1399 TVector3 p3Tmp;
1400
1401 for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1402 AliAODPWG4ParticleCorrelation* particlecorr = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1403 clusterID = particlecorr->GetCaloLabel(0) ;
1404 if(clusterID < 0) continue;
1405 cluster = FindCluster(clusters,clusterID,iclustmp);
1406 fhPhotonPt->Fill(particlecorr->Pt());
1407 fhPhotonPtCorrected->Fill(particlecorr->Pt() - cluster->GetNCells() * medianPhotonRho);
1408 fhPhotonPtDiff->Fill(cluster->GetNCells() * medianPhotonRho);
1409 fhPhotonPtDiffVsCentrality->Fill(GetEventCentrality(),cluster->GetNCells() * medianPhotonRho);
1410 fhPhotonPtDiffVsNcells->Fill(numberOfcells,cluster->GetNCells() * medianPhotonRho);
1411 fhPhotonPtDiffVsNtracks->Fill(GetCTSTracks()->GetEntriesFast(),cluster->GetNCells() * medianPhotonRho);
1412 fhPhotonPtDiffVsNclusters->Fill(ntrig,cluster->GetNCells() * medianPhotonRho);
1413
1414 fhPhotonPtCorrectedZoom->Fill(particlecorr->Pt() - cluster->GetNCells() * medianPhotonRho);
1415
1416 //test: sum_pt in the cone 0.3 for each photon
1417 //should be: random fake gamma from MB
1418 //is: each gamma for EMCEGA
1419 sumPtTmp=0.;
1420 sumPtCorrectTmp=0.;
1421
1422 for(Int_t iaod2 = 0; iaod2 < ntrig ; iaod2++){
1423 if(iaod==iaod2) continue;
1424 AliAODPWG4ParticleCorrelation* particlecorr2 = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod2));
1425 clusterID = particlecorr2->GetCaloLabel(0) ;
1426 if(clusterID < 0) continue;
1427 cluster2 = FindCluster(clusters,clusterID,iclustmp);
1428 photon2Corrected = particlecorr2->Pt() - cluster2->GetNCells() * medianPhotonRho;
1429
1430 //if(Pt()<0.5) continue; //<<hardcoded here //FIXME
1431 if( TMath::Sqrt((particlecorr->Eta()-particlecorr2->Eta())*(particlecorr->Eta()-particlecorr2->Eta()) +
1432 (particlecorr->Phi()-particlecorr2->Phi())*(particlecorr->Phi()-particlecorr2->Phi()) )<fGammaConeSize ){//if(/*cone is correct*/){
1433 sumPtTmp+= particlecorr2->Pt();
1434 sumPtCorrectTmp+=photon2Corrected;
1435 }
1436 }
1437 fhPhotonSumPtInCone->Fill(sumPtTmp);
1438 fhPhotonSumPtCorrectInCone->Fill(sumPtCorrectTmp);
1439
1440 //test: sum_pt in the cone 0.3 for each track
1441 //should be: random fake gamma from MB
1442 //is: each gamma for EMCEGA
1443 sumPtTmp=0.;
1444 for(Int_t ipr = 0;ipr < GetCTSTracks()->GetEntriesFast() ; ipr ++){
1445 trackTmp = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
1446 p3Tmp.SetXYZ(trackTmp->Px(),trackTmp->Py(),trackTmp->Pz());
1447 if( TMath::Sqrt((particlecorr->Eta()-p3Tmp.Eta())*(particlecorr->Eta()-p3Tmp.Eta()) +
1448 (particlecorr->Phi()-p3Tmp.Phi())*(particlecorr->Phi()-p3Tmp.Phi()) )<fGammaConeSize ){
1449 sumPtTmp+=p3Tmp.Pt();
1450 }
1451 }//end of loop over tracks
1452 fhPhotonSumPtChargedInCone->Fill(sumPtTmp);
1453 }
1454
1455 //End of Fill temporary photon histograms
1456
1457 //
1458 // Apply background subtraction for photons
1459 //
1460 fGamRho = medianPhotonRho;
1461 if(!IsBackgroundSubtractionGamma()) medianPhotonRho=0;
477d6cee 1462
290c1e1b 1463
1464 //
1465 //Get vertex for cluster momentum calculation <<----new here
1466 //
1467 Double_t vertex[] = {0,0,0} ; //vertex ;
1468 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
1469 GetReader()->GetVertex(vertex);
1470 fZvertex = vertex[2];
1471
1472 //
477d6cee 1473 //Loop on stored AOD particles, trigger
290c1e1b 1474 //
477d6cee 1475 for(Int_t iaod = 0; iaod < ntrig ; iaod++){
290c1e1b 1476 AliAODPWG4ParticleCorrelation* particlecorr = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1477 fhCuts->Fill(0);
1478 fhCuts2->Fill(0.,(Double_t)nJets);
1479 if(GetDebug() > 5) printf("OnlyIsolated %d !particlecorr->IsIsolated() %d \n",OnlyIsolated(), !particlecorr->IsIsolated());
477d6cee 1480
a3aebfff 1481 if(OnlyIsolated() && !particlecorr->IsIsolated()) continue;
290c1e1b 1482 fhCuts->Fill(1);
1483 fhCuts2->Fill(1.,nJets);
1484
1485 if(nJets>0) {
1486 fhCuts->Fill(2);
1487 }
477d6cee 1488
1489 //Recover the jet correlated, found previously.
a3aebfff 1490 AliAODJet * jet = particlecorr->GetJet();
477d6cee 1491 //If correlation not made before, do it now.
1492 if(fMakeCorrelationInHistoMaker){
1493 //Correlate with jets
290c1e1b 1494 Int_t ijet = SelectJet(particlecorr,aodRecJets);//input for jets is TClonesArray
477d6cee 1495 if(ijet > -1){
f37fa8d2 1496 if(GetDebug() > 2) printf ("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - Jet with index %d selected \n",ijet);
290c1e1b 1497 //jet = event->GetJet(ijet);
1498 jet = dynamic_cast<AliAODJet*>(aodRecJets-> At(ijet));
1499
8c8f7d7a 1500 if(jet) particlecorr->SetRefJet(jet);
290c1e1b 1501
477d6cee 1502 }
1503 }
1504
1505 if (!jet) continue ;
290c1e1b 1506 fhCuts->Fill(3);
1507 fhCuts2->Fill(7.,1.);
477d6cee 1508
290c1e1b 1509 //
1510 //Fill Correlation Histograms
1511 //
1512 clusterID = particlecorr->GetCaloLabel(0) ;
1513 if(!(clusterID<0)){
1514 cluster = FindCluster(clusters,clusterID,iclustmp);
1515 //fill tree variables
1516 fGamNcells = cluster->GetNCells();
1517 }
1518 Double_t ptTrig = particlecorr->Pt() - medianPhotonRho * fGamNcells;//<<---changed here
1519 Double_t ptJet = jet->Pt() - rhoEvent * jet->EffectiveAreaCharged();//<<---changed here
a3aebfff 1520 Double_t phiTrig = particlecorr->Phi();
477d6cee 1521 Double_t phiJet = jet->Phi();
a3aebfff 1522 Double_t etaTrig = particlecorr->Eta();
477d6cee 1523 Double_t etaJet = jet->Eta();
290c1e1b 1524 Double_t deltaPhi=phiTrig-phiJet;
1525 if(deltaPhi<0)deltaPhi+=(TMath::Pi()*2.);
477d6cee 1526 //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",
1527 // ptTrig,ptJet, phiJet-phiTrig, etaJet-etaTrig, ptTrig-ptJet, ptJet/ptTrig);
1528 fhDeltaPt ->Fill(ptTrig, ptTrig-ptJet);
290c1e1b 1529 // fhDeltaPhi->Fill(ptTrig, phiTrig-phiJet);//need to be shifted by 2pi
1530
1531 fhDeltaPhiCorrect->Fill(ptTrig, deltaPhi);//correct
1532
1533 Double_t deltaPhiCorrect = TMath::Abs( particlecorr->Phi() - jet->Phi() );
1534 if ( deltaPhiCorrect > TMath::Pi() ) deltaPhiCorrect = 2. * TMath::Pi() - deltaPhiCorrect ;
1535 fhDeltaPhi0PiCorrect->Fill(ptTrig, deltaPhiCorrect);
1536
1537 fhDeltaEta->Fill(ptTrig, etaTrig-etaJet);
477d6cee 1538 fhPtRatio ->Fill(ptTrig, ptJet/ptTrig);
1539 fhPt ->Fill(ptTrig, ptJet);
1540
290c1e1b 1541 fhSelectedJetPhiVsEta->Fill(phiJet,etaJet);
1542 fhSelectedJetChBkgEnergyVsPtJet->Fill(ptJet,(IsBackgroundJetFromReader()?rhoEvent * jet->EffectiveAreaCharged():jet->ChargedBgEnergy()) );
1543 fhSelectedJetChAreaVsPtJet->Fill(ptJet,jet->EffectiveAreaCharged());
1544 fhSelectedJetNjet->Fill(nJets);
1545 fhSelectedNtracks->Fill(GetCTSTracks()->GetEntriesFast());//to be checked
1546 fhSelectedPhotonNLMVsPt->Fill(ptTrig,particlecorr->GetFiducialArea());
1547
1548
1549 if(clusterID < 0 ){
1550 fhSelectedPhotonLambda0VsPt->Fill(ptTrig,-1);
1551 //fill tree variables
1552 fGamLambda0 = -1;
1553 fGamTime = -1;
1554 fGamNcells = 0;
1555 fGamSumPtNeu=0;
1556 }
1557 else{
1558 //Int_t iclus = -1;
1559 // TObjArray* clusters = GetEMCALClusters();
1560 //cluster = FindCluster(clusters,clusterID,iclustmp);
1561 Double_t lambda0=cluster->GetM02();
1562 fhSelectedPhotonLambda0VsPt->Fill(ptTrig,lambda0);
1563 //fill tree variables
1564 fGamLambda0 = lambda0;
1565 fGamTime = cluster->GetTOF();
1566 //fGamNcells = cluster->GetNCells();
1567
1568 fGamSumPtNeu=0;
1569 fGamNclusters=0;
1570 TLorentzVector mom ;
1571 //TVector3 p3Tmp;
1572 //Double_t scalarProduct=0;
1573 //Double_t vectorLength=particlecorr->P();
1574 for(Int_t icalo=0; icalo <clusters->GetEntriesFast(); icalo++){
1575 AliVCluster* calo = (AliVCluster *) clusters->At(icalo);
1576 if(clusterID==calo->GetID()) continue;//the same cluster as trigger
1577 calo->GetMomentum(mom,vertex) ;//Assume that come from vertex in straight line
1578 //printf("min pt %f\n",GetMinPt());
1579 if(mom.Pt()<GetMinPt()) continue; //<<hardcoded here //FIXME 0.5 check if correct
1580 p3Tmp.SetXYZ(mom.Px(),mom.Py(),mom.Pz());
1581 //calculate sum pt in the cone
1582 if( TMath::Sqrt((particlecorr->Eta()-p3Tmp.Eta())*(particlecorr->Eta()-p3Tmp.Eta()) +
1583 (particlecorr->Phi()-p3Tmp.Phi())*(particlecorr->Phi()-p3Tmp.Phi()) )<fGammaConeSize ){
1584 //scalarProduct = particlecorr->Px()*mom.Px() + particlecorr->Py()*mom.Py() + particlecorr->Pz()*mom.Pz();
1585 //scalarProduct/=mom.P();
1586 //scalarProduct/=vectorLength;
1587 //if(scalarProduct>TMath::Cos(0.3)) {//FIXME photon radius
1588 fGamSumPtNeu+=mom.Pt();
1589 fGamNclusters++;
1590 }
1591 }
1592 }
1593
1594 //sum pt of charged tracks in the gamma isolation cone
1595 //starts here
1596 fGamSumPtCh=0;
1597 fGamNtracks=0;
1598 for(itrack = 0; itrack < nCTSTracks ; itrack++){
1599 aodtrack = dynamic_cast <AliAODTrack*>(GetCTSTracks()->At(itrack));
8c8f7d7a 1600 if(!aodtrack) continue;
290c1e1b 1601 fhSelectedTrackPhiVsEta->Fill(aodtrack->Phi(),aodtrack->Eta());//fill histogram here
1602 // if(aodtrack->Pt()<0.15) continue;//hardcoded
1603 if(aodtrack->Pt()<fPtThresholdInCone) continue;
1604 if(!aodtrack->IsHybridGlobalConstrainedGlobal()) continue;
1605 if(TMath::Sqrt((particlecorr->Phi() - aodtrack->Phi())*(particlecorr->Phi() - aodtrack->Phi()) +
1606 (particlecorr->Eta() - aodtrack->Eta())*(particlecorr->Eta() - aodtrack->Eta()) ) <fGammaConeSize ) {
1607 fGamSumPtCh+=aodtrack->Pt();
1608 fGamNtracks++;
1609 }
1610 }
1611 //ends here
1612
1613 // for(Int_t itrack = 0; itrack < nCTSTracks ; itrack++){
1614 // aodtrack = dynamic_cast <AliAODTrack*>(GetCTSTracks()->At(itrack));
1615 // fhSelectedTrackPhiVsEta->Fill(aodtrack->Phi(),aodtrack->Eta());
1616 // }
1617
1618 //
1619 // Background Fragmentation function
1620 //
1621 TVector3 gammaVector,jetVector;
1622 gammaVector.SetXYZ(particlecorr->Px(),particlecorr->Py(),particlecorr->Pz());
1623 jetVector.SetXYZ(jet->Px(),jet->Py(),jet->Pz());
1624 CalculateBkg(gammaVector,jetVector,vertex,1);//jet perp
1625 CalculateBkg(gammaVector,jetVector,vertex,2);//RC
1626 CalculateBkg(gammaVector,jetVector,vertex,3);//mid point
1627 CalculateBkg(gammaVector,jetVector,vertex,4);//gamma perp
1628 //CalculateBkg(gammaVector,jetVector,vertex,5);/test
1629 Double_t angleJetGam = gammaVector.Angle(jetVector);
1630 //printf("angleJetGam %f\n",angleJetGam*180/TMath::Pi());
1631
1632 //
1633 // Fragmentation function
1634 //
477d6cee 1635 Float_t rad = 0, pt = 0, eta = 0, phi = 0;
1636 Int_t npartcone = 0;
1637 TVector3 p3;
477d6cee 1638
1639 Int_t ntracks = 0;
290c1e1b 1640 if(GetDebug()>3){
1641 printf ("fUseJetRefTracks %d\n",fUseJetRefTracks );
1642 printf ("jet->GetRefTracks() %p\n",jet->GetRefTracks());
1643 printf ("GetCTSTracks() %p\n",GetCTSTracks() );
1644 }
1645
477d6cee 1646 if(!fUseJetRefTracks)
be518ab0 1647 ntracks =GetCTSTracks()->GetEntriesFast();
477d6cee 1648 else //If you want to use jet tracks from JETAN
1649 ntracks = (jet->GetRefTracks())->GetEntriesFast();
1650
290c1e1b 1651 if(GetDebug()>3) printf ("ntracks %d\n",ntracks);
88f9563f 1652 AliVTrack* track = 0x0 ;
477d6cee 1653 for(Int_t ipr = 0;ipr < ntracks ; ipr ++ ){
1654 if(!fUseJetRefTracks)
290c1e1b 1655 track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
477d6cee 1656 else //If you want to use jet tracks from JETAN
88f9563f 1657 track = (AliVTrack *) ((jet->GetRefTracks())->At(ipr));
477d6cee 1658
1659 p3.SetXYZ(track->Px(),track->Py(),track->Pz());
1660 pt = p3.Pt();
1661 eta = p3.Eta();
1662 phi = p3.Phi() ;
1663 if(phi < 0) phi+=TMath::TwoPi();
1664
1665 //Check if there is any particle inside cone with pt larger than fPtThreshold
1666 rad = TMath::Sqrt((eta-etaJet)*(eta-etaJet)+ (phi-phiJet)*(phi-phiJet));
290c1e1b 1667 if(rad < fConeSize && pt > fPtThresholdInCone){
f37fa8d2 1668 //printf("charged in jet cone pt %f, phi %f, eta %f, R %f \n",pt,phi,eta,rad);
1669 npartcone++;
1670 fhFFz ->Fill(ptTrig, pt/ptTrig);
1671 fhFFxi->Fill(ptTrig, TMath::Log(ptTrig/pt));
1672 fhFFpt->Fill(ptTrig, pt);
290c1e1b 1673
1674 //according to jet axis
1675 fhJetFFz ->Fill(ptJet, pt/ptJet);
1676 fhJetFFxi->Fill(ptJet, TMath::Log(ptJet/pt));
1677 fhJetFFpt->Fill(ptJet, pt);
1678
1679
1680 if(TMath::Cos(angleJetGam)<0 && ptJet!=0 && pt!=0 ){
1681 fhJetFFzCor ->Fill(ptJet, -pt*TMath::Cos(angleJetGam)/ptJet);
1682 fhJetFFxiCor->Fill(ptJet, TMath::Log(ptJet/(-pt*TMath::Cos(angleJetGam))));
1683 }
477d6cee 1684 }
1685 }//Tracks
1686 fhNTracksInCone->Fill(ptTrig, npartcone);
290c1e1b 1687 //fill tree here for each photon-jet (isolation required usually)
1688
1689 fGamPt = ptTrig;
1690 //fGamLambda0 = ;//filled earlier
1691 fGamNLM = particlecorr->GetFiducialArea();
1692 //fGamSumPtCh = ;//filled earlier
1693 //fGamTime = particlecorr->GetTOF();//filled earlier
1694 //fGamNcells = particlecorr->GetNCells();//filled earlier
1695 fGamEta = etaTrig;
1696 fGamPhi = phiTrig;
1697 //fGamSumPtNeu = ;//filled earlier
1698 //fGamNtracks = ;//filled earlier
1699 //fGamNclusters= ;//filled earlier
1700 //fGamAvEne = ;//filled earlier
1701 fJetPhi = phiJet;
1702 fJetEta = etaJet;
1703 fJetPt = ptJet;
1704 fJetBkgChEne = (IsBackgroundJetFromReader()?rhoEvent * jet->EffectiveAreaCharged():jet->ChargedBgEnergy());
1705 fJetArea = jet->EffectiveAreaCharged();
1706 fJetNtracks = (jet->GetRefTracks())->GetEntriesFast();
1707 fEventNumber = 0;
1708 fNtracks = GetCTSTracks()->GetEntriesFast();
1709 fCentrality = GetEventCentrality();
1710 fIso = particlecorr->IsIsolated();
1711
1712 Int_t nTrk1GeV=0;
1713 Int_t nTrk2GeV=0;
1714 for(itrack=0;itrack < fJetNtracks;itrack++){
1715 track = (AliVTrack *) ((jet->GetRefTracks())->At(itrack));
1716 if(track->Pt()>1.) nTrk1GeV++;
1717 if(track->Pt()>2.) nTrk2GeV++;
1718 }
1719
1720 fJetNtracks1 = nTrk1GeV;
1721 fJetNtracks2 = nTrk2GeV;
1722
1723 if(fSaveGJTree) fTreeGJ->Fill();
477d6cee 1724 }//AOD trigger particle loop
a3aebfff 1725 if(GetDebug() > 1) printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - End fill histograms \n");
290c1e1b 1726}
1c5acb87 1727
1728
1729//__________________________________________________________________
1730void AliAnaParticleJetFinderCorrelation::Print(const Option_t * opt) const
1731{
1c5acb87 1732
477d6cee 1733 //Print some relevant parameters set for the analysis
1734 if(! opt)
1735 return;
1736
a3aebfff 1737 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
745913ae 1738 AliAnaCaloTrackCorrBaseClass::Print(" ");
a3aebfff 1739
477d6cee 1740 printf("Phi trigger-jet < %3.2f\n", fDeltaPhiMaxCut) ;
1741 printf("Phi trigger-jet > %3.2f\n", fDeltaPhiMinCut) ;
1742 printf("pT Ratio trigger/jet < %3.2f\n", fRatioMaxCut) ;
1743 printf("pT Ratio trigger/jet > %3.2f\n", fRatioMinCut) ;
1744 printf("fConeSize = %3.2f\n", fConeSize) ;
1745 printf("fPtThresholdInCone = %3.2f\n", fPtThresholdInCone) ;
1746 printf("fUseJetRefTracks = %d\n", fUseJetRefTracks) ;
1747 printf("fMakeCorrelationInHistoMaker = %d\n", fMakeCorrelationInHistoMaker) ;
1748 printf("Isolated Trigger? %d\n", fSelectIsolated) ;
290c1e1b 1749 printf("Reconstructed jet cone size = %3.2f\n", fJetConeSize) ;
1750 printf("Reconstructed jet minimum pt = %3.2f\n", fJetMinPt) ;
1751 printf("Reconstructed jet minimum area fraction = %3.2f\n", fJetAreaFraction) ;
1752
9717eff6 1753// if(!fNonStandardJetFromReader){
290c1e1b 1754 printf("fJetBranchName = %s\n", fJetBranchName.Data()) ;
9717eff6 1755// }
290c1e1b 1756 if(!fBackgroundJetFromReader){
1757 printf("fBkgJetBranchName = %s\n", fBkgJetBranchName.Data()) ;
1758 }
1759
1760 printf("Isolation cone size = %3.2f\n", fGammaConeSize) ;
1761 printf("fUseBackgroundSubtractionGamma = %d\n",fUseBackgroundSubtractionGamma);
1762 printf("fSaveGJTree = %d\n",fSaveGJTree);
1763 printf("fMostEnergetic = %d\n",fMostEnergetic);
1764 printf("fMostOpposite = %d\n",fMostOpposite);
1765
1766 printf("fUseHistogramJetBkg = %d\n",fUseHistogramJetBkg);
1767 printf("fUseHistogramTracks = %d\n",fUseHistogramTracks);
1768 printf("fUseHistogramJetTracks = %d\n",fUseHistogramJetTracks);
1769
1c5acb87 1770}
1771
290c1e1b 1772//__________________________________________________________________
1773void AliAnaParticleJetFinderCorrelation::CalculateBkg(TVector3 gamma, TVector3 jet,Double_t vertex[3],Int_t type=2){
1774 //
1775 // calculate background for fragmentation function and fill histograms
1776 // 1. 90 degrees from jet axis in random place = perpendicular cone
1777 // 2. Random cone not belonging to jet cone nor photon cone
1778 // 3. In the middle point from jet and photon momentum vectors
1779 // 4. 90 degrees from photon direction in random place = perpendicular cone 2
1780
1781 //
1782 // implementation of 2 works, 1 and 4 works
1783 //
1784 Double_t gammaPt = gamma.Pt();
1785 Double_t gammaEta = gamma.Eta();
1786 Double_t gammaPhi = gamma.Phi();
1787 Double_t jetEta = jet.Eta();
1788 Double_t jetPhi = jet.Phi();
1789
1790 //refference direction of background
1791 Double_t refEta=0.,refPhi=0.;
1792 Double_t rad = 0,rad2 = 0.;
1793 if(type==1){//perpendicular to jet axis
1794 //printf("vertex: %f %f %f \n",vertex[0],vertex[1],vertex[2]);
1795
1796 Double_t xVar;
1797 Double_t yVar;
1798 Double_t newX=0.;
1799 Double_t newY=0.;
1800 Double_t newZ=0.;
1801 //Z axis vector: [jet.Px(),jet.Py(),jet.Pz()]
1802 Double_t jx=jet.Px();
1803 Double_t jy=jet.Py();
1804 Double_t jz=jet.Pz();
1805 //if(jz==0)printf("problem\n");
1806 //X axis
1807 Double_t Xx=1.0-vertex[0];
1808 Double_t Xy=-1.0*vertex[1];
1809 Double_t Xz=jx/jz*(vertex[0]-1.)+vertex[1]*jy/jz;
1810 //Y axis
1811 Double_t Yx=jy*Xz-jz*Xy;
1812 Double_t Yy=jz*Xx-jx*Xz;
1813 Double_t Yz=jx*Xy-jy*Xx;
1814 //Determinant
1815 Double_t det = Xx*Yy*jz + Xy*Yz*jx + Xz*Yx*jy - Xx*Yz*jy - Xy*Yx*jz - Xz*Yy*jx;
1816 if(det==0)printf("problem det==0\n");
1817 Double_t detX = 0.;
1818 Double_t detY = 0.;
1819 Double_t detZ = 0.;
1820
1821// Double_t tmpScalar = jx*Xx+jy*Xy+jz*Xz;
1822// printf("scalar jet.P o X: %f\n",tmpScalar);
1823// tmpScalar = jet.Px()*Yx+jet.Py()*Yy+jet.Pz()*Yz;
1824// printf("scalar jet.P o Y: %f\n",tmpScalar);
1825// tmpScalar = Xx*Yx+Xy*Yy+Xz*Yz;
1826// printf("scalar X o Y: %f\n",tmpScalar);
1827
1828 TVector3 perp;
1829 //randomise
1830 do{
1831 refPhi=fGenerator->Rndm()*TMath::Pi()*2.;
1832 //refPhi=fGenerator->Uniform(-1,1)*TMath::Pi();
1833 xVar=TMath::Cos(refPhi);
1834 yVar=TMath::Sin(refPhi);
1835 //yVar=TMath::Cos(-TMath::Pi()/2.+refPhi);
1836 //zVar=0 in new surface frame
1837 detX = xVar*Yy*jz + Xz*yVar*jy - xVar*Yz*jy - Xy*yVar*jz;
1838 detY = Xx*yVar*jz + xVar*Yz*jx - xVar*Yx*jz - Xz*yVar*jx;
1839 detZ = Xy*yVar*jx + xVar*Yx*jy - Xx*yVar*jy - xVar*Yy*jx;
1840
1841 newX=detX/det;
1842 newY=detY/det;
1843 newZ=detZ/det;
1844
1845 perp.SetXYZ(newX,newY,newZ);
1846 refEta = perp.Eta();
1847 refPhi = perp.Phi();//output in <-pi, pi> range
1848 if(refPhi<0)refPhi+=2*TMath::Pi();
1849 rad = TMath::Sqrt((gammaEta-refEta)*(gammaEta-refEta) + (gammaPhi-refPhi)*(gammaPhi-refPhi));
1850 rad2 = TMath::Sqrt((jetEta-refEta)*(jetEta-refEta) + (jetPhi-refPhi)*(jetPhi-refPhi));
1851 //printf("refEta,refPhi,rad,rad2: %f %f %f %f\n",refEta,refPhi,rad,rad2);
9717eff6 1852 } while (rad<fJetConeSize+fGammaConeSize || rad2<2.*fJetConeSize || TMath::Abs(refEta)>0.9-fJetConeSize);
290c1e1b 1853 fhRandomPhiEta[2]->Fill(refPhi,refEta);
1854
1855 }
1856 else if(type==2){//random cone
1857 //randomise
1858 do{
1859 refPhi=fGenerator->Rndm()*TMath::Pi()*2.;
1860 refEta=fGenerator->Uniform(-(0.9-fJetConeSize),0.9-fJetConeSize);
1861 rad = TMath::Sqrt((gammaEta-refEta)*(gammaEta-refEta) + (gammaPhi-refPhi)*(gammaPhi-refPhi));
1862 rad2 = TMath::Sqrt((jetEta-refEta)*(jetEta-refEta) + (jetPhi-refPhi)*(jetPhi-refPhi));
1863 //check if reference is not within jet cone or gamma cone +0.4
1864 //example: FF fConeSize=0.4, fJetConeSize=0.4, fIsoGammaConeSize=0.3
9717eff6 1865 } while (rad<fJetConeSize+fGammaConeSize || rad2<2.*fJetConeSize);
290c1e1b 1866 //photon:0.7=0.4+0.3; jets:0.8=0.4 +0.4 //rad<fConeSize+fJetConeSize rad2<fConeSize+0.3
1867 fhRandomPhiEta[0]->Fill(refPhi,refEta);
1868 }
1869 else if(type==4){//perpendicular to photon axis
1870 Double_t xVar;
1871 Double_t yVar;
1872 Double_t newX=0.;
1873 Double_t newY=0.;
1874 Double_t newZ=0.;
1875 //Z axis vector: [jet.Px(),jet.Py(),jet.Pz()]
1876 Double_t jx=gamma.Px();
1877 Double_t jy=gamma.Py();
1878 Double_t jz=gamma.Pz();
1879 //if(jz==0)printf("problem\n");
1880 //X axis
1881 Double_t Xx=1.0-vertex[0];
1882 Double_t Xy=-1.0*vertex[1];
1883 Double_t Xz=jx/jz*(vertex[0]-1.)+vertex[1]*jy/jz;
1884 //Y axis
1885 Double_t Yx=jy*Xz-jz*Xy;
1886 Double_t Yy=jz*Xx-jx*Xz;
1887 Double_t Yz=jx*Xy-jy*Xx;
1888 //Determinant
1889 Double_t det = Xx*Yy*jz + Xy*Yz*jx + Xz*Yx*jy - Xx*Yz*jy - Xy*Yx*jz - Xz*Yy*jx;
1890 if(det==0)printf("problem det==0\n");
1891 Double_t detX = 0.;
1892 Double_t detY = 0.;
1893 Double_t detZ = 0.;
1894
1895// Double_t tmpScalar = jx*Xx+jy*Xy+jz*Xz;
1896// printf("scalar jet.P o X: %f\n",tmpScalar);
1897// tmpScalar = jet.Px()*Yx+jet.Py()*Yy+jet.Pz()*Yz;
1898// printf("scalar jet.P o Y: %f\n",tmpScalar);
1899// tmpScalar = Xx*Yx+Xy*Yy+Xz*Yz;
1900// printf("scalar X o Y: %f\n",tmpScalar);
1901
1902 TVector3 perp;
1903 //randomise
1904 do{
1905 refPhi=fGenerator->Rndm()*TMath::Pi()*2.;
1906 //refPhi=fGenerator->Uniform(-1,1)*TMath::Pi();
1907 xVar=TMath::Cos(refPhi);
1908 yVar=TMath::Sin(refPhi);
1909 //yVar=TMath::Cos(-TMath::Pi()/2.+refPhi);
1910 //zVar=0 in new surface frame
1911 detX = xVar*Yy*jz + Xz*yVar*jy - xVar*Yz*jy - Xy*yVar*jz;
1912 detY = Xx*yVar*jz + xVar*Yz*jx - xVar*Yx*jz - Xz*yVar*jx;
1913 detZ = Xy*yVar*jx + xVar*Yx*jy - Xx*yVar*jy - xVar*Yy*jx;
1914
1915 newX=detX/det;
1916 newY=detY/det;
1917 newZ=detZ/det;
1918
1919 perp.SetXYZ(newX,newY,newZ);
1920 refEta = perp.Eta();
1921 refPhi = perp.Phi();//output in <-pi, pi> range
1922 if(refPhi<0)refPhi+=2*TMath::Pi();
1923 rad = TMath::Sqrt((gammaEta-refEta)*(gammaEta-refEta) + (gammaPhi-refPhi)*(gammaPhi-refPhi));
1924 rad2 = TMath::Sqrt((jetEta-refEta)*(jetEta-refEta) + (jetPhi-refPhi)*(jetPhi-refPhi));
1925 //printf("refEta,refPhi,rad,rad2: %f %f %f %f\n",refEta,refPhi,rad,rad2);
9717eff6 1926 } while (rad<fJetConeSize+fGammaConeSize || rad2<2.*fJetConeSize || TMath::Abs(refEta)>0.9-fJetConeSize);
290c1e1b 1927 fhRandomPhiEta[1]->Fill(refPhi,refEta);
1928
1929 }
1930 else if(type==3){//mid point
1931
1932 Double_t jx=jet.Px();
1933 Double_t jy=jet.Py();
1934 Double_t jz=jet.Pz();
1935 // if(jz==0)printf("problem\n");
1936 Double_t gx=gamma.Px();
1937 Double_t gy=gamma.Py();
1938 Double_t gz=gamma.Pz();
1939
1940 Double_t cosAlpha=(jx*gx+jy*gy+jz*gz)/(jet.Mag()*gamma.Mag());
1941 Double_t cosinus=TMath::Sqrt((cosAlpha+1.)/2.);
1942 //perpendicular axis
1943 Double_t Zx=gy*jz-gz*jy;
1944 Double_t Zy=gz*jx-gx*jz;
1945 Double_t Zz=gx*jy-gy*jx;
1946
1947 //Determinant
1948 Double_t det = Zx*gy*jz + Zy*gz*jx + Zz*gx*jy - Zz*gy*jx - Zy*gx*jz - Zx*gz*jy;
1949
1950 Double_t newX=0.;
1951 Double_t newY=0.;
1952 Double_t newZ=0.;
1953 if(det!=0) {
1954 Double_t detX = -Zy*gz*cosinus +Zz*cosinus*jy + Zz*gy*cosinus - Zy*cosinus*jz;
1955 Double_t detY = Zx*cosinus*jz - Zz*gx*cosinus - Zz*cosinus*jx + Zx*gz*cosinus;
1956 Double_t detZ = -Zx*gy*cosinus + Zy*cosinus*jx + Zy*gx*cosinus - Zx*cosinus*jy;
1957
1958 newX=detX/det;
1959 newY=detY/det;
1960 newZ=detZ/det;
1961 }
1962
1963 TVector3 perp;
1964 perp.SetXYZ(newX,newY,newZ);
1965 refEta = perp.Eta();
1966 refPhi = perp.Phi();//output in <-pi, pi> range
1967 if(refPhi<0)refPhi+=2*TMath::Pi();
1968 rad = TMath::Sqrt((gammaEta-refEta)*(gammaEta-refEta) + (gammaPhi-refPhi)*(gammaPhi-refPhi));
1969 rad2 = TMath::Sqrt((jetEta-refEta)*(jetEta-refEta) + (jetPhi-refPhi)*(jetPhi-refPhi));
1970 //printf("refEta,refPhi,rad,rad2: %f %f %f %f\n",refEta,refPhi,rad,rad2);
1971
9717eff6 1972 if (rad<fJetConeSize+fGammaConeSize || rad2<2.*fJetConeSize || TMath::Abs(refEta)>0.9-fJetConeSize) fhRandomPhiEta[3]->Fill(refPhi,refEta);
290c1e1b 1973 }
1974 else if(type==5){//tmp
1975 //printf("vertex: %f %f %f \n",vertex[0],vertex[1],vertex[2]);
1976
1977 Double_t xVar;
1978 Double_t newX=0.;
1979 Double_t newY=0.;
1980 Double_t newZ=0.;
1981 //Z axis vector: [jet.Px(),jet.Py(),jet.Pz()]
1982 Double_t jx=jet.Px();
1983 Double_t jy=jet.Py();
1984 Double_t jz=jet.Pz();
1985 // if(jz==0)printf("problem\n");
1986 //X axis
1987 Double_t Xx=1.0-vertex[0];
1988 Double_t Xy=-1.0*vertex[1];
1989 Double_t Xz=jx/jz*(vertex[0]-1.)+vertex[1]*jy/jz;
1990 //Y axis
1991 Double_t Yx=jy*Xz-jz*Xy;
1992 Double_t Yy=jz*Xx-jx*Xz;
1993 Double_t Yz=jx*Xy-jy*Xx;
1994
1995 // X and Y length
1996 Double_t Xlength=TMath::Sqrt(Xx*Xx+Xy*Xy+Xz*Xz);
1997 Double_t Ylength=TMath::Sqrt(Yx*Yx+Yy*Yy+Yz*Yz);
1998 Double_t ratio=Ylength/Xlength;
1999
2000 TVector3 perp;
2001 //randomise
2002 do{
2003 refPhi=fGenerator->Rndm()*TMath::Pi()*2.;
2004 xVar=TMath::Tan(refPhi)/ratio;
2005 newX=xVar*Yx+Xx;
2006 newY=xVar*Yy+Xy;
2007 newZ=xVar*Yz+Xz;
2008
2009 perp.SetXYZ(newX,newY,newZ);
2010 refEta = perp.Eta();
2011 refPhi = perp.Phi();//output in <-pi, pi> range
2012 if(refPhi<0)refPhi+=2*TMath::Pi();
2013 rad = TMath::Sqrt((gammaEta-refEta)*(gammaEta-refEta) + (gammaPhi-refPhi)*(gammaPhi-refPhi));
2014 rad2 = TMath::Sqrt((jetEta-refEta)*(jetEta-refEta) + (jetPhi-refPhi)*(jetPhi-refPhi));
2015 //printf("refEta,refPhi,rad,rad2: %f %f %f %f\n",refEta,refPhi,rad,rad2);
9717eff6 2016 } while (rad<fJetConeSize+fGammaConeSize || rad2<2.*fJetConeSize || TMath::Abs(refEta)>0.9-fJetConeSize);
290c1e1b 2017 fhRandomPhiEta[4]->Fill(refPhi,refEta);
2018 }
2019
2020
2021
2022 //calculate FF in background
2023 Int_t ntracks = 0;
2024 ntracks =GetCTSTracks()->GetEntriesFast();
2025 AliVTrack* track = 0x0 ;
2026 TVector3 p3;
2027
2028 Double_t pt = 0, eta = 0, phi = 0;
2029 Int_t npartcone = 0;
2030 Double_t sumPt=0.;
2031 for(Int_t ipr = 0;ipr < ntracks ; ipr ++ ){
2032 track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
2033 p3.SetXYZ(track->Px(),track->Py(),track->Pz());
2034 pt = p3.Pt();
2035 if(pt<fPtThresholdInCone) {//0.150
2036 //printf("problem: track pt < %f MeV/c \n",fPtThresholdInCone);
2037 continue;
2038 }
2039 eta = p3.Eta() ;
2040 phi = p3.Phi() ;
2041 if(phi < 0) phi+=TMath::TwoPi();
2042 //Check if there is any particle inside cone with pt larger than fPtThreshold
2043 rad = TMath::Sqrt((eta-refEta)*(eta-refEta) + (phi-refPhi)*(phi-refPhi));
2044 if(rad < fConeSize && pt > fPtThresholdInCone){
2045 //printf("charged in jet cone pt %f, phi %f, eta %f, R %f \n",pt,phi,eta,rad);
2046 npartcone++;
2047 sumPt+=pt;
2048 if(type==1){//perp jet
2049 fhBkgFFz[1] ->Fill(gammaPt, pt/gammaPt);
2050 fhBkgFFxi[1]->Fill(gammaPt, TMath::Log(gammaPt/pt));
2051 fhBkgFFpt[1]->Fill(gammaPt, pt);
2052 }
2053 else if(type==2){//RC
2054 fhBkgFFz[0] ->Fill(gammaPt, pt/gammaPt);
2055 fhBkgFFxi[0]->Fill(gammaPt, TMath::Log(gammaPt/pt));
2056 fhBkgFFpt[0]->Fill(gammaPt, pt);
2057 }
2058 else if(type==3){//mid point
2059 fhBkgFFz[3] ->Fill(gammaPt, pt/gammaPt);
2060 fhBkgFFxi[3]->Fill(gammaPt, TMath::Log(gammaPt/pt));
2061 fhBkgFFpt[3]->Fill(gammaPt, pt);
2062 }
2063 else if(type==4){//perp jet
2064 fhBkgFFz[2] ->Fill(gammaPt, pt/gammaPt);
2065 fhBkgFFxi[2]->Fill(gammaPt, TMath::Log(gammaPt/pt));
2066 fhBkgFFpt[2]->Fill(gammaPt, pt);
2067 }
2068 else if(type==5){//test
2069 fhBkgFFz[4] ->Fill(gammaPt, pt/gammaPt);
2070 fhBkgFFxi[4]->Fill(gammaPt, TMath::Log(gammaPt/pt));
2071 fhBkgFFpt[4]->Fill(gammaPt, pt);
2072 }
2073
2074
2075 }
2076 }//end of loop over tracks
2077 Double_t sumOverTracks=0.;
2078 if(npartcone!=0) sumOverTracks = sumPt/npartcone;
2079 if(type==1) {
2080 fhBkgNTracksInCone[1]->Fill(gammaPt, npartcone);
2081 fhBkgSumPtInCone[1]->Fill(gammaPt,sumPt);
2082 fhBkgSumPtnTracksInCone[1]->Fill(gammaPt,sumOverTracks);
2083 }
2084 else if(type==2) {
2085 fhBkgNTracksInCone[0]->Fill(gammaPt, npartcone);
2086 fhBkgSumPtInCone[0]->Fill(gammaPt,sumPt);
2087 fhBkgSumPtnTracksInCone[0]->Fill(gammaPt,sumOverTracks);
2088 }
2089 else if(type==3) {
2090 fhBkgNTracksInCone[3]->Fill(gammaPt, npartcone);
2091 fhBkgSumPtInCone[3]->Fill(gammaPt,sumPt);
2092 fhBkgSumPtnTracksInCone[3]->Fill(gammaPt,sumOverTracks);
2093 }
2094 else if(type==4) {
2095 fhBkgNTracksInCone[2]->Fill(gammaPt, npartcone);
2096 fhBkgSumPtInCone[2]->Fill(gammaPt,sumPt);
2097 fhBkgSumPtnTracksInCone[2]->Fill(gammaPt,sumOverTracks);
2098 }
2099 else if(type==5) {
2100 fhBkgNTracksInCone[4]->Fill(gammaPt, npartcone);
2101 fhBkgSumPtInCone[4]->Fill(gammaPt,sumPt);
2102 fhBkgSumPtnTracksInCone[4]->Fill(gammaPt,sumOverTracks);
2103 }
2104}
2105
2106