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