]> 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
preparations for embedding (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"
b16bb001 19#include "AliLog.h"
2949a2ec 20#include "AliMCEvent.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),
37 fMatchingPar(0),
38 fJet2MinEta(-999),
39 fJet2MaxEta(-999),
40 fJet2MinPhi(-999),
41 fJet2MaxPhi(-999),
4643d2e8 42 fSelectPtHardBin(-999),
1ee1b5b8 43 fPythiaHeader(0),
1ee1b5b8 44 fPtHardBin(0),
45 fNTrials(0),
cd6431de 46 fTracks2(0),
47 fCaloClusters2(0),
48 fJets2(0),
49 fRho2(0),
50 fRho2Val(0),
1ee1b5b8 51 fHistNTrials(0),
1ee1b5b8 52 fHistEvents(0),
cd6431de 53 fHistJets1PhiEta(0),
54 fHistJets1PtArea(0),
55 fHistJets1CorrPtArea(0),
56 fHistJets2PhiEta(0),
57 fHistJets2PtArea(0),
58 fHistJets2CorrPtArea(0),
59 fHistMatchingLevelvsJet2Pt(0),
60 fHistClosestDeltaEtaPhivsJet2Pt(0),
61 fHistClosestDeltaPtvsJet2Pt(0),
62 fHistClosestDeltaCorrPtvsJet2Pt(0),
63 fHistNonMatchedJets1PtArea(0),
64 fHistNonMatchedJets2PtArea(0),
65 fHistNonMatchedJets1CorrPtArea(0),
66 fHistNonMatchedJets2CorrPtArea(0),
67 fHistJet1PtvsJet2Pt(0),
68 fHistJet1CorrPtvsJet2CorrPt(0),
69 fHistMissedJets2PtArea(0)
2949a2ec 70{
71 // Default constructor.
4643d2e8 72
a487deae 73 SetMakeGeneralHistograms(kTRUE);
2949a2ec 74}
75
76//________________________________________________________________________
77AliJetResponseMaker::AliJetResponseMaker(const char *name) :
e44e8726 78 AliAnalysisTaskEmcalJet(name, kTRUE),
cd6431de 79 fTracks2Name("MCParticles"),
80 fCalo2Name(""),
81 fJets2Name("MCJets"),
82 fRho2Name(""),
83 fPtBiasJet2Track(0),
84 fPtBiasJet2Clus(0),
85 fAreCollections1MC(kFALSE),
86 fAreCollections2MC(kTRUE),
87 fMatching(kNoMatching),
88 fMatchingPar(0.25),
89 fJet2MinEta(-999),
90 fJet2MaxEta(-999),
91 fJet2MinPhi(-999),
92 fJet2MaxPhi(-999),
4643d2e8 93 fSelectPtHardBin(-999),
1ee1b5b8 94 fPythiaHeader(0),
1ee1b5b8 95 fPtHardBin(0),
96 fNTrials(0),
cd6431de 97 fTracks2(0),
98 fCaloClusters2(0),
99 fJets2(0),
100 fRho2(0),
101 fRho2Val(0),
1ee1b5b8 102 fHistNTrials(0),
1ee1b5b8 103 fHistEvents(0),
cd6431de 104 fHistJets1PhiEta(0),
105 fHistJets1PtArea(0),
106 fHistJets1CorrPtArea(0),
107 fHistJets2PhiEta(0),
108 fHistJets2PtArea(0),
109 fHistJets2CorrPtArea(0),
110 fHistMatchingLevelvsJet2Pt(0),
111 fHistClosestDeltaEtaPhivsJet2Pt(0),
112 fHistClosestDeltaPtvsJet2Pt(0),
113 fHistClosestDeltaCorrPtvsJet2Pt(0),
114 fHistNonMatchedJets1PtArea(0),
115 fHistNonMatchedJets2PtArea(0),
116 fHistNonMatchedJets1CorrPtArea(0),
117 fHistNonMatchedJets2CorrPtArea(0),
118 fHistJet1PtvsJet2Pt(0),
119 fHistJet1CorrPtvsJet2CorrPt(0),
120 fHistMissedJets2PtArea(0)
2949a2ec 121{
122 // Standard constructor.
4643d2e8 123
a487deae 124 SetMakeGeneralHistograms(kTRUE);
2949a2ec 125}
126
127//________________________________________________________________________
128AliJetResponseMaker::~AliJetResponseMaker()
129{
130 // Destructor
131}
132
133//________________________________________________________________________
134void AliJetResponseMaker::UserCreateOutputObjects()
135{
136 // Create user objects.
137
a487deae 138 AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
2949a2ec 139
1ee1b5b8 140 const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
141 const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
142
143 fHistNTrials = new TH1F("fHistNTrials", "fHistNTrials", 11, 0, 11);
144 fHistNTrials->GetXaxis()->SetTitle("p_{T} hard bin");
145 fHistNTrials->GetYaxis()->SetTitle("trials");
146 fOutput->Add(fHistNTrials);
147
1ee1b5b8 148 fHistEvents = new TH1F("fHistEvents", "fHistEvents", 11, 0, 11);
149 fHistEvents->GetXaxis()->SetTitle("p_{T} hard bin");
150 fHistEvents->GetYaxis()->SetTitle("total events");
151 fOutput->Add(fHistEvents);
152
153 for (Int_t i = 1; i < 12; i++) {
154 fHistNTrials->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
1ee1b5b8 155 fHistEvents->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
156 }
157
cd6431de 158 fHistJets1PhiEta = new TH2F("fHistJets1PhiEta", "fHistJets1PhiEta", 20, -2, 2, 32, 0, 6.4);
159 fHistJets1PhiEta->GetXaxis()->SetTitle("#eta");
160 fHistJets1PhiEta->GetYaxis()->SetTitle("#phi");
161 fOutput->Add(fHistJets1PhiEta);
2949a2ec 162
cd6431de 163 fHistJets1PtArea = new TH2F("fHistJets1PtArea", "fHistJets1PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
164 fHistJets1PtArea->GetXaxis()->SetTitle("area");
165 fHistJets1PtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
166 fHistJets1PtArea->GetZaxis()->SetTitle("counts");
167 fOutput->Add(fHistJets1PtArea);
168
169 if (!fRhoName.IsNull()) {
170 fHistJets1CorrPtArea = new TH2F("fHistJets1CorrPtArea", "fHistJets1CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
171 fHistJets1CorrPtArea->GetXaxis()->SetTitle("area");
172 fHistJets1CorrPtArea->GetYaxis()->SetTitle("p_{T,1}^{corr} (GeV/c)");
173 fHistJets1CorrPtArea->GetZaxis()->SetTitle("counts");
174 fOutput->Add(fHistJets1CorrPtArea);
175 }
176
177 fHistJets2PhiEta = new TH2F("fHistJets2PhiEta", "fHistJets2PhiEta", 20, -2, 2, 32, 0, 6.4);
178 fHistJets2PhiEta->GetXaxis()->SetTitle("#eta");
179 fHistJets2PhiEta->GetYaxis()->SetTitle("#phi");
180 fOutput->Add(fHistJets2PhiEta);
7549c451 181
cd6431de 182 fHistJets2PtArea = new TH2F("fHistJets2PtArea", "fHistJets2PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
183 fHistJets2PtArea->GetXaxis()->SetTitle("area");
184 fHistJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
185 fHistJets2PtArea->GetZaxis()->SetTitle("counts");
186 fOutput->Add(fHistJets2PtArea);
187
188 if (!fRho2Name.IsNull()) {
189 fHistJets2CorrPtArea = new TH2F("fHistJets2CorrPtArea", "fHistJets2CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
190 fHistJets2CorrPtArea->GetXaxis()->SetTitle("area");
191 fHistJets2CorrPtArea->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)");
192 fHistJets2CorrPtArea->GetZaxis()->SetTitle("counts");
193 fOutput->Add(fHistJets2CorrPtArea);
194 }
195
196 fHistMatchingLevelvsJet2Pt = new TH2F("fHistMatchingLevelvsJet2Pt", "fHistMatchingLevelvsJet2Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt);
197 fHistMatchingLevelvsJet2Pt->GetXaxis()->SetTitle("Matching level");
198 fHistMatchingLevelvsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");
199 fHistMatchingLevelvsJet2Pt->GetZaxis()->SetTitle("counts");
200 fOutput->Add(fHistMatchingLevelvsJet2Pt);
201
202 fHistClosestDeltaEtaPhivsJet2Pt = new TH3F("fHistClosestDeltaEtaPhivsJet2Pt", "fHistClosestDeltaEtaPhivsJet2Pt", 40, -1, 1, 128, -1.6, 4.8, fNbins/2, fMinBinPt, fMaxBinPt);
203 fHistClosestDeltaEtaPhivsJet2Pt->GetXaxis()->SetTitle("#Delta#eta");
204 fHistClosestDeltaEtaPhivsJet2Pt->GetYaxis()->SetTitle("#Delta#phi");
205 fHistClosestDeltaEtaPhivsJet2Pt->GetZaxis()->SetTitle("p_{T,2}");
206 fOutput->Add(fHistClosestDeltaEtaPhivsJet2Pt);
207
208 fHistClosestDeltaPtvsJet2Pt = new TH2F("fHistClosestDeltaPtvsJet2Pt", "fHistClosestDeltaPtvsJet2Pt", fNbins/2, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
209 fHistClosestDeltaPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");
210 fHistClosestDeltaPtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
211 fHistClosestDeltaPtvsJet2Pt->GetZaxis()->SetTitle("counts");
212 fOutput->Add(fHistClosestDeltaPtvsJet2Pt);
213
214 if (!fRhoName.IsNull() || !fRho2Name.IsNull()) {
215 fHistClosestDeltaCorrPtvsJet2Pt = new TH2F("fHistClosestDeltaCorrPtvsJet2Pt", "fHistClosestDeltaCorrPtvsJet2Pt", fNbins/2, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
216 fHistClosestDeltaCorrPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");
217 fHistClosestDeltaCorrPtvsJet2Pt->GetYaxis()->SetTitle("#Deltap_{T}^{corr} (GeV/c)");
218 fHistClosestDeltaCorrPtvsJet2Pt->GetZaxis()->SetTitle("counts");
219 fOutput->Add(fHistClosestDeltaCorrPtvsJet2Pt);
220 }
221
222 fHistNonMatchedJets1PtArea = new TH2F("fHistNonMatchedJets1PtArea", "fHistNonMatchedJets1PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
223 fHistNonMatchedJets1PtArea->GetXaxis()->SetTitle("area");
224 fHistNonMatchedJets1PtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
225 fHistNonMatchedJets1PtArea->GetZaxis()->SetTitle("counts");
226 fOutput->Add(fHistNonMatchedJets1PtArea);
227
228 fHistNonMatchedJets2PtArea = new TH2F("fHistNonMatchedJets2PtArea", "fHistNonMatchedJets2PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
229 fHistNonMatchedJets2PtArea->GetXaxis()->SetTitle("area");
230 fHistNonMatchedJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
231 fHistNonMatchedJets2PtArea->GetZaxis()->SetTitle("counts");
232 fOutput->Add(fHistNonMatchedJets2PtArea);
233
234 if (!fRhoName.IsNull()) {
235 fHistNonMatchedJets1CorrPtArea = new TH2F("fHistNonMatchedJets1CorrPtArea", "fHistNonMatchedJets1CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
236 fHistNonMatchedJets1CorrPtArea->GetXaxis()->SetTitle("area");
237 fHistNonMatchedJets1CorrPtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
238 fHistNonMatchedJets1CorrPtArea->GetZaxis()->SetTitle("counts");
239 fOutput->Add(fHistNonMatchedJets1CorrPtArea);
240 }
241
242 if (!fRho2Name.IsNull()) {
243 fHistNonMatchedJets2CorrPtArea = new TH2F("fHistNonMatchedJets2CorrPtArea", "fHistNonMatchedJets2CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
244 fHistNonMatchedJets2CorrPtArea->GetXaxis()->SetTitle("area");
245 fHistNonMatchedJets2CorrPtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
246 fHistNonMatchedJets2CorrPtArea->GetZaxis()->SetTitle("counts");
247 fOutput->Add(fHistNonMatchedJets2CorrPtArea);
248 }
249
250 fHistJet1PtvsJet2Pt = new TH2F("fHistJet1PtvsJet2Pt", "fHistJet1PtvsJet2Pt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
251 fHistJet1PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,1}");
252 fHistJet1PtvsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");
253 fHistJet1PtvsJet2Pt->GetZaxis()->SetTitle("counts");
254 fOutput->Add(fHistJet1PtvsJet2Pt);
255
256 if (!fRhoName.IsNull() || !fRho2Name.IsNull()) {
257 if (fRhoName.IsNull())
258 fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt", fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
259 else if (fRho2Name.IsNull())
260 fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt", 2*fNbins, -fMaxBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
261 else
262 fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt", 2*fNbins, -fMaxBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
263 fHistJet1CorrPtvsJet2CorrPt->GetXaxis()->SetTitle("p_{T,1}^{corr}");
264 fHistJet1CorrPtvsJet2CorrPt->GetYaxis()->SetTitle("p_{T,2}^{corr}");
265 fHistJet1CorrPtvsJet2CorrPt->GetZaxis()->SetTitle("counts");
266 fOutput->Add(fHistJet1CorrPtvsJet2CorrPt);
267 }
268
269 fHistMissedJets2PtArea = new TH2F("fHistMissedJets2PtArea", "fHistMissedJets2PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
270 fHistMissedJets2PtArea->GetXaxis()->SetTitle("area");
271 fHistMissedJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
272 fHistMissedJets2PtArea->GetZaxis()->SetTitle("counts");
273 fOutput->Add(fHistMissedJets2PtArea);
2bddb6ae 274
99c67c1b 275 PostData(1, fOutput); // Post data for ALL output slots > 0 here, to get at least an empty histogram
2949a2ec 276}
277
ceefbfbc 278//________________________________________________________________________
a487deae 279Bool_t AliJetResponseMaker::AcceptJet(AliEmcalJet *jet) const
ceefbfbc 280{
281 // Return true if jet is accepted.
282
283 if (jet->Pt() <= fJetPtCut)
284 return kFALSE;
285 if (jet->Area() <= fJetAreaCut)
286 return kFALSE;
65bb5510 287
288 return kTRUE;
ceefbfbc 289}
290
cd6431de 291//________________________________________________________________________
292Bool_t AliJetResponseMaker::AcceptBiasJet2(AliEmcalJet *jet) const
293{
294 // Accept jet with a bias.
295
296 if (fLeadingHadronType == 0) {
297 if (jet->MaxTrackPt() < fPtBiasJet2Track) return kFALSE;
298 }
299 else if (fLeadingHadronType == 1) {
300 if (jet->MaxClusterPt() < fPtBiasJet2Clus) return kFALSE;
301 }
302 else {
303 if (jet->MaxTrackPt() < fPtBiasJet2Track && jet->MaxClusterPt() < fPtBiasJet2Clus) return kFALSE;
304 }
305
306 return kTRUE;
307}
308
b3ccdfe2 309//________________________________________________________________________
310void AliJetResponseMaker::ExecOnce()
311{
312 // Execute once.
313
cd6431de 314 if (!fJets2Name.IsNull() && !fJets2) {
315 fJets2 = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJets2Name));
316 if (!fJets2) {
317 AliError(Form("%s: Could not retrieve jets2 %s!", GetName(), fJets2Name.Data()));
b3ccdfe2 318 return;
319 }
cd6431de 320 else if (!fJets2->GetClass()->GetBaseClass("AliEmcalJet")) {
321 AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fJets2Name.Data()));
322 fJets2 = 0;
b3ccdfe2 323 return;
324 }
325 }
326
cd6431de 327 if (!fTracks2Name.IsNull() && !fTracks2) {
328 fTracks2 = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracks2Name));
329 if (!fTracks2) {
330 AliError(Form("%s: Could not retrieve tracks2 %s!", GetName(), fTracks2Name.Data()));
b3ccdfe2 331 return;
332 }
333 else {
cd6431de 334 TClass *cl = fTracks2->GetClass();
b3ccdfe2 335 if (!cl->GetBaseClass("AliVParticle") && !cl->GetBaseClass("AliEmcalParticle")) {
cd6431de 336 AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fTracks2Name.Data()));
337 fTracks2 = 0;
b3ccdfe2 338 return;
339 }
340 }
341 }
342
cd6431de 343 if (!fCalo2Name.IsNull() && !fCaloClusters2) {
344 fCaloClusters2 = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCalo2Name));
345 if (!fCaloClusters2) {
346 AliError(Form("%s: Could not retrieve clusters %s!", GetName(), fCalo2Name.Data()));
347 return;
348 } else {
349 TClass *cl = fCaloClusters2->GetClass();
350 if (!cl->GetBaseClass("AliVCluster") && !cl->GetBaseClass("AliEmcalParticle")) {
351 AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fCalo2Name.Data()));
352 fCaloClusters2 = 0;
353 return;
354 }
355 }
91bee8bc 356 }
cd6431de 357
358 if (!fRho2Name.IsNull() && !fRho2) {
359 fRho2 = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRho2Name));
360 if (!fRho2) {
361 AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRho2Name.Data()));
362 fInitialized = kFALSE;
363 return;
364 }
91bee8bc 365 }
cd6431de 366
367 if (fJet2MinEta == -999)
368 fJet2MinEta = fJetMinEta - fJetRadius;
369 if (fJet2MaxEta == -999)
370 fJet2MaxEta = fJetMaxEta + fJetRadius;
371 if (fJet2MinPhi == -999)
372 fJet2MinPhi = fJetMinPhi - fJetRadius;
373 if (fJet2MaxPhi == -999)
374 fJet2MaxPhi = fJetMaxPhi + fJetRadius;
375
376 AliAnalysisTaskEmcalJet::ExecOnce();
b3ccdfe2 377}
378
4643d2e8 379//________________________________________________________________________
380Bool_t AliJetResponseMaker::IsEventSelected()
381{
382 // Check if event is selected
383
384 if (fSelectPtHardBin != -999 && fSelectPtHardBin != fPtHardBin)
385 return kFALSE;
386
387 return AliAnalysisTaskEmcalJet::IsEventSelected();
388}
389
b3ccdfe2 390//________________________________________________________________________
391Bool_t AliJetResponseMaker::RetrieveEventObjects()
392{
393 // Retrieve event objects.
394
395 if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
396 return kFALSE;
cd6431de 397
398 if (fRho2)
399 fRho2Val = fRho2->GetVal();
b3ccdfe2 400
401 fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
402
403 if (!fPythiaHeader)
404 return kFALSE;
405
b3ccdfe2 406 const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
407 const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
408
409 Double_t pthard = fPythiaHeader->GetPtHard();
91bee8bc 410
9f52c61f 411 for (fPtHardBin = 0; fPtHardBin < 11; fPtHardBin++) {
b3ccdfe2 412 if (pthard >= ptHardLo[fPtHardBin] && pthard < ptHardHi[fPtHardBin])
413 break;
414 }
b3ccdfe2 415
416 fNTrials = fPythiaHeader->Trials();
417
418 return kTRUE;
419}
420
421//________________________________________________________________________
422Bool_t AliJetResponseMaker::Run()
423{
424 // Find the closest jets
cd6431de 425
426 if (fMatching == kNoMatching)
aa4d701c 427 return kTRUE;
428
cd6431de 429 DoJetLoop(fJets, fJets2, kFALSE);
430 DoJetLoop(fJets2, fJets, kTRUE);
b3ccdfe2 431
cd6431de 432 const Int_t nJets2 = fJets2->GetEntriesFast();
aa4d701c 433
cd6431de 434 for (Int_t i = 0; i < nJets2; i++) {
aa4d701c 435
cd6431de 436 AliEmcalJet* jet2 = static_cast<AliEmcalJet*>(fJets2->At(i));
aa4d701c 437
cd6431de 438 if (!jet2) {
aa4d701c 439 AliError(Form("Could not receive jet %d", i));
440 continue;
441 }
442
cd6431de 443 if (!AcceptJet(jet2))
aa4d701c 444 continue;
445
cd6431de 446 if (jet2->Eta() < fJet2MinEta || jet2->Eta() > fJet2MaxEta || jet2->Phi() < fJet2MinPhi || jet2->Phi() > fJet2MaxPhi)
aa4d701c 447 continue;
448
cd6431de 449 if (jet2->Pt() > fMaxBinPt)
aa4d701c 450 continue;
451
cd6431de 452 if (jet2->ClosestJet() && jet2->ClosestJet()->ClosestJet() == jet2 &&
453 jet2->ClosestJetDistance() < fMatchingPar) { // Matched jet found
454 jet2->SetMatchedToClosest();
455 jet2->ClosestJet()->SetMatchedToClosest();
aa4d701c 456 }
457 }
458
b3ccdfe2 459 return kTRUE;
460}
461
2949a2ec 462//________________________________________________________________________
cd6431de 463void AliJetResponseMaker::DoJetLoop(TClonesArray *jets1, TClonesArray *jets2, Bool_t first)
2949a2ec 464{
44476be7 465 // Do the jet loop.
1ee1b5b8 466
44476be7 467 Int_t nJets1 = jets1->GetEntriesFast();
468 Int_t nJets2 = jets2->GetEntriesFast();
1ee1b5b8 469
44476be7 470 for (Int_t i = 0; i < nJets1; i++) {
1ee1b5b8 471
44476be7 472 AliEmcalJet* jet1 = static_cast<AliEmcalJet*>(jets1->At(i));
2949a2ec 473
44476be7 474 if (!jet1) {
475 AliError(Form("Could not receive jet %d", i));
476 continue;
477 }
478
ceefbfbc 479 if (!AcceptJet(jet1))
44476be7 480 continue;
481
cd6431de 482 if (first) {
483 if (jet1->Eta() < fJet2MinEta || jet1->Eta() > fJet2MaxEta || jet1->Phi() < fJet2MinPhi || jet1->Phi() > fJet2MaxPhi)
ceefbfbc 484 continue;
485 }
486 else {
cd6431de 487 if (jet1->Eta() < fJetMinEta || jet1->Eta() > fJetMaxEta || jet1->Phi() < fJetMinPhi || jet1->Phi() > fJetMaxPhi)
ceefbfbc 488 continue;
489 }
490
44476be7 491 for (Int_t j = 0; j < nJets2; j++) {
492
493 AliEmcalJet* jet2 = static_cast<AliEmcalJet*>(jets2->At(j));
494
495 if (!jet2) {
496 AliError(Form("Could not receive jet %d", j));
497 continue;
498 }
499
ceefbfbc 500 if (!AcceptJet(jet2))
44476be7 501 continue;
ceefbfbc 502
cd6431de 503 if (first) {
a487deae 504 if (jet2->Eta() < fJetMinEta || jet2->Eta() > fJetMaxEta || jet2->Phi() < fJetMinPhi || jet2->Phi() > fJetMaxPhi)
ceefbfbc 505 continue;
506 }
507 else {
cd6431de 508 if (jet1->Eta() < fJet2MinEta || jet1->Eta() > fJet2MaxEta || jet1->Phi() < fJet2MinPhi || jet1->Phi() > fJet2MaxPhi)
ceefbfbc 509 continue;
510 }
cd6431de 511
512 Double_t d = GetMatchingLevel(jet1, jet2);
513
514 if (d < 0)
515 continue;
44476be7 516
517 if (d < jet1->ClosestJetDistance()) {
518 jet1->SetSecondClosestJet(jet1->ClosestJet(), jet1->ClosestJetDistance());
519 jet1->SetClosestJet(jet2, d);
520 }
521 else if (d < jet1->SecondClosestJetDistance()) {
522 jet1->SetSecondClosestJet(jet2, d);
523 }
524 }
525 }
1ee1b5b8 526}
527
cd6431de 528//________________________________________________________________________
529Double_t AliJetResponseMaker::GetMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2) const
530{
531 Double_t r = -1;
532
533 switch (fMatching) {
534 case kGeometrical:
535 {
536 Double_t deta = jet2->Eta() - jet1->Eta();
537 Double_t dphi = jet2->Phi() - jet1->Phi();
538 r = TMath::Sqrt(deta * deta + dphi * dphi);
539 }
540 break;
541 case kMCLabel:
542 AliError("MC label matching not implemented!");
543 break;
544 default:
545 ;
546 }
547
548 return r;
549}
550
1ee1b5b8 551//________________________________________________________________________
552Bool_t AliJetResponseMaker::FillHistograms()
553{
554 // Fill histograms.
555
4643d2e8 556 fHistEvents->SetBinContent(fPtHardBin + 1, fHistEvents->GetBinContent(fPtHardBin + 1) + 1);
557 fHistNTrials->SetBinContent(fPtHardBin + 1, fHistNTrials->GetBinContent(fPtHardBin + 1) + fNTrials);
4643d2e8 558
cd6431de 559 const Int_t nJets2 = fJets2->GetEntriesFast();
2949a2ec 560
cd6431de 561 for (Int_t i = 0; i < nJets2; i++) {
2949a2ec 562
cd6431de 563 AliEmcalJet* jet2 = static_cast<AliEmcalJet*>(fJets2->At(i));
2949a2ec 564
cd6431de 565 if (!jet2) {
566 AliError(Form("Could not receive jet2 %d", i));
2949a2ec 567 continue;
568 }
569
cd6431de 570 if (!AcceptJet(jet2))
571 continue;
572
573 if (!AcceptBiasJet2(jet2))
ceefbfbc 574 continue;
575
cd6431de 576 if (jet2->Eta() < fJet2MinEta || jet2->Eta() > fJet2MaxEta || jet2->Phi() < fJet2MinPhi || jet2->Phi() > fJet2MaxPhi)
2bddb6ae 577 continue;
578
cd6431de 579 if (jet2->Pt() > fMaxBinPt)
2949a2ec 580 continue;
581
cd6431de 582 if (jet2->MatchedJet()) {
aa4d701c 583
cd6431de 584 if (!AcceptBiasJet(jet2->MatchedJet()) ||
585 jet2->MatchedJet()->MaxTrackPt() > fMaxTrackPt || jet2->MatchedJet()->MaxClusterPt() > fMaxClusterPt ||
586 jet2->MatchedJet()->Pt() > fMaxBinPt) {
587 fHistMissedJets2PtArea->Fill(jet2->Area(), jet2->Pt());
2bddb6ae 588 }
589 else {
cd6431de 590 fHistMatchingLevelvsJet2Pt->Fill(jet2->ClosestJetDistance(), jet2->Pt());
591
592 Double_t deta = jet2->MatchedJet()->Eta() - jet2->Eta();
593 Double_t dphi = jet2->MatchedJet()->Phi() - jet2->Phi();
594 fHistClosestDeltaEtaPhivsJet2Pt->Fill(deta, dphi, jet2->Pt());
595
596 Double_t dpt = jet2->MatchedJet()->Pt() - jet2->Pt();
597 fHistClosestDeltaPtvsJet2Pt->Fill(jet2->Pt(), dpt);
598
599 fHistJet1PtvsJet2Pt->Fill(jet2->MatchedJet()->Pt(), jet2->Pt());
600
601 if (!fRhoName.IsNull() || !fRho2Name.IsNull()) {
602 dpt -= fRhoVal * jet2->MatchedJet()->Area() - fRho2Val * jet2->Area();
603 fHistClosestDeltaCorrPtvsJet2Pt->Fill(jet2->Pt(), dpt);
604 fHistJet1CorrPtvsJet2CorrPt->Fill(jet2->MatchedJet()->Pt() - fRhoVal * jet2->MatchedJet()->Area(), jet2->Pt() - fRho2Val * jet2->Area());
605 }
2bddb6ae 606 }
2949a2ec 607 }
99c67c1b 608 else {
cd6431de 609 fHistNonMatchedJets2PtArea->Fill(jet2->Area(), jet2->Pt());
610 fHistMissedJets2PtArea->Fill(jet2->Area(), jet2->Pt());
2949a2ec 611
cd6431de 612 if (!fRho2Name.IsNull())
613 fHistNonMatchedJets2CorrPtArea->Fill(jet2->Area(), jet2->Pt() - fRhoVal * jet2->Area());
1b76c28f 614 }
ceefbfbc 615
cd6431de 616 fHistJets2PtArea->Fill(jet2->Area(), jet2->Pt());
617 fHistJets2PhiEta->Fill(jet2->Eta(), jet2->Phi());
618
619 if (!fRho2Name.IsNull())
620 fHistJets2CorrPtArea->Fill(jet2->Area(), jet2->Pt() - fRho2Val * jet2->Area());
1b76c28f 621 }
622
cd6431de 623 const Int_t nJets1 = fJets->GetEntriesFast();
1b76c28f 624
cd6431de 625 for (Int_t i = 0; i < nJets1; i++) {
1b76c28f 626
cd6431de 627 AliEmcalJet* jet1 = static_cast<AliEmcalJet*>(fJets->At(i));
1b76c28f 628
cd6431de 629 if (!jet1) {
630 AliError(Form("Could not receive jet %d", i));
2949a2ec 631 continue;
1b76c28f 632 }
99c67c1b 633
cd6431de 634 if (!AcceptJet(jet1))
1b76c28f 635 continue;
cd6431de 636
637 if (!AcceptBiasJet(jet1))
91bee8bc 638 continue;
cd6431de 639
640 if (jet1->MaxTrackPt() > fMaxTrackPt || jet1->MaxClusterPt() > fMaxClusterPt)
91bee8bc 641 continue;
cd6431de 642
643 if (jet1->Eta() < fJetMinEta || jet1->Eta() > fJetMaxEta || jet1->Phi() < fJetMinPhi || jet1->Phi() > fJetMaxPhi)
91bee8bc 644 continue;
2949a2ec 645
cd6431de 646 if (!jet1->MatchedJet()) {
647 fHistNonMatchedJets1PtArea->Fill(jet1->Area(), jet1->Pt());
648 if (!fRhoName.IsNull())
649 fHistNonMatchedJets1CorrPtArea->Fill(jet1->Area(), jet1->Pt() - fRhoVal * jet1->Area());
99c67c1b 650 }
1b76c28f 651
cd6431de 652 fHistJets1PtArea->Fill(jet1->Area(), jet1->Pt());
653 fHistJets1PhiEta->Fill(jet1->Eta(), jet1->Phi());
1b76c28f 654
cd6431de 655 if (!fRhoName.IsNull())
656 fHistJets1CorrPtArea->Fill(jet1->Area(), jet1->Pt() - fRhoVal * jet1->Area());
1b76c28f 657 }
6fd5039f 658
659 return kTRUE;
1b76c28f 660}