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