3 // Emcal jet response matrix maker task.
7 #include "AliJetResponseMaker.h"
9 #include <TClonesArray.h>
11 #include <THnSparse.h>
12 #include <TLorentzVector.h>
14 #include "AliAnalysisManager.h"
15 #include "AliVCluster.h"
16 #include "AliVTrack.h"
17 #include "AliEmcalJet.h"
19 #include "AliRhoParameter.h"
20 #include "AliNamedArrayI.h"
21 #include "AliJetContainer.h"
22 #include "AliParticleContainer.h"
23 #include "AliClusterContainer.h"
25 ClassImp(AliJetResponseMaker)
27 //________________________________________________________________________
28 AliJetResponseMaker::AliJetResponseMaker() :
29 AliAnalysisTaskEmcalJet("AliJetResponseMaker", kTRUE),
30 fMatching(kNoMatching),
33 fUseCellsToMatch(kFALSE),
37 fDeltaEtaDeltaPhiAxis(0),
47 fHistJets1CorrPtArea(0),
49 fHistJets1CEFvsCEFPt(0),
53 fHistJets2CorrPtArea(0),
55 fHistJets2CEFvsCEFPt(0),
57 fHistCommonEnergy1vsJet1Pt(0),
58 fHistCommonEnergy2vsJet2Pt(0),
59 fHistDistancevsJet1Pt(0),
60 fHistDistancevsJet2Pt(0),
61 fHistDistancevsCommonEnergy1(0),
62 fHistDistancevsCommonEnergy2(0),
63 fHistCommonEnergy1vsCommonEnergy2(0),
64 fHistDeltaEtaDeltaPhi(0),
65 fHistJet2PtOverJet1PtvsJet2Pt(0),
66 fHistJet1PtOverJet2PtvsJet1Pt(0),
67 fHistDeltaPtvsJet1Pt(0),
68 fHistDeltaPtvsJet2Pt(0),
69 fHistDeltaPtOverJet1PtvsJet1Pt(0),
70 fHistDeltaPtOverJet2PtvsJet2Pt(0),
71 fHistDeltaPtvsDistance(0),
72 fHistDeltaPtvsCommonEnergy1(0),
73 fHistDeltaPtvsCommonEnergy2(0),
74 fHistDeltaPtvsArea1(0),
75 fHistDeltaPtvsArea2(0),
76 fHistDeltaPtvsDeltaArea(0),
77 fHistJet1PtvsJet2Pt(0),
78 fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt(0),
79 fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt(0),
80 fHistDeltaCorrPtvsJet1CorrPt(0),
81 fHistDeltaCorrPtvsJet2CorrPt(0),
82 fHistDeltaCorrPtvsDistance(0),
83 fHistDeltaCorrPtvsCommonEnergy1(0),
84 fHistDeltaCorrPtvsCommonEnergy2(0),
85 fHistDeltaCorrPtvsArea1(0),
86 fHistDeltaCorrPtvsArea2(0),
87 fHistDeltaCorrPtvsDeltaArea(0),
88 fHistJet1CorrPtvsJet2CorrPt(0),
89 fHistDeltaMCPtOverJet1MCPtvsJet1MCPt(0),
90 fHistDeltaMCPtOverJet2PtvsJet2Pt(0),
91 fHistDeltaMCPtvsJet1MCPt(0),
92 fHistDeltaMCPtvsJet2Pt(0),
93 fHistDeltaMCPtvsDistance(0),
94 fHistDeltaMCPtvsCommonEnergy1(0),
95 fHistDeltaMCPtvsCommonEnergy2(0),
96 fHistDeltaMCPtvsArea1(0),
97 fHistDeltaMCPtvsArea2(0),
98 fHistDeltaMCPtvsDeltaArea(0),
99 fHistJet1MCPtvsJet2Pt(0)
101 // Default constructor.
103 SetMakeGeneralHistograms(kTRUE);
106 //________________________________________________________________________
107 AliJetResponseMaker::AliJetResponseMaker(const char *name) :
108 AliAnalysisTaskEmcalJet(name, kTRUE),
109 fMatching(kNoMatching),
112 fUseCellsToMatch(kFALSE),
116 fDeltaEtaDeltaPhiAxis(0),
126 fHistJets1CorrPtArea(0),
127 fHistJets1NEFvsPt(0),
128 fHistJets1CEFvsCEFPt(0),
132 fHistJets2CorrPtArea(0),
133 fHistJets2NEFvsPt(0),
134 fHistJets2CEFvsCEFPt(0),
136 fHistCommonEnergy1vsJet1Pt(0),
137 fHistCommonEnergy2vsJet2Pt(0),
138 fHistDistancevsJet1Pt(0),
139 fHistDistancevsJet2Pt(0),
140 fHistDistancevsCommonEnergy1(0),
141 fHistDistancevsCommonEnergy2(0),
142 fHistCommonEnergy1vsCommonEnergy2(0),
143 fHistDeltaEtaDeltaPhi(0),
144 fHistJet2PtOverJet1PtvsJet2Pt(0),
145 fHistJet1PtOverJet2PtvsJet1Pt(0),
146 fHistDeltaPtvsJet1Pt(0),
147 fHistDeltaPtvsJet2Pt(0),
148 fHistDeltaPtOverJet1PtvsJet1Pt(0),
149 fHistDeltaPtOverJet2PtvsJet2Pt(0),
150 fHistDeltaPtvsDistance(0),
151 fHistDeltaPtvsCommonEnergy1(0),
152 fHistDeltaPtvsCommonEnergy2(0),
153 fHistDeltaPtvsArea1(0),
154 fHistDeltaPtvsArea2(0),
155 fHistDeltaPtvsDeltaArea(0),
156 fHistJet1PtvsJet2Pt(0),
157 fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt(0),
158 fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt(0),
159 fHistDeltaCorrPtvsJet1CorrPt(0),
160 fHistDeltaCorrPtvsJet2CorrPt(0),
161 fHistDeltaCorrPtvsDistance(0),
162 fHistDeltaCorrPtvsCommonEnergy1(0),
163 fHistDeltaCorrPtvsCommonEnergy2(0),
164 fHistDeltaCorrPtvsArea1(0),
165 fHistDeltaCorrPtvsArea2(0),
166 fHistDeltaCorrPtvsDeltaArea(0),
167 fHistJet1CorrPtvsJet2CorrPt(0),
168 fHistDeltaMCPtOverJet1MCPtvsJet1MCPt(0),
169 fHistDeltaMCPtOverJet2PtvsJet2Pt(0),
170 fHistDeltaMCPtvsJet1MCPt(0),
171 fHistDeltaMCPtvsJet2Pt(0),
172 fHistDeltaMCPtvsDistance(0),
173 fHistDeltaMCPtvsCommonEnergy1(0),
174 fHistDeltaMCPtvsCommonEnergy2(0),
175 fHistDeltaMCPtvsArea1(0),
176 fHistDeltaMCPtvsArea2(0),
177 fHistDeltaMCPtvsDeltaArea(0),
178 fHistJet1MCPtvsJet2Pt(0)
180 // Standard constructor.
182 SetMakeGeneralHistograms(kTRUE);
185 //________________________________________________________________________
186 AliJetResponseMaker::~AliJetResponseMaker()
192 //________________________________________________________________________
193 void AliJetResponseMaker::AllocateTH2()
195 // Allocate TH2 histograms.
197 fHistJets1PhiEta = new TH2F("fHistJets1PhiEta", "fHistJets1PhiEta", 40, -1, 1, 40, 0, TMath::Pi()*2);
198 fHistJets1PhiEta->GetXaxis()->SetTitle("#eta");
199 fHistJets1PhiEta->GetYaxis()->SetTitle("#phi");
200 fOutput->Add(fHistJets1PhiEta);
202 fHistJets1PtArea = new TH2F("fHistJets1PtArea", "fHistJets1PtArea", fNbins/2, 0, 1.5, fNbins, fMinBinPt, fMaxBinPt);
203 fHistJets1PtArea->GetXaxis()->SetTitle("area");
204 fHistJets1PtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
205 fHistJets1PtArea->GetZaxis()->SetTitle("counts");
206 fOutput->Add(fHistJets1PtArea);
209 fHistJets1CorrPtArea = new TH2F("fHistJets1CorrPtArea", "fHistJets1CorrPtArea",
210 fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
211 fHistJets1CorrPtArea->GetXaxis()->SetTitle("area");
212 fHistJets1CorrPtArea->GetYaxis()->SetTitle("p_{T,1}^{corr} (GeV/c)");
213 fHistJets1CorrPtArea->GetZaxis()->SetTitle("counts");
214 fOutput->Add(fHistJets1CorrPtArea);
217 fHistJets1ZvsPt = new TH2F("fHistJets1ZvsPt", "fHistJets1ZvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
218 fHistJets1ZvsPt->GetXaxis()->SetTitle("Z");
219 fHistJets1ZvsPt->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
220 fHistJets1ZvsPt->GetZaxis()->SetTitle("counts");
221 fOutput->Add(fHistJets1ZvsPt);
223 fHistJets1NEFvsPt = new TH2F("fHistJets1NEFvsPt", "fHistJets1NEFvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
224 fHistJets1NEFvsPt->GetXaxis()->SetTitle("NEF");
225 fHistJets1NEFvsPt->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
226 fHistJets1NEFvsPt->GetZaxis()->SetTitle("counts");
227 fOutput->Add(fHistJets1NEFvsPt);
229 fHistJets1CEFvsCEFPt = new TH2F("fHistJets1CEFvsCEFPt", "fHistJets1CEFvsCEFPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
230 fHistJets1CEFvsCEFPt->GetXaxis()->SetTitle("1-NEF");
231 fHistJets1CEFvsCEFPt->GetYaxis()->SetTitle("(1-NEF)*p_{T,1} (GeV/c)");
232 fHistJets1CEFvsCEFPt->GetZaxis()->SetTitle("counts");
233 fOutput->Add(fHistJets1CEFvsCEFPt);
237 fHistJets2PhiEta = new TH2F("fHistJets2PhiEta", "fHistJets2PhiEta", 40, -1, 1, 40, 0, TMath::Pi()*2);
238 fHistJets2PhiEta->GetXaxis()->SetTitle("#eta");
239 fHistJets2PhiEta->GetYaxis()->SetTitle("#phi");
240 fOutput->Add(fHistJets2PhiEta);
242 fHistJets2PtArea = new TH2F("fHistJets2PtArea", "fHistJets2PtArea", fNbins/2, 0, 1.5, fNbins, fMinBinPt, fMaxBinPt);
243 fHistJets2PtArea->GetXaxis()->SetTitle("area");
244 fHistJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
245 fHistJets2PtArea->GetZaxis()->SetTitle("counts");
246 fOutput->Add(fHistJets2PtArea);
249 fHistJets2CorrPtArea = new TH2F("fHistJets2CorrPtArea", "fHistJets2CorrPtArea",
250 fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
251 fHistJets2CorrPtArea->GetXaxis()->SetTitle("area");
252 fHistJets2CorrPtArea->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)");
253 fHistJets2CorrPtArea->GetZaxis()->SetTitle("counts");
254 fOutput->Add(fHistJets2CorrPtArea);
257 fHistJets2ZvsPt = new TH2F("fHistJets2ZvsPt", "fHistJets2ZvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
258 fHistJets2ZvsPt->GetXaxis()->SetTitle("Z");
259 fHistJets2ZvsPt->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
260 fHistJets2ZvsPt->GetZaxis()->SetTitle("counts");
261 fOutput->Add(fHistJets2ZvsPt);
263 fHistJets2NEFvsPt = new TH2F("fHistJets2NEFvsPt", "fHistJets2NEFvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
264 fHistJets2NEFvsPt->GetXaxis()->SetTitle("NEF");
265 fHistJets2NEFvsPt->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
266 fHistJets2NEFvsPt->GetZaxis()->SetTitle("counts");
267 fOutput->Add(fHistJets2NEFvsPt);
269 fHistJets2CEFvsCEFPt = new TH2F("fHistJets2CEFvsCEFPt", "fHistJets2CEFvsCEFPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
270 fHistJets2CEFvsCEFPt->GetXaxis()->SetTitle("1-NEF");
271 fHistJets2CEFvsCEFPt->GetYaxis()->SetTitle("(1-NEF)*p_{T,2} (GeV/c)");
272 fHistJets2CEFvsCEFPt->GetZaxis()->SetTitle("counts");
273 fOutput->Add(fHistJets2CEFvsCEFPt);
275 // Matching histograms
277 fHistCommonEnergy1vsJet1Pt = new TH2F("fHistCommonEnergy1vsJet1Pt", "fHistCommonEnergy1vsJet1Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt);
278 fHistCommonEnergy1vsJet1Pt->GetXaxis()->SetTitle("Common energy 1 (%)");
279 fHistCommonEnergy1vsJet1Pt->GetYaxis()->SetTitle("p_{T,1}");
280 fHistCommonEnergy1vsJet1Pt->GetZaxis()->SetTitle("counts");
281 fOutput->Add(fHistCommonEnergy1vsJet1Pt);
283 fHistCommonEnergy2vsJet2Pt = new TH2F("fHistCommonEnergy2vsJet2Pt", "fHistCommonEnergy2vsJet2Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt);
284 fHistCommonEnergy2vsJet2Pt->GetXaxis()->SetTitle("Common energy 2 (%)");
285 fHistCommonEnergy2vsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");
286 fHistCommonEnergy2vsJet2Pt->GetZaxis()->SetTitle("counts");
287 fOutput->Add(fHistCommonEnergy2vsJet2Pt);
289 fHistDistancevsJet1Pt = new TH2F("fHistDistancevsJet1Pt", "fHistDistancevsJet1Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt);
290 fHistDistancevsJet1Pt->GetXaxis()->SetTitle("Distance");
291 fHistDistancevsJet1Pt->GetYaxis()->SetTitle("p_{T,1}");
292 fHistDistancevsJet1Pt->GetZaxis()->SetTitle("counts");
293 fOutput->Add(fHistDistancevsJet1Pt);
295 fHistDistancevsJet2Pt = new TH2F("fHistDistancevsJet2Pt", "fHistDistancevsJet2Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt);
296 fHistDistancevsJet2Pt->GetXaxis()->SetTitle("Distance");
297 fHistDistancevsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");
298 fHistDistancevsJet2Pt->GetZaxis()->SetTitle("counts");
299 fOutput->Add(fHistDistancevsJet2Pt);
301 fHistDistancevsCommonEnergy1 = new TH2F("fHistDistancevsCommonEnergy1", "fHistDistancevsCommonEnergy1", fNbins/2, 0, 1.2, fNbins/2, 0, 1.2);
302 fHistDistancevsCommonEnergy1->GetXaxis()->SetTitle("Distance");
303 fHistDistancevsCommonEnergy1->GetYaxis()->SetTitle("Common energy 1 (%)");
304 fHistDistancevsCommonEnergy1->GetZaxis()->SetTitle("counts");
305 fOutput->Add(fHistDistancevsCommonEnergy1);
307 fHistDistancevsCommonEnergy2 = new TH2F("fHistDistancevsCommonEnergy2", "fHistDistancevsCommonEnergy2", fNbins/2, 0, 1.2, fNbins/2, 0, 1.2);
308 fHistDistancevsCommonEnergy2->GetXaxis()->SetTitle("Distance");
309 fHistDistancevsCommonEnergy2->GetYaxis()->SetTitle("Common energy 2 (%)");
310 fHistDistancevsCommonEnergy2->GetZaxis()->SetTitle("counts");
311 fOutput->Add(fHistDistancevsCommonEnergy2);
313 fHistCommonEnergy1vsCommonEnergy2 = new TH2F("fHistCommonEnergy1vsCommonEnergy2", "fHistCommonEnergy1vsCommonEnergy2", fNbins/2, 0, 1.2, fNbins/2, 0, 1.2);
314 fHistCommonEnergy1vsCommonEnergy2->GetXaxis()->SetTitle("Common energy 1 (%)");
315 fHistCommonEnergy1vsCommonEnergy2->GetYaxis()->SetTitle("Common energy 2 (%)");
316 fHistCommonEnergy1vsCommonEnergy2->GetZaxis()->SetTitle("counts");
317 fOutput->Add(fHistCommonEnergy1vsCommonEnergy2);
319 fHistDeltaEtaDeltaPhi = new TH2F("fHistDeltaEtaDeltaPhi", "fHistDeltaEtaDeltaPhi", fNbins/4, -1, 1, fNbins/4, -TMath::Pi()/2, TMath::Pi()*3/2);
320 fHistDeltaEtaDeltaPhi->GetXaxis()->SetTitle("Common energy 1 (%)");
321 fHistDeltaEtaDeltaPhi->GetYaxis()->SetTitle("Common energy 2 (%)");
322 fHistDeltaEtaDeltaPhi->GetZaxis()->SetTitle("counts");
323 fOutput->Add(fHistDeltaEtaDeltaPhi);
325 fHistJet2PtOverJet1PtvsJet2Pt = new TH2F("fHistJet2PtOverJet1PtvsJet2Pt", "fHistJet2PtOverJet1PtvsJet2Pt", fNbins, fMinBinPt, fMaxBinPt, 300, 0, 1.5);
326 fHistJet2PtOverJet1PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");
327 fHistJet2PtOverJet1PtvsJet2Pt->GetYaxis()->SetTitle("p_{T,2} / p_{T,1}");
328 fHistJet2PtOverJet1PtvsJet2Pt->GetZaxis()->SetTitle("counts");
329 fOutput->Add(fHistJet2PtOverJet1PtvsJet2Pt);
331 fHistJet1PtOverJet2PtvsJet1Pt = new TH2F("fHistJet1PtOverJet2PtvsJet1Pt", "fHistJet1PtOverJet2PtvsJet1Pt", fNbins, fMinBinPt, fMaxBinPt, 300, 0, 1.5);
332 fHistJet1PtOverJet2PtvsJet1Pt->GetXaxis()->SetTitle("p_{T,1}");
333 fHistJet1PtOverJet2PtvsJet1Pt->GetYaxis()->SetTitle("p_{T,1} / p_{T,2}");
334 fHistJet1PtOverJet2PtvsJet1Pt->GetZaxis()->SetTitle("counts");
335 fOutput->Add(fHistJet1PtOverJet2PtvsJet1Pt);
337 fHistDeltaPtvsJet1Pt = new TH2F("fHistDeltaPtvsJet1Pt", "fHistDeltaPtvsJet1Pt",
338 fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
339 fHistDeltaPtvsJet1Pt->GetXaxis()->SetTitle("p_{T,1}");
340 fHistDeltaPtvsJet1Pt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
341 fHistDeltaPtvsJet1Pt->GetZaxis()->SetTitle("counts");
342 fOutput->Add(fHistDeltaPtvsJet1Pt);
344 fHistDeltaPtvsJet2Pt = new TH2F("fHistDeltaPtvsJet2Pt", "fHistDeltaPtvsJet2Pt",
345 fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
346 fHistDeltaPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");
347 fHistDeltaPtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
348 fHistDeltaPtvsJet2Pt->GetZaxis()->SetTitle("counts");
349 fOutput->Add(fHistDeltaPtvsJet2Pt);
351 fHistDeltaPtOverJet1PtvsJet1Pt = new TH2F("fHistDeltaPtOverJet1PtvsJet1Pt", "fHistDeltaPtOverJet1PtvsJet1Pt",
352 fNbins, fMinBinPt, fMaxBinPt, fNbins, -5, 5);
353 fHistDeltaPtOverJet1PtvsJet1Pt->GetXaxis()->SetTitle("p_{T,1}");
354 fHistDeltaPtOverJet1PtvsJet1Pt->GetYaxis()->SetTitle("#deltap_{T} / p_{T,1}");
355 fHistDeltaPtOverJet1PtvsJet1Pt->GetZaxis()->SetTitle("counts");
356 fOutput->Add(fHistDeltaPtOverJet1PtvsJet1Pt);
358 fHistDeltaPtOverJet2PtvsJet2Pt = new TH2F("fHistDeltaPtOverJet2PtvsJet2Pt", "fHistDeltaPtOverJet2PtvsJet2Pt",
359 fNbins, fMinBinPt, fMaxBinPt, fNbins, -5, 5);
360 fHistDeltaPtOverJet2PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");
361 fHistDeltaPtOverJet2PtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} / p_{T,2}");
362 fHistDeltaPtOverJet2PtvsJet2Pt->GetZaxis()->SetTitle("counts");
363 fOutput->Add(fHistDeltaPtOverJet2PtvsJet2Pt);
365 fHistDeltaPtvsDistance = new TH2F("fHistDeltaPtvsDistance", "fHistDeltaPtvsDistance",
366 fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
367 fHistDeltaPtvsDistance->GetXaxis()->SetTitle("Distance");
368 fHistDeltaPtvsDistance->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
369 fHistDeltaPtvsDistance->GetZaxis()->SetTitle("counts");
370 fOutput->Add(fHistDeltaPtvsDistance);
372 fHistDeltaPtvsCommonEnergy1 = new TH2F("fHistDeltaPtvsCommonEnergy1", "fHistDeltaPtvsCommonEnergy1",
373 fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
374 fHistDeltaPtvsCommonEnergy1->GetXaxis()->SetTitle("Common energy 1 (%)");
375 fHistDeltaPtvsCommonEnergy1->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
376 fHistDeltaPtvsCommonEnergy1->GetZaxis()->SetTitle("counts");
377 fOutput->Add(fHistDeltaPtvsCommonEnergy1);
379 fHistDeltaPtvsCommonEnergy2 = new TH2F("fHistDeltaPtvsCommonEnergy2", "fHistDeltaPtvsCommonEnergy2",
380 fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
381 fHistDeltaPtvsCommonEnergy2->GetXaxis()->SetTitle("Common energy 2 (%)");
382 fHistDeltaPtvsCommonEnergy2->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
383 fHistDeltaPtvsCommonEnergy2->GetZaxis()->SetTitle("counts");
384 fOutput->Add(fHistDeltaPtvsCommonEnergy2);
386 fHistDeltaPtvsArea1 = new TH2F("fHistDeltaPtvsArea1", "fHistDeltaPtvsArea1",
387 fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
388 fHistDeltaPtvsArea1->GetXaxis()->SetTitle("A_{jet,1}");
389 fHistDeltaPtvsArea1->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
390 fHistDeltaPtvsArea1->GetZaxis()->SetTitle("counts");
391 fOutput->Add(fHistDeltaPtvsArea1);
393 fHistDeltaPtvsArea2 = new TH2F("fHistDeltaPtvsArea2", "fHistDeltaPtvsArea2",
394 fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
395 fHistDeltaPtvsArea2->GetXaxis()->SetTitle("A_{jet,2}");
396 fHistDeltaPtvsArea2->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
397 fHistDeltaPtvsArea2->GetZaxis()->SetTitle("counts");
398 fOutput->Add(fHistDeltaPtvsArea2);
400 fHistDeltaPtvsDeltaArea = new TH2F("fHistDeltaPtvsDeltaArea", "fHistDeltaPtvsDeltaArea",
401 fNbins, -1.+1./fNbins, 1.+1./fNbins, 2*fNbins, -fMaxBinPt, fMaxBinPt);
402 fHistDeltaPtvsDeltaArea->GetXaxis()->SetTitle("#deltaA_{jet}");
403 fHistDeltaPtvsDeltaArea->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
404 fHistDeltaPtvsDeltaArea->GetZaxis()->SetTitle("counts");
405 fOutput->Add(fHistDeltaPtvsDeltaArea);
407 fHistJet1PtvsJet2Pt = new TH2F("fHistJet1PtvsJet2Pt", "fHistJet1PtvsJet2Pt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
408 fHistJet1PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,1}");
409 fHistJet1PtvsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");
410 fHistJet1PtvsJet2Pt->GetZaxis()->SetTitle("counts");
411 fOutput->Add(fHistJet1PtvsJet2Pt);
413 if (fIsJet1Rho || fIsJet2Rho) {
415 fHistDeltaCorrPtvsJet1CorrPt = new TH2F("fHistDeltaCorrPtvsJet1CorrPt", "fHistDeltaCorrPtvsJet1CorrPt",
416 fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
418 fHistDeltaCorrPtvsJet1CorrPt = new TH2F("fHistDeltaCorrPtvsJet1CorrPt", "fHistDeltaCorrPtvsJet1CorrPt",
419 2*fNbins, -fMaxBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
421 fHistDeltaCorrPtvsJet1CorrPt->GetXaxis()->SetTitle("p_{T,1}^{corr}");
422 fHistDeltaCorrPtvsJet1CorrPt->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
423 fHistDeltaCorrPtvsJet1CorrPt->GetZaxis()->SetTitle("counts");
424 fOutput->Add(fHistDeltaCorrPtvsJet1CorrPt);
427 fHistDeltaCorrPtvsJet2CorrPt = new TH2F("fHistDeltaCorrPtvsJet2CorrPt", "fHistDeltaCorrPtvsJet2CorrPt",
428 fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
430 fHistDeltaCorrPtvsJet2CorrPt = new TH2F("fHistDeltaCorrPtvsJet2CorrPt", "fHistDeltaCorrPtvsJet2CorrPt",
431 2*fNbins, -fMaxBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
433 fHistDeltaCorrPtvsJet2CorrPt->GetXaxis()->SetTitle("p_{T,2}^{corr}");
434 fHistDeltaCorrPtvsJet2CorrPt->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
435 fHistDeltaCorrPtvsJet2CorrPt->GetZaxis()->SetTitle("counts");
436 fOutput->Add(fHistDeltaCorrPtvsJet2CorrPt);
438 fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt = new TH2F("fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt", "fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt",
439 2*fNbins, -fMaxBinPt, fMaxBinPt, fNbins, -5, 5);
440 fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt->GetXaxis()->SetTitle("p_{T,1}^{corr}");
441 fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt->GetYaxis()->SetTitle("#deltap_{T}^{corr} / p_{T,1}^{corr}");
442 fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt->GetZaxis()->SetTitle("counts");
443 fOutput->Add(fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt);
445 fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt = new TH2F("fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt", "fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt",
446 2*fNbins, -fMaxBinPt, fMaxBinPt, fNbins, -5, 5);
447 fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt->GetXaxis()->SetTitle("p_{T,2}^{corr}");
448 fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt->GetYaxis()->SetTitle("#deltap_{T}^{corr} / p_{T,2}^{corr}");
449 fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt->GetZaxis()->SetTitle("counts");
450 fOutput->Add(fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt);
452 fHistDeltaCorrPtvsDistance = new TH2F("fHistDeltaCorrPtvsDistance", "fHistDeltaCorrPtvsDistance",
453 fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
454 fHistDeltaCorrPtvsDistance->GetXaxis()->SetTitle("Distance");
455 fHistDeltaCorrPtvsDistance->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
456 fHistDeltaCorrPtvsDistance->GetZaxis()->SetTitle("counts");
457 fOutput->Add(fHistDeltaCorrPtvsDistance);
459 fHistDeltaCorrPtvsCommonEnergy1 = new TH2F("fHistDeltaCorrPtvsCommonEnergy1", "fHistDeltaCorrPtvsCommonEnergy1",
460 fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
461 fHistDeltaCorrPtvsCommonEnergy1->GetXaxis()->SetTitle("Common energy 1 (%)");
462 fHistDeltaCorrPtvsCommonEnergy1->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
463 fHistDeltaCorrPtvsCommonEnergy1->GetZaxis()->SetTitle("counts");
464 fOutput->Add(fHistDeltaCorrPtvsCommonEnergy1);
466 fHistDeltaCorrPtvsCommonEnergy2 = new TH2F("fHistDeltaCorrPtvsCommonEnergy2", "fHistDeltaCorrPtvsCommonEnergy2",
467 fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
468 fHistDeltaCorrPtvsCommonEnergy2->GetXaxis()->SetTitle("Common energy 2 (%)");
469 fHistDeltaCorrPtvsCommonEnergy2->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
470 fHistDeltaCorrPtvsCommonEnergy2->GetZaxis()->SetTitle("counts");
471 fOutput->Add(fHistDeltaCorrPtvsCommonEnergy2);
473 fHistDeltaCorrPtvsArea1 = new TH2F("fHistDeltaCorrPtvsArea1", "fHistDeltaCorrPtvsArea1",
474 fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
475 fHistDeltaCorrPtvsArea1->GetXaxis()->SetTitle("A_{jet,1}");
476 fHistDeltaCorrPtvsArea1->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
477 fHistDeltaCorrPtvsArea1->GetZaxis()->SetTitle("counts");
478 fOutput->Add(fHistDeltaCorrPtvsArea1);
480 fHistDeltaCorrPtvsArea2 = new TH2F("fHistDeltaCorrPtvsArea2", "fHistDeltaCorrPtvsArea2",
481 fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
482 fHistDeltaCorrPtvsArea2->GetXaxis()->SetTitle("A_{jet,2}");
483 fHistDeltaCorrPtvsArea2->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
484 fHistDeltaCorrPtvsArea2->GetZaxis()->SetTitle("counts");
485 fOutput->Add(fHistDeltaCorrPtvsArea2);
487 fHistDeltaCorrPtvsDeltaArea = new TH2F("fHistDeltaCorrPtvsDeltaArea", "fHistDeltaCorrPtvsDeltaArea",
488 fNbins, -1.+1./fNbins, 1.+1./fNbins, 2*fNbins, -fMaxBinPt, fMaxBinPt);
489 fHistDeltaCorrPtvsDeltaArea->GetXaxis()->SetTitle("#deltaA_{jet}");
490 fHistDeltaCorrPtvsDeltaArea->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
491 fHistDeltaCorrPtvsDeltaArea->GetZaxis()->SetTitle("counts");
492 fOutput->Add(fHistDeltaCorrPtvsDeltaArea);
495 fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt",
496 fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
497 else if (!fIsJet2Rho)
498 fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt",
499 2*fNbins, -fMaxBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
501 fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt",
502 2*fNbins, -fMaxBinPt, fMaxBinPt,
503 2*fNbins, -fMaxBinPt, fMaxBinPt);
504 fHistJet1CorrPtvsJet2CorrPt->GetXaxis()->SetTitle("p_{T,1}^{corr}");
505 fHistJet1CorrPtvsJet2CorrPt->GetYaxis()->SetTitle("p_{T,2}^{corr}");
506 fHistJet1CorrPtvsJet2CorrPt->GetZaxis()->SetTitle("counts");
507 fOutput->Add(fHistJet1CorrPtvsJet2CorrPt);
511 fHistDeltaMCPtvsJet1MCPt = new TH2F("fHistDeltaMCPtvsJet1MCPt", "fHistDeltaMCPtvsJet1MCPt",
512 fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
513 fHistDeltaMCPtvsJet1MCPt->GetXaxis()->SetTitle("p_{T,1}^{MC}");
514 fHistDeltaMCPtvsJet1MCPt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
515 fHistDeltaMCPtvsJet1MCPt->GetZaxis()->SetTitle("counts");
516 fOutput->Add(fHistDeltaMCPtvsJet1MCPt);
518 fHistDeltaMCPtvsJet2Pt = new TH2F("fHistDeltaMCPtvsJet2Pt", "fHistDeltaMCPtvsJet2Pt",
519 fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
520 fHistDeltaMCPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");
521 fHistDeltaMCPtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
522 fHistDeltaMCPtvsJet2Pt->GetZaxis()->SetTitle("counts");
523 fOutput->Add(fHistDeltaMCPtvsJet2Pt);
525 fHistDeltaMCPtOverJet1MCPtvsJet1MCPt = new TH2F("fHistDeltaMCPtOverJet1MCPtvsJet1MCPt", "fHistDeltaMCPtOverJet1MCPtvsJet1MCPt",
526 fNbins, fMinBinPt, fMaxBinPt, fNbins, -5, 5);
527 fHistDeltaMCPtOverJet1MCPtvsJet1MCPt->GetXaxis()->SetTitle("p_{T,1}^{MC}");
528 fHistDeltaMCPtOverJet1MCPtvsJet1MCPt->GetYaxis()->SetTitle("#deltap_{T} / p_{T,1}^{MC}");
529 fHistDeltaMCPtOverJet1MCPtvsJet1MCPt->GetZaxis()->SetTitle("counts");
530 fOutput->Add(fHistDeltaMCPtOverJet1MCPtvsJet1MCPt);
532 fHistDeltaMCPtOverJet2PtvsJet2Pt = new TH2F("fHistDeltaMCPtOverJet2PtvsJet2Pt", "fHistDeltaMCPtOverJet2PtvsJet2Pt",
533 fNbins, fMinBinPt, fMaxBinPt, fNbins, -5, 5);
534 fHistDeltaMCPtOverJet2PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");
535 fHistDeltaMCPtOverJet2PtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} / p_{T,2}");
536 fHistDeltaMCPtOverJet2PtvsJet2Pt->GetZaxis()->SetTitle("counts");
537 fOutput->Add(fHistDeltaMCPtOverJet2PtvsJet2Pt);
539 fHistDeltaMCPtvsDistance = new TH2F("fHistDeltaMCPtvsDistance", "fHistDeltaMCPtvsDistance",
540 fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
541 fHistDeltaMCPtvsDistance->GetXaxis()->SetTitle("Distance");
542 fHistDeltaMCPtvsDistance->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
543 fHistDeltaMCPtvsDistance->GetZaxis()->SetTitle("counts");
544 fOutput->Add(fHistDeltaMCPtvsDistance);
546 fHistDeltaMCPtvsCommonEnergy1 = new TH2F("fHistDeltaMCPtvsCommonEnergy1", "fHistDeltaMCPtvsCommonEnergy1",
547 fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
548 fHistDeltaMCPtvsCommonEnergy1->GetXaxis()->SetTitle("Common energy 1 (%)");
549 fHistDeltaMCPtvsCommonEnergy1->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
550 fHistDeltaMCPtvsCommonEnergy1->GetZaxis()->SetTitle("counts");
551 fOutput->Add(fHistDeltaMCPtvsCommonEnergy1);
553 fHistDeltaMCPtvsCommonEnergy2 = new TH2F("fHistDeltaMCPtvsCommonEnergy2", "fHistDeltaMCPtvsCommonEnergy2",
554 fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
555 fHistDeltaMCPtvsCommonEnergy2->GetXaxis()->SetTitle("Common energy 2 (%)");
556 fHistDeltaMCPtvsCommonEnergy2->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
557 fHistDeltaMCPtvsCommonEnergy2->GetZaxis()->SetTitle("counts");
558 fOutput->Add(fHistDeltaMCPtvsCommonEnergy2);
560 fHistDeltaMCPtvsArea1 = new TH2F("fHistDeltaMCPtvsArea1", "fHistDeltaMCPtvsArea1",
561 fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
562 fHistDeltaMCPtvsArea1->GetXaxis()->SetTitle("A_{jet,1}");
563 fHistDeltaMCPtvsArea1->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
564 fHistDeltaMCPtvsArea1->GetZaxis()->SetTitle("counts");
565 fOutput->Add(fHistDeltaMCPtvsArea1);
567 fHistDeltaMCPtvsArea2 = new TH2F("fHistDeltaMCPtvsArea2", "fHistDeltaMCPtvsArea2",
568 fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
569 fHistDeltaMCPtvsArea2->GetXaxis()->SetTitle("A_{jet,2}");
570 fHistDeltaMCPtvsArea2->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
571 fHistDeltaMCPtvsArea2->GetZaxis()->SetTitle("counts");
572 fOutput->Add(fHistDeltaMCPtvsArea2);
574 fHistDeltaMCPtvsDeltaArea = new TH2F("fHistDeltaMCPtvsDeltaArea", "fHistDeltaMCPtvsDeltaArea",
575 fNbins, -1.+1./fNbins, 1.+1./fNbins, 2*fNbins, -fMaxBinPt, fMaxBinPt);
576 fHistDeltaMCPtvsDeltaArea->GetXaxis()->SetTitle("#deltaA_{jet}");
577 fHistDeltaMCPtvsDeltaArea->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
578 fHistDeltaMCPtvsDeltaArea->GetZaxis()->SetTitle("counts");
579 fOutput->Add(fHistDeltaMCPtvsDeltaArea);
581 fHistJet1MCPtvsJet2Pt = new TH2F("fHistJet1MCPtvsJet2Pt", "fHistJet1MCPtvsJet2Pt",
582 fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
583 fHistJet1MCPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,1}^{MC}");
584 fHistJet1MCPtvsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");
585 fHistJet1MCPtvsJet2Pt->GetZaxis()->SetTitle("counts");
586 fOutput->Add(fHistJet1MCPtvsJet2Pt);
590 //________________________________________________________________________
591 void AliJetResponseMaker::AllocateTHnSparse()
593 // Allocate THnSparse histograms.
595 TString title[25]= {""};
596 Int_t nbins[25] = {0};
597 Double_t min[25] = {0.};
598 Double_t max[25] = {0.};
602 nbins[dim] = fNbins/4;
604 max[dim] = 2*TMath::Pi()*(1 + 1./(nbins[dim]-1));
608 nbins[dim] = fNbins/4;
613 title[dim] = "p_{T}";
619 title[dim] = "A_{jet}";
620 nbins[dim] = fNbins/4;
627 nbins[dim] = fNbins/4;
635 nbins[dim] = fNbins/4;
641 title[dim] = "p_{T,particle}^{leading} (GeV/c)";
647 Int_t dim1 = dim, dim2 = dim;
650 title[dim1] = "p_{T}^{corr}";
651 nbins[dim1] = fNbins*2;
658 title[dim1] = "p_{T}^{MC}";
659 nbins[dim1] = fNbins;
665 fHistJets1 = new THnSparseD("fHistJets1","fHistJets1",dim1,nbins,min,max);
666 for (Int_t i = 0; i < dim1; i++)
667 fHistJets1->GetAxis(i)->SetTitle(title[i]);
668 fOutput->Add(fHistJets1);
671 title[dim2] = "p_{T}^{corr}";
672 nbins[dim2] = fNbins*2;
678 fHistJets2 = new THnSparseD("fHistJets2","fHistJets2",dim2,nbins,min,max);
679 for (Int_t i = 0; i < dim2; i++)
680 fHistJets2->GetAxis(i)->SetTitle(title[i]);
681 fOutput->Add(fHistJets2);
687 title[dim] = "p_{T,1}";
693 title[dim] = "p_{T,2}";
699 title[dim] = "A_{jet,1}";
700 nbins[dim] = fNbins/4;
705 title[dim] = "A_{jet,2}";
706 nbins[dim] = fNbins/4;
711 title[dim] = "distance";
712 nbins[dim] = fNbins/2;
718 nbins[dim] = fNbins/2;
724 nbins[dim] = fNbins/2;
729 title[dim] = "p_{T,particle,1}^{leading} (GeV/c)";
735 title[dim] = "p_{T,particle,2}^{leading} (GeV/c)";
742 title[dim] = "#deltaA_{jet}";
743 nbins[dim] = fNbins/2;
748 title[dim] = "#deltap_{T}";
749 nbins[dim] = fNbins*2;
755 title[dim] = "p_{T,1}^{corr}";
756 nbins[dim] = fNbins*2;
762 title[dim] = "p_{T,2}^{corr}";
763 nbins[dim] = fNbins*2;
768 if (fDeltaPtAxis && (fIsJet1Rho || fIsJet2Rho)) {
769 title[dim] = "#deltap_{T}^{corr}";
770 nbins[dim] = fNbins*2;
775 if (fDeltaEtaDeltaPhiAxis) {
776 title[dim] = "#delta#eta";
777 nbins[dim] = fNbins/2;
782 title[dim] = "#delta#phi";
783 nbins[dim] = fNbins/2;
784 min[dim] = -TMath::Pi()/2;
785 max[dim] = TMath::Pi()*3/2;
789 title[dim] = "p_{T,1}^{MC}";
796 title[dim] = "#deltap_{T}^{MC}";
797 nbins[dim] = fNbins*2;
805 title[dim] = "NEF_{1}";
806 nbins[dim] = fNbins/4;
811 title[dim] = "NEF_{2}";
812 nbins[dim] = fNbins/4;
819 title[dim] = "Z_{1}";
820 nbins[dim] = fNbins/4;
825 title[dim] = "Z_{2}";
826 nbins[dim] = fNbins/4;
832 fHistMatching = new THnSparseD("fHistMatching","fHistMatching",dim,nbins,min,max);
834 for (Int_t i = 0; i < dim; i++)
835 fHistMatching->GetAxis(i)->SetTitle(title[i]);
837 fOutput->Add(fHistMatching);
840 //________________________________________________________________________
841 void AliJetResponseMaker::UserCreateOutputObjects()
843 // Create user objects.
845 AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
847 AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
848 AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
850 if (!jets1 || !jets2) return;
852 if (jets1->GetRhoName().IsNull()) fIsJet1Rho = kFALSE;
853 else fIsJet1Rho = kTRUE;
854 if (jets2->GetRhoName().IsNull()) fIsJet2Rho = kFALSE;
855 else fIsJet2Rho = kTRUE;
862 PostData(1, fOutput); // Post data for ALL output slots > 0 here, to get at least an empty histogram
865 //________________________________________________________________________
866 void AliJetResponseMaker::FillJetHisto(Double_t Phi, Double_t Eta, Double_t Pt, Double_t A, Double_t NEF, Double_t LeadingPt, Double_t CorrPt, Double_t MCPt, Int_t Set)
869 THnSparse *histo = 0;
878 Double_t contents[20]={0};
880 for (Int_t i = 0; i < histo->GetNdimensions(); i++) {
881 TString title(histo->GetAxis(i)->GetTitle());
884 else if (title=="#eta")
886 else if (title=="p_{T}")
888 else if (title=="A_{jet}")
890 else if (title=="NEF")
893 contents[i] = LeadingPt/Pt;
894 else if (title=="p_{T}^{corr}")
895 contents[i] = CorrPt;
896 else if (title=="p_{T}^{MC}")
898 else if (title=="p_{T,particle}^{leading} (GeV/c)")
899 contents[i] = LeadingPt;
901 AliWarning(Form("Unable to fill dimension %s!",title.Data()));
904 histo->Fill(contents);
908 fHistJets1PtArea->Fill(A, Pt);
909 fHistJets1PhiEta->Fill(Eta, Phi);
910 fHistJets1ZvsPt->Fill(LeadingPt/Pt, Pt);
911 fHistJets1NEFvsPt->Fill(NEF, Pt);
912 fHistJets1CEFvsCEFPt->Fill(1-NEF, (1-NEF)*Pt);
914 fHistJets1CorrPtArea->Fill(A, CorrPt);
917 fHistJets2PtArea->Fill(A, Pt);
918 fHistJets2PhiEta->Fill(Eta, Phi);
919 fHistJets2ZvsPt->Fill(LeadingPt/Pt, Pt);
920 fHistJets2NEFvsPt->Fill(NEF, Pt);
921 fHistJets2CEFvsCEFPt->Fill(1-NEF, (1-NEF)*Pt);
923 fHistJets2CorrPtArea->Fill(A, CorrPt);
928 //________________________________________________________________________
929 void AliJetResponseMaker::FillMatchingHistos(Double_t Pt1, Double_t Pt2, Double_t Eta1, Double_t Eta2, Double_t Phi1, Double_t Phi2,
930 Double_t A1, Double_t A2, Double_t d, Double_t CE1, Double_t CE2, Double_t CorrPt1, Double_t CorrPt2,
931 Double_t MCPt1, Double_t NEF1, Double_t NEF2, Double_t LeadingPt1, Double_t LeadingPt2)
934 Double_t contents[20]={0};
936 for (Int_t i = 0; i < fHistMatching->GetNdimensions(); i++) {
937 TString title(fHistMatching->GetAxis(i)->GetTitle());
938 if (title=="p_{T,1}")
940 else if (title=="p_{T,2}")
942 else if (title=="A_{jet,1}")
944 else if (title=="A_{jet,2}")
946 else if (title=="distance")
948 else if (title=="CE1")
950 else if (title=="CE2")
952 else if (title=="#deltaA_{jet}")
954 else if (title=="#deltap_{T}")
955 contents[i] = Pt1-Pt2;
956 else if (title=="#delta#eta")
957 contents[i] = Eta1-Eta2;
958 else if (title=="#delta#phi")
959 contents[i] = Phi1-Phi2;
960 else if (title=="p_{T,1}^{corr}")
961 contents[i] = CorrPt1;
962 else if (title=="p_{T,2}^{corr}")
963 contents[i] = CorrPt2;
964 else if (title=="#deltap_{T}^{corr}")
965 contents[i] = CorrPt1-CorrPt2;
966 else if (title=="p_{T,1}^{MC}")
968 else if (title=="#deltap_{T}^{MC}")
969 contents[i] = MCPt1-Pt2;
970 else if (title=="NEF_{1}")
972 else if (title=="NEF_{2}")
974 else if (title=="Z_{1}")
975 contents[i] = LeadingPt1/Pt1;
976 else if (title=="Z_{2}")
977 contents[i] = LeadingPt2/Pt2;
978 else if (title=="p_{T,particle,1}^{leading} (GeV/c)")
979 contents[i] = LeadingPt1;
980 else if (title=="p_{T,particle,2}^{leading} (GeV/c)")
981 contents[i] = LeadingPt2;
983 AliWarning(Form("Unable to fill dimension %s!",title.Data()));
986 fHistMatching->Fill(contents);
989 fHistCommonEnergy1vsJet1Pt->Fill(CE1, Pt1);
990 fHistCommonEnergy2vsJet2Pt->Fill(CE2, Pt2);
992 fHistDistancevsJet1Pt->Fill(d, Pt1);
993 fHistDistancevsJet2Pt->Fill(d, Pt2);
995 fHistDistancevsCommonEnergy1->Fill(d, CE1);
996 fHistDistancevsCommonEnergy2->Fill(d, CE2);
997 fHistCommonEnergy1vsCommonEnergy2->Fill(CE1, CE2);
999 fHistDeltaEtaDeltaPhi->Fill(Eta1-Eta2,Phi1-Phi2);
1001 fHistJet2PtOverJet1PtvsJet2Pt->Fill(Pt2, Pt2 / Pt1);
1002 fHistJet1PtOverJet2PtvsJet1Pt->Fill(Pt1, Pt1 / Pt2);
1004 Double_t dpt = Pt1 - Pt2;
1005 fHistDeltaPtvsJet1Pt->Fill(Pt1, dpt);
1006 fHistDeltaPtvsJet2Pt->Fill(Pt2, dpt);
1007 fHistDeltaPtOverJet1PtvsJet1Pt->Fill(Pt1, dpt/Pt1);
1008 fHistDeltaPtOverJet2PtvsJet2Pt->Fill(Pt2, dpt/Pt2);
1010 fHistDeltaPtvsDistance->Fill(d, dpt);
1011 fHistDeltaPtvsCommonEnergy1->Fill(CE1, dpt);
1012 fHistDeltaPtvsCommonEnergy2->Fill(CE2, dpt);
1014 fHistDeltaPtvsArea1->Fill(A1, dpt);
1015 fHistDeltaPtvsArea2->Fill(A2, dpt);
1017 Double_t darea = A1 - A2;
1018 fHistDeltaPtvsDeltaArea->Fill(darea, dpt);
1020 fHistJet1PtvsJet2Pt->Fill(Pt1, Pt2);
1022 if (fIsJet1Rho || fIsJet2Rho) {
1023 Double_t dcorrpt = CorrPt1 - CorrPt2;
1024 fHistDeltaCorrPtvsJet1CorrPt->Fill(CorrPt1, dcorrpt);
1025 fHistDeltaCorrPtvsJet2CorrPt->Fill(CorrPt2, dcorrpt);
1026 fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt->Fill(CorrPt1, dcorrpt/CorrPt1);
1027 fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt->Fill(CorrPt2, dcorrpt/CorrPt2);
1028 fHistDeltaCorrPtvsDistance->Fill(d, dcorrpt);
1029 fHistDeltaCorrPtvsCommonEnergy1->Fill(CE1, dcorrpt);
1030 fHistDeltaCorrPtvsCommonEnergy2->Fill(CE2, dcorrpt);
1031 fHistDeltaCorrPtvsArea1->Fill(A1, dcorrpt);
1032 fHistDeltaCorrPtvsArea2->Fill(A2, dcorrpt);
1033 fHistDeltaCorrPtvsDeltaArea->Fill(darea, dcorrpt);
1034 fHistJet1CorrPtvsJet2CorrPt->Fill(CorrPt1, CorrPt2);
1038 Double_t dmcpt = MCPt1 - Pt2;
1039 fHistDeltaMCPtvsJet1MCPt->Fill(MCPt1, dmcpt);
1040 fHistDeltaMCPtvsJet2Pt->Fill(Pt2, dmcpt);
1041 fHistDeltaMCPtOverJet1MCPtvsJet1MCPt->Fill(MCPt1, dmcpt/MCPt1);
1042 fHistDeltaMCPtOverJet2PtvsJet2Pt->Fill(Pt2, dmcpt/Pt2);
1043 fHistDeltaMCPtvsDistance->Fill(d, dmcpt);
1044 fHistDeltaMCPtvsCommonEnergy1->Fill(CE1, dmcpt);
1045 fHistDeltaMCPtvsCommonEnergy2->Fill(CE2, dmcpt);
1046 fHistDeltaMCPtvsArea1->Fill(A1, dmcpt);
1047 fHistDeltaMCPtvsArea2->Fill(A2, dmcpt);
1048 fHistDeltaMCPtvsDeltaArea->Fill(darea, dmcpt);
1049 fHistJet1MCPtvsJet2Pt->Fill(MCPt1, Pt2);
1054 //________________________________________________________________________
1055 void AliJetResponseMaker::ExecOnce()
1059 AliAnalysisTaskEmcalJet::ExecOnce();
1061 AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
1062 AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
1064 if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return;
1066 if (fMatching == kMCLabel && (!jets2->GetIsParticleLevel() || jets1->GetIsParticleLevel())) {
1067 if (jets1->GetIsParticleLevel() == jets2->GetIsParticleLevel()) {
1068 AliWarning("Changing matching type from MC label to same collection...");
1069 fMatching = kSameCollections;
1072 AliWarning("Changing matching type from MC label to geometrical...");
1073 fMatching = kGeometrical;
1076 else if (fMatching == kSameCollections && jets1->GetIsParticleLevel() != jets2->GetIsParticleLevel()) {
1077 if (jets2->GetIsParticleLevel() && !jets1->GetIsParticleLevel()) {
1078 AliWarning("Changing matching type from same collection to MC label...");
1079 fMatching = kMCLabel;
1082 AliWarning("Changing matching type from same collection to geometrical...");
1083 fMatching = kGeometrical;
1088 //________________________________________________________________________
1089 Bool_t AliJetResponseMaker::Run()
1091 // Find the closest jets
1093 if (fMatching == kNoMatching)
1096 return DoJetMatching();
1099 //________________________________________________________________________
1100 Bool_t AliJetResponseMaker::DoJetMatching()
1102 AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
1103 AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
1105 if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return kFALSE;
1109 AliEmcalJet* jet1 = 0;
1111 jets1->ResetCurrentID();
1112 while ((jet1 = jets1->GetNextJet())) {
1114 AliEmcalJet *jet2 = jet1->ClosestJet();
1116 if (!jet2) continue;
1117 if (jet2->ClosestJet() != jet1) continue;
1118 if (jet1->ClosestJetDistance() > fMatchingPar1 || jet2->ClosestJetDistance() > fMatchingPar2) continue;
1120 // Matched jet found
1121 jet1->SetMatchedToClosest(fMatching);
1122 jet2->SetMatchedToClosest(fMatching);
1123 AliDebug(2,Form("Found matching: jet1 pt = %f, eta = %f, phi = %f, jet2 pt = %f, eta = %f, phi = %f",
1124 jet1->Pt(), jet1->Eta(), jet1->Phi(),
1125 jet2->Pt(), jet2->Eta(), jet2->Phi()));
1131 //________________________________________________________________________
1132 void AliJetResponseMaker::DoJetLoop()
1136 AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
1137 AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
1139 if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return;
1141 AliEmcalJet* jet1 = 0;
1142 AliEmcalJet* jet2 = 0;
1144 jets2->ResetCurrentID();
1145 while ((jet2 = jets2->GetNextJet())) jet2->ResetMatching();
1147 jets1->ResetCurrentID();
1148 while ((jet1 = jets1->GetNextJet())) {
1149 jet1->ResetMatching();
1151 if (jet1->MCPt() < fMinJetMCPt) continue;
1153 jets2->ResetCurrentID();
1154 while ((jet2 = jets2->GetNextJet())) {
1155 SetMatchingLevel(jet1, jet2, fMatching);
1160 //________________________________________________________________________
1161 void AliJetResponseMaker::GetGeometricalMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d) const
1163 Double_t deta = jet2->Eta() - jet1->Eta();
1164 Double_t dphi = jet2->Phi() - jet1->Phi();
1165 d = TMath::Sqrt(deta * deta + dphi * dphi);
1168 //________________________________________________________________________
1169 void AliJetResponseMaker::GetMCLabelMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d1, Double_t &d2) const
1171 AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
1172 AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
1174 if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return;
1176 AliParticleContainer *tracks1 = jets1->GetParticleContainer();
1177 AliClusterContainer *clusters1 = jets1->GetClusterContainer();
1178 AliParticleContainer *tracks2 = jets2->GetParticleContainer();
1180 // d1 and d2 represent the matching level: 0 = maximum level of matching, 1 = the two jets are completely unrelated
1183 Double_t totalPt1 = d1; // the total pt of the reconstructed jet will be cleaned from the background
1185 // remove completely tracks that are not MC particles (label == 0)
1186 if (tracks1 && tracks1->GetArray()) {
1187 for (Int_t iTrack = 0; iTrack < jet1->GetNumberOfTracks(); iTrack++) {
1188 AliVParticle *track = jet1->TrackAt(iTrack,tracks1->GetArray());
1190 AliWarning(Form("Could not find track %d!", iTrack));
1194 Int_t MClabel = TMath::Abs(track->GetLabel());
1195 MClabel -= fMCLabelShift;
1196 if (MClabel != 0) continue;
1198 // this is not a MC particle; remove it completely
1199 AliDebug(3,Form("Track %d (pT = %f) is not a MC particle (MClabel = %d)!",iTrack,track->Pt(),MClabel));
1200 totalPt1 -= track->Pt();
1205 // remove completely clusters that are not MC particles (label == 0)
1206 if (fUseCellsToMatch && fCaloCells) {
1207 for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
1208 AliVCluster *clus = jet1->ClusterAt(iClus,clusters1->GetArray());
1210 AliWarning(Form("Could not find cluster %d!", iClus));
1213 TLorentzVector part;
1214 clus->GetMomentum(part, const_cast<Double_t*>(fVertex));
1216 for (Int_t iCell = 0; iCell < clus->GetNCells(); iCell++) {
1217 Int_t cellId = clus->GetCellAbsId(iCell);
1218 Double_t cellFrac = clus->GetCellAmplitudeFraction(iCell);
1220 Int_t MClabel = TMath::Abs(fCaloCells->GetCellMCLabel(cellId));
1221 MClabel -= fMCLabelShift;
1222 if (MClabel != 0) continue;
1224 // this is not a MC particle; remove it completely
1225 AliDebug(3,Form("Cell %d (frac = %f) is not a MC particle (MClabel = %d)!",iCell,cellFrac,MClabel));
1226 totalPt1 -= part.Pt() * cellFrac;
1227 d1 -= part.Pt() * cellFrac;
1232 for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
1233 AliVCluster *clus = jet1->ClusterAt(iClus,clusters1->GetArray());
1235 AliWarning(Form("Could not find cluster %d!", iClus));
1238 TLorentzVector part;
1239 clus->GetMomentum(part, const_cast<Double_t*>(fVertex));
1241 Int_t MClabel = TMath::Abs(clus->GetLabel());
1242 MClabel -= fMCLabelShift;
1243 if (MClabel != 0) continue;
1245 // this is not a MC particle; remove it completely
1246 AliDebug(3,Form("Cluster %d (pT = %f) is not a MC particle (MClabel = %d)!",iClus,part.Pt(),MClabel));
1247 totalPt1 -= part.Pt();
1252 for (Int_t iTrack2 = 0; iTrack2 < jet2->GetNumberOfTracks(); iTrack2++) {
1253 Bool_t track2Found = kFALSE;
1254 Int_t index2 = jet2->TrackAt(iTrack2);
1256 // now look for common particles in the track array
1257 for (Int_t iTrack = 0; iTrack < jet1->GetNumberOfTracks(); iTrack++) {
1258 AliVParticle *track = jet1->TrackAt(iTrack,tracks1->GetArray());
1260 AliWarning(Form("Could not find track %d!", iTrack));
1263 Int_t MClabel = TMath::Abs(track->GetLabel());
1264 MClabel -= fMCLabelShift;
1265 if (MClabel <= 0) continue;
1268 index = tracks2->GetIndexFromLabel(MClabel);
1270 AliDebug(2,Form("Track %d (pT = %f) does not have an associated MC particle (MClabel = %d)!",iTrack,track->Pt(),MClabel));
1274 if (index2 != index) continue;
1276 // found common particle
1280 AliVParticle *MCpart = tracks2->GetParticle(index2);
1281 AliDebug(3,Form("Track %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!",
1282 iTrack,track->Pt(),track->Eta(),track->Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
1286 track2Found = kTRUE;
1289 // now look for common particles in the cluster array
1290 if (fUseCellsToMatch && fCaloCells) { // if the cell colection is available, look for cells with a matched MC particle
1291 for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
1292 AliVCluster *clus = jet1->ClusterAt(iClus,fCaloClusters);
1294 AliWarning(Form("Could not find cluster %d!", iClus));
1297 TLorentzVector part;
1298 clus->GetMomentum(part, const_cast<Double_t*>(fVertex));
1300 for (Int_t iCell = 0; iCell < clus->GetNCells(); iCell++) {
1301 Int_t cellId = clus->GetCellAbsId(iCell);
1302 Double_t cellFrac = clus->GetCellAmplitudeFraction(iCell);
1304 Int_t MClabel = TMath::Abs(fCaloCells->GetCellMCLabel(cellId));
1305 MClabel -= fMCLabelShift;
1306 if (MClabel <= 0) continue;
1309 index1 = tracks2->GetIndexFromLabel(MClabel);
1311 AliDebug(3,Form("Cell %d (frac = %f) does not have an associated MC particle (MClabel = %d)!",iCell,cellFrac,MClabel));
1315 if (index2 != index1) continue;
1317 // found common particle
1318 d1 -= part.Pt() * cellFrac;
1320 if (!track2Found) { // only if it is not already found among charged tracks (charged particles are most likely already found)
1321 AliVParticle *MCpart = tracks2->GetParticle(index2);
1322 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)!",
1323 iCell,iClus,part.Pt(),part.Eta(),part.Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
1324 d2 -= MCpart->Pt() * cellFrac;
1327 track2Found = kTRUE;
1331 else { //otherwise look for the first contributor to the cluster, and if matched to a MC label remove it
1332 for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
1333 AliVCluster *clus = jet1->ClusterAt(iClus,fCaloClusters);
1335 AliWarning(Form("Could not find cluster %d!", iClus));
1338 TLorentzVector part;
1339 clus->GetMomentum(part, const_cast<Double_t*>(fVertex));
1341 Int_t MClabel = TMath::Abs(clus->GetLabel());
1342 MClabel -= fMCLabelShift;
1343 if (MClabel <= 0) continue;
1346 index = tracks2->GetIndexFromLabel(MClabel);
1349 AliDebug(3,Form("Cluster %d (pT = %f) does not have an associated MC particle (MClabel = %d)!",iClus,part.Pt(),MClabel));
1353 if (index2 != index) continue;
1355 // found common particle
1358 if (!track2Found) { // only if it is not already found among charged tracks (charged particles are most likely already found)
1359 AliVParticle *MCpart = tracks2->GetParticle(index2);
1360 AliDebug(3,Form("Cluster %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!",
1361 iClus,part.Pt(),part.Eta(),part.Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
1366 track2Found = kTRUE;
1388 //________________________________________________________________________
1389 void AliJetResponseMaker::GetSameCollectionsMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d1, Double_t &d2) const
1391 AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
1392 AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
1394 if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return;
1396 AliParticleContainer *tracks1 = jets1->GetParticleContainer();
1397 AliClusterContainer *clusters1 = jets1->GetClusterContainer();
1398 AliParticleContainer *tracks2 = jets2->GetParticleContainer();
1399 AliClusterContainer *clusters2 = jets2->GetClusterContainer();
1401 // d1 and d2 represent the matching level: 0 = maximum level of matching, 1 = the two jets are completely unrelated
1405 if (tracks1 && tracks2) {
1407 for (Int_t iTrack2 = 0; iTrack2 < jet2->GetNumberOfTracks(); iTrack2++) {
1408 Int_t index2 = jet2->TrackAt(iTrack2);
1409 for (Int_t iTrack1 = 0; iTrack1 < jet1->GetNumberOfTracks(); iTrack1++) {
1410 Int_t index1 = jet1->TrackAt(iTrack1);
1411 if (index2 == index1) { // found common particle
1412 AliVParticle *part1 = tracks1->GetParticle(index1);
1414 AliWarning(Form("Could not find track %d!", index1));
1417 AliVParticle *part2 = tracks2->GetParticle(index2);
1419 AliWarning(Form("Could not find track %d!", index2));
1432 if (clusters1 && clusters2) {
1434 if (fUseCellsToMatch) {
1435 const Int_t nClus1 = jet1->GetNumberOfClusters();
1437 Int_t ncells1[nClus1];
1438 UShort_t *cellsId1[nClus1];
1439 Double_t *cellsFrac1[nClus1];
1440 Int_t *sortedIndexes1[nClus1];
1441 Double_t ptClus1[nClus1];
1442 for (Int_t iClus1 = 0; iClus1 < nClus1; iClus1++) {
1443 Int_t index1 = jet1->ClusterAt(iClus1);
1444 AliVCluster *clus1 = clusters1->GetCluster(index1);
1446 AliWarning(Form("Could not find cluster %d!", index1));
1447 ncells1[iClus1] = 0;
1448 cellsId1[iClus1] = 0;
1449 cellsFrac1[iClus1] = 0;
1450 sortedIndexes1[iClus1] = 0;
1451 ptClus1[iClus1] = 0;
1454 TLorentzVector part1;
1455 clus1->GetMomentum(part1, const_cast<Double_t*>(fVertex));
1457 ncells1[iClus1] = clus1->GetNCells();
1458 cellsId1[iClus1] = clus1->GetCellsAbsId();
1459 cellsFrac1[iClus1] = clus1->GetCellsAmplitudeFraction();
1460 sortedIndexes1[iClus1] = new Int_t[ncells1[iClus1]];
1461 ptClus1[iClus1] = part1.Pt();
1463 TMath::Sort(ncells1[iClus1], cellsId1[iClus1], sortedIndexes1[iClus1], kFALSE);
1466 const Int_t nClus2 = jet2->GetNumberOfClusters();
1468 const Int_t maxNcells2 = 11520;
1469 Int_t sortedIndexes2[maxNcells2];
1470 for (Int_t iClus2 = 0; iClus2 < nClus2; iClus2++) {
1471 Int_t index2 = jet2->ClusterAt(iClus2);
1472 AliVCluster *clus2 = clusters2->GetCluster(index2);
1474 AliWarning(Form("Could not find cluster %d!", index2));
1477 Int_t ncells2 = clus2->GetNCells();
1478 if (ncells2 >= maxNcells2) {
1479 AliError(Form("Number of cells in the cluster %d >= %d",ncells2,maxNcells2));
1482 UShort_t *cellsId2 = clus2->GetCellsAbsId();
1483 Double_t *cellsFrac2 = clus2->GetCellsAmplitudeFraction();
1485 TLorentzVector part2;
1486 clus2->GetMomentum(part2, const_cast<Double_t*>(fVertex));
1487 Double_t ptClus2 = part2.Pt();
1489 TMath::Sort(ncells2, cellsId2, sortedIndexes2, kFALSE);
1491 for (Int_t iClus1 = 0; iClus1 < nClus1; iClus1++) {
1492 if (sortedIndexes1[iClus1] == 0)
1494 Int_t iCell1 = 0, iCell2 = 0;
1495 while (iCell1 < ncells1[iClus1] && iCell2 < ncells2) {
1496 if (cellsId1[iClus1][sortedIndexes1[iClus1][iCell1]] == cellsId2[sortedIndexes2[iCell2]]) { // found a common cell
1497 d1 -= cellsFrac1[iClus1][sortedIndexes1[iClus1][iCell1]] * ptClus1[iClus1];
1498 d2 -= cellsFrac2[sortedIndexes2[iCell2]] * ptClus2;
1502 else if (cellsId1[iClus1][sortedIndexes1[iClus1][iCell1]] > cellsId2[sortedIndexes2[iCell2]]) {
1513 for (Int_t iClus2 = 0; iClus2 < jet2->GetNumberOfClusters(); iClus2++) {
1514 Int_t index2 = jet2->ClusterAt(iClus2);
1515 for (Int_t iClus1 = 0; iClus1 < jet1->GetNumberOfClusters(); iClus1++) {
1516 Int_t index1 = jet1->ClusterAt(iClus1);
1517 if (index2 == index1) { // found common particle
1518 AliVCluster *clus1 = clusters1->GetCluster(index1);
1520 AliWarning(Form("Could not find cluster %d!", index1));
1523 AliVCluster *clus2 = clusters2->GetCluster(index2);
1525 AliWarning(Form("Could not find cluster %d!", index2));
1528 TLorentzVector part1, part2;
1529 clus1->GetMomentum(part1, const_cast<Double_t*>(fVertex));
1530 clus2->GetMomentum(part2, const_cast<Double_t*>(fVertex));
1558 //________________________________________________________________________
1559 void AliJetResponseMaker::SetMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, MatchingType matching)
1566 GetGeometricalMatchingLevel(jet1,jet2,d1);
1569 case kMCLabel: // jet1 = detector level and jet2 = particle level!
1570 GetMCLabelMatchingLevel(jet1,jet2,d1,d2);
1572 case kSameCollections:
1573 GetSameCollectionsMatchingLevel(jet1,jet2,d1,d2);
1581 if (d1 < jet1->ClosestJetDistance()) {
1582 jet1->SetSecondClosestJet(jet1->ClosestJet(), jet1->ClosestJetDistance());
1583 jet1->SetClosestJet(jet2, d1);
1585 else if (d1 < jet1->SecondClosestJetDistance()) {
1586 jet1->SetSecondClosestJet(jet2, d1);
1592 if (d2 < jet2->ClosestJetDistance()) {
1593 jet2->SetSecondClosestJet(jet2->ClosestJet(), jet2->ClosestJetDistance());
1594 jet2->SetClosestJet(jet1, d2);
1596 else if (d2 < jet2->SecondClosestJetDistance()) {
1597 jet2->SetSecondClosestJet(jet1, d2);
1602 //________________________________________________________________________
1603 Bool_t AliJetResponseMaker::FillHistograms()
1607 AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
1608 AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
1610 if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return kFALSE;
1612 AliEmcalJet* jet1 = 0;
1613 AliEmcalJet* jet2 = 0;
1615 jets2->ResetCurrentID();
1616 while ((jet2 = jets2->GetNextJet())) {
1618 AliDebug(2,Form("Processing jet (2) %d", jets2->GetCurrentID()));
1620 Double_t ptLeading2 = jets2->GetLeadingHadronPt(jet2);
1621 Double_t corrpt2 = jet2->Pt() - jets2->GetRhoVal() * jet2->Area();
1623 if (jets2->AcceptJet(jet2))
1624 FillJetHisto(jet2->Phi(), jet2->Eta(), jet2->Pt(), jet2->Area(), jet2->NEF(), ptLeading2,
1625 corrpt2, jet2->MCPt(), 2);
1627 jet1 = jet2->MatchedJet();
1629 if (!jet1) continue;
1630 if (!jets1->AcceptJet(jet1)) continue;
1631 if (jet1->MCPt() < fMinJetMCPt) continue;
1633 Double_t d=-1, ce1=-1, ce2=-1;
1634 if (jet2->GetMatchingType() == kGeometrical) {
1635 if (jets2->GetIsParticleLevel() && !jets1->GetIsParticleLevel()) // the other way around is not supported
1636 GetMCLabelMatchingLevel(jet1, jet2, ce1, ce2);
1637 else if (jets1->GetIsParticleLevel() == jets2->GetIsParticleLevel())
1638 GetSameCollectionsMatchingLevel(jet1, jet2, ce1, ce2);
1640 d = jet2->ClosestJetDistance();
1642 else if (jet2->GetMatchingType() == kMCLabel || jet2->GetMatchingType() == kSameCollections) {
1643 GetGeometricalMatchingLevel(jet1, jet2, d);
1645 ce1 = jet1->ClosestJetDistance();
1646 ce2 = jet2->ClosestJetDistance();
1649 Double_t ptLeading1 = jets1->GetLeadingHadronPt(jet1);
1650 Double_t corrpt1 = jet1->Pt() - jets1->GetRhoVal() * jet1->Area();
1652 FillMatchingHistos(jet1->Pt(), jet2->Pt(), jet1->Eta(), jet2->Eta(), jet1->Phi(), jet2->Phi(),
1653 jet1->Area(), jet2->Area(), d, ce1, ce2, corrpt1, corrpt2, jet1->MCPt(),
1654 jet1->NEF(), jet2->NEF(), ptLeading1, ptLeading2);
1657 jets1->ResetCurrentID();
1658 while ((jet1 = jets1->GetNextAcceptJet())) {
1659 if (jet1->MCPt() < fMinJetMCPt) continue;
1660 AliDebug(2,Form("Processing jet (1) %d", jets1->GetCurrentID()));
1662 Double_t ptLeading1 = jets1->GetLeadingHadronPt(jet1);
1663 Double_t corrpt1 = jet1->Pt() - jets1->GetRhoVal() * jet1->Area();
1665 FillJetHisto(jet1->Phi(), jet1->Eta(), jet1->Pt(), jet1->Area(), jet1->NEF(), ptLeading1,
1666 corrpt1, jet1->MCPt(), 1);