Possibility to set the first event.
[u/mrichter/AliRoot.git] / PWGJE / AliAnalysisTaskJetCore.cxx
CommitLineData
a9e585a7 1
568f8fa2 2// ******************************************
3// This task computes several jet observables like
4// the fraction of energy in inner and outer coronnas,
da93bb11 5// jet-track correlations,triggered jet shapes and
568f8fa2 6// correlation strength distribution of particles inside jets.
7// Author: lcunquei@cern.ch
8// *******************************************
9
10
11/**************************************************************************
12 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
13 * *
14 * Author: The ALICE Off-line Project. *
15 * Contributors are mentioned in the code where appropriate. *
16 * *
17 * Permission to use, copy, modify and distribute this software and its *
18 * documentation strictly for non-commercial purposes is hereby granted *
19 * without fee, provided that the above copyright notice appears in all *
20 * copies and that both the copyright notice and this permission notice *
21 * appear in the supporting documentation. The authors make no claims *
22 * about the suitability of this software for any purpose. It is *
23 * provided "as is" without express or implied warranty. *
24 **************************************************************************/
25
26
75bf77e3 27#include "TChain.h"
28#include "TTree.h"
29#include "TMath.h"
30#include "TH1F.h"
31#include "TH2F.h"
32#include "TH3F.h"
33#include "THnSparse.h"
34#include "TCanvas.h"
35
36#include "AliLog.h"
37
38#include "AliAnalysisTask.h"
39#include "AliAnalysisManager.h"
40
41#include "AliVEvent.h"
42#include "AliESDEvent.h"
43#include "AliESDInputHandler.h"
44#include "AliCentrality.h"
45#include "AliAnalysisHelperJetTasks.h"
46#include "AliInputEventHandler.h"
47#include "AliAODJetEventBackground.h"
48#include "AliAnalysisTaskFastEmbedding.h"
75bf77e3 49#include "AliAODEvent.h"
ea693273 50#include "AliAODHandler.h"
75bf77e3 51#include "AliAODJet.h"
52
53#include "AliAnalysisTaskJetCore.h"
54
55ClassImp(AliAnalysisTaskJetCore)
56
57AliAnalysisTaskJetCore::AliAnalysisTaskJetCore() :
58AliAnalysisTaskSE(),
59fESD(0x0),
5bd732d4 60fAODIn(0x0),
61fAODOut(0x0),
ea693273 62fAODExtension(0x0),
75bf77e3 63fBackgroundBranch(""),
ea693273 64fNonStdFile(""),
75bf77e3 65fIsPbPb(kTRUE),
66fOfflineTrgMask(AliVEvent::kAny),
67fMinContribVtx(1),
a9e585a7 68fVtxZMin(-10.),
69fVtxZMax(10.),
75bf77e3 70fEvtClassMin(0),
71fEvtClassMax(4),
8b47ec90 72fFilterMask(0),
ea693273 73fRadioFrac(0.2),
74fMinDist(0.1),
75bf77e3 75fCentMin(0.),
76fCentMax(100.),
77fNInputTracksMin(0),
78fNInputTracksMax(-1),
ea693273 79fAngStructCloseTracks(0),
da93bb11 80fCheckMethods(0),
529e2916 81fDoEventMixing(0),
75bf77e3 82fJetEtaMin(-.5),
83fJetEtaMax(.5),
529e2916 84fNevents(0x0),
85fTindex(0x0),
86fTrigBufferIndex(0x0),
75bf77e3 87fJetPtMin(20.),
88fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
89fJetPtFractionMin(0.5),
90fNMatchJets(4),
91fMatchMaxDist(0.8),
92fKeepJets(kFALSE),
75bf77e3 93fkNbranches(2),
94fkEvtClasses(12),
95fOutputList(0x0),
96fbEvent(kTRUE),
97fHistEvtSelection(0x0),
da93bb11 98fhnDeltaR(0x0),
529e2916 99fhnMixedEvents(0x0),
75bf77e3 100fh2JetCoreMethod1C10(0x0),
101fh2JetCoreMethod2C10(0x0),
ea693273 102fh2JetCoreMethod1C20(0x0),
103fh2JetCoreMethod2C20(0x0),
75bf77e3 104fh2JetCoreMethod1C30(0x0),
105fh2JetCoreMethod2C30(0x0),
75bf77e3 106fh2JetCoreMethod1C60(0x0),
107fh2JetCoreMethod2C60(0x0),
ea693273 108fh2AngStructpt1C10(0x0),
109fh2AngStructpt2C10(0x0),
110fh2AngStructpt3C10(0x0),
111fh2AngStructpt4C10(0x0),
112fh2AngStructpt1C20(0x0),
113fh2AngStructpt2C20(0x0),
114fh2AngStructpt3C20(0x0),
115fh2AngStructpt4C20(0x0),
116fh2AngStructpt1C30(0x0),
117fh2AngStructpt2C30(0x0),
118fh2AngStructpt3C30(0x0),
119fh2AngStructpt4C30(0x0),
120fh2AngStructpt1C60(0x0),
121fh2AngStructpt2C60(0x0),
122fh2AngStructpt3C60(0x0),
da93bb11 123fh2AngStructpt4C60(0x0),
5bd732d4 124fh2JetsumHT3R2(0x0),
125fh2JetsumHT3R2gap(0x0),
126fh2JetsumHT3R4(0x0),
127fh2JetsumHT3R4gap(0x0),
128fh2JetsumHT3R6(0x0),
129fh2JetsumHT3R6gap(0x0),
130fh2JetsumHT3R8(0x0),
131fh2JetsumHT3R8gap(0x0),
132fh2JetsumHT3R10(0x0),
133fh2JetsumHT3R10gap(0x0),
da93bb11 134fh3spectriggered(0x0),
135fh3specbiased(0x0),
20dcc500 136fh3spectot(0x0),
137fh3spectotb(0x0)
da93bb11 138
75bf77e3 139{
140 // default Constructor
141
529e2916 142
143 // Trigger buffer.
144 for(Int_t i=0; i<10; i++) {
145 for(Int_t j=0; j<7; j++) {
146 fTrigBuffer[i][j]=0;
147 }
148 }
149
150
151
152
153
75bf77e3 154 fJetBranchName[0] = "";
155 fJetBranchName[1] = "";
156
157 fListJets[0] = new TList;
158 fListJets[1] = new TList;
159}
160
161AliAnalysisTaskJetCore::AliAnalysisTaskJetCore(const char *name) :
162AliAnalysisTaskSE(name),
163fESD(0x0),
5bd732d4 164fAODIn(0x0),
165fAODOut(0x0),
ea693273 166fAODExtension(0x0),
75bf77e3 167fBackgroundBranch(""),
ea693273 168fNonStdFile(""),
75bf77e3 169fIsPbPb(kTRUE),
170fOfflineTrgMask(AliVEvent::kAny),
171fMinContribVtx(1),
a9e585a7 172fVtxZMin(-10.),
173fVtxZMax(10.),
75bf77e3 174fEvtClassMin(0),
175fEvtClassMax(4),
8b47ec90 176fFilterMask(0),
ea693273 177fRadioFrac(0.2),
178fMinDist(0.1),
75bf77e3 179fCentMin(0.),
180fCentMax(100.),
181fNInputTracksMin(0),
182fNInputTracksMax(-1),
ea693273 183fAngStructCloseTracks(0),
da93bb11 184fCheckMethods(0),
529e2916 185fDoEventMixing(0),
75bf77e3 186fJetEtaMin(-.5),
187fJetEtaMax(.5),
529e2916 188fNevents(0x0),
189fTindex(0x0),
190fTrigBufferIndex(0x0),
75bf77e3 191fJetPtMin(20.),
192fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
193fJetPtFractionMin(0.5),
194fNMatchJets(4),
195fMatchMaxDist(0.8),
196fKeepJets(kFALSE),
75bf77e3 197fkNbranches(2),
198fkEvtClasses(12),
199fOutputList(0x0),
200fbEvent(kTRUE),
201fHistEvtSelection(0x0),
da93bb11 202fhnDeltaR(0x0),
529e2916 203fhnMixedEvents(0x0),
75bf77e3 204fh2JetCoreMethod1C10(0x0),
205fh2JetCoreMethod2C10(0x0),
ea693273 206fh2JetCoreMethod1C20(0x0),
207fh2JetCoreMethod2C20(0x0),
75bf77e3 208fh2JetCoreMethod1C30(0x0),
209fh2JetCoreMethod2C30(0x0),
75bf77e3 210fh2JetCoreMethod1C60(0x0),
211fh2JetCoreMethod2C60(0x0),
ea693273 212fh2AngStructpt1C10(0x0),
213fh2AngStructpt2C10(0x0),
214fh2AngStructpt3C10(0x0),
215fh2AngStructpt4C10(0x0),
216fh2AngStructpt1C20(0x0),
217fh2AngStructpt2C20(0x0),
218fh2AngStructpt3C20(0x0),
219fh2AngStructpt4C20(0x0),
220fh2AngStructpt1C30(0x0),
221fh2AngStructpt2C30(0x0),
222fh2AngStructpt3C30(0x0),
223fh2AngStructpt4C30(0x0),
224fh2AngStructpt1C60(0x0),
225fh2AngStructpt2C60(0x0),
226fh2AngStructpt3C60(0x0),
da93bb11 227fh2AngStructpt4C60(0x0),
5bd732d4 228fh2JetsumHT3R2(0x0),
229fh2JetsumHT3R2gap(0x0),
230fh2JetsumHT3R4(0x0),
231fh2JetsumHT3R4gap(0x0),
232fh2JetsumHT3R6(0x0),
233fh2JetsumHT3R6gap(0x0),
234fh2JetsumHT3R8(0x0),
235fh2JetsumHT3R8gap(0x0),
236fh2JetsumHT3R10(0x0),
237fh2JetsumHT3R10gap(0x0),
da93bb11 238fh3spectriggered(0x0),
239fh3specbiased(0x0),
20dcc500 240fh3spectot(0x0),
241fh3spectotb(0x0)
75bf77e3 242 {
243 // Constructor
244
529e2916 245
246 for(Int_t i=0; i<10; i++) {
247 for(Int_t j=0; j<7; j++) {
248 fTrigBuffer[i][j]=0;
249 }
250 }
251
252
253
75bf77e3 254 fJetBranchName[0] = "";
255 fJetBranchName[1] = "";
256
257 fListJets[0] = new TList;
258 fListJets[1] = new TList;
259
260 DefineOutput(1, TList::Class());
261}
262
263AliAnalysisTaskJetCore::~AliAnalysisTaskJetCore()
264{
265 delete fListJets[0];
266 delete fListJets[1];
267}
268
269void AliAnalysisTaskJetCore::SetBranchNames(const TString &branch1, const TString &branch2)
270{
271 fJetBranchName[0] = branch1;
272 fJetBranchName[1] = branch2;
273}
274
275void AliAnalysisTaskJetCore::Init()
276{
277
278 // check for jet branches
279 if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
280 AliError("Jet branch name not set.");
281 }
282
283}
284
285void AliAnalysisTaskJetCore::UserCreateOutputObjects()
286{
287 // Create histograms
288 // Called once
289 OpenFile(1);
290 if(!fOutputList) fOutputList = new TList;
291 fOutputList->SetOwner(kTRUE);
292
293 Bool_t oldStatus = TH1::AddDirectoryStatus();
294 TH1::AddDirectory(kFALSE);
295
296
297 fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
298 fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
299 fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
300 fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
301 fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
302 fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
303 fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
304
529e2916 305 UInt_t entries = 0; // bit coded, see GetDimParams() below
306 entries = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 |1<<7;
da93bb11 307 fhnDeltaR = NewTHnSparseF("fhnDeltaR", entries);
c42c6838 308
20dcc500 309 //change binning in pTtrack
310 Double_t *xPt3 = new Double_t[10];
c42c6838 311 xPt3[0] = 0.;
20dcc500 312 for(int i = 1; i<=9;i++){
c42c6838 313 if(xPt3[i-1]<1)xPt3[i] = xPt3[i-1] + 0.2; // 1 - 5
20dcc500 314 else if(xPt3[i-1]<10)xPt3[i] = xPt3[i-1] + 3; // 5 - 12
315 else xPt3[i] = xPt3[i-1] + 150.; // 18
c42c6838 316 }
317 fhnDeltaR->SetBinEdges(2,xPt3);
318 delete [] xPt3;
319
20dcc500 320 //change binning in HTI
321 Double_t *xPt4 = new Double_t[14];
322 xPt4[0] = 0.;
323 for(int i = 1; i<=13;i++){
324 if(xPt4[i-1]<10)xPt4[i] = xPt4[i-1] + 1; // 1 - 10
325 else if(xPt4[i-1]<20)xPt4[i] = xPt4[i-1] + 5; // 10 - 12
326 else xPt4[i] = xPt4[i-1] + 30.; // 13
327 }
328 fhnDeltaR->SetBinEdges(6,xPt4);
329 delete [] xPt4;
330
331
332
c42c6838 333
334
335
529e2916 336
c42c6838 337 if(fDoEventMixing){
529e2916 338 UInt_t cifras = 0; // bit coded, see GetDimParams() below
339 cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 |1<<7;
c42c6838 340 fhnMixedEvents = NewTHnSparseF("fhnMixedEvents", cifras);}
ea693273 341
da93bb11 342 if(fCheckMethods){
529e2916 343
75bf77e3 344
345 fh2JetCoreMethod1C10 = new TH2F("JetCoreMethod1C10","",150, 0., 150.,100, 0., 1.5);
346 fh2JetCoreMethod2C10 = new TH2F("JetCoreMethod2C10","",150, 0., 150.,100, 0., 1.5);
ea693273 347 fh2JetCoreMethod1C20 = new TH2F("JetCoreMethod1C20","",150, 0., 150.,100, 0., 1.5);
348 fh2JetCoreMethod2C20 = new TH2F("JetCoreMethod2C20","",150, 0., 150.,100, 0., 1.5);
75bf77e3 349 fh2JetCoreMethod1C30 = new TH2F("JetCoreMethod1C30","",150, 0., 150.,100, 0., 1.5);
350 fh2JetCoreMethod2C30 = new TH2F("JetCoreMethod2C30","",150, 0., 150.,100, 0., 1.5);
75bf77e3 351 fh2JetCoreMethod1C60 = new TH2F("JetCoreMethod1C60","",150, 0., 150.,100, 0., 1.5);
da93bb11 352 fh2JetCoreMethod2C60 = new TH2F("JetCoreMethod2C60","",150, 0., 150.,100, 0., 1.5);}
353
ea693273 354
da93bb11 355 if(fAngStructCloseTracks>0){
a9e585a7 356 fh2AngStructpt1C10 = new TH2F("Ang struct pt1 C10","",15,0.,1.5,150,0.,10.);
357 fh2AngStructpt2C10 = new TH2F("Ang struct pt2 C10","",15,0.,1.5,150,0.,10.);
358 fh2AngStructpt3C10 = new TH2F("Ang struct pt3 C10","",15,0.,1.5,150,0.,10.);
359 fh2AngStructpt4C10 = new TH2F("Ang struct pt4 C10","",15,0.,1.5,150,0.,10.);
360 fh2AngStructpt1C20 = new TH2F("Ang struct pt1 C20","",15,0.,1.5,150,0.,10.);
361 fh2AngStructpt2C20 = new TH2F("Ang struct pt2 C20","",15,0.,1.5,150,0.,10.);
362 fh2AngStructpt3C20 = new TH2F("Ang struct pt3 C20","",15,0.,1.5,150,0.,10.);
363 fh2AngStructpt4C20 = new TH2F("Ang struct pt4 C20","",15,0.,1.5,150,0.,10.);
364 fh2AngStructpt1C30 = new TH2F("Ang struct pt1 C30","",15,0.,1.5,150,0.,10.);
365 fh2AngStructpt2C30 = new TH2F("Ang struct pt2 C30","",15,0.,1.5,150,0.,10.);
366 fh2AngStructpt3C30 = new TH2F("Ang struct pt3 C30","",15,0.,1.5,150,0.,10.);
367 fh2AngStructpt4C30 = new TH2F("Ang struct pt4 C30","",15,0.,1.5,150,0.,10.);
368 fh2AngStructpt1C60 = new TH2F("Ang struct pt1 C60","",15,0.,1.5,150,0.,10.);
369 fh2AngStructpt2C60 = new TH2F("Ang struct pt2 C60","",15,0.,1.5,150,0.,10.);
370 fh2AngStructpt3C60 = new TH2F("Ang struct pt3 C60","",15,0.,1.5,150,0.,10.);
371 fh2AngStructpt4C60 = new TH2F("Ang struct pt4 C60","",15,0.,1.5,150,0.,10.); }
20dcc500 372
5bd732d4 373 fh2JetsumHT3R2 = new TH2F("Pt sum R02 HT3","",20,0.,200.,100,0.,10.);
374 fh2JetsumHT3R2gap = new TH2F("Pt sum R02 HT3 gap","",20,0.,200.,100,0.,10.);
375 fh2JetsumHT3R4 = new TH2F("Pt sum R04 HT3","",20,0.,200.,100,0.,10.);
376 fh2JetsumHT3R4gap = new TH2F("Pt sum R04 HT3 gap","",20,0.,200.,100,0.,10.);
377 fh2JetsumHT3R6 = new TH2F("Pt sum R06 HT3","",20,0.,200.,100,0.,10.);
378 fh2JetsumHT3R6gap = new TH2F("Pt sum R06 HT3 gap","",20,0.,200.,100,0.,10.);
379 fh2JetsumHT3R8 = new TH2F("Pt sum R08 HT3","",20,0.,200.,100,0.,10.);
380 fh2JetsumHT3R8gap = new TH2F("Pt sum R08 HT3 gap","",20,0.,200.,100,0.,10.);
381 fh2JetsumHT3R10 = new TH2F("Pt sum R10 HT3","",20,0.,200.,100,0.,10.);
382 fh2JetsumHT3R10gap = new TH2F("Pt sum R10 HT3 gap","",20,0.,200.,100,0.,10.);
da93bb11 383 fh3spectriggered = new TH3F("Triggered spectrum","",10,0,100,50,0.,200,50,0.,50.);
384 fh3specbiased = new TH3F("Biased spectrum","",10,0,100,50,0.,200.,50,0.,50.);
20dcc500 385 fh3spectot = new TH3F("Total spectrum 0-10","",50,0.,200.,50,0.,50.,50,0.,50.);
386 fh3spectotb = new TH3F("Total spectrum 0-20","",50,0.,200.,50,0.,50.,50,0.,50.);
75bf77e3 387 fOutputList->Add(fHistEvtSelection);
a9e585a7 388
da93bb11 389 fOutputList->Add(fhnDeltaR);
390
529e2916 391 fOutputList->Add(fhnMixedEvents);
da93bb11 392
393
394
395 if(fCheckMethods){
529e2916 396
75bf77e3 397 fOutputList->Add(fh2JetCoreMethod1C10);
398 fOutputList->Add(fh2JetCoreMethod2C10);
ea693273 399 fOutputList->Add(fh2JetCoreMethod1C20);
400 fOutputList->Add(fh2JetCoreMethod2C20);
75bf77e3 401 fOutputList->Add(fh2JetCoreMethod1C30);
402 fOutputList->Add(fh2JetCoreMethod2C30);
75bf77e3 403 fOutputList->Add(fh2JetCoreMethod1C60);
da93bb11 404 fOutputList->Add(fh2JetCoreMethod2C60);}
a9e585a7 405
ea693273 406
75bf77e3 407
a9e585a7 408
409
410 if(fAngStructCloseTracks>0){
ea693273 411 fOutputList->Add(fh2AngStructpt1C10);
412 fOutputList->Add(fh2AngStructpt2C10);
413 fOutputList->Add(fh2AngStructpt3C10);
414 fOutputList->Add(fh2AngStructpt4C10);
415 fOutputList->Add(fh2AngStructpt1C20);
416 fOutputList->Add(fh2AngStructpt2C20);
417 fOutputList->Add(fh2AngStructpt3C20);
418 fOutputList->Add(fh2AngStructpt4C20);
419 fOutputList->Add(fh2AngStructpt1C30);
420 fOutputList->Add(fh2AngStructpt2C30);
421 fOutputList->Add(fh2AngStructpt3C30);
422 fOutputList->Add(fh2AngStructpt4C30);
423 fOutputList->Add(fh2AngStructpt1C60);
424 fOutputList->Add(fh2AngStructpt2C60);
425 fOutputList->Add(fh2AngStructpt3C60);
a9e585a7 426 fOutputList->Add(fh2AngStructpt4C60);}
20dcc500 427
428
5bd732d4 429 fOutputList->Add(fh2JetsumHT3R2);
430 fOutputList->Add(fh2JetsumHT3R2gap);
20dcc500 431
20dcc500 432
5bd732d4 433 fOutputList->Add(fh2JetsumHT3R4);
434 fOutputList->Add(fh2JetsumHT3R4gap);
20dcc500 435
20dcc500 436
5bd732d4 437 fOutputList->Add(fh2JetsumHT3R6);
438 fOutputList->Add(fh2JetsumHT3R6gap);
439
440
441 fOutputList->Add(fh2JetsumHT3R8);
442 fOutputList->Add(fh2JetsumHT3R8gap);
443
444
445 fOutputList->Add(fh2JetsumHT3R10);
446 fOutputList->Add(fh2JetsumHT3R10gap);
447
20dcc500 448
da93bb11 449 fOutputList->Add(fh3spectriggered);
450 fOutputList->Add(fh3specbiased);
529e2916 451 fOutputList->Add(fh3spectot);
20dcc500 452 fOutputList->Add(fh3spectotb);
75bf77e3 453 // =========== Switch on Sumw2 for all histos ===========
454 for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
455 TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
456 if (h1){
457 h1->Sumw2();
458 continue;
459 }
ea693273 460 THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
461 if (hn){
462 hn->Sumw2();
463 }
75bf77e3 464 }
465 TH1::AddDirectory(oldStatus);
466
467 PostData(1, fOutputList);
468}
469
470void AliAnalysisTaskJetCore::UserExec(Option_t *)
471{
472
473
474 if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
475 AliError("Jet branch name not set.");
476 return;
477 }
478
479 fESD=dynamic_cast<AliESDEvent*>(InputEvent());
480 if (!fESD) {
481 AliError("ESD not available");
5bd732d4 482 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
483 }
484 fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
485
486 static AliAODEvent* aod = 0;
487 // take all other information from the aod we take the tracks from
488 if(!aod){
489 if(!fESD)aod = fAODIn;
490 else aod = fAODOut;}
491
492
ea693273 493
494 if(fNonStdFile.Length()!=0){
495 // case that we have an AOD extension we need can fetch the jets from the extended output
496 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
497 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
498 if(!fAODExtension){
499 if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
500 }}
501
502
503
504
75bf77e3 505
506 // -- event selection --
507 fHistEvtSelection->Fill(1); // number of events before event selection
508
509 // physics selection
510 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
511 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
da93bb11 512 cout<<inputHandler->IsEventSelected()<<" "<<fOfflineTrgMask<<endl;
75bf77e3 513 if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
514 if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
515 fHistEvtSelection->Fill(2);
516 PostData(1, fOutputList);
517 return;
518 }
519
520 // vertex selection
5bd732d4 521 if(!aod){
52721bee 522 if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
523 fHistEvtSelection->Fill(3);
524 PostData(1, fOutputList);
525 }
5bd732d4 526 AliAODVertex* primVtx = aod->GetPrimaryVertex();
52721bee 527
528 if(!primVtx){
529 if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
530 fHistEvtSelection->Fill(3);
d27719bf 531 PostData(1, fOutputList);
532 return;
52721bee 533 }
534
75bf77e3 535 Int_t nTracksPrim = primVtx->GetNContributors();
536 if ((nTracksPrim < fMinContribVtx) ||
537 (primVtx->GetZ() < fVtxZMin) ||
538 (primVtx->GetZ() > fVtxZMax) ){
539 if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
540 fHistEvtSelection->Fill(3);
541 PostData(1, fOutputList);
542 return;
543 }
544
545 // event class selection (from jet helper task)
546 Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
547 if(fDebug) Printf("Event class %d", eventClass);
548 if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
549 fHistEvtSelection->Fill(4);
550 PostData(1, fOutputList);
551 return;
552 }
553
554 // centrality selection
555 AliCentrality *cent = 0x0;
da93bb11 556 Double_t centValue = 0.;
20dcc500 557 if(fESD) {cent = fESD->GetCentrality();
558 if(cent) centValue = cent->GetCentralityPercentile("V0M");}
5bd732d4 559 else centValue=aod->GetHeader()->GetCentrality();
20dcc500 560
75bf77e3 561 if(fDebug) printf("centrality: %f\n", centValue);
562 if (centValue < fCentMin || centValue > fCentMax){
563 fHistEvtSelection->Fill(4);
564 PostData(1, fOutputList);
565 return;
566 }
567
568
568f8fa2 569 fHistEvtSelection->Fill(0);
570 // accepted events
75bf77e3 571 // -- end event selection --
ea693273 572
75bf77e3 573 // get background
574 AliAODJetEventBackground* externalBackground = 0;
5bd732d4 575 if(fAODOut&&!externalBackground&&fBackgroundBranch.Length()){
576 externalBackground = (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch.Data()));
75bf77e3 577 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
578 }
ea693273 579 if(fAODExtension&&!externalBackground&&fBackgroundBranch.Length()){
580 externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
581 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
582 }
5bd732d4 583
584 if(fAODIn&&!externalBackground&&fBackgroundBranch.Length()){
585 externalBackground = (AliAODJetEventBackground*)(fAODIn->FindListObject(fBackgroundBranch.Data()));
586 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
587 }
ea693273 588
75bf77e3 589 Float_t rho = 0;
590 if(externalBackground)rho = externalBackground->GetBackground(0);
591
592
593 // fetch jets
594 TClonesArray *aodJets[2];
ea693273 595 aodJets[0]=0;
5bd732d4 596 if(fAODOut&&!aodJets[0]){
597 aodJets[0] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[0].Data()));
598 aodJets[1] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[1].Data())); }
ea693273 599 if(fAODExtension && !aodJets[0]){
529e2916 600 aodJets[0] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[0].Data()));
601 aodJets[1] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[1].Data())); }
5bd732d4 602 if(fAODIn&&!aodJets[0]){
603 aodJets[0] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[0].Data()));
604 aodJets[1] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[1].Data())); }
605
ea693273 606
529e2916 607 //Double_t ptsub[aodJets[0]->GetEntriesFast()];
608 //Int_t inord[aodJets[0]->GetEntriesFast()];
609 //for(Int_t n=0;n<aodJets[0]->GetEntriesFast();n++){
610 // ptsub[n]=0;
611 // inord[n]=0;}
a9e585a7 612
ea693273 613 TList ParticleList;
614 Int_t nT = GetListOfTracks(&ParticleList);
615 for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
75bf77e3 616 fListJets[iJetType]->Clear();
617 if (!aodJets[iJetType]) continue;
529e2916 618
75bf77e3 619 if(fDebug) Printf("%s: %d jets",fJetBranchName[iJetType].Data(),aodJets[iJetType]->GetEntriesFast());
620
a9e585a7 621
75bf77e3 622 for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
623 AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
529e2916 624 if (jet) fListJets[iJetType]->Add(jet);
625 // if(iJetType==0){
626 // ptsub[iJet]=jet->Pt()-rho*jet->EffectiveAreaCharged();}
627 }}
75bf77e3 628
629 Double_t etabig=0;
630 Double_t ptbig=0;
631 Double_t areabig=0;
632 Double_t phibig=0.;
633 Double_t etasmall=0;
634 Double_t ptsmall=0;
635 Double_t areasmall=0;
da93bb11 636 //Double_t distr=0.;
75bf77e3 637 Double_t phismall=0.;
a9e585a7 638
da93bb11 639
529e2916 640 // Double_t up1[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
641 // Double_t up2[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
642 // Double_t up3[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
643 // Double_t up4[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
644 // Double_t down1[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
645 // Double_t down2[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
646 // Double_t down3[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
647 // Double_t down4[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
648 Int_t iCount=0;
649 Int_t trigJet=-1;
650 Int_t trigBBTrack=-1;
651 Int_t trigInTrack=-1;
652
75bf77e3 653 for(Int_t i=0; i<fListJets[0]->GetEntries(); ++i){
654 AliAODJet* jetbig = (AliAODJet*)(fListJets[0]->At(i));
655 etabig = jetbig->Eta();
656 phibig = jetbig->Phi();
657 ptbig = jetbig->Pt();
658 if(ptbig==0) continue;
659 areabig = jetbig->EffectiveAreaCharged();
ea693273 660 Double_t ptcorr=ptbig-rho*areabig;
661 if(ptcorr<=0) continue;
75bf77e3 662 if((etabig<fJetEtaMin)||(etabig>fJetEtaMax)) continue;
663 Double_t dismin=100.;
664 Double_t ptmax=-10.;
665 Int_t index1=-1;
666 Int_t index2=-1;
da93bb11 667
529e2916 668 Int_t point=GetHardestTrackBackToJet(jetbig);
da93bb11 669 AliVParticle *partback = (AliVParticle*)ParticleList.At(point);
529e2916 670 if(!partback) continue;
da93bb11 671 fh3spectriggered->Fill(centValue,ptcorr,partback->Pt());
6134f035 672 //if(partback->Pt()<6.) continue;
da93bb11 673 AliAODTrack* leadtrack;
674 Int_t ippt=0;
675 Double_t ppt=-10;
529e2916 676 TRefArray *genTrackList = jetbig->GetRefTracks();
677 Int_t nTracksGenJet = genTrackList->GetEntriesFast();
678 AliAODTrack* genTrack;
da93bb11 679 for(Int_t ir=0; ir<nTracksGenJet; ++ir){
680 genTrack = (AliAODTrack*)(genTrackList->At(ir));
681 if(genTrack->Pt()>ppt){ppt=genTrack->Pt();
529e2916 682 ippt=ir;}}
da93bb11 683 leadtrack=(AliAODTrack*)(genTrackList->At(ippt));
529e2916 684 if(!leadtrack) continue;
da93bb11 685 fh3specbiased->Fill(centValue,ptcorr,leadtrack->Pt());
529e2916 686 if(centValue<10)fh3spectot->Fill(ptcorr,leadtrack->Pt(),partback->Pt());
20dcc500 687 if(centValue<20)fh3spectotb->Fill(ptcorr,leadtrack->Pt(),partback->Pt());
529e2916 688 //store one trigger info
689 if((partback->Pt()>10.)&&(iCount==0)){
690 trigJet=i;
691 trigBBTrack=point;
692 trigInTrack=ippt;
693 iCount=iCount+1;}
694
695
696 if(fCheckMethods){
697 for(Int_t j=0; j<fListJets[1]->GetEntries(); ++j){
75bf77e3 698 AliAODJet* jetsmall = (AliAODJet*)(fListJets[1]->At(j));
699 etasmall = jetsmall->Eta();
700 phismall = jetsmall->Phi();
701 ptsmall = jetsmall->Pt();
702 areasmall = jetsmall->EffectiveAreaCharged();
ea693273 703 Double_t tmpDeltaR=(phismall-phibig)*(phismall-phibig)+(etasmall-etabig)*(etasmall-etabig);
704 tmpDeltaR=TMath::Sqrt(tmpDeltaR);
705 //Fraction in the jet core
706 if((ptsmall>ptmax)&&(tmpDeltaR<=fRadioFrac)){ptmax=ptsmall;
75bf77e3 707 index2=j;}
ea693273 708 if(tmpDeltaR<=dismin){ dismin=tmpDeltaR;
709 index1=j;}} //en of loop over R=0.2 jets
75bf77e3 710 //method1:most concentric jet=core
711 if(dismin<fMinDist){ AliAODJet* jetmethod1 = (AliAODJet*)(fListJets[1]->At(index1));
ea693273 712 if(centValue<10) fh2JetCoreMethod1C10->Fill(ptcorr,jetmethod1->Pt()/ptbig);
713 if((centValue>20)&&(centValue<40)) fh2JetCoreMethod1C20->Fill(ptcorr,jetmethod1->Pt()/ptbig);
714 if((centValue>30)&&(centValue<60)) fh2JetCoreMethod1C30->Fill(ptcorr,jetmethod1->Pt()/ptbig);
715 if(centValue>60) fh2JetCoreMethod1C60->Fill(ptcorr,jetmethod1->Pt()/ptbig); }
75bf77e3 716 //method2:hardest contained jet=core
717 if(index2!=-1){
718 AliAODJet* jetmethod2 = (AliAODJet*)(fListJets[1]->At(index2));
ea693273 719 if(centValue<10) fh2JetCoreMethod2C10->Fill(ptcorr,jetmethod2->Pt()/ptbig);
720 if((centValue>20)&&(centValue<40)) fh2JetCoreMethod2C20->Fill(ptcorr,jetmethod2->Pt()/ptbig);
721 if((centValue>30)&&(centValue<60)) fh2JetCoreMethod2C30->Fill(ptcorr,jetmethod2->Pt()/ptbig);
da93bb11 722 if(centValue>60) fh2JetCoreMethod2C60->Fill(ptcorr,jetmethod2->Pt()/ptbig); }}
20dcc500 723 Double_t sumpt2a=0.;
724 Double_t sumpt2b=0.;
5bd732d4 725
20dcc500 726 Double_t sumpt4a=0.;
727 Double_t sumpt4b=0.;
5bd732d4 728
20dcc500 729 Double_t sumpt6a=0.;
730 Double_t sumpt6b=0.;
5bd732d4 731
20dcc500 732 Double_t sumpt8a=0.;
733 Double_t sumpt8b=0.;
5bd732d4 734
20dcc500 735 Double_t sumpt10a=0.;
736 Double_t sumpt10b=0.;
5bd732d4 737
ea693273 738
da93bb11 739 for(int it = 0;it<nT;++it){
ea693273 740 AliVParticle *part = (AliVParticle*)ParticleList.At(it);
529e2916 741 Double_t deltaR = jetbig->DeltaR(part);
5bd732d4 742 Double_t deltaEta = etabig-part->Eta();
20dcc500 743 if(partback->Pt()>10.){
744 if(centValue<10){
20dcc500 745
5bd732d4 746 if(deltaR<0.2){if(leadtrack->Pt()>3.){ sumpt2a=sumpt2a+part->Pt();
747 if(TMath::Abs(deltaEta)>0.8) sumpt2b=sumpt2b+part->Pt();}}
748 if(deltaR>=0.2 && deltaR<0.4){if(leadtrack->Pt()>3.){ sumpt4a=sumpt4a+part->Pt();
749 if(TMath::Abs(deltaEta)>0.8) sumpt4b=sumpt4b+part->Pt();}}
750 if(deltaR>=0.4 && deltaR<0.6){if(leadtrack->Pt()>3.){ sumpt6a=sumpt6a+part->Pt();
751 if(TMath::Abs(deltaEta)>0.8) sumpt6b=sumpt6b+part->Pt();}}
752 if(deltaR>=0.6 && deltaR<0.8){if(leadtrack->Pt()>3.){ sumpt8a=sumpt8a+part->Pt();
753 if(TMath::Abs(deltaEta)>0.8) sumpt8b=sumpt8b+part->Pt();}}
754 if(deltaR>=0.8 && deltaR<1.2){if(leadtrack->Pt()>3.){ sumpt10a=sumpt10a+part->Pt();
755 if(TMath::Abs(deltaEta)>0.8) sumpt10b=sumpt10b+part->Pt();}}}}
756
20dcc500 757
da93bb11 758 Double_t deltaPhi=phibig-part->Phi();
759 if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
760 if(deltaPhi>3./2.*TMath::Pi()) deltaPhi-=2.*TMath::Pi();
529e2916 761 Double_t jetEntries[8] = {centValue,ptcorr,part->Pt(),deltaR,deltaEta,deltaPhi,leadtrack->Pt(),partback->Pt()}; fhnDeltaR->Fill(jetEntries);
762 }
763 //end of track loop
20dcc500 764 Double_t rhoin2=rho*TMath::Pi()*0.2*0.2;
765 Double_t rhoin4=rho*TMath::Pi()*(0.4*0.4-0.2*0.2);
766 Double_t rhoin6=rho*TMath::Pi()*(0.6*0.6-0.4*0.4);
767 Double_t rhoin8=rho*TMath::Pi()*(0.8*0.8-0.6*0.6);
768 Double_t rhoin10=rho*TMath::Pi()*(1.2*1.2-0.8*0.8);
769
770
771 if(rho!=0){
5bd732d4 772 if(partback->Pt()>10.){
773 fh2JetsumHT3R2->Fill(ptcorr,sumpt2a/rhoin2);
774 fh2JetsumHT3R2gap->Fill(ptcorr,sumpt2b/rhoin2);
775 fh2JetsumHT3R4->Fill(ptcorr,sumpt4a/rhoin4);
776 fh2JetsumHT3R4gap->Fill(ptcorr,sumpt4b/rhoin4);
777 fh2JetsumHT3R6->Fill(ptcorr,sumpt6a/rhoin6);
778 fh2JetsumHT3R6gap->Fill(ptcorr,sumpt6b/rhoin6);
779 fh2JetsumHT3R8->Fill(ptcorr,sumpt8a/rhoin8);
780 fh2JetsumHT3R8gap->Fill(ptcorr,sumpt8b/rhoin8);
781 fh2JetsumHT3R10->Fill(ptcorr,sumpt10a/rhoin10);
782 fh2JetsumHT3R10gap->Fill(ptcorr,sumpt10b/rhoin10);
783 }}
20dcc500 784
785 }
529e2916 786 //end of jet loop
787
788
789
790
791 if(fDoEventMixing){
792 //check before if the trigger exists
793 // fTrigBuffer[i][0] = zvtx
794 // fTrigBuffer[i][1] = phi
795 // fTrigBuffer[i][2] = eta
796 // fTrigBuffer[i][3] = pt_jet
797 // fTrigBuffer[i][4] = pt_trig
798 // fTrigBuffer[i][5]= pt_track_in
799 // fTrigBuffer[i][6]= centrality
800 if(fTindex==11) fTindex=0;
801 if(fTrigBuffer[fTindex][3]>0){
802 if (TMath::Abs(fTrigBuffer[fTindex][0]-primVtx->GetZ()<2.)){
803 if (TMath::Abs(fTrigBuffer[fTindex][6]-centValue<10)){
804
805 for(int it = 0;it<nT;++it){
20dcc500 806 AliVParticle *part = (AliVParticle*)ParticleList.At(it);
807 Double_t DPhi = fTrigBuffer[fTindex][1] - part->Phi();
529e2916 808 Double_t DEta = fTrigBuffer[fTindex][2] - part->Eta();
809 Double_t DR=TMath::Sqrt(DPhi*DPhi+DEta*DEta);
810 if(DPhi<-0.5*TMath::Pi()) DPhi+=2.*TMath::Pi();
811 if(DPhi>3./2.*TMath::Pi()) DPhi-=2.*TMath::Pi();
812 Double_t triggerEntries[8] = {centValue,fTrigBuffer[fTindex][3],part->Pt(),DR,DEta,DPhi,fTrigBuffer[fTindex][4],fTrigBuffer[fTindex][5]};
813 fhnMixedEvents->Fill(triggerEntries);
814 }
815 fNevents=fNevents+1;
816 if(fNevents==9) {fTindex=fTindex+1;
817 fNevents=0;}
818 }}}
819
820
821 // Copy the triggers from the current event into the buffer.
822 //again, only if the trigger exists:
823 if(trigJet>-1){
824 AliAODJet* jetT = (AliAODJet*)(fListJets[0]->At(trigJet));
825 AliVParticle *partL = (AliVParticle*)ParticleList.At(trigInTrack);
826 AliVParticle *partT = (AliVParticle*)ParticleList.At(trigBBTrack);
827 fTrigBuffer[fTrigBufferIndex][0] = primVtx->GetZ();
828 fTrigBuffer[fTrigBufferIndex][1] = jetT->Phi();
829 fTrigBuffer[fTrigBufferIndex][2] = jetT->Eta();
830 fTrigBuffer[fTrigBufferIndex][3] = jetT->Pt()-rho*jetT->EffectiveAreaCharged();
831 fTrigBuffer[fTrigBufferIndex][4] = partT->Pt();
832 fTrigBuffer[fTrigBufferIndex][5] = partL->Pt();
833 fTrigBuffer[fTrigBufferIndex][6] = centValue;
834 fTrigBufferIndex++;
835 if(fTrigBufferIndex==9) fTrigBufferIndex=0;
836 }
837 }
838
839
840
841
da93bb11 842
ea693273 843 //////////////////ANGULAR STRUCTURE//////////////////////////////////////
844
845 //tracks up to R=0.8 distant from the jet axis
529e2916 846 // if(fAngStructCloseTracks==1){
847 // TList CloseTrackList;
848 // Int_t nn=GetListOfTracksCloseToJet(&CloseTrackList,jetbig);
849 // Double_t difR=0.04;
850 // for(Int_t l=0;l<15;l++){
851 // Double_t rr=l*0.1+0.1;
852 // for(int it = 0;it<nn;++it){
853 // AliVParticle *part1 = (AliVParticle*)CloseTrackList.At(it);
854 // for(int itu=it+1;itu<CloseTrackList.GetEntries();itu++){
855 // AliVParticle *part2 = (AliVParticle*)CloseTrackList.At(itu);
856 // Double_t ptm=part1->Pt();
857 // Double_t ptn=part2->Pt();
858 // Double_t Rnm = (part1->Eta()-part2->Eta())*(part1->Eta()-part2->Eta())+(part1->Phi()-part2->Phi())*(part1->Phi()-part2->Phi());
859 // Rnm=TMath::Sqrt(Rnm);
860 // Double_t deltag=(1./(TMath::Sqrt(2*TMath::Pi())*difR))*TMath::Exp(-1.*(rr-Rnm)*(rr-Rnm)/(2.*difR*difR));
861 // Double_t stepf=0.5*(1.+TMath::Erf((rr-Rnm)/(TMath::Sqrt(2.)*difR)));
862 // if((ptcorr<85.) && (ptcorr>=70.)){up1[l]=up1[l]+ptm*ptn*Rnm*Rnm*deltag;
863 // down1[l]=down1[l]+ptm*ptn*Rnm*Rnm*stepf;}
864 // if((ptcorr<100.) && (ptcorr>=85.)){up2[l]=up2[l]+ptm*ptn*Rnm*Rnm*deltag;
865 // down2[l]=down2[l]+ptm*ptn*Rnm*Rnm*stepf;}
866 // if((ptcorr<120.) && (ptcorr>=100.)){up3[l]=up3[l]+ptm*ptn*Rnm*Rnm*deltag;
867 // down3[l]=down3[l]+ptm*ptn*Rnm*Rnm*stepf;}
868 // if((ptcorr<140.) && (ptcorr>=120.)){up4[l]=up4[l]+ptm*ptn*Rnm*Rnm*deltag;
869 // down4[l]=down4[l]+ptm*ptn*Rnm*Rnm*stepf;}}}}
870 // }
ea693273 871
529e2916 872 // //only jet constituents
873 // if(fAngStructCloseTracks==2){
ea693273 874
529e2916 875 // Double_t difR=0.04;
876 // for(Int_t l=0;l<15;l++){
877 // Double_t rr=l*0.1+0.1;
ea693273 878
879
529e2916 880 // AliAODTrack* part1;
881 // AliAODTrack* part2;
da93bb11 882
529e2916 883 // TRefArray *genTrackListb = jetbig->GetRefTracks();
884 // Int_t nTracksGenJetb = genTrackListb->GetEntriesFast();
da93bb11 885
886
887
529e2916 888 // for(Int_t it=0; it<nTracksGenJetb; ++it){
889 // part1 = (AliAODTrack*)(genTrackListb->At(it));
890 // for(Int_t itu=0; itu<nTracksGenJetb; ++itu){
891 // part2 = (AliAODTrack*)(genTrackListb->At(itu));
892 // Double_t ptm=part1->Pt();
893 // Double_t ptn=part2->Pt();
894 // Double_t Rnm = (part1->Eta()-part2->Eta())*(part1->Eta()-part2->Eta())+(part1->Phi()-part2->Phi())*(part1->Phi()-part2->Phi());
895 // Rnm=TMath::Sqrt(Rnm);
896 // Double_t deltag=(1./(TMath::Sqrt(2*TMath::Pi())*difR))*TMath::Exp(-1.*(rr-Rnm)*(rr-Rnm)/(2.*difR*difR));
897 // Double_t stepf=0.5*(1.+TMath::Erf((rr-Rnm)/(TMath::Sqrt(2.)*difR)));
898 // if((ptcorr<85.) && (ptcorr>=70.)){up1[l]=up1[l]+ptm*ptn*Rnm*Rnm*deltag;
899 // down1[l]=down1[l]+ptm*ptn*Rnm*Rnm*stepf;}
900 // if((ptcorr<100.) && (ptcorr>=85.)){up2[l]=up2[l]+ptm*ptn*Rnm*Rnm*deltag;
901 // down2[l]=down2[l]+ptm*ptn*Rnm*Rnm*stepf;}
902 // if((ptcorr<120.) && (ptcorr>=100.)){up3[l]=up3[l]+ptm*ptn*Rnm*Rnm*deltag;
903 // down3[l]=down3[l]+ptm*ptn*Rnm*Rnm*stepf;}
904 // if((ptcorr<140.) && (ptcorr>=120.)){up4[l]=up4[l]+ptm*ptn*Rnm*Rnm*deltag;
905 // down4[l]=down4[l]+ptm*ptn*Rnm*Rnm*stepf;}}}}}
906 // }
907 // //end loop over R=0.4 jets
908 // if(fAngStructCloseTracks>0){
909 // for(Int_t l=0;l<15;l++){
910 // Double_t rr=l*0.1+0.1;
911 // if(down1[l]!=0){
912 // if(centValue<10.)fh2AngStructpt1C10->Fill(rr,rr*up1[l]/down1[l]);
913 // if(centValue>20. && centValue<40.) fh2AngStructpt1C20->Fill(rr,rr*up1[l]/down1[l]);
914 // if(centValue>30. && centValue<60.) fh2AngStructpt1C30->Fill(rr,rr*up1[l]/down1[l]);
915 // if(centValue>60.) fh2AngStructpt1C60->Fill(rr,rr*up1[l]/down1[l]);}
916 // if(down2[l]!=0){
917 // if(centValue<10.) fh2AngStructpt2C10->Fill(rr,rr*up2[l]/down2[l]);
918 // if(centValue>20. && centValue<40.) fh2AngStructpt2C20->Fill(rr,rr*up2[l]/down2[l]);
919 // if(centValue>30. && centValue<60.) fh2AngStructpt2C30->Fill(rr,rr*up2[l]/down2[l]);
920 // if(centValue>60.) fh2AngStructpt2C60->Fill(rr,rr*up2[l]/down2[l]);}
921 // if(down3[l]!=0){
922 // if(centValue<10.) fh2AngStructpt3C10->Fill(rr,rr*up3[l]/down3[l]);
923 // if(centValue>20. && centValue<40.) fh2AngStructpt3C20->Fill(rr,rr*up3[l]/down3[l]);
924 // if(centValue>30. && centValue<60.) fh2AngStructpt3C30->Fill(rr,rr*up3[l]/down3[l]);
925 // if(centValue>60.) fh2AngStructpt3C60->Fill(rr,rr*up3[l]/down3[l]);}
926 // if(down4[l]!=0){
927 // if(centValue<10.) fh2AngStructpt4C10->Fill(rr,rr*up4[l]/down4[l]);
928 // if(centValue>20. && centValue<40.) fh2AngStructpt4C20->Fill(rr,rr*up4[l]/down4[l]);
929 // if(centValue>30. && centValue<60.) fh2AngStructpt4C30->Fill(rr,rr*up4[l]/down4[l]);
930 // if(centValue>60.) fh2AngStructpt4C60->Fill(rr,rr*up4[l]/down4[l]);}}}
ea693273 931
932
933
934
935
75bf77e3 936
937
938 PostData(1, fOutputList);
da93bb11 939}
75bf77e3 940
941void AliAnalysisTaskJetCore::Terminate(const Option_t *)
942{
943 // Draw result to the screen
944 // Called once at the end of the query
945
946 if (!GetOutputData(1))
947 return;
948}
949
ea693273 950
951
5bd732d4 952
ea693273 953
954
955Int_t AliAnalysisTaskJetCore::GetListOfTracks(TList *list){
956
5bd732d4 957 Int_t iCount = 0;
958 AliAODEvent *aod = 0;
542c9952 959 if(!fESD){
960 if(!fAODIn)return iCount;
961 aod = fAODIn;
962 }
5bd732d4 963 else aod = fAODOut;
542c9952 964
5bd732d4 965 for(int it = 0;it < aod->GetNumberOfTracks();++it){
966 AliAODTrack *tr = aod->GetTrack(it);
ea693273 967 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
968 if(TMath::Abs(tr->Eta())>0.9)continue;
969 if(tr->Pt()<0.15)continue;
970 list->Add(tr);
971 //cout<<fAOD->GetNumberOfTracks()<<" "<<tr->Pt()<<endl;
972 iCount++;
973 }
974
529e2916 975
ea693273 976 return iCount;
977
978}
979
da93bb11 980 Int_t AliAnalysisTaskJetCore::GetHardestTrackBackToJet(AliAODJet *jetbig){
5bd732d4 981
982 AliAODEvent *aod = 0;
983 if(!fESD)aod = fAODIn;
984 else aod = fAODOut;
da93bb11 985 Int_t index=-1;
986 Double_t ptmax=-10;
987 Double_t dphi=0;
988 Double_t dif=0;
989 Int_t iCount=0;
5bd732d4 990 for(int it = 0;it < aod->GetNumberOfTracks();++it){
991 AliAODTrack *tr = aod->GetTrack(it);
da93bb11 992 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
993 if(TMath::Abs(tr->Eta())>0.9)continue;
994 if(tr->Pt()<0.15)continue;
995 iCount=iCount+1;
529e2916 996 dphi=RelativePhi(tr->Phi(),jetbig->Phi());
da93bb11 997 if(TMath::Abs(dphi)<TMath::Pi()-0.2) continue;
998 if(tr->Pt()>ptmax){ ptmax=tr->Pt();
b6844e9e 999 index=iCount-1;
1000 dif=dphi; }}
da93bb11 1001
1002 return index;
1003
1004 }
1005
1006
1007
1008
1009
1010
1011
1012
1013
ea693273 1014 Int_t AliAnalysisTaskJetCore::GetListOfTracksCloseToJet(TList *list,AliAODJet *jetbig){
1015
1016 Int_t iCount = 0;
5bd732d4 1017 AliAODEvent *aod = 0;
1018 if(!fESD)aod = fAODIn;
1019 else aod = fAODOut;
8b47ec90 1020
5bd732d4 1021 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1022 AliAODTrack *tr = aod->GetTrack(it);
ea693273 1023 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1024 if(TMath::Abs(tr->Eta())>0.9)continue;
1025 if(tr->Pt()<0.15)continue;
1026 Double_t disR=jetbig->DeltaR(tr);
1027 if(disR>0.8) continue;
1028 list->Add(tr);
1029 //cout<<fAOD->GetNumberOfTracks()<<" "<<tr->Pt()<<endl;
1030 iCount++;
1031 }
1032
1033 list->Sort();
1034 return iCount;
1035
1036}
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
75bf77e3 1048Int_t AliAnalysisTaskJetCore::GetNInputTracks()
1049{
1050
1051 Int_t nInputTracks = 0;
5bd732d4 1052 AliAODEvent *aod = 0;
1053 if(!fESD)aod = fAODIn;
1054 else aod = fAODOut;
75bf77e3 1055 TString jbname(fJetBranchName[1]);
1056 //needs complete event, use jets without background subtraction
1057 for(Int_t i=1; i<=3; ++i){
1058 if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0");
1059 }
1060 // use only HI event
1061 if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
1062 if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
1063
1064 if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
5bd732d4 1065 TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(aod->FindListObject(jbname.Data()));
75bf77e3 1066 if(!tmpAODjets){
1067 Printf("Jet branch %s not found", jbname.Data());
1068 Printf("AliAnalysisTaskJetCore::GetNInputTracks FAILED");
1069 return -1;
1070 }
1071
1072 for (Int_t iJet=0; iJet<tmpAODjets->GetEntriesFast(); iJet++){
1073 AliAODJet *jet = dynamic_cast<AliAODJet*>((*tmpAODjets)[iJet]);
1074 if(!jet) continue;
1075 TRefArray *trackList = jet->GetRefTracks();
1076 Int_t nTracks = trackList->GetEntriesFast();
1077 nInputTracks += nTracks;
1078 if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks);
1079 }
1080 if(fDebug) Printf("---> input tracks: %d", nInputTracks);
1081
1082 return nInputTracks;
1083}
1084
1085
1086
ea693273 1087Double_t AliAnalysisTaskJetCore::RelativePhi(Double_t mphi,Double_t vphi){
1088
1089 if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
1090 else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
1091 if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
1092 else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
1093 double dphi = mphi-vphi;
1094 if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
1095 else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
1096 return dphi;//dphi in [-Pi, Pi]
1097}
1098
75bf77e3 1099
1100
da93bb11 1101THnSparse* AliAnalysisTaskJetCore::NewTHnSparseF(const char* name, UInt_t entries)
1102{
1103 // generate new THnSparseF, axes are defined in GetDimParams()
1104
1105 Int_t count = 0;
1106 UInt_t tmp = entries;
1107 while(tmp!=0){
1108 count++;
1109 tmp = tmp &~ -tmp; // clear lowest bit
1110 }
1111
1112 TString hnTitle(name);
1113 const Int_t dim = count;
1114 Int_t nbins[dim];
1115 Double_t xmin[dim];
1116 Double_t xmax[dim];
1117
1118 Int_t i=0;
1119 Int_t c=0;
1120 while(c<dim && i<32){
1121 if(entries&(1<<i)){
1122
1123 TString label("");
1124 GetDimParams(i, label, nbins[c], xmin[c], xmax[c]);
1125 hnTitle += Form(";%s",label.Data());
1126 c++;
1127 }
1128
1129 i++;
1130 }
1131 hnTitle += ";";
1132
1133 return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
1134}
1135
1136void AliAnalysisTaskJetCore::GetDimParams(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
1137{
1138 // stores label and binning of axis for THnSparse
1139
1140 const Double_t pi = TMath::Pi();
1141
1142 switch(iEntry){
1143
1144 case 0:
1145 label = "V0 centrality (%)";
1146
1147 nbins = 10;
1148 xmin = 0.;
1149 xmax = 100.;
1150 break;
1151
1152
1153 case 1:
1154 label = "corrected jet pt";
c42c6838 1155 nbins = 20;
da93bb11 1156 xmin = 0.;
1157 xmax = 200.;
1158 break;
1159
1160
1161 case 2:
1162 label = "track pT";
1163
5bd732d4 1164 nbins = 10;
da93bb11 1165 xmin = 0.;
5bd732d4 1166 xmax = 150;
da93bb11 1167 break;
1168
1169
529e2916 1170 case 3:
da93bb11 1171 label = "deltaR";
1172 nbins = 15;
1173 xmin = 0.;
1174 xmax = 1.5;
1175 break;
529e2916 1176
1177
1178
da93bb11 1179 case 4:
1180 label = "deltaEta";
5bd732d4 1181 nbins = 15;
1182 xmin = -1.5;
1183 xmax = 1.5;
da93bb11 1184 break;
1185
1186
5bd732d4 1187 case 5:
da93bb11 1188 label = "deltaPhi";
1189 nbins = 90;
1190 xmin = -0.5*pi;
1191 xmax = 1.5*pi;
1192 break;
1193
1194
529e2916 1195
1196 case 6:
da93bb11 1197 label = "leading track";
5bd732d4 1198 nbins = 14;
da93bb11 1199 xmin = 0;
20dcc500 1200 xmax = 50;
da93bb11 1201 break;
1202
529e2916 1203 case 7:
da93bb11 1204
1205 label = "trigger track";
c42c6838 1206 nbins =5;
da93bb11 1207 xmin = 0;
1208 xmax = 50;
1209 break;
529e2916 1210
1211
1212
1213
1214
1215
1216
1217
da93bb11 1218 }
1219
1220}
1221