]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/EMCALJetTasks/AliJetResponseMaker.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliJetResponseMaker.cxx
CommitLineData
7549c451 1// $Id$
2949a2ec 2//
3// Emcal jet response matrix maker task.
4//
5// Author: S. Aiola
6
7#include "AliJetResponseMaker.h"
8
2949a2ec 9#include <TClonesArray.h>
2949a2ec 10#include <TH2F.h>
43032ce2 11#include <THnSparse.h>
2949a2ec 12#include <TLorentzVector.h>
13
ca5c29fa 14#include "AliAnalysisManager.h"
2949a2ec 15#include "AliVCluster.h"
16#include "AliVTrack.h"
17#include "AliEmcalJet.h"
787a3c4f 18#include "AliLog.h"
cd6431de 19#include "AliRhoParameter.h"
5be3857d 20#include "AliNamedArrayI.h"
7cd832c7 21#include "AliJetContainer.h"
22#include "AliParticleContainer.h"
23#include "AliClusterContainer.h"
2949a2ec 24
25ClassImp(AliJetResponseMaker)
26
27//________________________________________________________________________
28AliJetResponseMaker::AliJetResponseMaker() :
9239b066 29 AliAnalysisTaskEmcalJet("AliJetResponseMaker", kTRUE),
cd6431de 30 fMatching(kNoMatching),
7f76e479 31 fMatchingPar1(0),
32 fMatchingPar2(0),
a7477843 33 fUseCellsToMatch(kFALSE),
8159f4cd 34 fMinJetMCPt(1),
43032ce2 35 fHistoType(0),
36 fDeltaPtAxis(0),
37 fDeltaEtaDeltaPhiAxis(0),
faf21244 38 fNEFAxis(0),
39 fZAxis(0),
7cd832c7 40 fIsJet1Rho(kFALSE),
41 fIsJet2Rho(kFALSE),
68e03fc3 42 fHistRejectionReason1(0),
43 fHistRejectionReason2(0),
43032ce2 44 fHistJets1(0),
45 fHistJets2(0),
43032ce2 46 fHistMatching(0),
cd6431de 47 fHistJets1PhiEta(0),
48 fHistJets1PtArea(0),
49 fHistJets1CorrPtArea(0),
63fac07f 50 fHistJets1NEFvsPt(0),
51 fHistJets1CEFvsCEFPt(0),
52 fHistJets1ZvsPt(0),
cd6431de 53 fHistJets2PhiEta(0),
54 fHistJets2PtArea(0),
55 fHistJets2CorrPtArea(0),
63fac07f 56 fHistJets2NEFvsPt(0),
57 fHistJets2CEFvsCEFPt(0),
58 fHistJets2ZvsPt(0),
ca5c29fa 59 fHistCommonEnergy1vsJet1Pt(0),
60 fHistCommonEnergy2vsJet2Pt(0),
61 fHistDistancevsJet1Pt(0),
62 fHistDistancevsJet2Pt(0),
6a20534a 63 fHistDistancevsCommonEnergy1(0),
64 fHistDistancevsCommonEnergy2(0),
05077f28 65 fHistCommonEnergy1vsCommonEnergy2(0),
43032ce2 66 fHistDeltaEtaDeltaPhi(0),
c560b734 67 fHistJet2PtOverJet1PtvsJet2Pt(0),
7030f36f 68 fHistJet1PtOverJet2PtvsJet1Pt(0),
7030f36f 69 fHistDeltaPtvsJet1Pt(0),
cfc2ac24 70 fHistDeltaPtvsJet2Pt(0),
4c00a7af 71 fHistDeltaPtOverJet1PtvsJet1Pt(0),
72 fHistDeltaPtOverJet2PtvsJet2Pt(0),
05077f28 73 fHistDeltaPtvsDistance(0),
74 fHistDeltaPtvsCommonEnergy1(0),
75 fHistDeltaPtvsCommonEnergy2(0),
76 fHistDeltaPtvsArea1(0),
77 fHistDeltaPtvsArea2(0),
78 fHistDeltaPtvsDeltaArea(0),
79 fHistJet1PtvsJet2Pt(0),
4c00a7af 80 fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt(0),
81 fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt(0),
43032ce2 82 fHistDeltaCorrPtvsJet1CorrPt(0),
83 fHistDeltaCorrPtvsJet2CorrPt(0),
05077f28 84 fHistDeltaCorrPtvsDistance(0),
85 fHistDeltaCorrPtvsCommonEnergy1(0),
86 fHistDeltaCorrPtvsCommonEnergy2(0),
87 fHistDeltaCorrPtvsArea1(0),
88 fHistDeltaCorrPtvsArea2(0),
89 fHistDeltaCorrPtvsDeltaArea(0),
cd6431de 90 fHistJet1CorrPtvsJet2CorrPt(0),
4c00a7af 91 fHistDeltaMCPtOverJet1MCPtvsJet1MCPt(0),
92 fHistDeltaMCPtOverJet2PtvsJet2Pt(0),
43032ce2 93 fHistDeltaMCPtvsJet1MCPt(0),
94 fHistDeltaMCPtvsJet2Pt(0),
05077f28 95 fHistDeltaMCPtvsDistance(0),
96 fHistDeltaMCPtvsCommonEnergy1(0),
97 fHistDeltaMCPtvsCommonEnergy2(0),
98 fHistDeltaMCPtvsArea1(0),
99 fHistDeltaMCPtvsArea2(0),
100 fHistDeltaMCPtvsDeltaArea(0),
101 fHistJet1MCPtvsJet2Pt(0)
2949a2ec 102{
103 // Default constructor.
4643d2e8 104
a487deae 105 SetMakeGeneralHistograms(kTRUE);
2949a2ec 106}
107
108//________________________________________________________________________
109AliJetResponseMaker::AliJetResponseMaker(const char *name) :
9239b066 110 AliAnalysisTaskEmcalJet(name, kTRUE),
cd6431de 111 fMatching(kNoMatching),
7f76e479 112 fMatchingPar1(0),
113 fMatchingPar2(0),
a7477843 114 fUseCellsToMatch(kFALSE),
8159f4cd 115 fMinJetMCPt(1),
43032ce2 116 fHistoType(0),
117 fDeltaPtAxis(0),
118 fDeltaEtaDeltaPhiAxis(0),
faf21244 119 fNEFAxis(0),
120 fZAxis(0),
7cd832c7 121 fIsJet1Rho(kFALSE),
122 fIsJet2Rho(kFALSE),
68e03fc3 123 fHistRejectionReason1(0),
124 fHistRejectionReason2(0),
43032ce2 125 fHistJets1(0),
126 fHistJets2(0),
43032ce2 127 fHistMatching(0),
cd6431de 128 fHistJets1PhiEta(0),
129 fHistJets1PtArea(0),
130 fHistJets1CorrPtArea(0),
63fac07f 131 fHistJets1NEFvsPt(0),
132 fHistJets1CEFvsCEFPt(0),
133 fHistJets1ZvsPt(0),
cd6431de 134 fHistJets2PhiEta(0),
135 fHistJets2PtArea(0),
136 fHistJets2CorrPtArea(0),
63fac07f 137 fHistJets2NEFvsPt(0),
138 fHistJets2CEFvsCEFPt(0),
139 fHistJets2ZvsPt(0),
ca5c29fa 140 fHistCommonEnergy1vsJet1Pt(0),
141 fHistCommonEnergy2vsJet2Pt(0),
142 fHistDistancevsJet1Pt(0),
143 fHistDistancevsJet2Pt(0),
6a20534a 144 fHistDistancevsCommonEnergy1(0),
145 fHistDistancevsCommonEnergy2(0),
05077f28 146 fHistCommonEnergy1vsCommonEnergy2(0),
43032ce2 147 fHistDeltaEtaDeltaPhi(0),
c560b734 148 fHistJet2PtOverJet1PtvsJet2Pt(0),
7030f36f 149 fHistJet1PtOverJet2PtvsJet1Pt(0),
7030f36f 150 fHistDeltaPtvsJet1Pt(0),
cfc2ac24 151 fHistDeltaPtvsJet2Pt(0),
4c00a7af 152 fHistDeltaPtOverJet1PtvsJet1Pt(0),
153 fHistDeltaPtOverJet2PtvsJet2Pt(0),
05077f28 154 fHistDeltaPtvsDistance(0),
155 fHistDeltaPtvsCommonEnergy1(0),
156 fHistDeltaPtvsCommonEnergy2(0),
157 fHistDeltaPtvsArea1(0),
158 fHistDeltaPtvsArea2(0),
159 fHistDeltaPtvsDeltaArea(0),
160 fHistJet1PtvsJet2Pt(0),
4c00a7af 161 fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt(0),
162 fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt(0),
43032ce2 163 fHistDeltaCorrPtvsJet1CorrPt(0),
164 fHistDeltaCorrPtvsJet2CorrPt(0),
05077f28 165 fHistDeltaCorrPtvsDistance(0),
166 fHistDeltaCorrPtvsCommonEnergy1(0),
167 fHistDeltaCorrPtvsCommonEnergy2(0),
168 fHistDeltaCorrPtvsArea1(0),
169 fHistDeltaCorrPtvsArea2(0),
170 fHistDeltaCorrPtvsDeltaArea(0),
cd6431de 171 fHistJet1CorrPtvsJet2CorrPt(0),
4c00a7af 172 fHistDeltaMCPtOverJet1MCPtvsJet1MCPt(0),
173 fHistDeltaMCPtOverJet2PtvsJet2Pt(0),
43032ce2 174 fHistDeltaMCPtvsJet1MCPt(0),
175 fHistDeltaMCPtvsJet2Pt(0),
05077f28 176 fHistDeltaMCPtvsDistance(0),
177 fHistDeltaMCPtvsCommonEnergy1(0),
178 fHistDeltaMCPtvsCommonEnergy2(0),
179 fHistDeltaMCPtvsArea1(0),
180 fHistDeltaMCPtvsArea2(0),
181 fHistDeltaMCPtvsDeltaArea(0),
182 fHistJet1MCPtvsJet2Pt(0)
2949a2ec 183{
184 // Standard constructor.
4643d2e8 185
a487deae 186 SetMakeGeneralHistograms(kTRUE);
2949a2ec 187}
188
189//________________________________________________________________________
190AliJetResponseMaker::~AliJetResponseMaker()
191{
192 // Destructor
193}
194
43032ce2 195
2949a2ec 196//________________________________________________________________________
43032ce2 197void AliJetResponseMaker::AllocateTH2()
2949a2ec 198{
43032ce2 199 // Allocate TH2 histograms.
05077f28 200
cfc2ac24 201 fHistJets1PhiEta = new TH2F("fHistJets1PhiEta", "fHistJets1PhiEta", 40, -1, 1, 40, 0, TMath::Pi()*2);
cd6431de 202 fHistJets1PhiEta->GetXaxis()->SetTitle("#eta");
203 fHistJets1PhiEta->GetYaxis()->SetTitle("#phi");
204 fOutput->Add(fHistJets1PhiEta);
2949a2ec 205
faf21244 206 fHistJets1PtArea = new TH2F("fHistJets1PtArea", "fHistJets1PtArea", fNbins/2, 0, 1.5, fNbins, fMinBinPt, fMaxBinPt);
cd6431de 207 fHistJets1PtArea->GetXaxis()->SetTitle("area");
208 fHistJets1PtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
209 fHistJets1PtArea->GetZaxis()->SetTitle("counts");
210 fOutput->Add(fHistJets1PtArea);
43032ce2 211
325e2b7b 212 if (fIsJet1Rho) {
9adcb46d 213 fHistJets1CorrPtArea = new TH2F("fHistJets1CorrPtArea", "fHistJets1CorrPtArea",
faf21244 214 fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
05077f28 215 fHistJets1CorrPtArea->GetXaxis()->SetTitle("area");
216 fHistJets1CorrPtArea->GetYaxis()->SetTitle("p_{T,1}^{corr} (GeV/c)");
217 fHistJets1CorrPtArea->GetZaxis()->SetTitle("counts");
218 fOutput->Add(fHistJets1CorrPtArea);
cd6431de 219 }
220
63fac07f 221 fHistJets1ZvsPt = new TH2F("fHistJets1ZvsPt", "fHistJets1ZvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
222 fHistJets1ZvsPt->GetXaxis()->SetTitle("Z");
223 fHistJets1ZvsPt->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
224 fHistJets1ZvsPt->GetZaxis()->SetTitle("counts");
225 fOutput->Add(fHistJets1ZvsPt);
226
227 fHistJets1NEFvsPt = new TH2F("fHistJets1NEFvsPt", "fHistJets1NEFvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
228 fHistJets1NEFvsPt->GetXaxis()->SetTitle("NEF");
229 fHistJets1NEFvsPt->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
230 fHistJets1NEFvsPt->GetZaxis()->SetTitle("counts");
231 fOutput->Add(fHistJets1NEFvsPt);
232
233 fHistJets1CEFvsCEFPt = new TH2F("fHistJets1CEFvsCEFPt", "fHistJets1CEFvsCEFPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
234 fHistJets1CEFvsCEFPt->GetXaxis()->SetTitle("1-NEF");
235 fHistJets1CEFvsCEFPt->GetYaxis()->SetTitle("(1-NEF)*p_{T,1} (GeV/c)");
236 fHistJets1CEFvsCEFPt->GetZaxis()->SetTitle("counts");
237 fOutput->Add(fHistJets1CEFvsCEFPt);
238
05077f28 239 // Jets 2 histograms
cfc2ac24 240
43032ce2 241 fHistJets2PhiEta = new TH2F("fHistJets2PhiEta", "fHistJets2PhiEta", 40, -1, 1, 40, 0, TMath::Pi()*2);
242 fHistJets2PhiEta->GetXaxis()->SetTitle("#eta");
243 fHistJets2PhiEta->GetYaxis()->SetTitle("#phi");
244 fOutput->Add(fHistJets2PhiEta);
cd6431de 245
faf21244 246 fHistJets2PtArea = new TH2F("fHistJets2PtArea", "fHistJets2PtArea", fNbins/2, 0, 1.5, fNbins, fMinBinPt, fMaxBinPt);
43032ce2 247 fHistJets2PtArea->GetXaxis()->SetTitle("area");
248 fHistJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
249 fHistJets2PtArea->GetZaxis()->SetTitle("counts");
250 fOutput->Add(fHistJets2PtArea);
251
325e2b7b 252 if (fIsJet2Rho) {
43032ce2 253 fHistJets2CorrPtArea = new TH2F("fHistJets2CorrPtArea", "fHistJets2CorrPtArea",
faf21244 254 fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
43032ce2 255 fHistJets2CorrPtArea->GetXaxis()->SetTitle("area");
256 fHistJets2CorrPtArea->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)");
257 fHistJets2CorrPtArea->GetZaxis()->SetTitle("counts");
258 fOutput->Add(fHistJets2CorrPtArea);
259 }
ca5c29fa 260
63fac07f 261 fHistJets2ZvsPt = new TH2F("fHistJets2ZvsPt", "fHistJets2ZvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
262 fHistJets2ZvsPt->GetXaxis()->SetTitle("Z");
263 fHistJets2ZvsPt->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
264 fHistJets2ZvsPt->GetZaxis()->SetTitle("counts");
265 fOutput->Add(fHistJets2ZvsPt);
266
267 fHistJets2NEFvsPt = new TH2F("fHistJets2NEFvsPt", "fHistJets2NEFvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
268 fHistJets2NEFvsPt->GetXaxis()->SetTitle("NEF");
269 fHistJets2NEFvsPt->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
270 fHistJets2NEFvsPt->GetZaxis()->SetTitle("counts");
271 fOutput->Add(fHistJets2NEFvsPt);
272
273 fHistJets2CEFvsCEFPt = new TH2F("fHistJets2CEFvsCEFPt", "fHistJets2CEFvsCEFPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
274 fHistJets2CEFvsCEFPt->GetXaxis()->SetTitle("1-NEF");
275 fHistJets2CEFvsCEFPt->GetYaxis()->SetTitle("(1-NEF)*p_{T,2} (GeV/c)");
276 fHistJets2CEFvsCEFPt->GetZaxis()->SetTitle("counts");
277 fOutput->Add(fHistJets2CEFvsCEFPt);
278
05077f28 279 // Matching histograms
280
ca5c29fa 281 fHistCommonEnergy1vsJet1Pt = new TH2F("fHistCommonEnergy1vsJet1Pt", "fHistCommonEnergy1vsJet1Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt);
282 fHistCommonEnergy1vsJet1Pt->GetXaxis()->SetTitle("Common energy 1 (%)");
283 fHistCommonEnergy1vsJet1Pt->GetYaxis()->SetTitle("p_{T,1}");
284 fHistCommonEnergy1vsJet1Pt->GetZaxis()->SetTitle("counts");
285 fOutput->Add(fHistCommonEnergy1vsJet1Pt);
286
287 fHistCommonEnergy2vsJet2Pt = new TH2F("fHistCommonEnergy2vsJet2Pt", "fHistCommonEnergy2vsJet2Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt);
288 fHistCommonEnergy2vsJet2Pt->GetXaxis()->SetTitle("Common energy 2 (%)");
289 fHistCommonEnergy2vsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");
290 fHistCommonEnergy2vsJet2Pt->GetZaxis()->SetTitle("counts");
291 fOutput->Add(fHistCommonEnergy2vsJet2Pt);
292
293 fHistDistancevsJet1Pt = new TH2F("fHistDistancevsJet1Pt", "fHistDistancevsJet1Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt);
294 fHistDistancevsJet1Pt->GetXaxis()->SetTitle("Distance");
295 fHistDistancevsJet1Pt->GetYaxis()->SetTitle("p_{T,1}");
296 fHistDistancevsJet1Pt->GetZaxis()->SetTitle("counts");
297 fOutput->Add(fHistDistancevsJet1Pt);
298
299 fHistDistancevsJet2Pt = new TH2F("fHistDistancevsJet2Pt", "fHistDistancevsJet2Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt);
300 fHistDistancevsJet2Pt->GetXaxis()->SetTitle("Distance");
301 fHistDistancevsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");
302 fHistDistancevsJet2Pt->GetZaxis()->SetTitle("counts");
303 fOutput->Add(fHistDistancevsJet2Pt);
cd6431de 304
6a20534a 305 fHistDistancevsCommonEnergy1 = new TH2F("fHistDistancevsCommonEnergy1", "fHistDistancevsCommonEnergy1", fNbins/2, 0, 1.2, fNbins/2, 0, 1.2);
306 fHistDistancevsCommonEnergy1->GetXaxis()->SetTitle("Distance");
307 fHistDistancevsCommonEnergy1->GetYaxis()->SetTitle("Common energy 1 (%)");
308 fHistDistancevsCommonEnergy1->GetZaxis()->SetTitle("counts");
309 fOutput->Add(fHistDistancevsCommonEnergy1);
310
311 fHistDistancevsCommonEnergy2 = new TH2F("fHistDistancevsCommonEnergy2", "fHistDistancevsCommonEnergy2", fNbins/2, 0, 1.2, fNbins/2, 0, 1.2);
312 fHistDistancevsCommonEnergy2->GetXaxis()->SetTitle("Distance");
313 fHistDistancevsCommonEnergy2->GetYaxis()->SetTitle("Common energy 2 (%)");
314 fHistDistancevsCommonEnergy2->GetZaxis()->SetTitle("counts");
315 fOutput->Add(fHistDistancevsCommonEnergy2);
cfc2ac24 316
05077f28 317 fHistCommonEnergy1vsCommonEnergy2 = new TH2F("fHistCommonEnergy1vsCommonEnergy2", "fHistCommonEnergy1vsCommonEnergy2", fNbins/2, 0, 1.2, fNbins/2, 0, 1.2);
318 fHistCommonEnergy1vsCommonEnergy2->GetXaxis()->SetTitle("Common energy 1 (%)");
319 fHistCommonEnergy1vsCommonEnergy2->GetYaxis()->SetTitle("Common energy 2 (%)");
320 fHistCommonEnergy1vsCommonEnergy2->GetZaxis()->SetTitle("counts");
321 fOutput->Add(fHistCommonEnergy1vsCommonEnergy2);
322
43032ce2 323 fHistDeltaEtaDeltaPhi = new TH2F("fHistDeltaEtaDeltaPhi", "fHistDeltaEtaDeltaPhi", fNbins/4, -1, 1, fNbins/4, -TMath::Pi()/2, TMath::Pi()*3/2);
324 fHistDeltaEtaDeltaPhi->GetXaxis()->SetTitle("Common energy 1 (%)");
325 fHistDeltaEtaDeltaPhi->GetYaxis()->SetTitle("Common energy 2 (%)");
326 fHistDeltaEtaDeltaPhi->GetZaxis()->SetTitle("counts");
327 fOutput->Add(fHistDeltaEtaDeltaPhi);
328
7030f36f 329 fHistJet2PtOverJet1PtvsJet2Pt = new TH2F("fHistJet2PtOverJet1PtvsJet2Pt", "fHistJet2PtOverJet1PtvsJet2Pt", fNbins, fMinBinPt, fMaxBinPt, 300, 0, 1.5);
c560b734 330 fHistJet2PtOverJet1PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");
331 fHistJet2PtOverJet1PtvsJet2Pt->GetYaxis()->SetTitle("p_{T,2} / p_{T,1}");
332 fHistJet2PtOverJet1PtvsJet2Pt->GetZaxis()->SetTitle("counts");
333 fOutput->Add(fHistJet2PtOverJet1PtvsJet2Pt);
334
7030f36f 335 fHistJet1PtOverJet2PtvsJet1Pt = new TH2F("fHistJet1PtOverJet2PtvsJet1Pt", "fHistJet1PtOverJet2PtvsJet1Pt", fNbins, fMinBinPt, fMaxBinPt, 300, 0, 1.5);
336 fHistJet1PtOverJet2PtvsJet1Pt->GetXaxis()->SetTitle("p_{T,1}");
337 fHistJet1PtOverJet2PtvsJet1Pt->GetYaxis()->SetTitle("p_{T,1} / p_{T,2}");
338 fHistJet1PtOverJet2PtvsJet1Pt->GetZaxis()->SetTitle("counts");
339 fOutput->Add(fHistJet1PtOverJet2PtvsJet1Pt);
340
9adcb46d 341 fHistDeltaPtvsJet1Pt = new TH2F("fHistDeltaPtvsJet1Pt", "fHistDeltaPtvsJet1Pt",
342 fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
7030f36f 343 fHistDeltaPtvsJet1Pt->GetXaxis()->SetTitle("p_{T,1}");
344 fHistDeltaPtvsJet1Pt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
345 fHistDeltaPtvsJet1Pt->GetZaxis()->SetTitle("counts");
346 fOutput->Add(fHistDeltaPtvsJet1Pt);
347
9adcb46d 348 fHistDeltaPtvsJet2Pt = new TH2F("fHistDeltaPtvsJet2Pt", "fHistDeltaPtvsJet2Pt",
349 fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
cfc2ac24 350 fHistDeltaPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");
351 fHistDeltaPtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
352 fHistDeltaPtvsJet2Pt->GetZaxis()->SetTitle("counts");
353 fOutput->Add(fHistDeltaPtvsJet2Pt);
354
4c00a7af 355 fHistDeltaPtOverJet1PtvsJet1Pt = new TH2F("fHistDeltaPtOverJet1PtvsJet1Pt", "fHistDeltaPtOverJet1PtvsJet1Pt",
356 fNbins, fMinBinPt, fMaxBinPt, fNbins, -5, 5);
357 fHistDeltaPtOverJet1PtvsJet1Pt->GetXaxis()->SetTitle("p_{T,1}");
358 fHistDeltaPtOverJet1PtvsJet1Pt->GetYaxis()->SetTitle("#deltap_{T} / p_{T,1}");
359 fHistDeltaPtOverJet1PtvsJet1Pt->GetZaxis()->SetTitle("counts");
360 fOutput->Add(fHistDeltaPtOverJet1PtvsJet1Pt);
361
362 fHistDeltaPtOverJet2PtvsJet2Pt = new TH2F("fHistDeltaPtOverJet2PtvsJet2Pt", "fHistDeltaPtOverJet2PtvsJet2Pt",
363 fNbins, fMinBinPt, fMaxBinPt, fNbins, -5, 5);
364 fHistDeltaPtOverJet2PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");
365 fHistDeltaPtOverJet2PtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} / p_{T,2}");
366 fHistDeltaPtOverJet2PtvsJet2Pt->GetZaxis()->SetTitle("counts");
367 fOutput->Add(fHistDeltaPtOverJet2PtvsJet2Pt);
368
9adcb46d 369 fHistDeltaPtvsDistance = new TH2F("fHistDeltaPtvsDistance", "fHistDeltaPtvsDistance",
370 fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
05077f28 371 fHistDeltaPtvsDistance->GetXaxis()->SetTitle("Distance");
372 fHistDeltaPtvsDistance->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
373 fHistDeltaPtvsDistance->GetZaxis()->SetTitle("counts");
374 fOutput->Add(fHistDeltaPtvsDistance);
375
9adcb46d 376 fHistDeltaPtvsCommonEnergy1 = new TH2F("fHistDeltaPtvsCommonEnergy1", "fHistDeltaPtvsCommonEnergy1",
377 fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
05077f28 378 fHistDeltaPtvsCommonEnergy1->GetXaxis()->SetTitle("Common energy 1 (%)");
379 fHistDeltaPtvsCommonEnergy1->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
380 fHistDeltaPtvsCommonEnergy1->GetZaxis()->SetTitle("counts");
381 fOutput->Add(fHistDeltaPtvsCommonEnergy1);
382
9adcb46d 383 fHistDeltaPtvsCommonEnergy2 = new TH2F("fHistDeltaPtvsCommonEnergy2", "fHistDeltaPtvsCommonEnergy2",
384 fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
05077f28 385 fHistDeltaPtvsCommonEnergy2->GetXaxis()->SetTitle("Common energy 2 (%)");
386 fHistDeltaPtvsCommonEnergy2->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
387 fHistDeltaPtvsCommonEnergy2->GetZaxis()->SetTitle("counts");
388 fOutput->Add(fHistDeltaPtvsCommonEnergy2);
389
9adcb46d 390 fHistDeltaPtvsArea1 = new TH2F("fHistDeltaPtvsArea1", "fHistDeltaPtvsArea1",
faf21244 391 fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
05077f28 392 fHistDeltaPtvsArea1->GetXaxis()->SetTitle("A_{jet,1}");
393 fHistDeltaPtvsArea1->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
394 fHistDeltaPtvsArea1->GetZaxis()->SetTitle("counts");
395 fOutput->Add(fHistDeltaPtvsArea1);
396
9adcb46d 397 fHistDeltaPtvsArea2 = new TH2F("fHistDeltaPtvsArea2", "fHistDeltaPtvsArea2",
faf21244 398 fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
05077f28 399 fHistDeltaPtvsArea2->GetXaxis()->SetTitle("A_{jet,2}");
400 fHistDeltaPtvsArea2->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
401 fHistDeltaPtvsArea2->GetZaxis()->SetTitle("counts");
402 fOutput->Add(fHistDeltaPtvsArea2);
403
404 fHistDeltaPtvsDeltaArea = new TH2F("fHistDeltaPtvsDeltaArea", "fHistDeltaPtvsDeltaArea",
43032ce2 405 fNbins, -1.+1./fNbins, 1.+1./fNbins, 2*fNbins, -fMaxBinPt, fMaxBinPt);
05077f28 406 fHistDeltaPtvsDeltaArea->GetXaxis()->SetTitle("#deltaA_{jet}");
407 fHistDeltaPtvsDeltaArea->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
408 fHistDeltaPtvsDeltaArea->GetZaxis()->SetTitle("counts");
409 fOutput->Add(fHistDeltaPtvsDeltaArea);
410
411 fHistJet1PtvsJet2Pt = new TH2F("fHistJet1PtvsJet2Pt", "fHistJet1PtvsJet2Pt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
412 fHistJet1PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,1}");
413 fHistJet1PtvsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");
414 fHistJet1PtvsJet2Pt->GetZaxis()->SetTitle("counts");
415 fOutput->Add(fHistJet1PtvsJet2Pt);
cd6431de 416
7cd832c7 417 if (fIsJet1Rho || fIsJet2Rho) {
418 if (!fIsJet1Rho)
5c045cde 419 fHistDeltaCorrPtvsJet1CorrPt = new TH2F("fHistDeltaCorrPtvsJet1CorrPt", "fHistDeltaCorrPtvsJet1CorrPt",
420 fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
421 else
422 fHistDeltaCorrPtvsJet1CorrPt = new TH2F("fHistDeltaCorrPtvsJet1CorrPt", "fHistDeltaCorrPtvsJet1CorrPt",
423 2*fNbins, -fMaxBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
424
4c00a7af 425 fHistDeltaCorrPtvsJet1CorrPt->GetXaxis()->SetTitle("p_{T,1}^{corr}");
426 fHistDeltaCorrPtvsJet1CorrPt->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
427 fHistDeltaCorrPtvsJet1CorrPt->GetZaxis()->SetTitle("counts");
428 fOutput->Add(fHistDeltaCorrPtvsJet1CorrPt);
429
7cd832c7 430 if (!fIsJet2Rho)
5c045cde 431 fHistDeltaCorrPtvsJet2CorrPt = new TH2F("fHistDeltaCorrPtvsJet2CorrPt", "fHistDeltaCorrPtvsJet2CorrPt",
432 fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
433 else
434 fHistDeltaCorrPtvsJet2CorrPt = new TH2F("fHistDeltaCorrPtvsJet2CorrPt", "fHistDeltaCorrPtvsJet2CorrPt",
435 2*fNbins, -fMaxBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
436
4c00a7af 437 fHistDeltaCorrPtvsJet2CorrPt->GetXaxis()->SetTitle("p_{T,2}^{corr}");
438 fHistDeltaCorrPtvsJet2CorrPt->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
439 fHistDeltaCorrPtvsJet2CorrPt->GetZaxis()->SetTitle("counts");
440 fOutput->Add(fHistDeltaCorrPtvsJet2CorrPt);
441
442 fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt = new TH2F("fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt", "fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt",
443 2*fNbins, -fMaxBinPt, fMaxBinPt, fNbins, -5, 5);
444 fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt->GetXaxis()->SetTitle("p_{T,1}^{corr}");
445 fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt->GetYaxis()->SetTitle("#deltap_{T}^{corr} / p_{T,1}^{corr}");
446 fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt->GetZaxis()->SetTitle("counts");
447 fOutput->Add(fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt);
448
449 fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt = new TH2F("fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt", "fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt",
450 2*fNbins, -fMaxBinPt, fMaxBinPt, fNbins, -5, 5);
451 fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt->GetXaxis()->SetTitle("p_{T,2}^{corr}");
452 fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt->GetYaxis()->SetTitle("#deltap_{T}^{corr} / p_{T,2}^{corr}");
453 fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt->GetZaxis()->SetTitle("counts");
454 fOutput->Add(fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt);
cfc2ac24 455
9adcb46d 456 fHistDeltaCorrPtvsDistance = new TH2F("fHistDeltaCorrPtvsDistance", "fHistDeltaCorrPtvsDistance",
457 fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
05077f28 458 fHistDeltaCorrPtvsDistance->GetXaxis()->SetTitle("Distance");
459 fHistDeltaCorrPtvsDistance->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
460 fHistDeltaCorrPtvsDistance->GetZaxis()->SetTitle("counts");
461 fOutput->Add(fHistDeltaCorrPtvsDistance);
462
9adcb46d 463 fHistDeltaCorrPtvsCommonEnergy1 = new TH2F("fHistDeltaCorrPtvsCommonEnergy1", "fHistDeltaCorrPtvsCommonEnergy1",
464 fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
05077f28 465 fHistDeltaCorrPtvsCommonEnergy1->GetXaxis()->SetTitle("Common energy 1 (%)");
466 fHistDeltaCorrPtvsCommonEnergy1->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
467 fHistDeltaCorrPtvsCommonEnergy1->GetZaxis()->SetTitle("counts");
468 fOutput->Add(fHistDeltaCorrPtvsCommonEnergy1);
469
9adcb46d 470 fHistDeltaCorrPtvsCommonEnergy2 = new TH2F("fHistDeltaCorrPtvsCommonEnergy2", "fHistDeltaCorrPtvsCommonEnergy2",
471 fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
05077f28 472 fHistDeltaCorrPtvsCommonEnergy2->GetXaxis()->SetTitle("Common energy 2 (%)");
473 fHistDeltaCorrPtvsCommonEnergy2->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
474 fHistDeltaCorrPtvsCommonEnergy2->GetZaxis()->SetTitle("counts");
475 fOutput->Add(fHistDeltaCorrPtvsCommonEnergy2);
476
9adcb46d 477 fHistDeltaCorrPtvsArea1 = new TH2F("fHistDeltaCorrPtvsArea1", "fHistDeltaCorrPtvsArea1",
faf21244 478 fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
05077f28 479 fHistDeltaCorrPtvsArea1->GetXaxis()->SetTitle("A_{jet,1}");
480 fHistDeltaCorrPtvsArea1->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
481 fHistDeltaCorrPtvsArea1->GetZaxis()->SetTitle("counts");
482 fOutput->Add(fHistDeltaCorrPtvsArea1);
483
9adcb46d 484 fHistDeltaCorrPtvsArea2 = new TH2F("fHistDeltaCorrPtvsArea2", "fHistDeltaCorrPtvsArea2",
faf21244 485 fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
05077f28 486 fHistDeltaCorrPtvsArea2->GetXaxis()->SetTitle("A_{jet,2}");
487 fHistDeltaCorrPtvsArea2->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
488 fHistDeltaCorrPtvsArea2->GetZaxis()->SetTitle("counts");
489 fOutput->Add(fHistDeltaCorrPtvsArea2);
490
491 fHistDeltaCorrPtvsDeltaArea = new TH2F("fHistDeltaCorrPtvsDeltaArea", "fHistDeltaCorrPtvsDeltaArea",
43032ce2 492 fNbins, -1.+1./fNbins, 1.+1./fNbins, 2*fNbins, -fMaxBinPt, fMaxBinPt);
05077f28 493 fHistDeltaCorrPtvsDeltaArea->GetXaxis()->SetTitle("#deltaA_{jet}");
494 fHistDeltaCorrPtvsDeltaArea->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
495 fHistDeltaCorrPtvsDeltaArea->GetZaxis()->SetTitle("counts");
496 fOutput->Add(fHistDeltaCorrPtvsDeltaArea);
497
7cd832c7 498 if (!fIsJet1Rho)
9adcb46d 499 fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt",
500 fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
7cd832c7 501 else if (!fIsJet2Rho)
9adcb46d 502 fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt",
503 2*fNbins, -fMaxBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
cd6431de 504 else
9adcb46d 505 fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt",
506 2*fNbins, -fMaxBinPt, fMaxBinPt,
507 2*fNbins, -fMaxBinPt, fMaxBinPt);
cd6431de 508 fHistJet1CorrPtvsJet2CorrPt->GetXaxis()->SetTitle("p_{T,1}^{corr}");
509 fHistJet1CorrPtvsJet2CorrPt->GetYaxis()->SetTitle("p_{T,2}^{corr}");
510 fHistJet1CorrPtvsJet2CorrPt->GetZaxis()->SetTitle("counts");
511 fOutput->Add(fHistJet1CorrPtvsJet2CorrPt);
512 }
513
4358e58a 514 if (fIsEmbedded) {
4c00a7af 515 fHistDeltaMCPtvsJet1MCPt = new TH2F("fHistDeltaMCPtvsJet1MCPt", "fHistDeltaMCPtvsJet1MCPt",
43032ce2 516 fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
4c00a7af 517 fHistDeltaMCPtvsJet1MCPt->GetXaxis()->SetTitle("p_{T,1}^{MC}");
518 fHistDeltaMCPtvsJet1MCPt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
519 fHistDeltaMCPtvsJet1MCPt->GetZaxis()->SetTitle("counts");
520 fOutput->Add(fHistDeltaMCPtvsJet1MCPt);
05077f28 521
9adcb46d 522 fHistDeltaMCPtvsJet2Pt = new TH2F("fHistDeltaMCPtvsJet2Pt", "fHistDeltaMCPtvsJet2Pt",
523 fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
05077f28 524 fHistDeltaMCPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");
525 fHistDeltaMCPtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
526 fHistDeltaMCPtvsJet2Pt->GetZaxis()->SetTitle("counts");
527 fOutput->Add(fHistDeltaMCPtvsJet2Pt);
528
4c00a7af 529 fHistDeltaMCPtOverJet1MCPtvsJet1MCPt = new TH2F("fHistDeltaMCPtOverJet1MCPtvsJet1MCPt", "fHistDeltaMCPtOverJet1MCPtvsJet1MCPt",
530 fNbins, fMinBinPt, fMaxBinPt, fNbins, -5, 5);
531 fHistDeltaMCPtOverJet1MCPtvsJet1MCPt->GetXaxis()->SetTitle("p_{T,1}^{MC}");
532 fHistDeltaMCPtOverJet1MCPtvsJet1MCPt->GetYaxis()->SetTitle("#deltap_{T} / p_{T,1}^{MC}");
533 fHistDeltaMCPtOverJet1MCPtvsJet1MCPt->GetZaxis()->SetTitle("counts");
534 fOutput->Add(fHistDeltaMCPtOverJet1MCPtvsJet1MCPt);
535
536 fHistDeltaMCPtOverJet2PtvsJet2Pt = new TH2F("fHistDeltaMCPtOverJet2PtvsJet2Pt", "fHistDeltaMCPtOverJet2PtvsJet2Pt",
537 fNbins, fMinBinPt, fMaxBinPt, fNbins, -5, 5);
538 fHistDeltaMCPtOverJet2PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");
539 fHistDeltaMCPtOverJet2PtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} / p_{T,2}");
540 fHistDeltaMCPtOverJet2PtvsJet2Pt->GetZaxis()->SetTitle("counts");
541 fOutput->Add(fHistDeltaMCPtOverJet2PtvsJet2Pt);
542
9adcb46d 543 fHistDeltaMCPtvsDistance = new TH2F("fHistDeltaMCPtvsDistance", "fHistDeltaMCPtvsDistance",
544 fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
05077f28 545 fHistDeltaMCPtvsDistance->GetXaxis()->SetTitle("Distance");
546 fHistDeltaMCPtvsDistance->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
547 fHistDeltaMCPtvsDistance->GetZaxis()->SetTitle("counts");
548 fOutput->Add(fHistDeltaMCPtvsDistance);
549
9adcb46d 550 fHistDeltaMCPtvsCommonEnergy1 = new TH2F("fHistDeltaMCPtvsCommonEnergy1", "fHistDeltaMCPtvsCommonEnergy1",
551 fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
05077f28 552 fHistDeltaMCPtvsCommonEnergy1->GetXaxis()->SetTitle("Common energy 1 (%)");
553 fHistDeltaMCPtvsCommonEnergy1->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
554 fHistDeltaMCPtvsCommonEnergy1->GetZaxis()->SetTitle("counts");
555 fOutput->Add(fHistDeltaMCPtvsCommonEnergy1);
556
9adcb46d 557 fHistDeltaMCPtvsCommonEnergy2 = new TH2F("fHistDeltaMCPtvsCommonEnergy2", "fHistDeltaMCPtvsCommonEnergy2",
558 fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
05077f28 559 fHistDeltaMCPtvsCommonEnergy2->GetXaxis()->SetTitle("Common energy 2 (%)");
560 fHistDeltaMCPtvsCommonEnergy2->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
561 fHistDeltaMCPtvsCommonEnergy2->GetZaxis()->SetTitle("counts");
562 fOutput->Add(fHistDeltaMCPtvsCommonEnergy2);
563
9adcb46d 564 fHistDeltaMCPtvsArea1 = new TH2F("fHistDeltaMCPtvsArea1", "fHistDeltaMCPtvsArea1",
faf21244 565 fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
05077f28 566 fHistDeltaMCPtvsArea1->GetXaxis()->SetTitle("A_{jet,1}");
567 fHistDeltaMCPtvsArea1->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
568 fHistDeltaMCPtvsArea1->GetZaxis()->SetTitle("counts");
569 fOutput->Add(fHistDeltaMCPtvsArea1);
570
9adcb46d 571 fHistDeltaMCPtvsArea2 = new TH2F("fHistDeltaMCPtvsArea2", "fHistDeltaMCPtvsArea2",
faf21244 572 fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
05077f28 573 fHistDeltaMCPtvsArea2->GetXaxis()->SetTitle("A_{jet,2}");
574 fHistDeltaMCPtvsArea2->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
575 fHistDeltaMCPtvsArea2->GetZaxis()->SetTitle("counts");
576 fOutput->Add(fHistDeltaMCPtvsArea2);
577
578 fHistDeltaMCPtvsDeltaArea = new TH2F("fHistDeltaMCPtvsDeltaArea", "fHistDeltaMCPtvsDeltaArea",
43032ce2 579 fNbins, -1.+1./fNbins, 1.+1./fNbins, 2*fNbins, -fMaxBinPt, fMaxBinPt);
05077f28 580 fHistDeltaMCPtvsDeltaArea->GetXaxis()->SetTitle("#deltaA_{jet}");
581 fHistDeltaMCPtvsDeltaArea->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
582 fHistDeltaMCPtvsDeltaArea->GetZaxis()->SetTitle("counts");
583 fOutput->Add(fHistDeltaMCPtvsDeltaArea);
584
9adcb46d 585 fHistJet1MCPtvsJet2Pt = new TH2F("fHistJet1MCPtvsJet2Pt", "fHistJet1MCPtvsJet2Pt",
586 fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
4358e58a 587 fHistJet1MCPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,1}^{MC}");
588 fHistJet1MCPtvsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");
589 fHistJet1MCPtvsJet2Pt->GetZaxis()->SetTitle("counts");
590 fOutput->Add(fHistJet1MCPtvsJet2Pt);
591 }
43032ce2 592}
593
594//________________________________________________________________________
595void AliJetResponseMaker::AllocateTHnSparse()
596{
597 // Allocate THnSparse histograms.
598
69a6bda0 599 TString title[25]= {""};
600 Int_t nbins[25] = {0};
601 Double_t min[25] = {0.};
602 Double_t max[25] = {0.};
43032ce2 603 Int_t dim = 0;
604
605 title[dim] = "#phi";
606 nbins[dim] = fNbins/4;
607 min[dim] = 0;
608 max[dim] = 2*TMath::Pi()*(1 + 1./(nbins[dim]-1));
609 dim++;
610
611 title[dim] = "#eta";
612 nbins[dim] = fNbins/4;
613 min[dim] = -1;
614 max[dim] = 1;
615 dim++;
616
617 title[dim] = "p_{T}";
618 nbins[dim] = fNbins;
619 min[dim] = 0;
620 max[dim] = 250;
621 dim++;
622
623 title[dim] = "A_{jet}";
624 nbins[dim] = fNbins/4;
625 min[dim] = 0;
faf21244 626 max[dim] = 1.5;
43032ce2 627 dim++;
628
d3c0e12a 629 if (fNEFAxis) {
630 title[dim] = "NEF";
631 nbins[dim] = fNbins/4;
632 min[dim] = 0;
633 max[dim] = 1.2;
634 dim++;
635 }
43032ce2 636
d3c0e12a 637 if (fZAxis) {
638 title[dim] = "Z";
639 nbins[dim] = fNbins/4;
640 min[dim] = 0;
641 max[dim] = 1.2;
642 dim++;
643 }
644
645 title[dim] = "p_{T,particle}^{leading} (GeV/c)";
646 nbins[dim] = 120;
43032ce2 647 min[dim] = 0;
d3c0e12a 648 max[dim] = 120;
43032ce2 649 dim++;
650
651 Int_t dim1 = dim, dim2 = dim;
652
7cd832c7 653 if (fIsJet1Rho) {
43032ce2 654 title[dim1] = "p_{T}^{corr}";
655 nbins[dim1] = fNbins*2;
656 min[dim1] = -250;
657 max[dim1] = 250;
658 dim1++;
659 }
660
661 if (fIsEmbedded) {
662 title[dim1] = "p_{T}^{MC}";
663 nbins[dim1] = fNbins;
664 min[dim1] = 0;
665 max[dim1] = 250;
666 dim1++;
667 }
668
669 fHistJets1 = new THnSparseD("fHistJets1","fHistJets1",dim1,nbins,min,max);
670 for (Int_t i = 0; i < dim1; i++)
671 fHistJets1->GetAxis(i)->SetTitle(title[i]);
672 fOutput->Add(fHistJets1);
673
7cd832c7 674 if (fIsJet2Rho) {
43032ce2 675 title[dim2] = "p_{T}^{corr}";
676 nbins[dim2] = fNbins*2;
677 min[dim2] = -250;
678 max[dim2] = 250;
679 dim2++;
680 }
681
d3c0e12a 682 fHistJets2 = new THnSparseD("fHistJets2","fHistJets2",dim2,nbins,min,max);
43032ce2 683 for (Int_t i = 0; i < dim2; i++)
d3c0e12a 684 fHistJets2->GetAxis(i)->SetTitle(title[i]);
685 fOutput->Add(fHistJets2);
43032ce2 686
687 // Matching
688
689 dim = 0;
690
691 title[dim] = "p_{T,1}";
692 nbins[dim] = fNbins;
693 min[dim] = 0;
694 max[dim] = 250;
695 dim++;
696
697 title[dim] = "p_{T,2}";
698 nbins[dim] = fNbins;
699 min[dim] = 0;
700 max[dim] = 250;
701 dim++;
702
703 title[dim] = "A_{jet,1}";
704 nbins[dim] = fNbins/4;
705 min[dim] = 0;
faf21244 706 max[dim] = 1.5;
43032ce2 707 dim++;
708
709 title[dim] = "A_{jet,2}";
710 nbins[dim] = fNbins/4;
711 min[dim] = 0;
faf21244 712 max[dim] = 1.5;
43032ce2 713 dim++;
714
715 title[dim] = "distance";
716 nbins[dim] = fNbins/2;
717 min[dim] = 0;
718 max[dim] = 1.2;
719 dim++;
720
721 title[dim] = "CE1";
722 nbins[dim] = fNbins/2;
723 min[dim] = 0;
724 max[dim] = 1.2;
725 dim++;
726
727 title[dim] = "CE2";
728 nbins[dim] = fNbins/2;
729 min[dim] = 0;
730 max[dim] = 1.2;
731 dim++;
732
d3c0e12a 733 title[dim] = "p_{T,particle,1}^{leading} (GeV/c)";
734 nbins[dim] = 120;
735 min[dim] = 0;
736 max[dim] = 120;
737 dim++;
738
739 title[dim] = "p_{T,particle,2}^{leading} (GeV/c)";
740 nbins[dim] = 120;
741 min[dim] = 0;
742 max[dim] = 120;
743 dim++;
744
43032ce2 745 if (fDeltaPtAxis) {
746 title[dim] = "#deltaA_{jet}";
747 nbins[dim] = fNbins/2;
748 min[dim] = -1;
749 max[dim] = 1;
750 dim++;
751
752 title[dim] = "#deltap_{T}";
753 nbins[dim] = fNbins*2;
754 min[dim] = -250;
755 max[dim] = 250;
756 dim++;
757 }
7cd832c7 758 if (fIsJet1Rho) {
43032ce2 759 title[dim] = "p_{T,1}^{corr}";
760 nbins[dim] = fNbins*2;
761 min[dim] = -250;
762 max[dim] = 250;
763 dim++;
764 }
7cd832c7 765 if (fIsJet2Rho) {
43032ce2 766 title[dim] = "p_{T,2}^{corr}";
767 nbins[dim] = fNbins*2;
768 min[dim] = -250;
769 max[dim] = 250;
770 dim++;
771 }
7cd832c7 772 if (fDeltaPtAxis && (fIsJet1Rho || fIsJet2Rho)) {
43032ce2 773 title[dim] = "#deltap_{T}^{corr}";
774 nbins[dim] = fNbins*2;
775 min[dim] = -250;
776 max[dim] = 250;
777 dim++;
778 }
779 if (fDeltaEtaDeltaPhiAxis) {
780 title[dim] = "#delta#eta";
781 nbins[dim] = fNbins/2;
782 min[dim] = -1;
783 max[dim] = 1;
784 dim++;
785
786 title[dim] = "#delta#phi";
787 nbins[dim] = fNbins/2;
788 min[dim] = -TMath::Pi()/2;
789 max[dim] = TMath::Pi()*3/2;
790 dim++;
791 }
792 if (fIsEmbedded) {
793 title[dim] = "p_{T,1}^{MC}";
794 nbins[dim] = fNbins;
795 min[dim] = 0;
796 max[dim] = 250;
797 dim++;
798
799 if (fDeltaPtAxis) {
800 title[dim] = "#deltap_{T}^{MC}";
801 nbins[dim] = fNbins*2;
802 min[dim] = -250;
803 max[dim] = 250;
804 dim++;
805 }
806 }
faf21244 807
808 if (fNEFAxis) {
809 title[dim] = "NEF_{1}";
810 nbins[dim] = fNbins/4;
811 min[dim] = 0;
812 max[dim] = 1.2;
813 dim++;
814
815 title[dim] = "NEF_{2}";
816 nbins[dim] = fNbins/4;
817 min[dim] = 0;
818 max[dim] = 1.2;
819 dim++;
820 }
821
822 if (fZAxis) {
823 title[dim] = "Z_{1}";
824 nbins[dim] = fNbins/4;
825 min[dim] = 0;
826 max[dim] = 1.2;
827 dim++;
828
829 title[dim] = "Z_{2}";
830 nbins[dim] = fNbins/4;
831 min[dim] = 0;
832 max[dim] = 1.2;
833 dim++;
834 }
43032ce2 835
836 fHistMatching = new THnSparseD("fHistMatching","fHistMatching",dim,nbins,min,max);
837
838 for (Int_t i = 0; i < dim; i++)
839 fHistMatching->GetAxis(i)->SetTitle(title[i]);
840
841 fOutput->Add(fHistMatching);
842}
843
844//________________________________________________________________________
845void AliJetResponseMaker::UserCreateOutputObjects()
846{
847 // Create user objects.
848
9239b066 849 AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
43032ce2 850
7cd832c7 851 AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
852 AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
853
854 if (!jets1 || !jets2) return;
855
856 if (jets1->GetRhoName().IsNull()) fIsJet1Rho = kFALSE;
857 else fIsJet1Rho = kTRUE;
858 if (jets2->GetRhoName().IsNull()) fIsJet2Rho = kFALSE;
859 else fIsJet2Rho = kTRUE;
43032ce2 860
68e03fc3 861 fHistRejectionReason1 = new TH2F("fHistRejectionReason1", "fHistRejectionReason1", 32, 0, 32, 100, 0, 250);
862 fHistRejectionReason1->GetXaxis()->SetTitle("Rejection reason");
863 fHistRejectionReason1->GetYaxis()->SetTitle("p_{T,jet} (GeV/c)");
864 fHistRejectionReason1->GetZaxis()->SetTitle("counts");
865 SetRejectionReasonLabels(fHistRejectionReason1->GetXaxis());
866 fOutput->Add(fHistRejectionReason1);
867
868 fHistRejectionReason2 = new TH2F("fHistRejectionReason2", "fHistRejectionReason2", 32, 0, 32, 100, 0, 250);
869 fHistRejectionReason2->GetXaxis()->SetTitle("Rejection reason");
870 fHistRejectionReason2->GetYaxis()->SetTitle("p_{T,jet} (GeV/c)");
871 fHistRejectionReason2->GetZaxis()->SetTitle("counts");
872 SetRejectionReasonLabels(fHistRejectionReason2->GetXaxis());
873 fOutput->Add(fHistRejectionReason2);
874
43032ce2 875 if (fHistoType==0)
876 AllocateTH2();
877 else
878 AllocateTHnSparse();
4358e58a 879
99c67c1b 880 PostData(1, fOutput); // Post data for ALL output slots > 0 here, to get at least an empty histogram
2949a2ec 881}
882
43032ce2 883//________________________________________________________________________
d3c0e12a 884void AliJetResponseMaker::FillJetHisto(Double_t Phi, Double_t Eta, Double_t Pt, Double_t A, Double_t NEF, Double_t LeadingPt, Double_t CorrPt, Double_t MCPt, Int_t Set)
43032ce2 885{
886 if (fHistoType==1) {
887 THnSparse *histo = 0;
888 if (Set==1)
889 histo = fHistJets1;
890 else if (Set==2)
891 histo = fHistJets2;
43032ce2 892
893 if (!histo)
894 return;
895
896 Double_t contents[20]={0};
897
898 for (Int_t i = 0; i < histo->GetNdimensions(); i++) {
899 TString title(histo->GetAxis(i)->GetTitle());
900 if (title=="#phi")
901 contents[i] = Phi;
902 else if (title=="#eta")
903 contents[i] = Eta;
904 else if (title=="p_{T}")
905 contents[i] = Pt;
906 else if (title=="A_{jet}")
907 contents[i] = A;
908 else if (title=="NEF")
909 contents[i] = NEF;
910 else if (title=="Z")
d3c0e12a 911 contents[i] = LeadingPt/Pt;
43032ce2 912 else if (title=="p_{T}^{corr}")
913 contents[i] = CorrPt;
914 else if (title=="p_{T}^{MC}")
915 contents[i] = MCPt;
d3c0e12a 916 else if (title=="p_{T,particle}^{leading} (GeV/c)")
917 contents[i] = LeadingPt;
43032ce2 918 else
919 AliWarning(Form("Unable to fill dimension %s!",title.Data()));
920 }
921
922 histo->Fill(contents);
923 }
924 else {
925 if (Set == 1) {
926 fHistJets1PtArea->Fill(A, Pt);
927 fHistJets1PhiEta->Fill(Eta, Phi);
d3c0e12a 928 fHistJets1ZvsPt->Fill(LeadingPt/Pt, Pt);
43032ce2 929 fHistJets1NEFvsPt->Fill(NEF, Pt);
930 fHistJets1CEFvsCEFPt->Fill(1-NEF, (1-NEF)*Pt);
7cd832c7 931 if (fIsJet1Rho)
43032ce2 932 fHistJets1CorrPtArea->Fill(A, CorrPt);
933 }
934 else if (Set == 2) {
935 fHistJets2PtArea->Fill(A, Pt);
936 fHistJets2PhiEta->Fill(Eta, Phi);
d3c0e12a 937 fHistJets2ZvsPt->Fill(LeadingPt/Pt, Pt);
43032ce2 938 fHistJets2NEFvsPt->Fill(NEF, Pt);
939 fHistJets2CEFvsCEFPt->Fill(1-NEF, (1-NEF)*Pt);
7cd832c7 940 if (fIsJet2Rho)
e8751ce9 941 fHistJets2CorrPtArea->Fill(A, CorrPt);
43032ce2 942 }
943 }
944}
945
946//________________________________________________________________________
faf21244 947void AliJetResponseMaker::FillMatchingHistos(Double_t Pt1, Double_t Pt2, Double_t Eta1, Double_t Eta2, Double_t Phi1, Double_t Phi2,
948 Double_t A1, Double_t A2, Double_t d, Double_t CE1, Double_t CE2, Double_t CorrPt1, Double_t CorrPt2,
d3c0e12a 949 Double_t MCPt1, Double_t NEF1, Double_t NEF2, Double_t LeadingPt1, Double_t LeadingPt2)
43032ce2 950{
951 if (fHistoType==1) {
952 Double_t contents[20]={0};
953
954 for (Int_t i = 0; i < fHistMatching->GetNdimensions(); i++) {
955 TString title(fHistMatching->GetAxis(i)->GetTitle());
956 if (title=="p_{T,1}")
957 contents[i] = Pt1;
958 else if (title=="p_{T,2}")
959 contents[i] = Pt2;
960 else if (title=="A_{jet,1}")
961 contents[i] = A1;
962 else if (title=="A_{jet,2}")
963 contents[i] = A2;
964 else if (title=="distance")
965 contents[i] = d;
966 else if (title=="CE1")
967 contents[i] = CE1;
968 else if (title=="CE2")
969 contents[i] = CE2;
970 else if (title=="#deltaA_{jet}")
971 contents[i] = A1-A2;
972 else if (title=="#deltap_{T}")
973 contents[i] = Pt1-Pt2;
974 else if (title=="#delta#eta")
975 contents[i] = Eta1-Eta2;
976 else if (title=="#delta#phi")
977 contents[i] = Phi1-Phi2;
978 else if (title=="p_{T,1}^{corr}")
979 contents[i] = CorrPt1;
980 else if (title=="p_{T,2}^{corr}")
981 contents[i] = CorrPt2;
982 else if (title=="#deltap_{T}^{corr}")
983 contents[i] = CorrPt1-CorrPt2;
984 else if (title=="p_{T,1}^{MC}")
985 contents[i] = MCPt1;
986 else if (title=="#deltap_{T}^{MC}")
987 contents[i] = MCPt1-Pt2;
faf21244 988 else if (title=="NEF_{1}")
989 contents[i] = NEF1;
990 else if (title=="NEF_{2}")
991 contents[i] = NEF2;
992 else if (title=="Z_{1}")
d3c0e12a 993 contents[i] = LeadingPt1/Pt1;
faf21244 994 else if (title=="Z_{2}")
d3c0e12a 995 contents[i] = LeadingPt2/Pt2;
996 else if (title=="p_{T,particle,1}^{leading} (GeV/c)")
997 contents[i] = LeadingPt1;
998 else if (title=="p_{T,particle,2}^{leading} (GeV/c)")
999 contents[i] = LeadingPt2;
43032ce2 1000 else
1001 AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1002 }
1003
1004 fHistMatching->Fill(contents);
1005 }
1006 else {
1007 fHistCommonEnergy1vsJet1Pt->Fill(CE1, Pt1);
1008 fHistCommonEnergy2vsJet2Pt->Fill(CE2, Pt2);
1009
1010 fHistDistancevsJet1Pt->Fill(d, Pt1);
1011 fHistDistancevsJet2Pt->Fill(d, Pt2);
1012
1013 fHistDistancevsCommonEnergy1->Fill(d, CE1);
1014 fHistDistancevsCommonEnergy2->Fill(d, CE2);
1015 fHistCommonEnergy1vsCommonEnergy2->Fill(CE1, CE2);
1016
1017 fHistDeltaEtaDeltaPhi->Fill(Eta1-Eta2,Phi1-Phi2);
1018
1019 fHistJet2PtOverJet1PtvsJet2Pt->Fill(Pt2, Pt2 / Pt1);
1020 fHistJet1PtOverJet2PtvsJet1Pt->Fill(Pt1, Pt1 / Pt2);
1021
1022 Double_t dpt = Pt1 - Pt2;
1023 fHistDeltaPtvsJet1Pt->Fill(Pt1, dpt);
1024 fHistDeltaPtvsJet2Pt->Fill(Pt2, dpt);
1025 fHistDeltaPtOverJet1PtvsJet1Pt->Fill(Pt1, dpt/Pt1);
1026 fHistDeltaPtOverJet2PtvsJet2Pt->Fill(Pt2, dpt/Pt2);
1027
1028 fHistDeltaPtvsDistance->Fill(d, dpt);
1029 fHistDeltaPtvsCommonEnergy1->Fill(CE1, dpt);
1030 fHistDeltaPtvsCommonEnergy2->Fill(CE2, dpt);
1031
1032 fHistDeltaPtvsArea1->Fill(A1, dpt);
1033 fHistDeltaPtvsArea2->Fill(A2, dpt);
1034
1035 Double_t darea = A1 - A2;
1036 fHistDeltaPtvsDeltaArea->Fill(darea, dpt);
1037
1038 fHistJet1PtvsJet2Pt->Fill(Pt1, Pt2);
1039
7cd832c7 1040 if (fIsJet1Rho || fIsJet2Rho) {
43032ce2 1041 Double_t dcorrpt = CorrPt1 - CorrPt2;
1042 fHistDeltaCorrPtvsJet1CorrPt->Fill(CorrPt1, dcorrpt);
1043 fHistDeltaCorrPtvsJet2CorrPt->Fill(CorrPt2, dcorrpt);
1044 fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt->Fill(CorrPt1, dcorrpt/CorrPt1);
1045 fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt->Fill(CorrPt2, dcorrpt/CorrPt2);
1046 fHistDeltaCorrPtvsDistance->Fill(d, dcorrpt);
1047 fHistDeltaCorrPtvsCommonEnergy1->Fill(CE1, dcorrpt);
1048 fHistDeltaCorrPtvsCommonEnergy2->Fill(CE2, dcorrpt);
1049 fHistDeltaCorrPtvsArea1->Fill(A1, dcorrpt);
1050 fHistDeltaCorrPtvsArea2->Fill(A2, dcorrpt);
1051 fHistDeltaCorrPtvsDeltaArea->Fill(darea, dcorrpt);
1052 fHistJet1CorrPtvsJet2CorrPt->Fill(CorrPt1, CorrPt2);
1053 }
1054
1055 if (fIsEmbedded) {
1056 Double_t dmcpt = MCPt1 - Pt2;
1057 fHistDeltaMCPtvsJet1MCPt->Fill(MCPt1, dmcpt);
1058 fHistDeltaMCPtvsJet2Pt->Fill(Pt2, dmcpt);
1059 fHistDeltaMCPtOverJet1MCPtvsJet1MCPt->Fill(MCPt1, dmcpt/MCPt1);
1060 fHistDeltaMCPtOverJet2PtvsJet2Pt->Fill(Pt2, dmcpt/Pt2);
1061 fHistDeltaMCPtvsDistance->Fill(d, dmcpt);
1062 fHistDeltaMCPtvsCommonEnergy1->Fill(CE1, dmcpt);
1063 fHistDeltaMCPtvsCommonEnergy2->Fill(CE2, dmcpt);
1064 fHistDeltaMCPtvsArea1->Fill(A1, dmcpt);
1065 fHistDeltaMCPtvsArea2->Fill(A2, dmcpt);
1066 fHistDeltaMCPtvsDeltaArea->Fill(darea, dmcpt);
1067 fHistJet1MCPtvsJet2Pt->Fill(MCPt1, Pt2);
1068 }
1069 }
1070}
1071
b3ccdfe2 1072//________________________________________________________________________
1073void AliJetResponseMaker::ExecOnce()
1074{
1075 // Execute once.
1076
9239b066 1077 AliAnalysisTaskEmcalJet::ExecOnce();
cfc2ac24 1078
7cd832c7 1079 AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
1080 AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
cd6431de 1081
7cd832c7 1082 if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return;
cd6431de 1083
7cd832c7 1084 if (fMatching == kMCLabel && (!jets2->GetIsParticleLevel() || jets1->GetIsParticleLevel())) {
1085 if (jets1->GetIsParticleLevel() == jets2->GetIsParticleLevel()) {
43032ce2 1086 AliWarning("Changing matching type from MC label to same collection...");
1087 fMatching = kSameCollections;
1088 }
1089 else {
1090 AliWarning("Changing matching type from MC label to geometrical...");
1091 fMatching = kGeometrical;
1092 }
1093 }
7cd832c7 1094 else if (fMatching == kSameCollections && jets1->GetIsParticleLevel() != jets2->GetIsParticleLevel()) {
1095 if (jets2->GetIsParticleLevel() && !jets1->GetIsParticleLevel()) {
43032ce2 1096 AliWarning("Changing matching type from same collection to MC label...");
1097 fMatching = kMCLabel;
1098 }
1099 else {
1100 AliWarning("Changing matching type from same collection to geometrical...");
1101 fMatching = kGeometrical;
1102 }
1103 }
b3ccdfe2 1104}
1105
1106//________________________________________________________________________
1107Bool_t AliJetResponseMaker::Run()
1108{
1109 // Find the closest jets
cfc2ac24 1110
f660c2d6 1111 if (fMatching == kNoMatching)
aa4d701c 1112 return kTRUE;
f660c2d6 1113 else
1114 return DoJetMatching();
cfc2ac24 1115}
aa4d701c 1116
cfc2ac24 1117//________________________________________________________________________
f660c2d6 1118Bool_t AliJetResponseMaker::DoJetMatching()
cfc2ac24 1119{
7cd832c7 1120 AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
1121 AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
cfc2ac24 1122
7cd832c7 1123 if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return kFALSE;
cfc2ac24 1124
7cd832c7 1125 DoJetLoop();
cfc2ac24 1126
7cd832c7 1127 AliEmcalJet* jet1 = 0;
cfc2ac24 1128
7cd832c7 1129 jets1->ResetCurrentID();
1130 while ((jet1 = jets1->GetNextJet())) {
cfc2ac24 1131
7cd832c7 1132 AliEmcalJet *jet2 = jet1->ClosestJet();
cfc2ac24 1133
7cd832c7 1134 if (!jet2) continue;
1135 if (jet2->ClosestJet() != jet1) continue;
1136 if (jet1->ClosestJetDistance() > fMatchingPar1 || jet2->ClosestJetDistance() > fMatchingPar2) continue;
cfc2ac24 1137
7cd832c7 1138 // Matched jet found
1139 jet1->SetMatchedToClosest(fMatching);
1140 jet2->SetMatchedToClosest(fMatching);
1141 AliDebug(2,Form("Found matching: jet1 pt = %f, eta = %f, phi = %f, jet2 pt = %f, eta = %f, phi = %f",
1142 jet1->Pt(), jet1->Eta(), jet1->Phi(),
1143 jet2->Pt(), jet2->Eta(), jet2->Phi()));
aa4d701c 1144 }
1145
b3ccdfe2 1146 return kTRUE;
1147}
1148
2949a2ec 1149//________________________________________________________________________
7cd832c7 1150void AliJetResponseMaker::DoJetLoop()
2949a2ec 1151{
44476be7 1152 // Do the jet loop.
1ee1b5b8 1153
7cd832c7 1154 AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
1155 AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
fde82e42 1156
7cd832c7 1157 if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return;
fde82e42 1158
7cd832c7 1159 AliEmcalJet* jet1 = 0;
1160 AliEmcalJet* jet2 = 0;
fde82e42 1161
7cd832c7 1162 jets2->ResetCurrentID();
1163 while ((jet2 = jets2->GetNextJet())) jet2->ResetMatching();
05077f28 1164
7cd832c7 1165 jets1->ResetCurrentID();
1166 while ((jet1 = jets1->GetNextJet())) {
2103dc6a 1167 jet1->ResetMatching();
7cd832c7 1168
1169 if (jet1->MCPt() < fMinJetMCPt) continue;
44476be7 1170
7cd832c7 1171 jets2->ResetCurrentID();
1172 while ((jet2 = jets2->GetNextJet())) {
7f76e479 1173 SetMatchingLevel(jet1, jet2, fMatching);
2103dc6a 1174 } // jet2 loop
2103dc6a 1175 } // jet1 loop
1ee1b5b8 1176}
1177
cd6431de 1178//________________________________________________________________________
7f76e479 1179void AliJetResponseMaker::GetGeometricalMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d) const
cd6431de 1180{
7f76e479 1181 Double_t deta = jet2->Eta() - jet1->Eta();
1182 Double_t dphi = jet2->Phi() - jet1->Phi();
1183 d = TMath::Sqrt(deta * deta + dphi * dphi);
1184}
cd6431de 1185
7f76e479 1186//________________________________________________________________________
1187void AliJetResponseMaker::GetMCLabelMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d1, Double_t &d2) const
1188{
7cd832c7 1189 AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
1190 AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
1191
1192 if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return;
1193
1194 AliParticleContainer *tracks1 = jets1->GetParticleContainer();
1195 AliClusterContainer *clusters1 = jets1->GetClusterContainer();
1196 AliParticleContainer *tracks2 = jets2->GetParticleContainer();
1197
7f76e479 1198 // d1 and d2 represent the matching level: 0 = maximum level of matching, 1 = the two jets are completely unrelated
1199 d1 = jet1->Pt();
1200 d2 = jet2->Pt();
1201 Double_t totalPt1 = d1; // the total pt of the reconstructed jet will be cleaned from the background
1202
7cd832c7 1203 // remove completely tracks that are not MC particles (label == 0)
1204 if (tracks1 && tracks1->GetArray()) {
7f76e479 1205 for (Int_t iTrack = 0; iTrack < jet1->GetNumberOfTracks(); iTrack++) {
7cd832c7 1206 AliVParticle *track = jet1->TrackAt(iTrack,tracks1->GetArray());
7f76e479 1207 if (!track) {
1208 AliWarning(Form("Could not find track %d!", iTrack));
1209 continue;
1210 }
7cd832c7 1211
787a3c4f 1212 Int_t MClabel = TMath::Abs(track->GetLabel());
68faaa50 1213 MClabel -= fMCLabelShift;
7cd832c7 1214 if (MClabel != 0) continue;
1215
1216 // this is not a MC particle; remove it completely
1217 AliDebug(3,Form("Track %d (pT = %f) is not a MC particle (MClabel = %d)!",iTrack,track->Pt(),MClabel));
1218 totalPt1 -= track->Pt();
1219 d1 -= track->Pt();
1220 }
1221 }
1222
1223 // remove completely clusters that are not MC particles (label == 0)
1224 if (fUseCellsToMatch && fCaloCells) {
1225 for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
1226 AliVCluster *clus = jet1->ClusterAt(iClus,clusters1->GetArray());
1227 if (!clus) {
1228 AliWarning(Form("Could not find cluster %d!", iClus));
7f76e479 1229 continue;
1230 }
7cd832c7 1231 TLorentzVector part;
1232 clus->GetMomentum(part, const_cast<Double_t*>(fVertex));
1233
1234 for (Int_t iCell = 0; iCell < clus->GetNCells(); iCell++) {
1235 Int_t cellId = clus->GetCellAbsId(iCell);
1236 Double_t cellFrac = clus->GetCellAmplitudeFraction(iCell);
1237
1238 Int_t MClabel = TMath::Abs(fCaloCells->GetCellMCLabel(cellId));
68faaa50 1239 MClabel -= fMCLabelShift;
7cd832c7 1240 if (MClabel != 0) continue;
1241
1242 // this is not a MC particle; remove it completely
1243 AliDebug(3,Form("Cell %d (frac = %f) is not a MC particle (MClabel = %d)!",iCell,cellFrac,MClabel));
1244 totalPt1 -= part.Pt() * cellFrac;
1245 d1 -= part.Pt() * cellFrac;
1246 }
1247 }
1248 }
1249 else {
1250 for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
1251 AliVCluster *clus = jet1->ClusterAt(iClus,clusters1->GetArray());
1252 if (!clus) {
1253 AliWarning(Form("Could not find cluster %d!", iClus));
1254 continue;
7f76e479 1255 }
7cd832c7 1256 TLorentzVector part;
1257 clus->GetMomentum(part, const_cast<Double_t*>(fVertex));
7f76e479 1258
7cd832c7 1259 Int_t MClabel = TMath::Abs(clus->GetLabel());
68faaa50 1260 MClabel -= fMCLabelShift;
7cd832c7 1261 if (MClabel != 0) continue;
1262
1263 // this is not a MC particle; remove it completely
1264 AliDebug(3,Form("Cluster %d (pT = %f) is not a MC particle (MClabel = %d)!",iClus,part.Pt(),MClabel));
1265 totalPt1 -= part.Pt();
1266 d1 -= part.Pt();
1267 }
1268 }
1269
1270 for (Int_t iTrack2 = 0; iTrack2 < jet2->GetNumberOfTracks(); iTrack2++) {
1271 Bool_t track2Found = kFALSE;
1272 Int_t index2 = jet2->TrackAt(iTrack2);
1273
1274 // now look for common particles in the track array
1275 for (Int_t iTrack = 0; iTrack < jet1->GetNumberOfTracks(); iTrack++) {
1276 AliVParticle *track = jet1->TrackAt(iTrack,tracks1->GetArray());
1277 if (!track) {
1278 AliWarning(Form("Could not find track %d!", iTrack));
1279 continue;
1280 }
1281 Int_t MClabel = TMath::Abs(track->GetLabel());
68faaa50 1282 MClabel -= fMCLabelShift;
7cd832c7 1283 if (MClabel <= 0) continue;
1284
1285 Int_t index = -1;
1286 index = tracks2->GetIndexFromLabel(MClabel);
7f76e479 1287 if (index < 0) {
1288 AliDebug(2,Form("Track %d (pT = %f) does not have an associated MC particle (MClabel = %d)!",iTrack,track->Pt(),MClabel));
1289 continue;
1290 }
7cd832c7 1291
1292 if (index2 != index) continue;
1293
1294 // found common particle
1295 d1 -= track->Pt();
7f76e479 1296
7cd832c7 1297 if (!track2Found) {
1298 AliVParticle *MCpart = tracks2->GetParticle(index2);
7f76e479 1299 AliDebug(3,Form("Track %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!",
1300 iTrack,track->Pt(),track->Eta(),track->Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
1301 d2 -= MCpart->Pt();
7f76e479 1302 }
7cd832c7 1303
1304 track2Found = kTRUE;
cd6431de 1305 }
7cd832c7 1306
1307 // now look for common particles in the cluster array
1308 if (fUseCellsToMatch && fCaloCells) { // if the cell colection is available, look for cells with a matched MC particle
1309 for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
1310 AliVCluster *clus = jet1->ClusterAt(iClus,fCaloClusters);
1311 if (!clus) {
1312 AliWarning(Form("Could not find cluster %d!", iClus));
1313 continue;
1314 }
1315 TLorentzVector part;
1316 clus->GetMomentum(part, const_cast<Double_t*>(fVertex));
1317
7f76e479 1318 for (Int_t iCell = 0; iCell < clus->GetNCells(); iCell++) {
1319 Int_t cellId = clus->GetCellAbsId(iCell);
1320 Double_t cellFrac = clus->GetCellAmplitudeFraction(iCell);
7cd832c7 1321
787a3c4f 1322 Int_t MClabel = TMath::Abs(fCaloCells->GetCellMCLabel(cellId));
68faaa50 1323 MClabel -= fMCLabelShift;
7cd832c7 1324 if (MClabel <= 0) continue;
7f76e479 1325
7cd832c7 1326 Int_t index1 = -1;
1327 index1 = tracks2->GetIndexFromLabel(MClabel);
1328 if (index1 < 0) {
7f76e479 1329 AliDebug(3,Form("Cell %d (frac = %f) does not have an associated MC particle (MClabel = %d)!",iCell,cellFrac,MClabel));
1330 continue;
1331 }
7cd832c7 1332
1333 if (index2 != index1) continue;
1334
1335 // found common particle
1336 d1 -= part.Pt() * cellFrac;
1337
1338 if (!track2Found) { // only if it is not already found among charged tracks (charged particles are most likely already found)
1339 AliVParticle *MCpart = tracks2->GetParticle(index2);
1340 AliDebug(3,Form("Cell %d belonging to cluster %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!",
1341 iCell,iClus,part.Pt(),part.Eta(),part.Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
1342 d2 -= MCpart->Pt() * cellFrac;
7f76e479 1343 }
7cd832c7 1344
1345 track2Found = kTRUE;
cfc2ac24 1346 }
1347 }
7cd832c7 1348 }
1349 else { //otherwise look for the first contributor to the cluster, and if matched to a MC label remove it
1350 for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
1351 AliVCluster *clus = jet1->ClusterAt(iClus,fCaloClusters);
1352 if (!clus) {
1353 AliWarning(Form("Could not find cluster %d!", iClus));
cfc2ac24 1354 continue;
1355 }
7cd832c7 1356 TLorentzVector part;
1357 clus->GetMomentum(part, const_cast<Double_t*>(fVertex));
1358
1359 Int_t MClabel = TMath::Abs(clus->GetLabel());
68faaa50 1360 MClabel -= fMCLabelShift;
7cd832c7 1361 if (MClabel <= 0) continue;
1362
1363 Int_t index = -1;
1364 index = tracks2->GetIndexFromLabel(MClabel);
1365
cfc2ac24 1366 if (index < 0) {
7f76e479 1367 AliDebug(3,Form("Cluster %d (pT = %f) does not have an associated MC particle (MClabel = %d)!",iClus,part.Pt(),MClabel));
cfc2ac24 1368 continue;
1369 }
7cd832c7 1370
1371 if (index2 != index) continue;
1372
1373 // found common particle
1374 d1 -= part.Pt();
1375
1376 if (!track2Found) { // only if it is not already found among charged tracks (charged particles are most likely already found)
1377 AliVParticle *MCpart = tracks2->GetParticle(index2);
1378 AliDebug(3,Form("Cluster %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!",
1379 iClus,part.Pt(),part.Eta(),part.Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
1380
1381 d2 -= MCpart->Pt();
cfc2ac24 1382 }
7cd832c7 1383
1384 track2Found = kTRUE;
cfc2ac24 1385 }
7f76e479 1386 }
1387 }
05077f28 1388
1389 if (d1 < 0)
7f76e479 1390 d1 = 0;
05077f28 1391
1392 if (d2 < 0)
1393 d2 = 0;
1394
1395 if (totalPt1 < 1)
1396 d1 = -1;
7f76e479 1397 else
1398 d1 /= totalPt1;
f660c2d6 1399
05077f28 1400 if (jet2->Pt() < 1)
1401 d2 = -1;
7f76e479 1402 else
05077f28 1403 d2 /= jet2->Pt();
7f76e479 1404}
1405
1406//________________________________________________________________________
1407void AliJetResponseMaker::GetSameCollectionsMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d1, Double_t &d2) const
1408{
7cd832c7 1409 AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
1410 AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
1411
1412 if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return;
1413
1414 AliParticleContainer *tracks1 = jets1->GetParticleContainer();
1415 AliClusterContainer *clusters1 = jets1->GetClusterContainer();
1416 AliParticleContainer *tracks2 = jets2->GetParticleContainer();
1417 AliClusterContainer *clusters2 = jets2->GetClusterContainer();
1418
7f76e479 1419 // d1 and d2 represent the matching level: 0 = maximum level of matching, 1 = the two jets are completely unrelated
1420 d1 = jet1->Pt();
1421 d2 = jet2->Pt();
1422
7cd832c7 1423 if (tracks1 && tracks2) {
7f76e479 1424
1425 for (Int_t iTrack2 = 0; iTrack2 < jet2->GetNumberOfTracks(); iTrack2++) {
1426 Int_t index2 = jet2->TrackAt(iTrack2);
7cd832c7 1427 for (Int_t iTrack1 = 0; iTrack1 < jet1->GetNumberOfTracks(); iTrack1++) {
1428 Int_t index1 = jet1->TrackAt(iTrack1);
1429 if (index2 == index1) { // found common particle
1430 AliVParticle *part1 = tracks1->GetParticle(index1);
1431 if (!part1) {
1432 AliWarning(Form("Could not find track %d!", index1));
f660c2d6 1433 continue;
1434 }
7cd832c7 1435 AliVParticle *part2 = tracks2->GetParticle(index2);
7f76e479 1436 if (!part2) {
1437 AliWarning(Form("Could not find track %d!", index2));
cfc2ac24 1438 continue;
1439 }
7f76e479 1440
7cd832c7 1441 d1 -= part1->Pt();
7f76e479 1442 d2 -= part2->Pt();
1443 break;
cfc2ac24 1444 }
1445 }
cfc2ac24 1446 }
7f76e479 1447
1448 }
1449
7cd832c7 1450 if (clusters1 && clusters2) {
7f76e479 1451
a7477843 1452 if (fUseCellsToMatch) {
1453 const Int_t nClus1 = jet1->GetNumberOfClusters();
1454
1455 Int_t ncells1[nClus1];
1456 UShort_t *cellsId1[nClus1];
1457 Double_t *cellsFrac1[nClus1];
1458 Int_t *sortedIndexes1[nClus1];
1459 Double_t ptClus1[nClus1];
1460 for (Int_t iClus1 = 0; iClus1 < nClus1; iClus1++) {
1461 Int_t index1 = jet1->ClusterAt(iClus1);
7cd832c7 1462 AliVCluster *clus1 = clusters1->GetCluster(index1);
a7477843 1463 if (!clus1) {
1464 AliWarning(Form("Could not find cluster %d!", index1));
9adcb46d 1465 ncells1[iClus1] = 0;
1466 cellsId1[iClus1] = 0;
1467 cellsFrac1[iClus1] = 0;
1468 sortedIndexes1[iClus1] = 0;
1469 ptClus1[iClus1] = 0;
a7477843 1470 continue;
1471 }
9adcb46d 1472 TLorentzVector part1;
1473 clus1->GetMomentum(part1, const_cast<Double_t*>(fVertex));
1474
a7477843 1475 ncells1[iClus1] = clus1->GetNCells();
1476 cellsId1[iClus1] = clus1->GetCellsAbsId();
1477 cellsFrac1[iClus1] = clus1->GetCellsAmplitudeFraction();
1478 sortedIndexes1[iClus1] = new Int_t[ncells1[iClus1]];
a7477843 1479 ptClus1[iClus1] = part1.Pt();
1480
1481 TMath::Sort(ncells1[iClus1], cellsId1[iClus1], sortedIndexes1[iClus1], kFALSE);
1482 }
1483
1484 const Int_t nClus2 = jet2->GetNumberOfClusters();
1485
1486 const Int_t maxNcells2 = 11520;
1487 Int_t sortedIndexes2[maxNcells2];
1488 for (Int_t iClus2 = 0; iClus2 < nClus2; iClus2++) {
1489 Int_t index2 = jet2->ClusterAt(iClus2);
7cd832c7 1490 AliVCluster *clus2 = clusters2->GetCluster(index2);
a7477843 1491 if (!clus2) {
1492 AliWarning(Form("Could not find cluster %d!", index2));
1493 continue;
1494 }
1495 Int_t ncells2 = clus2->GetNCells();
9adcb46d 1496 if (ncells2 >= maxNcells2) {
1497 AliError(Form("Number of cells in the cluster %d >= %d",ncells2,maxNcells2));
a7477843 1498 continue;
1499 }
1500 UShort_t *cellsId2 = clus2->GetCellsAbsId();
1501 Double_t *cellsFrac2 = clus2->GetCellsAmplitudeFraction();
1502
1503 TLorentzVector part2;
1504 clus2->GetMomentum(part2, const_cast<Double_t*>(fVertex));
1505 Double_t ptClus2 = part2.Pt();
1506
1507 TMath::Sort(ncells2, cellsId2, sortedIndexes2, kFALSE);
1508
1509 for (Int_t iClus1 = 0; iClus1 < nClus1; iClus1++) {
9adcb46d 1510 if (sortedIndexes1[iClus1] == 0)
1511 continue;
a7477843 1512 Int_t iCell1 = 0, iCell2 = 0;
a7477843 1513 while (iCell1 < ncells1[iClus1] && iCell2 < ncells2) {
1514 if (cellsId1[iClus1][sortedIndexes1[iClus1][iCell1]] == cellsId2[sortedIndexes2[iCell2]]) { // found a common cell
1515 d1 -= cellsFrac1[iClus1][sortedIndexes1[iClus1][iCell1]] * ptClus1[iClus1];
1516 d2 -= cellsFrac2[sortedIndexes2[iCell2]] * ptClus2;
1517 iCell1++;
1518 iCell2++;
a7477843 1519 }
1520 else if (cellsId1[iClus1][sortedIndexes1[iClus1][iCell1]] > cellsId2[sortedIndexes2[iCell2]]) {
1521 iCell2++;
1522 }
1523 else {
1524 iCell1++;
1525 }
cac0e10b 1526 }
a7477843 1527 }
1528 }
1529 }
1530 else {
1531 for (Int_t iClus2 = 0; iClus2 < jet2->GetNumberOfClusters(); iClus2++) {
1532 Int_t index2 = jet2->ClusterAt(iClus2);
7cd832c7 1533 for (Int_t iClus1 = 0; iClus1 < jet1->GetNumberOfClusters(); iClus1++) {
1534 Int_t index1 = jet1->ClusterAt(iClus1);
1535 if (index2 == index1) { // found common particle
1536 AliVCluster *clus1 = clusters1->GetCluster(index1);
1537 if (!clus1) {
1538 AliWarning(Form("Could not find cluster %d!", index1));
a7477843 1539 continue;
1540 }
7cd832c7 1541 AliVCluster *clus2 = clusters2->GetCluster(index2);
a7477843 1542 if (!clus2) {
1543 AliWarning(Form("Could not find cluster %d!", index2));
1544 continue;
1545 }
7cd832c7 1546 TLorentzVector part1, part2;
1547 clus1->GetMomentum(part1, const_cast<Double_t*>(fVertex));
a7477843 1548 clus2->GetMomentum(part2, const_cast<Double_t*>(fVertex));
1549
7cd832c7 1550 d1 -= part1.Pt();
a7477843 1551 d2 -= part2.Pt();
1552 break;
cac0e10b 1553 }
7f76e479 1554 }
7f76e479 1555 }
1556 }
7f76e479 1557 }
1558
05077f28 1559 if (d1 < 0)
1560 d1 = 0;
1561
1562 if (d2 < 0)
1563 d2 = 0;
1564
1565 if (jet1->Pt() > 0)
7f76e479 1566 d1 /= jet1->Pt();
1567 else
05077f28 1568 d1 = -1;
7f76e479 1569
05077f28 1570 if (jet2->Pt() > 0)
7f76e479 1571 d2 /= jet2->Pt();
1572 else
05077f28 1573 d2 = -1;
7f76e479 1574}
1575
1576//________________________________________________________________________
1577void AliJetResponseMaker::SetMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, MatchingType matching)
1578{
1579 Double_t d1 = -1;
1580 Double_t d2 = -1;
1581
1582 switch (matching) {
1583 case kGeometrical:
1584 GetGeometricalMatchingLevel(jet1,jet2,d1);
1585 d2 = d1;
1586 break;
1587 case kMCLabel: // jet1 = detector level and jet2 = particle level!
1588 GetMCLabelMatchingLevel(jet1,jet2,d1,d2);
1589 break;
1590 case kSameCollections:
1591 GetSameCollectionsMatchingLevel(jet1,jet2,d1,d2);
cd6431de 1592 break;
1593 default:
1594 ;
1595 }
f660c2d6 1596
05077f28 1597 if (d1 >= 0) {
f660c2d6 1598
1599 if (d1 < jet1->ClosestJetDistance()) {
1600 jet1->SetSecondClosestJet(jet1->ClosestJet(), jet1->ClosestJetDistance());
1601 jet1->SetClosestJet(jet2, d1);
1602 }
1603 else if (d1 < jet1->SecondClosestJetDistance()) {
1604 jet1->SetSecondClosestJet(jet2, d1);
1605 }
1606 }
1607
05077f28 1608 if (d2 >= 0) {
f660c2d6 1609
1610 if (d2 < jet2->ClosestJetDistance()) {
1611 jet2->SetSecondClosestJet(jet2->ClosestJet(), jet2->ClosestJetDistance());
1612 jet2->SetClosestJet(jet1, d2);
1613 }
1614 else if (d2 < jet2->SecondClosestJetDistance()) {
1615 jet2->SetSecondClosestJet(jet1, d2);
1616 }
1617 }
cd6431de 1618}
1619
1ee1b5b8 1620//________________________________________________________________________
1621Bool_t AliJetResponseMaker::FillHistograms()
1622{
1623 // Fill histograms.
1624
7cd832c7 1625 AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
1626 AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
2949a2ec 1627
7cd832c7 1628 if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return kFALSE;
ca5c29fa 1629
7cd832c7 1630 AliEmcalJet* jet1 = 0;
1631 AliEmcalJet* jet2 = 0;
1632
1633 jets2->ResetCurrentID();
1634 while ((jet2 = jets2->GetNextJet())) {
d3c0e12a 1635
7cd832c7 1636 AliDebug(2,Form("Processing jet (2) %d", jets2->GetCurrentID()));
ca5c29fa 1637
b48964d0 1638 if (jet2->Pt() < jets2->GetJetPtCut()) continue;
1639
d3c0e12a 1640 Double_t ptLeading2 = jets2->GetLeadingHadronPt(jet2);
1641 Double_t corrpt2 = jet2->Pt() - jets2->GetRhoVal() * jet2->Area();
7cd832c7 1642
7cd832c7 1643 if (jets2->AcceptJet(jet2))
d3c0e12a 1644 FillJetHisto(jet2->Phi(), jet2->Eta(), jet2->Pt(), jet2->Area(), jet2->NEF(), ptLeading2,
1645 corrpt2, jet2->MCPt(), 2);
68e03fc3 1646 else
1647 fHistRejectionReason2->Fill(jets2->GetRejectionReasonBitPosition(), jet2->Pt());
7cd832c7 1648
1649 jet1 = jet2->MatchedJet();
1650
1651 if (!jet1) continue;
1652 if (!jets1->AcceptJet(jet1)) continue;
1653 if (jet1->MCPt() < fMinJetMCPt) continue;
1654
1655 Double_t d=-1, ce1=-1, ce2=-1;
1656 if (jet2->GetMatchingType() == kGeometrical) {
1657 if (jets2->GetIsParticleLevel() && !jets1->GetIsParticleLevel()) // the other way around is not supported
1658 GetMCLabelMatchingLevel(jet1, jet2, ce1, ce2);
1659 else if (jets1->GetIsParticleLevel() == jets2->GetIsParticleLevel())
1660 GetSameCollectionsMatchingLevel(jet1, jet2, ce1, ce2);
43032ce2 1661
7cd832c7 1662 d = jet2->ClosestJetDistance();
43032ce2 1663 }
7cd832c7 1664 else if (jet2->GetMatchingType() == kMCLabel || jet2->GetMatchingType() == kSameCollections) {
1665 GetGeometricalMatchingLevel(jet1, jet2, d);
43032ce2 1666
7cd832c7 1667 ce1 = jet1->ClosestJetDistance();
1668 ce2 = jet2->ClosestJetDistance();
2949a2ec 1669 }
1b76c28f 1670
d3c0e12a 1671 Double_t ptLeading1 = jets1->GetLeadingHadronPt(jet1);
7cd832c7 1672 Double_t corrpt1 = jet1->Pt() - jets1->GetRhoVal() * jet1->Area();
1b76c28f 1673
7cd832c7 1674 FillMatchingHistos(jet1->Pt(), jet2->Pt(), jet1->Eta(), jet2->Eta(), jet1->Phi(), jet2->Phi(),
1675 jet1->Area(), jet2->Area(), d, ce1, ce2, corrpt1, corrpt2, jet1->MCPt(),
d3c0e12a 1676 jet1->NEF(), jet2->NEF(), ptLeading1, ptLeading2);
7cd832c7 1677 }
2949a2ec 1678
7cd832c7 1679 jets1->ResetCurrentID();
68e03fc3 1680 while ((jet1 = jets1->GetNextJet())) {
1681 if (!jets1->AcceptJet(jet1)) {
1682 fHistRejectionReason1->Fill(jets1->GetRejectionReasonBitPosition(), jet1->Pt());
1683 continue;
1684 }
7cd832c7 1685 if (jet1->MCPt() < fMinJetMCPt) continue;
1686 AliDebug(2,Form("Processing jet (1) %d", jets1->GetCurrentID()));
ca5c29fa 1687
d3c0e12a 1688 Double_t ptLeading1 = jets1->GetLeadingHadronPt(jet1);
1689 Double_t corrpt1 = jet1->Pt() - jets1->GetRhoVal() * jet1->Area();
1690
1691 FillJetHisto(jet1->Phi(), jet1->Eta(), jet1->Pt(), jet1->Area(), jet1->NEF(), ptLeading1,
1692 corrpt1, jet1->MCPt(), 1);
1b76c28f 1693 }
6fd5039f 1694 return kTRUE;
1b76c28f 1695}