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