]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/AliAnalysisTaskJetResponseV2.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
- Set Sumw2 for all histos in the histo container
[u/mrichter/AliRoot.git] / PWGJE / AliAnalysisTaskJetResponseV2.cxx
CommitLineData
7a88ef30 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16//
d5ec3aef 17// task compares jets in two or three branches,
7a88ef30 18// written for analysis of jet embedding in HI events
19//
20
21
3c6a60f7 22#include "TChain.h"
23#include "TTree.h"
24#include "TMath.h"
25#include "TH1F.h"
8c483a6b 26#include "TH1I.h"
3c6a60f7 27#include "TH2F.h"
28#include "TH3F.h"
29#include "THnSparse.h"
30#include "TCanvas.h"
31
32#include "AliLog.h"
33
34#include "AliAnalysisTask.h"
35#include "AliAnalysisManager.h"
36
37#include "AliVEvent.h"
38#include "AliESDEvent.h"
39#include "AliESDInputHandler.h"
40#include "AliCentrality.h"
41#include "AliAnalysisHelperJetTasks.h"
42#include "AliInputEventHandler.h"
d314de7b 43#include "AliAODJetEventBackground.h"
31b9d515 44#include "AliAnalysisTaskFastEmbedding.h"
3c6a60f7 45
46#include "AliAODEvent.h"
47#include "AliAODJet.h"
53f9653b 48#include "AliAODHandler.h"
3c6a60f7 49
50#include "AliAnalysisTaskJetResponseV2.h"
51
c64cb1f6 52using std::cout;
53using std::endl;
54
3c6a60f7 55ClassImp(AliAnalysisTaskJetResponseV2)
56
57AliAnalysisTaskJetResponseV2::AliAnalysisTaskJetResponseV2() :
31b9d515 58AliAnalysisTaskSE(),
d5ec3aef 59 fESD(0x0),
60 fAOD(0x0),
f34e09e6 61 fAODOut(0x0),
53f9653b 62 fAODExtension(0x0),
63 fNonStdFile(""),
d5ec3aef 64 fBackgroundBranch(""),
65 fIsPbPb(kTRUE),
66 fOfflineTrgMask(AliVEvent::kAny),
67 fMinContribVtx(1),
68 fVtxZMin(-8.),
69 fVtxZMax(8.),
70 fEvtClassMin(0),
71 fEvtClassMax(4),
72 fCentMin(0.),
73 fCentMax(100.),
74 fNInputTracksMin(0),
75 fNInputTracksMax(-1),
76 fJetEtaMin(-.5),
77 fJetEtaMax(.5),
78 fJetPtMin(20.),
79 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
80 fJetPtFractionMin(0.5),
81 fNMatchJets(4),
82 fMatchMaxDist(0.8),
83 fKeepJets(kFALSE),
84 fkNbranches(2),
85 fkEvtClasses(12),
86 fOutputList(0x0),
87 fbEvent(kTRUE),
88 fbJetsMismatch1(kTRUE),
89 fbJetsMismatch2(kTRUE),
90 fbJetsRp(kTRUE),
91 fbJetsDeltaPt(kTRUE),
92 fbJetsEta(kTRUE),
93 fbJetsPhi(kTRUE),
94 fbJetsArea(kTRUE),
95 fbJets3Branches(kFALSE),
96 fbJetsBeforeCut1(kTRUE),
97 fbJetsBeforeCut2(kTRUE),
98 fHistEvtSelection(0x0),
99 fHistJetSelection(0x0),
100 fh2JetSelection(0x0),
101 fhnEvent(0x0),
102 fhnJetsMismatch1(0x0),
103 fhnJetsMismatch2(0x0),
104 fhnJetsRp(0x0),
105 fhnJetsDeltaPt(0x0),
106 fhnJetsEta(0x0),
107 fhnJetsPhi(0x0),
108 fhnJetsArea(0x0),
109 fhnJets3Branches(0x0),
110 fhnJetsBeforeCut1(0x0),
111 fhnJetsBeforeCut2(0x0)
3c6a60f7 112{
d5ec3aef 113 // default Constructor
3c6a60f7 114
d5ec3aef 115 fJetBranchName[0] = "";
116 fJetBranchName[1] = "";
117 fJetBranchName[2] = "";
3c6a60f7 118
d5ec3aef 119 fListJets[0] = new TList;
120 fListJets[1] = new TList;
121 fListJets[2] = new TList;
3c6a60f7 122}
123
124AliAnalysisTaskJetResponseV2::AliAnalysisTaskJetResponseV2(const char *name) :
d5ec3aef 125 AliAnalysisTaskSE(name),
126 fESD(0x0),
127 fAOD(0x0),
f34e09e6 128 fAODOut(0x0),
53f9653b 129 fAODExtension(0x0),
130 fNonStdFile(""),
d5ec3aef 131 fBackgroundBranch(""),
132 fIsPbPb(kTRUE),
133 fOfflineTrgMask(AliVEvent::kAny),
134 fMinContribVtx(1),
135 fVtxZMin(-8.),
136 fVtxZMax(8.),
137 fEvtClassMin(0),
138 fEvtClassMax(4),
139 fCentMin(0.),
140 fCentMax(100.),
141 fNInputTracksMin(0),
142 fNInputTracksMax(-1),
143 fJetEtaMin(-.5),
144 fJetEtaMax(.5),
145 fJetPtMin(20.),
146 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
147 fJetPtFractionMin(0.5),
148 fNMatchJets(4),
149 fMatchMaxDist(0.8),
150 fKeepJets(kFALSE),
151 fkNbranches(2),
152 fkEvtClasses(12),
153 fOutputList(0x0),
154 fbEvent(kTRUE),
155 fbJetsMismatch1(kTRUE),
156 fbJetsMismatch2(kTRUE),
157 fbJetsRp(kTRUE),
158 fbJetsDeltaPt(kTRUE),
159 fbJetsEta(kTRUE),
160 fbJetsPhi(kTRUE),
161 fbJetsArea(kTRUE),
162 fbJets3Branches(kFALSE),
163 fbJetsBeforeCut1(kTRUE),
164 fbJetsBeforeCut2(kTRUE),
165 fHistEvtSelection(0x0),
166 fHistJetSelection(0x0),
167 fh2JetSelection(0x0),
168 fhnEvent(0x0),
169 fhnJetsMismatch1(0x0),
170 fhnJetsMismatch2(0x0),
171 fhnJetsRp(0x0),
172 fhnJetsDeltaPt(0x0),
173 fhnJetsEta(0x0),
174 fhnJetsPhi(0x0),
175 fhnJetsArea(0x0),
176 fhnJets3Branches(0x0),
177 fhnJetsBeforeCut1(0x0),
178 fhnJetsBeforeCut2(0x0)
3c6a60f7 179{
d5ec3aef 180 // Constructor
3c6a60f7 181
d5ec3aef 182 fJetBranchName[0] = "";
183 fJetBranchName[1] = "";
184 fJetBranchName[2] = "";
3c6a60f7 185
d5ec3aef 186 fListJets[0] = new TList;
187 fListJets[1] = new TList;
188 fListJets[2] = new TList;
3c6a60f7 189
d5ec3aef 190 DefineOutput(1, TList::Class());
3c6a60f7 191}
192
193AliAnalysisTaskJetResponseV2::~AliAnalysisTaskJetResponseV2()
194{
d5ec3aef 195 delete fListJets[0];
196 delete fListJets[1];
197 delete fListJets[2];
3c6a60f7 198}
199
d5ec3aef 200void AliAnalysisTaskJetResponseV2::SetBranchNames(const TString &branch1, const TString &branch2, const TString &branch3)
3c6a60f7 201{
d5ec3aef 202 fJetBranchName[0] = branch1;
203 fJetBranchName[1] = branch2;
204 fJetBranchName[2] = branch3;
205
34d5296f 206 if(strlen(fJetBranchName[2].Data()) ) {
d5ec3aef 207 fkNbranches = 3;
34d5296f 208 fbJets3Branches = kTRUE;
209 }
3c6a60f7 210}
211
212void AliAnalysisTaskJetResponseV2::Init()
213{
214
d5ec3aef 215 // check for jet branches
216 if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
217 AliError("Jet branch name not set.");
218 }
3c6a60f7 219
220}
221
222void AliAnalysisTaskJetResponseV2::UserCreateOutputObjects()
223{
d5ec3aef 224 // Create histograms
225 // Called once
226 OpenFile(1);
227 if(!fOutputList) fOutputList = new TList;
228 fOutputList->SetOwner(kTRUE);
229
230 Bool_t oldStatus = TH1::AddDirectoryStatus();
231 TH1::AddDirectory(kFALSE);
232
233
234 fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
235 fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
236 fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
237 fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
238 fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
239 fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
240 fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
241
242 fHistJetSelection = new TH1I("fHistJetSelection", "jet selection", 8, -0.5, 7.5);
243 fHistJetSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
244 fHistJetSelection->GetXaxis()->SetBinLabel(2,"probes IN");
245 fHistJetSelection->GetXaxis()->SetBinLabel(3,"no matching jet");
246 fHistJetSelection->GetXaxis()->SetBinLabel(4,"not in list");
247 fHistJetSelection->GetXaxis()->SetBinLabel(5,"fraction cut");
248 fHistJetSelection->GetXaxis()->SetBinLabel(6,"acceptance cut");
249 fHistJetSelection->GetXaxis()->SetBinLabel(7,"p_{T} cut");
250 fHistJetSelection->GetXaxis()->SetBinLabel(8,"trigger exclude mask");
251
252 fh2JetSelection = new TH2F("fh2JetSelection", "jet selection", 8, -0.5, 7.5,100,0.,200.);
253 fh2JetSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
254 fh2JetSelection->GetXaxis()->SetBinLabel(2,"probes IN");
255 fh2JetSelection->GetXaxis()->SetBinLabel(3,"no matching jet");
256 fh2JetSelection->GetXaxis()->SetBinLabel(4,"not in list");
257 fh2JetSelection->GetXaxis()->SetBinLabel(5,"fraction cut");
258 fh2JetSelection->GetXaxis()->SetBinLabel(6,"acceptance cut");
259 fh2JetSelection->GetXaxis()->SetBinLabel(7,"p_{T} cut");
260 fh2JetSelection->GetXaxis()->SetBinLabel(8,"trigger exclude mask");
261
262
263 UInt_t entries = 0; // bit coded, see GetDimParams() below
264 UInt_t opt = 0; // bit coded, default (0) or high resolution (1)
265
266 if(fbEvent){
267 entries = 1<<0 | 1<<1 | 1<<2 | 1<<26; // cent : nInpTrks : rp psi : pT hard bin
268 opt = 1<<0 | 1<<1; // centrality and nInpTrks in high resolution
269 fhnEvent = NewTHnSparseF("fhnEvent", entries, opt);
270 }
31b9d515 271
d5ec3aef 272 if(fbJetsMismatch1){ // real mismatch (no related rec jet found)
273 // cent : nInpTrks : rp bins : probe pt : probe eta : probe phi : probe area : pT hard bin
274 entries = 1<<0 | 1<<1 | 1<<3 | 1<<6 | 1<<8 | 1<<10 | 1<<12 | 1<<26;
275 opt = 1<<6 | 1<<8 | 1<<10;
276 fhnJetsMismatch1 = NewTHnSparseF("fhnJetsMismatch1", entries, opt);
277 }
278
279 if(fbJetsMismatch2){ // acceptance + fraction cut
280 // cent : nInpTrks : rp bins : jetPt(2x) : jetEta(2x) : deltaEta : deltaR : fraction : pT hard bin
281 entries = 1<<0 | 1<<1 | 1<<3 | 1<<6 | 1<<7 | 1<<8 | 1<<9 | 1<<15 | 1<<17 | 1<<19 | 1<<26;
282 opt = 1<<6 | 1<<7 | 1<<8 | 1<<9;
283 fhnJetsMismatch2 = NewTHnSparseF("fhnJetsMismatch2", entries, opt);
284
285 }
d314de7b 286
d314de7b 287
d5ec3aef 288 if(fbJetsRp){
289 // cent : nInpTrks : rp bins : rp wrt jet(2x) : probe pT : probe area: deltaPt : rp phi : rho : correction for RP : local rho
290 /*
291 entries = 1<<0 | 1<<1 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<12 | 1<<14 | 1<<2 | 1<<22 | 1<<23 | 1<<24 | 1<<25;
292 opt = 1<<4 | 1<<5;*/
293 // cent : nInpTrks : rp bins : rp wrt jet(2x) : probe pT : deltaPt : rp phi : rho : correction for RP | pT hard bin
294 entries = 1<<0 | 1<<1 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<14 | 1<<2 | 1<<22 | 1<<23 | 1<<26;
295 opt = 1<<4 | 1<<5;
296 fhnJetsRp = NewTHnSparseF("fhnJetsRp", entries, opt);
297 }
298
299 // cent : nInpTrks : rp bins: deltaPt : jetPt(2x) : deltaArea : pT hard bin (hr delta pt)
300 if(fbJetsDeltaPt){
301 entries = 1<<0 | 1<<1 | 1<<3 | 1<<14 | 1<<6 | 1<<7 | 1<<18 | 1<<26;
302 opt = 1<<1 | 1<<14 | 1<<6 | 1<<7 | 1<<18;
303 fhnJetsDeltaPt = NewTHnSparseF("fhnJetsDeltaPt", entries, opt);
304 }
305
306 // cent : nInpTrks : rp bins : deltaPt : jetPt(2x) : deltaR : deltaEta : jetEta(2x) : pT hard bin (hr for eta)
307 if(fbJetsEta){
308 entries = 1<<0 | 1<<1 | 1<<3 | 1<<14 | 1<<6 | 1<<7 | 1<<17 | 1<<15 | 1<<8 | 1<<9 | 1<<26;
309 opt = 1<<15 | 1<<8 | 1<<9;
310 fhnJetsEta = NewTHnSparseF("fhnJetsEta", entries, opt);
311 }
312
313 // cent : nInpTrks : rp bins : jetPt(2x) : jetPhi(2x) : deltaPt : deltaPhi : pT hard bin
314 if(fbJetsPhi){
315 entries = 1<<0 | 1<<1 | 1<<3 | 1<<6 | 1<<7 | 1<<10 | 1<<11 | 1<<14 | 1<<16 | 1<<26;
316 opt = 1<<10 | 1<<11;
317 fhnJetsPhi = NewTHnSparseF("fhnJetsPhi", entries, opt);
318 }
319
320 // cent : nInpTrks : rp bins : deltaArea : jetArea(2x) : deltaR : fraction : distance next rec jet : pT next jet : deltaPt : jetPt(2x) : pT hard bin (hr for area)
321 if(fbJetsArea){
322 entries = 1<<0 | 1<<1 | 1<<3 | 1<<18 | 1<<12 | 1<<13 | 1<<17 | 1<<19 | 1<<20 | 1<<21 | 1<<14 | 1<<6 | 1<<7 | 1<<26;
323 opt = 1<<18 | 1<<12 | 1<<13;
324 fhnJetsArea = NewTHnSparseF("fhnJetsArea", entries, opt);
325 }
326
fedb5a47 327 // cent : nInpTrks : jetPt(3x) : deltaPt : delta : jetArea(3x) : fraction(2x) : deltaR(1x) : pT hard bin
d5ec3aef 328 if(fbJets3Branches){
fedb5a47 329 entries = 1<<0 | 1<<1 | 1<<6 | 1<<7 | 1<<27 | 1<<14 | 1<<28 | 1<<12 | 1<<13 | 1<<29 | 1<<19 | 1<<30 | 1<<17 | 1<<26;
b5fee0f9 330 opt = 1<<6 | 1<<7 | 1<<27 | 1<<14 | 1<<28;
d5ec3aef 331 fhnJets3Branches = NewTHnSparseF("fhnJets3Branches", entries, opt);
332 }
333
334
335
336
337 //before cut
338
339 // cent : nInpTrks : rp bins : fraction : jetPt(2x) : jetEta(2x) : jetPhi(2x) (low resolution) (with fraction, eta, phi, pt cuts possible)
340 if(fbJetsBeforeCut1){
341 entries = 1<<0 | 1<<1 | 1<<3 | 1<<19 | 1<<6 | 1<<7 | 1<<8 | 1<<9 | 1<<10 | 1<<11;
342 opt = 0;
343 fhnJetsBeforeCut1 = NewTHnSparseF("fhnJetsBeforeCut1", entries, opt);
344 }
345
346 // cent : nInpTrks : rp bins : deltaPt : jetPt(2x) : deltaR : deltaEta : jetEta(2x) (low resolution)
347 if(fbJetsBeforeCut2){
348 entries = 1<<0 | 1<<1 | 1<<3 | 1<<14 | 1<<6 | 1<<7 | 1<<17 | 1<<15 | 1<<8 | 1<<9;
349 opt = 0;
350 fhnJetsBeforeCut2 = NewTHnSparseF("fhnJetsBeforeCut2", entries, opt);
351 }
352
353 fOutputList->Add(fHistEvtSelection);
354 fOutputList->Add(fHistJetSelection);
355 fOutputList->Add(fh2JetSelection);
356 if(fbEvent) fOutputList->Add(fhnEvent);
357 if(fbJetsMismatch1) fOutputList->Add(fhnJetsMismatch1);
358 if(fbJetsMismatch2) fOutputList->Add(fhnJetsMismatch2);
359 if(fbJetsRp) fOutputList->Add(fhnJetsRp);
360 if(fbJetsDeltaPt) fOutputList->Add(fhnJetsDeltaPt);
361 if(fbJetsEta) fOutputList->Add(fhnJetsEta);
362 if(fbJetsPhi) fOutputList->Add(fhnJetsPhi);
363 if(fbJetsArea) fOutputList->Add(fhnJetsArea);
364 if(fbJets3Branches) fOutputList->Add(fhnJets3Branches);
365 if(fbJetsBeforeCut1) fOutputList->Add(fhnJetsBeforeCut1);
366 if(fbJetsBeforeCut2) fOutputList->Add(fhnJetsBeforeCut2);
367
368 // =========== Switch on Sumw2 for all histos ===========
369 for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
370 TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
371 if (h1){
372 h1->Sumw2();
373 continue;
374 }
375 THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
376 if (hn){
377 hn->Sumw2();
378 }
379 }
380 TH1::AddDirectory(oldStatus);
381
382 PostData(1, fOutputList);
3c6a60f7 383}
384
385void AliAnalysisTaskJetResponseV2::UserExec(Option_t *)
386{
d5ec3aef 387 // load events, apply event cuts, then compare leading jets
388
389 if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
390 AliError("Jet branch name not set.");
391 return;
392 }
393
394 fESD=dynamic_cast<AliESDEvent*>(InputEvent());
395 if (!fESD) {
396 AliError("ESD not available");
f34e09e6 397
d5ec3aef 398 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
f34e09e6 399 // assume that the AOD is in the general output...
400 fAODOut = AODEvent();
d5ec3aef 401 } else {
402 fAOD = dynamic_cast<AliAODEvent*>(AODEvent());
403 }
404 if (!fAOD) {
405 AliError("AOD not available");
406 return;
407 }
408
53f9653b 409 if(fNonStdFile.Length()!=0){
410 // case that we have an AOD extension we need can fetch the jets from the extended output
411 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
412 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
413 if(!fAODExtension){
414 if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
415 }
416 }
417
d5ec3aef 418 // -- event selection --
419 fHistEvtSelection->Fill(1); // number of events before event selection
420
421 // physics selection
422 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
423 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
424 if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
425 if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
426 fHistEvtSelection->Fill(2);
427 PostData(1, fOutputList);
428 return;
429 }
430
431 // vertex selection
432 AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
433 Int_t nTracksPrim = primVtx->GetNContributors();
434 if ((nTracksPrim < fMinContribVtx) ||
435 (primVtx->GetZ() < fVtxZMin) ||
436 (primVtx->GetZ() > fVtxZMax) ){
437 if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
438 fHistEvtSelection->Fill(3);
439 PostData(1, fOutputList);
440 return;
441 }
442
443 // event class selection (from jet helper task)
444 Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
445 if(fDebug) Printf("Event class %d", eventClass);
446 if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
447 fHistEvtSelection->Fill(4);
448 PostData(1, fOutputList);
449 return;
450 }
451
452 // centrality selection
453 AliCentrality *cent = 0x0;
454 Float_t centValue = 0.;
455 if(fESD) cent = fESD->GetCentrality();
456 if(cent) centValue = cent->GetCentralityPercentile("V0M");
457 if(fDebug) printf("centrality: %f\n", centValue);
458 if (centValue < fCentMin || centValue > fCentMax){
459 fHistEvtSelection->Fill(4);
460 PostData(1, fOutputList);
461 return;
462 }
463
464
465 // multiplicity due to input tracks
466 Int_t nInputTracks = GetNInputTracks();
467
468 if (nInputTracks < fNInputTracksMin || (fNInputTracksMax > -1 && nInputTracks > fNInputTracksMax)){
469 fHistEvtSelection->Fill(5);
470 PostData(1, fOutputList);
471 return;
472 }
31b9d515 473
474
d5ec3aef 475 fHistEvtSelection->Fill(0); // accepted events
476 // -- end event selection --
477
478 // pt hard
479 Double_t pthard = AliAnalysisTaskFastEmbedding::GetPtHard();
480 Int_t pthardbin = GetPtHardBin(pthard);
481
482 // reaction plane
483 Double_t rp = AliAnalysisHelperJetTasks::ReactionPlane(kFALSE);
484 if(fbEvent){
485 Double_t eventEntries[3] = { (Double_t)centValue, (Double_t)nInputTracks, rp };
486 fhnEvent->Fill(eventEntries);
487 }
488
489
490 // get background
491 AliAODJetEventBackground* externalBackground = 0;
492 if(!externalBackground&&fBackgroundBranch.Length()){
493 externalBackground = (AliAODJetEventBackground*)(fAOD->FindListObject(fBackgroundBranch.Data()));
f34e09e6 494 if(!externalBackground && fAODOut) externalBackground = (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch.Data()));
53f9653b 495 if(!externalBackground && fAODExtension) externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
d5ec3aef 496 //if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
497 }
498 Float_t rho = 0;
53f9653b 499 if(externalBackground) rho = externalBackground->GetBackground(0);
d5ec3aef 500
501
502 // fetch jets
503 TClonesArray *aodJets[3];
504 aodJets[0] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[0].Data())); // in general: embedded jet
f34e09e6 505 if(!aodJets[0] && fAODOut) aodJets[0] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[0].Data())); // in general: embedded jet
53f9653b 506 if(!aodJets[0] && fAODExtension) aodJets[0] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[0].Data()));
d5ec3aef 507 aodJets[1] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[1].Data())); // in general: embedded jet + UE version1
f34e09e6 508 if(!aodJets[1] && fAODOut) aodJets[1] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[1].Data())); // in general: embedded jet
53f9653b 509 if(!aodJets[1] && fAODExtension) aodJets[1] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[1].Data()));
b5fee0f9 510 if( strlen(fJetBranchName[2].Data()) ) {
d5ec3aef 511 aodJets[2] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[2].Data())); // in general: embedded jet + UE version2
f34e09e6 512 if(!aodJets[2] && fAODOut) aodJets[2] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[2].Data())); // in general: embedded jet
53f9653b 513 if(!aodJets[2] && fAODExtension) aodJets[2] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[2].Data()));
514 if(aodJets[2]) fkNbranches=3;
515 if(fDebug>10) printf("3rd branch: %s\n",fJetBranchName[2].Data());
b5fee0f9 516 }
d5ec3aef 517
53f9653b 518 if(fDebug>10) printf("fkNbranches %d\n",fkNbranches);
d5ec3aef 519 for (Int_t iJetType = 0; iJetType < fkNbranches; iJetType++) {
520 fListJets[iJetType]->Clear();
53f9653b 521 if (!aodJets[iJetType]) {
522 if(fDebug) Printf("%s: no jet branch\n",fJetBranchName[iJetType].Data());
523 continue;
524 }
d5ec3aef 525 if(fDebug) Printf("%s: %d jets",fJetBranchName[iJetType].Data(),aodJets[iJetType]->GetEntriesFast());
526
527 for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
528 AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
529 if (jet) fListJets[iJetType]->Add(jet);
530 }
531 fListJets[iJetType]->Sort();
532 }
31b9d515 533
534
d5ec3aef 535 // jet matching
536 static TArrayI aMatchIndex(fListJets[0]->GetEntries());
537 static TArrayF aPtFraction(fListJets[0]->GetEntries());
538 if(aMatchIndex.GetSize()<fListJets[0]->GetEntries()) aMatchIndex.Set(fListJets[0]->GetEntries());
539 if(aPtFraction.GetSize()<fListJets[0]->GetEntries()) aPtFraction.Set(fListJets[0]->GetEntries());
540
541 // stores matched jets in 'aMatchIndex' and fraction of pT in 'aPtFraction'
542 AliAnalysisHelperJetTasks::GetJetMatching(fListJets[0], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[0]->GetEntries()),
543 fListJets[1], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[1]->GetEntries()),
544 aMatchIndex, aPtFraction, fDebug, fMatchMaxDist, fIsPbPb?1:2);
545 static TArrayI aMatchIndexv2(fListJets[0]->GetEntries());
546 static TArrayF aPtFractionv2(fListJets[0]->GetEntries());
b5fee0f9 547 if(aMatchIndexv2.GetSize()<fListJets[0]->GetEntries()) aMatchIndexv2.Set(fListJets[0]->GetEntries());
548 if(aPtFractionv2.GetSize()<fListJets[0]->GetEntries()) aPtFractionv2.Set(fListJets[0]->GetEntries());
d5ec3aef 549 if( strlen(fJetBranchName[2].Data()) ) {
b5fee0f9 550 AliAnalysisHelperJetTasks::GetJetMatching(fListJets[0], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[0]->GetEntries()),
551 fListJets[2], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[2]->GetEntries()),
552 aMatchIndexv2, aPtFractionv2, fDebug, fMatchMaxDist, fIsPbPb?1:2);
d5ec3aef 553 }
3c6a60f7 554
d5ec3aef 555 // loop over matched jets
556 Int_t ir = -1; // index of matched reconstruced jet
557 Int_t ir2 = -1; // index of matched reconstruced jet
558 Float_t fraction = -1.;
559 Float_t fraction2 = -1.;
560 AliAODJet *jet[3] = { 0x0, 0x0, 0x0 };
561 Float_t jetEta[3] = { -990., -990., -990. };
562 Float_t jetPhi[3] = { -990., -990., -990. };
563 Float_t jetPt[3] = { -990., -990., -990. };
564 Float_t jetArea[3] = { -990., -990., -990. };
565 Float_t rpJet[3] = { -990., -990., -990. };
566 Int_t rpBin = -990;
3c6a60f7 567
d5ec3aef 568 for(Int_t ig=0; ig<fListJets[0]->GetEntries(); ++ig){
569 ir = aMatchIndex[ig];
b5fee0f9 570 if(aMatchIndexv2.GetSize()>ig) ir2 = aMatchIndexv2[ig];
571 else ir2 =0;
d5ec3aef 572
573 //fetch jets
574 jet[0] = (AliAODJet*)(fListJets[0]->At(ig));
575 if(ir>=0) jet[1] = (AliAODJet*)(fListJets[1]->At(ir));
576 else jet[1] = 0x0;
577 if(ir2>=0) jet[2] = (AliAODJet*)(fListJets[2]->At(ir2));
578 else jet[2] = 0x0;
579
580 for(Int_t i=0; i<fkNbranches; ++i){
581 if(!jet[i]){
582 jetEta[i] = -990;
583 jetPhi[i] = -990.;
584 jetPt[i] = -990.;
585 jetArea[i] = -990.;
586 rpJet[i] = -990.;
587 if(i==1) rpBin = -990;
588 }
589 else {
590 jetEta[i] = jet[i]->Eta();
591 jetPhi[i] = jet[i]->Phi();
592 jetPt[i] = GetPt(jet[i], i);
593 jetArea[i] = jet[i]->EffectiveAreaCharged();
594 rpJet[i] = TVector2::Phi_mpi_pi(rp-jetPhi[i]);
595 if(i==1) rpBin = AliAnalysisHelperJetTasks::GetPhiBin(TVector2::Phi_mpi_pi(rp-jetPhi[i]), 3);
d314de7b 596 }
d5ec3aef 597 }
598 fraction = aPtFraction[ig];
b5fee0f9 599 if(aPtFractionv2.GetSize()>ig) fraction2 = aPtFractionv2[ig];
600 else fraction2 = 0.;
31b9d515 601
d5ec3aef 602 // jet statistics
603 fHistJetSelection->Fill(1); // all probe jets
604 if(jet[0]) fh2JetSelection->Fill(1.,jetPt[0]); // all probe jets
31b9d515 605
d5ec3aef 606 if(ir<0){
607 fHistJetSelection->Fill(2);
608 if(jet[0]) fh2JetSelection->Fill(2.,jetPt[0]);
31b9d515 609
d5ec3aef 610 if(fbJetsMismatch1){
611 if(!jet[0]) continue;
612 Double_t jetEntriesMismatch1[7] = {
613 (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
614 (Double_t)jetPt[0], (Double_t)jetEta[0], (Double_t)jetPhi[0], (Double_t)jetArea[0]
615 };
616 fhnJetsMismatch1->Fill(jetEntriesMismatch1);
31b9d515 617 }
d5ec3aef 618
619 continue;
620 }
31b9d515 621
d5ec3aef 622 if(!jet[0] || !jet[1]){
623 fHistJetSelection->Fill(3);
624 if(jet[0]) fh2JetSelection->Fill(3.,jetPt[0]);
625 continue;
626 }
31b9d515 627
d5ec3aef 628 // look for distance to next rec jet
629 Float_t distNextJet = -0.01; // no neighbor
630 Float_t ptNextJet = -1.;
631 for(Int_t r=0; r<fListJets[1]->GetEntries(); ++r){
632 if(r==ir) continue;
633 Float_t tmpDeltaR = jet[1]->DeltaR((AliAODJet*)fListJets[1]->At(r));
634 if(distNextJet<0. || distNextJet>tmpDeltaR){
635 distNextJet = tmpDeltaR;
636 if(fKeepJets) ptNextJet = ((AliAODJet*)fListJets[1]->At(r))->GetPtSubtracted(0);
637 else ptNextJet = ((AliAODJet*)fListJets[1]->At(r))->Pt();
d314de7b 638 }
d5ec3aef 639 }
31b9d515 640
641
642
d5ec3aef 643 //Float_t localRho = jetArea[1]>0. ? (jetPt[1]+rho*jetArea[1] - jetPt[0]) / jetArea[1] : 0.;
644 //Float_t relRho = rho>0. ? localRho / rho : 0.;
31b9d515 645
646
d5ec3aef 647 // calculate parameters of associated jets
648 /* from Leticia, kT clusterizer
649 Float_t par0[4] = { 0.00409, 0.01229, 0.05, 0.26 };
650 Float_t par1[4] = { -2.97035e-03, -2.03182e-03, -1.25702e-03, -9.95107e-04 };
651 Float_t par2[4] = { 1.02865e-01, 1.49039e-01, 1.53910e-01, 1.51109e-01 };
652 */
653 // own, from embedded tracks
654 Float_t par0[4] = { 0.02841, 0.05039, 0.09092, 0.24089 };
655 Float_t par1[4] = { -4.26725e-04, -1.15273e-03, -1.56827e-03, -3.08003e-03 };
656 Float_t par2[4] = { 4.95415e-02, 9.79538e-02, 1.32814e-01, 1.71743e-01 };
31b9d515 657
d5ec3aef 658 Float_t rpCorr = 0.;
31b9d515 659
d5ec3aef 660 if(eventClass>0&&eventClass<4){
661 rpCorr = par0[eventClass-1] + par1[eventClass-1] * 2*TMath::Cos(rpJet[1]) + par2[eventClass-1] * 2*TMath::Cos(2*rpJet[1]);
662 rpCorr *= rho * jetArea[1];
663 }
664
665 Float_t delta = jetPt[2]-jetPt[1];
666 Float_t deltaPt = jetPt[1]-jetPt[0];
667 Float_t deltaEta = jetEta[1]-jetEta[0];
668 Float_t deltaPhi = TVector2::Phi_mpi_pi(jetPhi[1]-jetPhi[0]);
669 Float_t deltaR = TMath::Sqrt(deltaEta*deltaEta+deltaPhi*deltaPhi);
670 Float_t deltaArea = jetArea[1]-jetArea[0];
31b9d515 671
672
d5ec3aef 673 // fill thnsparse before acceptance cut
674 if(fbJetsBeforeCut1){
675 Double_t jetBeforeCutEntries1[10] = {
676 (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
677 (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetEta[0], (Double_t)jetEta[1],
678 (Double_t)jetPhi[0], (Double_t)jetPhi[1], (Double_t)fraction
679 };
680 fhnJetsBeforeCut1->Fill(jetBeforeCutEntries1);
681 }
31b9d515 682
d5ec3aef 683 if(fbJetsBeforeCut2){
684 Double_t jetBeforeCutEntries2[10] = {
685 (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
686 (Double_t)jetPt[0], (Double_t)jetPt[1],
687 (Double_t)jetEta[0], (Double_t)jetEta[1],
688 (Double_t)deltaPt, (Double_t)deltaEta, (Double_t)deltaR
689 };
690 fhnJetsBeforeCut2->Fill(jetBeforeCutEntries2);
691 }
31b9d515 692
d5ec3aef 693 // jet selection
694 Bool_t jetAccepted = kTRUE;
695 // minimum fraction required
696 if(fraction<fJetPtFractionMin){
697 fHistJetSelection->Fill(4);
698 fh2JetSelection->Fill(4.,jetPt[0]);
699 jetAccepted = kFALSE;
700 }
31b9d515 701
d5ec3aef 702 if(jetAccepted){
703 // jet acceptance + minimum pT check
704 if(jetEta[0]>fJetEtaMax || jetEta[0]<fJetEtaMin ||
705 jetEta[1]>fJetEtaMax || jetEta[1]<fJetEtaMin){
31b9d515 706
d5ec3aef 707 if(fDebug){
708 Printf("Jet not in eta acceptance.");
709 Printf("[0]: jet %d eta %.2f", ig, jetEta[0]);
710 Printf("[1]: jet %d eta %.2f", ir, jetEta[1]);
711 }
712 fHistJetSelection->Fill(5);
713 fh2JetSelection->Fill(5.,jetPt[0]);
714 jetAccepted = kFALSE;
31b9d515 715 }
d5ec3aef 716 }
717 if(jetAccepted){
718 if(jetPt[0] < fJetPtMin){
719 if(fDebug) Printf("Jet %d (pT %.1f GeV/c) has less than required pT.", ir, jetPt[0]);
720 fHistJetSelection->Fill(6);
721 fh2JetSelection->Fill(6.,jetPt[0]);
722 jetAccepted = kFALSE;
31b9d515 723 }
d5ec3aef 724 }
31b9d515 725
d5ec3aef 726 if(jetAccepted){
727 if(jet[1]->Trigger()&fJetTriggerExcludeMask){
728 fHistJetSelection->Fill(7);
729 fh2JetSelection->Fill(7.,jetPt[0]);
730 jetAccepted = kFALSE;
31b9d515 731 }
d5ec3aef 732
733 }
31b9d515 734
d5ec3aef 735 if(!jetAccepted){
736 if(fbJetsMismatch2){
737 Double_t jetEntriesMismatch2[10] = {
738 (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
739 (Double_t)jetPt[0], (Double_t)jetPt[1],
740 (Double_t)jetEta[0], (Double_t)jetEta[1],
741 (Double_t)deltaEta, (Double_t)deltaR,
742 (Double_t)fraction
743 };
744 fhnJetsMismatch2->Fill(jetEntriesMismatch2);
31b9d515 745 }
d5ec3aef 746 continue;
747 }
31b9d515 748
d5ec3aef 749 // all accepted jets
750 fHistJetSelection->Fill(0);
751 fh2JetSelection->Fill(0.,jetPt[0]);
31b9d515 752
d5ec3aef 753 // fill thnsparse
754 if(fbJetsRp){
755 Double_t jetEntriesRp[11] = {
756 (Double_t)centValue, (Double_t)nInputTracks, (Double_t) rp,
757 (Double_t)rpBin, (Double_t)rpJet[0], (Double_t)rpJet[1],
758 (Double_t)jetPt[0], (Double_t)deltaPt, (Double_t)rho, (Double_t)rpCorr,
759 (Double_t)pthardbin
760 };
761 fhnJetsRp->Fill(jetEntriesRp);
762 }
31b9d515 763
d5ec3aef 764 if(fbJetsDeltaPt){
765 Double_t jetEntriesDeltaPt[8] = {
766 (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
767 (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)deltaPt, (Double_t)deltaArea,
768 (Double_t)pthardbin
769 };
770 fhnJetsDeltaPt->Fill(jetEntriesDeltaPt);
771 }
31b9d515 772
d5ec3aef 773 if(fbJetsEta){
774 Double_t jetEntriesEta[11] = {
775 (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
776 (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetEta[0], (Double_t)jetEta[1],
777 (Double_t)deltaPt, (Double_t)deltaEta, (Double_t)deltaR,
778 (Double_t)pthardbin
779 };
780 fhnJetsEta->Fill(jetEntriesEta);
781 }
31b9d515 782
d5ec3aef 783 if(fbJetsPhi){
784 Double_t jetEntriesPhi[10] = {
785 (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
786 (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetPhi[0], (Double_t)jetPhi[1],
787 (Double_t)deltaPt, (Double_t)deltaPhi,
788 (Double_t)pthardbin
789 };
790 fhnJetsPhi->Fill(jetEntriesPhi);
791 }
31b9d515 792
d5ec3aef 793 if(fbJetsArea){
794 Double_t jetEntriesArea[14] = {
795 (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
796 (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetArea[0], (Double_t)jetArea[1],
797 (Double_t)deltaPt, (Double_t)deltaR, (Double_t)deltaArea,
798 (Double_t)fraction, (Double_t)distNextJet, (Double_t)ptNextJet,
799 (Double_t)pthardbin
800 };
801 fhnJetsArea->Fill(jetEntriesArea);
802 }
803
804 if(fbJets3Branches){
fedb5a47 805 Double_t jetEntries3Branches[14] = {
d5ec3aef 806 (Double_t)centValue, (Double_t)nInputTracks,
b5fee0f9 807 (Double_t)jetPt[0], (Double_t)jetPt[1],
808 (Double_t)jetArea[0], (Double_t)jetArea[1],
809 (Double_t)deltaPt, (Double_t)fraction, (Double_t)pthardbin,
fedb5a47 810 (Double_t)jetPt[2],(Double_t)delta,(Double_t)jetArea[2], (Double_t)fraction2, (Double_t)deltaR
d5ec3aef 811 };
812 fhnJets3Branches->Fill(jetEntries3Branches);
813 }
31b9d515 814
d5ec3aef 815 }
31b9d515 816
d5ec3aef 817 PostData(1, fOutputList);
3c6a60f7 818}
819
820void AliAnalysisTaskJetResponseV2::Terminate(const Option_t *)
821{
d5ec3aef 822 // Draw result to the screen
823 // Called once at the end of the query
3c6a60f7 824
d5ec3aef 825 if (!GetOutputData(1))
826 return;
3c6a60f7 827}
828
829Int_t AliAnalysisTaskJetResponseV2::GetNInputTracks()
830{
d5ec3aef 831 // number of tracks in the event, obtained via jet finder
832
833 Int_t nInputTracks = 0;
834
835 TString jbname(fJetBranchName[1]);
836 //needs complete event, use jets without background subtraction
837 for(Int_t i=1; i<=3; ++i){
838 if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0");
839 }
840 // use only HI event
841 if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
9993bed0 842 else if(jbname.Contains("AODMCextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
d5ec3aef 843 if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
9993bed0 844 else if(jbname.Contains("AODMCextra")) jbname.ReplaceAll("AODextra","AOD");
d5ec3aef 845
846 if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
847 TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(jbname.Data()));
f34e09e6 848 if(!tmpAODjets && fAODOut) tmpAODjets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(jbname.Data()));
53f9653b 849 if(!tmpAODjets && fAODExtension) tmpAODjets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(jbname.Data()));
d5ec3aef 850 if(!tmpAODjets){
851 Printf("Jet branch %s not found", jbname.Data());
852 Printf("AliAnalysisTaskJetResponseV2::GetNInputTracks FAILED");
853 return -1;
854 }
31b9d515 855
d5ec3aef 856 for (Int_t iJet=0; iJet<tmpAODjets->GetEntriesFast(); iJet++){
857 AliAODJet *jet = dynamic_cast<AliAODJet*>((*tmpAODjets)[iJet]);
858 if(!jet) continue;
859 TRefArray *trackList = jet->GetRefTracks();
860 Int_t nTracks = trackList->GetEntriesFast();
861 nInputTracks += nTracks;
862 if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks);
863 }
864 if(fDebug) Printf("---> input tracks: %d", nInputTracks);
865
866 return nInputTracks;
3c6a60f7 867}
d314de7b 868
869THnSparse* AliAnalysisTaskJetResponseV2::NewTHnSparseF(const char* name, UInt_t entries, UInt_t opt)
870{
d5ec3aef 871 // generate new THnSparseF, axes are defined in GetDimParams()
872
873 Int_t count = 0;
874 UInt_t tmp = entries;
875 while(tmp!=0){
876 count++;
877 tmp = tmp &~ -tmp; // clear lowest bit
878 }
879
880 TString hnTitle(name);
881 const Int_t dim = count;
882 Int_t nbins[dim];
883 Double_t xmin[dim];
884 Double_t xmax[dim];
885
886 Int_t i=0;
887 Int_t c=0;
888 while(c<dim && i<32){
889 if(entries&(1<<i)){
890 Bool_t highres = opt&(1<<i);
891 TString label("");
892 GetDimParams(i, highres, label, nbins[c], xmin[c], xmax[c]);
893 hnTitle += Form(";%s",label.Data());
894 c++;
895 }
31b9d515 896
d5ec3aef 897 i++;
898 }
899 hnTitle += ";";
31b9d515 900
d5ec3aef 901 return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
d314de7b 902}
903
904void AliAnalysisTaskJetResponseV2::GetDimParams(Int_t iEntry, Bool_t hr, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
905{
d5ec3aef 906 // stores label and binning of axis for THnSparse
d314de7b 907
d5ec3aef 908 const Double_t pi = TMath::Pi();
d314de7b 909
d5ec3aef 910 switch(iEntry){
31b9d515 911
d5ec3aef 912 case 0:
913 label = "V0 centrality (%)";
914 if(hr){
915 nbins = 100;
916 xmin = 0.;
917 xmax = 100.;
918 } else {
919 nbins = 8;
920 xmin = 0.;
921 xmax = 80.;
922 }
923 break;
d314de7b 924
925
d5ec3aef 926 case 1:
927 label = "nb. of input tracks";
928 if(fIsPbPb){
929 if(hr){
930 nbins = 400;
931 xmin = 0.;
932 xmax = 4000.;
31b9d515 933 } else {
d5ec3aef 934 nbins = 40;
935 xmin = 0.;
936 xmax = 4000.;
31b9d515 937 }
d5ec3aef 938 } else {
939 nbins = 40;
940 xmin = 0.;
941 xmax = 400.;
942 }
943 break;
d314de7b 944
945
d5ec3aef 946 case 2:
947 label = "event plane #psi";
948 if(hr){
949 nbins = 30;
950 xmin = 0.;
951 xmax = pi;
952 } else {
953 nbins = 30;
954 xmin = 0.;
955 xmax = pi;
956 }
957 break;
d314de7b 958
959
d5ec3aef 960 case 3:
961 label = "event plane bin";
962 nbins = 3;
963 xmin = -.5;
964 xmax = 2.5;
965 break;
d314de7b 966
967
d5ec3aef 968 case 4:
969 case 5:
970 if(iEntry==4)label = "#Delta#phi(RP-jet) (probe)";
971 if(iEntry==5)label = "#Delta#phi(RP-jet) (rec)";
972 nbins = 48;
973 xmin = -pi;
974 xmax = pi;
975 break;
d314de7b 976
977
d5ec3aef 978 case 6:
979 case 7:
b5fee0f9 980 case 27:
d5ec3aef 981 if(iEntry==6)label = "probe p_{T} (GeV/c)";
b5fee0f9 982 if(iEntry==7 || iEntry==27)label = "rec p_{T} (GeV/c)";
d5ec3aef 983 if(hr){
984 nbins = 300;
985 xmin = -50.;
986 xmax = 250.;
987 } else {
988 nbins = 50;
989 xmin = 0.;
990 xmax = 250.;
991 }
992 break;
d314de7b 993
994
d5ec3aef 995 case 8:
996 case 9:
997 if(iEntry==8)label = "probe #eta";
998 if(iEntry==9)label = "rec #eta";
999 if(hr){
1000 nbins = 56;
1001 xmin = -.7;
1002 xmax = .7;
1003 } else {
1004 nbins = 28;
1005 xmin = -.7;
1006 xmax = .7;
1007 }
1008 break;
d314de7b 1009
1010
d5ec3aef 1011 case 10:
1012 case 11:
1013 if(iEntry==10)label = "probe #phi";
1014 if(iEntry==11)label = "rec #phi";
1015 if(hr){
1016 nbins = 90; // modulo 18 (sectors)
1017 xmin = 0.;
1018 xmax = 2*pi;
1019 } else {
1020 nbins = 25;
1021 xmin = 0.;
1022 xmax = 2*pi;
1023 }
1024 break;
d314de7b 1025
1026
d5ec3aef 1027 case 12:
1028 case 13:
b5fee0f9 1029 case 29:
d5ec3aef 1030 if(iEntry==12)label = "probe area";
b5fee0f9 1031 if(iEntry==13 || iEntry==29)label = "rec area";
d5ec3aef 1032 if(hr){
1033 nbins = 100;
1034 xmin = 0.;
1035 xmax = 1.;
1036 } else {
1037 nbins = 25;
1038 xmin = 0.;
1039 xmax = 1.;
1040 }
1041 break;
d314de7b 1042
d5ec3aef 1043 case 14:
b5fee0f9 1044 case 28:
1045 if(iEntry==14) label = "#delta p_{T}";
1046 if(iEntry==28) label = "#delta";
d5ec3aef 1047 if(hr){
1048 nbins = 241;
1049 xmin = -120.5;
1050 xmax = 120.5;
1051 } else {
1052 nbins = 101;
1053 xmin = -101.;
1054 xmax = 101.;
1055 }
1056 break;
d314de7b 1057
d5ec3aef 1058 case 15:
1059 label = "#Delta#eta";
1060 if(hr){
1061 nbins = 51;
1062 xmin = -1.02;
1063 xmax = 1.02;
1064 } else {
1065 nbins = 51;
1066 xmin = -1.02;
1067 xmax = 1.02;
1068 }
1069 break;
d314de7b 1070
1071
d5ec3aef 1072 case 16:
1073 label = "#Delta#phi";
1074 if(hr){
1075 nbins = 45;
1076 xmin = -pi;
1077 xmax = pi;
1078 } else {
1079 nbins = 45;
1080 xmin = -pi;
1081 xmax = pi;
1082 }
1083 break;
d314de7b 1084
1085
d5ec3aef 1086 case 17:
1087 label = "#DeltaR";
1088 if(hr){
1089 nbins = 50;
1090 xmin = 0.;
1091 xmax = 1.;
1092 } else {
1093 nbins = 50;
1094 xmin = 0.;
1095 xmax = 1.;
1096 }
1097 break;
d314de7b 1098
1099
d5ec3aef 1100 case 18:
1101 label = "#Deltaarea";
1102 if(hr){
1103 nbins = 81;
1104 xmin = -.81;
1105 xmax = .81;
1106 } else {
1107 nbins = 33;
1108 xmin = -.825;
1109 xmax = .825;
1110 }
1111 break;
d314de7b 1112
1113
d5ec3aef 1114 case 19:
b5fee0f9 1115 case 30:
d5ec3aef 1116 label = "fraction";
1117 if(hr){
1118 nbins = 52;
1119 xmin = 0.;
1120 xmax = 1.04;
1121 } else {
1122 nbins = 52;
1123 xmin = 0.;
1124 xmax = 1.04;
1125 }
1126 break;
d314de7b 1127
1128
d5ec3aef 1129 case 20:
1130 label = "distance to closest rec jet";
1131 if(hr){
1132 nbins = 51;
1133 xmin = -0.02;
1134 xmax = 1.;
1135 } else {
1136 nbins = 51;
1137 xmin = -0.02;
1138 xmax = 1.;
1139 }
1140 break;
d314de7b 1141
1142
d5ec3aef 1143 case 21:
1144 label = "p_{T} of closest rec jet";
1145 nbins = 100;
1146 xmin = 0.;
1147 xmax = 200.;
1148 break;
d314de7b 1149
d5ec3aef 1150 case 22:
1151 label = "#rho";
1152 nbins = 125;
1153 xmin = 0.;
1154 xmax = 250.;
1155 break;
d314de7b 1156
d5ec3aef 1157 case 23:
1158 label = "abs. correction of #rho for RP";
1159 nbins = 51;
1160 xmin = -51.;
1161 xmax = 51.;
1162 break;
d314de7b 1163
d5ec3aef 1164 case 24:
1165 label = "local #rho";
1166 nbins = 125;
1167 xmin = 0.;
1168 xmax = 250.;
1169 break;
d314de7b 1170
d5ec3aef 1171 case 25:
1172 label = "local #rho / #rho";
1173 nbins = 500;
1174 xmin = 0.;
1175 xmax = 5.;
1176 break;
31b9d515 1177
d5ec3aef 1178 case 26:
1179 label = "p_{T,hard} bin";
1180 nbins = 10;
1181 xmin = -.5;
1182 xmax = 9.5;
1183 break;
31b9d515 1184
d5ec3aef 1185 }
31b9d515 1186
1187}
1188
1189//____________________________________________________________________________
1190Int_t AliAnalysisTaskJetResponseV2::GetPtHardBin(Double_t ptHard){
1191
d5ec3aef 1192 // returns pt hard bin (for LHC10e14, LHC11a1x, LHC11a2x simulations)
7a88ef30 1193
d5ec3aef 1194 const Int_t nBins = 10;
1195 Double_t binLimits[nBins] = { 5., 11., 21., 36., 57., 84., 117., 156., 200., 249. }; // lower limits
d314de7b 1196
d5ec3aef 1197 Int_t bin = -1;
1198 while(bin<nBins-1 && binLimits[bin+1]<ptHard){
1199 bin++;
1200 }
31b9d515 1201
d5ec3aef 1202 return bin;
31b9d515 1203}
1204
1205//____________________________________________________________________________
1206Double_t AliAnalysisTaskJetResponseV2::GetPt(AliAODJet *j, Int_t mode=0){
1207
d5ec3aef 1208 // returns jet pt, also negative pt after background subtraction if available
7a88ef30 1209
d5ec3aef 1210 Double_t pt = 0.;
31b9d515 1211
d5ec3aef 1212 if(fKeepJets && mode>0){ // background subtracted pt, can be negative
1213 pt = j->GetPtSubtracted(0);
1214 }
1215 else{
1216 pt = j->Pt(); // if negative, pt=0.01
1217 }
d314de7b 1218
d5ec3aef 1219 return pt;
d314de7b 1220}