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