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