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