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