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