AliPhysics  cdeda5a (cdeda5a)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliJetResponseMaker.cxx
Go to the documentation of this file.
1 /*************************************************************************
2  * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 //
16 // Emcal jet response matrix maker task.
17 //
18 // Author: Salvatore Aiola, Yale University, salvatore.aiola@cern.ch
19 
20 #include "AliJetResponseMaker.h"
21 
22 #include <TClonesArray.h>
23 #include <TH2F.h>
24 #include <THnSparse.h>
25 
26 #include "AliTLorentzVector.h"
27 #include "AliAnalysisManager.h"
28 #include "AliVCluster.h"
29 #include "AliVTrack.h"
30 #include "AliEmcalJet.h"
31 #include "AliLog.h"
32 #include "AliRhoParameter.h"
33 #include "AliNamedArrayI.h"
34 #include "AliJetContainer.h"
35 #include "AliParticleContainer.h"
36 #include "AliClusterContainer.h"
37 
39 
40 //________________________________________________________________________
43  fMatching(kNoMatching),
44  fMatchingPar1(0),
45  fMatchingPar2(0),
46  fUseCellsToMatch(kFALSE),
47  fMinJetMCPt(1),
48  fHistoType(0),
49  fDeltaPtAxis(0),
50  fDeltaEtaDeltaPhiAxis(0),
51  fNEFAxis(0),
52  fZAxis(0),
53  fZgAxis(0),
54  fdRAxis(0),
55  fPtgAxis(0),
56  fFlavourZAxis(0),
57  fFlavourPtAxis(0),
58  fJetRelativeEPAngle(0),
59  fIsJet1Rho(kFALSE),
60  fIsJet2Rho(kFALSE),
61  fHistRejectionReason1(0),
62  fHistRejectionReason2(0),
63  fHistJets1(0),
64  fHistJets2(0),
65  fHistMatching(0),
66  fHistJets1PhiEta(0),
67  fHistJets1PtArea(0),
68  fHistJets1CorrPtArea(0),
69  fHistJets1NEFvsPt(0),
70  fHistJets1CEFvsCEFPt(0),
71  fHistJets1ZvsPt(0),
72  fHistJets2PhiEta(0),
73  fHistJets2PtArea(0),
74  fHistJets2CorrPtArea(0),
75  fHistJets2NEFvsPt(0),
76  fHistJets2CEFvsCEFPt(0),
77  fHistJets2ZvsPt(0),
78  fHistCommonEnergy1vsJet1Pt(0),
79  fHistCommonEnergy2vsJet2Pt(0),
80  fHistDistancevsJet1Pt(0),
81  fHistDistancevsJet2Pt(0),
82  fHistDistancevsCommonEnergy1(0),
83  fHistDistancevsCommonEnergy2(0),
84  fHistCommonEnergy1vsCommonEnergy2(0),
85  fHistDeltaEtaDeltaPhi(0),
86  fHistJet2PtOverJet1PtvsJet2Pt(0),
87  fHistJet1PtOverJet2PtvsJet1Pt(0),
88  fHistDeltaPtvsJet1Pt(0),
89  fHistDeltaPtvsJet2Pt(0),
90  fHistDeltaPtOverJet1PtvsJet1Pt(0),
91  fHistDeltaPtOverJet2PtvsJet2Pt(0),
92  fHistDeltaPtvsDistance(0),
93  fHistDeltaPtvsCommonEnergy1(0),
94  fHistDeltaPtvsCommonEnergy2(0),
95  fHistDeltaPtvsArea1(0),
96  fHistDeltaPtvsArea2(0),
97  fHistDeltaPtvsDeltaArea(0),
98  fHistJet1PtvsJet2Pt(0),
99  fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt(0),
100  fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt(0),
101  fHistDeltaCorrPtvsJet1CorrPt(0),
102  fHistDeltaCorrPtvsJet2CorrPt(0),
103  fHistDeltaCorrPtvsDistance(0),
104  fHistDeltaCorrPtvsCommonEnergy1(0),
105  fHistDeltaCorrPtvsCommonEnergy2(0),
106  fHistDeltaCorrPtvsArea1(0),
107  fHistDeltaCorrPtvsArea2(0),
108  fHistDeltaCorrPtvsDeltaArea(0),
109  fHistJet1CorrPtvsJet2CorrPt(0),
110  fHistDeltaMCPtOverJet1MCPtvsJet1MCPt(0),
111  fHistDeltaMCPtOverJet2PtvsJet2Pt(0),
112  fHistDeltaMCPtvsJet1MCPt(0),
113  fHistDeltaMCPtvsJet2Pt(0),
114  fHistDeltaMCPtvsDistance(0),
115  fHistDeltaMCPtvsCommonEnergy1(0),
116  fHistDeltaMCPtvsCommonEnergy2(0),
117  fHistDeltaMCPtvsArea1(0),
118  fHistDeltaMCPtvsArea2(0),
119  fHistDeltaMCPtvsDeltaArea(0),
120  fHistJet1MCPtvsJet2Pt(0)
121 {
122  // Default constructor.
123 
124  SetMakeGeneralHistograms(kTRUE);
125 }
126 
127 //________________________________________________________________________
129  AliAnalysisTaskEmcalJet(name, kTRUE),
130  fMatching(kNoMatching),
131  fMatchingPar1(0),
132  fMatchingPar2(0),
133  fUseCellsToMatch(kFALSE),
134  fMinJetMCPt(1),
135  fHistoType(0),
136  fDeltaPtAxis(0),
137  fDeltaEtaDeltaPhiAxis(0),
138  fNEFAxis(0),
139  fZAxis(0),
140  fZgAxis(0),
141  fdRAxis(0),
142  fPtgAxis(0),
143  fFlavourZAxis(0),
144  fFlavourPtAxis(0),
145  fJetRelativeEPAngle(0),
146  fIsJet1Rho(kFALSE),
147  fIsJet2Rho(kFALSE),
148  fHistRejectionReason1(0),
149  fHistRejectionReason2(0),
150  fHistJets1(0),
151  fHistJets2(0),
152  fHistMatching(0),
153  fHistJets1PhiEta(0),
154  fHistJets1PtArea(0),
155  fHistJets1CorrPtArea(0),
156  fHistJets1NEFvsPt(0),
157  fHistJets1CEFvsCEFPt(0),
158  fHistJets1ZvsPt(0),
159  fHistJets2PhiEta(0),
160  fHistJets2PtArea(0),
161  fHistJets2CorrPtArea(0),
162  fHistJets2NEFvsPt(0),
163  fHistJets2CEFvsCEFPt(0),
164  fHistJets2ZvsPt(0),
165  fHistCommonEnergy1vsJet1Pt(0),
166  fHistCommonEnergy2vsJet2Pt(0),
167  fHistDistancevsJet1Pt(0),
168  fHistDistancevsJet2Pt(0),
169  fHistDistancevsCommonEnergy1(0),
170  fHistDistancevsCommonEnergy2(0),
171  fHistCommonEnergy1vsCommonEnergy2(0),
172  fHistDeltaEtaDeltaPhi(0),
173  fHistJet2PtOverJet1PtvsJet2Pt(0),
174  fHistJet1PtOverJet2PtvsJet1Pt(0),
175  fHistDeltaPtvsJet1Pt(0),
176  fHistDeltaPtvsJet2Pt(0),
177  fHistDeltaPtOverJet1PtvsJet1Pt(0),
178  fHistDeltaPtOverJet2PtvsJet2Pt(0),
179  fHistDeltaPtvsDistance(0),
180  fHistDeltaPtvsCommonEnergy1(0),
181  fHistDeltaPtvsCommonEnergy2(0),
182  fHistDeltaPtvsArea1(0),
183  fHistDeltaPtvsArea2(0),
184  fHistDeltaPtvsDeltaArea(0),
185  fHistJet1PtvsJet2Pt(0),
186  fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt(0),
187  fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt(0),
188  fHistDeltaCorrPtvsJet1CorrPt(0),
189  fHistDeltaCorrPtvsJet2CorrPt(0),
190  fHistDeltaCorrPtvsDistance(0),
191  fHistDeltaCorrPtvsCommonEnergy1(0),
192  fHistDeltaCorrPtvsCommonEnergy2(0),
193  fHistDeltaCorrPtvsArea1(0),
194  fHistDeltaCorrPtvsArea2(0),
195  fHistDeltaCorrPtvsDeltaArea(0),
196  fHistJet1CorrPtvsJet2CorrPt(0),
197  fHistDeltaMCPtOverJet1MCPtvsJet1MCPt(0),
198  fHistDeltaMCPtOverJet2PtvsJet2Pt(0),
199  fHistDeltaMCPtvsJet1MCPt(0),
200  fHistDeltaMCPtvsJet2Pt(0),
201  fHistDeltaMCPtvsDistance(0),
202  fHistDeltaMCPtvsCommonEnergy1(0),
203  fHistDeltaMCPtvsCommonEnergy2(0),
204  fHistDeltaMCPtvsArea1(0),
205  fHistDeltaMCPtvsArea2(0),
206  fHistDeltaMCPtvsDeltaArea(0),
207  fHistJet1MCPtvsJet2Pt(0)
208 {
209  // Standard constructor.
210 
212 }
213 
214 //________________________________________________________________________
216 {
217  // Destructor
218 }
219 
220 
221 //________________________________________________________________________
223 {
224  // Allocate TH2 histograms.
225 
226  fHistJets1PhiEta = new TH2F("fHistJets1PhiEta", "fHistJets1PhiEta", 40, -1, 1, 40, 0, TMath::Pi()*2);
227  fHistJets1PhiEta->GetXaxis()->SetTitle("#eta");
228  fHistJets1PhiEta->GetYaxis()->SetTitle("#phi");
230 
231  fHistJets1PtArea = new TH2F("fHistJets1PtArea", "fHistJets1PtArea", fNbins/2, 0, 1.5, fNbins, fMinBinPt, fMaxBinPt);
232  fHistJets1PtArea->GetXaxis()->SetTitle("area");
233  fHistJets1PtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
234  fHistJets1PtArea->GetZaxis()->SetTitle("counts");
236 
237  if (fIsJet1Rho) {
238  fHistJets1CorrPtArea = new TH2F("fHistJets1CorrPtArea", "fHistJets1CorrPtArea",
239  fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
240  fHistJets1CorrPtArea->GetXaxis()->SetTitle("area");
241  fHistJets1CorrPtArea->GetYaxis()->SetTitle("p_{T,1}^{corr} (GeV/c)");
242  fHistJets1CorrPtArea->GetZaxis()->SetTitle("counts");
244  }
245 
246  fHistJets1ZvsPt = new TH2F("fHistJets1ZvsPt", "fHistJets1ZvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
247  fHistJets1ZvsPt->GetXaxis()->SetTitle("Z");
248  fHistJets1ZvsPt->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
249  fHistJets1ZvsPt->GetZaxis()->SetTitle("counts");
250  fOutput->Add(fHistJets1ZvsPt);
251 
252  fHistJets1NEFvsPt = new TH2F("fHistJets1NEFvsPt", "fHistJets1NEFvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
253  fHistJets1NEFvsPt->GetXaxis()->SetTitle("NEF");
254  fHistJets1NEFvsPt->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
255  fHistJets1NEFvsPt->GetZaxis()->SetTitle("counts");
257 
258  fHistJets1CEFvsCEFPt = new TH2F("fHistJets1CEFvsCEFPt", "fHistJets1CEFvsCEFPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
259  fHistJets1CEFvsCEFPt->GetXaxis()->SetTitle("1-NEF");
260  fHistJets1CEFvsCEFPt->GetYaxis()->SetTitle("(1-NEF)*p_{T,1} (GeV/c)");
261  fHistJets1CEFvsCEFPt->GetZaxis()->SetTitle("counts");
263 
264  // Jets 2 histograms
265 
266  fHistJets2PhiEta = new TH2F("fHistJets2PhiEta", "fHistJets2PhiEta", 40, -1, 1, 40, 0, TMath::Pi()*2);
267  fHistJets2PhiEta->GetXaxis()->SetTitle("#eta");
268  fHistJets2PhiEta->GetYaxis()->SetTitle("#phi");
270 
271  fHistJets2PtArea = new TH2F("fHistJets2PtArea", "fHistJets2PtArea", fNbins/2, 0, 1.5, fNbins, fMinBinPt, fMaxBinPt);
272  fHistJets2PtArea->GetXaxis()->SetTitle("area");
273  fHistJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
274  fHistJets2PtArea->GetZaxis()->SetTitle("counts");
276 
277  if (fIsJet2Rho) {
278  fHistJets2CorrPtArea = new TH2F("fHistJets2CorrPtArea", "fHistJets2CorrPtArea",
279  fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
280  fHistJets2CorrPtArea->GetXaxis()->SetTitle("area");
281  fHistJets2CorrPtArea->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)");
282  fHistJets2CorrPtArea->GetZaxis()->SetTitle("counts");
284  }
285 
286  fHistJets2ZvsPt = new TH2F("fHistJets2ZvsPt", "fHistJets2ZvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
287  fHistJets2ZvsPt->GetXaxis()->SetTitle("Z");
288  fHistJets2ZvsPt->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
289  fHistJets2ZvsPt->GetZaxis()->SetTitle("counts");
290  fOutput->Add(fHistJets2ZvsPt);
291 
292  fHistJets2NEFvsPt = new TH2F("fHistJets2NEFvsPt", "fHistJets2NEFvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
293  fHistJets2NEFvsPt->GetXaxis()->SetTitle("NEF");
294  fHistJets2NEFvsPt->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
295  fHistJets2NEFvsPt->GetZaxis()->SetTitle("counts");
297 
298  fHistJets2CEFvsCEFPt = new TH2F("fHistJets2CEFvsCEFPt", "fHistJets2CEFvsCEFPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
299  fHistJets2CEFvsCEFPt->GetXaxis()->SetTitle("1-NEF");
300  fHistJets2CEFvsCEFPt->GetYaxis()->SetTitle("(1-NEF)*p_{T,2} (GeV/c)");
301  fHistJets2CEFvsCEFPt->GetZaxis()->SetTitle("counts");
303 
304  // Matching histograms
305 
306  fHistCommonEnergy1vsJet1Pt = new TH2F("fHistCommonEnergy1vsJet1Pt", "fHistCommonEnergy1vsJet1Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt);
307  fHistCommonEnergy1vsJet1Pt->GetXaxis()->SetTitle("Common energy 1 (%)");
308  fHistCommonEnergy1vsJet1Pt->GetYaxis()->SetTitle("p_{T,1}");
309  fHistCommonEnergy1vsJet1Pt->GetZaxis()->SetTitle("counts");
311 
312  fHistCommonEnergy2vsJet2Pt = new TH2F("fHistCommonEnergy2vsJet2Pt", "fHistCommonEnergy2vsJet2Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt);
313  fHistCommonEnergy2vsJet2Pt->GetXaxis()->SetTitle("Common energy 2 (%)");
314  fHistCommonEnergy2vsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");
315  fHistCommonEnergy2vsJet2Pt->GetZaxis()->SetTitle("counts");
317 
318  fHistDistancevsJet1Pt = new TH2F("fHistDistancevsJet1Pt", "fHistDistancevsJet1Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt);
319  fHistDistancevsJet1Pt->GetXaxis()->SetTitle("Distance");
320  fHistDistancevsJet1Pt->GetYaxis()->SetTitle("p_{T,1}");
321  fHistDistancevsJet1Pt->GetZaxis()->SetTitle("counts");
323 
324  fHistDistancevsJet2Pt = new TH2F("fHistDistancevsJet2Pt", "fHistDistancevsJet2Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt);
325  fHistDistancevsJet2Pt->GetXaxis()->SetTitle("Distance");
326  fHistDistancevsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");
327  fHistDistancevsJet2Pt->GetZaxis()->SetTitle("counts");
329 
330  fHistDistancevsCommonEnergy1 = new TH2F("fHistDistancevsCommonEnergy1", "fHistDistancevsCommonEnergy1", fNbins/2, 0, 1.2, fNbins/2, 0, 1.2);
331  fHistDistancevsCommonEnergy1->GetXaxis()->SetTitle("Distance");
332  fHistDistancevsCommonEnergy1->GetYaxis()->SetTitle("Common energy 1 (%)");
333  fHistDistancevsCommonEnergy1->GetZaxis()->SetTitle("counts");
335 
336  fHistDistancevsCommonEnergy2 = new TH2F("fHistDistancevsCommonEnergy2", "fHistDistancevsCommonEnergy2", fNbins/2, 0, 1.2, fNbins/2, 0, 1.2);
337  fHistDistancevsCommonEnergy2->GetXaxis()->SetTitle("Distance");
338  fHistDistancevsCommonEnergy2->GetYaxis()->SetTitle("Common energy 2 (%)");
339  fHistDistancevsCommonEnergy2->GetZaxis()->SetTitle("counts");
341 
342  fHistCommonEnergy1vsCommonEnergy2 = new TH2F("fHistCommonEnergy1vsCommonEnergy2", "fHistCommonEnergy1vsCommonEnergy2", fNbins/2, 0, 1.2, fNbins/2, 0, 1.2);
343  fHistCommonEnergy1vsCommonEnergy2->GetXaxis()->SetTitle("Common energy 1 (%)");
344  fHistCommonEnergy1vsCommonEnergy2->GetYaxis()->SetTitle("Common energy 2 (%)");
345  fHistCommonEnergy1vsCommonEnergy2->GetZaxis()->SetTitle("counts");
347 
348  fHistDeltaEtaDeltaPhi = new TH2F("fHistDeltaEtaDeltaPhi", "fHistDeltaEtaDeltaPhi", fNbins/4, -1, 1, fNbins/4, -TMath::Pi()/2, TMath::Pi()*3/2);
349  fHistDeltaEtaDeltaPhi->GetXaxis()->SetTitle("Common energy 1 (%)");
350  fHistDeltaEtaDeltaPhi->GetYaxis()->SetTitle("Common energy 2 (%)");
351  fHistDeltaEtaDeltaPhi->GetZaxis()->SetTitle("counts");
353 
354  fHistJet2PtOverJet1PtvsJet2Pt = new TH2F("fHistJet2PtOverJet1PtvsJet2Pt", "fHistJet2PtOverJet1PtvsJet2Pt", fNbins, fMinBinPt, fMaxBinPt, 300, 0, 1.5);
355  fHistJet2PtOverJet1PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");
356  fHistJet2PtOverJet1PtvsJet2Pt->GetYaxis()->SetTitle("p_{T,2} / p_{T,1}");
357  fHistJet2PtOverJet1PtvsJet2Pt->GetZaxis()->SetTitle("counts");
359 
360  fHistJet1PtOverJet2PtvsJet1Pt = new TH2F("fHistJet1PtOverJet2PtvsJet1Pt", "fHistJet1PtOverJet2PtvsJet1Pt", fNbins, fMinBinPt, fMaxBinPt, 300, 0, 1.5);
361  fHistJet1PtOverJet2PtvsJet1Pt->GetXaxis()->SetTitle("p_{T,1}");
362  fHistJet1PtOverJet2PtvsJet1Pt->GetYaxis()->SetTitle("p_{T,1} / p_{T,2}");
363  fHistJet1PtOverJet2PtvsJet1Pt->GetZaxis()->SetTitle("counts");
365 
366  fHistDeltaPtvsJet1Pt = new TH2F("fHistDeltaPtvsJet1Pt", "fHistDeltaPtvsJet1Pt",
368  fHistDeltaPtvsJet1Pt->GetXaxis()->SetTitle("p_{T,1}");
369  fHistDeltaPtvsJet1Pt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
370  fHistDeltaPtvsJet1Pt->GetZaxis()->SetTitle("counts");
372 
373  fHistDeltaPtvsJet2Pt = new TH2F("fHistDeltaPtvsJet2Pt", "fHistDeltaPtvsJet2Pt",
375  fHistDeltaPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");
376  fHistDeltaPtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
377  fHistDeltaPtvsJet2Pt->GetZaxis()->SetTitle("counts");
379 
380  fHistDeltaPtOverJet1PtvsJet1Pt = new TH2F("fHistDeltaPtOverJet1PtvsJet1Pt", "fHistDeltaPtOverJet1PtvsJet1Pt",
381  fNbins, fMinBinPt, fMaxBinPt, fNbins, -5, 5);
382  fHistDeltaPtOverJet1PtvsJet1Pt->GetXaxis()->SetTitle("p_{T,1}");
383  fHistDeltaPtOverJet1PtvsJet1Pt->GetYaxis()->SetTitle("#deltap_{T} / p_{T,1}");
384  fHistDeltaPtOverJet1PtvsJet1Pt->GetZaxis()->SetTitle("counts");
386 
387  fHistDeltaPtOverJet2PtvsJet2Pt = new TH2F("fHistDeltaPtOverJet2PtvsJet2Pt", "fHistDeltaPtOverJet2PtvsJet2Pt",
388  fNbins, fMinBinPt, fMaxBinPt, fNbins, -5, 5);
389  fHistDeltaPtOverJet2PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");
390  fHistDeltaPtOverJet2PtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} / p_{T,2}");
391  fHistDeltaPtOverJet2PtvsJet2Pt->GetZaxis()->SetTitle("counts");
393 
394  fHistDeltaPtvsDistance = new TH2F("fHistDeltaPtvsDistance", "fHistDeltaPtvsDistance",
395  fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
396  fHistDeltaPtvsDistance->GetXaxis()->SetTitle("Distance");
397  fHistDeltaPtvsDistance->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
398  fHistDeltaPtvsDistance->GetZaxis()->SetTitle("counts");
400 
401  fHistDeltaPtvsCommonEnergy1 = new TH2F("fHistDeltaPtvsCommonEnergy1", "fHistDeltaPtvsCommonEnergy1",
402  fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
403  fHistDeltaPtvsCommonEnergy1->GetXaxis()->SetTitle("Common energy 1 (%)");
404  fHistDeltaPtvsCommonEnergy1->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
405  fHistDeltaPtvsCommonEnergy1->GetZaxis()->SetTitle("counts");
407 
408  fHistDeltaPtvsCommonEnergy2 = new TH2F("fHistDeltaPtvsCommonEnergy2", "fHistDeltaPtvsCommonEnergy2",
409  fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
410  fHistDeltaPtvsCommonEnergy2->GetXaxis()->SetTitle("Common energy 2 (%)");
411  fHistDeltaPtvsCommonEnergy2->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
412  fHistDeltaPtvsCommonEnergy2->GetZaxis()->SetTitle("counts");
414 
415  fHistDeltaPtvsArea1 = new TH2F("fHistDeltaPtvsArea1", "fHistDeltaPtvsArea1",
416  fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
417  fHistDeltaPtvsArea1->GetXaxis()->SetTitle("A_{jet,1}");
418  fHistDeltaPtvsArea1->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
419  fHistDeltaPtvsArea1->GetZaxis()->SetTitle("counts");
421 
422  fHistDeltaPtvsArea2 = new TH2F("fHistDeltaPtvsArea2", "fHistDeltaPtvsArea2",
423  fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
424  fHistDeltaPtvsArea2->GetXaxis()->SetTitle("A_{jet,2}");
425  fHistDeltaPtvsArea2->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
426  fHistDeltaPtvsArea2->GetZaxis()->SetTitle("counts");
428 
429  fHistDeltaPtvsDeltaArea = new TH2F("fHistDeltaPtvsDeltaArea", "fHistDeltaPtvsDeltaArea",
430  fNbins, -1.+1./fNbins, 1.+1./fNbins, 2*fNbins, -fMaxBinPt, fMaxBinPt);
431  fHistDeltaPtvsDeltaArea->GetXaxis()->SetTitle("#deltaA_{jet}");
432  fHistDeltaPtvsDeltaArea->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
433  fHistDeltaPtvsDeltaArea->GetZaxis()->SetTitle("counts");
435 
436  fHistJet1PtvsJet2Pt = new TH2F("fHistJet1PtvsJet2Pt", "fHistJet1PtvsJet2Pt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
437  fHistJet1PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,1}");
438  fHistJet1PtvsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");
439  fHistJet1PtvsJet2Pt->GetZaxis()->SetTitle("counts");
441 
442  if (fIsJet1Rho || fIsJet2Rho) {
443  if (!fIsJet1Rho)
444  fHistDeltaCorrPtvsJet1CorrPt = new TH2F("fHistDeltaCorrPtvsJet1CorrPt", "fHistDeltaCorrPtvsJet1CorrPt",
446  else
447  fHistDeltaCorrPtvsJet1CorrPt = new TH2F("fHistDeltaCorrPtvsJet1CorrPt", "fHistDeltaCorrPtvsJet1CorrPt",
449 
450  fHistDeltaCorrPtvsJet1CorrPt->GetXaxis()->SetTitle("p_{T,1}^{corr}");
451  fHistDeltaCorrPtvsJet1CorrPt->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
452  fHistDeltaCorrPtvsJet1CorrPt->GetZaxis()->SetTitle("counts");
454 
455  if (!fIsJet2Rho)
456  fHistDeltaCorrPtvsJet2CorrPt = new TH2F("fHistDeltaCorrPtvsJet2CorrPt", "fHistDeltaCorrPtvsJet2CorrPt",
458  else
459  fHistDeltaCorrPtvsJet2CorrPt = new TH2F("fHistDeltaCorrPtvsJet2CorrPt", "fHistDeltaCorrPtvsJet2CorrPt",
461 
462  fHistDeltaCorrPtvsJet2CorrPt->GetXaxis()->SetTitle("p_{T,2}^{corr}");
463  fHistDeltaCorrPtvsJet2CorrPt->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
464  fHistDeltaCorrPtvsJet2CorrPt->GetZaxis()->SetTitle("counts");
466 
467  fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt = new TH2F("fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt", "fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt",
468  2*fNbins, -fMaxBinPt, fMaxBinPt, fNbins, -5, 5);
469  fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt->GetXaxis()->SetTitle("p_{T,1}^{corr}");
470  fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt->GetYaxis()->SetTitle("#deltap_{T}^{corr} / p_{T,1}^{corr}");
471  fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt->GetZaxis()->SetTitle("counts");
473 
474  fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt = new TH2F("fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt", "fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt",
475  2*fNbins, -fMaxBinPt, fMaxBinPt, fNbins, -5, 5);
476  fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt->GetXaxis()->SetTitle("p_{T,2}^{corr}");
477  fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt->GetYaxis()->SetTitle("#deltap_{T}^{corr} / p_{T,2}^{corr}");
478  fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt->GetZaxis()->SetTitle("counts");
480 
481  fHistDeltaCorrPtvsDistance = new TH2F("fHistDeltaCorrPtvsDistance", "fHistDeltaCorrPtvsDistance",
482  fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
483  fHistDeltaCorrPtvsDistance->GetXaxis()->SetTitle("Distance");
484  fHistDeltaCorrPtvsDistance->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
485  fHistDeltaCorrPtvsDistance->GetZaxis()->SetTitle("counts");
487 
488  fHistDeltaCorrPtvsCommonEnergy1 = new TH2F("fHistDeltaCorrPtvsCommonEnergy1", "fHistDeltaCorrPtvsCommonEnergy1",
489  fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
490  fHistDeltaCorrPtvsCommonEnergy1->GetXaxis()->SetTitle("Common energy 1 (%)");
491  fHistDeltaCorrPtvsCommonEnergy1->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
492  fHistDeltaCorrPtvsCommonEnergy1->GetZaxis()->SetTitle("counts");
494 
495  fHistDeltaCorrPtvsCommonEnergy2 = new TH2F("fHistDeltaCorrPtvsCommonEnergy2", "fHistDeltaCorrPtvsCommonEnergy2",
496  fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
497  fHistDeltaCorrPtvsCommonEnergy2->GetXaxis()->SetTitle("Common energy 2 (%)");
498  fHistDeltaCorrPtvsCommonEnergy2->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
499  fHistDeltaCorrPtvsCommonEnergy2->GetZaxis()->SetTitle("counts");
501 
502  fHistDeltaCorrPtvsArea1 = new TH2F("fHistDeltaCorrPtvsArea1", "fHistDeltaCorrPtvsArea1",
503  fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
504  fHistDeltaCorrPtvsArea1->GetXaxis()->SetTitle("A_{jet,1}");
505  fHistDeltaCorrPtvsArea1->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
506  fHistDeltaCorrPtvsArea1->GetZaxis()->SetTitle("counts");
508 
509  fHistDeltaCorrPtvsArea2 = new TH2F("fHistDeltaCorrPtvsArea2", "fHistDeltaCorrPtvsArea2",
510  fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
511  fHistDeltaCorrPtvsArea2->GetXaxis()->SetTitle("A_{jet,2}");
512  fHistDeltaCorrPtvsArea2->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
513  fHistDeltaCorrPtvsArea2->GetZaxis()->SetTitle("counts");
515 
516  fHistDeltaCorrPtvsDeltaArea = new TH2F("fHistDeltaCorrPtvsDeltaArea", "fHistDeltaCorrPtvsDeltaArea",
517  fNbins, -1.+1./fNbins, 1.+1./fNbins, 2*fNbins, -fMaxBinPt, fMaxBinPt);
518  fHistDeltaCorrPtvsDeltaArea->GetXaxis()->SetTitle("#deltaA_{jet}");
519  fHistDeltaCorrPtvsDeltaArea->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
520  fHistDeltaCorrPtvsDeltaArea->GetZaxis()->SetTitle("counts");
522 
523  if (!fIsJet1Rho)
524  fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt",
526  else if (!fIsJet2Rho)
527  fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt",
529  else
530  fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt",
532  2*fNbins, -fMaxBinPt, fMaxBinPt);
533  fHistJet1CorrPtvsJet2CorrPt->GetXaxis()->SetTitle("p_{T,1}^{corr}");
534  fHistJet1CorrPtvsJet2CorrPt->GetYaxis()->SetTitle("p_{T,2}^{corr}");
535  fHistJet1CorrPtvsJet2CorrPt->GetZaxis()->SetTitle("counts");
537  }
538 
539  if (fIsEmbedded) {
540  fHistDeltaMCPtvsJet1MCPt = new TH2F("fHistDeltaMCPtvsJet1MCPt", "fHistDeltaMCPtvsJet1MCPt",
542  fHistDeltaMCPtvsJet1MCPt->GetXaxis()->SetTitle("p_{T,1}^{MC}");
543  fHistDeltaMCPtvsJet1MCPt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
544  fHistDeltaMCPtvsJet1MCPt->GetZaxis()->SetTitle("counts");
546 
547  fHistDeltaMCPtvsJet2Pt = new TH2F("fHistDeltaMCPtvsJet2Pt", "fHistDeltaMCPtvsJet2Pt",
549  fHistDeltaMCPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");
550  fHistDeltaMCPtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
551  fHistDeltaMCPtvsJet2Pt->GetZaxis()->SetTitle("counts");
553 
554  fHistDeltaMCPtOverJet1MCPtvsJet1MCPt = new TH2F("fHistDeltaMCPtOverJet1MCPtvsJet1MCPt", "fHistDeltaMCPtOverJet1MCPtvsJet1MCPt",
555  fNbins, fMinBinPt, fMaxBinPt, fNbins, -5, 5);
556  fHistDeltaMCPtOverJet1MCPtvsJet1MCPt->GetXaxis()->SetTitle("p_{T,1}^{MC}");
557  fHistDeltaMCPtOverJet1MCPtvsJet1MCPt->GetYaxis()->SetTitle("#deltap_{T} / p_{T,1}^{MC}");
558  fHistDeltaMCPtOverJet1MCPtvsJet1MCPt->GetZaxis()->SetTitle("counts");
560 
561  fHistDeltaMCPtOverJet2PtvsJet2Pt = new TH2F("fHistDeltaMCPtOverJet2PtvsJet2Pt", "fHistDeltaMCPtOverJet2PtvsJet2Pt",
562  fNbins, fMinBinPt, fMaxBinPt, fNbins, -5, 5);
563  fHistDeltaMCPtOverJet2PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");
564  fHistDeltaMCPtOverJet2PtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} / p_{T,2}");
565  fHistDeltaMCPtOverJet2PtvsJet2Pt->GetZaxis()->SetTitle("counts");
567 
568  fHistDeltaMCPtvsDistance = new TH2F("fHistDeltaMCPtvsDistance", "fHistDeltaMCPtvsDistance",
569  fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
570  fHistDeltaMCPtvsDistance->GetXaxis()->SetTitle("Distance");
571  fHistDeltaMCPtvsDistance->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
572  fHistDeltaMCPtvsDistance->GetZaxis()->SetTitle("counts");
574 
575  fHistDeltaMCPtvsCommonEnergy1 = new TH2F("fHistDeltaMCPtvsCommonEnergy1", "fHistDeltaMCPtvsCommonEnergy1",
576  fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
577  fHistDeltaMCPtvsCommonEnergy1->GetXaxis()->SetTitle("Common energy 1 (%)");
578  fHistDeltaMCPtvsCommonEnergy1->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
579  fHistDeltaMCPtvsCommonEnergy1->GetZaxis()->SetTitle("counts");
581 
582  fHistDeltaMCPtvsCommonEnergy2 = new TH2F("fHistDeltaMCPtvsCommonEnergy2", "fHistDeltaMCPtvsCommonEnergy2",
583  fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
584  fHistDeltaMCPtvsCommonEnergy2->GetXaxis()->SetTitle("Common energy 2 (%)");
585  fHistDeltaMCPtvsCommonEnergy2->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
586  fHistDeltaMCPtvsCommonEnergy2->GetZaxis()->SetTitle("counts");
588 
589  fHistDeltaMCPtvsArea1 = new TH2F("fHistDeltaMCPtvsArea1", "fHistDeltaMCPtvsArea1",
590  fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
591  fHistDeltaMCPtvsArea1->GetXaxis()->SetTitle("A_{jet,1}");
592  fHistDeltaMCPtvsArea1->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
593  fHistDeltaMCPtvsArea1->GetZaxis()->SetTitle("counts");
595 
596  fHistDeltaMCPtvsArea2 = new TH2F("fHistDeltaMCPtvsArea2", "fHistDeltaMCPtvsArea2",
597  fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
598  fHistDeltaMCPtvsArea2->GetXaxis()->SetTitle("A_{jet,2}");
599  fHistDeltaMCPtvsArea2->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
600  fHistDeltaMCPtvsArea2->GetZaxis()->SetTitle("counts");
602 
603  fHistDeltaMCPtvsDeltaArea = new TH2F("fHistDeltaMCPtvsDeltaArea", "fHistDeltaMCPtvsDeltaArea",
604  fNbins, -1.+1./fNbins, 1.+1./fNbins, 2*fNbins, -fMaxBinPt, fMaxBinPt);
605  fHistDeltaMCPtvsDeltaArea->GetXaxis()->SetTitle("#deltaA_{jet}");
606  fHistDeltaMCPtvsDeltaArea->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
607  fHistDeltaMCPtvsDeltaArea->GetZaxis()->SetTitle("counts");
609 
610  fHistJet1MCPtvsJet2Pt = new TH2F("fHistJet1MCPtvsJet2Pt", "fHistJet1MCPtvsJet2Pt",
612  fHistJet1MCPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,1}^{MC}");
613  fHistJet1MCPtvsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");
614  fHistJet1MCPtvsJet2Pt->GetZaxis()->SetTitle("counts");
616  }
617 }
618 
619 //________________________________________________________________________
621 {
622  // Allocate THnSparse histograms.
623 
624  TString title[25]= {""};
625  Int_t nbins[25] = {0};
626  Double_t min[25] = {0.};
627  Double_t max[25] = {0.};
628  Int_t dim = 0;
629 
630  title[dim] = "#phi";
631  nbins[dim] = fNbins/4;
632  min[dim] = 0;
633  max[dim] = 2*TMath::Pi()*(1 + 1./(nbins[dim]-1));
634  dim++;
635 
636  title[dim] = "#eta";
637  nbins[dim] = fNbins/4;
638  min[dim] = -1;
639  max[dim] = 1;
640  dim++;
641 
642  title[dim] = "p_{T}";
643  nbins[dim] = fNbins;
644  min[dim] = 0;
645  max[dim] = 250;
646  dim++;
647 
648  title[dim] = "A_{jet}";
649  nbins[dim] = fNbins/4;
650  min[dim] = 0;
651  max[dim] = 1.5;
652  dim++;
653 
654  if (fNEFAxis) {
655  title[dim] = "NEF";
656  nbins[dim] = 120;
657  min[dim] = 0.0;
658  max[dim] = 1.2;
659  dim++;
660  }
661 
662  if (fZAxis) {
663  title[dim] = "Z";
664  nbins[dim] = 120;
665  min[dim] = 0.0;
666  max[dim] = 1.2;
667  dim++;
668  }
669 
670  if (fFlavourZAxis) {
671  title[dim] = "z_{flavour}";
672  nbins[dim] = 100;
673  min[dim] = 0.0;
674  max[dim] = 2.0;
675  dim++;
676  }
677 
678  if (fFlavourPtAxis) {
679  title[dim] = "p_{T}^{D}";
680  nbins[dim] = fNbins/2;
681  min[dim] = 0;
682  max[dim] = 125;
683  dim++;
684  }
685 
686  if (fJetRelativeEPAngle) {
687  title[dim] = "#theta_{jet}^{EP}";
688  nbins[dim] = 3;
689  min[dim] = 0;
690  max[dim] = TMath::Pi()/2;
691  dim++;
692  }
693 
694  title[dim] = "p_{T,particle}^{leading} (GeV/c)";
695  nbins[dim] = 120;
696  min[dim] = 0;
697  max[dim] = 120;
698  dim++;
699 
700  Int_t dim1 = dim, dim2 = dim;
701 
702  if (fIsJet1Rho) {
703  title[dim1] = "p_{T}^{corr}";
704  nbins[dim1] = fNbins*2;
705  min[dim1] = -250;
706  max[dim1] = 250;
707  dim1++;
708  }
709 
710  if (fIsEmbedded) {
711  title[dim1] = "p_{T}^{MC}";
712  nbins[dim1] = fNbins;
713  min[dim1] = 0;
714  max[dim1] = 250;
715  dim1++;
716  }
717 
718  fHistJets1 = new THnSparseD("fHistJets1","fHistJets1",dim1,nbins,min,max);
719  for (Int_t i = 0; i < dim1; i++)
720  fHistJets1->GetAxis(i)->SetTitle(title[i]);
721  fOutput->Add(fHistJets1);
722 
723  if (fIsJet2Rho) {
724  title[dim2] = "p_{T}^{corr}";
725  nbins[dim2] = fNbins*2;
726  min[dim2] = -250;
727  max[dim2] = 250;
728  dim2++;
729  }
730 
731  fHistJets2 = new THnSparseD("fHistJets2","fHistJets2",dim2,nbins,min,max);
732  for (Int_t i = 0; i < dim2; i++)
733  fHistJets2->GetAxis(i)->SetTitle(title[i]);
734  fOutput->Add(fHistJets2);
735 
736  // Matching
737 
738  dim = 0;
739 
740  title[dim] = "p_{T,1}";
741  nbins[dim] = fNbins;
742  min[dim] = 0;
743  max[dim] = 250;
744  dim++;
745 
746  title[dim] = "p_{T,2}";
747  nbins[dim] = fNbins;
748  min[dim] = 0;
749  max[dim] = 250;
750  dim++;
751 
752  title[dim] = "A_{jet,1}";
753  nbins[dim] = fNbins/4;
754  min[dim] = 0;
755  max[dim] = 1.5;
756  dim++;
757 
758  title[dim] = "A_{jet,2}";
759  nbins[dim] = fNbins/4;
760  min[dim] = 0;
761  max[dim] = 1.5;
762  dim++;
763 
764  title[dim] = "distance";
765  nbins[dim] = fNbins/2;
766  min[dim] = 0;
767  max[dim] = 1.2;
768  dim++;
769 
770  title[dim] = "CE1";
771  nbins[dim] = fNbins/2;
772  min[dim] = 0;
773  max[dim] = 1.2;
774  dim++;
775 
776  title[dim] = "CE2";
777  nbins[dim] = fNbins/2;
778  min[dim] = 0;
779  max[dim] = 1.2;
780  dim++;
781 
782  title[dim] = "p_{T,particle,1}^{leading} (GeV/c)";
783  nbins[dim] = 120;
784  min[dim] = 0;
785  max[dim] = 120;
786  dim++;
787 
788  title[dim] = "p_{T,particle,2}^{leading} (GeV/c)";
789  nbins[dim] = 120;
790  min[dim] = 0;
791  max[dim] = 120;
792  dim++;
793 
794  if (fDeltaPtAxis) {
795  title[dim] = "#deltaA_{jet}";
796  nbins[dim] = fNbins/2;
797  min[dim] = -1;
798  max[dim] = 1;
799  dim++;
800 
801  title[dim] = "#deltap_{T}";
802  nbins[dim] = fNbins*2;
803  min[dim] = -250;
804  max[dim] = 250;
805  dim++;
806  }
807  if (fIsJet1Rho) {
808  title[dim] = "p_{T,1}^{corr}";
809  nbins[dim] = fNbins*2;
810  min[dim] = -250;
811  max[dim] = 250;
812  dim++;
813  }
814  if (fIsJet2Rho) {
815  title[dim] = "p_{T,2}^{corr}";
816  nbins[dim] = fNbins*2;
817  min[dim] = -250;
818  max[dim] = 250;
819  dim++;
820  }
821  if (fDeltaPtAxis && (fIsJet1Rho || fIsJet2Rho)) {
822  title[dim] = "#deltap_{T}^{corr}";
823  nbins[dim] = fNbins*2;
824  min[dim] = -250;
825  max[dim] = 250;
826  dim++;
827  }
828  if (fDeltaEtaDeltaPhiAxis) {
829  title[dim] = "#delta#eta";
830  nbins[dim] = fNbins/2;
831  min[dim] = -1;
832  max[dim] = 1;
833  dim++;
834 
835  title[dim] = "#delta#phi";
836  nbins[dim] = fNbins/2;
837  min[dim] = -TMath::Pi()/2;
838  max[dim] = TMath::Pi()*3/2;
839  dim++;
840  }
841  if (fIsEmbedded) {
842  title[dim] = "p_{T,1}^{MC}";
843  nbins[dim] = fNbins;
844  min[dim] = 0;
845  max[dim] = 250;
846  dim++;
847 
848  if (fDeltaPtAxis) {
849  title[dim] = "#deltap_{T}^{MC}";
850  nbins[dim] = fNbins*2;
851  min[dim] = -250;
852  max[dim] = 250;
853  dim++;
854  }
855  }
856 
857  if (fNEFAxis) {
858  title[dim] = "NEF_{1}";
859  nbins[dim] = 120;
860  min[dim] = 0.0;
861  max[dim] = 1.2;
862  dim++;
863 
864  title[dim] = "NEF_{2}";
865  nbins[dim] = 120;
866  min[dim] = 0.0;
867  max[dim] = 1.2;
868  dim++;
869  }
870 
871  if (fZAxis) {
872  title[dim] = "Z_{1}";
873  nbins[dim] = 120;
874  min[dim] = 0.0;
875  max[dim] = 1.2;
876  dim++;
877 
878  title[dim] = "Z_{2}";
879  nbins[dim] = 120;
880  min[dim] = 0.0;
881  max[dim] = 1.2;
882  dim++;
883  }
884 
885  if (fFlavourZAxis) {
886  title[dim] = "z_{flavour,1}";
887  nbins[dim] = 100;
888  min[dim] = 0.0;
889  max[dim] = 2.0;
890  dim++;
891 
892  title[dim] = "z_{flavour,2}";
893  nbins[dim] = 100;
894  min[dim] = 0.0;
895  max[dim] = 2.0;
896  dim++;
897  }
898 
899  if (fFlavourPtAxis) {
900  title[dim] = "p_{T,1}^{D}";
901  nbins[dim] = fNbins/2;
902  min[dim] = 0;
903  max[dim] = 125;
904  dim++;
905 
906  title[dim] = "p_{T,2}^{D}";
907  nbins[dim] = fNbins/2;
908  min[dim] = 0;
909  max[dim] = 125;
910  dim++;
911  }
912 
913  if (fZgAxis) {
914  title[dim] = "Z_{g,1}";
915  nbins[dim] = 20;
916  min[dim] = 0.0;
917  max[dim] = 1.0;
918  dim++;
919  title[dim] = "Z_{g,2}";
920  nbins[dim] = 20;
921  min[dim] = 0.0;
922  max[dim] = 1.0;
923  dim++;
924  }
925 
926  if (fdRAxis) {
927  title[dim] = "dR_{1}";
928  nbins[dim] = 40;
929  min[dim] = 0.0;
930  max[dim] = 0.5;
931  dim++;
932  title[dim] = "dR_{2}";
933  nbins[dim] = 40;
934  min[dim] = 0.0;
935  max[dim] = 0.5;
936  dim++;
937  }
938 
939  if (fPtgAxis) {
940  title[dim] = "p_{T,g,1}";
941  nbins[dim] = 16;
942  min[dim] = 0.0;
943  max[dim] = 160.0;
944  dim++;
945  title[dim] = "p_{T,g,2}";
946  nbins[dim] = 16;
947  min[dim] = 0.0;
948  max[dim] = 160.0;
949  dim++;
950  }
951 
952  if (fJetRelativeEPAngle) {
953  title[dim] = "#theta_{jet,1}^{EP}";
954  nbins[dim] = 3;
955  min[dim] = 0;
956  max[dim] = TMath::Pi()/2;
957  dim++;
958 
959  title[dim] = "#theta_{jet,2}^{EP}";
960  nbins[dim] = 3;
961  min[dim] = 0;
962  max[dim] = TMath::Pi()/2;
963  dim++;
964  }
965 
966  fHistMatching = new THnSparseD("fHistMatching","fHistMatching",dim,nbins,min,max);
967 
968  for (Int_t i = 0; i < dim; i++)
969  fHistMatching->GetAxis(i)->SetTitle(title[i]);
970 
971  fOutput->Add(fHistMatching);
972 }
973 
974 //________________________________________________________________________
976 {
977  // Create user objects.
978 
980 
981  AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
982  AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
983 
984  if (!jets1 || !jets2) return;
985 
986  if (jets1->GetRhoName().IsNull()) fIsJet1Rho = kFALSE;
987  else fIsJet1Rho = kTRUE;
988  if (jets2->GetRhoName().IsNull()) fIsJet2Rho = kFALSE;
989  else fIsJet2Rho = kTRUE;
990 
991  fHistRejectionReason1 = new TH2F("fHistRejectionReason1", "fHistRejectionReason1", 32, 0, 32, 100, 0, 250);
992  fHistRejectionReason1->GetXaxis()->SetTitle("Rejection reason");
993  fHistRejectionReason1->GetYaxis()->SetTitle("p_{T,jet} (GeV/c)");
994  fHistRejectionReason1->GetZaxis()->SetTitle("counts");
997 
998  fHistRejectionReason2 = new TH2F("fHistRejectionReason2", "fHistRejectionReason2", 32, 0, 32, 100, 0, 250);
999  fHistRejectionReason2->GetXaxis()->SetTitle("Rejection reason");
1000  fHistRejectionReason2->GetYaxis()->SetTitle("p_{T,jet} (GeV/c)");
1001  fHistRejectionReason2->GetZaxis()->SetTitle("counts");
1004 
1005  if (fHistoType==0)
1006  AllocateTH2();
1007  else
1009 
1010  PostData(1, fOutput); // Post data for ALL output slots > 0 here, to get at least an empty histogram
1011 }
1012 
1013 //________________________________________________________________________
1015 {
1016  AliJetContainer* jets = GetJetContainer(Set-1);
1017 
1018  AliTLorentzVector leadPart;
1019  jets->GetLeadingHadronMomentum(leadPart, jet);
1020  Double_t zleading = GetParallelFraction(leadPart.Vect(), jet);
1021  if (zleading == 1 || (zleading > 1 && zleading - 1 < 1e-3)) zleading = 0.999; // so that it will contribute to the bin 0.9-1 rather than 1-1.1
1022 
1023  Double_t corrpt = jet->Pt() - jets->GetRhoVal() * jet->Area();
1024  Double_t zflavour = 0;
1025  Double_t ptflavour = 0;
1026  AliVParticle* hftrack = jet->GetFlavourTrack();
1027  if (hftrack) {
1028  zflavour = GetParallelFraction(hftrack, jet);
1029  ptflavour = hftrack->Pt();
1030 
1031  if (zflavour == 1 || (zflavour > 1 && zflavour - 1 < 1e-3)) zflavour = 0.999; // so that it will contribute to the bin 0.9-1 rather than 1-1.1
1032  }
1033  Double_t jetRelativeEPAngle = GetRelativeEPAngle(jet->Phi(), fEPV0);
1034 
1035  if (fHistoType==1) {
1036  THnSparse *histo = 0;
1037  if (Set==1) {
1038  histo = fHistJets1;
1039  }
1040  else if (Set==2) {
1041  histo = fHistJets2;
1042  }
1043 
1044  if (!histo) return;
1045 
1046  Double_t contents[20]={0};
1047 
1048  for (Int_t i = 0; i < histo->GetNdimensions(); i++) {
1049  TString title(histo->GetAxis(i)->GetTitle());
1050  if (title=="#phi")
1051  contents[i] = jet->Phi();
1052  else if (title=="#eta")
1053  contents[i] = jet->Eta();
1054  else if (title=="p_{T}")
1055  contents[i] = jet->Pt();
1056  else if (title=="A_{jet}")
1057  contents[i] = jet->Area();
1058  else if (title=="NEF")
1059  contents[i] = jet->NEF();
1060  else if (title=="Z")
1061  contents[i] = zleading;
1062  else if (title=="p_{T}^{corr}")
1063  contents[i] = corrpt;
1064  else if (title=="p_{T}^{MC}")
1065  contents[i] = jet->MCPt();
1066  else if (title=="p_{T,particle}^{leading} (GeV/c)")
1067  contents[i] = leadPart.Pt();
1068  else if (title=="z_{flavour}")
1069  contents[i] = zflavour;
1070  else if (title=="p_{T}^{D}")
1071  contents[i] = ptflavour;
1072  else if (title=="#theta_{jet}^{EP}")
1073  contents[i] = jetRelativeEPAngle;
1074  else
1075  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1076  }
1077 
1078  histo->Fill(contents);
1079  }
1080  else {
1081  if (Set == 1) {
1082  fHistJets1PtArea->Fill(jet->Area(), jet->Pt());
1083  fHistJets1PhiEta->Fill(jet->Eta(), jet->Phi());
1084  fHistJets1ZvsPt->Fill(zleading, jet->Pt());
1085  fHistJets1NEFvsPt->Fill(jet->NEF(), jet->Pt());
1086  fHistJets1CEFvsCEFPt->Fill(1-jet->NEF(), (1-jet->NEF())*jet->Pt());
1087  if (fIsJet1Rho)
1088  fHistJets1CorrPtArea->Fill(jet->Area(), corrpt);
1089  }
1090  else if (Set == 2) {
1091  fHistJets2PtArea->Fill(jet->Area(), jet->Pt());
1092  fHistJets2PhiEta->Fill(jet->Eta(), jet->Phi());
1093  fHistJets2ZvsPt->Fill(zleading, jet->Pt());
1094  fHistJets2NEFvsPt->Fill(jet->NEF(), jet->Pt());
1095  fHistJets2CEFvsCEFPt->Fill(1-jet->NEF(), (1-jet->NEF())*jet->Pt());
1096  if (fIsJet2Rho)
1097  fHistJets2CorrPtArea->Fill(jet->Area(), corrpt);
1098  }
1099  }
1100 }
1101 
1102 //________________________________________________________________________
1104 {
1105  AliJetContainer* jets1 = GetJetContainer(0);
1106  AliJetContainer* jets2 = GetJetContainer(1);
1107 
1108  /*
1109  Printf("Detector level jet");
1110  jet1->Print();
1111 
1112  jet1->PrintConstituents(jets1->GetParticleContainer() != 0 ? jets1->GetParticleContainer()->GetArray() : 0,
1113  jets1->GetClusterContainer() != 0 ? jets1->GetClusterContainer()->GetArray() : 0);
1114 
1115  Printf("Matched with particle level jet");
1116  jet2->Print();
1117  jet2->PrintConstituents(jets2->GetParticleContainer() != 0 ? jets2->GetParticleContainer()->GetArray() : 0,
1118  jets2->GetClusterContainer() != 0 ? jets2->GetClusterContainer()->GetArray() : 0);
1119  */
1120 
1121  AliTLorentzVector leadPart1;
1122  jets1->GetLeadingHadronMomentum(leadPart1, jet1);
1123  Double_t zleading1 = GetParallelFraction(leadPart1.Vect(), jet1);
1124  if (zleading1 == 1 || (zleading1 > 1 && zleading1 - 1 < 1e-3)) zleading1 = 0.999; // so that it will contribute to the bin 0.9-1 rather than 1-1.1
1125 
1126  Double_t corrpt1 = jet1->Pt() - jets1->GetRhoVal() * jet1->Area();
1127  Double_t zflavour1 = 0;
1128  Double_t ptflavour1 = 0;
1129  AliVParticle* hftrack1 = jet1->GetFlavourTrack();
1130  if (hftrack1) {
1131  zflavour1 = GetParallelFraction(hftrack1, jet1);
1132  ptflavour1 = hftrack1->Pt();
1133 
1134  if (zflavour1 == 1 || (zflavour1 > 1 && zflavour1 - 1 < 1e-3)) zflavour1 = 0.999; // so that it will contribute to the bin 0.9-1 rather than 1-1.1
1135  }
1136  Double_t jetRelativeEPAngle1 = GetRelativeEPAngle(jet1->Phi(), fEPV0);
1137 
1138 
1139  AliTLorentzVector leadPart2;
1140  jets2->GetLeadingHadronMomentum(leadPart2, jet2);
1141  Double_t zleading2 = GetParallelFraction(leadPart2.Vect(), jet2);
1142  if (zleading2 == 1 || (zleading2 > 1 && zleading2 - 1 < 1e-3)) zleading2 = 0.999; // so that it will contribute to the bin 0.9-1 rather than 1-1.1
1143 
1144  Double_t corrpt2 = jet2->Pt() - jets2->GetRhoVal() * jet2->Area();
1145  Double_t zflavour2 = 0;
1146  Double_t ptflavour2 = 0;
1147  AliVParticle* hftrack2 = jet2->GetFlavourTrack();
1148  if (hftrack2) {
1149  zflavour2 = GetParallelFraction(hftrack2, jet2);
1150  ptflavour2 = hftrack2->Pt();
1151 
1152  if (zflavour2 == 1 || (zflavour2 > 1 && zflavour2 - 1 < 1e-3)) zflavour2 = 0.999; // so that it will contribute to the bin 0.9-1 rather than 1-1.1
1153  }
1154  Double_t jetRelativeEPAngle2 = GetRelativeEPAngle(jet2->Phi(), fEPV0);
1155 
1156  if (fHistoType==1) {
1157  Double_t contents[20]={0};
1158 
1159  for (Int_t i = 0; i < fHistMatching->GetNdimensions(); i++) {
1160  TString title(fHistMatching->GetAxis(i)->GetTitle());
1161  if (title=="p_{T,1}")
1162  contents[i] = jet1->Pt();
1163  else if (title=="p_{T,2}")
1164  contents[i] = jet2->Pt();
1165  else if (title=="A_{jet,1}")
1166  contents[i] = jet1->Area();
1167  else if (title=="A_{jet,2}")
1168  contents[i] = jet2->Area();
1169  else if (title=="distance")
1170  contents[i] = d;
1171  else if (title=="CE1")
1172  contents[i] = CE1;
1173  else if (title=="CE2")
1174  contents[i] = CE2;
1175  else if (title=="#deltaA_{jet}")
1176  contents[i] = jet1->Area()-jet2->Area();
1177  else if (title=="#deltap_{T}")
1178  contents[i] = jet1->Pt()-jet2->Pt();
1179  else if (title=="#delta#eta")
1180  contents[i] = jet1->Eta()-jet2->Eta();
1181  else if (title=="#delta#phi")
1182  contents[i] = jet1->Phi()-jet2->Phi();
1183  else if (title=="p_{T,1}^{corr}")
1184  contents[i] = corrpt1;
1185  else if (title=="p_{T,2}^{corr}")
1186  contents[i] = corrpt2;
1187  else if (title=="#deltap_{T}^{corr}")
1188  contents[i] = corrpt1-corrpt2;
1189  else if (title=="p_{T,1}^{MC}")
1190  contents[i] = jet1->MCPt();
1191  else if (title=="#deltap_{T}^{MC}")
1192  contents[i] = jet1->MCPt()-jet2->Pt();
1193  else if (title=="NEF_{1}")
1194  contents[i] = jet1->NEF();
1195  else if (title=="NEF_{2}")
1196  contents[i] = jet2->NEF();
1197  else if (title=="Z_{1}")
1198  contents[i] = zleading1;
1199  else if (title=="Z_{2}")
1200  contents[i] = zleading2;
1201  else if (title=="p_{T,particle,1}^{leading} (GeV/c)")
1202  contents[i] = leadPart1.Pt();
1203  else if (title=="p_{T,particle,2}^{leading} (GeV/c)")
1204  contents[i] = leadPart2.Pt();
1205  else if (title=="z_{flavour,1}")
1206  contents[i] = zflavour1;
1207  else if (title=="z_{flavour,2}")
1208  contents[i] = zflavour2;
1209  else if (title=="p_{T,1}^{D}")
1210  contents[i] = ptflavour1;
1211  else if (title=="p_{T,2}^{D}")
1212  contents[i] = ptflavour2;
1213  else if (title=="Z_{g,1}")
1214  contents[i] = jet1->GetShapeProperties()->GetSoftDropZg();
1215  else if (title=="Z_{g,2}")
1216  contents[i] = jet2->GetShapeProperties()->GetSoftDropZg();
1217  else if (title=="dR_{1}")
1218  contents[i] = jet1->GetShapeProperties()->GetSoftDropdR();
1219  else if (title=="dR_{2}")
1220  contents[i] = jet2->GetShapeProperties()->GetSoftDropdR();
1221  else if (title=="p_{T,g,1}")
1222  contents[i] = ( jet1->GetShapeProperties()->GetSoftDropPtfrac() )*( jet1->Pt() );
1223  else if (title=="p_{T,g,2}")
1224  contents[i] = ( jet2->GetShapeProperties()->GetSoftDropPtfrac() )*( jet2->Pt() );
1225  else if (title=="#theta_{jet,1}^{EP}")
1226  contents[i] = jetRelativeEPAngle1;
1227  else if (title=="#theta_{jet,2}^{EP}")
1228  contents[i] = jetRelativeEPAngle2;
1229  else
1230  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1231  }
1232 
1233  fHistMatching->Fill(contents);
1234  }
1235  else {
1236  fHistCommonEnergy1vsJet1Pt->Fill(CE1, jet1->Pt());
1237  fHistCommonEnergy2vsJet2Pt->Fill(CE2, jet2->Pt());
1238 
1239  fHistDistancevsJet1Pt->Fill(d, jet1->Pt());
1240  fHistDistancevsJet2Pt->Fill(d, jet2->Pt());
1241 
1242  fHistDistancevsCommonEnergy1->Fill(d, CE1);
1243  fHistDistancevsCommonEnergy2->Fill(d, CE2);
1244  fHistCommonEnergy1vsCommonEnergy2->Fill(CE1, CE2);
1245 
1246  fHistDeltaEtaDeltaPhi->Fill(jet1->Eta()-jet2->Eta(),jet1->Phi()-jet2->Phi());
1247 
1248  fHistJet2PtOverJet1PtvsJet2Pt->Fill(jet2->Pt(), jet2->Pt() / jet1->Pt());
1249  fHistJet1PtOverJet2PtvsJet1Pt->Fill(jet1->Pt(), jet1->Pt() / jet2->Pt());
1250 
1251  Double_t dpt = jet1->Pt() - jet2->Pt();
1252  fHistDeltaPtvsJet1Pt->Fill(jet1->Pt(), dpt);
1253  fHistDeltaPtvsJet2Pt->Fill(jet2->Pt(), dpt);
1254  fHistDeltaPtOverJet1PtvsJet1Pt->Fill(jet1->Pt(), dpt/jet1->Pt());
1255  fHistDeltaPtOverJet2PtvsJet2Pt->Fill(jet2->Pt(), dpt/jet2->Pt());
1256 
1257  fHistDeltaPtvsDistance->Fill(d, dpt);
1258  fHistDeltaPtvsCommonEnergy1->Fill(CE1, dpt);
1259  fHistDeltaPtvsCommonEnergy2->Fill(CE2, dpt);
1260 
1261  fHistDeltaPtvsArea1->Fill(jet1->Area(), dpt);
1262  fHistDeltaPtvsArea2->Fill(jet2->Area(), dpt);
1263 
1264  Double_t darea = jet1->Area() - jet2->Area();
1265  fHistDeltaPtvsDeltaArea->Fill(darea, dpt);
1266 
1267  fHistJet1PtvsJet2Pt->Fill(jet1->Pt(), jet2->Pt());
1268 
1269  if (fIsJet1Rho || fIsJet2Rho) {
1270  Double_t dcorrpt = corrpt1 - corrpt2;
1271  fHistDeltaCorrPtvsJet1CorrPt->Fill(corrpt1, dcorrpt);
1272  fHistDeltaCorrPtvsJet2CorrPt->Fill(corrpt2, dcorrpt);
1273  fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt->Fill(corrpt1, dcorrpt/corrpt1);
1274  fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt->Fill(corrpt2, dcorrpt/corrpt2);
1275  fHistDeltaCorrPtvsDistance->Fill(d, dcorrpt);
1276  fHistDeltaCorrPtvsCommonEnergy1->Fill(CE1, dcorrpt);
1277  fHistDeltaCorrPtvsCommonEnergy2->Fill(CE2, dcorrpt);
1278  fHistDeltaCorrPtvsArea1->Fill(jet1->Area(), dcorrpt);
1279  fHistDeltaCorrPtvsArea2->Fill(jet2->Area(), dcorrpt);
1280  fHistDeltaCorrPtvsDeltaArea->Fill(darea, dcorrpt);
1281  fHistJet1CorrPtvsJet2CorrPt->Fill(corrpt1, corrpt2);
1282  }
1283 
1284  if (fIsEmbedded) {
1285  Double_t dmcpt = jet1->MCPt() - jet2->Pt();
1286  fHistDeltaMCPtvsJet1MCPt->Fill(jet1->MCPt(), dmcpt);
1287  fHistDeltaMCPtvsJet2Pt->Fill(jet2->MCPt(), dmcpt);
1288  fHistDeltaMCPtOverJet1MCPtvsJet1MCPt->Fill(jet1->MCPt(), dmcpt/jet1->MCPt());
1289  fHistDeltaMCPtOverJet2PtvsJet2Pt->Fill(jet2->Pt(), dmcpt/jet2->Pt());
1290  fHistDeltaMCPtvsDistance->Fill(d, dmcpt);
1291  fHistDeltaMCPtvsCommonEnergy1->Fill(CE1, dmcpt);
1292  fHistDeltaMCPtvsCommonEnergy2->Fill(CE2, dmcpt);
1293  fHistDeltaMCPtvsArea1->Fill(jet1->Area(), dmcpt);
1294  fHistDeltaMCPtvsArea2->Fill(jet2->Area(), dmcpt);
1295  fHistDeltaMCPtvsDeltaArea->Fill(darea, dmcpt);
1296  fHistJet1MCPtvsJet2Pt->Fill(jet1->MCPt(), jet2->Pt());
1297  }
1298  }
1299 }
1300 
1301 //________________________________________________________________________
1303 {
1304  // Execute once.
1305 
1307 
1308  AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
1309  AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
1310 
1311  if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return;
1312 
1313  if (fMatching == kMCLabel && (!jets2->GetIsParticleLevel() || jets1->GetIsParticleLevel())) {
1314  if (jets1->GetIsParticleLevel() == jets2->GetIsParticleLevel()) {
1315  AliWarning("Changing matching type from MC label to same collection...");
1317  }
1318  else {
1319  AliWarning("Changing matching type from MC label to geometrical...");
1321  }
1322  }
1323  else if (fMatching == kSameCollections && jets1->GetIsParticleLevel() != jets2->GetIsParticleLevel()) {
1324  if (jets2->GetIsParticleLevel() && !jets1->GetIsParticleLevel()) {
1325  AliWarning("Changing matching type from same collection to MC label...");
1326  fMatching = kMCLabel;
1327  }
1328  else {
1329  AliWarning("Changing matching type from same collection to geometrical...");
1331  }
1332  }
1333 }
1334 
1335 //________________________________________________________________________
1337 {
1338  // Find the closest jets
1339 
1340  if (fMatching == kNoMatching)
1341  return kTRUE;
1342  else
1343  return DoJetMatching();
1344 }
1345 
1346 //________________________________________________________________________
1348 {
1349  AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
1350  AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
1351 
1352  if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return kFALSE;
1353 
1354  DoJetLoop();
1355 
1356  AliEmcalJet* jet1 = 0;
1357 
1358  jets1->ResetCurrentID();
1359  while ((jet1 = jets1->GetNextJet())) {
1360 
1361  AliEmcalJet *jet2 = jet1->ClosestJet();
1362 
1363  if (!jet2) continue;
1364  if (jet2->ClosestJet() != jet1) continue;
1365  if (jet1->ClosestJetDistance() > fMatchingPar1 || jet2->ClosestJetDistance() > fMatchingPar2) continue;
1366 
1367  // Matched jet found
1370  AliDebug(2,Form("Found matching: jet1 pt = %f, eta = %f, phi = %f, jet2 pt = %f, eta = %f, phi = %f",
1371  jet1->Pt(), jet1->Eta(), jet1->Phi(),
1372  jet2->Pt(), jet2->Eta(), jet2->Phi()));
1373  }
1374 
1375  return kTRUE;
1376 }
1377 
1378 //________________________________________________________________________
1380 {
1381  // Do the jet loop.
1382 
1383  AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
1384  AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
1385 
1386  if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return;
1387 
1388  AliEmcalJet* jet1 = 0;
1389  AliEmcalJet* jet2 = 0;
1390 
1391  jets2->ResetCurrentID();
1392  while ((jet2 = jets2->GetNextJet())) jet2->ResetMatching();
1393 
1394  jets1->ResetCurrentID();
1395  while ((jet1 = jets1->GetNextJet())) {
1396  jet1->ResetMatching();
1397 
1398  if (jet1->MCPt() < fMinJetMCPt) continue;
1399 
1400  jets2->ResetCurrentID();
1401  while ((jet2 = jets2->GetNextJet())) {
1402  SetMatchingLevel(jet1, jet2, fMatching);
1403  } // jet2 loop
1404  } // jet1 loop
1405 }
1406 
1407 //________________________________________________________________________
1409 {
1410  d = jet1->DeltaR(jet2);
1411 }
1412 
1413 //________________________________________________________________________
1415 {
1416  AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
1417  AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
1418 
1419  if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return;
1420 
1421  // tracks1 just serves as a proxy to ensure that tracks are in jets1
1422  AliParticleContainer *tracks1 = jets1->GetParticleContainer();
1423  // tracks2 is used to retrieve MC labels associated with tracks in the container
1424  // NOTE: For multiple containers, this would need to be generalized!
1425  AliParticleContainer *tracks2 = jets2->GetParticleContainer();
1426 
1427  // d1 and d2 represent the matching level: 0 = maximum level of matching, 1 = the two jets are completely unrelated
1428  d1 = jet1->Pt();
1429  d2 = jet2->Pt();
1430  Double_t totalPt1 = d1; // the total pt of the reconstructed jet will be cleaned from the background
1431 
1432  // remove completely tracks that are not MC particles (label == 0)
1433  if (tracks1 && tracks1->GetArray()) {
1434  for (Int_t iTrack = 0; iTrack < jet1->GetNumberOfTracks(); iTrack++) {
1435  AliVParticle *track = jet1->Track(iTrack);
1436  if (!track) {
1437  AliWarning(Form("Could not find track %d!", iTrack));
1438  continue;
1439  }
1440 
1441  Int_t MClabel = TMath::Abs(track->GetLabel());
1442  MClabel -= fMCLabelShift;
1443  if (MClabel != 0) continue;
1444 
1445  // this is not a MC particle; remove it completely
1446  AliDebug(3,Form("Track %d (pT = %f) is not a MC particle (MClabel = %d)!",iTrack,track->Pt(),MClabel));
1447  totalPt1 -= track->Pt();
1448  d1 -= track->Pt();
1449  }
1450  }
1451 
1452  // remove completely clusters that are not MC particles (label == 0)
1453  if (fUseCellsToMatch && fCaloCells) {
1454  for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
1455  AliVCluster *clus = jet1->Cluster(iClus);
1456  if (!clus) {
1457  AliWarning(Form("Could not find cluster %d!", iClus));
1458  continue;
1459  }
1460  AliTLorentzVector part;
1461  clus->GetMomentum(part, fVertex);
1462 
1463  for (Int_t iCell = 0; iCell < clus->GetNCells(); iCell++) {
1464  Int_t cellId = clus->GetCellAbsId(iCell);
1465  Double_t cellFrac = clus->GetCellAmplitudeFraction(iCell);
1466 
1467  Int_t MClabel = TMath::Abs(fCaloCells->GetCellMCLabel(cellId));
1468  MClabel -= fMCLabelShift;
1469  if (MClabel != 0) continue;
1470 
1471  // this is not a MC particle; remove it completely
1472  AliDebug(3,Form("Cell %d (frac = %f) is not a MC particle (MClabel = %d)!",iCell,cellFrac,MClabel));
1473  totalPt1 -= part.Pt() * cellFrac;
1474  d1 -= part.Pt() * cellFrac;
1475  }
1476  }
1477  }
1478  else {
1479  for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
1480  AliVCluster *clus = jet1->Cluster(iClus);
1481  if (!clus) {
1482  AliWarning(Form("Could not find cluster %d!", iClus));
1483  continue;
1484  }
1485  TLorentzVector part;
1486  clus->GetMomentum(part, fVertex);
1487 
1488  Int_t MClabel = TMath::Abs(clus->GetLabel());
1489  MClabel -= fMCLabelShift;
1490  if (MClabel != 0) continue;
1491 
1492  // this is not a MC particle; remove it completely
1493  AliDebug(3,Form("Cluster %d (pT = %f) is not a MC particle (MClabel = %d)!",iClus,part.Pt(),MClabel));
1494  totalPt1 -= part.Pt();
1495  d1 -= part.Pt();
1496  }
1497  }
1498 
1499  for (Int_t iTrack2 = 0; iTrack2 < jet2->GetNumberOfTracks(); iTrack2++) {
1500  Bool_t track2Found = kFALSE;
1501  Int_t index2 = jet2->TrackAt(iTrack2);
1502 
1503  // now look for common particles in the track array
1504  for (Int_t iTrack = 0; iTrack < jet1->GetNumberOfTracks(); iTrack++) {
1505  AliVParticle *track = jet1->Track(iTrack);
1506  if (!track) {
1507  AliWarning(Form("Could not find track %d!", iTrack));
1508  continue;
1509  }
1510  Int_t MClabel = TMath::Abs(track->GetLabel());
1511  MClabel -= fMCLabelShift;
1512  if (MClabel <= 0) continue;
1513 
1514  Int_t index = -1;
1515  index = tracks2->GetIndexFromLabel(MClabel);
1516  if (index < 0) {
1517  AliDebug(2,Form("Track %d (pT = %f) does not have an associated MC particle (MClabel = %d)!",iTrack,track->Pt(),MClabel));
1518  continue;
1519  }
1520 
1521  if (index2 != index) continue;
1522 
1523  // found common particle
1524  d1 -= track->Pt();
1525 
1526  if (!track2Found) {
1527  AliVParticle *MCpart = jet2->Track(index2);
1528  AliDebug(3,Form("Track %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!",
1529  iTrack,track->Pt(),track->Eta(),track->Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
1530  d2 -= MCpart->Pt();
1531  }
1532 
1533  track2Found = kTRUE;
1534  }
1535 
1536  // now look for common particles in the cluster array
1537  if (fUseCellsToMatch && fCaloCells) { // if the cell colection is available, look for cells with a matched MC particle
1538  for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
1539  AliVCluster *clus = jet1->Cluster(iClus);
1540  if (!clus) {
1541  AliWarning(Form("Could not find cluster %d!", iClus));
1542  continue;
1543  }
1544  AliTLorentzVector part;
1545  clus->GetMomentum(part, fVertex);
1546 
1547  for (Int_t iCell = 0; iCell < clus->GetNCells(); iCell++) {
1548  Int_t cellId = clus->GetCellAbsId(iCell);
1549  Double_t cellFrac = clus->GetCellAmplitudeFraction(iCell);
1550 
1551  Int_t MClabel = TMath::Abs(fCaloCells->GetCellMCLabel(cellId));
1552  MClabel -= fMCLabelShift;
1553  if (MClabel <= 0) continue;
1554 
1555  Int_t index1 = -1;
1556  index1 = tracks2->GetIndexFromLabel(MClabel);
1557  if (index1 < 0) {
1558  AliDebug(3,Form("Cell %d (frac = %f) does not have an associated MC particle (MClabel = %d)!",iCell,cellFrac,MClabel));
1559  continue;
1560  }
1561 
1562  if (index2 != index1) continue;
1563 
1564  // found common particle
1565  d1 -= part.Pt() * cellFrac;
1566 
1567  if (!track2Found) { // only if it is not already found among charged tracks (charged particles are most likely already found)
1568  AliVParticle *MCpart = jet2->Track(index2);
1569  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)!",
1570  iCell,iClus,part.Pt(),part.Eta(),part.Phi_0_2pi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
1571  d2 -= MCpart->Pt() * cellFrac;
1572  }
1573 
1574  track2Found = kTRUE;
1575  }
1576  }
1577  }
1578  else { //otherwise look for the first contributor to the cluster, and if matched to a MC label remove it
1579  for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
1580  AliVCluster *clus = jet1->Cluster(iClus);
1581  if (!clus) {
1582  AliWarning(Form("Could not find cluster %d!", iClus));
1583  continue;
1584  }
1585  AliTLorentzVector part;
1586  clus->GetMomentum(part, fVertex);
1587 
1588  Int_t MClabel = TMath::Abs(clus->GetLabel());
1589  MClabel -= fMCLabelShift;
1590  if (MClabel <= 0) continue;
1591 
1592  Int_t index = -1;
1593  index = tracks2->GetIndexFromLabel(MClabel);
1594 
1595  if (index < 0) {
1596  AliDebug(3,Form("Cluster %d (pT = %f) does not have an associated MC particle (MClabel = %d)!",iClus,part.Pt(),MClabel));
1597  continue;
1598  }
1599 
1600  if (index2 != index) continue;
1601 
1602  // found common particle
1603  d1 -= part.Pt();
1604 
1605  if (!track2Found) { // only if it is not already found among charged tracks (charged particles are most likely already found)
1606  AliVParticle *MCpart = jet2->Track(index2);
1607  AliDebug(3,Form("Cluster %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!",
1608  iClus,part.Pt(),part.Eta(),part.Phi_0_2pi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
1609 
1610  d2 -= MCpart->Pt();
1611  }
1612 
1613  track2Found = kTRUE;
1614  }
1615  }
1616  }
1617 
1618  if (d1 < 0)
1619  d1 = 0;
1620 
1621  if (d2 < 0)
1622  d2 = 0;
1623 
1624  if (totalPt1 < 1)
1625  d1 = -1;
1626  else
1627  d1 /= totalPt1;
1628 
1629  if (jet2->Pt() < 1)
1630  d2 = -1;
1631  else
1632  d2 /= jet2->Pt();
1633 }
1634 
1635 //________________________________________________________________________
1637 {
1638  AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
1639  AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
1640 
1641  if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return;
1642 
1643  // All of the containers are simply used as proxies for whether tracks or clusters are in a jet
1644  AliParticleContainer *tracks1 = jets1->GetParticleContainer();
1645  AliClusterContainer *clusters1 = jets1->GetClusterContainer();
1646  AliParticleContainer *tracks2 = jets2->GetParticleContainer();
1647  AliClusterContainer *clusters2 = jets2->GetClusterContainer();
1648 
1649  // d1 and d2 represent the matching level: 0 = maximum level of matching, 1 = the two jets are completely unrelated
1650  d1 = jet1->Pt();
1651  d2 = jet2->Pt();
1652 
1653  if (tracks1 && tracks2) {
1654 
1655  for (Int_t iTrack2 = 0; iTrack2 < jet2->GetNumberOfTracks(); iTrack2++) {
1656  Int_t index2 = jet2->TrackAt(iTrack2);
1657  for (Int_t iTrack1 = 0; iTrack1 < jet1->GetNumberOfTracks(); iTrack1++) {
1658  Int_t index1 = jet1->TrackAt(iTrack1);
1659  if (index2 == index1) { // found common particle
1660  AliVParticle *part1 = jet1->Track(iTrack1);
1661  if (!part1) {
1662  AliWarning(Form("Could not find track %d!", index1));
1663  continue;
1664  }
1665  AliVParticle *part2 = jet2->Track(iTrack2);
1666  if (!part2) {
1667  AliWarning(Form("Could not find track %d!", index2));
1668  continue;
1669  }
1670 
1671  d1 -= part1->Pt();
1672  d2 -= part2->Pt();
1673  break;
1674  }
1675  }
1676  }
1677 
1678  }
1679 
1680  if (clusters1 && clusters2) {
1681 
1682  if (fUseCellsToMatch && fCaloCells) {
1683  // Note: this section of the code needs to be revised and tested heavily
1684  // While fixing some inconsistencies in AliAnalysisTaskEMCALClusterizeFast
1685  // some issues came up, e.g. memory leaks (fixed) and inconsistent use of
1686  // fCaloCells. In principle the two cluster collections may use cells
1687  // from different sources (embedding / non-embedding). This is not handled
1688  // correctly in the current version of this code.
1689  AliWarning("ATTENTION ATTENTION ATTENTION: this section of the AliJetResponseMaker code needs to be revised and tested before using it for physics!!!");
1690  const Int_t nClus1 = jet1->GetNumberOfClusters();
1691 
1692  Int_t ncells1[nClus1];
1693  UShort_t *cellsId1[nClus1];
1694  Double_t *cellsFrac1[nClus1];
1695  Double_t *cellsClusFrac1[nClus1];
1696  Int_t *sortedIndexes1[nClus1];
1697  Double_t ptClus1[nClus1];
1698  for (Int_t iClus1 = 0; iClus1 < nClus1; iClus1++) {
1699  Int_t index1 = jet1->ClusterAt(iClus1);
1700  AliVCluster *clus1 = clusters1->GetCluster(index1);
1701  if (!clus1) {
1702  AliWarning(Form("Could not find cluster %d!", index1));
1703  ncells1[iClus1] = 0;
1704  cellsId1[iClus1] = 0;
1705  cellsFrac1[iClus1] = 0;
1706  cellsClusFrac1[iClus1] = 0;
1707  sortedIndexes1[iClus1] = 0;
1708  ptClus1[iClus1] = 0;
1709  continue;
1710  }
1711  TLorentzVector part1;
1712  clus1->GetMomentum(part1, fVertex);
1713 
1714  ncells1[iClus1] = clus1->GetNCells();
1715  cellsId1[iClus1] = clus1->GetCellsAbsId();
1716  cellsFrac1[iClus1] = clus1->GetCellsAmplitudeFraction();
1717  cellsClusFrac1[iClus1] = new Double_t[ncells1[iClus1]];
1718  sortedIndexes1[iClus1] = new Int_t[ncells1[iClus1]];
1719  ptClus1[iClus1] = part1.Pt();
1720 
1721  for (Int_t iCell = 0; iCell < ncells1[iClus1]; iCell++) {
1722  cellsClusFrac1[iClus1][iCell] = fCaloCells->GetCellAmplitude(cellsId1[iClus1][iCell]) / clus1->E();
1723  }
1724 
1725  TMath::Sort(ncells1[iClus1], cellsId1[iClus1], sortedIndexes1[iClus1], kFALSE);
1726  }
1727 
1728  const Int_t nClus2 = jet2->GetNumberOfClusters();
1729 
1730  const Int_t maxNcells2 = 11520;
1731  Int_t sortedIndexes2[maxNcells2];
1732  for (Int_t iClus2 = 0; iClus2 < nClus2; iClus2++) {
1733  Int_t index2 = jet2->ClusterAt(iClus2);
1734  AliVCluster *clus2 = clusters2->GetCluster(index2);
1735  if (!clus2) {
1736  AliWarning(Form("Could not find cluster %d!", index2));
1737  continue;
1738  }
1739  Int_t ncells2 = clus2->GetNCells();
1740  if (ncells2 >= maxNcells2) {
1741  AliError(Form("Number of cells in the cluster %d >= %d",ncells2,maxNcells2));
1742  continue;
1743  }
1744  UShort_t *cellsId2 = clus2->GetCellsAbsId();
1745  Double_t *cellsFrac2 = clus2->GetCellsAmplitudeFraction();
1746  Double_t *cellsClusFrac2 = new Double_t[ncells2];
1747 
1748  for (Int_t iCell = 0; iCell < ncells2; iCell++) {
1749  cellsClusFrac2[iCell] = fCaloCells->GetCellAmplitude(cellsId2[iCell]) / clus2->E();
1750  }
1751 
1752  TLorentzVector part2;
1753  clus2->GetMomentum(part2, fVertex);
1754  Double_t ptClus2 = part2.Pt();
1755 
1756  TMath::Sort(ncells2, cellsId2, sortedIndexes2, kFALSE);
1757 
1758  for (Int_t iClus1 = 0; iClus1 < nClus1; iClus1++) {
1759  if (sortedIndexes1[iClus1] == 0)
1760  continue;
1761  Int_t iCell1 = 0, iCell2 = 0;
1762  while (iCell1 < ncells1[iClus1] && iCell2 < ncells2) {
1763  if (cellsId1[iClus1][sortedIndexes1[iClus1][iCell1]] == cellsId2[sortedIndexes2[iCell2]]) { // found a common cell
1764  d1 -= cellsFrac1[iClus1][sortedIndexes1[iClus1][iCell1]] * cellsClusFrac1[iClus1][sortedIndexes1[iClus1][iCell1]] * ptClus1[iClus1];
1765  d2 -= cellsFrac2[sortedIndexes2[iCell2]] * cellsClusFrac2[sortedIndexes2[iCell2]] * ptClus2;
1766  iCell1++;
1767  iCell2++;
1768  }
1769  else if (cellsId1[iClus1][sortedIndexes1[iClus1][iCell1]] > cellsId2[sortedIndexes2[iCell2]]) {
1770  iCell2++;
1771  }
1772  else {
1773  iCell1++;
1774  }
1775  }
1776  }
1777  delete[] cellsClusFrac2;
1778  }
1779  for (Int_t iClus1 = 0; iClus1 < nClus1; iClus1++) {
1780  delete[] cellsClusFrac1[iClus1];
1781  delete[] sortedIndexes1[iClus1];
1782  }
1783  }
1784  else {
1785  for (Int_t iClus2 = 0; iClus2 < jet2->GetNumberOfClusters(); iClus2++) {
1786  Int_t index2 = jet2->ClusterAt(iClus2);
1787  for (Int_t iClus1 = 0; iClus1 < jet1->GetNumberOfClusters(); iClus1++) {
1788  Int_t index1 = jet1->ClusterAt(iClus1);
1789  if (index2 == index1) { // found common particle
1790  AliVCluster *clus1 = jet1->Cluster(iClus1);
1791  if (!clus1) {
1792  AliWarning(Form("Could not find cluster %d!", index1));
1793  continue;
1794  }
1795  AliVCluster *clus2 = jet2->Cluster(iClus2);
1796  if (!clus2) {
1797  AliWarning(Form("Could not find cluster %d!", index2));
1798  continue;
1799  }
1800  TLorentzVector part1, part2;
1801  clus1->GetMomentum(part1, fVertex);
1802  clus2->GetMomentum(part2, fVertex);
1803 
1804  d1 -= part1.Pt();
1805  d2 -= part2.Pt();
1806  break;
1807  }
1808  }
1809  }
1810  }
1811  }
1812 
1813  if (d1 < 0)
1814  d1 = 0;
1815 
1816  if (d2 < 0)
1817  d2 = 0;
1818 
1819  if (jet1->Pt() > 0)
1820  d1 /= jet1->Pt();
1821  else
1822  d1 = -1;
1823 
1824  if (jet2->Pt() > 0)
1825  d2 /= jet2->Pt();
1826  else
1827  d2 = -1;
1828 }
1829 
1830 //________________________________________________________________________
1832 {
1833  Double_t d1 = -1;
1834  Double_t d2 = -1;
1835 
1836  switch (matching) {
1837  case kGeometrical:
1838  GetGeometricalMatchingLevel(jet1,jet2,d1);
1839  d2 = d1;
1840  break;
1841  case kMCLabel: // jet1 = detector level and jet2 = particle level!
1842  GetMCLabelMatchingLevel(jet1,jet2,d1,d2);
1843  break;
1844  case kSameCollections:
1845  GetSameCollectionsMatchingLevel(jet1,jet2,d1,d2);
1846  break;
1847  default:
1848  ;
1849  }
1850 
1851  if (d1 >= 0) {
1852 
1853  if (d1 < jet1->ClosestJetDistance()) {
1854  jet1->SetSecondClosestJet(jet1->ClosestJet(), jet1->ClosestJetDistance());
1855  jet1->SetClosestJet(jet2, d1);
1856  }
1857  else if (d1 < jet1->SecondClosestJetDistance()) {
1858  jet1->SetSecondClosestJet(jet2, d1);
1859  }
1860  }
1861 
1862  if (d2 >= 0) {
1863 
1864  if (d2 < jet2->ClosestJetDistance()) {
1865  jet2->SetSecondClosestJet(jet2->ClosestJet(), jet2->ClosestJetDistance());
1866  jet2->SetClosestJet(jet1, d2);
1867  }
1868  else if (d2 < jet2->SecondClosestJetDistance()) {
1869  jet2->SetSecondClosestJet(jet1, d2);
1870  }
1871  }
1872 }
1873 
1874 //________________________________________________________________________
1876 {
1877  // Fill histograms.
1878 
1879  AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
1880  AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
1881 
1882  if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return kFALSE;
1883 
1884  AliEmcalJet* jet1 = 0;
1885  AliEmcalJet* jet2 = 0;
1886 
1887  jets2->ResetCurrentID();
1888  while ((jet2 = jets2->GetNextJet())) {
1889 
1890  AliDebug(2,Form("Processing jet (2) %d", jets2->GetCurrentID()));
1891 
1892  if (jet2->Pt() < jets2->GetJetPtCut()) continue;
1893 
1894  UInt_t rejectionReason = 0;
1895  if (jets2->AcceptJet(jet2, rejectionReason))
1896  FillJetHisto(jet2, 2);
1897  else
1898  fHistRejectionReason2->Fill(jets2->GetRejectionReasonBitPosition(rejectionReason), jet2->Pt());
1899 
1900  jet1 = jet2->MatchedJet();
1901 
1902  if (!jet1) continue;
1903  rejectionReason = 0;
1904  if (!jets1->AcceptJet(jet1, rejectionReason)) continue;
1905  if (jet1->MCPt() < fMinJetMCPt) continue;
1906 
1907  Double_t d=-1, ce1=-1, ce2=-1;
1908  if (jet2->GetMatchingType() == kGeometrical) {
1909  if (jets2->GetIsParticleLevel() && !jets1->GetIsParticleLevel()) // the other way around is not supported
1910  GetMCLabelMatchingLevel(jet1, jet2, ce1, ce2);
1911  else if (jets1->GetIsParticleLevel() == jets2->GetIsParticleLevel())
1912  GetSameCollectionsMatchingLevel(jet1, jet2, ce1, ce2);
1913 
1914  d = jet2->ClosestJetDistance();
1915  }
1916  else if (jet2->GetMatchingType() == kMCLabel || jet2->GetMatchingType() == kSameCollections) {
1917  GetGeometricalMatchingLevel(jet1, jet2, d);
1918 
1919  ce1 = jet1->ClosestJetDistance();
1920  ce2 = jet2->ClosestJetDistance();
1921  }
1922 
1923  FillMatchingHistos(jet1, jet2, d, ce1, ce2);
1924  }
1925 
1926  jets1->ResetCurrentID();
1927  while ((jet1 = jets1->GetNextJet())) {
1928  UInt_t rejectionReason = 0;
1929  if (!jets1->AcceptJet(jet1, rejectionReason)) {
1930  fHistRejectionReason1->Fill(jets1->GetRejectionReasonBitPosition(rejectionReason), jet1->Pt());
1931  continue;
1932  }
1933  if (jet1->MCPt() < fMinJetMCPt) continue;
1934  AliDebug(2,Form("Processing jet (1) %d", jets1->GetCurrentID()));
1935 
1936  FillJetHisto(jet1, 1);
1937  }
1938  return kTRUE;
1939 }
1940 
1951 {
1952  Double_t dphi = (epAngle - jetAngle);
1953 
1954  // ran into trouble with a few dEP<-Pi so trying this...
1955  if( dphi<-1*TMath::Pi() ){
1956  dphi = dphi + 1*TMath::Pi();
1957  } // this assumes we are doing full jets currently
1958 
1959  if( (dphi>0) && (dphi<1*TMath::Pi()/2) ){
1960  // Do nothing! we are in quadrant 1
1961  }else if( (dphi>1*TMath::Pi()/2) && (dphi<1*TMath::Pi()) ){
1962  dphi = 1*TMath::Pi() - dphi;
1963  }else if( (dphi<0) && (dphi>-1*TMath::Pi()/2) ){
1964  dphi = fabs(dphi);
1965  }else if( (dphi<-1*TMath::Pi()/2) && (dphi>-1*TMath::Pi()) ){
1966  dphi = dphi + 1*TMath::Pi();
1967  }
1968 
1969  // test
1970  if( dphi < 0 || dphi > TMath::Pi()/2 )
1971  AliWarning(Form("%s: dPHI not in range [0, 0.5*Pi]!", GetName()));
1972 
1973  return dphi; // dphi in [0, Pi/2]
1974 }
void SetSecondClosestJet(AliEmcalJet *j, Double_t d)
Definition: AliEmcalJet.h:212
TH2 * fHistJets2PtArea
phi-eta distribution of jets 2
Double_t Area() const
Definition: AliEmcalJet.h:117
TH2 * fHistRejectionReason1
whether the jet2 collection has to be average subtracted
TH2 * fHistJets2CorrPtArea
inclusive jet pt vs. area histogram 2
Double_t GetRhoVal() const
const TString & GetRhoName() const
double Double_t
Definition: External.C:58
TH2 * fHistDeltaPtvsCommonEnergy1
delta pt between matched jets vs distance
TH2 * fHistJet2PtOverJet1PtvsJet2Pt
delta eta vs delta phi of matched jets
Definition: External.C:236
const char * title
Definition: MakeQAPdf.C:26
Double_t MCPt() const
Definition: AliEmcalJet.h:140
TH2 * fHistDeltaPtvsDeltaArea
delta pt between matched jets vs jet 2 area
AliEmcalJet * ClosestJet() const
Definition: AliEmcalJet.h:215
AliJetContainer * GetJetContainer(Int_t i=0) const
Double_t ClosestJetDistance() const
Definition: AliEmcalJet.h:216
Double_t Eta() const
Definition: AliEmcalJet.h:108
TH2 * fHistJets2PhiEta
Constituent Pt over Jet Pt ratio vs. jet pt 1.
TH2 * fHistDeltaCorrPtvsArea2
delta pt corr between matched jets vs jet 1 area
void FillJetHisto(AliEmcalJet *jet, Int_t Set)
Double_t Phi() const
Definition: AliEmcalJet.h:104
Int_t fJetRelativeEPAngle
add jet angle relative to the EP in matching THnSparse (default=0)
AliEmcalJet * MatchedJet() const
Definition: AliEmcalJet.h:219
Declaration of class AliTLorentzVector.
void GetSameCollectionsMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d1, Double_t &d2) const
Double_t fMinBinPt
min pt in histograms
Double_t fEPV0
!event plane V0
Int_t ClusterAt(Int_t idx) const
Definition: AliEmcalJet.h:124
TH2 * fHistDeltaMCPtvsCommonEnergy2
jet 1 MC pt - jet2 pt vs common energy 1 (%)
AliClusterContainer * GetClusterContainer() const
TH2 * fHistJets1ZvsPt
Jet charged energy fraction vs. charged jet pt 1.
AliVParticle * Track(Int_t idx) const
virtual Bool_t AcceptJet(Int_t i, UInt_t &rejectionReason) const
Double_t GetRelativeEPAngle(Double_t jetAngle, Double_t epAngle) const
THnSparse * fHistMatching
jet2 THnSparse
TH2 * fHistDeltaEtaDeltaPhi
common energy 1 (%) vs common energy 2 (%)
ClassImp(AliJetResponseMaker) AliJetResponseMaker
TH2 * fHistDeltaMCPtvsCommonEnergy1
jet 1 MC pt - jet2 pt vs distance
TH2 * fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt
delta pt corr / jet 1 corr pt between matched jets vs jet 1 corr pt
TH2 * fHistCommonEnergy1vsCommonEnergy2
distance vs common energy 2 (%)
TH2 * fHistDeltaPtvsCommonEnergy2
delta pt between matched jets vs common energy 1 (%)
TH2 * fHistDeltaMCPtvsDistance
jet 1 MC pt - jet2 pt vs jet 2 pt
AliVParticle * GetFlavourTrack(Int_t i=0) const
Container for particles within the EMCAL framework.
TH2 * fHistJets1NEFvsPt
inclusive jet pt vs. area histogram 1
Int_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:147
TH2 * fHistJets2ZvsPt
Jet charged energy fraction vs. charged jet pt 2.
Bool_t fIsEmbedded
trigger, embedded signal
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:126
Bool_t fIsJet2Rho
whether the jet1 collection has to be average subtracted
void ResetMatching()
TH2 * fHistDeltaPtvsJet1Pt
jet 1 pt over jet 2 pt vs jet 1 pt
TH2 * fHistDeltaCorrPtvsCommonEnergy2
delta pt corr between matched jets vs common energy 1 (%)
AliParticleContainer * GetParticleContainer() const
void GetLeadingHadronMomentum(TLorentzVector &mom, const AliEmcalJet *jet) const
TH2 * fHistDeltaPtvsJet2Pt
delta pt between matched jets vs jet 1 pt
TH2 * fHistJets1PhiEta
matching THnSparse
int Int_t
Definition: External.C:63
UShort_t GetNumberOfClusters() const
Definition: AliEmcalJet.h:125
unsigned int UInt_t
Definition: External.C:33
THnSparse * fHistJets1
Rejection reason vs. jet pt.
UShort_t GetMatchingType() const
Definition: AliEmcalJet.h:220
TH2 * fHistJet1PtOverJet2PtvsJet1Pt
jet 2 pt over jet 1 pt vs jet 2 pt
Double_t Phi_0_2pi() const
TH2 * fHistDeltaCorrPtvsDistance
delta pt corr between matched jets vs jet 2 corr pt
void FillMatchingHistos(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t d, Double_t CE1, Double_t CE2)
TH2 * fHistJet1PtvsJet2Pt
delta pt between matched jets vs delta area
TH2 * fHistDeltaCorrPtvsJet1CorrPt
delta pt corr / jet 2 corr pt between matched jets vs jet 2 corr pt
TH2 * fHistDeltaMCPtOverJet2PtvsJet2Pt
jet 1 MC pt - jet2 pt / jet 1 MC pt vs jet 1 pt
TH2 * fHistDistancevsCommonEnergy2
distance vs common energy 1 (%)
Double_t GetJetPtCut() const
TH2 * fHistRejectionReason2
Rejection reason vs. jet pt.
Int_t fMCLabelShift
if MC label > fMCLabelShift, MC label -= fMCLabelShift
TH2 * fHistDistancevsJet1Pt
common energy 2 (%) vs jet 2 pt
TH2 * fHistDeltaMCPtOverJet1MCPtvsJet1MCPt
correlation jet 1 corr pt vs jet 2 corr pt
static Double_t GetParallelFraction(AliVParticle *part1, AliVParticle *part2)
AliVCluster * GetCluster(Int_t i) const
TH2 * fHistDeltaPtOverJet1PtvsJet1Pt
delta pt between matched jets vs jet 2 pt
TObjArray fJetCollArray
jet collection array
Double_t DeltaR(const AliVParticle *part) const
TH2 * fHistDeltaCorrPtvsArea1
delta pt corr between matched jets vs common energy 2 (%)
AliVCaloCells * fCaloCells
!cells
TH2 * fHistCommonEnergy1vsJet1Pt
Constituent Pt over Jet Pt ratio vs. jet pt 2.
void SetMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, MatchingType matching)
Double_t Pt() const
Definition: AliEmcalJet.h:96
TH2 * fHistJets1CorrPtArea
inclusive jet pt vs. area histogram 1
AliEmcalJet * GetNextJet()
AliEmcalList * fOutput
!output list
Double_t fMaxBinPt
max pt in histograms
TH2 * fHistDeltaMCPtvsJet1MCPt
jet 1 MC pt - jet2 pt / jet 2 pt vs jet 2 pt
void SetClosestJet(AliEmcalJet *j, Double_t d)
Definition: AliEmcalJet.h:211
TH2 * fHistDeltaMCPtvsArea2
jet 1 MC pt - jet2 pt vs jet 1 area
TH2 * fHistDeltaMCPtvsJet2Pt
jet 1 MC pt - jet2 pt vs jet 1 MC pt
TH2 * fHistDistancevsJet2Pt
distance vs jet 1 pt
TH2 * fHistDeltaPtvsArea1
delta pt between matched jets vs common energy 2 (%)
Double_t fVertex[3]
!event vertex
TH2 * fHistDeltaCorrPtvsJet2CorrPt
delta pt corr between matched jets vs jet 1 corr pt
void GetGeometricalMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d) const
THnSparse * fHistJets2
jet1 THnSparse
TH2 * fHistDeltaPtOverJet2PtvsJet2Pt
delta pt / jet 1 pt between matched jets vs jet 1 pt
AliVCluster * Cluster(Int_t idx) const
TH2 * fHistDeltaPtvsArea2
delta pt between matched jets vs jet 1 area
void SetMakeGeneralHistograms(Bool_t g)
TH2 * fHistJets2NEFvsPt
inclusive jet pt vs. area histogram 2
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:44
TH2 * fHistCommonEnergy2vsJet2Pt
common energy 1 (%) vs jet 1 pt
unsigned short UShort_t
Definition: External.C:28
TH2 * fHistDeltaPtvsDistance
delta pt / jet 2 pt between matched jets vs jet 2 pt
TH2 * fHistJet1MCPtvsJet2Pt
jet 1 MC pt - jet2 pt vs delta area
TH2 * fHistJets2CEFvsCEFPt
Jet neutral energy fraction vs. jet pt 2.
void SetRejectionReasonLabels(TAxis *axis)
TH2 * fHistDistancevsCommonEnergy1
distance vs jet 2 pt
TH2 * fHistJets1CEFvsCEFPt
Jet neutral energy fraction vs. jet pt 1.
const Int_t nbins
bool Bool_t
Definition: External.C:53
TH2 * fHistDeltaMCPtvsArea1
jet 1 MC pt - jet2 pt vs common energy 2 (%)
AliEmcalJetShapeProperties * GetShapeProperties() const
Definition: AliEmcalJet.h:249
TH2 * fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt
correlation jet 1 pt vs jet 2 pt
TH2 * fHistJets1PtArea
phi-eta distribution of jets 1
Double_t NEF() const
Definition: AliEmcalJet.h:135
TH2 * fHistDeltaMCPtvsDeltaArea
jet 1 MC pt - jet2 pt vs jet 2 area
Container structure for EMCAL clusters.
TH2 * fHistDeltaCorrPtvsCommonEnergy1
delta pt corr between matched jets vs distance
void GetMCLabelMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d1, Double_t &d2) const
Container for jet within the EMCAL jet framework.
void SetMatchedToClosest(UShort_t m)
Definition: AliEmcalJet.h:213
Int_t fNbins
no. of pt bins
TH2 * fHistJet1CorrPtvsJet2CorrPt
delta pt corr between matched jets vs delta area
TH2 * fHistDeltaCorrPtvsDeltaArea
delta pt corr between matched jets vs jet 2 area