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