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