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