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