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