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