- fix init of trackfilter
[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
dfaeec0d 788
789
16a4c8e8 790 AliVParticle* leadtrackb=0;
791 if(fFlagJetHadron!=0){
792 Int_t nTb = GetHardestTrackBackToJet(jetbig);
793 leadtrackb = (AliVParticle*)ParticleList.At(nTb);
794 if(!leadtrackb) continue;
795 }
796
797
798
799
d5e7475c 800
0bf6381d 801 //store one trigger info
802 if(iCount==0){
803 trigJet=i;
804 trigBBTrack=nT;
805 trigInTrack=ippt;
806 iCount=iCount+1;}
807
d5e7475c 808
0bf6381d 809 if(fCheckMethods){
810 for(Int_t j=0; j<fListJets[1]->GetEntries(); ++j){
811 AliAODJet* jetsmall = (AliAODJet*)(fListJets[1]->At(j));
812 etasmall = jetsmall->Eta();
813 phismall = jetsmall->Phi();
814 ptsmall = jetsmall->Pt();
815 areasmall = jetsmall->EffectiveAreaCharged();
816 Double_t tmpDeltaR=(phismall-phibig)*(phismall-phibig)+(etasmall-etabig)*(etasmall-etabig);
817 tmpDeltaR=TMath::Sqrt(tmpDeltaR);
818 //Fraction in the jet core
819 if((ptsmall>ptmax)&&(tmpDeltaR<=fRadioFrac)){ptmax=ptsmall;
820 index2=j;}
ea693273 821 if(tmpDeltaR<=dismin){ dismin=tmpDeltaR;
822 index1=j;}} //en of loop over R=0.2 jets
75bf77e3 823 //method1:most concentric jet=core
824 if(dismin<fMinDist){ AliAODJet* jetmethod1 = (AliAODJet*)(fListJets[1]->At(index1));
ea693273 825 if(centValue<10) fh2JetCoreMethod1C10->Fill(ptcorr,jetmethod1->Pt()/ptbig);
826 if((centValue>20)&&(centValue<40)) fh2JetCoreMethod1C20->Fill(ptcorr,jetmethod1->Pt()/ptbig);
827 if((centValue>30)&&(centValue<60)) fh2JetCoreMethod1C30->Fill(ptcorr,jetmethod1->Pt()/ptbig);
828 if(centValue>60) fh2JetCoreMethod1C60->Fill(ptcorr,jetmethod1->Pt()/ptbig); }
75bf77e3 829 //method2:hardest contained jet=core
830 if(index2!=-1){
831 AliAODJet* jetmethod2 = (AliAODJet*)(fListJets[1]->At(index2));
ea693273 832 if(centValue<10) fh2JetCoreMethod2C10->Fill(ptcorr,jetmethod2->Pt()/ptbig);
833 if((centValue>20)&&(centValue<40)) fh2JetCoreMethod2C20->Fill(ptcorr,jetmethod2->Pt()/ptbig);
834 if((centValue>30)&&(centValue<60)) fh2JetCoreMethod2C30->Fill(ptcorr,jetmethod2->Pt()/ptbig);
da93bb11 835 if(centValue>60) fh2JetCoreMethod2C60->Fill(ptcorr,jetmethod2->Pt()/ptbig); }}
dfaeec0d 836 if(centValue<10&&leadtrack) fh2Ntriggers2C10->Fill(leadtrack->Pt(),partback->Pt());
837 if(centValue<20&&leadtrack) fh2Ntriggers2C20->Fill(leadtrack->Pt(),partback->Pt());
d90d5d75 838 if(fDoEventMixing==0 && fFlagOnlyRecoil==0){
41bcb7e7 839 for(int it = 0;it<ParticleList.GetEntries();++it){
840 AliVParticle *part = (AliVParticle*)ParticleList.At(it);
841 Double_t deltaR = jetbig->DeltaR(part);
842 Double_t deltaEta = etabig-part->Eta();
d90d5d75 843
41bcb7e7 844 Double_t deltaPhi=phibig-part->Phi();
845 if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
846 if(deltaPhi>3./2.*TMath::Pi()) deltaPhi-=2.*TMath::Pi();
847 Double_t pTcont=0;
848 if(fFlagJetHadron==0) pTcont=leadtrack->Pt();
849 if(fFlagJetHadron!=0) pTcont=leadtrackb->Pt();
850 Double_t jetEntries[8] = {centValue,ptcorr,part->Pt(),deltaR,deltaEta,deltaPhi,pTcont,partback->Pt()};
851 fhnDeltaR->Fill(jetEntries);}
852
853
d90d5d75 854 }
85b5b73e 855
d5e7475c 856
d5e7475c 857 //end of track loop, we only do it if EM is switched off
858
ba143e7f 859
860
861
862
863
864
865
20dcc500 866
867 }
ae1e07c1 868 if(injet>0) fh3JetDensity->Fill(ParticleList.GetEntries(),injet/accep,partback->Pt());
869 if(injet4>0)fh3JetDensityA4->Fill(ParticleList.GetEntries(),injet4/accep,partback->Pt());
529e2916 870 //end of jet loop
871
41bcb7e7 872 //}
529e2916 873
874
447d2a5b 875 if(fDoEventMixing>0){
529e2916 876 //check before if the trigger exists
877 // fTrigBuffer[i][0] = zvtx
878 // fTrigBuffer[i][1] = phi
879 // fTrigBuffer[i][2] = eta
880 // fTrigBuffer[i][3] = pt_jet
881 // fTrigBuffer[i][4] = pt_trig
447d2a5b 882 // fTrigBuffer[i][5]= centrality
883 if(fTindex==10) fTindex=0;
529e2916 884 if(fTrigBuffer[fTindex][3]>0){
885 if (TMath::Abs(fTrigBuffer[fTindex][0]-primVtx->GetZ()<2.)){
447d2a5b 886 if (TMath::Abs(fTrigBuffer[fTindex][5]-centValue<5)){
529e2916 887
888 for(int it = 0;it<nT;++it){
20dcc500 889 AliVParticle *part = (AliVParticle*)ParticleList.At(it);
890 Double_t DPhi = fTrigBuffer[fTindex][1] - part->Phi();
529e2916 891 Double_t DEta = fTrigBuffer[fTindex][2] - part->Eta();
892 Double_t DR=TMath::Sqrt(DPhi*DPhi+DEta*DEta);
893 if(DPhi<-0.5*TMath::Pi()) DPhi+=2.*TMath::Pi();
894 if(DPhi>3./2.*TMath::Pi()) DPhi-=2.*TMath::Pi();
447d2a5b 895 Double_t triggerEntries[7] = {centValue,fTrigBuffer[fTindex][3],part->Pt(),DR,DEta,DPhi,fTrigBuffer[fTindex][4]};
529e2916 896 fhnMixedEvents->Fill(triggerEntries);
897 }
898 fNevents=fNevents+1;
447d2a5b 899 if(fNevents==10) fTindex=fTindex+1;
529e2916 900 }}}
ae1e07c1 901
447d2a5b 902 if(fTindex==10&&fNevents==10) fCountAgain=0;
529e2916 903
904 // Copy the triggers from the current event into the buffer.
905 //again, only if the trigger exists:
447d2a5b 906 if(fCountAgain==0){
529e2916 907 if(trigJet>-1){
447d2a5b 908 AliAODJet* jetT = (AliAODJet*)(fListJets[0]->At(trigJet)); AliVParticle *partT = (AliVParticle*)ParticleList.At(trigBBTrack);
529e2916 909 fTrigBuffer[fTrigBufferIndex][0] = primVtx->GetZ();
910 fTrigBuffer[fTrigBufferIndex][1] = jetT->Phi();
911 fTrigBuffer[fTrigBufferIndex][2] = jetT->Eta();
912 fTrigBuffer[fTrigBufferIndex][3] = jetT->Pt()-rho*jetT->EffectiveAreaCharged();
913 fTrigBuffer[fTrigBufferIndex][4] = partT->Pt();
447d2a5b 914 fTrigBuffer[fTrigBufferIndex][5] = centValue;
529e2916 915 fTrigBufferIndex++;
447d2a5b 916 if(fTrigBufferIndex==9) {fTrigBufferIndex=0;
917 fCountAgain=1;}
529e2916 918 }
447d2a5b 919 }
529e2916 920
447d2a5b 921 }
529e2916 922
4f6d1cc1 923 /////////////////////////////////////////////////////////////////////////////
924 ////////////////////// Rongrong's analysis //////////////////////////////////
925 if(fRunAnaAzimuthalCorrelation)
926 {
927 fhTTPt->Fill(centValue,partback->Pt());
928 for(Int_t ij=0; ij<fListJets[0]->GetEntries(); ij++)
929 {
930 AliAODJet* jet = (AliAODJet*)(fListJets[0]->At(ij));
931 Double_t jetPt = jet->Pt();
932 Double_t jetEta = jet->Eta();
933 Double_t jetPhi = jet->Phi();
934 if(jetPt==0) continue;
935 if((jetEta<fJetEtaMin)||(jetEta>fJetEtaMax)) continue;
936 Double_t jetArea = jet->EffectiveAreaCharged();
937 Double_t jetPtCorr=jetPt-rho*jetArea;
830b5724 938 Double_t dPhi=jetPhi-partback->Phi();
4f6d1cc1 939 if(dPhi>2*TMath::Pi()) dPhi -= 2*TMath::Pi();
830b5724 940 if(dPhi<-2*TMath::Pi()) dPhi += 2*TMath::Pi();
941 if(dPhi<-0.5*TMath::Pi()) dPhi += 2*TMath::Pi();
942 if(dPhi>1.5*TMath::Pi()) dPhi -= 2*TMath::Pi();
943
4f6d1cc1 944 Double_t fill[] = {partback->Pt(),jetPtCorr,dPhi,jetArea,centValue};
945 fHJetPhiCorr->Fill(fill);
946 }
947 }
948 /////////////////////////////////////////////////////////////////////////////
949 /////////////////////////////////////////////////////////////////////////////
529e2916 950
da93bb11 951
ea693273 952 //////////////////ANGULAR STRUCTURE//////////////////////////////////////
953
954 //tracks up to R=0.8 distant from the jet axis
529e2916 955 // if(fAngStructCloseTracks==1){
956 // TList CloseTrackList;
957 // Int_t nn=GetListOfTracksCloseToJet(&CloseTrackList,jetbig);
958 // Double_t difR=0.04;
959 // for(Int_t l=0;l<15;l++){
960 // Double_t rr=l*0.1+0.1;
961 // for(int it = 0;it<nn;++it){
962 // AliVParticle *part1 = (AliVParticle*)CloseTrackList.At(it);
963 // for(int itu=it+1;itu<CloseTrackList.GetEntries();itu++){
964 // AliVParticle *part2 = (AliVParticle*)CloseTrackList.At(itu);
965 // Double_t ptm=part1->Pt();
966 // Double_t ptn=part2->Pt();
967 // Double_t Rnm = (part1->Eta()-part2->Eta())*(part1->Eta()-part2->Eta())+(part1->Phi()-part2->Phi())*(part1->Phi()-part2->Phi());
968 // Rnm=TMath::Sqrt(Rnm);
969 // Double_t deltag=(1./(TMath::Sqrt(2*TMath::Pi())*difR))*TMath::Exp(-1.*(rr-Rnm)*(rr-Rnm)/(2.*difR*difR));
970 // Double_t stepf=0.5*(1.+TMath::Erf((rr-Rnm)/(TMath::Sqrt(2.)*difR)));
971 // if((ptcorr<85.) && (ptcorr>=70.)){up1[l]=up1[l]+ptm*ptn*Rnm*Rnm*deltag;
972 // down1[l]=down1[l]+ptm*ptn*Rnm*Rnm*stepf;}
973 // if((ptcorr<100.) && (ptcorr>=85.)){up2[l]=up2[l]+ptm*ptn*Rnm*Rnm*deltag;
974 // down2[l]=down2[l]+ptm*ptn*Rnm*Rnm*stepf;}
975 // if((ptcorr<120.) && (ptcorr>=100.)){up3[l]=up3[l]+ptm*ptn*Rnm*Rnm*deltag;
976 // down3[l]=down3[l]+ptm*ptn*Rnm*Rnm*stepf;}
977 // if((ptcorr<140.) && (ptcorr>=120.)){up4[l]=up4[l]+ptm*ptn*Rnm*Rnm*deltag;
978 // down4[l]=down4[l]+ptm*ptn*Rnm*Rnm*stepf;}}}}
979 // }
ea693273 980
529e2916 981 // //only jet constituents
982 // if(fAngStructCloseTracks==2){
ea693273 983
529e2916 984 // Double_t difR=0.04;
985 // for(Int_t l=0;l<15;l++){
986 // Double_t rr=l*0.1+0.1;
ea693273 987
988
529e2916 989 // AliAODTrack* part1;
990 // AliAODTrack* part2;
da93bb11 991
529e2916 992 // TRefArray *genTrackListb = jetbig->GetRefTracks();
993 // Int_t nTracksGenJetb = genTrackListb->GetEntriesFast();
da93bb11 994
995
996
529e2916 997 // for(Int_t it=0; it<nTracksGenJetb; ++it){
998 // part1 = (AliAODTrack*)(genTrackListb->At(it));
999 // for(Int_t itu=0; itu<nTracksGenJetb; ++itu){
1000 // part2 = (AliAODTrack*)(genTrackListb->At(itu));
1001 // Double_t ptm=part1->Pt();
1002 // Double_t ptn=part2->Pt();
1003 // Double_t Rnm = (part1->Eta()-part2->Eta())*(part1->Eta()-part2->Eta())+(part1->Phi()-part2->Phi())*(part1->Phi()-part2->Phi());
1004 // Rnm=TMath::Sqrt(Rnm);
1005 // Double_t deltag=(1./(TMath::Sqrt(2*TMath::Pi())*difR))*TMath::Exp(-1.*(rr-Rnm)*(rr-Rnm)/(2.*difR*difR));
1006 // Double_t stepf=0.5*(1.+TMath::Erf((rr-Rnm)/(TMath::Sqrt(2.)*difR)));
1007 // if((ptcorr<85.) && (ptcorr>=70.)){up1[l]=up1[l]+ptm*ptn*Rnm*Rnm*deltag;
1008 // down1[l]=down1[l]+ptm*ptn*Rnm*Rnm*stepf;}
1009 // if((ptcorr<100.) && (ptcorr>=85.)){up2[l]=up2[l]+ptm*ptn*Rnm*Rnm*deltag;
1010 // down2[l]=down2[l]+ptm*ptn*Rnm*Rnm*stepf;}
1011 // if((ptcorr<120.) && (ptcorr>=100.)){up3[l]=up3[l]+ptm*ptn*Rnm*Rnm*deltag;
1012 // down3[l]=down3[l]+ptm*ptn*Rnm*Rnm*stepf;}
1013 // if((ptcorr<140.) && (ptcorr>=120.)){up4[l]=up4[l]+ptm*ptn*Rnm*Rnm*deltag;
1014 // down4[l]=down4[l]+ptm*ptn*Rnm*Rnm*stepf;}}}}}
1015 // }
1016 // //end loop over R=0.4 jets
1017 // if(fAngStructCloseTracks>0){
1018 // for(Int_t l=0;l<15;l++){
1019 // Double_t rr=l*0.1+0.1;
1020 // if(down1[l]!=0){
1021 // if(centValue<10.)fh2AngStructpt1C10->Fill(rr,rr*up1[l]/down1[l]);
1022 // if(centValue>20. && centValue<40.) fh2AngStructpt1C20->Fill(rr,rr*up1[l]/down1[l]);
1023 // if(centValue>30. && centValue<60.) fh2AngStructpt1C30->Fill(rr,rr*up1[l]/down1[l]);
1024 // if(centValue>60.) fh2AngStructpt1C60->Fill(rr,rr*up1[l]/down1[l]);}
1025 // if(down2[l]!=0){
1026 // if(centValue<10.) fh2AngStructpt2C10->Fill(rr,rr*up2[l]/down2[l]);
1027 // if(centValue>20. && centValue<40.) fh2AngStructpt2C20->Fill(rr,rr*up2[l]/down2[l]);
1028 // if(centValue>30. && centValue<60.) fh2AngStructpt2C30->Fill(rr,rr*up2[l]/down2[l]);
1029 // if(centValue>60.) fh2AngStructpt2C60->Fill(rr,rr*up2[l]/down2[l]);}
1030 // if(down3[l]!=0){
1031 // if(centValue<10.) fh2AngStructpt3C10->Fill(rr,rr*up3[l]/down3[l]);
1032 // if(centValue>20. && centValue<40.) fh2AngStructpt3C20->Fill(rr,rr*up3[l]/down3[l]);
1033 // if(centValue>30. && centValue<60.) fh2AngStructpt3C30->Fill(rr,rr*up3[l]/down3[l]);
1034 // if(centValue>60.) fh2AngStructpt3C60->Fill(rr,rr*up3[l]/down3[l]);}
1035 // if(down4[l]!=0){
1036 // if(centValue<10.) fh2AngStructpt4C10->Fill(rr,rr*up4[l]/down4[l]);
1037 // if(centValue>20. && centValue<40.) fh2AngStructpt4C20->Fill(rr,rr*up4[l]/down4[l]);
1038 // if(centValue>30. && centValue<60.) fh2AngStructpt4C30->Fill(rr,rr*up4[l]/down4[l]);
1039 // if(centValue>60.) fh2AngStructpt4C60->Fill(rr,rr*up4[l]/down4[l]);}}}
ea693273 1040
1041
1042
1043
1044
75bf77e3 1045
1046
1047 PostData(1, fOutputList);
da93bb11 1048}
75bf77e3 1049
1050void AliAnalysisTaskJetCore::Terminate(const Option_t *)
1051{
1052 // Draw result to the screen
1053 // Called once at the end of the query
1054
1055 if (!GetOutputData(1))
1056 return;
1057}
1058
ea693273 1059
1060
5bd732d4 1061
ea693273 1062
1063
1064Int_t AliAnalysisTaskJetCore::GetListOfTracks(TList *list){
1065
85b5b73e 1066 Int_t iCount = 0;
1067 AliAODEvent *aod = 0;
1068
ba143e7f 1069 if(!fESD)aod = fAODIn;
5bd732d4 1070 else aod = fAODOut;
0bf6381d 1071
1072 if(!aod)return 0;
1073
1074 Int_t index=-1;
8205e054 1075 Double_t ptmax=-10;
85b5b73e 1076
1077
1078
1079 for(int it = 0;it < aod->GetNumberOfTracks();++it){
5bd732d4 1080 AliAODTrack *tr = aod->GetTrack(it);
ae73f4c9 1081 Bool_t bGood = false;
1082 if(fFilterType == 0)bGood = true;
1083 else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
d90d5d75 1084 else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
1085 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
ae73f4c9 1086 if(bGood==false) continue;
ea693273 1087 if(TMath::Abs(tr->Eta())>0.9)continue;
1088 if(tr->Pt()<0.15)continue;
1089 list->Add(tr);
ea693273 1090 iCount++;
d90d5d75 1091 if(fFilterType==2 && fFilterMaskBestPt>0){// only set the trigger track index for good quality tracks
873306cf 1092 if(tr->TestFilterBit(fFilterMaskBestPt)){
1093 if(tr->Pt()>ptmax){
1094 ptmax=tr->Pt();
1095 index=iCount-1;
1096 }
1097 }
1098 }
1099 else{
1100 if(tr->Pt()>ptmax){
1101 ptmax=tr->Pt();
1102 index=iCount-1;
1103 }
1104 }
85b5b73e 1105 }
ea693273 1106
529e2916 1107
85b5b73e 1108 // else if (type == kTrackAODMCCharged) {
1109 // TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1110 // if(!tca)return iCount;
1111 // for(int it = 0;it < tca->GetEntriesFast();++it){
1112 // AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1113 // if(!part)continue;
1114 // if(part->Pt()<0.15)continue;
1115 // if(!part->IsPhysicalPrimary())continue;
1116 // if(part->Charge()==0)continue;
1117 // if(TMath::Abs(part->Eta())>0.9)continue;
1118 // list->Add(part);
1119 // iCount++;
1120 // if(part->Pt()>ptmax){ ptmax=part->Pt();
1121 // index=iCount-1;}}}
1122 return index;
ea693273 1123
1124}
1125
da93bb11 1126 Int_t AliAnalysisTaskJetCore::GetHardestTrackBackToJet(AliAODJet *jetbig){
5bd732d4 1127
1128 AliAODEvent *aod = 0;
1129 if(!fESD)aod = fAODIn;
1130 else aod = fAODOut;
da93bb11 1131 Int_t index=-1;
1132 Double_t ptmax=-10;
1133 Double_t dphi=0;
1134 Double_t dif=0;
1135 Int_t iCount=0;
5bd732d4 1136 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1137 AliAODTrack *tr = aod->GetTrack(it);
da93bb11 1138 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1139 if(TMath::Abs(tr->Eta())>0.9)continue;
1140 if(tr->Pt()<0.15)continue;
1141 iCount=iCount+1;
529e2916 1142 dphi=RelativePhi(tr->Phi(),jetbig->Phi());
16a4c8e8 1143 if(TMath::Abs(dphi)<TMath::Pi()-0.6) continue;
da93bb11 1144 if(tr->Pt()>ptmax){ ptmax=tr->Pt();
b6844e9e 1145 index=iCount-1;
1146 dif=dphi; }}
da93bb11 1147
1148 return index;
1149
1150 }
1151
1152
1153
1154
1155
1156
1157
1158
1159
ea693273 1160 Int_t AliAnalysisTaskJetCore::GetListOfTracksCloseToJet(TList *list,AliAODJet *jetbig){
1161
1162 Int_t iCount = 0;
5bd732d4 1163 AliAODEvent *aod = 0;
1164 if(!fESD)aod = fAODIn;
1165 else aod = fAODOut;
8b47ec90 1166
16a4c8e8 1167 for(int it = 0;it < aod->GetNumberOfTracks();++it){
5bd732d4 1168 AliAODTrack *tr = aod->GetTrack(it);
ea693273 1169 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1170 if(TMath::Abs(tr->Eta())>0.9)continue;
1171 if(tr->Pt()<0.15)continue;
1172 Double_t disR=jetbig->DeltaR(tr);
1173 if(disR>0.8) continue;
1174 list->Add(tr);
1175 //cout<<fAOD->GetNumberOfTracks()<<" "<<tr->Pt()<<endl;
1176 iCount++;
1177 }
1178
1179 list->Sort();
1180 return iCount;
1181
1182}
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
75bf77e3 1194Int_t AliAnalysisTaskJetCore::GetNInputTracks()
1195{
1196
1197 Int_t nInputTracks = 0;
5bd732d4 1198 AliAODEvent *aod = 0;
1199 if(!fESD)aod = fAODIn;
1200 else aod = fAODOut;
75bf77e3 1201 TString jbname(fJetBranchName[1]);
1202 //needs complete event, use jets without background subtraction
1203 for(Int_t i=1; i<=3; ++i){
1204 if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0");
1205 }
1206 // use only HI event
1207 if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
1208 if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
1209
1210 if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
5bd732d4 1211 TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(aod->FindListObject(jbname.Data()));
75bf77e3 1212 if(!tmpAODjets){
1213 Printf("Jet branch %s not found", jbname.Data());
1214 Printf("AliAnalysisTaskJetCore::GetNInputTracks FAILED");
1215 return -1;
1216 }
1217
1218 for (Int_t iJet=0; iJet<tmpAODjets->GetEntriesFast(); iJet++){
1219 AliAODJet *jet = dynamic_cast<AliAODJet*>((*tmpAODjets)[iJet]);
1220 if(!jet) continue;
1221 TRefArray *trackList = jet->GetRefTracks();
1222 Int_t nTracks = trackList->GetEntriesFast();
1223 nInputTracks += nTracks;
1224 if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks);
1225 }
1226 if(fDebug) Printf("---> input tracks: %d", nInputTracks);
1227
1228 return nInputTracks;
1229}
1230
1231
1232
ea693273 1233Double_t AliAnalysisTaskJetCore::RelativePhi(Double_t mphi,Double_t vphi){
1234
1235 if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
1236 else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
1237 if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
1238 else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
1239 double dphi = mphi-vphi;
1240 if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
1241 else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
1242 return dphi;//dphi in [-Pi, Pi]
1243}
1244
d5e7475c 1245Int_t AliAnalysisTaskJetCore::GetPhiBin(Double_t phi)
1246{
1247 Int_t phibin=-1;
1248 if(!(TMath::Abs(phi)<=2*TMath::Pi())){AliError("phi w.r.t. RP out of defined range");return -1;}
1249 Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi)));
1250 phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi()));
1251 if(phibin<0||phibin>=fNRPBins){AliError("Phi Bin not defined");}
1252 return phibin;
1253}
1254
1255
75bf77e3 1256
1257
da93bb11 1258THnSparse* AliAnalysisTaskJetCore::NewTHnSparseF(const char* name, UInt_t entries)
1259{
1260 // generate new THnSparseF, axes are defined in GetDimParams()
1261
1262 Int_t count = 0;
1263 UInt_t tmp = entries;
1264 while(tmp!=0){
1265 count++;
1266 tmp = tmp &~ -tmp; // clear lowest bit
1267 }
1268
1269 TString hnTitle(name);
1270 const Int_t dim = count;
1271 Int_t nbins[dim];
1272 Double_t xmin[dim];
1273 Double_t xmax[dim];
1274
1275 Int_t i=0;
1276 Int_t c=0;
1277 while(c<dim && i<32){
1278 if(entries&(1<<i)){
1279
1280 TString label("");
1281 GetDimParams(i, label, nbins[c], xmin[c], xmax[c]);
1282 hnTitle += Form(";%s",label.Data());
1283 c++;
1284 }
1285
1286 i++;
1287 }
1288 hnTitle += ";";
1289
1290 return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
1291}
1292
1293void AliAnalysisTaskJetCore::GetDimParams(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
1294{
1295 // stores label and binning of axis for THnSparse
1296
1297 const Double_t pi = TMath::Pi();
1298
1299 switch(iEntry){
1300
1301 case 0:
1302 label = "V0 centrality (%)";
1303
1304 nbins = 10;
1305 xmin = 0.;
1306 xmax = 100.;
1307 break;
1308
1309
1310 case 1:
1311 label = "corrected jet pt";
c42c6838 1312 nbins = 20;
da93bb11 1313 xmin = 0.;
1314 xmax = 200.;
1315 break;
1316
1317
1318 case 2:
1319 label = "track pT";
1320
ba143e7f 1321 nbins = 9;
da93bb11 1322 xmin = 0.;
5bd732d4 1323 xmax = 150;
da93bb11 1324 break;
1325
1326
529e2916 1327 case 3:
da93bb11 1328 label = "deltaR";
1329 nbins = 15;
1330 xmin = 0.;
1331 xmax = 1.5;
1332 break;
529e2916 1333
1334
1335
da93bb11 1336 case 4:
1337 label = "deltaEta";
ba143e7f 1338 nbins = 8;
1339 xmin = -1.6;
1340 xmax = 1.6;
da93bb11 1341 break;
1342
1343
5bd732d4 1344 case 5:
da93bb11 1345 label = "deltaPhi";
1346 nbins = 90;
1347 xmin = -0.5*pi;
1348 xmax = 1.5*pi;
1349 break;
1350
1351
529e2916 1352
1353 case 6:
da93bb11 1354 label = "leading track";
ba143e7f 1355 nbins = 13;
da93bb11 1356 xmin = 0;
20dcc500 1357 xmax = 50;
da93bb11 1358 break;
1359
529e2916 1360 case 7:
da93bb11 1361
1362 label = "trigger track";
ba143e7f 1363 nbins =10;
da93bb11 1364 xmin = 0;
1365 xmax = 50;
1366 break;
529e2916 1367
1368
1369
1370
1371
1372
1373
1374
da93bb11 1375 }
1376
1377}
1378