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