]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/EMCALJetTasks/AliJetResponseMaker.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Updates
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliJetResponseMaker.cxx
CommitLineData
7549c451 1// $Id$
2949a2ec 2//
3// Emcal jet response matrix maker task.
4//
5// Author: S. Aiola
6
7#include "AliJetResponseMaker.h"
8
2949a2ec 9#include <TClonesArray.h>
10#include <TH1F.h>
11#include <TH2F.h>
b16bb001 12#include <TH3F.h>
2949a2ec 13#include <TLorentzVector.h>
14
2949a2ec 15#include "AliVCluster.h"
16#include "AliVTrack.h"
17#include "AliEmcalJet.h"
1ee1b5b8 18#include "AliGenPythiaEventHeader.h"
2949a2ec 19#include "AliMCEvent.h"
787a3c4f 20#include "AliLog.h"
cd6431de 21#include "AliRhoParameter.h"
2949a2ec 22
23ClassImp(AliJetResponseMaker)
24
25//________________________________________________________________________
26AliJetResponseMaker::AliJetResponseMaker() :
e44e8726 27 AliAnalysisTaskEmcalJet("AliJetResponseMaker", kTRUE),
cd6431de 28 fTracks2Name(""),
29 fCalo2Name(""),
30 fJets2Name(""),
31 fRho2Name(""),
32 fPtBiasJet2Track(0),
33 fPtBiasJet2Clus(0),
34 fAreCollections1MC(kFALSE),
35 fAreCollections2MC(kTRUE),
36 fMatching(kNoMatching),
7f76e479 37 fMatchingPar1(0),
38 fMatchingPar2(0),
cd6431de 39 fJet2MinEta(-999),
40 fJet2MaxEta(-999),
41 fJet2MinPhi(-999),
42 fJet2MaxPhi(-999),
4643d2e8 43 fSelectPtHardBin(-999),
1ee1b5b8 44 fPythiaHeader(0),
1ee1b5b8 45 fPtHardBin(0),
46 fNTrials(0),
cd6431de 47 fTracks2(0),
48 fCaloClusters2(0),
49 fJets2(0),
50 fRho2(0),
51 fRho2Val(0),
cfc2ac24 52 fTracks2Map(0),
1ee1b5b8 53 fHistNTrials(0),
1ee1b5b8 54 fHistEvents(0),
cd6431de 55 fHistJets1PhiEta(0),
56 fHistJets1PtArea(0),
57 fHistJets1CorrPtArea(0),
58 fHistJets2PhiEta(0),
59 fHistJets2PtArea(0),
60 fHistJets2CorrPtArea(0),
cfc2ac24 61 fHistJets2PhiEtaAcceptance(0),
62 fHistJets2PtAreaAcceptance(0),
63 fHistJets2CorrPtAreaAcceptance(0),
cd6431de 64 fHistMatchingLevelvsJet2Pt(0),
6a20534a 65 fHistDistancevsCommonEnergy1(0),
66 fHistDistancevsCommonEnergy2(0),
c560b734 67 fHistJet2PtOverJet1PtvsJet2Pt(0),
cfc2ac24 68 fHistDeltaEtaPhivsJet2Pt(0),
69 fHistDeltaPtvsJet2Pt(0),
70 fHistDeltaPtvsMatchingLevel(0),
71 fHistDeltaCorrPtvsJet2Pt(0),
72 fHistDeltaCorrPtvsMatchingLevel(0),
cd6431de 73 fHistNonMatchedJets1PtArea(0),
74 fHistNonMatchedJets2PtArea(0),
75 fHistNonMatchedJets1CorrPtArea(0),
76 fHistNonMatchedJets2CorrPtArea(0),
77 fHistJet1PtvsJet2Pt(0),
78 fHistJet1CorrPtvsJet2CorrPt(0),
79 fHistMissedJets2PtArea(0)
2949a2ec 80{
81 // Default constructor.
4643d2e8 82
a487deae 83 SetMakeGeneralHistograms(kTRUE);
2949a2ec 84}
85
86//________________________________________________________________________
87AliJetResponseMaker::AliJetResponseMaker(const char *name) :
e44e8726 88 AliAnalysisTaskEmcalJet(name, kTRUE),
cd6431de 89 fTracks2Name("MCParticles"),
90 fCalo2Name(""),
91 fJets2Name("MCJets"),
92 fRho2Name(""),
93 fPtBiasJet2Track(0),
94 fPtBiasJet2Clus(0),
95 fAreCollections1MC(kFALSE),
96 fAreCollections2MC(kTRUE),
97 fMatching(kNoMatching),
7f76e479 98 fMatchingPar1(0),
99 fMatchingPar2(0),
cd6431de 100 fJet2MinEta(-999),
101 fJet2MaxEta(-999),
102 fJet2MinPhi(-999),
103 fJet2MaxPhi(-999),
4643d2e8 104 fSelectPtHardBin(-999),
1ee1b5b8 105 fPythiaHeader(0),
1ee1b5b8 106 fPtHardBin(0),
107 fNTrials(0),
cd6431de 108 fTracks2(0),
109 fCaloClusters2(0),
110 fJets2(0),
111 fRho2(0),
112 fRho2Val(0),
cfc2ac24 113 fTracks2Map(0),
1ee1b5b8 114 fHistNTrials(0),
1ee1b5b8 115 fHistEvents(0),
cd6431de 116 fHistJets1PhiEta(0),
117 fHistJets1PtArea(0),
118 fHistJets1CorrPtArea(0),
119 fHistJets2PhiEta(0),
120 fHistJets2PtArea(0),
121 fHistJets2CorrPtArea(0),
cfc2ac24 122 fHistJets2PhiEtaAcceptance(0),
123 fHistJets2PtAreaAcceptance(0),
124 fHistJets2CorrPtAreaAcceptance(0),
cd6431de 125 fHistMatchingLevelvsJet2Pt(0),
6a20534a 126 fHistDistancevsCommonEnergy1(0),
127 fHistDistancevsCommonEnergy2(0),
c560b734 128 fHistJet2PtOverJet1PtvsJet2Pt(0),
cfc2ac24 129 fHistDeltaEtaPhivsJet2Pt(0),
130 fHistDeltaPtvsJet2Pt(0),
131 fHistDeltaPtvsMatchingLevel(0),
132 fHistDeltaCorrPtvsJet2Pt(0),
133 fHistDeltaCorrPtvsMatchingLevel(0),
cd6431de 134 fHistNonMatchedJets1PtArea(0),
135 fHistNonMatchedJets2PtArea(0),
136 fHistNonMatchedJets1CorrPtArea(0),
137 fHistNonMatchedJets2CorrPtArea(0),
138 fHistJet1PtvsJet2Pt(0),
139 fHistJet1CorrPtvsJet2CorrPt(0),
140 fHistMissedJets2PtArea(0)
2949a2ec 141{
142 // Standard constructor.
4643d2e8 143
a487deae 144 SetMakeGeneralHistograms(kTRUE);
2949a2ec 145}
146
147//________________________________________________________________________
148AliJetResponseMaker::~AliJetResponseMaker()
149{
150 // Destructor
151}
152
153//________________________________________________________________________
154void AliJetResponseMaker::UserCreateOutputObjects()
155{
156 // Create user objects.
157
a487deae 158 AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
2949a2ec 159
1ee1b5b8 160 const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
161 const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
162
163 fHistNTrials = new TH1F("fHistNTrials", "fHistNTrials", 11, 0, 11);
164 fHistNTrials->GetXaxis()->SetTitle("p_{T} hard bin");
165 fHistNTrials->GetYaxis()->SetTitle("trials");
166 fOutput->Add(fHistNTrials);
167
1ee1b5b8 168 fHistEvents = new TH1F("fHistEvents", "fHistEvents", 11, 0, 11);
169 fHistEvents->GetXaxis()->SetTitle("p_{T} hard bin");
170 fHistEvents->GetYaxis()->SetTitle("total events");
171 fOutput->Add(fHistEvents);
172
173 for (Int_t i = 1; i < 12; i++) {
174 fHistNTrials->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
1ee1b5b8 175 fHistEvents->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
176 }
177
cfc2ac24 178 fHistJets1PhiEta = new TH2F("fHistJets1PhiEta", "fHistJets1PhiEta", 40, -1, 1, 40, 0, TMath::Pi()*2);
cd6431de 179 fHistJets1PhiEta->GetXaxis()->SetTitle("#eta");
180 fHistJets1PhiEta->GetYaxis()->SetTitle("#phi");
181 fOutput->Add(fHistJets1PhiEta);
2949a2ec 182
cd6431de 183 fHistJets1PtArea = new TH2F("fHistJets1PtArea", "fHistJets1PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
184 fHistJets1PtArea->GetXaxis()->SetTitle("area");
185 fHistJets1PtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
186 fHistJets1PtArea->GetZaxis()->SetTitle("counts");
187 fOutput->Add(fHistJets1PtArea);
188
189 if (!fRhoName.IsNull()) {
190 fHistJets1CorrPtArea = new TH2F("fHistJets1CorrPtArea", "fHistJets1CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
191 fHistJets1CorrPtArea->GetXaxis()->SetTitle("area");
192 fHistJets1CorrPtArea->GetYaxis()->SetTitle("p_{T,1}^{corr} (GeV/c)");
193 fHistJets1CorrPtArea->GetZaxis()->SetTitle("counts");
194 fOutput->Add(fHistJets1CorrPtArea);
195 }
196
cfc2ac24 197 fHistJets2PhiEta = new TH2F("fHistJets2PhiEta", "fHistJets2PhiEta", 40, -1, 1, 40, 0, TMath::Pi()*2);
cd6431de 198 fHistJets2PhiEta->GetXaxis()->SetTitle("#eta");
199 fHistJets2PhiEta->GetYaxis()->SetTitle("#phi");
200 fOutput->Add(fHistJets2PhiEta);
cfc2ac24 201
202 fHistJets2PhiEtaAcceptance = new TH2F("fHistJets2PhiEtaAcceptance", "fHistJets2PhiEtaAcceptance", 40, -1, 1, 40, 0, TMath::Pi()*2);
203 fHistJets2PhiEtaAcceptance->GetXaxis()->SetTitle("#eta");
204 fHistJets2PhiEtaAcceptance->GetYaxis()->SetTitle("#phi");
205 fOutput->Add(fHistJets2PhiEtaAcceptance);
7549c451 206
cd6431de 207 fHistJets2PtArea = new TH2F("fHistJets2PtArea", "fHistJets2PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
208 fHistJets2PtArea->GetXaxis()->SetTitle("area");
209 fHistJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
210 fHistJets2PtArea->GetZaxis()->SetTitle("counts");
211 fOutput->Add(fHistJets2PtArea);
212
cfc2ac24 213 fHistJets2PtAreaAcceptance = new TH2F("fHistJets2PtAreaAcceptance", "fHistJets2PtAreaAcceptance", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
214 fHistJets2PtAreaAcceptance->GetXaxis()->SetTitle("area");
215 fHistJets2PtAreaAcceptance->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
216 fHistJets2PtAreaAcceptance->GetZaxis()->SetTitle("counts");
217 fOutput->Add(fHistJets2PtAreaAcceptance);
218
cd6431de 219 if (!fRho2Name.IsNull()) {
220 fHistJets2CorrPtArea = new TH2F("fHistJets2CorrPtArea", "fHistJets2CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
221 fHistJets2CorrPtArea->GetXaxis()->SetTitle("area");
222 fHistJets2CorrPtArea->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)");
223 fHistJets2CorrPtArea->GetZaxis()->SetTitle("counts");
224 fOutput->Add(fHistJets2CorrPtArea);
cfc2ac24 225
226 fHistJets2CorrPtAreaAcceptance = new TH2F("fHistJets2CorrPtAreaAcceptance", "fHistJets2CorrPtAreaAcceptance", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
227 fHistJets2CorrPtAreaAcceptance->GetXaxis()->SetTitle("area");
228 fHistJets2CorrPtAreaAcceptance->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)");
229 fHistJets2CorrPtAreaAcceptance->GetZaxis()->SetTitle("counts");
230 fOutput->Add(fHistJets2CorrPtAreaAcceptance);
cd6431de 231 }
232
233 fHistMatchingLevelvsJet2Pt = new TH2F("fHistMatchingLevelvsJet2Pt", "fHistMatchingLevelvsJet2Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt);
234 fHistMatchingLevelvsJet2Pt->GetXaxis()->SetTitle("Matching level");
235 fHistMatchingLevelvsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");
236 fHistMatchingLevelvsJet2Pt->GetZaxis()->SetTitle("counts");
237 fOutput->Add(fHistMatchingLevelvsJet2Pt);
238
6a20534a 239 fHistDistancevsCommonEnergy1 = new TH2F("fHistDistancevsCommonEnergy1", "fHistDistancevsCommonEnergy1", fNbins/2, 0, 1.2, fNbins/2, 0, 1.2);
240 fHistDistancevsCommonEnergy1->GetXaxis()->SetTitle("Distance");
241 fHistDistancevsCommonEnergy1->GetYaxis()->SetTitle("Common energy 1 (%)");
242 fHistDistancevsCommonEnergy1->GetZaxis()->SetTitle("counts");
243 fOutput->Add(fHistDistancevsCommonEnergy1);
244
245 fHistDistancevsCommonEnergy2 = new TH2F("fHistDistancevsCommonEnergy2", "fHistDistancevsCommonEnergy2", fNbins/2, 0, 1.2, fNbins/2, 0, 1.2);
246 fHistDistancevsCommonEnergy2->GetXaxis()->SetTitle("Distance");
247 fHistDistancevsCommonEnergy2->GetYaxis()->SetTitle("Common energy 2 (%)");
248 fHistDistancevsCommonEnergy2->GetZaxis()->SetTitle("counts");
249 fOutput->Add(fHistDistancevsCommonEnergy2);
cfc2ac24 250
c560b734 251 fHistJet2PtOverJet1PtvsJet2Pt = new TH2F("fHistJet2PtOverJet1PtvsJet2Pt", "fHistJet2PtOverJet1PtvsJet2Pt", fNbins, fMinBinPt, fMaxBinPt, 240, 0, 1.2);
252 fHistJet2PtOverJet1PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");
253 fHistJet2PtOverJet1PtvsJet2Pt->GetYaxis()->SetTitle("p_{T,2} / p_{T,1}");
254 fHistJet2PtOverJet1PtvsJet2Pt->GetZaxis()->SetTitle("counts");
255 fOutput->Add(fHistJet2PtOverJet1PtvsJet2Pt);
256
cfc2ac24 257 fHistDeltaEtaPhivsJet2Pt = new TH3F("fHistDeltaEtaPhivsJet2Pt", "fHistDeltaEtaPhivsJet2Pt", 40, -1, 1, 128, -1.6, 4.8, fNbins/2, fMinBinPt, fMaxBinPt);
258 fHistDeltaEtaPhivsJet2Pt->GetXaxis()->SetTitle("#Delta#eta");
259 fHistDeltaEtaPhivsJet2Pt->GetYaxis()->SetTitle("#Delta#phi");
260 fHistDeltaEtaPhivsJet2Pt->GetZaxis()->SetTitle("p_{T,2}");
261 fOutput->Add(fHistDeltaEtaPhivsJet2Pt);
262
c560b734 263 fHistDeltaPtvsJet2Pt = new TH2F("fHistDeltaPtvsJet2Pt", "fHistDeltaPtvsJet2Pt", fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
cfc2ac24 264 fHistDeltaPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");
265 fHistDeltaPtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
266 fHistDeltaPtvsJet2Pt->GetZaxis()->SetTitle("counts");
267 fOutput->Add(fHistDeltaPtvsJet2Pt);
268
269 fHistDeltaPtvsMatchingLevel = new TH2F("fHistDeltaPtvsMatchingLevel", "fHistDeltaPtvsMatchingLevel", fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
270 fHistDeltaPtvsMatchingLevel->GetXaxis()->SetTitle("Matching level");
271 fHistDeltaPtvsMatchingLevel->GetYaxis()->SetTitle("#Deltap_{T} (GeV/c)");
272 fHistDeltaPtvsMatchingLevel->GetZaxis()->SetTitle("counts");
273 fOutput->Add(fHistDeltaPtvsMatchingLevel);
cd6431de 274
275 if (!fRhoName.IsNull() || !fRho2Name.IsNull()) {
cfc2ac24 276 fHistDeltaCorrPtvsJet2Pt = new TH2F("fHistDeltaCorrPtvsJet2Pt", "fHistDeltaCorrPtvsJet2Pt", fNbins/2, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
277 fHistDeltaCorrPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");
278 fHistDeltaCorrPtvsJet2Pt->GetYaxis()->SetTitle("#Deltap_{T}^{corr} (GeV/c)");
279 fHistDeltaCorrPtvsJet2Pt->GetZaxis()->SetTitle("counts");
280 fOutput->Add(fHistDeltaCorrPtvsJet2Pt);
281
282 fHistDeltaCorrPtvsMatchingLevel = new TH2F("fHistDeltaCorrPtvsMatchingLevel", "fHistDeltaCorrPtvsMatchingLevel", fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
283 fHistDeltaCorrPtvsMatchingLevel->GetXaxis()->SetTitle("Matching level");
284 fHistDeltaCorrPtvsMatchingLevel->GetYaxis()->SetTitle("#Deltap_{T}^{corr} (GeV/c)");
285 fHistDeltaCorrPtvsMatchingLevel->GetZaxis()->SetTitle("counts");
787a3c4f 286 fOutput->Add(fHistDeltaCorrPtvsMatchingLevel);
cd6431de 287 }
288
289 fHistNonMatchedJets1PtArea = new TH2F("fHistNonMatchedJets1PtArea", "fHistNonMatchedJets1PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
290 fHistNonMatchedJets1PtArea->GetXaxis()->SetTitle("area");
291 fHistNonMatchedJets1PtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
292 fHistNonMatchedJets1PtArea->GetZaxis()->SetTitle("counts");
293 fOutput->Add(fHistNonMatchedJets1PtArea);
294
295 fHistNonMatchedJets2PtArea = new TH2F("fHistNonMatchedJets2PtArea", "fHistNonMatchedJets2PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
296 fHistNonMatchedJets2PtArea->GetXaxis()->SetTitle("area");
297 fHistNonMatchedJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
298 fHistNonMatchedJets2PtArea->GetZaxis()->SetTitle("counts");
299 fOutput->Add(fHistNonMatchedJets2PtArea);
300
301 if (!fRhoName.IsNull()) {
302 fHistNonMatchedJets1CorrPtArea = new TH2F("fHistNonMatchedJets1CorrPtArea", "fHistNonMatchedJets1CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
303 fHistNonMatchedJets1CorrPtArea->GetXaxis()->SetTitle("area");
304 fHistNonMatchedJets1CorrPtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
305 fHistNonMatchedJets1CorrPtArea->GetZaxis()->SetTitle("counts");
306 fOutput->Add(fHistNonMatchedJets1CorrPtArea);
307 }
308
309 if (!fRho2Name.IsNull()) {
310 fHistNonMatchedJets2CorrPtArea = new TH2F("fHistNonMatchedJets2CorrPtArea", "fHistNonMatchedJets2CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
311 fHistNonMatchedJets2CorrPtArea->GetXaxis()->SetTitle("area");
312 fHistNonMatchedJets2CorrPtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
313 fHistNonMatchedJets2CorrPtArea->GetZaxis()->SetTitle("counts");
314 fOutput->Add(fHistNonMatchedJets2CorrPtArea);
315 }
316
317 fHistJet1PtvsJet2Pt = new TH2F("fHistJet1PtvsJet2Pt", "fHistJet1PtvsJet2Pt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
318 fHistJet1PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,1}");
319 fHistJet1PtvsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");
320 fHistJet1PtvsJet2Pt->GetZaxis()->SetTitle("counts");
321 fOutput->Add(fHistJet1PtvsJet2Pt);
322
323 if (!fRhoName.IsNull() || !fRho2Name.IsNull()) {
324 if (fRhoName.IsNull())
325 fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt", fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
326 else if (fRho2Name.IsNull())
327 fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt", 2*fNbins, -fMaxBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
328 else
329 fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt", 2*fNbins, -fMaxBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
330 fHistJet1CorrPtvsJet2CorrPt->GetXaxis()->SetTitle("p_{T,1}^{corr}");
331 fHistJet1CorrPtvsJet2CorrPt->GetYaxis()->SetTitle("p_{T,2}^{corr}");
332 fHistJet1CorrPtvsJet2CorrPt->GetZaxis()->SetTitle("counts");
333 fOutput->Add(fHistJet1CorrPtvsJet2CorrPt);
334 }
335
336 fHistMissedJets2PtArea = new TH2F("fHistMissedJets2PtArea", "fHistMissedJets2PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
337 fHistMissedJets2PtArea->GetXaxis()->SetTitle("area");
338 fHistMissedJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
339 fHistMissedJets2PtArea->GetZaxis()->SetTitle("counts");
340 fOutput->Add(fHistMissedJets2PtArea);
2bddb6ae 341
99c67c1b 342 PostData(1, fOutput); // Post data for ALL output slots > 0 here, to get at least an empty histogram
2949a2ec 343}
344
ceefbfbc 345//________________________________________________________________________
a487deae 346Bool_t AliJetResponseMaker::AcceptJet(AliEmcalJet *jet) const
ceefbfbc 347{
348 // Return true if jet is accepted.
349
350 if (jet->Pt() <= fJetPtCut)
351 return kFALSE;
352 if (jet->Area() <= fJetAreaCut)
353 return kFALSE;
65bb5510 354
355 return kTRUE;
ceefbfbc 356}
357
cd6431de 358//________________________________________________________________________
359Bool_t AliJetResponseMaker::AcceptBiasJet2(AliEmcalJet *jet) const
360{
361 // Accept jet with a bias.
362
363 if (fLeadingHadronType == 0) {
364 if (jet->MaxTrackPt() < fPtBiasJet2Track) return kFALSE;
365 }
366 else if (fLeadingHadronType == 1) {
367 if (jet->MaxClusterPt() < fPtBiasJet2Clus) return kFALSE;
368 }
369 else {
370 if (jet->MaxTrackPt() < fPtBiasJet2Track && jet->MaxClusterPt() < fPtBiasJet2Clus) return kFALSE;
371 }
372
373 return kTRUE;
374}
375
b3ccdfe2 376//________________________________________________________________________
377void AliJetResponseMaker::ExecOnce()
378{
379 // Execute once.
380
cd6431de 381 if (!fJets2Name.IsNull() && !fJets2) {
382 fJets2 = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJets2Name));
383 if (!fJets2) {
384 AliError(Form("%s: Could not retrieve jets2 %s!", GetName(), fJets2Name.Data()));
b3ccdfe2 385 return;
386 }
cd6431de 387 else if (!fJets2->GetClass()->GetBaseClass("AliEmcalJet")) {
388 AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fJets2Name.Data()));
389 fJets2 = 0;
b3ccdfe2 390 return;
391 }
392 }
393
cd6431de 394 if (!fTracks2Name.IsNull() && !fTracks2) {
395 fTracks2 = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracks2Name));
396 if (!fTracks2) {
397 AliError(Form("%s: Could not retrieve tracks2 %s!", GetName(), fTracks2Name.Data()));
b3ccdfe2 398 return;
399 }
400 else {
cd6431de 401 TClass *cl = fTracks2->GetClass();
b3ccdfe2 402 if (!cl->GetBaseClass("AliVParticle") && !cl->GetBaseClass("AliEmcalParticle")) {
cd6431de 403 AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fTracks2Name.Data()));
404 fTracks2 = 0;
b3ccdfe2 405 return;
406 }
407 }
cfc2ac24 408
409 if (fAreCollections2MC) {
410 fTracks2Map = dynamic_cast<TH1*>(InputEvent()->FindListObject(fTracks2Name + "_Map"));
7f76e479 411 // this is needed to map the MC labels with the indexes of the MC particle collection
412 // if teh map is not given, the MC labels are assumed to be consistent with the indexes (which is not the case if AliEmcalMCTrackSelector is used)
cfc2ac24 413 if (!fTracks2Map) {
7f76e479 414 AliWarning(Form("%s: Could not retrieve map for tracks2 %s! Will assume MC labels consistent with indexes...", GetName(), fTracks2Name.Data()));
415 fTracks2Map = new TH1I("tracksMap","tracksMap",9999,0,1);
416 for (Int_t i = 0; i < 9999; i++) {
417 fTracks2Map->SetBinContent(i,i);
418 }
cfc2ac24 419 }
420 }
b3ccdfe2 421 }
422
cd6431de 423 if (!fCalo2Name.IsNull() && !fCaloClusters2) {
424 fCaloClusters2 = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCalo2Name));
425 if (!fCaloClusters2) {
426 AliError(Form("%s: Could not retrieve clusters %s!", GetName(), fCalo2Name.Data()));
427 return;
428 } else {
429 TClass *cl = fCaloClusters2->GetClass();
430 if (!cl->GetBaseClass("AliVCluster") && !cl->GetBaseClass("AliEmcalParticle")) {
431 AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fCalo2Name.Data()));
432 fCaloClusters2 = 0;
433 return;
434 }
435 }
91bee8bc 436 }
cd6431de 437
438 if (!fRho2Name.IsNull() && !fRho2) {
439 fRho2 = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRho2Name));
440 if (!fRho2) {
441 AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRho2Name.Data()));
442 fInitialized = kFALSE;
443 return;
444 }
91bee8bc 445 }
cd6431de 446
447 if (fJet2MinEta == -999)
448 fJet2MinEta = fJetMinEta - fJetRadius;
449 if (fJet2MaxEta == -999)
450 fJet2MaxEta = fJetMaxEta + fJetRadius;
451 if (fJet2MinPhi == -999)
452 fJet2MinPhi = fJetMinPhi - fJetRadius;
453 if (fJet2MaxPhi == -999)
454 fJet2MaxPhi = fJetMaxPhi + fJetRadius;
455
456 AliAnalysisTaskEmcalJet::ExecOnce();
b3ccdfe2 457}
458
4643d2e8 459//________________________________________________________________________
460Bool_t AliJetResponseMaker::IsEventSelected()
461{
462 // Check if event is selected
463
464 if (fSelectPtHardBin != -999 && fSelectPtHardBin != fPtHardBin)
465 return kFALSE;
466
467 return AliAnalysisTaskEmcalJet::IsEventSelected();
468}
469
b3ccdfe2 470//________________________________________________________________________
471Bool_t AliJetResponseMaker::RetrieveEventObjects()
472{
473 // Retrieve event objects.
474
475 if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
476 return kFALSE;
cd6431de 477
478 if (fRho2)
479 fRho2Val = fRho2->GetVal();
b3ccdfe2 480
b3ccdfe2 481 const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
482 const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
787a3c4f 483
484 if (MCEvent())
485 fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
486
487 if (fPythiaHeader) {
488 Double_t pthard = fPythiaHeader->GetPtHard();
489
490 for (fPtHardBin = 0; fPtHardBin < 11; fPtHardBin++) {
491 if (pthard >= ptHardLo[fPtHardBin] && pthard < ptHardHi[fPtHardBin])
492 break;
493 }
494
495 fNTrials = fPythiaHeader->Trials();
b3ccdfe2 496 }
b3ccdfe2 497
b3ccdfe2 498 return kTRUE;
499}
500
501//________________________________________________________________________
502Bool_t AliJetResponseMaker::Run()
503{
504 // Find the closest jets
cfc2ac24 505
f660c2d6 506 if (fMatching == kNoMatching)
aa4d701c 507 return kTRUE;
f660c2d6 508 else
509 return DoJetMatching();
cfc2ac24 510}
aa4d701c 511
cfc2ac24 512//________________________________________________________________________
f660c2d6 513Bool_t AliJetResponseMaker::DoJetMatching()
cfc2ac24 514{
515 DoJetLoop(kFALSE);
cfc2ac24 516
f660c2d6 517 const Int_t nJets = fJets->GetEntriesFast();
cfc2ac24 518
f660c2d6 519 for (Int_t i = 0; i < nJets; i++) {
cfc2ac24 520
521 AliEmcalJet* jet1 = static_cast<AliEmcalJet*>(fJets->At(i));
522
523 if (!jet1) {
524 AliError(Form("Could not receive jet %d", i));
525 continue;
526 }
527
528 if (!AcceptJet(jet1))
529 continue;
530
531 if (jet1->Eta() < fJetMinEta || jet1->Eta() > fJetMaxEta || jet1->Phi() < fJetMinPhi || jet1->Phi() > fJetMaxPhi)
532 continue;
533
f660c2d6 534 if (jet1->ClosestJet() && jet1->ClosestJet()->ClosestJet() == jet1 &&
7f76e479 535 jet1->ClosestJetDistance() < fMatchingPar1 && jet1->ClosestJet()->ClosestJetDistance() < fMatchingPar2) { // Matched jet found
cfc2ac24 536 jet1->SetMatchedToClosest(fMatching);
cfc2ac24 537 jet1->ClosestJet()->SetMatchedToClosest(fMatching);
aa4d701c 538 }
539 }
540
b3ccdfe2 541 return kTRUE;
542}
543
2949a2ec 544//________________________________________________________________________
cfc2ac24 545void AliJetResponseMaker::DoJetLoop(Bool_t order)
2949a2ec 546{
44476be7 547 // Do the jet loop.
1ee1b5b8 548
cfc2ac24 549 TClonesArray *jets1 = 0;
550 TClonesArray *jets2 = 0;
551
552 if (order) {
553 jets1 = fJets2;
554 jets2 = fJets;
555 }
556 else {
557 jets1 = fJets;
558 jets2 = fJets2;
559 }
560
44476be7 561 Int_t nJets1 = jets1->GetEntriesFast();
562 Int_t nJets2 = jets2->GetEntriesFast();
1ee1b5b8 563
44476be7 564 for (Int_t i = 0; i < nJets1; i++) {
1ee1b5b8 565
44476be7 566 AliEmcalJet* jet1 = static_cast<AliEmcalJet*>(jets1->At(i));
2949a2ec 567
44476be7 568 if (!jet1) {
569 AliError(Form("Could not receive jet %d", i));
570 continue;
571 }
572
ceefbfbc 573 if (!AcceptJet(jet1))
44476be7 574 continue;
575
cfc2ac24 576 if (order) {
cd6431de 577 if (jet1->Eta() < fJet2MinEta || jet1->Eta() > fJet2MaxEta || jet1->Phi() < fJet2MinPhi || jet1->Phi() > fJet2MaxPhi)
ceefbfbc 578 continue;
579 }
580 else {
cd6431de 581 if (jet1->Eta() < fJetMinEta || jet1->Eta() > fJetMaxEta || jet1->Phi() < fJetMinPhi || jet1->Phi() > fJetMaxPhi)
ceefbfbc 582 continue;
583 }
584
44476be7 585 for (Int_t j = 0; j < nJets2; j++) {
586
587 AliEmcalJet* jet2 = static_cast<AliEmcalJet*>(jets2->At(j));
588
589 if (!jet2) {
590 AliError(Form("Could not receive jet %d", j));
591 continue;
592 }
593
ceefbfbc 594 if (!AcceptJet(jet2))
44476be7 595 continue;
ceefbfbc 596
cfc2ac24 597 if (order) {
a487deae 598 if (jet2->Eta() < fJetMinEta || jet2->Eta() > fJetMaxEta || jet2->Phi() < fJetMinPhi || jet2->Phi() > fJetMaxPhi)
ceefbfbc 599 continue;
600 }
601 else {
cd6431de 602 if (jet1->Eta() < fJet2MinEta || jet1->Eta() > fJet2MaxEta || jet1->Phi() < fJet2MinPhi || jet1->Phi() > fJet2MaxPhi)
ceefbfbc 603 continue;
604 }
cd6431de 605
7f76e479 606 SetMatchingLevel(jet1, jet2, fMatching);
44476be7 607 }
608 }
1ee1b5b8 609}
610
cd6431de 611//________________________________________________________________________
7f76e479 612void AliJetResponseMaker::GetGeometricalMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d) const
cd6431de 613{
7f76e479 614 Double_t deta = jet2->Eta() - jet1->Eta();
615 Double_t dphi = jet2->Phi() - jet1->Phi();
616 d = TMath::Sqrt(deta * deta + dphi * dphi);
617}
cd6431de 618
7f76e479 619//________________________________________________________________________
620void AliJetResponseMaker::GetMCLabelMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d1, Double_t &d2) const
621{
622 // d1 and d2 represent the matching level: 0 = maximum level of matching, 1 = the two jets are completely unrelated
623 d1 = jet1->Pt();
624 d2 = jet2->Pt();
625 Double_t totalPt1 = d1; // the total pt of the reconstructed jet will be cleaned from the background
626
627 for (Int_t iTrack2 = 0; iTrack2 < jet2->GetNumberOfTracks(); iTrack2++) {
628 Bool_t track2Found = kFALSE;
629 Int_t index2 = jet2->TrackAt(iTrack2);
630 for (Int_t iTrack = 0; iTrack < jet1->GetNumberOfTracks(); iTrack++) {
631 AliVParticle *track = jet1->TrackAt(iTrack,fTracks);
632 if (!track) {
633 AliWarning(Form("Could not find track %d!", iTrack));
634 continue;
635 }
787a3c4f 636 Int_t MClabel = TMath::Abs(track->GetLabel());
7f76e479 637 Int_t index = -1;
638
787a3c4f 639 if (MClabel == 0) {// this is not a MC particle; remove it completely
7f76e479 640 AliDebug(3,Form("Track %d (pT = %f) is not a MC particle (MClabel = %d)!",iTrack,track->Pt(),MClabel));
641 totalPt1 -= track->Pt();
642 d1 -= track->Pt();
643 continue;
644 }
645 else if (MClabel < fTracks2Map->GetNbinsX()-2) {
646 index = fTracks2Map->GetBinContent(MClabel);
647 }
648
649 if (index < 0) {
650 AliDebug(2,Form("Track %d (pT = %f) does not have an associated MC particle (MClabel = %d)!",iTrack,track->Pt(),MClabel));
651 continue;
652 }
653
654 if (index2 == index) { // found common particle
655 track2Found = kTRUE;
656 d1 -= track->Pt();
657 AliVParticle *MCpart = static_cast<AliVParticle*>(fTracks2->At(index2));
658 AliDebug(3,Form("Track %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!",
659 iTrack,track->Pt(),track->Eta(),track->Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
660 d2 -= MCpart->Pt();
661 break;
662 }
cd6431de 663 }
7f76e479 664 for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
665 AliVCluster *clus = jet1->ClusterAt(iClus,fCaloClusters);
666 if (!clus) {
667 AliWarning(Form("Could not find cluster %d!", iClus));
668 continue;
669 }
670 TLorentzVector part;
671 clus->GetMomentum(part, const_cast<Double_t*>(fVertex));
672
673 if (fCaloCells) { // if the cell colection is available, look for cells with a matched MC particle
674 for (Int_t iCell = 0; iCell < clus->GetNCells(); iCell++) {
675 Int_t cellId = clus->GetCellAbsId(iCell);
676 Double_t cellFrac = clus->GetCellAmplitudeFraction(iCell);
677
787a3c4f 678 Int_t MClabel = TMath::Abs(fCaloCells->GetCellMCLabel(cellId));
7f76e479 679 Int_t index = -1;
680
787a3c4f 681 if (MClabel == 0) {// this is not a MC particle; remove it completely
7f76e479 682 AliDebug(3,Form("Cell %d (frac = %f) is not a MC particle (MClabel = %d)!",iCell,cellFrac,MClabel));
683 totalPt1 -= part.Pt() * cellFrac;
684 d1 -= part.Pt() * cellFrac;
685 continue;
686 }
687 else if (MClabel < fTracks2Map->GetNbinsX()-2) {
688 index = fTracks2Map->GetBinContent(MClabel);
689 }
690
691 if (index < 0) {
692 AliDebug(3,Form("Cell %d (frac = %f) does not have an associated MC particle (MClabel = %d)!",iCell,cellFrac,MClabel));
693 continue;
694 }
695 if (index2 == index) { // found common particle
696 d1 -= part.Pt() * cellFrac;
697
698 if (!track2Found) {// only if it is not already found among charged tracks (charged particles are most likely already found)
699 AliVParticle *MCpart = static_cast<AliVParticle*>(fTracks2->At(index2));
700 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)!",
701 iCell,iClus,part.Pt(),part.Eta(),part.Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
702 d2 -= MCpart->Pt() * cellFrac;
703 }
704 break;
705 }
cfc2ac24 706 }
707 }
7f76e479 708 else { //otherwise look for the first contributor to the cluster, and if matched to a MC label remove it
787a3c4f 709 Int_t MClabel = TMath::Abs(clus->GetLabel());
7f76e479 710 Int_t index = -1;
711
787a3c4f 712 if (MClabel == 0) {// this is not a MC particle; remove it completely
7f76e479 713 AliDebug(3,Form("Cluster %d (pT = %f) is not a MC particle (MClabel = %d)!",iClus,part.Pt(),MClabel));
714 totalPt1 -= part.Pt();
715 d1 -= part.Pt();
cfc2ac24 716 continue;
717 }
7f76e479 718 else if (MClabel < fTracks2Map->GetNbinsX()-2) {
719 index = fTracks2Map->GetBinContent(MClabel);
f660c2d6 720 }
7f76e479 721
cfc2ac24 722 if (index < 0) {
7f76e479 723 AliDebug(3,Form("Cluster %d (pT = %f) does not have an associated MC particle (MClabel = %d)!",iClus,part.Pt(),MClabel));
cfc2ac24 724 continue;
725 }
7f76e479 726 if (index2 == index) { // found common particle
727 d1 -= part.Pt();
728
729 if (!track2Found) {// only if it is not already found among charged tracks (charged particles are most likely already found)
f660c2d6 730 AliVParticle *MCpart = static_cast<AliVParticle*>(fTracks2->At(index2));
7f76e479 731 AliDebug(3,Form("Cluster %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!",
732 iClus,part.Pt(),part.Eta(),part.Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
733
f660c2d6 734 d2 -= MCpart->Pt();
cfc2ac24 735 }
7f76e479 736 break;
cfc2ac24 737 }
738 }
7f76e479 739 }
740 }
741 if (d1 <= 0 || totalPt1 < 1)
742 d1 = 0;
743 else
744 d1 /= totalPt1;
f660c2d6 745
7f76e479 746 if (jet2->Pt() > 0 && d2 > 0)
747 d2 /= jet2->Pt();
748 else
749 d2 = 0;
750}
751
752//________________________________________________________________________
753void AliJetResponseMaker::GetSameCollectionsMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d1, Double_t &d2) const
754{
755 // d1 and d2 represent the matching level: 0 = maximum level of matching, 1 = the two jets are completely unrelated
756 d1 = jet1->Pt();
757 d2 = jet2->Pt();
758
759 if (fTracks && fTracks2) {
760
761 for (Int_t iTrack2 = 0; iTrack2 < jet2->GetNumberOfTracks(); iTrack2++) {
762 Int_t index2 = jet2->TrackAt(iTrack2);
763 for (Int_t iTrack = 0; iTrack < jet1->GetNumberOfTracks(); iTrack++) {
764 Int_t index = jet1->TrackAt(iTrack);
765 if (index2 == index) { // found common particle
766 AliVParticle *part = static_cast<AliVParticle*>(fTracks->At(index));
767 if (!part) {
768 AliWarning(Form("Could not find track %d!", index));
f660c2d6 769 continue;
770 }
7f76e479 771 AliVParticle *part2 = static_cast<AliVParticle*>(fTracks2->At(index2));
772 if (!part2) {
773 AliWarning(Form("Could not find track %d!", index2));
cfc2ac24 774 continue;
775 }
7f76e479 776
777 d1 -= part->Pt();
778 d2 -= part2->Pt();
779 break;
cfc2ac24 780 }
781 }
cfc2ac24 782 }
7f76e479 783
784 }
785
786 if (fCaloClusters && fCaloClusters2) {
787
788 for (Int_t iClus2 = 0; iClus2 < jet2->GetNumberOfClusters(); iClus2++) {
789 Int_t index2 = jet2->ClusterAt(iClus2);
790 for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
791 Int_t index = jet1->ClusterAt(iClus);
792 AliVCluster *clus = static_cast<AliVCluster*>(fCaloClusters->At(index));
793 if (!clus) {
794 AliWarning(Form("Could not find cluster %d!", index));
795 continue;
796 }
797 AliVCluster *clus2 = static_cast<AliVCluster*>(fCaloClusters2->At(index2));
798 if (!clus2) {
799 AliWarning(Form("Could not find cluster %d!", index2));
800 continue;
801 }
802 TLorentzVector part, part2;
803 clus->GetMomentum(part, const_cast<Double_t*>(fVertex));
804 clus2->GetMomentum(part2, const_cast<Double_t*>(fVertex));
805
806 d1 -= part.Pt();
807 d2 -= part2.Pt();
808 break;
809 }
810 }
811
812 }
813
814 if (jet1->Pt() > 0 && d1 > 0)
815 d1 /= jet1->Pt();
816 else
817 d1 = 0;
818
819 if (jet2->Pt() > 0 && d2 > 0)
820 d2 /= jet2->Pt();
821 else
822 d2 = 0;
823}
824
825//________________________________________________________________________
826void AliJetResponseMaker::SetMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, MatchingType matching)
827{
828 Double_t d1 = -1;
829 Double_t d2 = -1;
830
831 switch (matching) {
832 case kGeometrical:
833 GetGeometricalMatchingLevel(jet1,jet2,d1);
834 d2 = d1;
835 break;
836 case kMCLabel: // jet1 = detector level and jet2 = particle level!
837 GetMCLabelMatchingLevel(jet1,jet2,d1,d2);
838 break;
839 case kSameCollections:
840 GetSameCollectionsMatchingLevel(jet1,jet2,d1,d2);
cd6431de 841 break;
842 default:
843 ;
844 }
f660c2d6 845
846 if (d1 > 0) {
847
848 if (d1 < jet1->ClosestJetDistance()) {
849 jet1->SetSecondClosestJet(jet1->ClosestJet(), jet1->ClosestJetDistance());
850 jet1->SetClosestJet(jet2, d1);
851 }
852 else if (d1 < jet1->SecondClosestJetDistance()) {
853 jet1->SetSecondClosestJet(jet2, d1);
854 }
855 }
856
857 if (d2 > 0) {
858
859 if (d2 < jet2->ClosestJetDistance()) {
860 jet2->SetSecondClosestJet(jet2->ClosestJet(), jet2->ClosestJetDistance());
861 jet2->SetClosestJet(jet1, d2);
862 }
863 else if (d2 < jet2->SecondClosestJetDistance()) {
864 jet2->SetSecondClosestJet(jet1, d2);
865 }
866 }
cd6431de 867}
868
1ee1b5b8 869//________________________________________________________________________
870Bool_t AliJetResponseMaker::FillHistograms()
871{
872 // Fill histograms.
873
4643d2e8 874 fHistEvents->SetBinContent(fPtHardBin + 1, fHistEvents->GetBinContent(fPtHardBin + 1) + 1);
875 fHistNTrials->SetBinContent(fPtHardBin + 1, fHistNTrials->GetBinContent(fPtHardBin + 1) + fNTrials);
4643d2e8 876
cd6431de 877 const Int_t nJets2 = fJets2->GetEntriesFast();
2949a2ec 878
cd6431de 879 for (Int_t i = 0; i < nJets2; i++) {
2949a2ec 880
cd6431de 881 AliEmcalJet* jet2 = static_cast<AliEmcalJet*>(fJets2->At(i));
2949a2ec 882
cd6431de 883 if (!jet2) {
884 AliError(Form("Could not receive jet2 %d", i));
2949a2ec 885 continue;
cfc2ac24 886 }
887
cd6431de 888 if (!AcceptJet(jet2))
889 continue;
890
cfc2ac24 891 if (AcceptBiasJet(jet2) &&
892 (jet2->Eta() > fJetMinEta && jet2->Eta() < fJetMaxEta && jet2->Phi() > fJetMinPhi && jet2->Phi() < fJetMaxPhi)) {
893
894 fHistJets2PtAreaAcceptance->Fill(jet2->Area(), jet2->Pt());
895 fHistJets2PhiEtaAcceptance->Fill(jet2->Eta(), jet2->Phi());
896
897 if (!fRho2Name.IsNull())
898 fHistJets2CorrPtAreaAcceptance->Fill(jet2->Area(), jet2->Pt() - fRho2Val * jet2->Area());
899 }
900
cd6431de 901 if (!AcceptBiasJet2(jet2))
ceefbfbc 902 continue;
903
cd6431de 904 if (jet2->Eta() < fJet2MinEta || jet2->Eta() > fJet2MaxEta || jet2->Phi() < fJet2MinPhi || jet2->Phi() > fJet2MaxPhi)
2bddb6ae 905 continue;
906
cfc2ac24 907 fHistJets2PtArea->Fill(jet2->Area(), jet2->Pt());
908 fHistJets2PhiEta->Fill(jet2->Eta(), jet2->Phi());
909
910 if (!fRho2Name.IsNull())
911 fHistJets2CorrPtArea->Fill(jet2->Area(), jet2->Pt() - fRho2Val * jet2->Area());
2949a2ec 912
cd6431de 913 if (jet2->MatchedJet()) {
aa4d701c 914
cd6431de 915 if (!AcceptBiasJet(jet2->MatchedJet()) ||
916 jet2->MatchedJet()->MaxTrackPt() > fMaxTrackPt || jet2->MatchedJet()->MaxClusterPt() > fMaxClusterPt ||
917 jet2->MatchedJet()->Pt() > fMaxBinPt) {
918 fHistMissedJets2PtArea->Fill(jet2->Area(), jet2->Pt());
2bddb6ae 919 }
920 else {
7f76e479 921 Double_t d1=-1, d2=-1;
922 if (jet2->GetMatchingType() == kGeometrical) {
923
924 if (fAreCollections2MC && !fAreCollections1MC)
925 GetMCLabelMatchingLevel(jet2->MatchedJet(), jet2, d1, d2);
926 else if ((fAreCollections1MC && fAreCollections2MC) || (!fAreCollections1MC && !fAreCollections2MC))
927 GetSameCollectionsMatchingLevel(jet2->MatchedJet(), jet2, d1, d2);
928
6a20534a 929 fHistDistancevsCommonEnergy1->Fill(jet2->ClosestJetDistance(), d1);
930 fHistDistancevsCommonEnergy2->Fill(jet2->ClosestJetDistance(), d2);
7f76e479 931 }
932 else if (jet2->GetMatchingType() == kMCLabel || jet2->GetMatchingType() == kSameCollections) {
933 GetGeometricalMatchingLevel(jet2->MatchedJet(), jet2, d1);
6a20534a 934 fHistDistancevsCommonEnergy1->Fill(d1, jet2->MatchedJet()->ClosestJetDistance());
935 fHistDistancevsCommonEnergy2->Fill(d1, jet2->ClosestJetDistance());
7f76e479 936 }
cfc2ac24 937
cd6431de 938 fHistMatchingLevelvsJet2Pt->Fill(jet2->ClosestJetDistance(), jet2->Pt());
939
940 Double_t deta = jet2->MatchedJet()->Eta() - jet2->Eta();
941 Double_t dphi = jet2->MatchedJet()->Phi() - jet2->Phi();
cfc2ac24 942 fHistDeltaEtaPhivsJet2Pt->Fill(deta, dphi, jet2->Pt());
cd6431de 943
944 Double_t dpt = jet2->MatchedJet()->Pt() - jet2->Pt();
cfc2ac24 945 fHistDeltaPtvsJet2Pt->Fill(jet2->Pt(), dpt);
946 fHistDeltaPtvsMatchingLevel->Fill(jet2->ClosestJetDistance(), dpt);
cd6431de 947
c560b734 948 fHistJet2PtOverJet1PtvsJet2Pt->Fill(jet2->Pt(),jet2->Pt() / jet2->MatchedJet()->Pt());
949
cd6431de 950 fHistJet1PtvsJet2Pt->Fill(jet2->MatchedJet()->Pt(), jet2->Pt());
951
952 if (!fRhoName.IsNull() || !fRho2Name.IsNull()) {
953 dpt -= fRhoVal * jet2->MatchedJet()->Area() - fRho2Val * jet2->Area();
cfc2ac24 954 fHistDeltaCorrPtvsJet2Pt->Fill(jet2->Pt(), dpt);
955 fHistDeltaCorrPtvsMatchingLevel->Fill(jet2->ClosestJetDistance(), dpt);
cd6431de 956 fHistJet1CorrPtvsJet2CorrPt->Fill(jet2->MatchedJet()->Pt() - fRhoVal * jet2->MatchedJet()->Area(), jet2->Pt() - fRho2Val * jet2->Area());
957 }
2bddb6ae 958 }
2949a2ec 959 }
99c67c1b 960 else {
cd6431de 961 fHistNonMatchedJets2PtArea->Fill(jet2->Area(), jet2->Pt());
962 fHistMissedJets2PtArea->Fill(jet2->Area(), jet2->Pt());
2949a2ec 963
cd6431de 964 if (!fRho2Name.IsNull())
965 fHistNonMatchedJets2CorrPtArea->Fill(jet2->Area(), jet2->Pt() - fRhoVal * jet2->Area());
1b76c28f 966 }
967 }
968
cd6431de 969 const Int_t nJets1 = fJets->GetEntriesFast();
1b76c28f 970
cd6431de 971 for (Int_t i = 0; i < nJets1; i++) {
1b76c28f 972
cd6431de 973 AliEmcalJet* jet1 = static_cast<AliEmcalJet*>(fJets->At(i));
1b76c28f 974
cd6431de 975 if (!jet1) {
976 AliError(Form("Could not receive jet %d", i));
2949a2ec 977 continue;
1b76c28f 978 }
99c67c1b 979
cd6431de 980 if (!AcceptJet(jet1))
1b76c28f 981 continue;
cd6431de 982
983 if (!AcceptBiasJet(jet1))
91bee8bc 984 continue;
cd6431de 985
986 if (jet1->MaxTrackPt() > fMaxTrackPt || jet1->MaxClusterPt() > fMaxClusterPt)
91bee8bc 987 continue;
cd6431de 988
989 if (jet1->Eta() < fJetMinEta || jet1->Eta() > fJetMaxEta || jet1->Phi() < fJetMinPhi || jet1->Phi() > fJetMaxPhi)
91bee8bc 990 continue;
2949a2ec 991
cd6431de 992 if (!jet1->MatchedJet()) {
993 fHistNonMatchedJets1PtArea->Fill(jet1->Area(), jet1->Pt());
994 if (!fRhoName.IsNull())
995 fHistNonMatchedJets1CorrPtArea->Fill(jet1->Area(), jet1->Pt() - fRhoVal * jet1->Area());
99c67c1b 996 }
1b76c28f 997
cd6431de 998 fHistJets1PtArea->Fill(jet1->Area(), jet1->Pt());
999 fHistJets1PhiEta->Fill(jet1->Eta(), jet1->Phi());
1b76c28f 1000
cd6431de 1001 if (!fRhoName.IsNull())
1002 fHistJets1CorrPtArea->Fill(jet1->Area(), jet1->Pt() - fRhoVal * jet1->Area());
1b76c28f 1003 }
6fd5039f 1004
1005 return kTRUE;
1b76c28f 1006}