1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //_________________________________________________________________________
17 // Class for the analysis of particle (direct gamma) -jet (jet found with finder) correlations
18 //*-- Author: Gustavo Conesa (LNF-INFN)
19 //*-- Modified: Adam Matyja (INP PAN Cracow)
20 //////////////////////////////////////////////////////////////////////////////
23 // --- ROOT system ---
25 #include "TClonesArray.h"
27 //#include "Riostream.h"
29 //---- AliRoot system ----
30 #include "AliCaloTrackReader.h"
31 #include "AliAODJet.h"
32 #include "AliAnaParticleJetFinderCorrelation.h"
33 #include "AliAODPWG4ParticleCorrelation.h"
34 #include "AliVTrack.h"
35 #include "AliAODCaloCluster.h"
36 #include "AliAODEvent.h"
39 #include "AliAODJetEventBackground.h"
43 #include "AliMCAnalysisUtils.h"
44 #include "AliAODMCParticle.h"
46 #include "AliEMCALGeometry.h"
48 ClassImp(AliAnaParticleJetFinderCorrelation)
50 //________________________________________________________________________
51 AliAnaParticleJetFinderCorrelation::AliAnaParticleJetFinderCorrelation() :
52 AliAnaCaloTrackCorrBaseClass(),
53 fDeltaPhiMaxCut(0.), fDeltaPhiMinCut(0.), fRatioMaxCut(0.), fRatioMinCut(0.),
54 fConeSize(0.), fPtThresholdInCone(0.),fUseJetRefTracks(kTRUE),
55 fMakeCorrelationInHistoMaker(kFALSE), fSelectIsolated(kTRUE),
56 fJetConeSize(0.4),fJetMinPt(5),fJetAreaFraction(0.6),
57 //fNonStandardJetFromReader(kTRUE),
58 fJetBranchName("jets"),
59 fBackgroundJetFromReader(kTRUE),
60 fBkgJetBranchName("jets"),
61 fGammaConeSize(0.3),fUseBackgroundSubtractionGamma(kFALSE),fSaveGJTree(kTRUE),
62 fMostEnergetic(kFALSE),fMostOpposite(kTRUE), fUseHistogramJetBkg(kTRUE),
63 fUseHistogramTracks(kTRUE),fUseHistogramJetTracks(kTRUE),fMCStudies(kFALSE),fGenerator(0),
64 fhDeltaEta(0), /*fhDeltaPhi(0),*/fhDeltaPhiCorrect(0),fhDeltaPhi0PiCorrect(0), fhDeltaPt(0), fhPtRatio(0), fhPt(0),
65 fhFFz(0),fhFFxi(0),fhFFpt(0),fhNTracksInCone(0),
66 fhJetFFz(0),fhJetFFxi(0),fhJetFFpt(0),fhJetFFzCor(0),fhJetFFxiCor(0),
67 fhBkgFFz(),fhBkgFFxi(),fhBkgFFpt(),fhBkgNTracksInCone(),fhBkgSumPtInCone(),fhBkgSumPtnTracksInCone(),//<<---new
68 fhNjetsNgammas(0),fhCuts(0),
69 fhDeltaEtaBefore(0),fhDeltaPhiBefore(0),fhDeltaPtBefore(0),fhPtRatioBefore(0),
70 fhPtBefore(0),fhDeltaPhi0PiCorrectBefore(0),
71 fhJetPtBefore(0),fhJetPt(0),fhJetPtMostEne(0),fhJetPhi(0),fhJetEta(0),fhJetEtaVsPt(0),
72 fhJetPhiVsEta(0),fhJetEtaVsNpartInJet(0),fhJetEtaVsNpartInJetBkg(0),fhJetChBkgEnergyVsPt(0),fhJetChAreaVsPt(0),/*fhJetNjet(0),*/
73 fhTrackPhiVsEta(0),fhTrackAveTrackPt(0),fhJetNjetOverPtCut(),
74 /*fhJetChBkgEnergyVsPtEtaGt05(0),fhJetChBkgEnergyVsPtEtaLe05(0),fhJetChAreaVsPtEtaGt05(0),fhJetChAreaVsPtEtaLe05(0),*/
75 fhJetChBkgEnergyVsArea(0),fhJetRhoVsPt(0),fhJetRhoVsCentrality(0),//fhJetBkgRho(0),
76 fhJetNparticlesInJet(0),fhJetDeltaEtaDeltaPhi(0),fhJetDeltaEtaDeltaPhiAllTracks(0),
77 fhJetAveTrackPt(0),fhJetNtracksInJetAboveThr(),fhJetRatioNTrkAboveToNTrk(),fhJetNtrackRatioMostEne(),
78 fhJetNtrackRatioJet5GeV(),fhJetNtrackRatioLead5GeV(),
79 fhBkgJetBackground(),fhBkgJetSigma(),fhBkgJetArea(),fhPhotonPtMostEne(0),
80 fhPhotonAverageEnergy(0),fhPhotonRatioAveEneToMostEne(0),fhPhotonAverageEnergyMinus1(0),fhPhotonRatioAveEneMinus1ToMostEne(0),
81 fhPhotonNgammaMoreAverageToNgamma(0),fhPhotonNgammaMoreAverageMinus1ToNgamma(0),fhPhotonNgammaOverPtCut(),
82 fhPhotonBkgRhoVsNtracks(0),fhPhotonBkgRhoVsNclusters(0),fhPhotonBkgRhoVsCentrality(0),
83 fhPhotonBkgRhoVsNcells(0),fhPhotonPt(0),fhPhotonPtCorrected(0),fhPhotonPtCorrectedZoom(0),fhPhotonPtDiff(0),
84 fhPhotonPtDiffVsCentrality(0),fhPhotonPtDiffVsNcells(0),fhPhotonPtDiffVsNtracks(0),fhPhotonPtDiffVsNclusters(0),
85 fhPhotonSumPtInCone(0),fhPhotonSumPtCorrectInCone(0),fhPhotonSumPtChargedInCone(0),
86 fhSelectedJetPhiVsEta(0),fhSelectedJetChBkgEnergyVsPtJet(0),fhSelectedJetChAreaVsPtJet(0),fhSelectedJetNjet(0),fhSelectedNtracks(0),
87 fhSelectedTrackPhiVsEta(0),fhCuts2(0),
88 fhSelectedPhotonNLMVsPt(0),fhSelectedPhotonLambda0VsPt(0), fhRandomPhiEta(),
89 fhMCPhotonCuts(0),fhMCPhotonPt(0),fhMCPhotonEtaPhi(0),fhMCJetOrigin(0),
90 fhMCJetNPartVsPt(0),fhMCJetChNPartVsPt(0),fhMCJetNPart150VsPt(0),fhMCJetChNPart150VsPt(0),fhMCJetChNPart150ConeVsPt(0),
91 fhMCJetRatioChFull(0),fhMCJetRatioCh150Ch(0),
92 fhMCJetEtaPhi(0),fhMCJetChEtaPhi(0),fhMCJet150EtaPhi(0),fhMCJetCh150EtaPhi(0),fhMCJetCh150ConeEtaPhi(0),
141 fMCJetCh150ConePt(0),
142 fMCJetCh150ConeNPart(0),
143 fMCJetCh150ConeEta(0),
144 fMCJetCh150ConePhi(0)
148 //printf("constructor\n");
150 //Initialize parameters
152 for(Int_t i=0;i<10;i++){
153 fhJetNjetOverPtCut[i] = 0;
154 fhPhotonNgammaOverPtCut[i] = 0;
156 fGenerator = new TRandom2();
157 fGenerator->SetSeed(0);
160 //___________________________________________________________________
161 AliAnaParticleJetFinderCorrelation::~AliAnaParticleJetFinderCorrelation(){
166 //___________________________________________________________________
167 TList * AliAnaParticleJetFinderCorrelation::GetCreateOutputObjects()
169 // Create histograms to be saved in output file and
170 // store them in fOutputContainer
171 //printf("GetCreateOutputObjects\n");
173 TList * outputContainer = new TList() ;
174 outputContainer->SetName("ParticleJetFinderHistos") ;
176 Int_t nptbins = GetHistogramRanges()->GetHistoPtBins();
177 // Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins();
178 // Int_t netabins = GetHistogramRanges()->GetHistoEtaBins();
179 Float_t ptmax = GetHistogramRanges()->GetHistoPtMax();
180 // Float_t phimax = GetHistogramRanges()->GetHistoPhiMax();
181 // Float_t etamax = GetHistogramRanges()->GetHistoEtaMax();
182 Float_t ptmin = GetHistogramRanges()->GetHistoPtMin();
183 // Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
184 // Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
186 // fhDeltaPhi = new TH2F("DeltaPhi","#phi_{trigger} - #phi_{jet} vs p_{T trigger}",nptbins,ptmin,ptmax,100,-6,4);
187 // fhDeltaPhi->SetYTitle("#Delta #phi");
188 // fhDeltaPhi->SetXTitle("p_{T trigger} (GeV/c)");
189 // outputContainer->Add(fhDeltaPhi);
191 fhDeltaPhiCorrect = new TH2F("DeltaPhiCorrect","#phi_{trigger} - #phi_{jet} vs p_{T trigger}",nptbins,ptmin,ptmax,100,0,6.5);
192 fhDeltaPhiCorrect->SetYTitle("#Delta #phi");
193 fhDeltaPhiCorrect->SetXTitle("p_{T trigger} (GeV/c)");
194 outputContainer->Add(fhDeltaPhiCorrect);
196 fhDeltaPhi0PiCorrect = new TH2F("DeltaPhi0PiCorrect","#phi_{trigger} - #phi_{jet} (0,#pi) vs p_{T trigger}",nptbins,ptmin,ptmax,100,0,3.5);
197 fhDeltaPhi0PiCorrect->SetYTitle("#Delta #phi");
198 fhDeltaPhi0PiCorrect->SetXTitle("p_{T trigger} (GeV/c)");
199 outputContainer->Add(fhDeltaPhi0PiCorrect);
202 fhDeltaEta = new TH2F("DeltaEta","#eta_{trigger} - #eta_{jet} vs p_{T trigger}",nptbins,ptmin,ptmax,100,-2,2);
203 fhDeltaEta->SetYTitle("#Delta #eta");
204 fhDeltaEta->SetXTitle("p_{T trigger} (GeV/c)");
205 outputContainer->Add(fhDeltaEta);
207 fhDeltaPt = new TH2F("DeltaPt","p_{T trigger} - p_{T jet} vs p_{T trigger}",nptbins,ptmin,ptmax,150,-50,100);
208 fhDeltaPt->SetYTitle("#Delta p_{T}");
209 fhDeltaPt->SetXTitle("p_{T trigger} (GeV/c)");
210 outputContainer->Add(fhDeltaPt);
212 fhPtRatio = new TH2F("PtRatio","p_{T jet} / p_{T trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,200,0,2.);
213 fhPtRatio->SetYTitle("ratio");
214 fhPtRatio->SetXTitle("p_{T trigger} (GeV/c)");
215 outputContainer->Add(fhPtRatio);
217 fhPt = new TH2F("Pt","p_{T jet} vs p_{T trigger}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
218 fhPt->SetYTitle("p_{T jet}(GeV/c)");
219 fhPt->SetXTitle("p_{T trigger} (GeV/c)");
220 outputContainer->Add(fhPt);
222 fhFFz = new TH2F("FFz","z = p_{T i charged}/p_{T trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,200,0.,2);
223 fhFFz->SetYTitle("z");
224 fhFFz->SetXTitle("p_{T trigger}");
225 outputContainer->Add(fhFFz) ;
227 fhFFxi = new TH2F("FFxi","#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger}", nptbins,ptmin,ptmax,100,0.,10.);
228 fhFFxi->SetYTitle("#xi");
229 fhFFxi->SetXTitle("p_{T trigger}");
230 outputContainer->Add(fhFFxi) ;
232 fhFFpt = new TH2F("FFpt","p_{T i charged} vs p_{T trigger}", nptbins,ptmin,ptmax,100,0.,50.);
233 fhFFpt->SetYTitle("p_{T charged hadron}");
234 fhFFpt->SetXTitle("p_{T trigger}");
235 outputContainer->Add(fhFFpt) ;
237 fhNTracksInCone = new TH2F("NTracksInCone","Number of tracks in cone vs p_{T trigger}", nptbins,ptmin,ptmax,100,0.,150.);
238 fhNTracksInCone->SetYTitle("p_{T charged hadron}");
239 fhNTracksInCone->SetXTitle("p_{T trigger}");
240 outputContainer->Add(fhNTracksInCone) ;
242 //FF according to jet axis
243 fhJetFFz = new TH2F("JetFFz","z = p_{T i charged}/p_{T jet} vs p_{T jet}",nptbins,ptmin,ptmax,200,0.,2);
244 fhJetFFz->SetYTitle("z");
245 fhJetFFz->SetXTitle("p_{T jet}");
246 outputContainer->Add(fhJetFFz) ;
248 fhJetFFxi = new TH2F("JetFFxi","#xi = ln(p_{T jet}/p_{T i charged}) vs p_{T jet}", nptbins,ptmin,ptmax,100,0.,10.);
249 fhJetFFxi->SetYTitle("#xi");
250 fhJetFFxi->SetXTitle("p_{T jet}");
251 outputContainer->Add(fhJetFFxi) ;
253 fhJetFFpt = new TH2F("JetFFpt","p_{T i charged} vs p_{T jet}", nptbins,ptmin,ptmax,100,0.,50.);
254 fhJetFFpt->SetYTitle("p_{T charged hadron}");
255 fhJetFFpt->SetXTitle("p_{T jet}");
256 outputContainer->Add(fhJetFFpt) ;
258 fhJetFFzCor = new TH2F("JetFFzCor","z = -cos(#alpha(jet,trig))*p_{T i charged}/p_{T jet} vs p_{T jet}",nptbins,ptmin,ptmax,200,0.,2);
259 fhJetFFzCor->SetYTitle("z");
260 fhJetFFzCor->SetXTitle("p_{T jet}");
261 outputContainer->Add(fhJetFFzCor) ;
263 fhJetFFxiCor = new TH2F("JetFFxiCor","#xi = ln(p_{T jet}/(-cos(#alpha(jet,trig))*p_{T i charged})) vs p_{T jet}", nptbins,ptmin,ptmax,100,0.,10.);
264 fhJetFFxiCor->SetYTitle("#xi");
265 fhJetFFxiCor->SetXTitle("p_{T jet}");
266 outputContainer->Add(fhJetFFxiCor) ;
270 fhBkgFFz[0] = new TH2F("BkgFFzRC", "z = p_{T i charged}/p_{T trigger} vs p_{T trigger} Bkg RC" ,nptbins,ptmin,ptmax,200,0.,2);
271 fhBkgFFz[1] = new TH2F("BkgFFzPCG", "z = p_{T i charged}/p_{T trigger} vs p_{T trigger} Bkg PCG" ,nptbins,ptmin,ptmax,200,0.,2);
272 fhBkgFFz[2] = new TH2F("BkgFFzPCJ", "z = p_{T i charged}/p_{T trigger} vs p_{T trigger} Bkg PCJ" ,nptbins,ptmin,ptmax,200,0.,2);
273 fhBkgFFz[3] = new TH2F("BkgFFzMP", "z = p_{T i charged}/p_{T trigger} vs p_{T trigger} Bkg MP" ,nptbins,ptmin,ptmax,200,0.,2);
274 fhBkgFFz[4] = new TH2F("BkgFFzTest","z = p_{T i charged}/p_{T trigger} vs p_{T trigger} Bkg Test",nptbins,ptmin,ptmax,200,0.,2);
275 for(Int_t i=0;i<5;i++){
276 fhBkgFFz[i]->SetYTitle("z");
277 fhBkgFFz[i]->SetXTitle("p_{T trigger}");
278 outputContainer->Add(fhBkgFFz[i]) ;
281 fhBkgFFxi[0] = new TH2F("BkgFFxiRC", "#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger} Bkg RC", nptbins,ptmin,ptmax,100,0.,10.);
282 fhBkgFFxi[1] = new TH2F("BkgFFxiPCG", "#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger} Bkg PCG", nptbins,ptmin,ptmax,100,0.,10.);
283 fhBkgFFxi[2] = new TH2F("BkgFFxiPCJ", "#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger} Bkg PCJ", nptbins,ptmin,ptmax,100,0.,10.);
284 fhBkgFFxi[3] = new TH2F("BkgFFxiMP", "#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger} Bkg MP", nptbins,ptmin,ptmax,100,0.,10.);
285 fhBkgFFxi[4] = new TH2F("BkgFFxiTest","#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger} Bkg Test",nptbins,ptmin,ptmax,100,0.,10.);
286 for(Int_t i=0;i<5;i++){
287 fhBkgFFxi[i]->SetYTitle("#xi");
288 fhBkgFFxi[i]->SetXTitle("p_{T trigger}");
289 outputContainer->Add(fhBkgFFxi[i]) ;
292 fhBkgFFpt[0] = new TH2F("BkgFFptRC", "p_{T i charged} vs p_{T trigger} Bkg RC", nptbins,ptmin,ptmax,100,0.,50.);
293 fhBkgFFpt[1] = new TH2F("BkgFFptPCG", "p_{T i charged} vs p_{T trigger} Bkg PCG", nptbins,ptmin,ptmax,100,0.,50.);
294 fhBkgFFpt[2] = new TH2F("BkgFFptPCJ", "p_{T i charged} vs p_{T trigger} Bkg PCJ", nptbins,ptmin,ptmax,100,0.,50.);
295 fhBkgFFpt[3] = new TH2F("BkgFFptMP", "p_{T i charged} vs p_{T trigger} Bkg MP", nptbins,ptmin,ptmax,100,0.,50.);
296 fhBkgFFpt[4] = new TH2F("BkgFFptTest","p_{T i charged} vs p_{T trigger} Bkg Test", nptbins,ptmin,ptmax,100,0.,50.);
297 for(Int_t i=0;i<5;i++){
298 fhBkgFFpt[i]->SetYTitle("p_{T charged hadron}");
299 fhBkgFFpt[i]->SetXTitle("p_{T trigger}");
300 outputContainer->Add(fhBkgFFpt[i]) ;
303 fhBkgNTracksInCone[0] = new TH2F("BkgNTracksInConeRC", "Number of tracks in cone vs p_{T trigger} Bkg RC", nptbins,ptmin,ptmax,100,0.,150.);
304 fhBkgNTracksInCone[1] = new TH2F("BkgNTracksInConePCG", "Number of tracks in cone vs p_{T trigger} Bkg PCG", nptbins,ptmin,ptmax,100,0.,150.);
305 fhBkgNTracksInCone[2] = new TH2F("BkgNTracksInConePCJ", "Number of tracks in cone vs p_{T trigger} Bkg PCJ", nptbins,ptmin,ptmax,100,0.,150.);
306 fhBkgNTracksInCone[3] = new TH2F("BkgNTracksInConeMP", "Number of tracks in cone vs p_{T trigger} Bkg MP", nptbins,ptmin,ptmax,100,0.,150.);
307 fhBkgNTracksInCone[4] = new TH2F("BkgNTracksInConeTest","Number of tracks in cone vs p_{T trigger} Bkg Test", nptbins,ptmin,ptmax,100,0.,150.);
308 for(Int_t i=0;i<5;i++){
309 fhBkgNTracksInCone[i]->SetYTitle("Number of tracks");
310 fhBkgNTracksInCone[i]->SetXTitle("p_{T trigger}");
311 outputContainer->Add(fhBkgNTracksInCone[i]) ;
314 fhBkgSumPtInCone[0] = new TH2F("BkgSumPtInConeRC", "Sum P_{T} in cone vs p_{T trigger} Bkg RC", nptbins,ptmin,ptmax,100,0.,100.);
315 fhBkgSumPtInCone[1] = new TH2F("BkgSumPtInConePCG", "Sum P_{T} in cone vs p_{T trigger} Bkg PCG", nptbins,ptmin,ptmax,100,0.,100.);
316 fhBkgSumPtInCone[2] = new TH2F("BkgSumPtInConePCJ", "Sum P_{T} in cone vs p_{T trigger} Bkg PCJ", nptbins,ptmin,ptmax,100,0.,100.);
317 fhBkgSumPtInCone[3] = new TH2F("BkgSumPtInConeMP", "Sum P_{T} in cone vs p_{T trigger} Bkg MP", nptbins,ptmin,ptmax,100,0.,100.);
318 fhBkgSumPtInCone[4] = new TH2F("BkgSumPtInConeTest","Sum P_{T} in cone vs p_{T trigger} Bkg Test", nptbins,ptmin,ptmax,100,0.,100.);
319 for(Int_t i=0;i<5;i++){
320 fhBkgSumPtInCone[i]->SetYTitle("Sum P_{T}");
321 fhBkgSumPtInCone[i]->SetXTitle("p_{T trigger}");
322 outputContainer->Add(fhBkgSumPtInCone[i]) ;
325 fhBkgSumPtnTracksInCone[0] = new TH2F("BkgSumPtnTracksInConeRC", "Sum p_{T} / Number of tracks in cone vs p_{T trigger} Bkg RC", nptbins,ptmin,ptmax,100,0.,20.);
326 fhBkgSumPtnTracksInCone[1] = new TH2F("BkgSumPtnTracksInConePCG", "Sum p_{T} / Number of tracks in cone vs p_{T trigger} Bkg PCG", nptbins,ptmin,ptmax,100,0.,20.);
327 fhBkgSumPtnTracksInCone[2] = new TH2F("BkgSumPtnTracksInConePCJ", "Sum p_{T} / Number of tracks in cone vs p_{T trigger} Bkg PCJ", nptbins,ptmin,ptmax,100,0.,20.);
328 fhBkgSumPtnTracksInCone[3] = new TH2F("BkgSumPtnTracksInConeMP", "Sum p_{T} / Number of tracks in cone vs p_{T trigger} Bkg MP", nptbins,ptmin,ptmax,100,0.,20.);
329 fhBkgSumPtnTracksInCone[4] = new TH2F("BkgSumPtnTracksInConeTest","Sum p_{T} / Number of tracks in cone vs p_{T trigger} Bkg Test", nptbins,ptmin,ptmax,100,0.,20.);
330 for(Int_t i=0;i<5;i++){
331 fhBkgSumPtnTracksInCone[i]->SetYTitle("Sum p_{T}/Number of tracks");
332 fhBkgSumPtnTracksInCone[i]->SetXTitle("p_{T trigger}");
333 outputContainer->Add(fhBkgSumPtnTracksInCone[i]) ;
337 //temporary histograms
338 fhNjetsNgammas = new TH2F("NjetsNgammas"," Number of jets vs number of gammas in event",20,0.,100.,10,0.,80.);
339 fhNjetsNgammas->SetYTitle("Number of gammas");
340 fhNjetsNgammas->SetXTitle("Number of jets");
341 outputContainer->Add(fhNjetsNgammas) ;
343 fhCuts = new TH1F("Cuts"," Cuts",10,0.,10.);
344 fhCuts->SetYTitle("Counts");
345 fhCuts->SetXTitle("Cut number");
346 outputContainer->Add(fhCuts) ;
348 fhDeltaPhiBefore = new TH2F("DeltaPhiBefore","#phi_{trigger} - #phi_{jet} vs p_{T trigger}",nptbins,ptmin,ptmax,100,0,6.5);
349 fhDeltaPhiBefore->SetYTitle("#Delta #phi");
350 fhDeltaPhiBefore->SetXTitle("p_{T trigger} (GeV/c)");
351 outputContainer->Add(fhDeltaPhiBefore);
353 fhDeltaEtaBefore = new TH2F("DeltaEtaBefore","#eta_{trigger} - #eta_{jet} vs p_{T trigger}",nptbins,ptmin,ptmax,100,-2,2);
354 fhDeltaEtaBefore->SetYTitle("#Delta #eta");
355 fhDeltaEtaBefore->SetXTitle("p_{T trigger} (GeV/c)");
356 outputContainer->Add(fhDeltaEtaBefore);
358 fhDeltaPtBefore = new TH2F("DeltaPtBefore","p_{T trigger} - p_{T jet} vs p_{T trigger}",nptbins,ptmin,ptmax,100,-50,50);
359 fhDeltaPtBefore->SetYTitle("#Delta p_{T}");
360 fhDeltaPtBefore->SetXTitle("p_{T trigger} (GeV/c)");
361 outputContainer->Add(fhDeltaPtBefore);
363 fhPtRatioBefore = new TH2F("PtRatioBefore","p_{T jet} / p_{T trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,200,0,2.);
364 fhPtRatioBefore->SetYTitle("ratio");
365 fhPtRatioBefore->SetXTitle("p_{T trigger} (GeV/c)");
366 outputContainer->Add(fhPtRatioBefore);
368 fhPtBefore = new TH2F("PtBefore","p_{T jet} vs p_{T trigger}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
369 fhPtBefore->SetYTitle("p_{T jet}(GeV/c)");
370 fhPtBefore->SetXTitle("p_{T trigger} (GeV/c)");
371 outputContainer->Add(fhPtBefore);
373 fhDeltaPhi0PiCorrectBefore = new TH2F("DeltaPhi0PiCorrectBefore","#phi_{trigger} - #phi_{jet} (0,#pi) vs p_{T trigger}",nptbins,ptmin,ptmax,100,0,3.5);
374 fhDeltaPhi0PiCorrectBefore->SetYTitle("#Delta #phi");
375 fhDeltaPhi0PiCorrectBefore->SetXTitle("p_{T trigger} (GeV/c)");
376 outputContainer->Add(fhDeltaPhi0PiCorrectBefore);
378 //temporary jet histograms
379 fhJetPtBefore = new TH1F("JetPtBefore","JetPtBefore",150,-50,100);
380 fhJetPtBefore->SetYTitle("Counts");
381 fhJetPtBefore->SetXTitle("p_{T jet}(GeV/c)");
382 outputContainer->Add(fhJetPtBefore) ;
384 fhJetPt = new TH1F("JetPt","JetPt",150,-50,100);
385 fhJetPt->SetYTitle("Counts");
386 fhJetPt->SetXTitle("p_{T jet}(GeV/c)");
387 outputContainer->Add(fhJetPt) ;
389 fhJetPtMostEne = new TH1F("JetPtMostEne","JetPtMostEne",150,0,150);
390 fhJetPtMostEne->SetYTitle("Counts");
391 fhJetPtMostEne->SetXTitle("p_{T jet}(GeV/c)");
392 outputContainer->Add(fhJetPtMostEne) ;
394 fhJetPhi = new TH1F("JetPhi","JetPhi",130,0,6.5);
395 fhJetPhi->SetYTitle("Counts");
396 fhJetPhi->SetXTitle("#phi_{jet}");
397 outputContainer->Add(fhJetPhi) ;
399 fhJetEta = new TH1F("JetEta","JetEta",100,-1,1);
400 fhJetEta->SetYTitle("Counts");
401 fhJetEta->SetXTitle("#eta_{jet}");
402 outputContainer->Add(fhJetEta) ;
404 fhJetEtaVsPt = new TH2F("JetEtaVsPt","JetEtaVsPt",100,0,100,50,-1,1);
405 fhJetEtaVsPt->SetYTitle("#eta_{jet}");
406 fhJetEtaVsPt->SetXTitle("p_{T,jet}(GeV/c)");
407 outputContainer->Add(fhJetEtaVsPt) ;
409 fhJetPhiVsEta = new TH2F("JetPhiVsEta","JetPhiVsEta",65,0,6.5,50,-1,1);
410 fhJetPhiVsEta->SetYTitle("#eta_{jet}");
411 fhJetPhiVsEta->SetXTitle("#phi_{jet}");
412 outputContainer->Add(fhJetPhiVsEta) ;
414 fhJetEtaVsNpartInJet= new TH2F("JetEtaVsNpartInJet","JetEtaVsNpartInJet",50,-1,1,100,0.,200.);
415 fhJetEtaVsNpartInJet->SetYTitle("N_{tracks-in-jet}");
416 fhJetEtaVsNpartInJet->SetXTitle("#eta_{jet}");
417 outputContainer->Add(fhJetEtaVsNpartInJet) ;
419 fhJetEtaVsNpartInJetBkg= new TH2F("JetEtaVsNpartInJetBkg","JetEtaVsNpartInJetBkg",50,-1,1,100,0.,200.);
420 fhJetEtaVsNpartInJetBkg->SetYTitle("N_{tracks-in-jet}");
421 fhJetEtaVsNpartInJetBkg->SetXTitle("#eta_{jet}");
422 outputContainer->Add(fhJetEtaVsNpartInJetBkg) ;
424 fhJetChBkgEnergyVsPt = new TH2F("JetBkgChEnergyVsPt","JetBkgChEnergyVsPt",100,0,100,90,0,90);
425 fhJetChBkgEnergyVsPt->SetYTitle("Jet Bkg Energy (GeV)");
426 fhJetChBkgEnergyVsPt->SetXTitle("p_{T,jet} (GeV/c)");
427 outputContainer->Add(fhJetChBkgEnergyVsPt);
429 fhJetChAreaVsPt = new TH2F("JetChAreaVsPt","JetChAreaVsPt",100,0,100,50,0,1);
430 fhJetChAreaVsPt->SetYTitle("Jet Area");
431 fhJetChAreaVsPt->SetXTitle("p_{T,jet} (GeV/c)");
432 outputContainer->Add(fhJetChAreaVsPt);
434 if(IsHistogramTracks()){
435 fhTrackPhiVsEta = new TH2F("TrackPhiVsEta","TrackPhiVsEta",65,0,6.5,50,-1,1);
436 fhTrackPhiVsEta->SetYTitle("#eta_{track}");
437 fhTrackPhiVsEta->SetXTitle("#phi_{track}");
438 outputContainer->Add(fhTrackPhiVsEta) ;
440 fhTrackAveTrackPt = new TH1F("TrackAveTrackPt","TrackAveTrackPt",45,0,1.5);
441 fhTrackAveTrackPt->SetYTitle("Counts");
442 fhTrackAveTrackPt->SetXTitle("Average p_{T,track} (GeV/c)");
443 outputContainer->Add(fhTrackAveTrackPt);
445 }//end of IsHistogramTracks()
447 for(Int_t i=0;i<10;i++){
448 fhJetNjetOverPtCut[i] = new TH1F(Form("JetNjetOverPtCut%d", i),Form("JetNjetOverPtCut%d", i),100,0,100);
449 fhJetNjetOverPtCut[i]->SetYTitle("Counts");
450 fhJetNjetOverPtCut[i]->SetXTitle("N_{jets} over threshold");
451 outputContainer->Add(fhJetNjetOverPtCut[i]);
454 fhJetChBkgEnergyVsArea = new TH2F("JetBkgChEnergyVsArea","JetBkgChEnergyVsArea",100,0,100,70,0,0.7);
455 fhJetChBkgEnergyVsArea->SetXTitle("Jet Bkg Energy (GeV)");
456 fhJetChBkgEnergyVsArea->SetYTitle("Area");
457 outputContainer->Add(fhJetChBkgEnergyVsArea);
459 fhJetRhoVsPt = new TH2F("JetRhoVsPt","JetRhoVsPt",100,0,100,100,0,150);
460 fhJetRhoVsPt->SetYTitle("Rho");
461 fhJetRhoVsPt->SetXTitle("p_{T,jet} (GeV/c)");
462 outputContainer->Add(fhJetRhoVsPt);
464 if(IsHistogramJetBkg()){
465 fhJetRhoVsCentrality = new TH2F("JetRhoVsCentrality","JetRhoVsCentrality",100,0,100,100,0,200);
466 fhJetRhoVsCentrality->SetYTitle("Rho");
467 fhJetRhoVsCentrality->SetXTitle("Centrality");
468 outputContainer->Add(fhJetRhoVsCentrality);
471 fhJetNparticlesInJet = new TH1F("JetNparticlesInJet","JetNparticlesInJet",100,0,200);
472 fhJetNparticlesInJet->SetXTitle("N^{particles}");
473 fhJetNparticlesInJet->SetYTitle("N^{jets}");
474 outputContainer->Add(fhJetNparticlesInJet);
476 fhJetDeltaEtaDeltaPhi = new TH2F("JetDeltaEtaDeltaPhi","#Delta #eta^{jet-track} vs. #Delta #phi^{jet-track} for jet tracks",100,-0.8,0.8,100,-0.8,0.8);
477 fhJetDeltaEtaDeltaPhi->SetXTitle("#Delta #eta^{jet-track}");
478 fhJetDeltaEtaDeltaPhi->SetYTitle("#Delta #phi^{jet-track}");
479 outputContainer->Add(fhJetDeltaEtaDeltaPhi );
482 fhJetDeltaEtaDeltaPhiAllTracks = new TH2F("JetDeltaEtaDeltaPhiAllTracks","#Delta #eta^{jet-track} vs. #Delta #phi^{jet-track} for all tracks",100,-3.2,3.2,100,-3.2,3.2);
483 fhJetDeltaEtaDeltaPhiAllTracks->SetXTitle("#Delta #eta^{jet-track}");
484 fhJetDeltaEtaDeltaPhiAllTracks->SetYTitle("#Delta #phi^{jet-track}");
485 outputContainer->Add(fhJetDeltaEtaDeltaPhiAllTracks);
488 if(IsHistogramJetTracks()){
489 fhJetAveTrackPt = new TH1F("JetAveTrackPt","JetAveTrackPt",45,0.,1.5);
490 fhJetAveTrackPt->SetXTitle("Average p_{T,track} (GeV/c)");
491 fhJetAveTrackPt->SetYTitle("Counts");
492 outputContainer->Add(fhJetAveTrackPt);
494 for(Int_t i=0;i<6;i++){
495 if(i==0) fhJetNtracksInJetAboveThr[i] = new TH2F(Form("JetNtracksInJetAboveThr%d", i),Form("JetNtracksInJetAboveThr%d", i),100,0,100,100,0,200);
496 else fhJetNtracksInJetAboveThr[i] = new TH2F(Form("JetNtracksInJetAboveThr%d", i),Form("JetNtracksInJetAboveThr%d", i),100,0,100,100,0,100);
497 fhJetNtracksInJetAboveThr[i]->SetXTitle("p_{T,jet} (GeV/c)");
498 fhJetNtracksInJetAboveThr[i]->SetYTitle("N_{tracks} over threshold");
499 outputContainer->Add(fhJetNtracksInJetAboveThr[i]);
502 for(Int_t i=0;i<5;i++){
503 fhJetRatioNTrkAboveToNTrk[i] = new TH2F(Form("JetRatioNTrkAboveToNTrk%d", i),Form("JetRatioNTrkAboveToNTrk%d", i),100,0,100,40,0,1);
504 fhJetRatioNTrkAboveToNTrk[i]->SetXTitle("p_{T,jet} (GeV/c)");
505 fhJetRatioNTrkAboveToNTrk[i]->SetYTitle("Ratio N_{tracks} over threshold to N_{tracks}");
506 outputContainer->Add(fhJetRatioNTrkAboveToNTrk[i]);
508 fhJetNtrackRatioMostEne[i] = new TH2F(Form("JetNtrackRatioMostEne%d", i),Form("JetNtrackRatioMostEne%d", i),100,0,100,40,0,1);
509 fhJetNtrackRatioMostEne[i]->SetXTitle("p_{T,jet} (GeV/c)");
510 fhJetNtrackRatioMostEne[i]->SetYTitle("Ratio N_{tracks} over threshold to N_{tracks}");
511 outputContainer->Add(fhJetNtrackRatioMostEne[i]);
513 fhJetNtrackRatioJet5GeV[i] = new TH2F(Form("JetNtrackRatioJet5GeV%d", i),Form("JetNtrackRatioJet5GeV%d", i),100,0,100,40,0,1);
514 fhJetNtrackRatioJet5GeV[i]->SetXTitle("p_{T,jet} (GeV/c)");
515 fhJetNtrackRatioJet5GeV[i]->SetYTitle("Ratio N_{tracks} over threshold to N_{tracks}");
516 outputContainer->Add(fhJetNtrackRatioJet5GeV[i]);
518 fhJetNtrackRatioLead5GeV[i] = new TH2F(Form("JetNtrackRatioLead5GeV%d", i),Form("JetNtrackRatioLead5GeV%d", i),100,0,100,40,0,1);
519 fhJetNtrackRatioLead5GeV[i]->SetXTitle("p_{T,jet} (GeV/c)");
520 fhJetNtrackRatioLead5GeV[i]->SetYTitle("Ratio N_{tracks} over threshold to N_{tracks}");
521 outputContainer->Add(fhJetNtrackRatioLead5GeV[i]);
523 }//end of if IsHistogramJetTracks
525 //temporary background jets histograms
526 if(IsHistogramJetBkg()){
527 for(Int_t i=0;i<4;i++){
528 fhBkgJetBackground[i] = new TH1F(Form("BkgJetBackground%d", i),Form("BkgJetBackground%d", i),100,0,200);
529 fhBkgJetBackground[i]->SetXTitle("<#rho> (GeV/c)");
530 fhBkgJetBackground[i]->SetYTitle("Counts");
531 outputContainer->Add(fhBkgJetBackground[i]);
533 fhBkgJetSigma[i] = new TH1F(Form("BkgJetSigma%d", i),Form("BkgJetSigma%d", i),100,0,50);
534 fhBkgJetSigma[i]->SetXTitle("#sigma (GeV/c)");
535 fhBkgJetSigma[i]->SetYTitle("Counts");
536 outputContainer->Add(fhBkgJetSigma[i]);
538 fhBkgJetArea[i] = new TH1F(Form("BkgJetArea%d", i),Form("BkgJetArea%d", i),100,0,1);
539 fhBkgJetArea[i]->SetXTitle("<A>");
540 fhBkgJetArea[i]->SetYTitle("Counts");
541 outputContainer->Add(fhBkgJetArea[i]);
545 //temporary photon histograms
546 fhPhotonPtMostEne = new TH1F("PhotonPtMostEne","PhotonPtMostEne",100,0,100);
547 fhPhotonPtMostEne->SetYTitle("Counts");
548 fhPhotonPtMostEne->SetXTitle("p_{T,#gamma} (GeV/c)");
549 outputContainer->Add(fhPhotonPtMostEne);
551 // fhPhotonIndexMostEne = new TH1F("PhotonIndexMostEne","PhotonIndexMostEne",100,0,100);
552 // fhPhotonIndexMostEne->SetYTitle("Counts");
553 // fhPhotonIndexMostEne->SetXTitle("Index");
554 // outputContainer->Add(fhPhotonIndexMostEne);
556 fhPhotonAverageEnergy = new TH1F("PhotonAverageEnergy","PhotonAverageEnergy",100,0,10);
557 fhPhotonAverageEnergy->SetYTitle("Counts");
558 fhPhotonAverageEnergy->SetXTitle("p_{T,#gamma} (GeV/c)");
559 outputContainer->Add(fhPhotonAverageEnergy);
561 fhPhotonRatioAveEneToMostEne = new TH1F("PhotonRatioAveEneToMostEne","PhotonRatioAveEneToMostEne",100,0,1);
562 fhPhotonRatioAveEneToMostEne->SetYTitle("Counts");
563 fhPhotonRatioAveEneToMostEne->SetXTitle("Ratio");
564 outputContainer->Add(fhPhotonRatioAveEneToMostEne);
566 fhPhotonAverageEnergyMinus1 = new TH1F("PhotonAverageEnergyMinus1","PhotonAverageEnergyMinus1",100,0,10);
567 fhPhotonAverageEnergyMinus1->SetYTitle("Counts");
568 fhPhotonAverageEnergyMinus1->SetXTitle("p_{T,#gamma} (GeV/c)");
569 outputContainer->Add(fhPhotonAverageEnergyMinus1);
571 fhPhotonRatioAveEneMinus1ToMostEne = new TH1F("PhotonRatioAveEneMinus1ToMostEne","PhotonRatioAveEneMinus1ToMostEne",100,0,1);
572 fhPhotonRatioAveEneMinus1ToMostEne->SetYTitle("Counts");
573 fhPhotonRatioAveEneMinus1ToMostEne->SetXTitle("Ratio");
574 outputContainer->Add(fhPhotonRatioAveEneMinus1ToMostEne);
576 fhPhotonNgammaMoreAverageToNgamma = new TH1F("PhotonNgammaMoreAverageToNgamma","PhotonNgammaMoreAverageToNgamma",100,0,1);
577 fhPhotonNgammaMoreAverageToNgamma->SetYTitle("Counts");
578 fhPhotonNgammaMoreAverageToNgamma->SetXTitle("Ratio");
579 outputContainer->Add(fhPhotonNgammaMoreAverageToNgamma);
581 fhPhotonNgammaMoreAverageMinus1ToNgamma = new TH1F("PhotonNgammaMoreAverageMinus1ToNgamma","PhotonNgammaMoreAverageMinus1ToNgamma",100,0,1);
582 fhPhotonNgammaMoreAverageMinus1ToNgamma->SetYTitle("Counts");
583 fhPhotonNgammaMoreAverageMinus1ToNgamma->SetXTitle("Ratio");
584 outputContainer->Add(fhPhotonNgammaMoreAverageMinus1ToNgamma);
586 for(Int_t i=0;i<10;i++){
587 fhPhotonNgammaOverPtCut[i] = new TH1F(Form("PhotonNgammaOverPtCut%d",i),Form("PhotonNgammaOverPtCut%d",i),100,0,100);
588 fhPhotonNgammaOverPtCut[i]->SetYTitle("Counts");
589 fhPhotonNgammaOverPtCut[i]->SetXTitle("N_{#gamma} over threshold");
590 outputContainer->Add(fhPhotonNgammaOverPtCut[i]);
593 fhPhotonBkgRhoVsNtracks = new TH2F("PhotonBkgRhoVsNtracks","PhotonBkgRhoVsNtracks",200,0,2500,75,0,1.5);
594 //fhPhotonBkgRhoVsNtracks->SetXTitle("Counts");
595 fhPhotonBkgRhoVsNtracks->SetXTitle("Ntracks");
596 fhPhotonBkgRhoVsNtracks->SetYTitle("Rho");
597 outputContainer->Add(fhPhotonBkgRhoVsNtracks);
599 fhPhotonBkgRhoVsNclusters = new TH2F("PhotonBkgRhoVsNclusters","PhotonBkgRhoVsNclusters",50,0,100,75,0,1.5);
600 fhPhotonBkgRhoVsNclusters->SetXTitle("Nclusters");
601 fhPhotonBkgRhoVsNclusters->SetYTitle("Rho");
602 outputContainer->Add(fhPhotonBkgRhoVsNclusters);
604 fhPhotonBkgRhoVsCentrality = new TH2F("PhotonBkgRhoVsCentrality","PhotonBkgRhoVsCentrality",100,0,100,75,0,1.5);
605 fhPhotonBkgRhoVsCentrality->SetXTitle("Centrality");
606 fhPhotonBkgRhoVsCentrality->SetYTitle("Rho");
607 outputContainer->Add(fhPhotonBkgRhoVsCentrality);
609 fhPhotonBkgRhoVsNcells = new TH2F("PhotonBkgRhoVsNcells","PhotonBkgRhoVsNcells",100,0,200,75,0,1.5);
610 fhPhotonBkgRhoVsNcells->SetXTitle("N_{cells}");
611 fhPhotonBkgRhoVsNcells->SetYTitle("Rho");
612 outputContainer->Add(fhPhotonBkgRhoVsNcells);
614 fhPhotonPt = new TH1F("PhotonPt","PhotonPt",220,-10,100);
615 fhPhotonPt->SetXTitle("p_{T,#gamma} (GeV/c)");
616 fhPhotonPt->SetYTitle("Counts");
617 outputContainer->Add(fhPhotonPt);
619 fhPhotonPtCorrected = new TH1F("PhotonPtCorrected","PhotonPtCorrected",220,-10,100);
620 fhPhotonPtCorrected->SetXTitle("p_{T,#gamma} (GeV/c)");
621 fhPhotonPtCorrected->SetYTitle("Counts");
622 outputContainer->Add(fhPhotonPtCorrected);
624 fhPhotonPtDiff = new TH1F("PhotonPtDiff","PhotonPtDiff",50,0,10);
625 fhPhotonPtDiff->SetXTitle("p_{T,#gamma} (GeV/c)");
626 fhPhotonPtDiff->SetYTitle("Counts");
627 outputContainer->Add(fhPhotonPtDiff);
629 fhPhotonPtDiffVsNtracks = new TH2F("PhotonPtDiffVsNtracks","PhotonPtDiffVsNtracks",200,0,2500,50,0,5);
630 fhPhotonPtDiffVsNtracks->SetXTitle("N_{tracks}");
631 fhPhotonPtDiffVsNtracks->SetYTitle("<#rho^{#gamma}>*N_{cells}");
632 outputContainer->Add(fhPhotonPtDiffVsNtracks);
634 fhPhotonPtDiffVsNclusters = new TH2F("PhotonPtDiffVsNclusters","PhotonPtDiffVsNclusters",50,0,100,50,0,5);
635 fhPhotonPtDiffVsNclusters->SetXTitle("N_{clusters}");
636 fhPhotonPtDiffVsNclusters->SetYTitle("<#rho^{#gamma}>*N_{cells}");
637 outputContainer->Add(fhPhotonPtDiffVsNclusters);
639 fhPhotonPtDiffVsCentrality = new TH2F("PhotonPtDiffVsCentrality","PhotonPtDiffVsCentrality",100,0,100,50,0,5);
640 fhPhotonPtDiffVsCentrality->SetXTitle("Centrality");
641 fhPhotonPtDiffVsCentrality->SetYTitle("<#rho^{#gamma}>*N_{cells}");
642 outputContainer->Add(fhPhotonPtDiffVsCentrality);
644 fhPhotonPtDiffVsNcells = new TH2F("PhotonPtDiffVsNcells","PhotonPtDiffVsNcells",100,0,200,50,0,5);
645 fhPhotonPtDiffVsNcells->SetXTitle("N_{cells}");
646 fhPhotonPtDiffVsNcells->SetYTitle("<#rho^{#gamma}>*N_{cells}");
647 outputContainer->Add(fhPhotonPtDiffVsNcells);
650 fhPhotonPtCorrectedZoom = new TH1F("PhotonPtCorrectedZoom","PhotonPtCorrectedZoom",200,-5,5);
651 fhPhotonPtCorrectedZoom->SetXTitle("p_{T,#gamma} (GeV/c)");
652 fhPhotonPtCorrectedZoom->SetYTitle("Counts");
653 outputContainer->Add(fhPhotonPtCorrectedZoom);
655 fhPhotonSumPtInCone = new TH1F("PhotonSumPtInCone","PhotonSumPtInCone",100,0,100);
656 fhPhotonSumPtInCone->SetXTitle("#Sigma p_{T,#gamma} (GeV/c)");
657 fhPhotonSumPtInCone->SetYTitle("Counts");
658 outputContainer->Add(fhPhotonSumPtInCone);
660 fhPhotonSumPtCorrectInCone = new TH1F("PhotonSumPtCorrectInCone","PhotonSumPtCorrectInCone",100,-20,80);
661 fhPhotonSumPtCorrectInCone->SetXTitle("#Sigma p_{T,#gamma} (GeV/c)");
662 fhPhotonSumPtCorrectInCone->SetYTitle("Counts");
663 outputContainer->Add(fhPhotonSumPtCorrectInCone);
665 fhPhotonSumPtChargedInCone = new TH1F("PhotonSumPtChargedInCone","PhotonSumPtChargedInCone",100,0,100);
666 fhPhotonSumPtChargedInCone->SetXTitle("#Sigma p_{T,#gamma}^{ch} (GeV/c)");
667 fhPhotonSumPtChargedInCone->SetYTitle("Counts");
668 outputContainer->Add(fhPhotonSumPtChargedInCone);
671 //temporary jet histograms after selection
672 fhSelectedJetPhiVsEta = new TH2F("SelectedJetSelectedPhiVsEta","SelectedJetPhiVsEta",65,0,6.5,50,-1,1);
673 fhSelectedJetPhiVsEta->SetYTitle("#eta_{jet}");
674 fhSelectedJetPhiVsEta->SetXTitle("#phi_{jet}");
675 outputContainer->Add(fhSelectedJetPhiVsEta) ;
677 fhSelectedJetChBkgEnergyVsPtJet = new TH2F("SelectedJetBkgChEnergyVsPtJet","SelectedJetBkgChEnergyVsPtJet",100,0,100,90,0,90);
678 fhSelectedJetChBkgEnergyVsPtJet->SetYTitle("Jet Bkg Energy (GeV)");
679 fhSelectedJetChBkgEnergyVsPtJet->SetXTitle("p_{T,jet} (GeV/c)");
680 outputContainer->Add(fhSelectedJetChBkgEnergyVsPtJet);
682 fhSelectedJetChAreaVsPtJet = new TH2F("SelectedJetChAreaVsPtJet","SelectedJetChAreaVsPtJet",100,0,100,50,0,1);
683 fhSelectedJetChAreaVsPtJet->SetYTitle("Jet Area");
684 fhSelectedJetChAreaVsPtJet->SetXTitle("p_{T,jet} (GeV/c)");
685 outputContainer->Add(fhSelectedJetChAreaVsPtJet);
687 fhSelectedJetNjet = new TH1F("SelectedJetNjet","SelectedJetNjet",100,0,100);
688 fhSelectedJetNjet->SetYTitle("Counts");
689 fhSelectedJetNjet->SetXTitle("N_{jets} per event");
690 outputContainer->Add(fhSelectedJetNjet);
692 fhSelectedNtracks = new TH1F("SelectedNtracks","SelectedNtracks",100,0,2000);
693 fhSelectedNtracks->SetYTitle("Counts");
694 fhSelectedNtracks->SetXTitle("N_{tracks} per event");
695 outputContainer->Add(fhSelectedNtracks);
697 fhSelectedTrackPhiVsEta = new TH2F("SelectedTrackPhiVsEta","SelectedTrackPhiVsEta",65,0,6.5,50,-1,1);
698 fhSelectedTrackPhiVsEta->SetYTitle("#eta_{track}");
699 fhSelectedTrackPhiVsEta->SetXTitle("#phi_{track}");
700 outputContainer->Add(fhSelectedTrackPhiVsEta) ;
702 fhCuts2 = new TH1F("Cuts2","Cuts2",10,0.,10.);
703 fhCuts2->SetYTitle("Counts");
704 fhCuts2->SetXTitle("Cut number");
705 outputContainer->Add(fhCuts2);
707 fhSelectedPhotonNLMVsPt = new TH2F("SelectedPhotonNLMVsPt","SelectedPhotonNLMVsPt",100,0,100,10,0,10);
708 fhSelectedPhotonNLMVsPt->SetYTitle("NLM");
709 fhSelectedPhotonNLMVsPt->SetXTitle("p_{T,#gamma} (GeV/c)");
710 outputContainer->Add(fhSelectedPhotonNLMVsPt);
712 fhSelectedPhotonLambda0VsPt = new TH2F("SelectedPhotonLambda0VsPt","SelectedPhotonLambda0VsPt",100,0,100,50,0,5);
713 fhSelectedPhotonLambda0VsPt->SetYTitle("#lambda_{0}");
714 fhSelectedPhotonLambda0VsPt->SetXTitle("p_{T,#gamma} (GeV/c)");
715 outputContainer->Add(fhSelectedPhotonLambda0VsPt);
718 fhRandomPhiEta[0] = new TH2F("RandomPhiEtaRC", "RandomPhiEtaRC", 100,0,6.5,100,-1.,1.);
719 fhRandomPhiEta[1] = new TH2F("RandomPhiEtaPCG", "RandomPhiEtaPerpConePhoton",100,0,6.5,100,-1.,1.);
720 fhRandomPhiEta[2] = new TH2F("RandomPhiEtaPCJ", "RandomPhiEtaPerpConeJet", 100,0,6.5,100,-1.,1.);
721 fhRandomPhiEta[3] = new TH2F("RandomPhiEtaMP", "RandomPhiEtaMidPoint", 100,0,6.5,100,-1.,1.);
722 fhRandomPhiEta[4] = new TH2F("RandomPhiEtaTest","RandomPhiEtaTest", 100,0,6.5,100,-1.,1.);
723 for(Int_t i=0;i<5;i++){
724 fhRandomPhiEta[i]->SetXTitle("#phi");
725 fhRandomPhiEta[i]->SetYTitle("#eta");
726 outputContainer->Add(fhRandomPhiEta[i]);
729 //MC generated photons and jets information
731 fhMCPhotonCuts = new TH1F("MCPhotonCuts","MCPhotonCuts",10,0.,10.);
732 fhMCPhotonCuts->SetYTitle("Counts");
733 fhMCPhotonCuts->SetXTitle("Cut number");
734 outputContainer->Add(fhMCPhotonCuts);
736 fhMCPhotonPt = new TH1F("MCPhotonPt","MCPhotonPt",100,0.,100.);
737 fhMCPhotonPt->SetYTitle("Counts");
738 fhMCPhotonPt->SetXTitle("p_{T,#gamma}^{gen} (GeV/c)");
739 outputContainer->Add(fhMCPhotonPt);
741 fhMCPhotonEtaPhi = new TH2F("MCPhotonEtaPhi","MCPhotonEtaPhi",65,0,6.5,50,-1,1);
742 fhMCPhotonEtaPhi->SetYTitle("#eta_{#gamma}");
743 fhMCPhotonEtaPhi->SetXTitle("#phi_{#gamma}");
744 outputContainer->Add(fhMCPhotonEtaPhi) ;
746 fhMCJetOrigin = new TH1F("MCJetOrigin","MCJetOrigin",35,-10.,25.);
747 fhMCJetOrigin->SetYTitle("Counts");
748 fhMCJetOrigin->SetXTitle("PDG number");
749 outputContainer->Add(fhMCJetOrigin);
751 fhMCJetNPartVsPt = new TH2F("MCJetNPartVsPt","MCJetNPartVsPt",100,0.,100.,100,0.,100.);
752 fhMCJetNPartVsPt->SetYTitle("N_{particles}");
753 fhMCJetNPartVsPt->SetXTitle("p_{T,full-jet}^{gen} (GeV/c)");
754 outputContainer->Add(fhMCJetNPartVsPt) ;
756 fhMCJetChNPartVsPt = new TH2F("MCJetChNPartVsPt","MCJetChNPartVsPt",100,0.,100.,100,0.,100.);
757 fhMCJetChNPartVsPt->SetYTitle("N_{particles}");
758 fhMCJetChNPartVsPt->SetXTitle("p_{T,charged-jet}^{gen} (GeV/c)");
759 outputContainer->Add(fhMCJetChNPartVsPt) ;
761 fhMCJetNPart150VsPt = new TH2F("MCJetNPart150VsPt","MCJetNPart150VsPt",100,0.,100.,100,0.,100.);
762 fhMCJetNPart150VsPt->SetYTitle("N_{particles (p_{T}>150 MeV/c)}");
763 fhMCJetNPart150VsPt->SetXTitle("p_{T,full-jet}^{gen} (GeV/c)");
764 outputContainer->Add(fhMCJetNPart150VsPt) ;
766 fhMCJetChNPart150VsPt = new TH2F("MCJetChNPart150VsPt","MCJetChNPart150VsPt",100,0.,100.,100,0.,100.);
767 fhMCJetChNPart150VsPt->SetYTitle("N_{particles (p_{T}>150 MeV/c)}");
768 fhMCJetChNPart150VsPt->SetXTitle("p_{T,charged-jet}^{gen} (GeV/c)");
769 outputContainer->Add(fhMCJetChNPart150VsPt) ;
771 fhMCJetChNPart150ConeVsPt = new TH2F("MCJetChNPart150ConeVsPt","MCJetChNPart150ConeVsPt",100,0.,100.,50,0.,50.);
772 fhMCJetChNPart150ConeVsPt->SetYTitle("N_{particles (p_{T}>150 MeV/c)}");
773 fhMCJetChNPart150ConeVsPt->SetXTitle("p_{T,charged-jet,R=0.4}^{gen} (GeV/c)");
774 outputContainer->Add(fhMCJetChNPart150ConeVsPt) ;
776 fhMCJetRatioChFull = new TH1F("MCJetRatioChFull","MCJetRatioChFull",100,0.,1.);
777 fhMCJetRatioChFull->SetYTitle("Counts");
778 fhMCJetRatioChFull->SetXTitle("p_{T,charged-jet}^{gen}/p_{T,full-jet}^{gen}");
779 outputContainer->Add(fhMCJetRatioChFull);
781 fhMCJetRatioCh150Ch = new TH1F("MCJetRatioCh150Ch","MCJetRatioCh150Ch",100,0.,1.);
782 fhMCJetRatioCh150Ch->SetYTitle("Counts");
783 fhMCJetRatioCh150Ch->SetXTitle("p_{T,charged-jet,pT>150MeV/c}^{gen}/p_{T,charged-jet}^{gen}");
784 outputContainer->Add(fhMCJetRatioCh150Ch);
786 fhMCJetEtaPhi = new TH2F("MCJetEtaPhi","MCJetEtaPhi",65,0,6.5,50,-1,1);
787 fhMCJetEtaPhi->SetYTitle("#eta_{jet}");
788 fhMCJetEtaPhi->SetXTitle("#phi_{jet}");
789 outputContainer->Add(fhMCJetEtaPhi) ;
791 fhMCJetChEtaPhi = new TH2F("MCJetChEtaPhi","MCJetChEtaPhi",65,0,6.5,50,-1,1);
792 fhMCJetChEtaPhi->SetYTitle("#eta_{jet}");
793 fhMCJetChEtaPhi->SetXTitle("#phi_{jet}");
794 outputContainer->Add(fhMCJetChEtaPhi) ;
796 fhMCJet150EtaPhi = new TH2F("MCJet150EtaPhi","MCJet150EtaPhi",65,0,6.5,50,-1,1);
797 fhMCJet150EtaPhi->SetYTitle("#eta_{jet}");
798 fhMCJet150EtaPhi->SetXTitle("#phi_{jet}");
799 outputContainer->Add(fhMCJet150EtaPhi) ;
801 fhMCJetCh150EtaPhi = new TH2F("MCJetCh150EtaPhi","MCJetCh150EtaPhi",65,0,6.5,50,-1,1);
802 fhMCJetCh150EtaPhi->SetYTitle("#eta_{jet}");
803 fhMCJetCh150EtaPhi->SetXTitle("#phi_{jet}");
804 outputContainer->Add(fhMCJetCh150EtaPhi) ;
806 fhMCJetCh150ConeEtaPhi = new TH2F("MCJetCh150ConeEtaPhi","MCJetCh150ConeEtaPhi",65,0,6.5,50,-1,1);
807 fhMCJetCh150ConeEtaPhi->SetYTitle("#eta_{jet}");
808 fhMCJetCh150ConeEtaPhi->SetXTitle("#phi_{jet}");
809 outputContainer->Add(fhMCJetCh150ConeEtaPhi) ;
814 fTreeGJ=new TTree("fTreeGJ","fTreeGJ");
815 fTreeGJ->Branch("fGamPt" ,&fGamPt ,"fGamPt/D");//! fGamPt
816 fTreeGJ->Branch("fGamLambda0" ,&fGamLambda0 ,"fGamLambda0/D");
817 fTreeGJ->Branch("fGamNLM" ,&fGamNLM ,"fGamNLM/I");
818 fTreeGJ->Branch("fGamSumPtCh" ,&fGamSumPtCh ,"fGamSumPtCh/D");
819 fTreeGJ->Branch("fGamNtracks" ,&fGamNtracks ,"fGamNtracks/I");
820 fTreeGJ->Branch("fGamTime" ,&fGamTime ,"fGamTime/D");
821 fTreeGJ->Branch("fGamNcells" ,&fGamNcells ,"fGamNcells/I");
822 fTreeGJ->Branch("fGamEta" ,&fGamEta ,"fGamEta/D");
823 fTreeGJ->Branch("fGamPhi" ,&fGamPhi ,"fGamPhi/D");
824 fTreeGJ->Branch("fGamSumPtNeu",&fGamSumPtNeu ,"fGamSumPtNeu/D");
825 fTreeGJ->Branch("fGamNclusters" ,&fGamNclusters ,"fGamNclusters/I");
826 fTreeGJ->Branch("fGamAvEne" ,&fGamAvEne ,"fGamAvEne/D");
827 fTreeGJ->Branch("fJetPhi" ,&fJetPhi ,"fJetPhi/D");
828 fTreeGJ->Branch("fJetEta" ,&fJetEta ,"fJetEta/D");
829 fTreeGJ->Branch("fJetPt" ,&fJetPt ,"fJetPt/D");
830 fTreeGJ->Branch("fJetBkgChEne",&fJetBkgChEne ,"fJetBkgChEne/D");
831 fTreeGJ->Branch("fJetArea" ,&fJetArea ,"fJetArea/D");
832 fTreeGJ->Branch("fJetNtracks" ,&fJetNtracks ,"fJetNtracks/I");
833 fTreeGJ->Branch("fJetNtracks1" ,&fJetNtracks1 ,"fJetNtracks1/I");
834 fTreeGJ->Branch("fJetNtracks2" ,&fJetNtracks2 ,"fJetNtracks2/I");
835 fTreeGJ->Branch("fJetRho" ,&fJetRho ,"fJetRho/D");
836 fTreeGJ->Branch("fEventNumber",&fEventNumber ,"fEventNumber/I");
837 fTreeGJ->Branch("fNtracks" ,&fNtracks ,"fNtracks/I");
838 fTreeGJ->Branch("fZvertex" ,&fZvertex ,"fZvertex/D");
839 fTreeGJ->Branch("fCentrality" ,&fCentrality ,"fCentrality/D");
840 fTreeGJ->Branch("fIso" ,&fIso ,"fIso/O");
841 fTreeGJ->Branch("fGamRho" ,&fGamRho ,"fGamRho/D");
843 fTreeGJ->Branch("fMCGamPt" ,&fMCGamPt ,"fMCGamPt/D" );
844 fTreeGJ->Branch("fMCGamEta" ,&fMCGamEta ,"fMCGamEta/D" );
845 fTreeGJ->Branch("fMCGamPhi" ,&fMCGamPhi ,"fMCGamPhi/D" );
846 fTreeGJ->Branch("fMCPartonType" ,&fMCPartonType ,"fMCPartonType/I" );
847 fTreeGJ->Branch("fMCJetPt" ,&fMCJetPt ,"fMCJetPt/D" );
848 fTreeGJ->Branch("fMCJetChPt" ,&fMCJetChPt ,"fMCJetChPt/D" );
849 fTreeGJ->Branch("fMCJet150Pt" ,&fMCJet150Pt ,"fMCJet150Pt/D" );
850 fTreeGJ->Branch("fMCJetCh150Pt" ,&fMCJetCh150Pt ,"fMCJetCh150Pt/D" );
851 fTreeGJ->Branch("fMCJetNPart" ,&fMCJetNPart ,"fMCJetNPart/I" );
852 fTreeGJ->Branch("fMCJetChNPart" ,&fMCJetChNPart ,"fMCJetChNPart/I" );
853 fTreeGJ->Branch("fMCJet150NPart" ,&fMCJet150NPart ,"fMCJet150NPart/I" );
854 fTreeGJ->Branch("fMCJetCh150NPart" ,&fMCJetCh150NPart ,"fMCJetCh150NPart/I");
855 fTreeGJ->Branch("fMCJetEta" ,&fMCJetEta ,"fMCJetEta/D" );
856 fTreeGJ->Branch("fMCJetPhi" ,&fMCJetPhi ,"fMCJetPhi/D" );
857 fTreeGJ->Branch("fMCJetChEta" ,&fMCJetChEta ,"fMCJetChEta/D" );
858 fTreeGJ->Branch("fMCJetChPhi" ,&fMCJetChPhi ,"fMCJetChPhi/D" );
859 fTreeGJ->Branch("fMCJet150Eta" ,&fMCJet150Eta ,"fMCJet150Eta/D" );
860 fTreeGJ->Branch("fMCJet150Phi" ,&fMCJet150Phi ,"fMCJet150Phi/D" );
861 fTreeGJ->Branch("fMCJetCh150Eta" ,&fMCJetCh150Eta ,"fMCJetCh150Eta/D" );
862 fTreeGJ->Branch("fMCJetCh150Phi" ,&fMCJetCh150Phi ,"fMCJetCh150Phi/D" );
864 fTreeGJ->Branch("fMCJetCh150ConePt" ,&fMCJetCh150ConePt ,"fMCJetCh150ConePt/D" );
865 fTreeGJ->Branch("fMCJetCh150ConeNPart" ,&fMCJetCh150ConeNPart ,"fMCJetCh150ConeNPart/I");
866 fTreeGJ->Branch("fMCJetCh150ConeEta" ,&fMCJetCh150ConeEta ,"fMCJetCh150ConeEta/D" );
867 fTreeGJ->Branch("fMCJetCh150ConePhi" ,&fMCJetCh150ConePhi ,"fMCJetCh150ConePhi/D" );
869 outputContainer->Add(fTreeGJ);
872 return outputContainer;
876 //_______________________________________________________
877 void AliAnaParticleJetFinderCorrelation::InitParameters()
879 //printf("InitParameters\n");
880 //Initialize the parameters of the analysis.
881 SetInputAODName("PWG4Particle");
882 AddToHistogramsName("AnaJetFinderCorr_");
884 fDeltaPhiMinCut = 2.6 ;
885 fDeltaPhiMaxCut = 3.7 ;
889 fPtThresholdInCone = 0.5 ;
890 fUseJetRefTracks = kFALSE ;
891 fMakeCorrelationInHistoMaker = kFALSE ;
892 fSelectIsolated = kTRUE;
894 fJetMinPt = 15. ; //GeV/c
895 fJetAreaFraction = 0.6 ;
896 fJetBranchName = "jets";
897 fBkgJetBranchName = "jets";
898 fGammaConeSize = 0.4;
899 fUseBackgroundSubtractionGamma = kFALSE;
901 fMostEnergetic = kFALSE;
902 fMostOpposite = kTRUE;
903 fUseHistogramJetBkg = kTRUE;
904 fUseHistogramTracks = kTRUE;
905 fUseHistogramJetTracks = kTRUE;
909 //__________________________________________________________________
910 Int_t AliAnaParticleJetFinderCorrelation::SelectJet(AliAODPWG4Particle * particle, TClonesArray *aodRecJets)
912 //Input for jets is TClonesArray *aodRecJets
913 //Returns the index of the jet that is opposite to the particle
914 //printf(" SelectJet ");
916 Double_t particlePt=particle->Pt();
917 if(fUseBackgroundSubtractionGamma) {
918 Int_t clusterIDtmp = particle->GetCaloLabel(0) ;
920 AliVCluster *cluster=0;
921 if(!(clusterIDtmp<0) ){
923 TObjArray* clusters = GetEMCALClusters();
924 cluster = FindCluster(clusters,clusterIDtmp,iclustmp);
925 nCells = cluster->GetNCells();
927 particlePt-=(fGamRho*nCells);
930 //printf("Particle with negative or 0 pt\n");
934 Int_t njets = aodRecJets->GetEntriesFast();
935 AliAODJet * jet = 0 ;
937 Double_t dphiprev= 10000;
938 Double_t particlePhi=particle->Phi();
939 Double_t deltaPhi=-10000.;// in the range (0; 2*pi)
942 for(Int_t ijet = 0; ijet < njets ; ijet++){
943 jet = dynamic_cast<AliAODJet*>(aodRecJets->At(ijet));
946 AliInfo("Jet not in container");
949 fhCuts2->Fill(2.,1.);
951 if(jetPt<fJetMinPt) continue;
952 fhCuts2->Fill(3.,1.);
953 if(fBackgroundJetFromReader ){
954 jetPt-= (fJetRho * jet->EffectiveAreaCharged() );
956 if(jetPt<0.) continue;
957 //put jet eta requirement here |eta_jet|<0.9-jet_cone_size
958 fhCuts2->Fill(4.,1.);
959 if(TMath::Abs(jet->Eta()) > (0.9 - fJetConeSize) ) continue;
960 fhCuts2->Fill(5.,1.);
961 if(jet->EffectiveAreaCharged()<fJetAreaFraction*TMath::Pi()*fJetConeSize*fJetConeSize) continue;
962 fhCuts2->Fill(6.,1.);
963 //printf("jet found\n");
964 Double_t deltaPhi0pi = TMath::Abs(particle->Phi()-jet->Phi());
965 Double_t ratio = jetPt/particlePt;
967 deltaPhi = particlePhi - jet->Phi() ;
968 if ( deltaPhi0pi > TMath::Pi() ) deltaPhi0pi = 2. * TMath::Pi() - deltaPhi0pi ;
969 if(deltaPhi<0) deltaPhi +=(TMath::Pi()*2.);
971 fhDeltaPtBefore ->Fill(particlePt, particlePt - jetPt);
972 fhDeltaPhiBefore->Fill(particlePt, deltaPhi);
973 fhDeltaEtaBefore->Fill(particlePt, particle->Eta() - jet->Eta());
974 fhPtRatioBefore ->Fill(particlePt, jetPt/particlePt);
975 fhPtBefore ->Fill(particlePt, jetPt);
977 fhDeltaPhi0PiCorrectBefore->Fill(particlePt, deltaPhi0pi);//good
980 printf("AliAnaParticleJetFinderCorrelation::SelectJet() - Jet %d, Ratio pT %2.3f, Delta phi %2.3f\n",ijet,ratio,deltaPhi);
982 if((deltaPhi > fDeltaPhiMinCut) && (deltaPhi < fDeltaPhiMaxCut) &&
983 (ratio > fRatioMinCut) && (ratio < fRatioMaxCut) &&
984 (TMath::Abs(deltaPhi-TMath::Pi()) < TMath::Abs(dphiprev-TMath::Pi()) )){//In case there is more than one jet select the most opposite.
994 //__________________________________________________________________
995 void AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD()
997 //Particle-Jet Correlation Analysis, fill AODs
998 // printf("I use MakeAnalysisFillAOD\n");
999 //Get the event, check if there are AODs available, if not, skip this analysis
1000 AliAODEvent * event = NULL;
1002 // cout<<"GetReader()->GetOutputEvent() "<<GetReader()->GetOutputEvent()<<endl;
1003 // cout<<"GetReader()->GetDataType() "<<GetReader()->GetDataType() <<endl;
1004 // cout<<"AliCaloTrackReader::kAOD "<<AliCaloTrackReader::kAOD<<endl;
1005 // cout<<"GetReader()->GetInputEvent() "<<GetReader()->GetInputEvent()<<endl;
1007 if(GetReader()->GetOutputEvent())
1009 event = dynamic_cast<AliAODEvent*>(GetReader()->GetOutputEvent());
1011 else if(GetReader()->GetDataType() == AliCaloTrackReader::kAOD)
1013 event = dynamic_cast<AliAODEvent*>(GetReader()->GetInputEvent());
1017 if(GetDebug() > 3) AliInfo("There are no jets available for this analysis");
1021 if(!GetInputAODBranch() || !event)
1023 AliFatal(Form("No input particles in AOD with name branch < %s > \n",
1024 GetInputAODName().Data()));
1025 return; // Trick coverity
1028 if(strcmp(GetInputAODBranch()->GetClass()->GetName(), "AliAODPWG4ParticleCorrelation"))
1030 AliFatal(Form("Wrong type of AOD object, change AOD class name in input AOD: It should be <AliAODPWG4ParticleCorrelation> and not <%s>",
1031 GetInputAODBranch()->GetClass()->GetName()));
1032 return; // Trick coverity
1036 // non-standard jets
1039 TClonesArray *aodRecJets = 0;
1040 //if(IsNonStandardJetFromReader()){//jet branch from reader
1041 if(GetDebug() > 3) AliInfo(Form("GetNonStandardJets function (from reader) is called"));
1042 aodRecJets = GetNonStandardJets();
1043 if(GetDebug() > 3) AliInfo(Form("aodRecJets %p",aodRecJets));
1046 if(GetDebug() > 3) event->Print();
1047 AliFatal("List of jets is null");
1050 nJets=aodRecJets->GetEntries();
1051 if(GetDebug() > 3) printf("nJets %d\n",nJets);
1056 //printf("Why number of jets = 0? Check what type of collision it is. If PbPb -problem.\n");
1063 AliAODJetEventBackground* aodBkgJets = 0;
1064 if(IsBackgroundJetFromReader()){//jet branch from reader
1065 if(GetDebug() > 3) printf("GetBackgroundJets function is called\n");
1066 aodBkgJets = GetBackgroundJets();
1067 if(GetDebug() > 3) printf("aodBkgJets %p\n",aodBkgJets);
1068 if(aodBkgJets==0x0){
1069 if(GetDebug() > 3) event->Print();
1070 AliFatal("No jet background found\n");
1071 return; // Trick coverity
1073 if(GetDebug() > 3) aodBkgJets->Print("c");
1076 Double_t rhoEvent=0.;
1077 if(aodBkgJets && IsBackgroundJetFromReader() ) {
1078 rhoEvent = aodBkgJets->GetBackground(2);//hardcoded
1083 Int_t ntrig = GetInputAODBranch()->GetEntriesFast() ;
1085 printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - Begin jet finder correlation analysis, fill AODs \n");
1086 printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - In particle branch aod entries %d\n", ntrig);
1087 printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - In standard jet branch aod entries %d\n", event->GetNJets());
1088 printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - In non standard jet branch aod entries %d\n", nJets);
1091 //if(nJets==0) return;//to speed up
1092 // cout<<"ntrig po return "<<ntrig<<endl;
1095 //calculate average cell energy without most energetic photon
1097 Double_t medianPhotonRho=0.;
1098 TObjArray* clusters = GetEMCALClusters();
1100 Int_t iclustmp = -1;
1101 AliVCluster *cluster=0;
1103 if(IsBackgroundSubtractionGamma()){
1105 // Find most energetic photon without background subtraction
1109 AliAODPWG4ParticleCorrelation* particlecorr =0;
1110 for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1111 particlecorr = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1112 if(particlecorr->Pt() > maxPt) {
1113 maxPt = particlecorr->Pt();
1119 // calculate background energy per cell
1121 Int_t numberOfcells=0;
1123 Double_t *photonRhoArr=new Double_t[ntrig-1];
1124 Int_t photonRhoArrayIndex=0;
1125 for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1126 particlecorr = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1127 if(iaod==maxIndex) continue;
1128 clusterIDtmp = particlecorr->GetCaloLabel(0) ;
1129 if(clusterIDtmp < 0) continue;
1130 cluster = FindCluster(clusters,clusterIDtmp,iclustmp);
1131 photonRhoArr[photonRhoArrayIndex]=particlecorr->Pt()/ cluster->GetNCells();
1132 numberOfcells+=cluster->GetNCells();
1133 photonRhoArrayIndex++;
1135 if(photonRhoArrayIndex>0) medianPhotonRho=TMath::Median(photonRhoArrayIndex,photonRhoArr);
1136 delete [] photonRhoArr;
1138 }//end of if background calculation for gamma
1139 fGamRho = medianPhotonRho;
1143 //take most energetic photon and most energetic jet and combine
1147 // most energetic photon with background subtraction
1149 Double_t mostEnePhotonPt=0.;
1150 Int_t indexMostEnePhoton=-1;
1151 AliAODPWG4ParticleCorrelation* particle =0;
1152 Double_t ptCorrect=0.;
1154 for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1155 particle = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1156 clusterIDtmp = particle->GetCaloLabel(0) ;
1157 if(!(clusterIDtmp<0)){
1158 cluster = FindCluster(clusters,clusterIDtmp,iclustmp);
1159 nCells = cluster->GetNCells();
1161 ptCorrect = particle->Pt() - medianPhotonRho * nCells;
1162 if( ptCorrect > mostEnePhotonPt ){
1163 mostEnePhotonPt = ptCorrect;
1164 indexMostEnePhoton = iaod ;
1167 // printf ("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - Photon with index %d selected \n",indexMostEnePhoton);
1169 // most energetic jet with background subtraction
1171 Double_t mostEneJetPt=0.;
1172 Int_t indexMostEneJet=-1;
1173 AliAODJet * jet = 0 ;
1174 //Double_t ptCorrect=0.;
1175 for(Int_t ijet = 0; ijet < nJets ; ijet++){
1176 jet = dynamic_cast<AliAODJet*>(aodRecJets->At(ijet));
1178 if(jet->Pt()<fJetMinPt) continue;
1179 if(TMath::Abs(jet->Eta()) > (0.9 - fJetConeSize) ) continue;
1180 if(jet->EffectiveAreaCharged()<fJetAreaFraction*TMath::Pi()*fJetConeSize*fJetConeSize) continue;
1181 ptCorrect = jet->Pt() - rhoEvent * jet->EffectiveAreaCharged();
1182 if(ptCorrect > mostEneJetPt){
1183 mostEneJetPt = ptCorrect;
1184 indexMostEneJet = ijet;
1187 // printf ("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - Jet with index %d selected \n",indexMostEneJet);
1190 // assign jet to photon
1192 if(indexMostEneJet>=0 && indexMostEnePhoton>=0){
1193 particle = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(indexMostEnePhoton));
1194 jet = dynamic_cast<AliAODJet*>(aodRecJets-> At(indexMostEneJet));
1195 if(jet)particle->SetRefJet(jet);
1197 }//end of take most energetic photon and most ene. jet after bkg subtraction
1200 //Bool_t isJetFound=kFALSE;//new
1201 //Loop on stored AOD particles, trigger
1202 for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1203 AliAODPWG4ParticleCorrelation* particle = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1205 //Correlate with jets
1206 Int_t ijet = SelectJet(particle,aodRecJets);//input for jets is TClonesArray
1209 if(GetDebug() > 2) printf ("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - Jet with index %d selected \n",ijet);
1210 AliAODJet *jet = dynamic_cast<AliAODJet*>(aodRecJets-> At(ijet));
1211 if(jet)particle->SetRefJet(jet);
1212 //printf("Most opposite found\n");
1215 // if(GetReader()->WriteDeltaAODToFile() && isJetFound) WriteJetsToOutputBranch(aodRecJets);
1216 }//end of take most opposite photon and jet after bkg subtraction
1219 if(GetDebug() > 1 ) printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - End fill AODs \n");
1222 //__________________________________________________________________
1223 void AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms()
1225 //Particle-Jet Correlation Analysis, fill histograms
1226 if(GetDebug() > 3 ) {
1227 printf("I use MakeAnalysisFillHistograms\n");
1228 printf("ntracks before iso %d\n",GetCTSTracks()->GetEntriesFast() );
1231 //Get the event, check if there are AODs available, if not, skip this analysis
1232 AliAODEvent * event = NULL;
1234 //printf("GetReader()->GetOutputEvent() %d\n",GetReader()->GetOutputEvent() );
1235 //printf("GetReader()->GetDataType() %d\n",GetReader()->GetDataType() );
1236 //printf("AliCaloTrackReader::kAOD %d\n",AliCaloTrackReader::kAOD );
1237 //printf("GetReader()->GetInputEvent() %d\n",GetReader()->GetInputEvent() );
1239 if(GetReader()->GetOutputEvent())
1241 event = dynamic_cast<AliAODEvent*>(GetReader()->GetOutputEvent());
1243 else if(GetReader()->GetDataType() == AliCaloTrackReader::kAOD)
1245 event = dynamic_cast<AliAODEvent*>(GetReader()->GetInputEvent());
1248 if(GetDebug() > 3) printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - There are no jets available for this analysis\n");
1252 if(!GetInputAODBranch() || !event){
1254 AliFatal(Form("No input particles in AOD with name branch < %s > \n",
1255 GetInputAODName().Data()));
1256 return; // Trick coverity
1260 TClonesArray *aodRecJets = 0;
1261 //if(IsNonStandardJetFromReader()){//branch read in reader from reader
1262 if (GetDebug() > 3) printf("GetNonStandardJets function (from reader) is called\n");
1263 aodRecJets = GetNonStandardJets();
1264 if(aodRecJets==0x0){
1265 if(GetDebug() > 3) event->Print();
1266 AliFatal("Jets container not found\n");
1267 return; // trick coverity
1269 nJets=aodRecJets->GetEntries();
1272 // printf("Why number of jets = 0? Check what type of collision it is. If PbPb -problem.\n");
1273 GetReader()->FillInputNonStandardJets();
1274 aodRecJets = GetNonStandardJets();
1275 if(aodRecJets) nJets=aodRecJets->GetEntries();
1276 // printf("nJets = %d\n",nJets);
1283 AliAODJetEventBackground* aodBkgJets = 0;
1284 if(IsBackgroundJetFromReader()){//jet branch from reader
1285 if(GetDebug() > 3) printf("GetBackgroundJets function is called\n");
1286 aodBkgJets = GetBackgroundJets();
1287 if(GetDebug() > 3) printf("aodBkgJets %p\n",aodBkgJets);
1288 if(aodBkgJets==0x0){
1289 if(GetDebug() > 3) event->Print();
1290 AliFatal("No jet background container found");
1291 return; // trick coverity
1293 if(GetDebug() > 3) aodBkgJets->Print("c");
1298 // only background jets informations
1300 // Float_t pTback = 0;
1301 Double_t rhoEvent=0.;
1303 if(IsBackgroundJetFromReader() ) rhoEvent = aodBkgJets->GetBackground(2);
1304 if(IsHistogramJetBkg()) {
1305 fhJetRhoVsCentrality->Fill(GetEventCentrality(),rhoEvent);
1306 for(Int_t i=0;i<4;i++){
1307 fhBkgJetBackground[i]->Fill(aodBkgJets->GetBackground(i));
1308 fhBkgJetSigma[i]->Fill(aodBkgJets->GetSigma(i));
1309 fhBkgJetArea[i]->Fill(aodBkgJets->GetMeanarea(i));
1311 }//end of if fill HistogramJetBkg
1312 }//end if aodBkgJets exists
1315 //only track information
1317 Int_t nCTSTracks = GetCTSTracks()->GetEntriesFast();
1318 AliAODTrack *aodtrack;
1320 if(IsHistogramTracks()) {
1321 Double_t sumTrackPt=0;
1322 for(itrack = 0; itrack < nCTSTracks ; itrack++){
1323 aodtrack = dynamic_cast <AliAODTrack*>(GetCTSTracks()->At(itrack));
1324 if(!aodtrack) continue;
1325 fhTrackPhiVsEta->Fill(aodtrack->Phi(),aodtrack->Eta());
1326 sumTrackPt+=aodtrack->Pt();
1329 fhTrackAveTrackPt->Fill(sumTrackPt/nCTSTracks);
1330 }//end if IsHistogramTracks
1333 //only jet informations
1335 AliAODJet * jettmp = 0 ;
1336 Double_t jetPttmp = 0.;
1337 Double_t var1tmp = 0.;
1338 Double_t var2tmp = 0.;
1339 // fhJetNjet->Fill(nJets);
1340 Double_t ptMostEne=0;
1341 // Int_t indexMostEne=-1;
1342 Int_t nJetsOverThreshold[10]={nJets,0,0,0,0,0,0,0,0,0};
1344 Double_t sumJetTrackPt=0.;
1345 Int_t sumNJetTrack=0;
1346 Int_t nTracksInJet=0;
1348 for(Int_t ijet = 0; ijet < nJets ; ijet++){
1349 jettmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ijet));
1350 if(!jettmp) continue;
1351 fhJetPtBefore->Fill(jettmp->Pt());
1352 jetPttmp = jettmp->Pt() - rhoEvent * jettmp->EffectiveAreaCharged();//<<---changed here
1354 //calculate number of tracks above 1,2,3,4,5 GeV in jet
1355 AliVTrack* jettrack = 0x0 ;
1357 Int_t nTrackThrGeV[5]={0,0,0,0,0};
1358 nTracksInJet=(jettmp->GetRefTracks())->GetEntriesFast();
1359 fhJetNparticlesInJet->Fill(nTracksInJet);
1360 if(nTracksInJet==0) continue;
1361 sumNJetTrack+=nTracksInJet;
1362 for(itrack=0;itrack<nTracksInJet;itrack++){
1363 jettrack=(AliVTrack *) ((jettmp->GetRefTracks())->At(itrack));
1364 if(!jettrack) continue;
1366 fhJetDeltaEtaDeltaPhi->Fill(jettmp->Eta()-jettrack->Eta(),jettmp->Phi()-jettrack->Phi());
1367 sumJetTrackPt+=jettrack->Pt();
1368 if(IsHistogramJetTracks()){
1369 if(jettrack->Pt()>1.) nTrackThrGeV[0]++;
1370 if(jettrack->Pt()>2.) nTrackThrGeV[1]++;
1371 if(jettrack->Pt()>3.) nTrackThrGeV[2]++;
1372 if(jettrack->Pt()>4.) nTrackThrGeV[3]++;
1373 if(jettrack->Pt()>5.) nTrackThrGeV[4]++;
1374 }//end of if IsHistogramJetTracks
1375 }//end of jet track loop
1376 //printf("jetPt %f ntrks %d ntrks>1,2,3,4,5GeV in jet %d, %d, %d, %d, %d\n",jetPttmp,nTracksInJet,nTrackThrGeV[0],nTrackThrGeV[1],nTrackThrGeV[2],nTrackThrGeV[3],nTrackThrGeV[4]);
1378 for(itrack = 0; itrack < nCTSTracks ; itrack++){
1379 aodtrack = dynamic_cast <AliAODTrack*>(GetCTSTracks()->At(itrack));
1380 if(aodtrack) fhJetDeltaEtaDeltaPhiAllTracks->Fill(jettmp->Eta()-aodtrack->Eta(),jettmp->Phi()-aodtrack->Phi());
1384 if(IsHistogramJetTracks()){
1385 fhJetNtracksInJetAboveThr[0]->Fill(jetPttmp,nTracksInJet);//all jets
1387 for(itrk=0;itrk<5;itrk++) {
1388 fhJetNtracksInJetAboveThr[itrk+1]->Fill(jetPttmp,nTrackThrGeV[itrk]);//all jets
1389 fhJetRatioNTrkAboveToNTrk[itrk]->Fill(jetPttmp,(Double_t)nTrackThrGeV[itrk]/(Double_t)nTracksInJet);//all jets
1391 if(ijet==0){//most ene jet
1392 for(itrk=0;itrk<5;itrk++)
1393 fhJetNtrackRatioMostEne[itrk]->Fill(jetPttmp,(Double_t)nTrackThrGeV[itrk]/(Double_t)nTracksInJet);
1395 if(jetPttmp>5){//jet with pt>5GeV/c
1396 for(itrk=0;itrk<5;itrk++)
1397 fhJetNtrackRatioJet5GeV[itrk]->Fill(jetPttmp,(Double_t)nTrackThrGeV[itrk]/(Double_t)nTracksInJet);
1399 if(nTrackThrGeV[4]>0){//jet with leading particle pt>5GeV/c
1400 for(itrk=0;itrk<5;itrk++)
1401 fhJetNtrackRatioLead5GeV[itrk]->Fill(jetPttmp,(Double_t)nTrackThrGeV[itrk]/(Double_t)nTracksInJet);
1403 }//end of if IsHistogramJetTracks
1406 Double_t effectiveChargedBgEnergy=(IsBackgroundJetFromReader()?rhoEvent * jettmp->EffectiveAreaCharged():jettmp->ChargedBgEnergy());
1409 fhJetChBkgEnergyVsArea->Fill(effectiveChargedBgEnergy,jettmp->EffectiveAreaCharged());
1410 //if(jettmp->EffectiveAreaCharged()>0)
1411 fhJetRhoVsPt->Fill(jetPttmp,jettmp->ChargedBgEnergy()*jettmp->EffectiveAreaCharged());
1413 if(jetPttmp>ptMostEne) {
1414 ptMostEne = jetPttmp;
1415 //indexMostEne=ijet;
1417 fhJetPt->Fill(jetPttmp);
1418 fhJetChBkgEnergyVsPt->Fill(jetPttmp,effectiveChargedBgEnergy);
1419 fhJetChAreaVsPt->Fill(jetPttmp,jettmp->EffectiveAreaCharged());
1420 if(GetDebug()>5) printf("ChargedBgEnergy %f EffectiveAreaCharged %f\n", jettmp->ChargedBgEnergy(),jettmp->EffectiveAreaCharged());
1421 for(iCounter=1;iCounter<10;iCounter++){
1422 if(jetPttmp>iCounter) nJetsOverThreshold[iCounter]++;
1425 var1tmp = jettmp->Phi();
1426 var2tmp = jettmp->Eta();
1427 fhJetPhi->Fill(var1tmp);
1428 fhJetEta->Fill(var2tmp);
1429 fhJetPhiVsEta->Fill(var1tmp,var2tmp);
1430 fhJetEtaVsPt->Fill(jetPttmp,var2tmp);
1431 fhJetEtaVsNpartInJet->Fill(var2tmp,nTracksInJet);
1433 fhJetEtaVsNpartInJetBkg->Fill(var2tmp,nTracksInJet);
1436 if(IsHistogramJetTracks()){
1438 //printf("average track pt %f\n",sumJetTrackPt/sumNJetTrack);
1439 fhJetAveTrackPt->Fill(sumJetTrackPt/sumNJetTrack);
1441 }//end of if IsHistogramJetTracks
1444 fhJetPtMostEne->Fill(ptMostEne);
1445 for(iCounter=0;iCounter<10;iCounter++){
1446 fhJetNjetOverPtCut[iCounter]->Fill(nJetsOverThreshold[iCounter]);
1447 nJetsOverThreshold[iCounter]=0;
1450 //end of only jet part
1455 Int_t ntrig = GetInputAODBranch()->GetEntriesFast() ;
1457 printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - Begin jet finder correlation analysis, fill histograms \n");
1458 printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - In particle branch aod entries %d\n", ntrig);
1459 printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - In jet output branch aod entries %d\n", event->GetNJets());
1461 fhNjetsNgammas->Fill(nJets,ntrig);
1463 //if(nJets==0) return;//to speed up
1464 //printf("ntrig %d, njets %d\n",ntrig,nJets);
1467 //Fill only temporary photon histograms
1473 nJetsOverThreshold[0]=ntrig;
1474 for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1475 AliAODPWG4ParticleCorrelation* particlecorr = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1476 tmpPt = particlecorr->Pt();
1481 for(iCounter=1;iCounter<10;iCounter++){
1482 if(tmpPt>iCounter) nJetsOverThreshold[iCounter]++;
1486 for(iCounter=0;iCounter<10;iCounter++){
1487 fhPhotonNgammaOverPtCut[iCounter]->Fill(nJetsOverThreshold[iCounter]);
1488 // nJetsOverThreshold[iCounter]=0;
1492 TObjArray* clusters = GetEMCALClusters();
1493 //printf("calculate median bkg energy for photons ");
1494 Double_t medianPhotonRho=0.;
1496 Int_t iclustmp = -1;
1497 Int_t numberOfcells=0;
1498 AliVCluster *cluster = 0;
1500 Double_t *photonRhoArr=new Double_t[ntrig-1];
1501 fhPhotonPtMostEne->Fill(maxPt);
1502 // fhPhotonIndexMostEne->Fill(indexMaxPt);
1503 fhPhotonAverageEnergy->Fill(sumPt/ntrig);
1504 fhPhotonRatioAveEneToMostEne->Fill(sumPt/(ntrig*maxPt));
1505 fhPhotonAverageEnergyMinus1->Fill((sumPt-maxPt)/(ntrig-1));
1506 fGamAvEne=(sumPt-maxPt)/(ntrig-1);
1507 fhPhotonRatioAveEneMinus1ToMostEne->Fill((sumPt-maxPt)/((ntrig-1)*maxPt));
1509 Int_t counterGamma=0;
1510 Int_t counterGammaMinus1=0;
1512 Int_t photonRhoArrayIndex=0;
1513 //TObjArray* clusterstmp = GetEMCALClusters();
1514 for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1515 AliAODPWG4ParticleCorrelation* particlecorr = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1516 if( particlecorr->Pt() > sumPt/ntrig ) counterGamma++;
1517 if( particlecorr->Pt() > (sumPt-maxPt)/(ntrig-1) ) counterGammaMinus1++;
1519 if(iaod==maxIndex) continue;
1520 clusterID = particlecorr->GetCaloLabel(0) ;
1521 if(clusterID < 0) continue;
1522 cluster = FindCluster(clusters,clusterID,iclustmp);
1523 photonRhoArr[photonRhoArrayIndex]=particlecorr->Pt()/ cluster->GetNCells();
1524 numberOfcells+=cluster->GetNCells();
1525 photonRhoArrayIndex++;
1527 if(photonRhoArrayIndex>0) medianPhotonRho=TMath::Median(photonRhoArrayIndex,photonRhoArr);
1528 delete [] photonRhoArr;
1529 fhPhotonNgammaMoreAverageToNgamma->Fill((Double_t)counterGamma / (Double_t)ntrig);
1530 fhPhotonNgammaMoreAverageMinus1ToNgamma->Fill((Double_t)counterGammaMinus1 / (Double_t)ntrig);
1532 //printf("median = %f\n",medianPhotonRho);
1533 fhPhotonBkgRhoVsNtracks->Fill(GetCTSTracks()->GetEntriesFast(),medianPhotonRho);
1534 fhPhotonBkgRhoVsNclusters->Fill(ntrig,medianPhotonRho);
1535 fhPhotonBkgRhoVsCentrality->Fill(GetEventCentrality(),medianPhotonRho);
1536 fhPhotonBkgRhoVsNcells->Fill(numberOfcells,medianPhotonRho);
1539 AliVCluster *cluster2 = 0;
1540 Double_t photon2Corrected=0;
1541 Double_t sumPtTmp=0.;
1542 Double_t sumPtCorrectTmp=0.;
1543 AliVTrack* trackTmp = 0x0 ;
1546 for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1547 AliAODPWG4ParticleCorrelation* particlecorr = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1548 clusterID = particlecorr->GetCaloLabel(0) ;
1549 if(clusterID < 0) continue;
1550 cluster = FindCluster(clusters,clusterID,iclustmp);
1551 fhPhotonPt->Fill(particlecorr->Pt());
1552 fhPhotonPtCorrected->Fill(particlecorr->Pt() - cluster->GetNCells() * medianPhotonRho);
1553 fhPhotonPtDiff->Fill(cluster->GetNCells() * medianPhotonRho);
1554 fhPhotonPtDiffVsCentrality->Fill(GetEventCentrality(),cluster->GetNCells() * medianPhotonRho);
1555 fhPhotonPtDiffVsNcells->Fill(numberOfcells,cluster->GetNCells() * medianPhotonRho);
1556 fhPhotonPtDiffVsNtracks->Fill(GetCTSTracks()->GetEntriesFast(),cluster->GetNCells() * medianPhotonRho);
1557 fhPhotonPtDiffVsNclusters->Fill(ntrig,cluster->GetNCells() * medianPhotonRho);
1559 fhPhotonPtCorrectedZoom->Fill(particlecorr->Pt() - cluster->GetNCells() * medianPhotonRho);
1561 //test: sum_pt in the cone 0.3 for each photon
1562 //should be: random fake gamma from MB
1563 //is: each gamma for EMCEGA
1567 for(Int_t iaod2 = 0; iaod2 < ntrig ; iaod2++){
1568 if(iaod==iaod2) continue;
1569 AliAODPWG4ParticleCorrelation* particlecorr2 = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod2));
1570 clusterID = particlecorr2->GetCaloLabel(0) ;
1571 if(clusterID < 0) continue;
1572 cluster2 = FindCluster(clusters,clusterID,iclustmp);
1573 photon2Corrected = particlecorr2->Pt() - cluster2->GetNCells() * medianPhotonRho;
1575 //if(Pt()<0.5) continue; //<<hardcoded here //FIXME
1576 if( TMath::Sqrt((particlecorr->Eta()-particlecorr2->Eta())*(particlecorr->Eta()-particlecorr2->Eta()) +
1577 (particlecorr->Phi()-particlecorr2->Phi())*(particlecorr->Phi()-particlecorr2->Phi()) )<fGammaConeSize ){//if(/*cone is correct*/){
1578 sumPtTmp+= particlecorr2->Pt();
1579 sumPtCorrectTmp+=photon2Corrected;
1582 fhPhotonSumPtInCone->Fill(sumPtTmp);
1583 fhPhotonSumPtCorrectInCone->Fill(sumPtCorrectTmp);
1585 //test: sum_pt in the cone 0.3 for each track
1586 //should be: random fake gamma from MB
1587 //is: each gamma for EMCEGA
1589 for(Int_t ipr = 0;ipr < GetCTSTracks()->GetEntriesFast() ; ipr ++){
1590 trackTmp = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
1591 p3Tmp.SetXYZ(trackTmp->Px(),trackTmp->Py(),trackTmp->Pz());
1592 if( TMath::Sqrt((particlecorr->Eta()-p3Tmp.Eta())*(particlecorr->Eta()-p3Tmp.Eta()) +
1593 (particlecorr->Phi()-p3Tmp.Phi())*(particlecorr->Phi()-p3Tmp.Phi()) )<fGammaConeSize ){
1594 sumPtTmp+=p3Tmp.Pt();
1596 }//end of loop over tracks
1597 fhPhotonSumPtChargedInCone->Fill(sumPtTmp);
1600 //End of Fill temporary photon histograms
1603 // Apply background subtraction for photons
1605 fGamRho = medianPhotonRho;
1606 if(!IsBackgroundSubtractionGamma()) medianPhotonRho=0;
1610 //Get vertex for cluster momentum calculation <<----new here
1612 Double_t vertex[] = {0,0,0} ; //vertex ;
1613 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
1614 GetReader()->GetVertex(vertex);
1615 fZvertex = vertex[2];
1618 //Loop on stored AOD particles, trigger
1620 for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1621 AliAODPWG4ParticleCorrelation* particlecorr = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1623 fhCuts2->Fill(0.,(Double_t)nJets);
1624 if(GetDebug() > 5) printf("OnlyIsolated %d !particlecorr->IsIsolated() %d \n",OnlyIsolated(), !particlecorr->IsIsolated());
1626 if(OnlyIsolated() && !particlecorr->IsIsolated()) continue;
1628 fhCuts2->Fill(1.,nJets);
1634 //Recover the jet correlated, found previously.
1635 AliAODJet * jet = particlecorr->GetJet();
1636 //If correlation not made before, do it now.
1637 if(fMakeCorrelationInHistoMaker){
1638 //Correlate with jets
1639 Int_t ijet = SelectJet(particlecorr,aodRecJets);//input for jets is TClonesArray
1641 if(GetDebug() > 2) printf ("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - Jet with index %d selected \n",ijet);
1642 //jet = event->GetJet(ijet);
1643 jet = dynamic_cast<AliAODJet*>(aodRecJets-> At(ijet));
1645 if(jet) particlecorr->SetRefJet(jet);
1650 if (!jet) continue ;
1652 fhCuts2->Fill(7.,1.);
1654 //check MC genereted information
1655 if(fMCStudies) FindMCgenInfo();
1658 //Fill Correlation Histograms
1660 clusterID = particlecorr->GetCaloLabel(0) ;
1662 cluster = FindCluster(clusters,clusterID,iclustmp);
1663 //fill tree variables
1664 fGamNcells = cluster->GetNCells();
1666 Double_t ptTrig = particlecorr->Pt() - medianPhotonRho * fGamNcells;//<<---changed here
1667 Double_t ptJet = jet->Pt() - rhoEvent * jet->EffectiveAreaCharged();//<<---changed here
1668 Double_t phiTrig = particlecorr->Phi();
1669 Double_t phiJet = jet->Phi();
1670 Double_t etaTrig = particlecorr->Eta();
1671 Double_t etaJet = jet->Eta();
1672 Double_t deltaPhi=phiTrig-phiJet;
1673 if(deltaPhi<0)deltaPhi+=(TMath::Pi()*2.);
1674 //printf("pT trigger %2.3f, pT jet %2.3f, Delta phi %2.3f, Delta eta %2.3f, Delta pT %2.3f, ratio %2.3f \n",
1675 // ptTrig,ptJet, phiJet-phiTrig, etaJet-etaTrig, ptTrig-ptJet, ptJet/ptTrig);
1676 fhDeltaPt ->Fill(ptTrig, ptTrig-ptJet);
1677 // fhDeltaPhi->Fill(ptTrig, phiTrig-phiJet);//need to be shifted by 2pi
1679 fhDeltaPhiCorrect->Fill(ptTrig, deltaPhi);//correct
1681 Double_t deltaPhiCorrect = TMath::Abs( particlecorr->Phi() - jet->Phi() );
1682 if ( deltaPhiCorrect > TMath::Pi() ) deltaPhiCorrect = 2. * TMath::Pi() - deltaPhiCorrect ;
1683 fhDeltaPhi0PiCorrect->Fill(ptTrig, deltaPhiCorrect);
1685 fhDeltaEta->Fill(ptTrig, etaTrig-etaJet);
1686 fhPtRatio ->Fill(ptTrig, ptJet/ptTrig);
1687 fhPt ->Fill(ptTrig, ptJet);
1689 fhSelectedJetPhiVsEta->Fill(phiJet,etaJet);
1690 fhSelectedJetChBkgEnergyVsPtJet->Fill(ptJet,(IsBackgroundJetFromReader()?rhoEvent * jet->EffectiveAreaCharged():jet->ChargedBgEnergy()) );
1691 fhSelectedJetChAreaVsPtJet->Fill(ptJet,jet->EffectiveAreaCharged());
1692 fhSelectedJetNjet->Fill(nJets);
1693 fhSelectedNtracks->Fill(GetCTSTracks()->GetEntriesFast());//to be checked
1694 fhSelectedPhotonNLMVsPt->Fill(ptTrig,particlecorr->GetFiducialArea());
1698 fhSelectedPhotonLambda0VsPt->Fill(ptTrig,-1);
1699 //fill tree variables
1707 // TObjArray* clusters = GetEMCALClusters();
1708 //cluster = FindCluster(clusters,clusterID,iclustmp);
1709 Double_t lambda0=cluster->GetM02();
1710 fhSelectedPhotonLambda0VsPt->Fill(ptTrig,lambda0);
1711 //fill tree variables
1712 fGamLambda0 = lambda0;
1713 fGamTime = cluster->GetTOF();
1714 //fGamNcells = cluster->GetNCells();
1718 TLorentzVector mom ;
1720 //Double_t scalarProduct=0;
1721 //Double_t vectorLength=particlecorr->P();
1722 for(Int_t icalo=0; icalo <clusters->GetEntriesFast(); icalo++){
1723 AliVCluster* calo = (AliVCluster *) clusters->At(icalo);
1724 if(clusterID==calo->GetID()) continue;//the same cluster as trigger
1725 calo->GetMomentum(mom,vertex) ;//Assume that come from vertex in straight line
1726 //printf("min pt %f\n",GetMinPt());
1727 if(mom.Pt()<GetMinPt()) continue; //<<hardcoded here //FIXME 0.5 check if correct
1728 p3Tmp.SetXYZ(mom.Px(),mom.Py(),mom.Pz());
1729 //calculate sum pt in the cone
1730 if( TMath::Sqrt((particlecorr->Eta()-p3Tmp.Eta())*(particlecorr->Eta()-p3Tmp.Eta()) +
1731 (particlecorr->Phi()-p3Tmp.Phi())*(particlecorr->Phi()-p3Tmp.Phi()) )<fGammaConeSize ){
1732 //scalarProduct = particlecorr->Px()*mom.Px() + particlecorr->Py()*mom.Py() + particlecorr->Pz()*mom.Pz();
1733 //scalarProduct/=mom.P();
1734 //scalarProduct/=vectorLength;
1735 //if(scalarProduct>TMath::Cos(0.3)) {//FIXME photon radius
1736 fGamSumPtNeu+=mom.Pt();
1742 //sum pt of charged tracks in the gamma isolation cone
1746 for(itrack = 0; itrack < nCTSTracks ; itrack++){
1747 aodtrack = dynamic_cast <AliAODTrack*>(GetCTSTracks()->At(itrack));
1748 if(!aodtrack) continue;
1749 fhSelectedTrackPhiVsEta->Fill(aodtrack->Phi(),aodtrack->Eta());//fill histogram here
1750 // if(aodtrack->Pt()<0.15) continue;//hardcoded
1751 if(aodtrack->Pt()<fPtThresholdInCone) continue;
1752 if(!aodtrack->IsHybridGlobalConstrainedGlobal()) continue;
1753 if(TMath::Sqrt((particlecorr->Phi() - aodtrack->Phi())*(particlecorr->Phi() - aodtrack->Phi()) +
1754 (particlecorr->Eta() - aodtrack->Eta())*(particlecorr->Eta() - aodtrack->Eta()) ) <fGammaConeSize ) {
1755 fGamSumPtCh+=aodtrack->Pt();
1761 // for(Int_t itrack = 0; itrack < nCTSTracks ; itrack++){
1762 // aodtrack = dynamic_cast <AliAODTrack*>(GetCTSTracks()->At(itrack));
1763 // fhSelectedTrackPhiVsEta->Fill(aodtrack->Phi(),aodtrack->Eta());
1767 // Background Fragmentation function
1769 TVector3 gammaVector,jetVector;
1770 gammaVector.SetXYZ(particlecorr->Px(),particlecorr->Py(),particlecorr->Pz());
1771 jetVector.SetXYZ(jet->Px(),jet->Py(),jet->Pz());
1772 CalculateBkg(gammaVector,jetVector,vertex,1);//jet perp
1773 CalculateBkg(gammaVector,jetVector,vertex,2);//RC
1774 CalculateBkg(gammaVector,jetVector,vertex,3);//mid point
1775 CalculateBkg(gammaVector,jetVector,vertex,4);//gamma perp
1776 //CalculateBkg(gammaVector,jetVector,vertex,5);/test
1777 Double_t angleJetGam = gammaVector.Angle(jetVector);
1778 //printf("angleJetGam %f\n",angleJetGam*180/TMath::Pi());
1781 // Fragmentation function
1783 Float_t rad = 0, pt = 0, eta = 0, phi = 0;
1784 Int_t npartcone = 0;
1789 printf ("fUseJetRefTracks %d\n",fUseJetRefTracks );
1790 printf ("jet->GetRefTracks() %p\n",jet->GetRefTracks());
1791 printf ("GetCTSTracks() %p\n",GetCTSTracks() );
1794 if(!fUseJetRefTracks)
1795 ntracks =GetCTSTracks()->GetEntriesFast();
1796 else //If you want to use jet tracks from JETAN
1797 ntracks = (jet->GetRefTracks())->GetEntriesFast();
1799 if(GetDebug()>3) printf ("ntracks %d\n",ntracks);
1800 AliVTrack* track = 0x0 ;
1801 for(Int_t ipr = 0;ipr < ntracks ; ipr ++ ){
1802 if(!fUseJetRefTracks)
1803 track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
1804 else //If you want to use jet tracks from JETAN
1805 track = (AliVTrack *) ((jet->GetRefTracks())->At(ipr));
1807 p3.SetXYZ(track->Px(),track->Py(),track->Pz());
1811 if(phi < 0) phi+=TMath::TwoPi();
1813 //Check if there is any particle inside cone with pt larger than fPtThreshold
1814 rad = TMath::Sqrt((eta-etaJet)*(eta-etaJet)+ (phi-phiJet)*(phi-phiJet));
1815 if(rad < fConeSize && pt > fPtThresholdInCone){
1816 //printf("charged in jet cone pt %f, phi %f, eta %f, R %f \n",pt,phi,eta,rad);
1818 fhFFz ->Fill(ptTrig, pt/ptTrig);
1819 fhFFxi->Fill(ptTrig, TMath::Log(ptTrig/pt));
1820 fhFFpt->Fill(ptTrig, pt);
1822 //according to jet axis
1823 fhJetFFz ->Fill(ptJet, pt/ptJet);
1824 fhJetFFxi->Fill(ptJet, TMath::Log(ptJet/pt));
1825 fhJetFFpt->Fill(ptJet, pt);
1828 if(TMath::Cos(angleJetGam)<0 && ptJet!=0 && pt!=0 ){
1829 fhJetFFzCor ->Fill(ptJet, -pt*TMath::Cos(angleJetGam)/ptJet);
1830 fhJetFFxiCor->Fill(ptJet, TMath::Log(ptJet/(-pt*TMath::Cos(angleJetGam))));
1834 fhNTracksInCone->Fill(ptTrig, npartcone);
1835 //fill tree here for each photon-jet (isolation required usually)
1838 //fGamLambda0 = ;//filled earlier
1839 fGamNLM = particlecorr->GetFiducialArea();
1840 //fGamSumPtCh = ;//filled earlier
1841 //fGamTime = particlecorr->GetTOF();//filled earlier
1842 //fGamNcells = particlecorr->GetNCells();//filled earlier
1845 //fGamSumPtNeu = ;//filled earlier
1846 //fGamNtracks = ;//filled earlier
1847 //fGamNclusters= ;//filled earlier
1848 //fGamAvEne = ;//filled earlier
1852 fJetBkgChEne = (IsBackgroundJetFromReader()?rhoEvent * jet->EffectiveAreaCharged():jet->ChargedBgEnergy());
1853 fJetArea = jet->EffectiveAreaCharged();
1854 fJetNtracks = (jet->GetRefTracks())->GetEntriesFast();
1856 fNtracks = GetCTSTracks()->GetEntriesFast();
1857 fCentrality = GetEventCentrality();
1858 fIso = particlecorr->IsIsolated();
1862 for(itrack=0;itrack < fJetNtracks;itrack++){
1863 track = (AliVTrack *) ((jet->GetRefTracks())->At(itrack));
1864 if(track->Pt()>1.) nTrk1GeV++;
1865 if(track->Pt()>2.) nTrk2GeV++;
1868 fJetNtracks1 = nTrk1GeV;
1869 fJetNtracks2 = nTrk2GeV;
1871 if(fSaveGJTree) fTreeGJ->Fill();
1872 }//AOD trigger particle loop
1873 if(GetDebug() > 1) printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - End fill histograms \n");
1877 //__________________________________________________________________
1878 void AliAnaParticleJetFinderCorrelation::Print(const Option_t * opt) const
1881 //Print some relevant parameters set for the analysis
1885 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1886 AliAnaCaloTrackCorrBaseClass::Print(" ");
1888 printf("Phi trigger-jet < %3.2f\n", fDeltaPhiMaxCut) ;
1889 printf("Phi trigger-jet > %3.2f\n", fDeltaPhiMinCut) ;
1890 printf("pT Ratio trigger/jet < %3.2f\n", fRatioMaxCut) ;
1891 printf("pT Ratio trigger/jet > %3.2f\n", fRatioMinCut) ;
1892 printf("fConeSize = %3.2f\n", fConeSize) ;
1893 printf("fPtThresholdInCone = %3.2f\n", fPtThresholdInCone) ;
1894 printf("fUseJetRefTracks = %d\n", fUseJetRefTracks) ;
1895 printf("fMakeCorrelationInHistoMaker = %d\n", fMakeCorrelationInHistoMaker) ;
1896 printf("Isolated Trigger? %d\n", fSelectIsolated) ;
1897 printf("Reconstructed jet cone size = %3.2f\n", fJetConeSize) ;
1898 printf("Reconstructed jet minimum pt before background subtraction = %3.2f\n", fJetMinPt) ;
1899 printf("Reconstructed jet minimum area fraction = %3.2f\n", fJetAreaFraction) ;
1901 //if(!fNonStandardJetFromReader){
1902 printf("fJetBranchName = %s\n", fJetBranchName.Data()) ;
1904 if(!fBackgroundJetFromReader){
1905 printf("fBkgJetBranchName = %s\n", fBkgJetBranchName.Data()) ;
1908 printf("Isolation cone size = %3.2f\n", fGammaConeSize) ;
1909 printf("fUseBackgroundSubtractionGamma = %d\n",fUseBackgroundSubtractionGamma);
1910 printf("fSaveGJTree = %d\n",fSaveGJTree);
1911 printf("fMostEnergetic = %d\n",fMostEnergetic);
1912 printf("fMostOpposite = %d\n",fMostOpposite);
1914 printf("fUseHistogramJetBkg = %d\n",fUseHistogramJetBkg);
1915 printf("fUseHistogramTracks = %d\n",fUseHistogramTracks);
1916 printf("fUseHistogramJetTracks = %d\n",fUseHistogramJetTracks);
1917 printf("fMCStudies = %d\n",fMCStudies);
1921 //__________________________________________________________________
1922 void AliAnaParticleJetFinderCorrelation::CalculateBkg(TVector3 gamma, TVector3 jet,Double_t vertex[3],Int_t type=2){
1924 // calculate background for fragmentation function and fill histograms
1925 // 1. 90 degrees from jet axis in random place = perpendicular cone
1926 // 2. Random cone not belonging to jet cone nor photon cone
1927 // 3. In the middle point from jet and photon momentum vectors
1928 // 4. 90 degrees from photon direction in random place = perpendicular cone 2
1931 // implementation of 2 works, 1 and 4 works
1933 Double_t gammaPt = gamma.Pt();
1934 Double_t gammaEta = gamma.Eta();
1935 Double_t gammaPhi = gamma.Phi();
1936 Double_t jetEta = jet.Eta();
1937 Double_t jetPhi = jet.Phi();
1939 //refference direction of background
1940 Double_t refEta=0.,refPhi=0.;
1941 Double_t rad = 0,rad2 = 0.;
1942 if(type==1){//perpendicular to jet axis
1943 //printf("vertex: %f %f %f \n",vertex[0],vertex[1],vertex[2]);
1950 //Z axis vector: [jet.Px(),jet.Py(),jet.Pz()]
1951 Double_t jx=jet.Px();
1952 Double_t jy=jet.Py();
1953 Double_t jz=jet.Pz();
1954 //if(jz==0)printf("problem\n");
1956 Double_t Xx=1.0-vertex[0];
1957 Double_t Xy=-1.0*vertex[1];
1958 Double_t Xz=jx/jz*(vertex[0]-1.)+vertex[1]*jy/jz;
1960 Double_t Yx=jy*Xz-jz*Xy;
1961 Double_t Yy=jz*Xx-jx*Xz;
1962 Double_t Yz=jx*Xy-jy*Xx;
1964 Double_t det = Xx*Yy*jz + Xy*Yz*jx + Xz*Yx*jy - Xx*Yz*jy - Xy*Yx*jz - Xz*Yy*jx;
1965 if(det==0)printf("problem det==0\n");
1970 // Double_t tmpScalar = jx*Xx+jy*Xy+jz*Xz;
1971 // printf("scalar jet.P o X: %f\n",tmpScalar);
1972 // tmpScalar = jet.Px()*Yx+jet.Py()*Yy+jet.Pz()*Yz;
1973 // printf("scalar jet.P o Y: %f\n",tmpScalar);
1974 // tmpScalar = Xx*Yx+Xy*Yy+Xz*Yz;
1975 // printf("scalar X o Y: %f\n",tmpScalar);
1980 refPhi=fGenerator->Rndm()*TMath::Pi()*2.;
1981 //refPhi=fGenerator->Uniform(-1,1)*TMath::Pi();
1982 xVar=TMath::Cos(refPhi);
1983 yVar=TMath::Sin(refPhi);
1984 //yVar=TMath::Cos(-TMath::Pi()/2.+refPhi);
1985 //zVar=0 in new surface frame
1986 detX = xVar*Yy*jz + Xz*yVar*jy - xVar*Yz*jy - Xy*yVar*jz;
1987 detY = Xx*yVar*jz + xVar*Yz*jx - xVar*Yx*jz - Xz*yVar*jx;
1988 detZ = Xy*yVar*jx + xVar*Yx*jy - Xx*yVar*jy - xVar*Yy*jx;
1994 perp.SetXYZ(newX,newY,newZ);
1995 refEta = perp.Eta();
1996 refPhi = perp.Phi();//output in <-pi, pi> range
1997 if(refPhi<0)refPhi+=2*TMath::Pi();
1998 rad = TMath::Sqrt((gammaEta-refEta)*(gammaEta-refEta) + (gammaPhi-refPhi)*(gammaPhi-refPhi));
1999 rad2 = TMath::Sqrt((jetEta-refEta)*(jetEta-refEta) + (jetPhi-refPhi)*(jetPhi-refPhi));
2000 //printf("refEta,refPhi,rad,rad2: %f %f %f %f\n",refEta,refPhi,rad,rad2);
2001 } while (rad<fJetConeSize+fGammaConeSize || rad2<2.*fJetConeSize || TMath::Abs(refEta)>0.9-fJetConeSize);
2002 fhRandomPhiEta[2]->Fill(refPhi,refEta);
2005 else if(type==2){//random cone
2008 refPhi=fGenerator->Rndm()*TMath::Pi()*2.;
2009 refEta=fGenerator->Uniform(-(0.9-fJetConeSize),0.9-fJetConeSize);
2010 rad = TMath::Sqrt((gammaEta-refEta)*(gammaEta-refEta) + (gammaPhi-refPhi)*(gammaPhi-refPhi));
2011 rad2 = TMath::Sqrt((jetEta-refEta)*(jetEta-refEta) + (jetPhi-refPhi)*(jetPhi-refPhi));
2012 //check if reference is not within jet cone or gamma cone +0.4
2013 //example: FF fConeSize=0.4, fJetConeSize=0.4, fIsoGammaConeSize=0.3
2014 } while (rad<fJetConeSize+fGammaConeSize || rad2<2.*fJetConeSize);
2015 //photon:0.7=0.4+0.3; jets:0.8=0.4 +0.4 //rad<fConeSize+fJetConeSize rad2<fConeSize+0.3
2016 fhRandomPhiEta[0]->Fill(refPhi,refEta);
2018 else if(type==4){//perpendicular to photon axis
2024 //Z axis vector: [jet.Px(),jet.Py(),jet.Pz()]
2025 Double_t jx=gamma.Px();
2026 Double_t jy=gamma.Py();
2027 Double_t jz=gamma.Pz();
2028 //if(jz==0)printf("problem\n");
2030 Double_t Xx=1.0-vertex[0];
2031 Double_t Xy=-1.0*vertex[1];
2032 Double_t Xz=jx/jz*(vertex[0]-1.)+vertex[1]*jy/jz;
2034 Double_t Yx=jy*Xz-jz*Xy;
2035 Double_t Yy=jz*Xx-jx*Xz;
2036 Double_t Yz=jx*Xy-jy*Xx;
2038 Double_t det = Xx*Yy*jz + Xy*Yz*jx + Xz*Yx*jy - Xx*Yz*jy - Xy*Yx*jz - Xz*Yy*jx;
2039 if(det==0)printf("problem det==0\n");
2044 // Double_t tmpScalar = jx*Xx+jy*Xy+jz*Xz;
2045 // printf("scalar jet.P o X: %f\n",tmpScalar);
2046 // tmpScalar = jet.Px()*Yx+jet.Py()*Yy+jet.Pz()*Yz;
2047 // printf("scalar jet.P o Y: %f\n",tmpScalar);
2048 // tmpScalar = Xx*Yx+Xy*Yy+Xz*Yz;
2049 // printf("scalar X o Y: %f\n",tmpScalar);
2054 refPhi=fGenerator->Rndm()*TMath::Pi()*2.;
2055 //refPhi=fGenerator->Uniform(-1,1)*TMath::Pi();
2056 xVar=TMath::Cos(refPhi);
2057 yVar=TMath::Sin(refPhi);
2058 //yVar=TMath::Cos(-TMath::Pi()/2.+refPhi);
2059 //zVar=0 in new surface frame
2060 detX = xVar*Yy*jz + Xz*yVar*jy - xVar*Yz*jy - Xy*yVar*jz;
2061 detY = Xx*yVar*jz + xVar*Yz*jx - xVar*Yx*jz - Xz*yVar*jx;
2062 detZ = Xy*yVar*jx + xVar*Yx*jy - Xx*yVar*jy - xVar*Yy*jx;
2068 perp.SetXYZ(newX,newY,newZ);
2069 refEta = perp.Eta();
2070 refPhi = perp.Phi();//output in <-pi, pi> range
2071 if(refPhi<0)refPhi+=2*TMath::Pi();
2072 rad = TMath::Sqrt((gammaEta-refEta)*(gammaEta-refEta) + (gammaPhi-refPhi)*(gammaPhi-refPhi));
2073 rad2 = TMath::Sqrt((jetEta-refEta)*(jetEta-refEta) + (jetPhi-refPhi)*(jetPhi-refPhi));
2074 //printf("refEta,refPhi,rad,rad2: %f %f %f %f\n",refEta,refPhi,rad,rad2);
2075 } while (rad<fJetConeSize+fGammaConeSize || rad2<2.*fJetConeSize || TMath::Abs(refEta)>0.9-fJetConeSize);
2076 fhRandomPhiEta[1]->Fill(refPhi,refEta);
2079 else if(type==3){//mid point
2081 Double_t jx=jet.Px();
2082 Double_t jy=jet.Py();
2083 Double_t jz=jet.Pz();
2084 // if(jz==0)printf("problem\n");
2085 Double_t gx=gamma.Px();
2086 Double_t gy=gamma.Py();
2087 Double_t gz=gamma.Pz();
2089 Double_t cosAlpha=(jx*gx+jy*gy+jz*gz)/(jet.Mag()*gamma.Mag());
2090 Double_t cosinus=TMath::Sqrt((cosAlpha+1.)/2.);
2091 //perpendicular axis
2092 Double_t Zx=gy*jz-gz*jy;
2093 Double_t Zy=gz*jx-gx*jz;
2094 Double_t Zz=gx*jy-gy*jx;
2097 Double_t det = Zx*gy*jz + Zy*gz*jx + Zz*gx*jy - Zz*gy*jx - Zy*gx*jz - Zx*gz*jy;
2103 Double_t detX = -Zy*gz*cosinus +Zz*cosinus*jy + Zz*gy*cosinus - Zy*cosinus*jz;
2104 Double_t detY = Zx*cosinus*jz - Zz*gx*cosinus - Zz*cosinus*jx + Zx*gz*cosinus;
2105 Double_t detZ = -Zx*gy*cosinus + Zy*cosinus*jx + Zy*gx*cosinus - Zx*cosinus*jy;
2113 perp.SetXYZ(newX,newY,newZ);
2114 refEta = perp.Eta();
2115 refPhi = perp.Phi();//output in <-pi, pi> range
2116 if(refPhi<0)refPhi+=2*TMath::Pi();
2117 rad = TMath::Sqrt((gammaEta-refEta)*(gammaEta-refEta) + (gammaPhi-refPhi)*(gammaPhi-refPhi));
2118 rad2 = TMath::Sqrt((jetEta-refEta)*(jetEta-refEta) + (jetPhi-refPhi)*(jetPhi-refPhi));
2119 //printf("refEta,refPhi,rad,rad2: %f %f %f %f\n",refEta,refPhi,rad,rad2);
2121 if (rad<fJetConeSize+fGammaConeSize || rad2<2.*fJetConeSize || TMath::Abs(refEta)>0.9-fJetConeSize) fhRandomPhiEta[3]->Fill(refPhi,refEta);
2123 else if(type==5){//tmp
2124 //printf("vertex: %f %f %f \n",vertex[0],vertex[1],vertex[2]);
2130 //Z axis vector: [jet.Px(),jet.Py(),jet.Pz()]
2131 Double_t jx=jet.Px();
2132 Double_t jy=jet.Py();
2133 Double_t jz=jet.Pz();
2134 // if(jz==0)printf("problem\n");
2136 Double_t Xx=1.0-vertex[0];
2137 Double_t Xy=-1.0*vertex[1];
2138 Double_t Xz=jx/jz*(vertex[0]-1.)+vertex[1]*jy/jz;
2140 Double_t Yx=jy*Xz-jz*Xy;
2141 Double_t Yy=jz*Xx-jx*Xz;
2142 Double_t Yz=jx*Xy-jy*Xx;
2145 Double_t Xlength=TMath::Sqrt(Xx*Xx+Xy*Xy+Xz*Xz);
2146 Double_t Ylength=TMath::Sqrt(Yx*Yx+Yy*Yy+Yz*Yz);
2147 Double_t ratio=Ylength/Xlength;
2152 refPhi=fGenerator->Rndm()*TMath::Pi()*2.;
2153 xVar=TMath::Tan(refPhi)/ratio;
2158 perp.SetXYZ(newX,newY,newZ);
2159 refEta = perp.Eta();
2160 refPhi = perp.Phi();//output in <-pi, pi> range
2161 if(refPhi<0)refPhi+=2*TMath::Pi();
2162 rad = TMath::Sqrt((gammaEta-refEta)*(gammaEta-refEta) + (gammaPhi-refPhi)*(gammaPhi-refPhi));
2163 rad2 = TMath::Sqrt((jetEta-refEta)*(jetEta-refEta) + (jetPhi-refPhi)*(jetPhi-refPhi));
2164 //printf("refEta,refPhi,rad,rad2: %f %f %f %f\n",refEta,refPhi,rad,rad2);
2165 } while (rad<fJetConeSize+fGammaConeSize || rad2<2.*fJetConeSize || TMath::Abs(refEta)>0.9-fJetConeSize);
2166 fhRandomPhiEta[4]->Fill(refPhi,refEta);
2171 //calculate FF in background
2173 ntracks =GetCTSTracks()->GetEntriesFast();
2174 AliVTrack* track = 0x0 ;
2177 Double_t pt = 0, eta = 0, phi = 0;
2178 Int_t npartcone = 0;
2180 for(Int_t ipr = 0;ipr < ntracks ; ipr ++ ){
2181 track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
2182 p3.SetXYZ(track->Px(),track->Py(),track->Pz());
2184 if(pt<fPtThresholdInCone) {//0.150
2185 //printf("problem: track pt < %f MeV/c \n",fPtThresholdInCone);
2190 if(phi < 0) phi+=TMath::TwoPi();
2191 //Check if there is any particle inside cone with pt larger than fPtThreshold
2192 rad = TMath::Sqrt((eta-refEta)*(eta-refEta) + (phi-refPhi)*(phi-refPhi));
2193 if(rad < fConeSize && pt > fPtThresholdInCone){
2194 //printf("charged in jet cone pt %f, phi %f, eta %f, R %f \n",pt,phi,eta,rad);
2197 if(type==1){//perp jet
2198 fhBkgFFz[1] ->Fill(gammaPt, pt/gammaPt);
2199 fhBkgFFxi[1]->Fill(gammaPt, TMath::Log(gammaPt/pt));
2200 fhBkgFFpt[1]->Fill(gammaPt, pt);
2202 else if(type==2){//RC
2203 fhBkgFFz[0] ->Fill(gammaPt, pt/gammaPt);
2204 fhBkgFFxi[0]->Fill(gammaPt, TMath::Log(gammaPt/pt));
2205 fhBkgFFpt[0]->Fill(gammaPt, pt);
2207 else if(type==3){//mid point
2208 fhBkgFFz[3] ->Fill(gammaPt, pt/gammaPt);
2209 fhBkgFFxi[3]->Fill(gammaPt, TMath::Log(gammaPt/pt));
2210 fhBkgFFpt[3]->Fill(gammaPt, pt);
2212 else if(type==4){//perp jet
2213 fhBkgFFz[2] ->Fill(gammaPt, pt/gammaPt);
2214 fhBkgFFxi[2]->Fill(gammaPt, TMath::Log(gammaPt/pt));
2215 fhBkgFFpt[2]->Fill(gammaPt, pt);
2217 else if(type==5){//test
2218 fhBkgFFz[4] ->Fill(gammaPt, pt/gammaPt);
2219 fhBkgFFxi[4]->Fill(gammaPt, TMath::Log(gammaPt/pt));
2220 fhBkgFFpt[4]->Fill(gammaPt, pt);
2225 }//end of loop over tracks
2226 Double_t sumOverTracks=0.;
2227 if(npartcone!=0) sumOverTracks = sumPt/npartcone;
2229 fhBkgNTracksInCone[1]->Fill(gammaPt, npartcone);
2230 fhBkgSumPtInCone[1]->Fill(gammaPt,sumPt);
2231 fhBkgSumPtnTracksInCone[1]->Fill(gammaPt,sumOverTracks);
2234 fhBkgNTracksInCone[0]->Fill(gammaPt, npartcone);
2235 fhBkgSumPtInCone[0]->Fill(gammaPt,sumPt);
2236 fhBkgSumPtnTracksInCone[0]->Fill(gammaPt,sumOverTracks);
2239 fhBkgNTracksInCone[3]->Fill(gammaPt, npartcone);
2240 fhBkgSumPtInCone[3]->Fill(gammaPt,sumPt);
2241 fhBkgSumPtnTracksInCone[3]->Fill(gammaPt,sumOverTracks);
2244 fhBkgNTracksInCone[2]->Fill(gammaPt, npartcone);
2245 fhBkgSumPtInCone[2]->Fill(gammaPt,sumPt);
2246 fhBkgSumPtnTracksInCone[2]->Fill(gammaPt,sumOverTracks);
2249 fhBkgNTracksInCone[4]->Fill(gammaPt, npartcone);
2250 fhBkgSumPtInCone[4]->Fill(gammaPt,sumPt);
2251 fhBkgSumPtnTracksInCone[4]->Fill(gammaPt,sumOverTracks);
2259 //__________________________________________________________________
2260 void AliAnaParticleJetFinderCorrelation::FindMCgenInfo(){
2262 // Find information about photon and (quark or gluon) on generated level
2265 //frequently used variables
2270 //Double_t photonY = -100 ;
2271 //Double_t photonE = -1 ;
2272 Double_t photonPt = -1 ;
2273 Double_t photonPhi = 100 ;
2274 Double_t photonEta = -1 ;
2275 Bool_t inacceptance = kFALSE;
2276 AliAODMCParticle * primTmp = NULL;
2279 Int_t nParticlesInJet=0;
2280 Int_t nChargedParticlesInJet=0;
2281 Int_t nParticlesInJet150=0;
2282 Int_t nChargedParticlesInJet150=0;
2283 Int_t nChargedParticlesInJet150Cone=0;
2285 Double_t eneParticlesInJet=0.;
2286 Double_t eneChargedParticlesInJet=0.;
2287 Double_t eneParticlesInJet150=0.;
2288 Double_t eneChargedParticlesInJet150=0.;
2289 Double_t eneChargedParticlesInJet150Cone=0.;
2291 Double_t pxParticlesInJet=0.;
2292 Double_t pxChargedParticlesInJet=0.;
2293 Double_t pxParticlesInJet150=0.;
2294 Double_t pxChargedParticlesInJet150=0.;
2295 Double_t pxChargedParticlesInJet150Cone=0.;
2297 Double_t pyParticlesInJet=0.;
2298 Double_t pyChargedParticlesInJet=0.;
2299 Double_t pyParticlesInJet150=0.;
2300 Double_t pyChargedParticlesInJet150=0.;
2301 Double_t pyChargedParticlesInJet150Cone=0.;
2303 Double_t etaParticlesInJet=0.;
2304 Double_t etaChargedParticlesInJet=0.;
2305 Double_t etaParticlesInJet150=0.;
2306 Double_t etaChargedParticlesInJet150=0.;
2307 Double_t etaChargedParticlesInJet150Cone=0.;
2309 Double_t phiParticlesInJet=0.;
2310 Double_t phiChargedParticlesInJet=0.;
2311 Double_t phiParticlesInJet150=0.;
2312 Double_t phiChargedParticlesInJet150=0.;
2313 Double_t phiChargedParticlesInJet150Cone=0.;
2315 Double_t ptParticlesInJet=0.;
2316 Double_t ptChargedParticlesInJet=0.;
2317 Double_t ptParticlesInJet150=0.;
2318 Double_t ptChargedParticlesInJet150=0.;
2319 Double_t ptChargedParticlesInJet150Cone=0.;
2321 Double_t coneJet=0.;
2322 Double_t coneChargedJet=0.;
2323 Double_t coneJet150=0.;
2324 Double_t coneChargedJet150=0.;
2326 std::vector<Int_t> jetParticleIndex;
2328 if(GetReader()->ReadStack()) {//ESD
2329 if(GetDebug()>3) printf("I use stack\n");
2331 else if(GetReader()->ReadAODMCParticles()) {//AOD
2332 TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
2335 //index =6 and 7 is hard scattering (jet-quark or photon)
2336 primTmp = (AliAODMCParticle *) mcparticles->At(6);
2337 pdg=primTmp->GetPdgCode();
2338 if(GetDebug()>3) printf("id 6 pdg %d, pt %f ",pdg,primTmp->Pt() );
2339 if(TMath::Abs(pdg)<=6 ||pdg==21) {
2340 fhMCJetOrigin->Fill(pdg);
2343 primTmp = (AliAODMCParticle *) mcparticles->At(7);
2344 pdg=primTmp->GetPdgCode();
2345 if(GetDebug()>3) printf("id 7 pdg %d, pt %f\n",pdg,primTmp->Pt() );
2346 if(TMath::Abs(pdg)<=6 ||pdg==21) {
2347 fhMCJetOrigin->Fill(pdg);
2352 Int_t nprim = mcparticles->GetEntriesFast();
2353 for(Int_t i=0; i < nprim; i++) {
2354 if(GetReader()->AcceptOnlyHIJINGLabels() && !GetReader()->IsHIJINGLabel(i)) continue ;
2355 AliAODMCParticle * prim = (AliAODMCParticle *) mcparticles->At(i);
2356 pdg = prim->GetPdgCode();
2357 mother=prim->GetMother();
2358 //photon=22, gluon=21, quarks=(1,...,6), antiquarks=(-1,...,-6)
2359 if(pdg == 22){//photon
2360 fhMCPhotonCuts->Fill(0);
2361 if(prim->GetStatus()!=1) continue;
2362 fhMCPhotonCuts->Fill(1);
2363 if(GetDebug()>5) printf("id %d, prim %d, physPrim %d, status %d\n",i,prim->IsPrimary(),prim->IsPhysicalPrimary(),prim->GetStatus());
2365 primTmp = (AliAODMCParticle *) mcparticles->At(mother);
2366 mother=primTmp->GetMother();
2368 if(mother<6)continue;
2369 fhMCPhotonCuts->Fill(2);
2370 primTmp = (AliAODMCParticle *) mcparticles->At(mother);
2371 if(primTmp->GetPdgCode()!=22)continue;
2372 fhMCPhotonCuts->Fill(3);
2374 //Get photon kinematics
2375 photonPt = prim->Pt() ;
2376 photonPhi = prim->Phi() ;
2377 if(photonPhi < 0) photonPhi+=TMath::TwoPi();
2378 photonEta = prim->Eta() ;
2379 fhMCPhotonPt->Fill(photonPt);
2380 fhMCPhotonEtaPhi->Fill(photonPhi,photonEta);
2382 //Check if photons hit the Calorimeter
2384 lv.SetPxPyPzE(prim->Px(),prim->Py(),prim->Pz(),prim->E());
2385 inacceptance = kFALSE;
2386 if(GetCaloUtils()->IsEMCALGeoMatrixSet()){
2387 fhMCPhotonCuts->Fill(4);
2389 if(GetEMCALGeometry()) {
2390 GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(prim->Eta(),prim->Phi(),absID);
2391 if(absID >= 0) inacceptance = kTRUE;
2392 if(GetDebug() > 3) printf("In EMCAL Real acceptance? %d\n",inacceptance);
2395 if(GetFiducialCut()->IsInFiducialCut(lv,"EMCAL")) inacceptance = kTRUE ;
2396 if(GetDebug() > 3) printf("In EMCAL fiducial cut acceptance? %d\n",inacceptance);
2398 }else{//no EMCAL nor EMCALGeoMatrixSet
2399 printf("not EMCALGeoMatrix set\n");
2400 }//end of check if EMCAL
2401 if(inacceptance)fhMCPhotonCuts->Fill(5);
2403 printf("Photon Energy %f, Pt %f\n",prim->E(),prim->Pt());
2405 fMCGamEta=photonEta;
2406 fMCGamPhi=photonPhi;
2407 }//end of check if photon
2409 if(prim->GetStatus()!=1) continue;
2411 printf("id %d, prim %d, physPrim %d, status %d, pdg %d, E %f ",i,prim->IsPrimary(),prim->IsPhysicalPrimary(),prim->GetStatus(),prim->GetPdgCode(),prim->E());
2413 primTmp = (AliAODMCParticle *) mcparticles->At(mother);
2414 mother=primTmp->GetMother();
2415 if(GetDebug() > 5) printf("next mother %d ",mother);
2417 if(GetDebug() > 5) printf("\n");
2418 if(mother<6)continue;//soft part
2419 primTmp = (AliAODMCParticle *) mcparticles->At(mother);
2420 pdg=primTmp->GetPdgCode();
2421 if( !(TMath::Abs(pdg)<=6 || pdg==21) ) continue;//origin not hard q or g
2423 //jetParticleIndex.Add(&i);
2424 jetParticleIndex.push_back(i);
2427 eneParticlesInJet+=prim->E();
2428 pxParticlesInJet+=prim->Px();
2429 pyParticlesInJet+=prim->Py();
2430 etaParticlesInJet+=(prim->E()*prim->Eta());
2431 photonPhi = prim->Phi() ;
2432 if(photonPhi < 0) photonPhi+=TMath::TwoPi();
2433 phiParticlesInJet+=(prim->E()*photonPhi);
2435 if(prim->Charge()!=0) {
2436 nChargedParticlesInJet++;
2437 eneChargedParticlesInJet+=prim->E();
2438 pxChargedParticlesInJet+=prim->Px();
2439 pyChargedParticlesInJet+=prim->Py();
2440 etaChargedParticlesInJet+=(prim->E()*prim->Eta());
2441 phiChargedParticlesInJet+=(prim->E()*photonPhi);
2443 if(prim->Pt()>0.150 && TMath::Abs(prim->Eta())<0.9 ) {//in TPC acceptance :)
2444 nParticlesInJet150++;
2445 eneParticlesInJet150+=prim->E();
2446 pxParticlesInJet150+=prim->Px();
2447 pyParticlesInJet150+=prim->Py();
2448 etaParticlesInJet150+=(prim->E()*prim->Eta());
2449 phiParticlesInJet150+=(prim->E()*photonPhi);
2451 if(prim->Charge()!=0 && prim->Pt()>0.150 && TMath::Abs(prim->Eta())<0.9 ) {//in TPC acceptance :)
2452 nChargedParticlesInJet150++;
2453 eneChargedParticlesInJet150+=prim->E();
2454 pxChargedParticlesInJet150+=prim->Px();
2455 pyChargedParticlesInJet150+=prim->Py();
2456 etaChargedParticlesInJet150+=(prim->E()*prim->Eta());
2457 phiChargedParticlesInJet150+=(prim->E()*photonPhi);
2461 }//end of loop over primaries
2463 if(eneParticlesInJet != 0.) {
2464 etaParticlesInJet/=eneParticlesInJet ;
2465 phiParticlesInJet/=eneParticlesInJet ;
2467 if(eneChargedParticlesInJet != 0) {
2468 etaChargedParticlesInJet/=eneChargedParticlesInJet;
2469 phiChargedParticlesInJet/=eneChargedParticlesInJet;
2471 if(eneParticlesInJet150 != 0) {
2472 etaParticlesInJet150/=eneParticlesInJet150;
2473 phiParticlesInJet150/=eneParticlesInJet150;
2475 if(eneChargedParticlesInJet150 != 0) {
2476 etaChargedParticlesInJet150/=eneChargedParticlesInJet150;
2477 phiChargedParticlesInJet150/=eneChargedParticlesInJet150;
2480 ptParticlesInJet=TMath::Sqrt(pxParticlesInJet*pxParticlesInJet+pyParticlesInJet*pyParticlesInJet);
2481 ptChargedParticlesInJet=TMath::Sqrt(pxChargedParticlesInJet*pxChargedParticlesInJet+pyChargedParticlesInJet*pyChargedParticlesInJet);
2482 ptParticlesInJet150=TMath::Sqrt(pxParticlesInJet150*pxParticlesInJet150+pyParticlesInJet150*pyParticlesInJet150);
2483 ptChargedParticlesInJet150=TMath::Sqrt(pxChargedParticlesInJet150*pxChargedParticlesInJet150+pyChargedParticlesInJet150*pyChargedParticlesInJet150);
2485 Double_t distance=0.;
2488 Double_t mostPtCharged=0.;
2489 Int_t mostmostPtChargedId=-1;
2490 std::vector<Int_t>::iterator it;
2491 for( it=jetParticleIndex.begin(); it!=jetParticleIndex.end(); ++it ){
2492 AliAODMCParticle * prim = (AliAODMCParticle *) mcparticles->At(*it);
2495 if(phi < 0) phi+=TMath::TwoPi();
2497 distance=TMath::Sqrt((eta-etaParticlesInJet)*(eta-etaParticlesInJet)+(phi-phiParticlesInJet)*(phi-phiParticlesInJet));
2498 if(distance>coneJet) coneJet=distance;
2500 distance=TMath::Sqrt((eta-etaChargedParticlesInJet)*(eta-etaChargedParticlesInJet)+(phi-phiChargedParticlesInJet)*(phi-phiChargedParticlesInJet));
2501 if(distance>coneChargedJet) coneChargedJet=distance;
2503 distance=TMath::Sqrt((eta-etaParticlesInJet150)*(eta-etaParticlesInJet150)+(phi-phiParticlesInJet150)*(phi-phiParticlesInJet150));
2504 if(distance>coneJet150 && TMath::Abs(eta)<0.9 ) coneJet150=distance;
2506 distance=TMath::Sqrt((eta-etaChargedParticlesInJet150)*(eta-etaChargedParticlesInJet150)+(phi-phiChargedParticlesInJet150)*(phi-phiChargedParticlesInJet150));
2507 if(distance>coneChargedJet150 && TMath::Abs(eta)<0.9) coneChargedJet150=distance;
2509 if(prim->Charge()!=0 && prim->Pt()>0.150 && TMath::Abs(eta)<0.9) {
2510 if(prim->Pt()>mostPtCharged) {
2511 mostPtCharged=prim->Pt();
2512 mostmostPtChargedId=(*it);
2517 if(prim->Charge()!=0 && prim->Pt()>0.150 && TMath::Abs(eta)<0.9) {
2518 nChargedParticlesInJet150Cone++;
2519 eneChargedParticlesInJet150Cone+=prim->E();
2520 pxChargedParticlesInJet150Cone+=prim->Px();
2521 pyChargedParticlesInJet150Cone+=prim->Py();
2522 etaChargedParticlesInJet150Cone+=(prim->E()*eta);
2523 phiChargedParticlesInJet150Cone+=(prim->E()*phi);
2528 }//end of loop over jet particle indexes
2529 if(eneChargedParticlesInJet150Cone != 0) {
2530 etaChargedParticlesInJet150Cone/=eneChargedParticlesInJet150Cone;
2531 phiChargedParticlesInJet150Cone/=eneChargedParticlesInJet150Cone;
2533 ptChargedParticlesInJet150Cone=TMath::Sqrt(pxChargedParticlesInJet150Cone*pxChargedParticlesInJet150Cone+pyChargedParticlesInJet150Cone*pyChargedParticlesInJet150Cone);
2534 if(nChargedParticlesInJet150>0 && nChargedParticlesInJet150Cone<1){//no particles in cone? take the most energetic one
2535 nChargedParticlesInJet150Cone=1;
2536 etaChargedParticlesInJet150Cone=((AliAODMCParticle *)mcparticles->At(mostmostPtChargedId))->Eta();
2537 phiChargedParticlesInJet150Cone=((AliAODMCParticle *)mcparticles->At(mostmostPtChargedId))->Phi();
2538 ptChargedParticlesInJet150Cone=((AliAODMCParticle *)mcparticles->At(mostmostPtChargedId))->Pt();
2542 }//mc array exists and data is MC
2545 jetParticleIndex.clear();
2549 if(GetDebug() > 3) {
2550 printf("cone full %f, charged %f, full150 %f, charged150 %f\n",coneJet,coneChargedJet,coneJet150,coneChargedJet150);
2551 printf("Npart %d, NchPart %d, Npart(pt>150M) %d, NchPart(pt>150M) %d, NchPart(pt>150M)Cone %d\n",nParticlesInJet,nChargedParticlesInJet,nParticlesInJet150,nChargedParticlesInJet150,nChargedParticlesInJet150Cone);
2552 printf("Etot %f, Ech %f, E(pt>150M) %f, Ech(pt>150M) %f\n",eneParticlesInJet,eneChargedParticlesInJet,eneParticlesInJet150,eneChargedParticlesInJet150);
2553 printf("pt %f, ptch %f, pt(pt>150M) %f,ptch(pt>150M) %f,ptch(pt>150M)Cone %f\n",ptParticlesInJet,ptChargedParticlesInJet,ptParticlesInJet150,ptChargedParticlesInJet150,ptChargedParticlesInJet150Cone);
2554 printf("eta/phi tot %f/%f, ch %f/%f, tot150 %f/%f, ch150 %f/%f, ch150cone %f/%f\n",etaParticlesInJet,phiParticlesInJet,etaChargedParticlesInJet,phiChargedParticlesInJet,etaParticlesInJet150,phiParticlesInJet150,etaChargedParticlesInJet150,phiChargedParticlesInJet150,etaChargedParticlesInJet150Cone,phiChargedParticlesInJet150Cone);
2557 if(ptParticlesInJet) fhMCJetRatioChFull->Fill(ptChargedParticlesInJet/ptParticlesInJet);
2558 if(ptChargedParticlesInJet) fhMCJetRatioCh150Ch->Fill(ptChargedParticlesInJet150/ptChargedParticlesInJet);
2560 fhMCJetNPartVsPt ->Fill(ptParticlesInJet,nParticlesInJet);
2561 fhMCJetChNPartVsPt ->Fill(ptChargedParticlesInJet,nChargedParticlesInJet);
2562 fhMCJetNPart150VsPt ->Fill(ptParticlesInJet150,nParticlesInJet150);
2563 fhMCJetChNPart150VsPt->Fill(ptChargedParticlesInJet150,nChargedParticlesInJet150);
2564 fhMCJetChNPart150ConeVsPt->Fill(ptChargedParticlesInJet150Cone,nChargedParticlesInJet150Cone);
2566 fhMCJetEtaPhi->Fill(phiParticlesInJet,etaParticlesInJet);
2567 fhMCJetChEtaPhi->Fill(phiChargedParticlesInJet,etaChargedParticlesInJet);
2568 fhMCJet150EtaPhi->Fill(phiParticlesInJet150,etaParticlesInJet150);
2569 fhMCJetCh150EtaPhi->Fill(phiChargedParticlesInJet150,etaChargedParticlesInJet150);
2570 fhMCJetCh150ConeEtaPhi->Fill(phiChargedParticlesInJet150Cone,etaChargedParticlesInJet150Cone);
2573 fMCJetPt = ptParticlesInJet;
2574 fMCJetChPt = ptChargedParticlesInJet;
2575 fMCJet150Pt = ptParticlesInJet150;
2576 fMCJetCh150Pt = ptChargedParticlesInJet150;
2577 fMCJetNPart = nParticlesInJet;
2578 fMCJetChNPart = nChargedParticlesInJet;
2579 fMCJet150NPart = nParticlesInJet150;
2580 fMCJetCh150NPart = nChargedParticlesInJet150;
2581 fMCJetEta = etaParticlesInJet ;
2582 fMCJetPhi = phiParticlesInJet ;
2583 fMCJetChEta = etaChargedParticlesInJet ;
2584 fMCJetChPhi = phiChargedParticlesInJet ;
2585 fMCJet150Eta = etaParticlesInJet150 ;
2586 fMCJet150Phi = phiParticlesInJet150 ;
2587 fMCJetCh150Eta = etaChargedParticlesInJet150;
2588 fMCJetCh150Phi = phiChargedParticlesInJet150;
2590 fMCJetCh150ConePt = ptChargedParticlesInJet150Cone;
2591 fMCJetCh150ConeNPart = nChargedParticlesInJet150Cone;
2592 fMCJetCh150ConeEta = etaChargedParticlesInJet150Cone;
2593 fMCJetCh150ConePhi = phiChargedParticlesInJet150Cone;
2595 }//end of method FindMCgenInfo