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