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