]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PWGJE/AliAnalysisTaskJetResponseV2.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Merge branch 'master' of https://git.cern.ch/reps/AliRoot into HEAD
[u/mrichter/AliRoot.git] / PWGJE / AliAnalysisTaskJetResponseV2.cxx
... / ...
CommitLineData
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//
17// task compares jets in two or three branches,
18// written for analysis of jet embedding in HI events
19//
20
21
22#include "TChain.h"
23#include "TTree.h"
24#include "TMath.h"
25#include "TH1F.h"
26#include "TH1I.h"
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"
43#include "AliAODJetEventBackground.h"
44#include "AliAnalysisTaskFastEmbedding.h"
45
46#include "AliAODEvent.h"
47#include "AliAODJet.h"
48#include "AliAODHandler.h"
49
50#include "AliAnalysisTaskJetResponseV2.h"
51
52using std::cout;
53using std::endl;
54
55ClassImp(AliAnalysisTaskJetResponseV2)
56
57AliAnalysisTaskJetResponseV2::AliAnalysisTaskJetResponseV2() :
58AliAnalysisTaskSE(),
59 fESD(0x0),
60 fAOD(0x0),
61 fAODOut(0x0),
62 fAODExtension(0x0),
63 fNonStdFile(""),
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)
112{
113 // default Constructor
114
115 fJetBranchName[0] = "";
116 fJetBranchName[1] = "";
117 fJetBranchName[2] = "";
118
119 fListJets[0] = new TList;
120 fListJets[1] = new TList;
121 fListJets[2] = new TList;
122}
123
124AliAnalysisTaskJetResponseV2::AliAnalysisTaskJetResponseV2(const char *name) :
125 AliAnalysisTaskSE(name),
126 fESD(0x0),
127 fAOD(0x0),
128 fAODOut(0x0),
129 fAODExtension(0x0),
130 fNonStdFile(""),
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)
179{
180 // Constructor
181
182 fJetBranchName[0] = "";
183 fJetBranchName[1] = "";
184 fJetBranchName[2] = "";
185
186 fListJets[0] = new TList;
187 fListJets[1] = new TList;
188 fListJets[2] = new TList;
189
190 DefineOutput(1, TList::Class());
191}
192
193AliAnalysisTaskJetResponseV2::~AliAnalysisTaskJetResponseV2()
194{
195 delete fListJets[0];
196 delete fListJets[1];
197 delete fListJets[2];
198}
199
200void AliAnalysisTaskJetResponseV2::SetBranchNames(const TString &branch1, const TString &branch2, const TString &branch3)
201{
202 fJetBranchName[0] = branch1;
203 fJetBranchName[1] = branch2;
204 fJetBranchName[2] = branch3;
205
206 if(strlen(fJetBranchName[2].Data()) ) {
207 fkNbranches = 3;
208 fbJets3Branches = kTRUE;
209 }
210}
211
212void AliAnalysisTaskJetResponseV2::Init()
213{
214
215 // check for jet branches
216 if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
217 AliError("Jet branch name not set.");
218 }
219
220}
221
222void AliAnalysisTaskJetResponseV2::UserCreateOutputObjects()
223{
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 }
271
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 }
286
287
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
327 // cent : nInpTrks : jetPt(3x) : deltaPt : delta : jetArea(3x) : fraction(2x) : deltaR(1x) : pT hard bin
328 if(fbJets3Branches){
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;
330 opt = 1<<6 | 1<<7 | 1<<27 | 1<<14 | 1<<28;
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);
383}
384
385void AliAnalysisTaskJetResponseV2::UserExec(Option_t *)
386{
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");
397
398 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
399 // assume that the AOD is in the general output...
400 fAODOut = AODEvent();
401 } else {
402 fAOD = dynamic_cast<AliAODEvent*>(AODEvent());
403 }
404 if (!fAOD) {
405 AliError("AOD not available");
406 return;
407 }
408
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
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 }
473
474
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()));
494 if(!externalBackground && fAODOut) externalBackground = (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch.Data()));
495 if(!externalBackground && fAODExtension) externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
496 //if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
497 }
498 Float_t rho = 0;
499 if(externalBackground) rho = externalBackground->GetBackground(0);
500
501
502 // fetch jets
503 TClonesArray *aodJets[3];
504 aodJets[0] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[0].Data())); // in general: embedded jet
505 if(!aodJets[0] && fAODOut) aodJets[0] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[0].Data())); // in general: embedded jet
506 if(!aodJets[0] && fAODExtension) aodJets[0] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[0].Data()));
507 aodJets[1] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[1].Data())); // in general: embedded jet + UE version1
508 if(!aodJets[1] && fAODOut) aodJets[1] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[1].Data())); // in general: embedded jet
509 if(!aodJets[1] && fAODExtension) aodJets[1] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[1].Data()));
510 if( strlen(fJetBranchName[2].Data()) ) {
511 aodJets[2] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[2].Data())); // in general: embedded jet + UE version2
512 if(!aodJets[2] && fAODOut) aodJets[2] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[2].Data())); // in general: embedded jet
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());
516 }
517
518 if(fDebug>10) printf("fkNbranches %d\n",fkNbranches);
519 for (Int_t iJetType = 0; iJetType < fkNbranches; iJetType++) {
520 fListJets[iJetType]->Clear();
521 if (!aodJets[iJetType]) {
522 if(fDebug) Printf("%s: no jet branch\n",fJetBranchName[iJetType].Data());
523 continue;
524 }
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 }
533
534
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());
547 if(aMatchIndexv2.GetSize()<fListJets[0]->GetEntries()) aMatchIndexv2.Set(fListJets[0]->GetEntries());
548 if(aPtFractionv2.GetSize()<fListJets[0]->GetEntries()) aPtFractionv2.Set(fListJets[0]->GetEntries());
549 if( strlen(fJetBranchName[2].Data()) ) {
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);
553 }
554
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;
567
568 for(Int_t ig=0; ig<fListJets[0]->GetEntries(); ++ig){
569 ir = aMatchIndex[ig];
570 if(aMatchIndexv2.GetSize()>ig) ir2 = aMatchIndexv2[ig];
571 else ir2 =0;
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);
596 }
597 }
598 fraction = aPtFraction[ig];
599 if(aPtFractionv2.GetSize()>ig) fraction2 = aPtFractionv2[ig];
600 else fraction2 = 0.;
601
602 // jet statistics
603 fHistJetSelection->Fill(1); // all probe jets
604 if(jet[0]) fh2JetSelection->Fill(1.,jetPt[0]); // all probe jets
605
606 if(ir<0){
607 fHistJetSelection->Fill(2);
608 if(jet[0]) fh2JetSelection->Fill(2.,jetPt[0]);
609
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);
617 }
618
619 continue;
620 }
621
622 if(!jet[0] || !jet[1]){
623 fHistJetSelection->Fill(3);
624 if(jet[0]) fh2JetSelection->Fill(3.,jetPt[0]);
625 continue;
626 }
627
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();
638 }
639 }
640
641
642
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.;
645
646
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 };
657
658 Float_t rpCorr = 0.;
659
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];
671
672
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 }
682
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 }
692
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 }
701
702 if(jetAccepted){
703 // jet acceptance + minimum pT check
704 if(jetEta[0]>fJetEtaMax || jetEta[0]<fJetEtaMin ||
705 jetEta[1]>fJetEtaMax || jetEta[1]<fJetEtaMin){
706
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;
715 }
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;
723 }
724 }
725
726 if(jetAccepted){
727 if(jet[1]->Trigger()&fJetTriggerExcludeMask){
728 fHistJetSelection->Fill(7);
729 fh2JetSelection->Fill(7.,jetPt[0]);
730 jetAccepted = kFALSE;
731 }
732
733 }
734
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);
745 }
746 continue;
747 }
748
749 // all accepted jets
750 fHistJetSelection->Fill(0);
751 fh2JetSelection->Fill(0.,jetPt[0]);
752
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 }
763
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 }
772
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 }
782
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 }
792
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){
805 Double_t jetEntries3Branches[14] = {
806 (Double_t)centValue, (Double_t)nInputTracks,
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,
810 (Double_t)jetPt[2],(Double_t)delta,(Double_t)jetArea[2], (Double_t)fraction2, (Double_t)deltaR
811 };
812 fhnJets3Branches->Fill(jetEntries3Branches);
813 }
814
815 }
816
817 PostData(1, fOutputList);
818}
819
820void AliAnalysisTaskJetResponseV2::Terminate(const Option_t *)
821{
822 // Draw result to the screen
823 // Called once at the end of the query
824
825 if (!GetOutputData(1))
826 return;
827}
828
829Int_t AliAnalysisTaskJetResponseV2::GetNInputTracks()
830{
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");
842 else if(jbname.Contains("AODMCextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
843 if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
844 else if(jbname.Contains("AODMCextra")) jbname.ReplaceAll("AODextra","AOD");
845
846 if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
847 TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(jbname.Data()));
848 if(!tmpAODjets && fAODOut) tmpAODjets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(jbname.Data()));
849 if(!tmpAODjets && fAODExtension) tmpAODjets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(jbname.Data()));
850 if(!tmpAODjets){
851 Printf("Jet branch %s not found", jbname.Data());
852 Printf("AliAnalysisTaskJetResponseV2::GetNInputTracks FAILED");
853 return -1;
854 }
855
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;
867}
868
869THnSparse* AliAnalysisTaskJetResponseV2::NewTHnSparseF(const char* name, UInt_t entries, UInt_t opt)
870{
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 }
896
897 i++;
898 }
899 hnTitle += ";";
900
901 return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
902}
903
904void AliAnalysisTaskJetResponseV2::GetDimParams(Int_t iEntry, Bool_t hr, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
905{
906 // stores label and binning of axis for THnSparse
907
908 const Double_t pi = TMath::Pi();
909
910 switch(iEntry){
911
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;
924
925
926 case 1:
927 label = "nb. of input tracks";
928 if(fIsPbPb){
929 if(hr){
930 nbins = 400;
931 xmin = 0.;
932 xmax = 4000.;
933 } else {
934 nbins = 40;
935 xmin = 0.;
936 xmax = 4000.;
937 }
938 } else {
939 nbins = 40;
940 xmin = 0.;
941 xmax = 400.;
942 }
943 break;
944
945
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;
958
959
960 case 3:
961 label = "event plane bin";
962 nbins = 3;
963 xmin = -.5;
964 xmax = 2.5;
965 break;
966
967
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;
976
977
978 case 6:
979 case 7:
980 case 27:
981 if(iEntry==6)label = "probe p_{T} (GeV/c)";
982 if(iEntry==7 || iEntry==27)label = "rec p_{T} (GeV/c)";
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;
993
994
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;
1009
1010
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;
1025
1026
1027 case 12:
1028 case 13:
1029 case 29:
1030 if(iEntry==12)label = "probe area";
1031 if(iEntry==13 || iEntry==29)label = "rec area";
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;
1042
1043 case 14:
1044 case 28:
1045 if(iEntry==14) label = "#delta p_{T}";
1046 if(iEntry==28) label = "#delta";
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;
1057
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;
1070
1071
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;
1084
1085
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;
1098
1099
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;
1112
1113
1114 case 19:
1115 case 30:
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;
1127
1128
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;
1141
1142
1143 case 21:
1144 label = "p_{T} of closest rec jet";
1145 nbins = 100;
1146 xmin = 0.;
1147 xmax = 200.;
1148 break;
1149
1150 case 22:
1151 label = "#rho";
1152 nbins = 125;
1153 xmin = 0.;
1154 xmax = 250.;
1155 break;
1156
1157 case 23:
1158 label = "abs. correction of #rho for RP";
1159 nbins = 51;
1160 xmin = -51.;
1161 xmax = 51.;
1162 break;
1163
1164 case 24:
1165 label = "local #rho";
1166 nbins = 125;
1167 xmin = 0.;
1168 xmax = 250.;
1169 break;
1170
1171 case 25:
1172 label = "local #rho / #rho";
1173 nbins = 500;
1174 xmin = 0.;
1175 xmax = 5.;
1176 break;
1177
1178 case 26:
1179 label = "p_{T,hard} bin";
1180 nbins = 10;
1181 xmin = -.5;
1182 xmax = 9.5;
1183 break;
1184
1185 }
1186
1187}
1188
1189//____________________________________________________________________________
1190Int_t AliAnalysisTaskJetResponseV2::GetPtHardBin(Double_t ptHard){
1191
1192 // returns pt hard bin (for LHC10e14, LHC11a1x, LHC11a2x simulations)
1193
1194 const Int_t nBins = 10;
1195 Double_t binLimits[nBins] = { 5., 11., 21., 36., 57., 84., 117., 156., 200., 249. }; // lower limits
1196
1197 Int_t bin = -1;
1198 while(bin<nBins-1 && binLimits[bin+1]<ptHard){
1199 bin++;
1200 }
1201
1202 return bin;
1203}
1204
1205//____________________________________________________________________________
1206Double_t AliAnalysisTaskJetResponseV2::GetPt(AliAODJet *j, Int_t mode=0){
1207
1208 // returns jet pt, also negative pt after background subtraction if available
1209
1210 Double_t pt = 0.;
1211
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 }
1218
1219 return pt;
1220}