AliPhysics  97344c9 (97344c9)
 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  fDBCAxis(0),
57  fFlavourZAxis(0),
58  fFlavourPtAxis(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  fDBCAxis(0),
144  fFlavourZAxis(0),
145  fFlavourPtAxis(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  title[dim] = "p_{T,particle}^{leading} (GeV/c)";
687  nbins[dim] = 120;
688  min[dim] = 0;
689  max[dim] = 120;
690  dim++;
691 
692  Int_t dim1 = dim, dim2 = dim;
693 
694  if (fIsJet1Rho) {
695  title[dim1] = "p_{T}^{corr}";
696  nbins[dim1] = fNbins*2;
697  min[dim1] = -250;
698  max[dim1] = 250;
699  dim1++;
700  }
701 
702  if (fIsEmbedded) {
703  title[dim1] = "p_{T}^{MC}";
704  nbins[dim1] = fNbins;
705  min[dim1] = 0;
706  max[dim1] = 250;
707  dim1++;
708  }
709 
710  fHistJets1 = new THnSparseD("fHistJets1","fHistJets1",dim1,nbins,min,max);
711  for (Int_t i = 0; i < dim1; i++)
712  fHistJets1->GetAxis(i)->SetTitle(title[i]);
713  fOutput->Add(fHistJets1);
714 
715  if (fIsJet2Rho) {
716  title[dim2] = "p_{T}^{corr}";
717  nbins[dim2] = fNbins*2;
718  min[dim2] = -250;
719  max[dim2] = 250;
720  dim2++;
721  }
722 
723  fHistJets2 = new THnSparseD("fHistJets2","fHistJets2",dim2,nbins,min,max);
724  for (Int_t i = 0; i < dim2; i++)
725  fHistJets2->GetAxis(i)->SetTitle(title[i]);
726  fOutput->Add(fHistJets2);
727 
728  // Matching
729 
730  dim = 0;
731 
732  title[dim] = "p_{T,1}";
733  nbins[dim] = fNbins;
734  min[dim] = 0;
735  max[dim] = 250;
736  dim++;
737 
738  title[dim] = "p_{T,2}";
739  nbins[dim] = fNbins;
740  min[dim] = 0;
741  max[dim] = 250;
742  dim++;
743 
744  title[dim] = "A_{jet,1}";
745  nbins[dim] = fNbins/4;
746  min[dim] = 0;
747  max[dim] = 1.5;
748  dim++;
749 
750  title[dim] = "A_{jet,2}";
751  nbins[dim] = fNbins/4;
752  min[dim] = 0;
753  max[dim] = 1.5;
754  dim++;
755 
756  title[dim] = "distance";
757  nbins[dim] = fNbins/2;
758  min[dim] = 0;
759  max[dim] = 1.2;
760  dim++;
761 
762  title[dim] = "CE1";
763  nbins[dim] = fNbins/2;
764  min[dim] = 0;
765  max[dim] = 1.2;
766  dim++;
767 
768  title[dim] = "CE2";
769  nbins[dim] = fNbins/2;
770  min[dim] = 0;
771  max[dim] = 1.2;
772  dim++;
773 
774  title[dim] = "p_{T,particle,1}^{leading} (GeV/c)";
775  nbins[dim] = 120;
776  min[dim] = 0;
777  max[dim] = 120;
778  dim++;
779 
780  title[dim] = "p_{T,particle,2}^{leading} (GeV/c)";
781  nbins[dim] = 120;
782  min[dim] = 0;
783  max[dim] = 120;
784  dim++;
785 
786  if (fDeltaPtAxis) {
787  title[dim] = "#deltaA_{jet}";
788  nbins[dim] = fNbins/2;
789  min[dim] = -1;
790  max[dim] = 1;
791  dim++;
792 
793  title[dim] = "#deltap_{T}";
794  nbins[dim] = fNbins*2;
795  min[dim] = -250;
796  max[dim] = 250;
797  dim++;
798  }
799  if (fIsJet1Rho) {
800  title[dim] = "p_{T,1}^{corr}";
801  nbins[dim] = fNbins*2;
802  min[dim] = -250;
803  max[dim] = 250;
804  dim++;
805  }
806  if (fIsJet2Rho) {
807  title[dim] = "p_{T,2}^{corr}";
808  nbins[dim] = fNbins*2;
809  min[dim] = -250;
810  max[dim] = 250;
811  dim++;
812  }
813  if (fDeltaPtAxis && (fIsJet1Rho || fIsJet2Rho)) {
814  title[dim] = "#deltap_{T}^{corr}";
815  nbins[dim] = fNbins*2;
816  min[dim] = -250;
817  max[dim] = 250;
818  dim++;
819  }
820  if (fDeltaEtaDeltaPhiAxis) {
821  title[dim] = "#delta#eta";
822  nbins[dim] = fNbins/2;
823  min[dim] = -1;
824  max[dim] = 1;
825  dim++;
826 
827  title[dim] = "#delta#phi";
828  nbins[dim] = fNbins/2;
829  min[dim] = -TMath::Pi()/2;
830  max[dim] = TMath::Pi()*3/2;
831  dim++;
832  }
833  if (fIsEmbedded) {
834  title[dim] = "p_{T,1}^{MC}";
835  nbins[dim] = fNbins;
836  min[dim] = 0;
837  max[dim] = 250;
838  dim++;
839 
840  if (fDeltaPtAxis) {
841  title[dim] = "#deltap_{T}^{MC}";
842  nbins[dim] = fNbins*2;
843  min[dim] = -250;
844  max[dim] = 250;
845  dim++;
846  }
847  }
848 
849  if (fNEFAxis) {
850  title[dim] = "NEF_{1}";
851  nbins[dim] = 120;
852  min[dim] = 0.0;
853  max[dim] = 1.2;
854  dim++;
855 
856  title[dim] = "NEF_{2}";
857  nbins[dim] = 120;
858  min[dim] = 0.0;
859  max[dim] = 1.2;
860  dim++;
861  }
862 
863  if (fZAxis) {
864  title[dim] = "Z_{1}";
865  nbins[dim] = 120;
866  min[dim] = 0.0;
867  max[dim] = 1.2;
868  dim++;
869 
870  title[dim] = "Z_{2}";
871  nbins[dim] = 120;
872  min[dim] = 0.0;
873  max[dim] = 1.2;
874  dim++;
875  }
876 
877  if (fFlavourZAxis) {
878  title[dim] = "z_{flavour,1}";
879  nbins[dim] = 100;
880  min[dim] = 0.0;
881  max[dim] = 2.0;
882  dim++;
883 
884  title[dim] = "z_{flavour,2}";
885  nbins[dim] = 100;
886  min[dim] = 0.0;
887  max[dim] = 2.0;
888  dim++;
889  }
890 
891  if (fFlavourPtAxis) {
892  title[dim] = "p_{T,1}^{D}";
893  nbins[dim] = fNbins/2;
894  min[dim] = 0;
895  max[dim] = 125;
896  dim++;
897 
898  title[dim] = "p_{T,2}^{D}";
899  nbins[dim] = fNbins/2;
900  min[dim] = 0;
901  max[dim] = 125;
902  dim++;
903  }
904 
905  if (fZgAxis) {
906  title[dim] = "Z_{g,1}";
907  nbins[dim] = 20;
908  min[dim] = 0.0;
909  max[dim] = 1.0;
910  dim++;
911  title[dim] = "Z_{g,2}";
912  nbins[dim] = 20;
913  min[dim] = 0.0;
914  max[dim] = 1.0;
915  dim++;
916  }
917 
918  if (fdRAxis) {
919  title[dim] = "dR_{1}";
920  nbins[dim] = 40;
921  min[dim] = 0.0;
922  max[dim] = 0.5;
923  dim++;
924  title[dim] = "dR_{2}";
925  nbins[dim] = 40;
926  min[dim] = 0.0;
927  max[dim] = 0.5;
928  dim++;
929  }
930 
931  if (fPtgAxis) {
932  title[dim] = "p_{T,g,1}";
933  nbins[dim] = 16;
934  min[dim] = 0.0;
935  max[dim] = 160.0;
936  dim++;
937  title[dim] = "p_{T,g,2}";
938  nbins[dim] = 16;
939  min[dim] = 0.0;
940  max[dim] = 160.0;
941  dim++;
942  }
943 
944  if (fDBCAxis) {
945  title[dim] = "DBC_{1}";
946  nbins[dim] = 20;
947  min[dim] = 0.0;
948  max[dim] = 20.0;
949  dim++;
950  title[dim] = "DBC_{2}";
951  nbins[dim] = 20;
952  min[dim] = 0.0;
953  max[dim] = 20.0;
954  dim++;
955  }
956 
957  fHistMatching = new THnSparseD("fHistMatching","fHistMatching",dim,nbins,min,max);
958 
959  for (Int_t i = 0; i < dim; i++)
960  fHistMatching->GetAxis(i)->SetTitle(title[i]);
961 
962  fOutput->Add(fHistMatching);
963 }
964 
965 //________________________________________________________________________
967 {
968  // Create user objects.
969 
971 
972  AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
973  AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
974 
975  if (!jets1 || !jets2) return;
976 
977  if (jets1->GetRhoName().IsNull()) fIsJet1Rho = kFALSE;
978  else fIsJet1Rho = kTRUE;
979  if (jets2->GetRhoName().IsNull()) fIsJet2Rho = kFALSE;
980  else fIsJet2Rho = kTRUE;
981 
982  fHistRejectionReason1 = new TH2F("fHistRejectionReason1", "fHistRejectionReason1", 32, 0, 32, 100, 0, 250);
983  fHistRejectionReason1->GetXaxis()->SetTitle("Rejection reason");
984  fHistRejectionReason1->GetYaxis()->SetTitle("p_{T,jet} (GeV/c)");
985  fHistRejectionReason1->GetZaxis()->SetTitle("counts");
988 
989  fHistRejectionReason2 = new TH2F("fHistRejectionReason2", "fHistRejectionReason2", 32, 0, 32, 100, 0, 250);
990  fHistRejectionReason2->GetXaxis()->SetTitle("Rejection reason");
991  fHistRejectionReason2->GetYaxis()->SetTitle("p_{T,jet} (GeV/c)");
992  fHistRejectionReason2->GetZaxis()->SetTitle("counts");
995 
996  if (fHistoType==0)
997  AllocateTH2();
998  else
1000 
1001  PostData(1, fOutput); // Post data for ALL output slots > 0 here, to get at least an empty histogram
1002 }
1003 
1004 //________________________________________________________________________
1006 {
1007  AliJetContainer* jets = GetJetContainer(Set-1);
1008 
1009  AliTLorentzVector leadPart;
1010  jets->GetLeadingHadronMomentum(leadPart, jet);
1011  Double_t zleading = GetParallelFraction(leadPart.Vect(), jet);
1012  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
1013 
1014  Double_t corrpt = jet->Pt() - jets->GetRhoVal() * jet->Area();
1015  Double_t zflavour = 0;
1016  Double_t ptflavour = 0;
1017  AliVParticle* hftrack = jet->GetFlavourTrack();
1018  if (hftrack) {
1019  zflavour = GetParallelFraction(hftrack, jet);
1020  ptflavour = hftrack->Pt();
1021 
1022  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
1023  }
1024 
1025  if (fHistoType==1) {
1026  THnSparse *histo = 0;
1027  if (Set==1) {
1028  histo = fHistJets1;
1029  }
1030  else if (Set==2) {
1031  histo = fHistJets2;
1032  }
1033 
1034  if (!histo) return;
1035 
1036  Double_t contents[20]={0};
1037 
1038  for (Int_t i = 0; i < histo->GetNdimensions(); i++) {
1039  TString title(histo->GetAxis(i)->GetTitle());
1040  if (title=="#phi")
1041  contents[i] = jet->Phi();
1042  else if (title=="#eta")
1043  contents[i] = jet->Eta();
1044  else if (title=="p_{T}")
1045  contents[i] = jet->Pt();
1046  else if (title=="A_{jet}")
1047  contents[i] = jet->Area();
1048  else if (title=="NEF")
1049  contents[i] = jet->NEF();
1050  else if (title=="Z")
1051  contents[i] = zleading;
1052  else if (title=="p_{T}^{corr}")
1053  contents[i] = corrpt;
1054  else if (title=="p_{T}^{MC}")
1055  contents[i] = jet->MCPt();
1056  else if (title=="p_{T,particle}^{leading} (GeV/c)")
1057  contents[i] = leadPart.Pt();
1058  else if (title=="z_{flavour}")
1059  contents[i] = zflavour;
1060  else if (title=="p_{T}^{D}")
1061  contents[i] = ptflavour;
1062  else
1063  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1064  }
1065 
1066  histo->Fill(contents);
1067  }
1068  else {
1069  if (Set == 1) {
1070  fHistJets1PtArea->Fill(jet->Area(), jet->Pt());
1071  fHistJets1PhiEta->Fill(jet->Eta(), jet->Phi());
1072  fHistJets1ZvsPt->Fill(zleading, jet->Pt());
1073  fHistJets1NEFvsPt->Fill(jet->NEF(), jet->Pt());
1074  fHistJets1CEFvsCEFPt->Fill(1-jet->NEF(), (1-jet->NEF())*jet->Pt());
1075  if (fIsJet1Rho)
1076  fHistJets1CorrPtArea->Fill(jet->Area(), corrpt);
1077  }
1078  else if (Set == 2) {
1079  fHistJets2PtArea->Fill(jet->Area(), jet->Pt());
1080  fHistJets2PhiEta->Fill(jet->Eta(), jet->Phi());
1081  fHistJets2ZvsPt->Fill(zleading, jet->Pt());
1082  fHistJets2NEFvsPt->Fill(jet->NEF(), jet->Pt());
1083  fHistJets2CEFvsCEFPt->Fill(1-jet->NEF(), (1-jet->NEF())*jet->Pt());
1084  if (fIsJet2Rho)
1085  fHistJets2CorrPtArea->Fill(jet->Area(), corrpt);
1086  }
1087  }
1088 }
1089 
1090 //________________________________________________________________________
1092 {
1093  AliJetContainer* jets1 = GetJetContainer(0);
1094  AliJetContainer* jets2 = GetJetContainer(1);
1095 
1096  /*
1097  Printf("Detector level jet");
1098  jet1->Print();
1099 
1100  jet1->PrintConstituents(jets1->GetParticleContainer() != 0 ? jets1->GetParticleContainer()->GetArray() : 0,
1101  jets1->GetClusterContainer() != 0 ? jets1->GetClusterContainer()->GetArray() : 0);
1102 
1103  Printf("Matched with particle level jet");
1104  jet2->Print();
1105  jet2->PrintConstituents(jets2->GetParticleContainer() != 0 ? jets2->GetParticleContainer()->GetArray() : 0,
1106  jets2->GetClusterContainer() != 0 ? jets2->GetClusterContainer()->GetArray() : 0);
1107  */
1108 
1109  AliTLorentzVector leadPart1;
1110  jets1->GetLeadingHadronMomentum(leadPart1, jet1);
1111  Double_t zleading1 = GetParallelFraction(leadPart1.Vect(), jet1);
1112  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
1113 
1114  Double_t corrpt1 = jet1->Pt() - jets1->GetRhoVal() * jet1->Area();
1115  Double_t zflavour1 = 0;
1116  Double_t ptflavour1 = 0;
1117  AliVParticle* hftrack1 = jet1->GetFlavourTrack();
1118  if (hftrack1) {
1119  zflavour1 = GetParallelFraction(hftrack1, jet1);
1120  ptflavour1 = hftrack1->Pt();
1121 
1122  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
1123  }
1124 
1125 
1126  AliTLorentzVector leadPart2;
1127  jets2->GetLeadingHadronMomentum(leadPart2, jet2);
1128  Double_t zleading2 = GetParallelFraction(leadPart2.Vect(), jet2);
1129  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
1130 
1131  Double_t corrpt2 = jet2->Pt() - jets2->GetRhoVal() * jet2->Area();
1132  Double_t zflavour2 = 0;
1133  Double_t ptflavour2 = 0;
1134  AliVParticle* hftrack2 = jet2->GetFlavourTrack();
1135  if (hftrack2) {
1136  zflavour2 = GetParallelFraction(hftrack2, jet2);
1137  ptflavour2 = hftrack2->Pt();
1138 
1139  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
1140  }
1141 
1142  if (fHistoType==1) {
1143  Double_t contents[20]={0};
1144 
1145  for (Int_t i = 0; i < fHistMatching->GetNdimensions(); i++) {
1146  TString title(fHistMatching->GetAxis(i)->GetTitle());
1147  if (title=="p_{T,1}")
1148  contents[i] = jet1->Pt();
1149  else if (title=="p_{T,2}")
1150  contents[i] = jet2->Pt();
1151  else if (title=="A_{jet,1}")
1152  contents[i] = jet1->Area();
1153  else if (title=="A_{jet,2}")
1154  contents[i] = jet2->Area();
1155  else if (title=="distance")
1156  contents[i] = d;
1157  else if (title=="CE1")
1158  contents[i] = CE1;
1159  else if (title=="CE2")
1160  contents[i] = CE2;
1161  else if (title=="#deltaA_{jet}")
1162  contents[i] = jet1->Area()-jet2->Area();
1163  else if (title=="#deltap_{T}")
1164  contents[i] = jet1->Pt()-jet2->Pt();
1165  else if (title=="#delta#eta")
1166  contents[i] = jet1->Eta()-jet2->Eta();
1167  else if (title=="#delta#phi")
1168  contents[i] = jet1->Phi()-jet2->Phi();
1169  else if (title=="p_{T,1}^{corr}")
1170  contents[i] = corrpt1;
1171  else if (title=="p_{T,2}^{corr}")
1172  contents[i] = corrpt2;
1173  else if (title=="#deltap_{T}^{corr}")
1174  contents[i] = corrpt1-corrpt2;
1175  else if (title=="p_{T,1}^{MC}")
1176  contents[i] = jet1->MCPt();
1177  else if (title=="#deltap_{T}^{MC}")
1178  contents[i] = jet1->MCPt()-jet2->Pt();
1179  else if (title=="NEF_{1}")
1180  contents[i] = jet1->NEF();
1181  else if (title=="NEF_{2}")
1182  contents[i] = jet2->NEF();
1183  else if (title=="Z_{1}")
1184  contents[i] = zleading1;
1185  else if (title=="Z_{2}")
1186  contents[i] = zleading2;
1187  else if (title=="p_{T,particle,1}^{leading} (GeV/c)")
1188  contents[i] = leadPart1.Pt();
1189  else if (title=="p_{T,particle,2}^{leading} (GeV/c)")
1190  contents[i] = leadPart2.Pt();
1191  else if (title=="z_{flavour,1}")
1192  contents[i] = zflavour1;
1193  else if (title=="z_{flavour,2}")
1194  contents[i] = zflavour2;
1195  else if (title=="p_{T,1}^{D}")
1196  contents[i] = ptflavour1;
1197  else if (title=="p_{T,2}^{D}")
1198  contents[i] = ptflavour2;
1199  else if (title=="Z_{g,1}")
1200  contents[i] = jet1->GetShapeProperties()->GetSoftDropZg();
1201  else if (title=="Z_{g,2}")
1202  contents[i] = jet2->GetShapeProperties()->GetSoftDropZg();
1203  else if (title=="dR_{1}")
1204  contents[i] = jet1->GetShapeProperties()->GetSoftDropdR();
1205  else if (title=="dR_{2}")
1206  contents[i] = jet2->GetShapeProperties()->GetSoftDropdR();
1207  else if (title=="p_{T,g,1}")
1208  contents[i] = ( jet1->GetShapeProperties()->GetSoftDropPtfrac() )*( jet1->Pt() );
1209  else if (title=="p_{T,g,2}")
1210  contents[i] = ( jet2->GetShapeProperties()->GetSoftDropPtfrac() )*( jet2->Pt() );
1211  else if (title=="DBC_{1}")
1212  contents[i] = ( jet1->GetShapeProperties()->GetSoftDropDropCount() );
1213  else if (title=="DBC_{2}")
1214  contents[i] = ( jet2->GetShapeProperties()->GetSoftDropDropCount() );
1215  else
1216  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1217  }
1218 
1219  fHistMatching->Fill(contents);
1220  }
1221  else {
1222  fHistCommonEnergy1vsJet1Pt->Fill(CE1, jet1->Pt());
1223  fHistCommonEnergy2vsJet2Pt->Fill(CE2, jet2->Pt());
1224 
1225  fHistDistancevsJet1Pt->Fill(d, jet1->Pt());
1226  fHistDistancevsJet2Pt->Fill(d, jet2->Pt());
1227 
1228  fHistDistancevsCommonEnergy1->Fill(d, CE1);
1229  fHistDistancevsCommonEnergy2->Fill(d, CE2);
1230  fHistCommonEnergy1vsCommonEnergy2->Fill(CE1, CE2);
1231 
1232  fHistDeltaEtaDeltaPhi->Fill(jet1->Eta()-jet2->Eta(),jet1->Phi()-jet2->Phi());
1233 
1234  fHistJet2PtOverJet1PtvsJet2Pt->Fill(jet2->Pt(), jet2->Pt() / jet1->Pt());
1235  fHistJet1PtOverJet2PtvsJet1Pt->Fill(jet1->Pt(), jet1->Pt() / jet2->Pt());
1236 
1237  Double_t dpt = jet1->Pt() - jet2->Pt();
1238  fHistDeltaPtvsJet1Pt->Fill(jet1->Pt(), dpt);
1239  fHistDeltaPtvsJet2Pt->Fill(jet2->Pt(), dpt);
1240  fHistDeltaPtOverJet1PtvsJet1Pt->Fill(jet1->Pt(), dpt/jet1->Pt());
1241  fHistDeltaPtOverJet2PtvsJet2Pt->Fill(jet2->Pt(), dpt/jet2->Pt());
1242 
1243  fHistDeltaPtvsDistance->Fill(d, dpt);
1244  fHistDeltaPtvsCommonEnergy1->Fill(CE1, dpt);
1245  fHistDeltaPtvsCommonEnergy2->Fill(CE2, dpt);
1246 
1247  fHistDeltaPtvsArea1->Fill(jet1->Area(), dpt);
1248  fHistDeltaPtvsArea2->Fill(jet2->Area(), dpt);
1249 
1250  Double_t darea = jet1->Area() - jet2->Area();
1251  fHistDeltaPtvsDeltaArea->Fill(darea, dpt);
1252 
1253  fHistJet1PtvsJet2Pt->Fill(jet1->Pt(), jet2->Pt());
1254 
1255  if (fIsJet1Rho || fIsJet2Rho) {
1256  Double_t dcorrpt = corrpt1 - corrpt2;
1257  fHistDeltaCorrPtvsJet1CorrPt->Fill(corrpt1, dcorrpt);
1258  fHistDeltaCorrPtvsJet2CorrPt->Fill(corrpt2, dcorrpt);
1259  fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt->Fill(corrpt1, dcorrpt/corrpt1);
1260  fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt->Fill(corrpt2, dcorrpt/corrpt2);
1261  fHistDeltaCorrPtvsDistance->Fill(d, dcorrpt);
1262  fHistDeltaCorrPtvsCommonEnergy1->Fill(CE1, dcorrpt);
1263  fHistDeltaCorrPtvsCommonEnergy2->Fill(CE2, dcorrpt);
1264  fHistDeltaCorrPtvsArea1->Fill(jet1->Area(), dcorrpt);
1265  fHistDeltaCorrPtvsArea2->Fill(jet2->Area(), dcorrpt);
1266  fHistDeltaCorrPtvsDeltaArea->Fill(darea, dcorrpt);
1267  fHistJet1CorrPtvsJet2CorrPt->Fill(corrpt1, corrpt2);
1268  }
1269 
1270  if (fIsEmbedded) {
1271  Double_t dmcpt = jet1->MCPt() - jet2->Pt();
1272  fHistDeltaMCPtvsJet1MCPt->Fill(jet1->MCPt(), dmcpt);
1273  fHistDeltaMCPtvsJet2Pt->Fill(jet2->MCPt(), dmcpt);
1274  fHistDeltaMCPtOverJet1MCPtvsJet1MCPt->Fill(jet1->MCPt(), dmcpt/jet1->MCPt());
1275  fHistDeltaMCPtOverJet2PtvsJet2Pt->Fill(jet2->Pt(), dmcpt/jet2->Pt());
1276  fHistDeltaMCPtvsDistance->Fill(d, dmcpt);
1277  fHistDeltaMCPtvsCommonEnergy1->Fill(CE1, dmcpt);
1278  fHistDeltaMCPtvsCommonEnergy2->Fill(CE2, dmcpt);
1279  fHistDeltaMCPtvsArea1->Fill(jet1->Area(), dmcpt);
1280  fHistDeltaMCPtvsArea2->Fill(jet2->Area(), dmcpt);
1281  fHistDeltaMCPtvsDeltaArea->Fill(darea, dmcpt);
1282  fHistJet1MCPtvsJet2Pt->Fill(jet1->MCPt(), jet2->Pt());
1283  }
1284  }
1285 }
1286 
1287 //________________________________________________________________________
1289 {
1290  // Execute once.
1291 
1293 
1294  AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
1295  AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
1296 
1297  if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return;
1298 
1299  if (fMatching == kMCLabel && (!jets2->GetIsParticleLevel() || jets1->GetIsParticleLevel())) {
1300  if (jets1->GetIsParticleLevel() == jets2->GetIsParticleLevel()) {
1301  AliWarning("Changing matching type from MC label to same collection...");
1303  }
1304  else {
1305  AliWarning("Changing matching type from MC label to geometrical...");
1307  }
1308  }
1309  else if (fMatching == kSameCollections && jets1->GetIsParticleLevel() != jets2->GetIsParticleLevel()) {
1310  if (jets2->GetIsParticleLevel() && !jets1->GetIsParticleLevel()) {
1311  AliWarning("Changing matching type from same collection to MC label...");
1312  fMatching = kMCLabel;
1313  }
1314  else {
1315  AliWarning("Changing matching type from same collection to geometrical...");
1317  }
1318  }
1319 }
1320 
1321 //________________________________________________________________________
1323 {
1324  // Find the closest jets
1325 
1326  if (fMatching == kNoMatching)
1327  return kTRUE;
1328  else
1329  return DoJetMatching();
1330 }
1331 
1332 //________________________________________________________________________
1334 {
1335  AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
1336  AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
1337 
1338  if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return kFALSE;
1339 
1340  DoJetLoop();
1341 
1342  AliEmcalJet* jet1 = 0;
1343 
1344  jets1->ResetCurrentID();
1345  while ((jet1 = jets1->GetNextJet())) {
1346 
1347  AliEmcalJet *jet2 = jet1->ClosestJet();
1348 
1349  if (!jet2) continue;
1350  if (jet2->ClosestJet() != jet1) continue;
1351  if (jet1->ClosestJetDistance() > fMatchingPar1 || jet2->ClosestJetDistance() > fMatchingPar2) continue;
1352 
1353  // Matched jet found
1356  AliDebug(2,Form("Found matching: jet1 pt = %f, eta = %f, phi = %f, jet2 pt = %f, eta = %f, phi = %f",
1357  jet1->Pt(), jet1->Eta(), jet1->Phi(),
1358  jet2->Pt(), jet2->Eta(), jet2->Phi()));
1359  }
1360 
1361  return kTRUE;
1362 }
1363 
1364 //________________________________________________________________________
1366 {
1367  // Do the jet loop.
1368 
1369  AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
1370  AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
1371 
1372  if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return;
1373 
1374  AliEmcalJet* jet1 = 0;
1375  AliEmcalJet* jet2 = 0;
1376 
1377  jets2->ResetCurrentID();
1378  while ((jet2 = jets2->GetNextJet())) jet2->ResetMatching();
1379 
1380  jets1->ResetCurrentID();
1381  while ((jet1 = jets1->GetNextJet())) {
1382  jet1->ResetMatching();
1383 
1384  if (jet1->MCPt() < fMinJetMCPt) continue;
1385 
1386  jets2->ResetCurrentID();
1387  while ((jet2 = jets2->GetNextJet())) {
1388  SetMatchingLevel(jet1, jet2, fMatching);
1389  } // jet2 loop
1390  } // jet1 loop
1391 }
1392 
1393 //________________________________________________________________________
1395 {
1396  d = jet1->DeltaR(jet2);
1397 }
1398 
1399 //________________________________________________________________________
1401 {
1402  AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
1403  AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
1404 
1405  if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return;
1406 
1407  AliParticleContainer *tracks1 = jets1->GetParticleContainer();
1408  AliClusterContainer *clusters1 = jets1->GetClusterContainer();
1409  AliParticleContainer *tracks2 = jets2->GetParticleContainer();
1410 
1411  // d1 and d2 represent the matching level: 0 = maximum level of matching, 1 = the two jets are completely unrelated
1412  d1 = jet1->Pt();
1413  d2 = jet2->Pt();
1414  Double_t totalPt1 = d1; // the total pt of the reconstructed jet will be cleaned from the background
1415 
1416  // remove completely tracks that are not MC particles (label == 0)
1417  if (tracks1 && tracks1->GetArray()) {
1418  for (Int_t iTrack = 0; iTrack < jet1->GetNumberOfTracks(); iTrack++) {
1419  AliVParticle *track = jet1->TrackAt(iTrack,tracks1->GetArray());
1420  if (!track) {
1421  AliWarning(Form("Could not find track %d!", iTrack));
1422  continue;
1423  }
1424 
1425  Int_t MClabel = TMath::Abs(track->GetLabel());
1426  MClabel -= fMCLabelShift;
1427  if (MClabel != 0) continue;
1428 
1429  // this is not a MC particle; remove it completely
1430  AliDebug(3,Form("Track %d (pT = %f) is not a MC particle (MClabel = %d)!",iTrack,track->Pt(),MClabel));
1431  totalPt1 -= track->Pt();
1432  d1 -= track->Pt();
1433  }
1434  }
1435 
1436  // remove completely clusters that are not MC particles (label == 0)
1437  if (fUseCellsToMatch && fCaloCells) {
1438  for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
1439  AliVCluster *clus = jet1->ClusterAt(iClus,clusters1->GetArray());
1440  if (!clus) {
1441  AliWarning(Form("Could not find cluster %d!", iClus));
1442  continue;
1443  }
1444  AliTLorentzVector part;
1445  clus->GetMomentum(part, fVertex);
1446 
1447  for (Int_t iCell = 0; iCell < clus->GetNCells(); iCell++) {
1448  Int_t cellId = clus->GetCellAbsId(iCell);
1449  Double_t cellFrac = clus->GetCellAmplitudeFraction(iCell);
1450 
1451  Int_t MClabel = TMath::Abs(fCaloCells->GetCellMCLabel(cellId));
1452  MClabel -= fMCLabelShift;
1453  if (MClabel != 0) continue;
1454 
1455  // this is not a MC particle; remove it completely
1456  AliDebug(3,Form("Cell %d (frac = %f) is not a MC particle (MClabel = %d)!",iCell,cellFrac,MClabel));
1457  totalPt1 -= part.Pt() * cellFrac;
1458  d1 -= part.Pt() * cellFrac;
1459  }
1460  }
1461  }
1462  else {
1463  for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
1464  AliVCluster *clus = jet1->ClusterAt(iClus,clusters1->GetArray());
1465  if (!clus) {
1466  AliWarning(Form("Could not find cluster %d!", iClus));
1467  continue;
1468  }
1469  TLorentzVector part;
1470  clus->GetMomentum(part, fVertex);
1471 
1472  Int_t MClabel = TMath::Abs(clus->GetLabel());
1473  MClabel -= fMCLabelShift;
1474  if (MClabel != 0) continue;
1475 
1476  // this is not a MC particle; remove it completely
1477  AliDebug(3,Form("Cluster %d (pT = %f) is not a MC particle (MClabel = %d)!",iClus,part.Pt(),MClabel));
1478  totalPt1 -= part.Pt();
1479  d1 -= part.Pt();
1480  }
1481  }
1482 
1483  for (Int_t iTrack2 = 0; iTrack2 < jet2->GetNumberOfTracks(); iTrack2++) {
1484  Bool_t track2Found = kFALSE;
1485  Int_t index2 = jet2->TrackAt(iTrack2);
1486 
1487  // now look for common particles in the track array
1488  for (Int_t iTrack = 0; iTrack < jet1->GetNumberOfTracks(); iTrack++) {
1489  AliVParticle *track = jet1->TrackAt(iTrack,tracks1->GetArray());
1490  if (!track) {
1491  AliWarning(Form("Could not find track %d!", iTrack));
1492  continue;
1493  }
1494  Int_t MClabel = TMath::Abs(track->GetLabel());
1495  MClabel -= fMCLabelShift;
1496  if (MClabel <= 0) continue;
1497 
1498  Int_t index = -1;
1499  index = tracks2->GetIndexFromLabel(MClabel);
1500  if (index < 0) {
1501  AliDebug(2,Form("Track %d (pT = %f) does not have an associated MC particle (MClabel = %d)!",iTrack,track->Pt(),MClabel));
1502  continue;
1503  }
1504 
1505  if (index2 != index) continue;
1506 
1507  // found common particle
1508  d1 -= track->Pt();
1509 
1510  if (!track2Found) {
1511  AliVParticle *MCpart = tracks2->GetParticle(index2);
1512  AliDebug(3,Form("Track %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!",
1513  iTrack,track->Pt(),track->Eta(),track->Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
1514  d2 -= MCpart->Pt();
1515  }
1516 
1517  track2Found = kTRUE;
1518  }
1519 
1520  // now look for common particles in the cluster array
1521  if (fUseCellsToMatch && fCaloCells) { // if the cell colection is available, look for cells with a matched MC particle
1522  for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
1523  AliVCluster *clus = jet1->ClusterAt(iClus,fCaloClusters);
1524  if (!clus) {
1525  AliWarning(Form("Could not find cluster %d!", iClus));
1526  continue;
1527  }
1528  AliTLorentzVector part;
1529  clus->GetMomentum(part, fVertex);
1530 
1531  for (Int_t iCell = 0; iCell < clus->GetNCells(); iCell++) {
1532  Int_t cellId = clus->GetCellAbsId(iCell);
1533  Double_t cellFrac = clus->GetCellAmplitudeFraction(iCell);
1534 
1535  Int_t MClabel = TMath::Abs(fCaloCells->GetCellMCLabel(cellId));
1536  MClabel -= fMCLabelShift;
1537  if (MClabel <= 0) continue;
1538 
1539  Int_t index1 = -1;
1540  index1 = tracks2->GetIndexFromLabel(MClabel);
1541  if (index1 < 0) {
1542  AliDebug(3,Form("Cell %d (frac = %f) does not have an associated MC particle (MClabel = %d)!",iCell,cellFrac,MClabel));
1543  continue;
1544  }
1545 
1546  if (index2 != index1) continue;
1547 
1548  // found common particle
1549  d1 -= part.Pt() * cellFrac;
1550 
1551  if (!track2Found) { // only if it is not already found among charged tracks (charged particles are most likely already found)
1552  AliVParticle *MCpart = tracks2->GetParticle(index2);
1553  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)!",
1554  iCell,iClus,part.Pt(),part.Eta(),part.Phi_0_2pi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
1555  d2 -= MCpart->Pt() * cellFrac;
1556  }
1557 
1558  track2Found = kTRUE;
1559  }
1560  }
1561  }
1562  else { //otherwise look for the first contributor to the cluster, and if matched to a MC label remove it
1563  for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
1564  AliVCluster *clus = jet1->ClusterAt(iClus,fCaloClusters);
1565  if (!clus) {
1566  AliWarning(Form("Could not find cluster %d!", iClus));
1567  continue;
1568  }
1569  AliTLorentzVector part;
1570  clus->GetMomentum(part, fVertex);
1571 
1572  Int_t MClabel = TMath::Abs(clus->GetLabel());
1573  MClabel -= fMCLabelShift;
1574  if (MClabel <= 0) continue;
1575 
1576  Int_t index = -1;
1577  index = tracks2->GetIndexFromLabel(MClabel);
1578 
1579  if (index < 0) {
1580  AliDebug(3,Form("Cluster %d (pT = %f) does not have an associated MC particle (MClabel = %d)!",iClus,part.Pt(),MClabel));
1581  continue;
1582  }
1583 
1584  if (index2 != index) continue;
1585 
1586  // found common particle
1587  d1 -= part.Pt();
1588 
1589  if (!track2Found) { // only if it is not already found among charged tracks (charged particles are most likely already found)
1590  AliVParticle *MCpart = tracks2->GetParticle(index2);
1591  AliDebug(3,Form("Cluster %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!",
1592  iClus,part.Pt(),part.Eta(),part.Phi_0_2pi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
1593 
1594  d2 -= MCpart->Pt();
1595  }
1596 
1597  track2Found = kTRUE;
1598  }
1599  }
1600  }
1601 
1602  if (d1 < 0)
1603  d1 = 0;
1604 
1605  if (d2 < 0)
1606  d2 = 0;
1607 
1608  if (totalPt1 < 1)
1609  d1 = -1;
1610  else
1611  d1 /= totalPt1;
1612 
1613  if (jet2->Pt() < 1)
1614  d2 = -1;
1615  else
1616  d2 /= jet2->Pt();
1617 }
1618 
1619 //________________________________________________________________________
1621 {
1622  AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
1623  AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
1624 
1625  if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return;
1626 
1627  AliParticleContainer *tracks1 = jets1->GetParticleContainer();
1628  AliClusterContainer *clusters1 = jets1->GetClusterContainer();
1629  AliParticleContainer *tracks2 = jets2->GetParticleContainer();
1630  AliClusterContainer *clusters2 = jets2->GetClusterContainer();
1631 
1632  // d1 and d2 represent the matching level: 0 = maximum level of matching, 1 = the two jets are completely unrelated
1633  d1 = jet1->Pt();
1634  d2 = jet2->Pt();
1635 
1636  if (tracks1 && tracks2) {
1637 
1638  for (Int_t iTrack2 = 0; iTrack2 < jet2->GetNumberOfTracks(); iTrack2++) {
1639  Int_t index2 = jet2->TrackAt(iTrack2);
1640  for (Int_t iTrack1 = 0; iTrack1 < jet1->GetNumberOfTracks(); iTrack1++) {
1641  Int_t index1 = jet1->TrackAt(iTrack1);
1642  if (index2 == index1) { // found common particle
1643  AliVParticle *part1 = tracks1->GetParticle(index1);
1644  if (!part1) {
1645  AliWarning(Form("Could not find track %d!", index1));
1646  continue;
1647  }
1648  AliVParticle *part2 = tracks2->GetParticle(index2);
1649  if (!part2) {
1650  AliWarning(Form("Could not find track %d!", index2));
1651  continue;
1652  }
1653 
1654  d1 -= part1->Pt();
1655  d2 -= part2->Pt();
1656  break;
1657  }
1658  }
1659  }
1660 
1661  }
1662 
1663  if (clusters1 && clusters2) {
1664 
1665  if (fUseCellsToMatch && fCaloCells) {
1666  // Note: this section of the code needs to be revised and tested heavily
1667  // While fixing some inconsistencies in AliAnalysisTaskEMCALClusterizeFast
1668  // some issues came up, e.g. memory leaks (fixed) and inconsistent use of
1669  // fCaloCells. In principle the two cluster collections may use cells
1670  // from different sources (embedding / non-embedding). This is not handled
1671  // correctly in the current version of this code.
1672  AliWarning("ATTENTION ATTENTION ATTENTION: this section of the AliJetResponseMaker code needs to be revised and tested before using it for physics!!!");
1673  const Int_t nClus1 = jet1->GetNumberOfClusters();
1674 
1675  Int_t ncells1[nClus1];
1676  UShort_t *cellsId1[nClus1];
1677  Double_t *cellsFrac1[nClus1];
1678  Double_t *cellsClusFrac1[nClus1];
1679  Int_t *sortedIndexes1[nClus1];
1680  Double_t ptClus1[nClus1];
1681  for (Int_t iClus1 = 0; iClus1 < nClus1; iClus1++) {
1682  Int_t index1 = jet1->ClusterAt(iClus1);
1683  AliVCluster *clus1 = clusters1->GetCluster(index1);
1684  if (!clus1) {
1685  AliWarning(Form("Could not find cluster %d!", index1));
1686  ncells1[iClus1] = 0;
1687  cellsId1[iClus1] = 0;
1688  cellsFrac1[iClus1] = 0;
1689  cellsClusFrac1[iClus1] = 0;
1690  sortedIndexes1[iClus1] = 0;
1691  ptClus1[iClus1] = 0;
1692  continue;
1693  }
1694  TLorentzVector part1;
1695  clus1->GetMomentum(part1, fVertex);
1696 
1697  ncells1[iClus1] = clus1->GetNCells();
1698  cellsId1[iClus1] = clus1->GetCellsAbsId();
1699  cellsFrac1[iClus1] = clus1->GetCellsAmplitudeFraction();
1700  cellsClusFrac1[iClus1] = new Double_t[ncells1[iClus1]];
1701  sortedIndexes1[iClus1] = new Int_t[ncells1[iClus1]];
1702  ptClus1[iClus1] = part1.Pt();
1703 
1704  for (Int_t iCell = 0; iCell < ncells1[iClus1]; iCell++) {
1705  cellsClusFrac1[iClus1][iCell] = fCaloCells->GetCellAmplitude(cellsId1[iClus1][iCell]) / clus1->E();
1706  }
1707 
1708  TMath::Sort(ncells1[iClus1], cellsId1[iClus1], sortedIndexes1[iClus1], kFALSE);
1709  }
1710 
1711  const Int_t nClus2 = jet2->GetNumberOfClusters();
1712 
1713  const Int_t maxNcells2 = 11520;
1714  Int_t sortedIndexes2[maxNcells2];
1715  for (Int_t iClus2 = 0; iClus2 < nClus2; iClus2++) {
1716  Int_t index2 = jet2->ClusterAt(iClus2);
1717  AliVCluster *clus2 = clusters2->GetCluster(index2);
1718  if (!clus2) {
1719  AliWarning(Form("Could not find cluster %d!", index2));
1720  continue;
1721  }
1722  Int_t ncells2 = clus2->GetNCells();
1723  if (ncells2 >= maxNcells2) {
1724  AliError(Form("Number of cells in the cluster %d >= %d",ncells2,maxNcells2));
1725  continue;
1726  }
1727  UShort_t *cellsId2 = clus2->GetCellsAbsId();
1728  Double_t *cellsFrac2 = clus2->GetCellsAmplitudeFraction();
1729  Double_t *cellsClusFrac2 = new Double_t[ncells2];
1730 
1731  for (Int_t iCell = 0; iCell < ncells2; iCell++) {
1732  cellsClusFrac2[iCell] = fCaloCells->GetCellAmplitude(cellsId2[iCell]) / clus2->E();
1733  }
1734 
1735  TLorentzVector part2;
1736  clus2->GetMomentum(part2, fVertex);
1737  Double_t ptClus2 = part2.Pt();
1738 
1739  TMath::Sort(ncells2, cellsId2, sortedIndexes2, kFALSE);
1740 
1741  for (Int_t iClus1 = 0; iClus1 < nClus1; iClus1++) {
1742  if (sortedIndexes1[iClus1] == 0)
1743  continue;
1744  Int_t iCell1 = 0, iCell2 = 0;
1745  while (iCell1 < ncells1[iClus1] && iCell2 < ncells2) {
1746  if (cellsId1[iClus1][sortedIndexes1[iClus1][iCell1]] == cellsId2[sortedIndexes2[iCell2]]) { // found a common cell
1747  d1 -= cellsFrac1[iClus1][sortedIndexes1[iClus1][iCell1]] * cellsClusFrac1[iClus1][sortedIndexes1[iClus1][iCell1]] * ptClus1[iClus1];
1748  d2 -= cellsFrac2[sortedIndexes2[iCell2]] * cellsClusFrac2[sortedIndexes2[iCell2]] * ptClus2;
1749  iCell1++;
1750  iCell2++;
1751  }
1752  else if (cellsId1[iClus1][sortedIndexes1[iClus1][iCell1]] > cellsId2[sortedIndexes2[iCell2]]) {
1753  iCell2++;
1754  }
1755  else {
1756  iCell1++;
1757  }
1758  }
1759  }
1760  delete[] cellsClusFrac2;
1761  }
1762  for (Int_t iClus1 = 0; iClus1 < nClus1; iClus1++) {
1763  delete[] cellsClusFrac1[iClus1];
1764  delete[] sortedIndexes1[iClus1];
1765  }
1766  }
1767  else {
1768  for (Int_t iClus2 = 0; iClus2 < jet2->GetNumberOfClusters(); iClus2++) {
1769  Int_t index2 = jet2->ClusterAt(iClus2);
1770  for (Int_t iClus1 = 0; iClus1 < jet1->GetNumberOfClusters(); iClus1++) {
1771  Int_t index1 = jet1->ClusterAt(iClus1);
1772  if (index2 == index1) { // found common particle
1773  AliVCluster *clus1 = clusters1->GetCluster(index1);
1774  if (!clus1) {
1775  AliWarning(Form("Could not find cluster %d!", index1));
1776  continue;
1777  }
1778  AliVCluster *clus2 = clusters2->GetCluster(index2);
1779  if (!clus2) {
1780  AliWarning(Form("Could not find cluster %d!", index2));
1781  continue;
1782  }
1783  TLorentzVector part1, part2;
1784  clus1->GetMomentum(part1, fVertex);
1785  clus2->GetMomentum(part2, fVertex);
1786 
1787  d1 -= part1.Pt();
1788  d2 -= part2.Pt();
1789  break;
1790  }
1791  }
1792  }
1793  }
1794  }
1795 
1796  if (d1 < 0)
1797  d1 = 0;
1798 
1799  if (d2 < 0)
1800  d2 = 0;
1801 
1802  if (jet1->Pt() > 0)
1803  d1 /= jet1->Pt();
1804  else
1805  d1 = -1;
1806 
1807  if (jet2->Pt() > 0)
1808  d2 /= jet2->Pt();
1809  else
1810  d2 = -1;
1811 }
1812 
1813 //________________________________________________________________________
1815 {
1816  Double_t d1 = -1;
1817  Double_t d2 = -1;
1818 
1819  switch (matching) {
1820  case kGeometrical:
1821  GetGeometricalMatchingLevel(jet1,jet2,d1);
1822  d2 = d1;
1823  break;
1824  case kMCLabel: // jet1 = detector level and jet2 = particle level!
1825  GetMCLabelMatchingLevel(jet1,jet2,d1,d2);
1826  break;
1827  case kSameCollections:
1828  GetSameCollectionsMatchingLevel(jet1,jet2,d1,d2);
1829  break;
1830  default:
1831  ;
1832  }
1833 
1834  if (d1 >= 0) {
1835 
1836  if (d1 < jet1->ClosestJetDistance()) {
1837  jet1->SetSecondClosestJet(jet1->ClosestJet(), jet1->ClosestJetDistance());
1838  jet1->SetClosestJet(jet2, d1);
1839  }
1840  else if (d1 < jet1->SecondClosestJetDistance()) {
1841  jet1->SetSecondClosestJet(jet2, d1);
1842  }
1843  }
1844 
1845  if (d2 >= 0) {
1846 
1847  if (d2 < jet2->ClosestJetDistance()) {
1848  jet2->SetSecondClosestJet(jet2->ClosestJet(), jet2->ClosestJetDistance());
1849  jet2->SetClosestJet(jet1, d2);
1850  }
1851  else if (d2 < jet2->SecondClosestJetDistance()) {
1852  jet2->SetSecondClosestJet(jet1, d2);
1853  }
1854  }
1855 }
1856 
1857 //________________________________________________________________________
1859 {
1860  // Fill histograms.
1861 
1862  AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
1863  AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
1864 
1865  if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return kFALSE;
1866 
1867  AliEmcalJet* jet1 = 0;
1868  AliEmcalJet* jet2 = 0;
1869 
1870  jets2->ResetCurrentID();
1871  while ((jet2 = jets2->GetNextJet())) {
1872 
1873  AliDebug(2,Form("Processing jet (2) %d", jets2->GetCurrentID()));
1874 
1875  if (jet2->Pt() < jets2->GetJetPtCut()) continue;
1876 
1877  UInt_t rejectionReason = 0;
1878  if (jets2->AcceptJet(jet2, rejectionReason))
1879  FillJetHisto(jet2, 2);
1880  else
1881  fHistRejectionReason2->Fill(jets2->GetRejectionReasonBitPosition(rejectionReason), jet2->Pt());
1882 
1883  jet1 = jet2->MatchedJet();
1884 
1885  if (!jet1) continue;
1886  rejectionReason = 0;
1887  if (!jets1->AcceptJet(jet1, rejectionReason)) continue;
1888  if (jet1->MCPt() < fMinJetMCPt) continue;
1889 
1890  Double_t d=-1, ce1=-1, ce2=-1;
1891  if (jet2->GetMatchingType() == kGeometrical) {
1892  if (jets2->GetIsParticleLevel() && !jets1->GetIsParticleLevel()) // the other way around is not supported
1893  GetMCLabelMatchingLevel(jet1, jet2, ce1, ce2);
1894  else if (jets1->GetIsParticleLevel() == jets2->GetIsParticleLevel())
1895  GetSameCollectionsMatchingLevel(jet1, jet2, ce1, ce2);
1896 
1897  d = jet2->ClosestJetDistance();
1898  }
1899  else if (jet2->GetMatchingType() == kMCLabel || jet2->GetMatchingType() == kSameCollections) {
1900  GetGeometricalMatchingLevel(jet1, jet2, d);
1901 
1902  ce1 = jet1->ClosestJetDistance();
1903  ce2 = jet2->ClosestJetDistance();
1904  }
1905 
1906  FillMatchingHistos(jet1, jet2, d, ce1, ce2);
1907  }
1908 
1909  jets1->ResetCurrentID();
1910  while ((jet1 = jets1->GetNextJet())) {
1911  UInt_t rejectionReason = 0;
1912  if (!jets1->AcceptJet(jet1, rejectionReason)) {
1913  fHistRejectionReason1->Fill(jets1->GetRejectionReasonBitPosition(rejectionReason), jet1->Pt());
1914  continue;
1915  }
1916  if (jet1->MCPt() < fMinJetMCPt) continue;
1917  AliDebug(2,Form("Processing jet (1) %d", jets1->GetCurrentID()));
1918 
1919  FillJetHisto(jet1, 1);
1920  }
1921  return kTRUE;
1922 }
void SetSecondClosestJet(AliEmcalJet *j, Double_t d)
Definition: AliEmcalJet.h:216
TH2 * fHistJets2PtArea
phi-eta distribution of jets 2
Double_t Area() const
Definition: AliEmcalJet.h:123
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:146
TH2 * fHistDeltaPtvsDeltaArea
delta pt between matched jets vs jet 2 area
AliEmcalJet * ClosestJet() const
Definition: AliEmcalJet.h:219
AliJetContainer * GetJetContainer(Int_t i=0) const
Double_t ClosestJetDistance() const
Definition: AliEmcalJet.h:220
Double_t Eta() const
Definition: AliEmcalJet.h:114
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:110
AliEmcalJet * MatchedJet() const
Definition: AliEmcalJet.h:223
Declaration of class AliTLorentzVector.
void GetSameCollectionsMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d1, Double_t &d2) const
Double_t fMinBinPt
min pt in histograms
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.
virtual Bool_t AcceptJet(Int_t i, UInt_t &rejectionReason) 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
TH2 * fHistJets2ZvsPt
Jet charged energy fraction vs. charged jet pt 2.
Bool_t fIsEmbedded
trigger, embedded signal
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:132
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
TClonesArray * fCaloClusters
!clusters
TH2 * fHistDeltaPtvsJet2Pt
delta pt between matched jets vs jet 1 pt
Short_t ClusterAt(Int_t idx) const
Definition: AliEmcalJet.h:130
TH2 * fHistJets1PhiEta
matching THnSparse
int Int_t
Definition: External.C:63
UShort_t GetNumberOfClusters() const
Definition: AliEmcalJet.h:131
unsigned int UInt_t
Definition: External.C:33
THnSparse * fHistJets1
Rejection reason vs. jet pt.
UShort_t GetMatchingType() const
Definition: AliEmcalJet.h:224
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
virtual AliVParticle * GetParticle(Int_t i=-1) const
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:102
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:215
TH2 * fHistDeltaMCPtvsArea2
jet 1 MC pt - jet2 pt vs jet 1 area
Short_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:153
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
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:253
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:141
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:217
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