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