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