]> 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
from Redmer
[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>
723f03e9 12#include <TProfile.h>
2949a2ec 13#include <TLorentzVector.h>
ca5c29fa 14#include <TSystem.h>
15#include <TFile.h>
723f03e9 16#include <TChain.h>
ca5c29fa 17#include <TKey.h>
18#include <TProfile.h>
2949a2ec 19
ca5c29fa 20#include "AliAnalysisManager.h"
2949a2ec 21#include "AliVCluster.h"
22#include "AliVTrack.h"
23#include "AliEmcalJet.h"
1ee1b5b8 24#include "AliGenPythiaEventHeader.h"
97f1b773 25#include "AliAODMCHeader.h"
2949a2ec 26#include "AliMCEvent.h"
787a3c4f 27#include "AliLog.h"
cd6431de 28#include "AliRhoParameter.h"
5be3857d 29#include "AliNamedArrayI.h"
2949a2ec 30
31ClassImp(AliJetResponseMaker)
32
33//________________________________________________________________________
34AliJetResponseMaker::AliJetResponseMaker() :
e44e8726 35 AliAnalysisTaskEmcalJet("AliJetResponseMaker", kTRUE),
cd6431de 36 fTracks2Name(""),
37 fCalo2Name(""),
38 fJets2Name(""),
39 fRho2Name(""),
40 fPtBiasJet2Track(0),
41 fPtBiasJet2Clus(0),
42 fAreCollections1MC(kFALSE),
43 fAreCollections2MC(kTRUE),
44 fMatching(kNoMatching),
7f76e479 45 fMatchingPar1(0),
46 fMatchingPar2(0),
cd6431de 47 fJet2MinEta(-999),
48 fJet2MaxEta(-999),
49 fJet2MinPhi(-999),
50 fJet2MaxPhi(-999),
4643d2e8 51 fSelectPtHardBin(-999),
4358e58a 52 fIsEmbedded(kFALSE),
53 fIsPythia(kTRUE),
1ee1b5b8 54 fPythiaHeader(0),
1ee1b5b8 55 fPtHardBin(0),
56 fNTrials(0),
cd6431de 57 fTracks2(0),
58 fCaloClusters2(0),
59 fJets2(0),
60 fRho2(0),
61 fRho2Val(0),
cfc2ac24 62 fTracks2Map(0),
ca5c29fa 63 fHistTrialsAfterSel(0),
64 fHistEventsAfterSel(0),
65 fHistTrials(0),
66 fHistXsection(0),
1ee1b5b8 67 fHistEvents(0),
fde82e42 68 fMCEnergy1vsEnergy2(0),
cd6431de 69 fHistJets1PhiEta(0),
70 fHistJets1PtArea(0),
71 fHistJets1CorrPtArea(0),
ca5c29fa 72 fHistLeadingJets1PtArea(0),
73 fHistLeadingJets1CorrPtArea(0),
63fac07f 74 fHistJets1NEFvsPt(0),
75 fHistJets1CEFvsCEFPt(0),
76 fHistJets1ZvsPt(0),
cd6431de 77 fHistJets2PhiEta(0),
78 fHistJets2PtArea(0),
79 fHistJets2CorrPtArea(0),
ca5c29fa 80 fHistLeadingJets2PtArea(0),
81 fHistLeadingJets2CorrPtArea(0),
cfc2ac24 82 fHistJets2PhiEtaAcceptance(0),
83 fHistJets2PtAreaAcceptance(0),
84 fHistJets2CorrPtAreaAcceptance(0),
ca5c29fa 85 fHistLeadingJets2PtAreaAcceptance(0),
86 fHistLeadingJets2CorrPtAreaAcceptance(0),
63fac07f 87 fHistJets2NEFvsPt(0),
88 fHistJets2CEFvsCEFPt(0),
89 fHistJets2ZvsPt(0),
05077f28 90 fHistNonMatchedJets1PtArea(0),
91 fHistNonMatchedJets2PtArea(0),
92 fHistNonMatchedJets1CorrPtArea(0),
93 fHistNonMatchedJets2CorrPtArea(0),
94 fHistMissedJets2PtArea(0),
ca5c29fa 95 fHistCommonEnergy1vsJet1Pt(0),
96 fHistCommonEnergy2vsJet2Pt(0),
97 fHistDistancevsJet1Pt(0),
98 fHistDistancevsJet2Pt(0),
6a20534a 99 fHistDistancevsCommonEnergy1(0),
100 fHistDistancevsCommonEnergy2(0),
05077f28 101 fHistCommonEnergy1vsCommonEnergy2(0),
c560b734 102 fHistJet2PtOverJet1PtvsJet2Pt(0),
7030f36f 103 fHistJet1PtOverJet2PtvsJet1Pt(0),
ca5c29fa 104 fHistDeltaEtaPhi(0),
7030f36f 105 fHistDeltaPtvsJet1Pt(0),
cfc2ac24 106 fHistDeltaPtvsJet2Pt(0),
05077f28 107 fHistDeltaPtvsDistance(0),
108 fHistDeltaPtvsCommonEnergy1(0),
109 fHistDeltaPtvsCommonEnergy2(0),
110 fHistDeltaPtvsArea1(0),
111 fHistDeltaPtvsArea2(0),
112 fHistDeltaPtvsDeltaArea(0),
113 fHistJet1PtvsJet2Pt(0),
7030f36f 114 fHistDeltaCorrPtvsJet1Pt(0),
cfc2ac24 115 fHistDeltaCorrPtvsJet2Pt(0),
05077f28 116 fHistDeltaCorrPtvsDistance(0),
117 fHistDeltaCorrPtvsCommonEnergy1(0),
118 fHistDeltaCorrPtvsCommonEnergy2(0),
119 fHistDeltaCorrPtvsArea1(0),
120 fHistDeltaCorrPtvsArea2(0),
121 fHistDeltaCorrPtvsDeltaArea(0),
cd6431de 122 fHistJet1CorrPtvsJet2CorrPt(0),
05077f28 123 fHistDeltaMCPtvsJet1Pt(0),
124 fHistDeltaMCPtvsJet2Pt(0),
125 fHistDeltaMCPtvsDistance(0),
126 fHistDeltaMCPtvsCommonEnergy1(0),
127 fHistDeltaMCPtvsCommonEnergy2(0),
128 fHistDeltaMCPtvsArea1(0),
129 fHistDeltaMCPtvsArea2(0),
130 fHistDeltaMCPtvsDeltaArea(0),
131 fHistJet1MCPtvsJet2Pt(0)
2949a2ec 132{
133 // Default constructor.
4643d2e8 134
a487deae 135 SetMakeGeneralHistograms(kTRUE);
2949a2ec 136}
137
138//________________________________________________________________________
139AliJetResponseMaker::AliJetResponseMaker(const char *name) :
e44e8726 140 AliAnalysisTaskEmcalJet(name, kTRUE),
cd6431de 141 fTracks2Name("MCParticles"),
142 fCalo2Name(""),
143 fJets2Name("MCJets"),
144 fRho2Name(""),
145 fPtBiasJet2Track(0),
146 fPtBiasJet2Clus(0),
147 fAreCollections1MC(kFALSE),
148 fAreCollections2MC(kTRUE),
149 fMatching(kNoMatching),
7f76e479 150 fMatchingPar1(0),
151 fMatchingPar2(0),
cd6431de 152 fJet2MinEta(-999),
153 fJet2MaxEta(-999),
154 fJet2MinPhi(-999),
155 fJet2MaxPhi(-999),
4643d2e8 156 fSelectPtHardBin(-999),
4358e58a 157 fIsEmbedded(kFALSE),
158 fIsPythia(kTRUE),
1ee1b5b8 159 fPythiaHeader(0),
1ee1b5b8 160 fPtHardBin(0),
161 fNTrials(0),
cd6431de 162 fTracks2(0),
163 fCaloClusters2(0),
164 fJets2(0),
165 fRho2(0),
166 fRho2Val(0),
cfc2ac24 167 fTracks2Map(0),
ca5c29fa 168 fHistTrialsAfterSel(0),
169 fHistEventsAfterSel(0),
170 fHistTrials(0),
171 fHistXsection(0),
1ee1b5b8 172 fHistEvents(0),
fde82e42 173 fMCEnergy1vsEnergy2(0),
cd6431de 174 fHistJets1PhiEta(0),
175 fHistJets1PtArea(0),
176 fHistJets1CorrPtArea(0),
ca5c29fa 177 fHistLeadingJets1PtArea(0),
178 fHistLeadingJets1CorrPtArea(0),
63fac07f 179 fHistJets1NEFvsPt(0),
180 fHistJets1CEFvsCEFPt(0),
181 fHistJets1ZvsPt(0),
cd6431de 182 fHistJets2PhiEta(0),
183 fHistJets2PtArea(0),
184 fHistJets2CorrPtArea(0),
ca5c29fa 185 fHistLeadingJets2PtArea(0),
186 fHistLeadingJets2CorrPtArea(0),
cfc2ac24 187 fHistJets2PhiEtaAcceptance(0),
188 fHistJets2PtAreaAcceptance(0),
189 fHistJets2CorrPtAreaAcceptance(0),
ca5c29fa 190 fHistLeadingJets2PtAreaAcceptance(0),
191 fHistLeadingJets2CorrPtAreaAcceptance(0),
63fac07f 192 fHistJets2NEFvsPt(0),
193 fHistJets2CEFvsCEFPt(0),
194 fHistJets2ZvsPt(0),
05077f28 195 fHistNonMatchedJets1PtArea(0),
196 fHistNonMatchedJets2PtArea(0),
197 fHistNonMatchedJets1CorrPtArea(0),
198 fHistNonMatchedJets2CorrPtArea(0),
199 fHistMissedJets2PtArea(0),
ca5c29fa 200 fHistCommonEnergy1vsJet1Pt(0),
201 fHistCommonEnergy2vsJet2Pt(0),
202 fHistDistancevsJet1Pt(0),
203 fHistDistancevsJet2Pt(0),
6a20534a 204 fHistDistancevsCommonEnergy1(0),
205 fHistDistancevsCommonEnergy2(0),
05077f28 206 fHistCommonEnergy1vsCommonEnergy2(0),
c560b734 207 fHistJet2PtOverJet1PtvsJet2Pt(0),
7030f36f 208 fHistJet1PtOverJet2PtvsJet1Pt(0),
ca5c29fa 209 fHistDeltaEtaPhi(0),
7030f36f 210 fHistDeltaPtvsJet1Pt(0),
cfc2ac24 211 fHistDeltaPtvsJet2Pt(0),
05077f28 212 fHistDeltaPtvsDistance(0),
213 fHistDeltaPtvsCommonEnergy1(0),
214 fHistDeltaPtvsCommonEnergy2(0),
215 fHistDeltaPtvsArea1(0),
216 fHistDeltaPtvsArea2(0),
217 fHistDeltaPtvsDeltaArea(0),
218 fHistJet1PtvsJet2Pt(0),
7030f36f 219 fHistDeltaCorrPtvsJet1Pt(0),
cfc2ac24 220 fHistDeltaCorrPtvsJet2Pt(0),
05077f28 221 fHistDeltaCorrPtvsDistance(0),
222 fHistDeltaCorrPtvsCommonEnergy1(0),
223 fHistDeltaCorrPtvsCommonEnergy2(0),
224 fHistDeltaCorrPtvsArea1(0),
225 fHistDeltaCorrPtvsArea2(0),
226 fHistDeltaCorrPtvsDeltaArea(0),
cd6431de 227 fHistJet1CorrPtvsJet2CorrPt(0),
05077f28 228 fHistDeltaMCPtvsJet1Pt(0),
229 fHistDeltaMCPtvsJet2Pt(0),
230 fHistDeltaMCPtvsDistance(0),
231 fHistDeltaMCPtvsCommonEnergy1(0),
232 fHistDeltaMCPtvsCommonEnergy2(0),
233 fHistDeltaMCPtvsArea1(0),
234 fHistDeltaMCPtvsArea2(0),
235 fHistDeltaMCPtvsDeltaArea(0),
236 fHistJet1MCPtvsJet2Pt(0)
2949a2ec 237{
238 // Standard constructor.
4643d2e8 239
a487deae 240 SetMakeGeneralHistograms(kTRUE);
2949a2ec 241}
242
243//________________________________________________________________________
244AliJetResponseMaker::~AliJetResponseMaker()
245{
246 // Destructor
247}
248
ca5c29fa 249//________________________________________________________________________
250Bool_t AliJetResponseMaker::PythiaInfoFromFile(const char* currFile, Float_t &fXsec, Float_t &fTrials, Int_t &pthard)
251{
252 //
253 // Get the cross section and the trails either from pyxsec.root or from pysec_hists.root
254 // Get the pt hard bin from the file path
255 // This is to called in Notify and should provide the path to the AOD/ESD file
256 // (Partially copied from AliAnalysisHelperJetTasks)
257
258 TString file(currFile);
259 fXsec = 0;
260 fTrials = 1;
261
2103dc6a 262 if(file.Contains(".zip#")){
ca5c29fa 263 Ssiz_t pos1 = file.Index("root_archive",12,0,TString::kExact);
264 Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
265 Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
266 file.Replace(pos+1,pos2-pos1,"");
267 }
268 else {
269 // not an archive take the basename....
270 file.ReplaceAll(gSystem->BaseName(file.Data()),"");
271 }
272 Printf("%s",file.Data());
273
274 // Get the pt hard bin
275 TString strPthard(file);
276 strPthard.Remove(strPthard.Last('/'));
277 strPthard.Remove(strPthard.Last('/'));
278 strPthard.Remove(0,strPthard.Last('/')+1);
279 if (strPthard.IsDec())
280 pthard = strPthard.Atoi();
281 else
282 AliWarning(Form("Could not extract file number from path %s", strPthard.Data()));
283
284 TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root")); // problem that we cannot really test the existance of a file in a archive so we have to lvie with open error message from root
285 if(!fxsec){
286 // next trial fetch the histgram file
287 fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
288 if(!fxsec){
289 // not a severe condition but inciate that we have no information
290 return kFALSE;
291 }
292 else{
293 // find the tlist we want to be independtent of the name so use the Tkey
294 TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
295 if(!key){
296 fxsec->Close();
297 return kFALSE;
298 }
299 TList *list = dynamic_cast<TList*>(key->ReadObj());
300 if(!list){
301 fxsec->Close();
302 return kFALSE;
303 }
304 fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
305 fTrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
306 fxsec->Close();
307 }
308 } // no tree pyxsec.root
309 else {
310 TTree *xtree = (TTree*)fxsec->Get("Xsection");
311 if(!xtree){
312 fxsec->Close();
313 return kFALSE;
314 }
315 UInt_t ntrials = 0;
316 Double_t xsection = 0;
317 xtree->SetBranchAddress("xsection",&xsection);
318 xtree->SetBranchAddress("ntrials",&ntrials);
319 xtree->GetEntry(0);
320 fTrials = ntrials;
321 fXsec = xsection;
322 fxsec->Close();
323 }
324 return kTRUE;
325}
326
327//________________________________________________________________________
328Bool_t AliJetResponseMaker::UserNotify()
329{
4358e58a 330 if (!fIsPythia)
331 return kTRUE;
332
ca5c29fa 333 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
334 if (!tree) {
335 AliError(Form("%s - UserNotify: No current tree!",GetName()));
336 return kFALSE;
337 }
338
339 Float_t xsection = 0;
340 Float_t trials = 0;
341 Int_t pthard = 0;
342
343 TFile *curfile = tree->GetCurrentFile();
344 if (!curfile) {
345 AliError(Form("%s - UserNotify: No current file!",GetName()));
346 return kFALSE;
347 }
348
723f03e9 349 TChain *chain = dynamic_cast<TChain*>(tree);
350 if (chain)
351 tree = chain->GetTree();
352
353 Int_t nevents = tree->GetEntriesFast();
354
ca5c29fa 355 PythiaInfoFromFile(curfile->GetName(), xsection, trials, pthard);
356
723f03e9 357 fHistTrials->Fill(pthard, trials);
358 fHistXsection->Fill(pthard, xsection);
359 fHistEvents->Fill(pthard, nevents);
ca5c29fa 360
361 return kTRUE;
362}
363
2949a2ec 364//________________________________________________________________________
365void AliJetResponseMaker::UserCreateOutputObjects()
366{
367 // Create user objects.
368
a487deae 369 AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
2949a2ec 370
05077f28 371 // General histograms
372 if (fIsPythia) {
373 fHistTrialsAfterSel = new TH1F("fHistTrialsAfterSel", "fHistTrialsAfterSel", 11, 0, 11);
374 fHistTrialsAfterSel->GetXaxis()->SetTitle("p_{T} hard bin");
375 fHistTrialsAfterSel->GetYaxis()->SetTitle("trials");
376 fOutput->Add(fHistTrialsAfterSel);
377
378 fHistEventsAfterSel = new TH1F("fHistEventsAfterSel", "fHistEventsAfterSel", 11, 0, 11);
379 fHistEventsAfterSel->GetXaxis()->SetTitle("p_{T} hard bin");
380 fHistEventsAfterSel->GetYaxis()->SetTitle("total events");
381 fOutput->Add(fHistEventsAfterSel);
382
383 fHistTrials = new TH1F("fHistTrials", "fHistTrials", 11, 0, 11);
384 fHistTrials->GetXaxis()->SetTitle("p_{T} hard bin");
385 fHistTrials->GetYaxis()->SetTitle("trials");
386 fOutput->Add(fHistTrials);
387
388 fHistXsection = new TProfile("fHistXsection", "fHistXsection", 11, 0, 11);
389 fHistXsection->GetXaxis()->SetTitle("p_{T} hard bin");
390 fHistXsection->GetYaxis()->SetTitle("xsection");
391 fOutput->Add(fHistXsection);
392
393 fHistEvents = new TH1F("fHistEvents", "fHistEvents", 11, 0, 11);
394 fHistEvents->GetXaxis()->SetTitle("p_{T} hard bin");
395 fHistEvents->GetYaxis()->SetTitle("total events");
396 fOutput->Add(fHistEvents);
397
398 const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
399 const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
400
401 for (Int_t i = 1; i < 12; i++) {
402 fHistTrialsAfterSel->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
403 fHistEventsAfterSel->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
404
405 fHistTrials->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
406 fHistXsection->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
407 fHistEvents->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
408 }
1ee1b5b8 409 }
410
fde82e42 411 if (fIsEmbedded) {
412 fMCEnergy1vsEnergy2 = new TH2F("fMCEnergy1vsEnergy2", "fMCEnergy1vsEnergy2", fNbins, fMinBinPt, fMaxBinPt*4, fNbins, fMinBinPt, fMaxBinPt*4);
413 fMCEnergy1vsEnergy2->GetXaxis()->SetTitle("#Sigmap_{T,1}^{MC}");
414 fMCEnergy1vsEnergy2->GetYaxis()->SetTitle("#Sigmap_{T,2}");
415 fMCEnergy1vsEnergy2->GetZaxis()->SetTitle("counts");
416 fOutput->Add(fMCEnergy1vsEnergy2);
417 }
418
05077f28 419 // Jets 1 histograms
420 fHistLeadingJets1PtArea = new TH2F("fHistLeadingJets1PtArea", "fHistLeadingJets1PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
421 fHistLeadingJets1PtArea->GetXaxis()->SetTitle("area");
422 fHistLeadingJets1PtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
423 fHistLeadingJets1PtArea->GetZaxis()->SetTitle("counts");
424 fOutput->Add(fHistLeadingJets1PtArea);
425
cfc2ac24 426 fHistJets1PhiEta = new TH2F("fHistJets1PhiEta", "fHistJets1PhiEta", 40, -1, 1, 40, 0, TMath::Pi()*2);
cd6431de 427 fHistJets1PhiEta->GetXaxis()->SetTitle("#eta");
428 fHistJets1PhiEta->GetYaxis()->SetTitle("#phi");
429 fOutput->Add(fHistJets1PhiEta);
2949a2ec 430
cd6431de 431 fHistJets1PtArea = new TH2F("fHistJets1PtArea", "fHistJets1PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
432 fHistJets1PtArea->GetXaxis()->SetTitle("area");
433 fHistJets1PtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
434 fHistJets1PtArea->GetZaxis()->SetTitle("counts");
435 fOutput->Add(fHistJets1PtArea);
436
437 if (!fRhoName.IsNull()) {
ca5c29fa 438 fHistLeadingJets1CorrPtArea = new TH2F("fHistLeadingJets1CorrPtArea", "fHistLeadingJets1CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
439 fHistLeadingJets1CorrPtArea->GetXaxis()->SetTitle("area");
440 fHistLeadingJets1CorrPtArea->GetYaxis()->SetTitle("p_{T,1}^{corr} (GeV/c)");
441 fHistLeadingJets1CorrPtArea->GetZaxis()->SetTitle("counts");
442 fOutput->Add(fHistLeadingJets1CorrPtArea);
05077f28 443
444 fHistJets1CorrPtArea = new TH2F("fHistJets1CorrPtArea", "fHistJets1CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
445 fHistJets1CorrPtArea->GetXaxis()->SetTitle("area");
446 fHistJets1CorrPtArea->GetYaxis()->SetTitle("p_{T,1}^{corr} (GeV/c)");
447 fHistJets1CorrPtArea->GetZaxis()->SetTitle("counts");
448 fOutput->Add(fHistJets1CorrPtArea);
cd6431de 449 }
450
63fac07f 451 fHistJets1ZvsPt = new TH2F("fHistJets1ZvsPt", "fHistJets1ZvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
452 fHistJets1ZvsPt->GetXaxis()->SetTitle("Z");
453 fHistJets1ZvsPt->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
454 fHistJets1ZvsPt->GetZaxis()->SetTitle("counts");
455 fOutput->Add(fHistJets1ZvsPt);
456
457 fHistJets1NEFvsPt = new TH2F("fHistJets1NEFvsPt", "fHistJets1NEFvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
458 fHistJets1NEFvsPt->GetXaxis()->SetTitle("NEF");
459 fHistJets1NEFvsPt->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
460 fHistJets1NEFvsPt->GetZaxis()->SetTitle("counts");
461 fOutput->Add(fHistJets1NEFvsPt);
462
463 fHistJets1CEFvsCEFPt = new TH2F("fHistJets1CEFvsCEFPt", "fHistJets1CEFvsCEFPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
464 fHistJets1CEFvsCEFPt->GetXaxis()->SetTitle("1-NEF");
465 fHistJets1CEFvsCEFPt->GetYaxis()->SetTitle("(1-NEF)*p_{T,1} (GeV/c)");
466 fHistJets1CEFvsCEFPt->GetZaxis()->SetTitle("counts");
467 fOutput->Add(fHistJets1CEFvsCEFPt);
468
05077f28 469 // Jets 2 histograms
cfc2ac24 470
05077f28 471 fHistLeadingJets2PtAreaAcceptance = new TH2F("fHistLeadingJets2PtAreaAcceptance", "fHistLeadingJets2PtAreaAcceptance", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
472 fHistLeadingJets2PtAreaAcceptance->GetXaxis()->SetTitle("area");
473 fHistLeadingJets2PtAreaAcceptance->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
474 fHistLeadingJets2PtAreaAcceptance->GetZaxis()->SetTitle("counts");
475 fOutput->Add(fHistLeadingJets2PtAreaAcceptance);
cd6431de 476
ca5c29fa 477 fHistLeadingJets2PtArea = new TH2F("fHistLeadingJets2PtArea", "fHistLeadingJets2PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
478 fHistLeadingJets2PtArea->GetXaxis()->SetTitle("area");
479 fHistLeadingJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
480 fHistLeadingJets2PtArea->GetZaxis()->SetTitle("counts");
481 fOutput->Add(fHistLeadingJets2PtArea);
482
05077f28 483 fHistJets2PhiEtaAcceptance = new TH2F("fHistJets2PhiEtaAcceptance", "fHistJets2PhiEtaAcceptance", 40, -1, 1, 40, 0, TMath::Pi()*2);
484 fHistJets2PhiEtaAcceptance->GetXaxis()->SetTitle("#eta");
485 fHistJets2PhiEtaAcceptance->GetYaxis()->SetTitle("#phi");
486 fOutput->Add(fHistJets2PhiEtaAcceptance);
487
cfc2ac24 488 fHistJets2PtAreaAcceptance = new TH2F("fHistJets2PtAreaAcceptance", "fHistJets2PtAreaAcceptance", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
489 fHistJets2PtAreaAcceptance->GetXaxis()->SetTitle("area");
490 fHistJets2PtAreaAcceptance->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
491 fHistJets2PtAreaAcceptance->GetZaxis()->SetTitle("counts");
492 fOutput->Add(fHistJets2PtAreaAcceptance);
493
05077f28 494 fHistJets2PhiEta = new TH2F("fHistJets2PhiEta", "fHistJets2PhiEta", 40, -1, 1, 40, 0, TMath::Pi()*2);
495 fHistJets2PhiEta->GetXaxis()->SetTitle("#eta");
496 fHistJets2PhiEta->GetYaxis()->SetTitle("#phi");
497 fOutput->Add(fHistJets2PhiEta);
cfc2ac24 498
05077f28 499 fHistJets2PtArea = new TH2F("fHistJets2PtArea", "fHistJets2PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
500 fHistJets2PtArea->GetXaxis()->SetTitle("area");
501 fHistJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
502 fHistJets2PtArea->GetZaxis()->SetTitle("counts");
503 fOutput->Add(fHistJets2PtArea);
ca5c29fa 504
05077f28 505 if (!fRho2Name.IsNull()) {
cfc2ac24 506 fHistJets2CorrPtAreaAcceptance = new TH2F("fHistJets2CorrPtAreaAcceptance", "fHistJets2CorrPtAreaAcceptance", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
507 fHistJets2CorrPtAreaAcceptance->GetXaxis()->SetTitle("area");
508 fHistJets2CorrPtAreaAcceptance->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)");
509 fHistJets2CorrPtAreaAcceptance->GetZaxis()->SetTitle("counts");
510 fOutput->Add(fHistJets2CorrPtAreaAcceptance);
ca5c29fa 511
512 fHistLeadingJets2CorrPtAreaAcceptance = new TH2F("fHistLeadingJets2CorrPtAreaAcceptance", "fHistLeadingJets2CorrPtAreaAcceptance", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
513 fHistLeadingJets2CorrPtAreaAcceptance->GetXaxis()->SetTitle("area");
514 fHistLeadingJets2CorrPtAreaAcceptance->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)");
515 fHistLeadingJets2CorrPtAreaAcceptance->GetZaxis()->SetTitle("counts");
516 fOutput->Add(fHistLeadingJets2CorrPtAreaAcceptance);
05077f28 517
518 fHistLeadingJets2CorrPtArea = new TH2F("fHistLeadingJets2CorrPtArea", "fHistLeadingJets2CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
519 fHistLeadingJets2CorrPtArea->GetXaxis()->SetTitle("area");
520 fHistLeadingJets2CorrPtArea->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)");
521 fHistLeadingJets2CorrPtArea->GetZaxis()->SetTitle("counts");
522 fOutput->Add(fHistLeadingJets2CorrPtArea);
523
524 fHistJets2CorrPtArea = new TH2F("fHistJets2CorrPtArea", "fHistJets2CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
525 fHistJets2CorrPtArea->GetXaxis()->SetTitle("area");
526 fHistJets2CorrPtArea->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)");
527 fHistJets2CorrPtArea->GetZaxis()->SetTitle("counts");
528 fOutput->Add(fHistJets2CorrPtArea);
cd6431de 529 }
530
63fac07f 531 fHistJets2ZvsPt = new TH2F("fHistJets2ZvsPt", "fHistJets2ZvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
532 fHistJets2ZvsPt->GetXaxis()->SetTitle("Z");
533 fHistJets2ZvsPt->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
534 fHistJets2ZvsPt->GetZaxis()->SetTitle("counts");
535 fOutput->Add(fHistJets2ZvsPt);
536
537 fHistJets2NEFvsPt = new TH2F("fHistJets2NEFvsPt", "fHistJets2NEFvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
538 fHistJets2NEFvsPt->GetXaxis()->SetTitle("NEF");
539 fHistJets2NEFvsPt->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
540 fHistJets2NEFvsPt->GetZaxis()->SetTitle("counts");
541 fOutput->Add(fHistJets2NEFvsPt);
542
543 fHistJets2CEFvsCEFPt = new TH2F("fHistJets2CEFvsCEFPt", "fHistJets2CEFvsCEFPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
544 fHistJets2CEFvsCEFPt->GetXaxis()->SetTitle("1-NEF");
545 fHistJets2CEFvsCEFPt->GetYaxis()->SetTitle("(1-NEF)*p_{T,2} (GeV/c)");
546 fHistJets2CEFvsCEFPt->GetZaxis()->SetTitle("counts");
547 fOutput->Add(fHistJets2CEFvsCEFPt);
548
05077f28 549 // Matching histograms
550
551 fHistNonMatchedJets1PtArea = new TH2F("fHistNonMatchedJets1PtArea", "fHistNonMatchedJets1PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
552 fHistNonMatchedJets1PtArea->GetXaxis()->SetTitle("area");
553 fHistNonMatchedJets1PtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
554 fHistNonMatchedJets1PtArea->GetZaxis()->SetTitle("counts");
555 fOutput->Add(fHistNonMatchedJets1PtArea);
556
557 fHistNonMatchedJets2PtArea = new TH2F("fHistNonMatchedJets2PtArea", "fHistNonMatchedJets2PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
558 fHistNonMatchedJets2PtArea->GetXaxis()->SetTitle("area");
559 fHistNonMatchedJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
560 fHistNonMatchedJets2PtArea->GetZaxis()->SetTitle("counts");
561 fOutput->Add(fHistNonMatchedJets2PtArea);
562
563 if (!fRhoName.IsNull()) {
564 fHistNonMatchedJets1CorrPtArea = new TH2F("fHistNonMatchedJets1CorrPtArea", "fHistNonMatchedJets1CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
565 fHistNonMatchedJets1CorrPtArea->GetXaxis()->SetTitle("area");
566 fHistNonMatchedJets1CorrPtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
567 fHistNonMatchedJets1CorrPtArea->GetZaxis()->SetTitle("counts");
568 fOutput->Add(fHistNonMatchedJets1CorrPtArea);
569 }
570
571 if (!fRho2Name.IsNull()) {
572 fHistNonMatchedJets2CorrPtArea = new TH2F("fHistNonMatchedJets2CorrPtArea", "fHistNonMatchedJets2CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
573 fHistNonMatchedJets2CorrPtArea->GetXaxis()->SetTitle("area");
574 fHistNonMatchedJets2CorrPtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
575 fHistNonMatchedJets2CorrPtArea->GetZaxis()->SetTitle("counts");
576 fOutput->Add(fHistNonMatchedJets2CorrPtArea);
577 }
578
579 fHistMissedJets2PtArea = new TH2F("fHistMissedJets2PtArea", "fHistMissedJets2PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
580 fHistMissedJets2PtArea->GetXaxis()->SetTitle("area");
581 fHistMissedJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
582 fHistMissedJets2PtArea->GetZaxis()->SetTitle("counts");
583 fOutput->Add(fHistMissedJets2PtArea);
584
ca5c29fa 585 fHistCommonEnergy1vsJet1Pt = new TH2F("fHistCommonEnergy1vsJet1Pt", "fHistCommonEnergy1vsJet1Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt);
586 fHistCommonEnergy1vsJet1Pt->GetXaxis()->SetTitle("Common energy 1 (%)");
587 fHistCommonEnergy1vsJet1Pt->GetYaxis()->SetTitle("p_{T,1}");
588 fHistCommonEnergy1vsJet1Pt->GetZaxis()->SetTitle("counts");
589 fOutput->Add(fHistCommonEnergy1vsJet1Pt);
590
591 fHistCommonEnergy2vsJet2Pt = new TH2F("fHistCommonEnergy2vsJet2Pt", "fHistCommonEnergy2vsJet2Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt);
592 fHistCommonEnergy2vsJet2Pt->GetXaxis()->SetTitle("Common energy 2 (%)");
593 fHistCommonEnergy2vsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");
594 fHistCommonEnergy2vsJet2Pt->GetZaxis()->SetTitle("counts");
595 fOutput->Add(fHistCommonEnergy2vsJet2Pt);
596
597 fHistDistancevsJet1Pt = new TH2F("fHistDistancevsJet1Pt", "fHistDistancevsJet1Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt);
598 fHistDistancevsJet1Pt->GetXaxis()->SetTitle("Distance");
599 fHistDistancevsJet1Pt->GetYaxis()->SetTitle("p_{T,1}");
600 fHistDistancevsJet1Pt->GetZaxis()->SetTitle("counts");
601 fOutput->Add(fHistDistancevsJet1Pt);
602
603 fHistDistancevsJet2Pt = new TH2F("fHistDistancevsJet2Pt", "fHistDistancevsJet2Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt);
604 fHistDistancevsJet2Pt->GetXaxis()->SetTitle("Distance");
605 fHistDistancevsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");
606 fHistDistancevsJet2Pt->GetZaxis()->SetTitle("counts");
607 fOutput->Add(fHistDistancevsJet2Pt);
cd6431de 608
6a20534a 609 fHistDistancevsCommonEnergy1 = new TH2F("fHistDistancevsCommonEnergy1", "fHistDistancevsCommonEnergy1", fNbins/2, 0, 1.2, fNbins/2, 0, 1.2);
610 fHistDistancevsCommonEnergy1->GetXaxis()->SetTitle("Distance");
611 fHistDistancevsCommonEnergy1->GetYaxis()->SetTitle("Common energy 1 (%)");
612 fHistDistancevsCommonEnergy1->GetZaxis()->SetTitle("counts");
613 fOutput->Add(fHistDistancevsCommonEnergy1);
614
615 fHistDistancevsCommonEnergy2 = new TH2F("fHistDistancevsCommonEnergy2", "fHistDistancevsCommonEnergy2", fNbins/2, 0, 1.2, fNbins/2, 0, 1.2);
616 fHistDistancevsCommonEnergy2->GetXaxis()->SetTitle("Distance");
617 fHistDistancevsCommonEnergy2->GetYaxis()->SetTitle("Common energy 2 (%)");
618 fHistDistancevsCommonEnergy2->GetZaxis()->SetTitle("counts");
619 fOutput->Add(fHistDistancevsCommonEnergy2);
cfc2ac24 620
05077f28 621 fHistCommonEnergy1vsCommonEnergy2 = new TH2F("fHistCommonEnergy1vsCommonEnergy2", "fHistCommonEnergy1vsCommonEnergy2", fNbins/2, 0, 1.2, fNbins/2, 0, 1.2);
622 fHistCommonEnergy1vsCommonEnergy2->GetXaxis()->SetTitle("Common energy 1 (%)");
623 fHistCommonEnergy1vsCommonEnergy2->GetYaxis()->SetTitle("Common energy 2 (%)");
624 fHistCommonEnergy1vsCommonEnergy2->GetZaxis()->SetTitle("counts");
625 fOutput->Add(fHistCommonEnergy1vsCommonEnergy2);
626
7030f36f 627 fHistJet2PtOverJet1PtvsJet2Pt = new TH2F("fHistJet2PtOverJet1PtvsJet2Pt", "fHistJet2PtOverJet1PtvsJet2Pt", fNbins, fMinBinPt, fMaxBinPt, 300, 0, 1.5);
c560b734 628 fHistJet2PtOverJet1PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");
629 fHistJet2PtOverJet1PtvsJet2Pt->GetYaxis()->SetTitle("p_{T,2} / p_{T,1}");
630 fHistJet2PtOverJet1PtvsJet2Pt->GetZaxis()->SetTitle("counts");
631 fOutput->Add(fHistJet2PtOverJet1PtvsJet2Pt);
632
7030f36f 633 fHistJet1PtOverJet2PtvsJet1Pt = new TH2F("fHistJet1PtOverJet2PtvsJet1Pt", "fHistJet1PtOverJet2PtvsJet1Pt", fNbins, fMinBinPt, fMaxBinPt, 300, 0, 1.5);
634 fHistJet1PtOverJet2PtvsJet1Pt->GetXaxis()->SetTitle("p_{T,1}");
635 fHistJet1PtOverJet2PtvsJet1Pt->GetYaxis()->SetTitle("p_{T,1} / p_{T,2}");
636 fHistJet1PtOverJet2PtvsJet1Pt->GetZaxis()->SetTitle("counts");
637 fOutput->Add(fHistJet1PtOverJet2PtvsJet1Pt);
638
ca5c29fa 639 fHistDeltaEtaPhi = new TH2F("fHistDeltaEtaPhi", "fHistDeltaEtaPhi", 200, -1, 1, 250, -1.6, 4.8);
05077f28 640 fHistDeltaEtaPhi->GetXaxis()->SetTitle("#delta#eta");
641 fHistDeltaEtaPhi->GetYaxis()->SetTitle("#delta#phi");
ca5c29fa 642 fHistDeltaEtaPhi->GetZaxis()->SetTitle("counts");
643 fOutput->Add(fHistDeltaEtaPhi);
cfc2ac24 644
7030f36f 645 fHistDeltaPtvsJet1Pt = new TH2F("fHistDeltaPtvsJet1Pt", "fHistDeltaPtvsJet1Pt", fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
646 fHistDeltaPtvsJet1Pt->GetXaxis()->SetTitle("p_{T,1}");
647 fHistDeltaPtvsJet1Pt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
648 fHistDeltaPtvsJet1Pt->GetZaxis()->SetTitle("counts");
649 fOutput->Add(fHistDeltaPtvsJet1Pt);
650
c560b734 651 fHistDeltaPtvsJet2Pt = new TH2F("fHistDeltaPtvsJet2Pt", "fHistDeltaPtvsJet2Pt", fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
cfc2ac24 652 fHistDeltaPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");
653 fHistDeltaPtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
654 fHistDeltaPtvsJet2Pt->GetZaxis()->SetTitle("counts");
655 fOutput->Add(fHistDeltaPtvsJet2Pt);
656
05077f28 657 fHistDeltaPtvsDistance = new TH2F("fHistDeltaPtvsDistance", "fHistDeltaPtvsDistance", fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
658 fHistDeltaPtvsDistance->GetXaxis()->SetTitle("Distance");
659 fHistDeltaPtvsDistance->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
660 fHistDeltaPtvsDistance->GetZaxis()->SetTitle("counts");
661 fOutput->Add(fHistDeltaPtvsDistance);
662
663 fHistDeltaPtvsCommonEnergy1 = new TH2F("fHistDeltaPtvsCommonEnergy1", "fHistDeltaPtvsCommonEnergy1", fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
664 fHistDeltaPtvsCommonEnergy1->GetXaxis()->SetTitle("Common energy 1 (%)");
665 fHistDeltaPtvsCommonEnergy1->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
666 fHistDeltaPtvsCommonEnergy1->GetZaxis()->SetTitle("counts");
667 fOutput->Add(fHistDeltaPtvsCommonEnergy1);
668
669 fHistDeltaPtvsCommonEnergy2 = new TH2F("fHistDeltaPtvsCommonEnergy2", "fHistDeltaPtvsCommonEnergy2", fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
670 fHistDeltaPtvsCommonEnergy2->GetXaxis()->SetTitle("Common energy 2 (%)");
671 fHistDeltaPtvsCommonEnergy2->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
672 fHistDeltaPtvsCommonEnergy2->GetZaxis()->SetTitle("counts");
673 fOutput->Add(fHistDeltaPtvsCommonEnergy2);
674
675 fHistDeltaPtvsArea1 = new TH2F("fHistDeltaPtvsArea1", "fHistDeltaPtvsArea1", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
676 fHistDeltaPtvsArea1->GetXaxis()->SetTitle("A_{jet,1}");
677 fHistDeltaPtvsArea1->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
678 fHistDeltaPtvsArea1->GetZaxis()->SetTitle("counts");
679 fOutput->Add(fHistDeltaPtvsArea1);
680
681 fHistDeltaPtvsArea2 = new TH2F("fHistDeltaPtvsArea2", "fHistDeltaPtvsArea2", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
682 fHistDeltaPtvsArea2->GetXaxis()->SetTitle("A_{jet,2}");
683 fHistDeltaPtvsArea2->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
684 fHistDeltaPtvsArea2->GetZaxis()->SetTitle("counts");
685 fOutput->Add(fHistDeltaPtvsArea2);
686
687 fHistDeltaPtvsDeltaArea = new TH2F("fHistDeltaPtvsDeltaArea", "fHistDeltaPtvsDeltaArea",
688 80, -fJetRadius * fJetRadius * TMath::Pi() * 3, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
689 fHistDeltaPtvsDeltaArea->GetXaxis()->SetTitle("#deltaA_{jet}");
690 fHistDeltaPtvsDeltaArea->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
691 fHistDeltaPtvsDeltaArea->GetZaxis()->SetTitle("counts");
692 fOutput->Add(fHistDeltaPtvsDeltaArea);
693
694 fHistJet1PtvsJet2Pt = new TH2F("fHistJet1PtvsJet2Pt", "fHistJet1PtvsJet2Pt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
695 fHistJet1PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,1}");
696 fHistJet1PtvsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");
697 fHistJet1PtvsJet2Pt->GetZaxis()->SetTitle("counts");
698 fOutput->Add(fHistJet1PtvsJet2Pt);
cd6431de 699
700 if (!fRhoName.IsNull() || !fRho2Name.IsNull()) {
7030f36f 701 fHistDeltaCorrPtvsJet1Pt = new TH2F("fHistDeltaCorrPtvsJet1Pt", "fHistDeltaCorrPtvsJet1Pt", fNbins/2, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
702 fHistDeltaCorrPtvsJet1Pt->GetXaxis()->SetTitle("p_{T,1}");
05077f28 703 fHistDeltaCorrPtvsJet1Pt->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
7030f36f 704 fHistDeltaCorrPtvsJet1Pt->GetZaxis()->SetTitle("counts");
705 fOutput->Add(fHistDeltaCorrPtvsJet1Pt);
706
cfc2ac24 707 fHistDeltaCorrPtvsJet2Pt = new TH2F("fHistDeltaCorrPtvsJet2Pt", "fHistDeltaCorrPtvsJet2Pt", fNbins/2, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
708 fHistDeltaCorrPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");
05077f28 709 fHistDeltaCorrPtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
cfc2ac24 710 fHistDeltaCorrPtvsJet2Pt->GetZaxis()->SetTitle("counts");
711 fOutput->Add(fHistDeltaCorrPtvsJet2Pt);
712
05077f28 713 fHistDeltaCorrPtvsDistance = new TH2F("fHistDeltaCorrPtvsDistance", "fHistDeltaCorrPtvsDistance", fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
714 fHistDeltaCorrPtvsDistance->GetXaxis()->SetTitle("Distance");
715 fHistDeltaCorrPtvsDistance->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
716 fHistDeltaCorrPtvsDistance->GetZaxis()->SetTitle("counts");
717 fOutput->Add(fHistDeltaCorrPtvsDistance);
718
719 fHistDeltaCorrPtvsCommonEnergy1 = new TH2F("fHistDeltaCorrPtvsCommonEnergy1", "fHistDeltaCorrPtvsCommonEnergy1", fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
720 fHistDeltaCorrPtvsCommonEnergy1->GetXaxis()->SetTitle("Common energy 1 (%)");
721 fHistDeltaCorrPtvsCommonEnergy1->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
722 fHistDeltaCorrPtvsCommonEnergy1->GetZaxis()->SetTitle("counts");
723 fOutput->Add(fHistDeltaCorrPtvsCommonEnergy1);
724
725 fHistDeltaCorrPtvsCommonEnergy2 = new TH2F("fHistDeltaCorrPtvsCommonEnergy2", "fHistDeltaCorrPtvsCommonEnergy2", fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
726 fHistDeltaCorrPtvsCommonEnergy2->GetXaxis()->SetTitle("Common energy 2 (%)");
727 fHistDeltaCorrPtvsCommonEnergy2->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
728 fHistDeltaCorrPtvsCommonEnergy2->GetZaxis()->SetTitle("counts");
729 fOutput->Add(fHistDeltaCorrPtvsCommonEnergy2);
730
731 fHistDeltaCorrPtvsArea1 = new TH2F("fHistDeltaCorrPtvsArea1", "fHistDeltaCorrPtvsArea1", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
732 fHistDeltaCorrPtvsArea1->GetXaxis()->SetTitle("A_{jet,1}");
733 fHistDeltaCorrPtvsArea1->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
734 fHistDeltaCorrPtvsArea1->GetZaxis()->SetTitle("counts");
735 fOutput->Add(fHistDeltaCorrPtvsArea1);
736
737 fHistDeltaCorrPtvsArea2 = new TH2F("fHistDeltaCorrPtvsArea2", "fHistDeltaCorrPtvsArea2", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
738 fHistDeltaCorrPtvsArea2->GetXaxis()->SetTitle("A_{jet,2}");
739 fHistDeltaCorrPtvsArea2->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
740 fHistDeltaCorrPtvsArea2->GetZaxis()->SetTitle("counts");
741 fOutput->Add(fHistDeltaCorrPtvsArea2);
742
743 fHistDeltaCorrPtvsDeltaArea = new TH2F("fHistDeltaCorrPtvsDeltaArea", "fHistDeltaCorrPtvsDeltaArea",
744 80, -fJetRadius * fJetRadius * TMath::Pi() * 3, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
745 fHistDeltaCorrPtvsDeltaArea->GetXaxis()->SetTitle("#deltaA_{jet}");
746 fHistDeltaCorrPtvsDeltaArea->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
747 fHistDeltaCorrPtvsDeltaArea->GetZaxis()->SetTitle("counts");
748 fOutput->Add(fHistDeltaCorrPtvsDeltaArea);
749
cd6431de 750 if (fRhoName.IsNull())
751 fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt", fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
752 else if (fRho2Name.IsNull())
753 fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt", 2*fNbins, -fMaxBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
754 else
755 fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt", 2*fNbins, -fMaxBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
756 fHistJet1CorrPtvsJet2CorrPt->GetXaxis()->SetTitle("p_{T,1}^{corr}");
757 fHistJet1CorrPtvsJet2CorrPt->GetYaxis()->SetTitle("p_{T,2}^{corr}");
758 fHistJet1CorrPtvsJet2CorrPt->GetZaxis()->SetTitle("counts");
759 fOutput->Add(fHistJet1CorrPtvsJet2CorrPt);
760 }
761
4358e58a 762 if (fIsEmbedded) {
05077f28 763 fHistDeltaMCPtvsJet1Pt = new TH2F("fHistDeltaMCPtvsJet1Pt", "fHistDeltaMCPtvsJet1Pt", fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
764 fHistDeltaMCPtvsJet1Pt->GetXaxis()->SetTitle("p_{T,1}");
765 fHistDeltaMCPtvsJet1Pt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
766 fHistDeltaMCPtvsJet1Pt->GetZaxis()->SetTitle("counts");
767 fOutput->Add(fHistDeltaMCPtvsJet1Pt);
768
769 fHistDeltaMCPtvsJet2Pt = new TH2F("fHistDeltaMCPtvsJet2Pt", "fHistDeltaMCPtvsJet2Pt", fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
770 fHistDeltaMCPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");
771 fHistDeltaMCPtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
772 fHistDeltaMCPtvsJet2Pt->GetZaxis()->SetTitle("counts");
773 fOutput->Add(fHistDeltaMCPtvsJet2Pt);
774
775 fHistDeltaMCPtvsDistance = new TH2F("fHistDeltaMCPtvsDistance", "fHistDeltaMCPtvsDistance", fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
776 fHistDeltaMCPtvsDistance->GetXaxis()->SetTitle("Distance");
777 fHistDeltaMCPtvsDistance->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
778 fHistDeltaMCPtvsDistance->GetZaxis()->SetTitle("counts");
779 fOutput->Add(fHistDeltaMCPtvsDistance);
780
781 fHistDeltaMCPtvsCommonEnergy1 = new TH2F("fHistDeltaMCPtvsCommonEnergy1", "fHistDeltaMCPtvsCommonEnergy1", fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
782 fHistDeltaMCPtvsCommonEnergy1->GetXaxis()->SetTitle("Common energy 1 (%)");
783 fHistDeltaMCPtvsCommonEnergy1->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
784 fHistDeltaMCPtvsCommonEnergy1->GetZaxis()->SetTitle("counts");
785 fOutput->Add(fHistDeltaMCPtvsCommonEnergy1);
786
787 fHistDeltaMCPtvsCommonEnergy2 = new TH2F("fHistDeltaMCPtvsCommonEnergy2", "fHistDeltaMCPtvsCommonEnergy2", fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
788 fHistDeltaMCPtvsCommonEnergy2->GetXaxis()->SetTitle("Common energy 2 (%)");
789 fHistDeltaMCPtvsCommonEnergy2->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
790 fHistDeltaMCPtvsCommonEnergy2->GetZaxis()->SetTitle("counts");
791 fOutput->Add(fHistDeltaMCPtvsCommonEnergy2);
792
793 fHistDeltaMCPtvsArea1 = new TH2F("fHistDeltaMCPtvsArea1", "fHistDeltaMCPtvsArea1", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
794 fHistDeltaMCPtvsArea1->GetXaxis()->SetTitle("A_{jet,1}");
795 fHistDeltaMCPtvsArea1->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
796 fHistDeltaMCPtvsArea1->GetZaxis()->SetTitle("counts");
797 fOutput->Add(fHistDeltaMCPtvsArea1);
798
799 fHistDeltaMCPtvsArea2 = new TH2F("fHistDeltaMCPtvsArea2", "fHistDeltaMCPtvsArea2", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
800 fHistDeltaMCPtvsArea2->GetXaxis()->SetTitle("A_{jet,2}");
801 fHistDeltaMCPtvsArea2->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
802 fHistDeltaMCPtvsArea2->GetZaxis()->SetTitle("counts");
803 fOutput->Add(fHistDeltaMCPtvsArea2);
804
805 fHistDeltaMCPtvsDeltaArea = new TH2F("fHistDeltaMCPtvsDeltaArea", "fHistDeltaMCPtvsDeltaArea",
806 80, -fJetRadius * fJetRadius * TMath::Pi() * 3, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
807 fHistDeltaMCPtvsDeltaArea->GetXaxis()->SetTitle("#deltaA_{jet}");
808 fHistDeltaMCPtvsDeltaArea->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
809 fHistDeltaMCPtvsDeltaArea->GetZaxis()->SetTitle("counts");
810 fOutput->Add(fHistDeltaMCPtvsDeltaArea);
811
4358e58a 812 fHistJet1MCPtvsJet2Pt = new TH2F("fHistJet1MCPtvsJet2Pt", "fHistJet1MCPtvsJet2Pt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
813 fHistJet1MCPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,1}^{MC}");
814 fHistJet1MCPtvsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");
815 fHistJet1MCPtvsJet2Pt->GetZaxis()->SetTitle("counts");
816 fOutput->Add(fHistJet1MCPtvsJet2Pt);
817 }
818
99c67c1b 819 PostData(1, fOutput); // Post data for ALL output slots > 0 here, to get at least an empty histogram
2949a2ec 820}
821
ceefbfbc 822//________________________________________________________________________
a487deae 823Bool_t AliJetResponseMaker::AcceptJet(AliEmcalJet *jet) const
ceefbfbc 824{
825 // Return true if jet is accepted.
826
827 if (jet->Pt() <= fJetPtCut)
828 return kFALSE;
829 if (jet->Area() <= fJetAreaCut)
830 return kFALSE;
65bb5510 831
832 return kTRUE;
ceefbfbc 833}
834
cd6431de 835//________________________________________________________________________
836Bool_t AliJetResponseMaker::AcceptBiasJet2(AliEmcalJet *jet) const
837{
838 // Accept jet with a bias.
839
840 if (fLeadingHadronType == 0) {
841 if (jet->MaxTrackPt() < fPtBiasJet2Track) return kFALSE;
842 }
843 else if (fLeadingHadronType == 1) {
844 if (jet->MaxClusterPt() < fPtBiasJet2Clus) return kFALSE;
845 }
846 else {
847 if (jet->MaxTrackPt() < fPtBiasJet2Track && jet->MaxClusterPt() < fPtBiasJet2Clus) return kFALSE;
848 }
849
850 return kTRUE;
851}
852
b3ccdfe2 853//________________________________________________________________________
854void AliJetResponseMaker::ExecOnce()
855{
856 // Execute once.
857
cd6431de 858 if (!fJets2Name.IsNull() && !fJets2) {
859 fJets2 = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJets2Name));
860 if (!fJets2) {
861 AliError(Form("%s: Could not retrieve jets2 %s!", GetName(), fJets2Name.Data()));
b3ccdfe2 862 return;
863 }
cd6431de 864 else if (!fJets2->GetClass()->GetBaseClass("AliEmcalJet")) {
865 AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fJets2Name.Data()));
866 fJets2 = 0;
b3ccdfe2 867 return;
868 }
869 }
870
cd6431de 871 if (!fTracks2Name.IsNull() && !fTracks2) {
872 fTracks2 = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracks2Name));
873 if (!fTracks2) {
874 AliError(Form("%s: Could not retrieve tracks2 %s!", GetName(), fTracks2Name.Data()));
b3ccdfe2 875 return;
876 }
877 else {
cd6431de 878 TClass *cl = fTracks2->GetClass();
b3ccdfe2 879 if (!cl->GetBaseClass("AliVParticle") && !cl->GetBaseClass("AliEmcalParticle")) {
cd6431de 880 AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fTracks2Name.Data()));
881 fTracks2 = 0;
b3ccdfe2 882 return;
883 }
884 }
cfc2ac24 885
886 if (fAreCollections2MC) {
5be3857d 887 fTracks2Map = dynamic_cast<AliNamedArrayI*>(InputEvent()->FindListObject(fTracks2Name + "_Map"));
7f76e479 888 // this is needed to map the MC labels with the indexes of the MC particle collection
889 // 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 890 if (!fTracks2Map) {
7f76e479 891 AliWarning(Form("%s: Could not retrieve map for tracks2 %s! Will assume MC labels consistent with indexes...", GetName(), fTracks2Name.Data()));
5be3857d 892 fTracks2Map = new AliNamedArrayI("tracksMap",9999);
7f76e479 893 for (Int_t i = 0; i < 9999; i++) {
5be3857d 894 fTracks2Map->AddAt(i,i);
7f76e479 895 }
cfc2ac24 896 }
897 }
b3ccdfe2 898 }
899
cd6431de 900 if (!fCalo2Name.IsNull() && !fCaloClusters2) {
901 fCaloClusters2 = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCalo2Name));
902 if (!fCaloClusters2) {
903 AliError(Form("%s: Could not retrieve clusters %s!", GetName(), fCalo2Name.Data()));
904 return;
905 } else {
906 TClass *cl = fCaloClusters2->GetClass();
907 if (!cl->GetBaseClass("AliVCluster") && !cl->GetBaseClass("AliEmcalParticle")) {
908 AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fCalo2Name.Data()));
909 fCaloClusters2 = 0;
910 return;
911 }
912 }
91bee8bc 913 }
cd6431de 914
915 if (!fRho2Name.IsNull() && !fRho2) {
916 fRho2 = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRho2Name));
917 if (!fRho2) {
918 AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRho2Name.Data()));
919 fInitialized = kFALSE;
920 return;
921 }
91bee8bc 922 }
cd6431de 923
924 if (fJet2MinEta == -999)
925 fJet2MinEta = fJetMinEta - fJetRadius;
926 if (fJet2MaxEta == -999)
927 fJet2MaxEta = fJetMaxEta + fJetRadius;
928 if (fJet2MinPhi == -999)
929 fJet2MinPhi = fJetMinPhi - fJetRadius;
930 if (fJet2MaxPhi == -999)
931 fJet2MaxPhi = fJetMaxPhi + fJetRadius;
932
933 AliAnalysisTaskEmcalJet::ExecOnce();
b3ccdfe2 934}
935
4643d2e8 936//________________________________________________________________________
937Bool_t AliJetResponseMaker::IsEventSelected()
938{
939 // Check if event is selected
940
941 if (fSelectPtHardBin != -999 && fSelectPtHardBin != fPtHardBin)
942 return kFALSE;
943
944 return AliAnalysisTaskEmcalJet::IsEventSelected();
945}
946
b3ccdfe2 947//________________________________________________________________________
948Bool_t AliJetResponseMaker::RetrieveEventObjects()
949{
950 // Retrieve event objects.
951
952 if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
953 return kFALSE;
cd6431de 954
955 if (fRho2)
956 fRho2Val = fRho2->GetVal();
b3ccdfe2 957
b3ccdfe2 958 const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
959 const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
787a3c4f 960
05077f28 961 if (MCEvent()) {
787a3c4f 962 fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
05077f28 963 if (!fPythiaHeader) {
97f1b773 964 // Check if AOD
965 AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
966
05077f28 967 if (aodMCH) {
968 for(UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++) {
97f1b773 969 fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
05077f28 970 if (fPythiaHeader) break;
97f1b773 971 }
05077f28 972 }
97f1b773 973 }
974 }
975
787a3c4f 976 if (fPythiaHeader) {
977 Double_t pthard = fPythiaHeader->GetPtHard();
978
979 for (fPtHardBin = 0; fPtHardBin < 11; fPtHardBin++) {
980 if (pthard >= ptHardLo[fPtHardBin] && pthard < ptHardHi[fPtHardBin])
981 break;
982 }
983
984 fNTrials = fPythiaHeader->Trials();
b3ccdfe2 985 }
b3ccdfe2 986
b3ccdfe2 987 return kTRUE;
988}
989
990//________________________________________________________________________
991Bool_t AliJetResponseMaker::Run()
992{
993 // Find the closest jets
cfc2ac24 994
f660c2d6 995 if (fMatching == kNoMatching)
aa4d701c 996 return kTRUE;
f660c2d6 997 else
998 return DoJetMatching();
cfc2ac24 999}
aa4d701c 1000
cfc2ac24 1001//________________________________________________________________________
f660c2d6 1002Bool_t AliJetResponseMaker::DoJetMatching()
cfc2ac24 1003{
1004 DoJetLoop(kFALSE);
cfc2ac24 1005
f660c2d6 1006 const Int_t nJets = fJets->GetEntriesFast();
cfc2ac24 1007
f660c2d6 1008 for (Int_t i = 0; i < nJets; i++) {
cfc2ac24 1009
1010 AliEmcalJet* jet1 = static_cast<AliEmcalJet*>(fJets->At(i));
1011
1012 if (!jet1) {
1013 AliError(Form("Could not receive jet %d", i));
1014 continue;
1015 }
1016
1017 if (!AcceptJet(jet1))
1018 continue;
1019
1020 if (jet1->Eta() < fJetMinEta || jet1->Eta() > fJetMaxEta || jet1->Phi() < fJetMinPhi || jet1->Phi() > fJetMaxPhi)
1021 continue;
1022
f660c2d6 1023 if (jet1->ClosestJet() && jet1->ClosestJet()->ClosestJet() == jet1 &&
7f76e479 1024 jet1->ClosestJetDistance() < fMatchingPar1 && jet1->ClosestJet()->ClosestJetDistance() < fMatchingPar2) { // Matched jet found
cfc2ac24 1025 jet1->SetMatchedToClosest(fMatching);
cfc2ac24 1026 jet1->ClosestJet()->SetMatchedToClosest(fMatching);
aa4d701c 1027 }
1028 }
1029
b3ccdfe2 1030 return kTRUE;
1031}
1032
2949a2ec 1033//________________________________________________________________________
cfc2ac24 1034void AliJetResponseMaker::DoJetLoop(Bool_t order)
2949a2ec 1035{
44476be7 1036 // Do the jet loop.
1ee1b5b8 1037
cfc2ac24 1038 TClonesArray *jets1 = 0;
1039 TClonesArray *jets2 = 0;
1040
1041 if (order) {
1042 jets1 = fJets2;
1043 jets2 = fJets;
1044 }
1045 else {
1046 jets1 = fJets;
1047 jets2 = fJets2;
1048 }
1049
44476be7 1050 Int_t nJets1 = jets1->GetEntriesFast();
1051 Int_t nJets2 = jets2->GetEntriesFast();
1ee1b5b8 1052
05077f28 1053 for (Int_t j = 0; j < nJets2; j++) {
1054
1055 AliEmcalJet* jet2 = static_cast<AliEmcalJet*>(jets2->At(j));
1056
1057 if (!jet2) {
1058 AliError(Form("Could not receive jet %d", j));
1059 continue;
1060 }
2103dc6a 1061
05077f28 1062 jet2->ResetMatching();
fde82e42 1063
1064
1065 if (!AcceptJet(jet2))
1066 continue;
1067
1068 if (jet2->Eta() < fJet2MinEta || jet2->Eta() > fJet2MaxEta || jet2->Phi() < fJet2MinPhi || jet2->Phi() > fJet2MaxPhi)
1069 continue;
05077f28 1070 }
1071
44476be7 1072 for (Int_t i = 0; i < nJets1; i++) {
1ee1b5b8 1073
44476be7 1074 AliEmcalJet* jet1 = static_cast<AliEmcalJet*>(jets1->At(i));
2949a2ec 1075
44476be7 1076 if (!jet1) {
1077 AliError(Form("Could not receive jet %d", i));
1078 continue;
2103dc6a 1079 }
1080
1081 jet1->ResetMatching();
44476be7 1082
ceefbfbc 1083 if (!AcceptJet(jet1))
44476be7 1084 continue;
1085
fde82e42 1086 if (jet1->MCPt() < 1)
1087 continue;
1088
cfc2ac24 1089 if (order) {
cd6431de 1090 if (jet1->Eta() < fJet2MinEta || jet1->Eta() > fJet2MaxEta || jet1->Phi() < fJet2MinPhi || jet1->Phi() > fJet2MaxPhi)
ceefbfbc 1091 continue;
1092 }
1093 else {
cd6431de 1094 if (jet1->Eta() < fJetMinEta || jet1->Eta() > fJetMaxEta || jet1->Phi() < fJetMinPhi || jet1->Phi() > fJetMaxPhi)
ceefbfbc 1095 continue;
1096 }
1097
44476be7 1098 for (Int_t j = 0; j < nJets2; j++) {
1099
1100 AliEmcalJet* jet2 = static_cast<AliEmcalJet*>(jets2->At(j));
1101
1102 if (!jet2) {
1103 AliError(Form("Could not receive jet %d", j));
1104 continue;
2103dc6a 1105 }
44476be7 1106
ceefbfbc 1107 if (!AcceptJet(jet2))
44476be7 1108 continue;
ceefbfbc 1109
cfc2ac24 1110 if (order) {
a487deae 1111 if (jet2->Eta() < fJetMinEta || jet2->Eta() > fJetMaxEta || jet2->Phi() < fJetMinPhi || jet2->Phi() > fJetMaxPhi)
ceefbfbc 1112 continue;
1113 }
1114 else {
fde82e42 1115 if (jet2->Eta() < fJet2MinEta || jet2->Eta() > fJet2MaxEta || jet2->Phi() < fJet2MinPhi || jet2->Phi() > fJet2MaxPhi)
ceefbfbc 1116 continue;
1117 }
cd6431de 1118
7f76e479 1119 SetMatchingLevel(jet1, jet2, fMatching);
2103dc6a 1120 } // jet2 loop
1121
2103dc6a 1122 } // jet1 loop
1ee1b5b8 1123}
1124
cd6431de 1125//________________________________________________________________________
7f76e479 1126void AliJetResponseMaker::GetGeometricalMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d) const
cd6431de 1127{
7f76e479 1128 Double_t deta = jet2->Eta() - jet1->Eta();
1129 Double_t dphi = jet2->Phi() - jet1->Phi();
1130 d = TMath::Sqrt(deta * deta + dphi * dphi);
1131}
cd6431de 1132
7f76e479 1133//________________________________________________________________________
1134void AliJetResponseMaker::GetMCLabelMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d1, Double_t &d2) const
1135{
1136 // d1 and d2 represent the matching level: 0 = maximum level of matching, 1 = the two jets are completely unrelated
1137 d1 = jet1->Pt();
1138 d2 = jet2->Pt();
1139 Double_t totalPt1 = d1; // the total pt of the reconstructed jet will be cleaned from the background
1140
1141 for (Int_t iTrack2 = 0; iTrack2 < jet2->GetNumberOfTracks(); iTrack2++) {
1142 Bool_t track2Found = kFALSE;
1143 Int_t index2 = jet2->TrackAt(iTrack2);
1144 for (Int_t iTrack = 0; iTrack < jet1->GetNumberOfTracks(); iTrack++) {
1145 AliVParticle *track = jet1->TrackAt(iTrack,fTracks);
1146 if (!track) {
1147 AliWarning(Form("Could not find track %d!", iTrack));
1148 continue;
1149 }
787a3c4f 1150 Int_t MClabel = TMath::Abs(track->GetLabel());
7f76e479 1151 Int_t index = -1;
1152
787a3c4f 1153 if (MClabel == 0) {// this is not a MC particle; remove it completely
7f76e479 1154 AliDebug(3,Form("Track %d (pT = %f) is not a MC particle (MClabel = %d)!",iTrack,track->Pt(),MClabel));
1155 totalPt1 -= track->Pt();
1156 d1 -= track->Pt();
1157 continue;
1158 }
5be3857d 1159 else if (MClabel < fTracks2Map->GetSize()) {
1160 index = fTracks2Map->At(MClabel);
7f76e479 1161 }
1162
1163 if (index < 0) {
1164 AliDebug(2,Form("Track %d (pT = %f) does not have an associated MC particle (MClabel = %d)!",iTrack,track->Pt(),MClabel));
1165 continue;
1166 }
1167
1168 if (index2 == index) { // found common particle
1169 track2Found = kTRUE;
1170 d1 -= track->Pt();
1171 AliVParticle *MCpart = static_cast<AliVParticle*>(fTracks2->At(index2));
1172 AliDebug(3,Form("Track %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!",
1173 iTrack,track->Pt(),track->Eta(),track->Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
1174 d2 -= MCpart->Pt();
1175 break;
1176 }
cd6431de 1177 }
7f76e479 1178 for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
1179 AliVCluster *clus = jet1->ClusterAt(iClus,fCaloClusters);
1180 if (!clus) {
1181 AliWarning(Form("Could not find cluster %d!", iClus));
1182 continue;
1183 }
1184 TLorentzVector part;
1185 clus->GetMomentum(part, const_cast<Double_t*>(fVertex));
1186
1187 if (fCaloCells) { // if the cell colection is available, look for cells with a matched MC particle
1188 for (Int_t iCell = 0; iCell < clus->GetNCells(); iCell++) {
1189 Int_t cellId = clus->GetCellAbsId(iCell);
1190 Double_t cellFrac = clus->GetCellAmplitudeFraction(iCell);
1191
787a3c4f 1192 Int_t MClabel = TMath::Abs(fCaloCells->GetCellMCLabel(cellId));
7f76e479 1193 Int_t index = -1;
1194
787a3c4f 1195 if (MClabel == 0) {// this is not a MC particle; remove it completely
7f76e479 1196 AliDebug(3,Form("Cell %d (frac = %f) is not a MC particle (MClabel = %d)!",iCell,cellFrac,MClabel));
1197 totalPt1 -= part.Pt() * cellFrac;
1198 d1 -= part.Pt() * cellFrac;
1199 continue;
1200 }
5be3857d 1201 else if (MClabel < fTracks2Map->GetSize()) {
1202 index = fTracks2Map->At(MClabel);
7f76e479 1203 }
1204
1205 if (index < 0) {
1206 AliDebug(3,Form("Cell %d (frac = %f) does not have an associated MC particle (MClabel = %d)!",iCell,cellFrac,MClabel));
1207 continue;
1208 }
1209 if (index2 == index) { // found common particle
1210 d1 -= part.Pt() * cellFrac;
1211
1212 if (!track2Found) {// only if it is not already found among charged tracks (charged particles are most likely already found)
1213 AliVParticle *MCpart = static_cast<AliVParticle*>(fTracks2->At(index2));
1214 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)!",
1215 iCell,iClus,part.Pt(),part.Eta(),part.Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
1216 d2 -= MCpart->Pt() * cellFrac;
1217 }
1218 break;
1219 }
cfc2ac24 1220 }
1221 }
7f76e479 1222 else { //otherwise look for the first contributor to the cluster, and if matched to a MC label remove it
787a3c4f 1223 Int_t MClabel = TMath::Abs(clus->GetLabel());
7f76e479 1224 Int_t index = -1;
1225
787a3c4f 1226 if (MClabel == 0) {// this is not a MC particle; remove it completely
7f76e479 1227 AliDebug(3,Form("Cluster %d (pT = %f) is not a MC particle (MClabel = %d)!",iClus,part.Pt(),MClabel));
1228 totalPt1 -= part.Pt();
1229 d1 -= part.Pt();
cfc2ac24 1230 continue;
1231 }
5be3857d 1232 else if (MClabel < fTracks2Map->GetSize()) {
1233 index = fTracks2Map->At(MClabel);
f660c2d6 1234 }
7f76e479 1235
cfc2ac24 1236 if (index < 0) {
7f76e479 1237 AliDebug(3,Form("Cluster %d (pT = %f) does not have an associated MC particle (MClabel = %d)!",iClus,part.Pt(),MClabel));
cfc2ac24 1238 continue;
1239 }
7f76e479 1240 if (index2 == index) { // found common particle
1241 d1 -= part.Pt();
1242
1243 if (!track2Found) {// only if it is not already found among charged tracks (charged particles are most likely already found)
f660c2d6 1244 AliVParticle *MCpart = static_cast<AliVParticle*>(fTracks2->At(index2));
7f76e479 1245 AliDebug(3,Form("Cluster %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!",
1246 iClus,part.Pt(),part.Eta(),part.Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
1247
f660c2d6 1248 d2 -= MCpart->Pt();
cfc2ac24 1249 }
7f76e479 1250 break;
cfc2ac24 1251 }
1252 }
7f76e479 1253 }
1254 }
05077f28 1255
1256 if (d1 < 0)
7f76e479 1257 d1 = 0;
05077f28 1258
1259 if (d2 < 0)
1260 d2 = 0;
1261
1262 if (totalPt1 < 1)
1263 d1 = -1;
7f76e479 1264 else
1265 d1 /= totalPt1;
f660c2d6 1266
05077f28 1267 if (jet2->Pt() < 1)
1268 d2 = -1;
7f76e479 1269 else
05077f28 1270 d2 /= jet2->Pt();
7f76e479 1271}
1272
1273//________________________________________________________________________
1274void AliJetResponseMaker::GetSameCollectionsMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d1, Double_t &d2) const
1275{
1276 // d1 and d2 represent the matching level: 0 = maximum level of matching, 1 = the two jets are completely unrelated
1277 d1 = jet1->Pt();
1278 d2 = jet2->Pt();
1279
1280 if (fTracks && fTracks2) {
1281
1282 for (Int_t iTrack2 = 0; iTrack2 < jet2->GetNumberOfTracks(); iTrack2++) {
1283 Int_t index2 = jet2->TrackAt(iTrack2);
1284 for (Int_t iTrack = 0; iTrack < jet1->GetNumberOfTracks(); iTrack++) {
1285 Int_t index = jet1->TrackAt(iTrack);
1286 if (index2 == index) { // found common particle
1287 AliVParticle *part = static_cast<AliVParticle*>(fTracks->At(index));
1288 if (!part) {
1289 AliWarning(Form("Could not find track %d!", index));
f660c2d6 1290 continue;
1291 }
7f76e479 1292 AliVParticle *part2 = static_cast<AliVParticle*>(fTracks2->At(index2));
1293 if (!part2) {
1294 AliWarning(Form("Could not find track %d!", index2));
cfc2ac24 1295 continue;
1296 }
7f76e479 1297
1298 d1 -= part->Pt();
1299 d2 -= part2->Pt();
1300 break;
cfc2ac24 1301 }
1302 }
cfc2ac24 1303 }
7f76e479 1304
1305 }
1306
1307 if (fCaloClusters && fCaloClusters2) {
1308
1309 for (Int_t iClus2 = 0; iClus2 < jet2->GetNumberOfClusters(); iClus2++) {
1310 Int_t index2 = jet2->ClusterAt(iClus2);
1311 for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
1312 Int_t index = jet1->ClusterAt(iClus);
cac0e10b 1313 if (index2 == index) { // found common particle
1314 AliVCluster *clus = static_cast<AliVCluster*>(fCaloClusters->At(index));
1315 if (!clus) {
1316 AliWarning(Form("Could not find cluster %d!", index));
1317 continue;
1318 }
1319 AliVCluster *clus2 = static_cast<AliVCluster*>(fCaloClusters2->At(index2));
1320 if (!clus2) {
1321 AliWarning(Form("Could not find cluster %d!", index2));
1322 continue;
1323 }
1324 TLorentzVector part, part2;
1325 clus->GetMomentum(part, const_cast<Double_t*>(fVertex));
1326 clus2->GetMomentum(part2, const_cast<Double_t*>(fVertex));
1327
1328 d1 -= part.Pt();
1329 d2 -= part2.Pt();
1330 break;
7f76e479 1331 }
7f76e479 1332 }
1333 }
1334
1335 }
1336
05077f28 1337 if (d1 < 0)
1338 d1 = 0;
1339
1340 if (d2 < 0)
1341 d2 = 0;
1342
1343 if (jet1->Pt() > 0)
7f76e479 1344 d1 /= jet1->Pt();
1345 else
05077f28 1346 d1 = -1;
7f76e479 1347
05077f28 1348 if (jet2->Pt() > 0)
7f76e479 1349 d2 /= jet2->Pt();
1350 else
05077f28 1351 d2 = -1;
7f76e479 1352}
1353
1354//________________________________________________________________________
1355void AliJetResponseMaker::SetMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, MatchingType matching)
1356{
1357 Double_t d1 = -1;
1358 Double_t d2 = -1;
1359
1360 switch (matching) {
1361 case kGeometrical:
1362 GetGeometricalMatchingLevel(jet1,jet2,d1);
1363 d2 = d1;
1364 break;
1365 case kMCLabel: // jet1 = detector level and jet2 = particle level!
1366 GetMCLabelMatchingLevel(jet1,jet2,d1,d2);
1367 break;
1368 case kSameCollections:
1369 GetSameCollectionsMatchingLevel(jet1,jet2,d1,d2);
cd6431de 1370 break;
1371 default:
1372 ;
1373 }
f660c2d6 1374
05077f28 1375 if (d1 >= 0) {
f660c2d6 1376
1377 if (d1 < jet1->ClosestJetDistance()) {
1378 jet1->SetSecondClosestJet(jet1->ClosestJet(), jet1->ClosestJetDistance());
1379 jet1->SetClosestJet(jet2, d1);
1380 }
1381 else if (d1 < jet1->SecondClosestJetDistance()) {
1382 jet1->SetSecondClosestJet(jet2, d1);
1383 }
1384 }
1385
05077f28 1386 if (d2 >= 0) {
f660c2d6 1387
1388 if (d2 < jet2->ClosestJetDistance()) {
1389 jet2->SetSecondClosestJet(jet2->ClosestJet(), jet2->ClosestJetDistance());
1390 jet2->SetClosestJet(jet1, d2);
1391 }
1392 else if (d2 < jet2->SecondClosestJetDistance()) {
1393 jet2->SetSecondClosestJet(jet1, d2);
1394 }
1395 }
cd6431de 1396}
1397
1ee1b5b8 1398//________________________________________________________________________
1399Bool_t AliJetResponseMaker::FillHistograms()
1400{
1401 // Fill histograms.
1402
ca5c29fa 1403 static Int_t indexes[9999] = {-1};
1404
05077f28 1405 if (fHistEventsAfterSel)
1406 fHistEventsAfterSel->SetBinContent(fPtHardBin + 1, fHistEventsAfterSel->GetBinContent(fPtHardBin + 1) + 1);
1407 if (fHistTrialsAfterSel)
1408 fHistTrialsAfterSel->SetBinContent(fPtHardBin + 1, fHistTrialsAfterSel->GetBinContent(fPtHardBin + 1) + fNTrials);
ca5c29fa 1409
1410 GetSortedArray(indexes, fJets2, fRho2Val);
4643d2e8 1411
cd6431de 1412 const Int_t nJets2 = fJets2->GetEntriesFast();
2949a2ec 1413
ca5c29fa 1414 Int_t naccJets2 = 0;
1415 Int_t naccJets2Acceptance = 0;
1416
fde82e42 1417 Double_t totalPt2 = 0;
1418
cd6431de 1419 for (Int_t i = 0; i < nJets2; i++) {
2949a2ec 1420
ca5c29fa 1421 AliEmcalJet* jet2 = static_cast<AliEmcalJet*>(fJets2->At(indexes[i]));
2949a2ec 1422
cd6431de 1423 if (!jet2) {
1424 AliError(Form("Could not receive jet2 %d", i));
2949a2ec 1425 continue;
cfc2ac24 1426 }
1427
cd6431de 1428 if (!AcceptJet(jet2))
1429 continue;
1430
cfc2ac24 1431 if (AcceptBiasJet(jet2) &&
1432 (jet2->Eta() > fJetMinEta && jet2->Eta() < fJetMaxEta && jet2->Phi() > fJetMinPhi && jet2->Phi() < fJetMaxPhi)) {
1433
1434 fHistJets2PtAreaAcceptance->Fill(jet2->Area(), jet2->Pt());
1435 fHistJets2PhiEtaAcceptance->Fill(jet2->Eta(), jet2->Phi());
fde82e42 1436
1437 totalPt2 += jet2->Pt();
cfc2ac24 1438
ca5c29fa 1439 if (naccJets2Acceptance < fNLeadingJets)
1440 fHistLeadingJets2PtAreaAcceptance->Fill(jet2->Area(), jet2->Pt());
1441
1442 if (!fRho2Name.IsNull()) {
cfc2ac24 1443 fHistJets2CorrPtAreaAcceptance->Fill(jet2->Area(), jet2->Pt() - fRho2Val * jet2->Area());
ca5c29fa 1444 if (naccJets2Acceptance < fNLeadingJets)
1445 fHistLeadingJets2CorrPtAreaAcceptance->Fill(jet2->Area(), jet2->Pt() - fRho2Val * jet2->Area());
1446 }
63fac07f 1447
1448 if (fTracks2) {
1449 for (Int_t it = 0; it < jet2->GetNumberOfTracks(); it++) {
1450 AliVParticle *track2 = jet2->TrackAt(it, fTracks2);
1451 if (track2)
1452 fHistJets2ZvsPt->Fill(track2->Pt() / jet2->Pt(), jet2->Pt());
1453 }
1454 }
1455
1456 if (fCaloClusters2) {
1457 for (Int_t ic = 0; ic < jet2->GetNumberOfClusters(); ic++) {
1458 AliVCluster *cluster2 = jet2->ClusterAt(ic, fCaloClusters2);
1459
1460 if (cluster2) {
1461 TLorentzVector nPart2;
1462 cluster2->GetMomentum(nPart2, fVertex);
1463 fHistJets2ZvsPt->Fill(nPart2.Et() / jet2->Pt(), jet2->Pt());
1464 }
1465 }
1466 }
1467
1468 fHistJets2NEFvsPt->Fill(jet2->NEF(), jet2->Pt());
1469 fHistJets2CEFvsCEFPt->Fill(1-jet2->NEF(), (1-jet2->NEF())*jet2->Pt());
ca5c29fa 1470
1471 naccJets2Acceptance++;
cfc2ac24 1472 }
1473
cd6431de 1474 if (!AcceptBiasJet2(jet2))
ceefbfbc 1475 continue;
1476
cd6431de 1477 if (jet2->Eta() < fJet2MinEta || jet2->Eta() > fJet2MaxEta || jet2->Phi() < fJet2MinPhi || jet2->Phi() > fJet2MaxPhi)
2bddb6ae 1478 continue;
ca5c29fa 1479
cfc2ac24 1480 fHistJets2PtArea->Fill(jet2->Area(), jet2->Pt());
1481 fHistJets2PhiEta->Fill(jet2->Eta(), jet2->Phi());
ca5c29fa 1482
1483 if (naccJets2 < fNLeadingJets)
1484 fHistLeadingJets2PtArea->Fill(jet2->Area(), jet2->Pt());
cfc2ac24 1485
ca5c29fa 1486 if (!fRho2Name.IsNull()) {
cfc2ac24 1487 fHistJets2CorrPtArea->Fill(jet2->Area(), jet2->Pt() - fRho2Val * jet2->Area());
ca5c29fa 1488 if (naccJets2 < fNLeadingJets)
1489 fHistLeadingJets2CorrPtArea->Fill(jet2->Area(), jet2->Pt() - fRho2Val * jet2->Area());
1490 }
1491
1492 naccJets2++;
2949a2ec 1493
cd6431de 1494 if (jet2->MatchedJet()) {
aa4d701c 1495
cd6431de 1496 if (!AcceptBiasJet(jet2->MatchedJet()) ||
ca5c29fa 1497 jet2->MatchedJet()->MaxTrackPt() > fMaxTrackPt || jet2->MatchedJet()->MaxClusterPt() > fMaxClusterPt) {
cd6431de 1498 fHistMissedJets2PtArea->Fill(jet2->Area(), jet2->Pt());
2bddb6ae 1499 }
1500 else {
ca5c29fa 1501 if (jet2->MatchedJet()->Pt() > fMaxBinPt)
1502 fHistMissedJets2PtArea->Fill(jet2->Area(), jet2->Pt());
1503
05077f28 1504 Double_t d=-1, ce1=-1, ce2=-1;
7f76e479 1505 if (jet2->GetMatchingType() == kGeometrical) {
05077f28 1506 if (fAreCollections2MC && !fAreCollections1MC) // the other way around is not supported
1507 GetMCLabelMatchingLevel(jet2->MatchedJet(), jet2, ce1, ce2);
7f76e479 1508 else if ((fAreCollections1MC && fAreCollections2MC) || (!fAreCollections1MC && !fAreCollections2MC))
05077f28 1509 GetSameCollectionsMatchingLevel(jet2->MatchedJet(), jet2, ce1, ce2);
ca5c29fa 1510
05077f28 1511 d = jet2->ClosestJetDistance();
7f76e479 1512 }
1513 else if (jet2->GetMatchingType() == kMCLabel || jet2->GetMatchingType() == kSameCollections) {
05077f28 1514 GetGeometricalMatchingLevel(jet2->MatchedJet(), jet2, d);
1515
1516 ce1 = jet2->MatchedJet()->ClosestJetDistance();
1517 ce2 = jet2->ClosestJetDistance();
1518 }
ca5c29fa 1519
05077f28 1520 fHistCommonEnergy1vsJet1Pt->Fill(ce1, jet2->MatchedJet()->Pt());
1521 fHistCommonEnergy2vsJet2Pt->Fill(ce2, jet2->Pt());
ca5c29fa 1522
05077f28 1523 fHistDistancevsJet1Pt->Fill(d, jet2->MatchedJet()->Pt());
1524 fHistDistancevsJet2Pt->Fill(d, jet2->Pt());
ca5c29fa 1525
05077f28 1526 fHistDistancevsCommonEnergy1->Fill(d, ce1);
507f74bc 1527 fHistDistancevsCommonEnergy2->Fill(d, ce2);
05077f28 1528 fHistCommonEnergy1vsCommonEnergy2->Fill(ce1, ce2);
1529
1530 fHistJet2PtOverJet1PtvsJet2Pt->Fill(jet2->Pt(), jet2->Pt() / jet2->MatchedJet()->Pt());
1531 fHistJet1PtOverJet2PtvsJet1Pt->Fill(jet2->MatchedJet()->Pt(), jet2->MatchedJet()->Pt() / jet2->Pt());
cd6431de 1532
1533 Double_t deta = jet2->MatchedJet()->Eta() - jet2->Eta();
1534 Double_t dphi = jet2->MatchedJet()->Phi() - jet2->Phi();
ca5c29fa 1535 fHistDeltaEtaPhi->Fill(deta, dphi, jet2->Pt());
cd6431de 1536
1537 Double_t dpt = jet2->MatchedJet()->Pt() - jet2->Pt();
7030f36f 1538 fHistDeltaPtvsJet1Pt->Fill(jet2->MatchedJet()->Pt(), dpt);
cfc2ac24 1539 fHistDeltaPtvsJet2Pt->Fill(jet2->Pt(), dpt);
cd6431de 1540
05077f28 1541 fHistDeltaPtvsDistance->Fill(d, dpt);
1542 fHistDeltaPtvsCommonEnergy1->Fill(ce1, dpt);
1543 fHistDeltaPtvsCommonEnergy2->Fill(ce2, dpt);
1544
1545 fHistDeltaPtvsArea1->Fill(jet2->MatchedJet()->Area(), dpt);
1546 fHistDeltaPtvsArea2->Fill(jet2->Area(), dpt);
1547
1548 Double_t darea = jet2->MatchedJet()->Area() - jet2->Area();
1549 fHistDeltaPtvsDeltaArea->Fill(darea, dpt);
c560b734 1550
cd6431de 1551 fHistJet1PtvsJet2Pt->Fill(jet2->MatchedJet()->Pt(), jet2->Pt());
4358e58a 1552
cd6431de 1553 if (!fRhoName.IsNull() || !fRho2Name.IsNull()) {
05077f28 1554 Double_t corrpt1 = jet2->MatchedJet()->Pt() - fRhoVal * jet2->MatchedJet()->Area();
1555 Double_t corrpt2 = jet2->Pt() - fRho2Val * jet2->Area();
1556 Double_t dcorrpt = corrpt1 - corrpt2;
1557 fHistDeltaCorrPtvsJet1Pt->Fill(jet2->MatchedJet()->Pt(), dcorrpt);
1558 fHistDeltaCorrPtvsJet2Pt->Fill(jet2->Pt(), dcorrpt);
1559 fHistDeltaCorrPtvsDistance->Fill(d, dcorrpt);
1560 fHistDeltaCorrPtvsCommonEnergy1->Fill(ce1, dcorrpt);
1561 fHistDeltaCorrPtvsCommonEnergy2->Fill(ce2, dcorrpt);
1562 fHistDeltaCorrPtvsArea1->Fill(jet2->MatchedJet()->Area(), dcorrpt);
1563 fHistDeltaCorrPtvsArea2->Fill(jet2->Area(), dcorrpt);
1564 fHistDeltaCorrPtvsDeltaArea->Fill(darea, dcorrpt);
1565 fHistJet1CorrPtvsJet2CorrPt->Fill(corrpt1, corrpt2);
1566 }
1567
1568 if (fIsEmbedded) {
1569 Double_t dmcpt = jet2->MatchedJet()->MCPt() - jet2->Pt();
1570 fHistDeltaMCPtvsJet1Pt->Fill(jet2->MatchedJet()->MCPt(), dmcpt);
1571 fHistDeltaMCPtvsJet2Pt->Fill(jet2->Pt(), dmcpt);
1572 fHistDeltaMCPtvsDistance->Fill(d, dmcpt);
1573 fHistDeltaMCPtvsCommonEnergy1->Fill(ce1, dmcpt);
1574 fHistDeltaMCPtvsCommonEnergy2->Fill(ce2, dmcpt);
1575 fHistDeltaMCPtvsArea1->Fill(jet2->MatchedJet()->Area(), dmcpt);
1576 fHistDeltaMCPtvsArea2->Fill(jet2->Area(), dmcpt);
1577 fHistDeltaMCPtvsDeltaArea->Fill(darea, dmcpt);
1578 fHistJet1MCPtvsJet2Pt->Fill(jet2->MatchedJet()->MCPt(), jet2->Pt());
cd6431de 1579 }
2bddb6ae 1580 }
2949a2ec 1581 }
99c67c1b 1582 else {
cd6431de 1583 fHistNonMatchedJets2PtArea->Fill(jet2->Area(), jet2->Pt());
1584 fHistMissedJets2PtArea->Fill(jet2->Area(), jet2->Pt());
2949a2ec 1585
cd6431de 1586 if (!fRho2Name.IsNull())
1587 fHistNonMatchedJets2CorrPtArea->Fill(jet2->Area(), jet2->Pt() - fRhoVal * jet2->Area());
1b76c28f 1588 }
1589 }
1590
ca5c29fa 1591 GetSortedArray(indexes, fJets, fRhoVal);
1592
cd6431de 1593 const Int_t nJets1 = fJets->GetEntriesFast();
ca5c29fa 1594 Int_t naccJets1 = 0;
fde82e42 1595 Double_t totalMCPt1 = 0;
1b76c28f 1596
cd6431de 1597 for (Int_t i = 0; i < nJets1; i++) {
1b76c28f 1598
ca5c29fa 1599 AliDebug(2,Form("Processing jet %d", i));
1600 AliEmcalJet* jet1 = static_cast<AliEmcalJet*>(fJets->At(indexes[i]));
1b76c28f 1601
cd6431de 1602 if (!jet1) {
1603 AliError(Form("Could not receive jet %d", i));
2949a2ec 1604 continue;
1b76c28f 1605 }
99c67c1b 1606
cd6431de 1607 if (!AcceptJet(jet1))
1b76c28f 1608 continue;
cd6431de 1609
1610 if (!AcceptBiasJet(jet1))
91bee8bc 1611 continue;
cd6431de 1612
1613 if (jet1->MaxTrackPt() > fMaxTrackPt || jet1->MaxClusterPt() > fMaxClusterPt)
91bee8bc 1614 continue;
cd6431de 1615
1616 if (jet1->Eta() < fJetMinEta || jet1->Eta() > fJetMaxEta || jet1->Phi() < fJetMinPhi || jet1->Phi() > fJetMaxPhi)
91bee8bc 1617 continue;
2949a2ec 1618
cd6431de 1619 if (!jet1->MatchedJet()) {
1620 fHistNonMatchedJets1PtArea->Fill(jet1->Area(), jet1->Pt());
1621 if (!fRhoName.IsNull())
1622 fHistNonMatchedJets1CorrPtArea->Fill(jet1->Area(), jet1->Pt() - fRhoVal * jet1->Area());
99c67c1b 1623 }
1b76c28f 1624
cd6431de 1625 fHistJets1PtArea->Fill(jet1->Area(), jet1->Pt());
1626 fHistJets1PhiEta->Fill(jet1->Eta(), jet1->Phi());
1b76c28f 1627
fde82e42 1628 totalMCPt1 += jet1->MCPt();
1629
ca5c29fa 1630 if (naccJets1 < fNLeadingJets)
1631 fHistLeadingJets1PtArea->Fill(jet1->Area(), jet1->Pt());
1632
1633 if (!fRhoName.IsNull()) {
cd6431de 1634 fHistJets1CorrPtArea->Fill(jet1->Area(), jet1->Pt() - fRhoVal * jet1->Area());
ca5c29fa 1635
1636 if (naccJets1 < fNLeadingJets)
1637 fHistLeadingJets1CorrPtArea->Fill(jet1->Area(), jet1->Pt() - fRhoVal * jet1->Area());
1638 }
1639
63fac07f 1640 if (fTracks) {
1641 for (Int_t it = 0; it < jet1->GetNumberOfTracks(); it++) {
1642 AliVParticle *track1 = jet1->TrackAt(it, fTracks2);
1643 if (track1)
1644 fHistJets1ZvsPt->Fill(track1->Pt() / jet1->Pt(), jet1->Pt());
1645 }
1646 }
1647
1648 if (fCaloClusters) {
1649 for (Int_t ic = 0; ic < jet1->GetNumberOfClusters(); ic++) {
1650 AliVCluster *cluster1 = jet1->ClusterAt(ic, fCaloClusters);
1651
1652 if (cluster1) {
1653 TLorentzVector nPart1;
1654 cluster1->GetMomentum(nPart1, fVertex);
1655 fHistJets2ZvsPt->Fill(nPart1.Et() / jet1->Pt(), jet1->Pt());
1656 }
1657 }
1658 }
1659
1660 fHistJets1NEFvsPt->Fill(jet1->NEF(), jet1->Pt());
1661 fHistJets1CEFvsCEFPt->Fill(1-jet1->NEF(), (1-jet1->NEF())*jet1->Pt());
1662
ca5c29fa 1663 naccJets1++;
1b76c28f 1664 }
6fd5039f 1665
fde82e42 1666 if (fIsEmbedded)
1667 fMCEnergy1vsEnergy2->Fill(totalMCPt1, totalPt2);
1668
6fd5039f 1669 return kTRUE;
1b76c28f 1670}