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