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