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