fzhou's changes
[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),
b9bc208b 70fDebug(0),
75bf77e3 71fOfflineTrgMask(AliVEvent::kAny),
72fMinContribVtx(1),
a9e585a7 73fVtxZMin(-10.),
74fVtxZMax(10.),
75bf77e3 75fEvtClassMin(0),
76fEvtClassMax(4),
8b47ec90 77fFilterMask(0),
41bcb7e7 78fFilterMaskBestPt(0),
d90d5d75 79fFilterType(0),
ea693273 80fRadioFrac(0.2),
81fMinDist(0.1),
75bf77e3 82fCentMin(0.),
83fCentMax(100.),
84fNInputTracksMin(0),
85fNInputTracksMax(-1),
ea693273 86fAngStructCloseTracks(0),
da93bb11 87fCheckMethods(0),
529e2916 88fDoEventMixing(0),
447d2a5b 89fFlagPhiBkg(0),
ca1009f1 90fFlagEtaBkg(0),
16a4c8e8 91fFlagJetHadron(0),
d5e7475c 92fFlagRandom(0),
ae1e07c1 93fFlagOnlyRecoil(0),
41bcb7e7 94fFlagOnlyHardest(1),
85b5b73e 95fTrackTypeRec(kTrackUndef),
d5e7475c 96fRPAngle(0),
41bcb7e7 97fNRPBins(50),
75bf77e3 98fJetEtaMin(-.5),
99fJetEtaMax(.5),
447d2a5b 100fNevents(0),
101fTindex(0),
102fTrigBufferIndex(0),
103fCountAgain(0),
75bf77e3 104fJetPtMin(20.),
105fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
106fJetPtFractionMin(0.5),
107fNMatchJets(4),
108fMatchMaxDist(0.8),
109fKeepJets(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),
c15cc27b 152fh3spectriggeredC10(0x0),
3353f803 153fh3spectriggeredC20(0x0),
d5e7475c 154fh3spectriggeredC3060(0x0)
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),
b9bc208b 188fDebug(0),
75bf77e3 189fOfflineTrgMask(AliVEvent::kAny),
190fMinContribVtx(1),
a9e585a7 191fVtxZMin(-10.),
192fVtxZMax(10.),
75bf77e3 193fEvtClassMin(0),
194fEvtClassMax(4),
8b47ec90 195fFilterMask(0),
41bcb7e7 196fFilterMaskBestPt(0),
d90d5d75 197fFilterType(0),
ea693273 198fRadioFrac(0.2),
199fMinDist(0.1),
75bf77e3 200fCentMin(0.),
201fCentMax(100.),
202fNInputTracksMin(0),
203fNInputTracksMax(-1),
ea693273 204fAngStructCloseTracks(0),
da93bb11 205fCheckMethods(0),
529e2916 206fDoEventMixing(0),
447d2a5b 207fFlagPhiBkg(0),
ca1009f1 208fFlagEtaBkg(0),
16a4c8e8 209fFlagJetHadron(0),
d5e7475c 210fFlagRandom(0),
ae1e07c1 211fFlagOnlyRecoil(0),
41bcb7e7 212fFlagOnlyHardest(1),
85b5b73e 213fTrackTypeRec(kTrackUndef),
d5e7475c 214fRPAngle(0),
41bcb7e7 215fNRPBins(50),
75bf77e3 216fJetEtaMin(-.5),
217fJetEtaMax(.5),
447d2a5b 218fNevents(0),
219fTindex(0),
220fTrigBufferIndex(0),
221fCountAgain(0),
75bf77e3 222fJetPtMin(20.),
223fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
224fJetPtFractionMin(0.5),
225fNMatchJets(4),
226fMatchMaxDist(0.8),
227fKeepJets(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),
c15cc27b 270fh3spectriggeredC10(0x0),
3353f803 271fh3spectriggeredC20(0x0),
d5e7475c 272fh3spectriggeredC3060(0x0)
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.);
b9bc208b 411 fh2RPJetsC10=new TH2F("RPJetC10","",35,0.,3.5,100,0.,100.);
412 fh2RPJetsC20=new TH2F("RPJetC20","",35,0.,3.5,100,0.,100.);
413 fh2RPTC10=new TH2F("RPTriggerC10","",35,0.,3.5,50,0.,50.);
414 fh2RPTC20=new TH2F("RPTriggerC20","",35,0.,3.5,50,0.,50.);
d90d5d75 415 fh3spectriggeredC10 = new TH3F("Triggered spectrumC10","",100,0.,1.,140,-80.,200.,50,0.,50.);
416 fh3spectriggeredC20 = new TH3F("Triggered spectrumC20","",100,0.,1.,140,-80.,200.,50,0.,50.);
ae1e07c1 417 fh3spectriggeredC3060 = new TH3F("Triggered spectrumC3060","",100,0.,1.,140,-80.,200.,10,0.,50.);
3353f803 418
d5e7475c 419
420
75bf77e3 421 fOutputList->Add(fHistEvtSelection);
a9e585a7 422
da93bb11 423 fOutputList->Add(fhnDeltaR);
424
529e2916 425 fOutputList->Add(fhnMixedEvents);
da93bb11 426
427
428
429 if(fCheckMethods){
529e2916 430
75bf77e3 431 fOutputList->Add(fh2JetCoreMethod1C10);
432 fOutputList->Add(fh2JetCoreMethod2C10);
ea693273 433 fOutputList->Add(fh2JetCoreMethod1C20);
434 fOutputList->Add(fh2JetCoreMethod2C20);
75bf77e3 435 fOutputList->Add(fh2JetCoreMethod1C30);
436 fOutputList->Add(fh2JetCoreMethod2C30);
75bf77e3 437 fOutputList->Add(fh2JetCoreMethod1C60);
da93bb11 438 fOutputList->Add(fh2JetCoreMethod2C60);}
a9e585a7 439
d5e7475c 440 fOutputList->Add(fh3JetTrackC3060);
e03c221a 441 fOutputList->Add(fh3JetTrackC20);
85b5b73e 442
4e90d948 443
75bf77e3 444
a9e585a7 445
446
447 if(fAngStructCloseTracks>0){
ea693273 448 fOutputList->Add(fh2AngStructpt1C10);
449 fOutputList->Add(fh2AngStructpt2C10);
450 fOutputList->Add(fh2AngStructpt3C10);
451 fOutputList->Add(fh2AngStructpt4C10);
452 fOutputList->Add(fh2AngStructpt1C20);
453 fOutputList->Add(fh2AngStructpt2C20);
454 fOutputList->Add(fh2AngStructpt3C20);
455 fOutputList->Add(fh2AngStructpt4C20);
456 fOutputList->Add(fh2AngStructpt1C30);
457 fOutputList->Add(fh2AngStructpt2C30);
458 fOutputList->Add(fh2AngStructpt3C30);
459 fOutputList->Add(fh2AngStructpt4C30);
460 fOutputList->Add(fh2AngStructpt1C60);
461 fOutputList->Add(fh2AngStructpt2C60);
462 fOutputList->Add(fh2AngStructpt3C60);
a9e585a7 463 fOutputList->Add(fh2AngStructpt4C60);}
20dcc500 464
465
d5e7475c 466
5bd732d4 467
ba143e7f 468
8205e054 469 fOutputList->Add(fh2Ntriggers);
d90d5d75 470 fOutputList->Add(fh2Ntriggers2C10);
471 fOutputList->Add(fh2Ntriggers2C20);
ae1e07c1 472 fOutputList->Add(fh3JetDensity);
473 fOutputList->Add(fh3JetDensityA4);
d90d5d75 474 fOutputList->Add(fh2RPJetsC10);
475 fOutputList->Add(fh2RPJetsC20);
476 fOutputList->Add(fh2RPTC10);
477 fOutputList->Add(fh2RPTC20);
c15cc27b 478 fOutputList->Add(fh3spectriggeredC10);
41bcb7e7 479 fOutputList->Add(fh3spectriggeredC20);
480 fOutputList->Add(fh3spectriggeredC3060);
3353f803 481
d5e7475c 482
75bf77e3 483 // =========== Switch on Sumw2 for all histos ===========
484 for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
485 TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
486 if (h1){
487 h1->Sumw2();
488 continue;
489 }
ea693273 490 THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
491 if (hn){
492 hn->Sumw2();
493 }
75bf77e3 494 }
495 TH1::AddDirectory(oldStatus);
496
497 PostData(1, fOutputList);
498}
499
500void AliAnalysisTaskJetCore::UserExec(Option_t *)
501{
502
503
504 if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
505 AliError("Jet branch name not set.");
506 return;
507 }
508
509 fESD=dynamic_cast<AliESDEvent*>(InputEvent());
510 if (!fESD) {
511 AliError("ESD not available");
5bd732d4 512 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
513 }
514 fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
515
516 static AliAODEvent* aod = 0;
517 // take all other information from the aod we take the tracks from
518 if(!aod){
519 if(!fESD)aod = fAODIn;
520 else aod = fAODOut;}
521
522
ea693273 523
524 if(fNonStdFile.Length()!=0){
525 // case that we have an AOD extension we need can fetch the jets from the extended output
526 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
527 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
528 if(!fAODExtension){
529 if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
530 }}
531
532
75bf77e3 533 // -- event selection --
534 fHistEvtSelection->Fill(1); // number of events before event selection
535
536 // physics selection
537 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
538 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
0a1ffd97 539 std::cout<<inputHandler->IsEventSelected()<<" "<<fOfflineTrgMask<<std::endl;
75bf77e3 540 if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
541 if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
542 fHistEvtSelection->Fill(2);
543 PostData(1, fOutputList);
544 return;
545 }
546
547 // vertex selection
5bd732d4 548 if(!aod){
52721bee 549 if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
550 fHistEvtSelection->Fill(3);
551 PostData(1, fOutputList);
552 }
5bd732d4 553 AliAODVertex* primVtx = aod->GetPrimaryVertex();
52721bee 554
555 if(!primVtx){
556 if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
557 fHistEvtSelection->Fill(3);
d27719bf 558 PostData(1, fOutputList);
559 return;
52721bee 560 }
561
75bf77e3 562 Int_t nTracksPrim = primVtx->GetNContributors();
563 if ((nTracksPrim < fMinContribVtx) ||
564 (primVtx->GetZ() < fVtxZMin) ||
565 (primVtx->GetZ() > fVtxZMax) ){
566 if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
567 fHistEvtSelection->Fill(3);
568 PostData(1, fOutputList);
569 return;
570 }
571
572 // event class selection (from jet helper task)
573 Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
574 if(fDebug) Printf("Event class %d", eventClass);
575 if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
576 fHistEvtSelection->Fill(4);
577 PostData(1, fOutputList);
578 return;
579 }
580
581 // centrality selection
582 AliCentrality *cent = 0x0;
da93bb11 583 Double_t centValue = 0.;
8205e054 584 if(fIsPbPb){
20dcc500 585 if(fESD) {cent = fESD->GetCentrality();
586 if(cent) centValue = cent->GetCentralityPercentile("V0M");}
5bd732d4 587 else centValue=aod->GetHeader()->GetCentrality();
20dcc500 588
75bf77e3 589 if(fDebug) printf("centrality: %f\n", centValue);
3353f803 590 if (centValue < fCentMin || centValue > fCentMax){
591 fHistEvtSelection->Fill(4);
592 PostData(1, fOutputList);
593 return;
8205e054 594 }}
75bf77e3 595
596
568f8fa2 597 fHistEvtSelection->Fill(0);
598 // accepted events
75bf77e3 599 // -- end event selection --
ea693273 600
75bf77e3 601 // get background
602 AliAODJetEventBackground* externalBackground = 0;
5bd732d4 603 if(fAODOut&&!externalBackground&&fBackgroundBranch.Length()){
604 externalBackground = (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch.Data()));
75bf77e3 605 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
606 }
ea693273 607 if(fAODExtension&&!externalBackground&&fBackgroundBranch.Length()){
608 externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
609 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
610 }
5bd732d4 611
612 if(fAODIn&&!externalBackground&&fBackgroundBranch.Length()){
613 externalBackground = (AliAODJetEventBackground*)(fAODIn->FindListObject(fBackgroundBranch.Data()));
614 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
615 }
ea693273 616
75bf77e3 617 Float_t rho = 0;
75bf77e3 618
d5e7475c 619 if(fFlagRandom==0){
620 if(externalBackground)rho = externalBackground->GetBackground(0);}
621 if(fFlagRandom==1){
622 if(externalBackground)rho = externalBackground->GetBackground(2);}
75bf77e3 623
624 // fetch jets
625 TClonesArray *aodJets[2];
ea693273 626 aodJets[0]=0;
5bd732d4 627 if(fAODOut&&!aodJets[0]){
628 aodJets[0] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[0].Data()));
629 aodJets[1] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[1].Data())); }
ea693273 630 if(fAODExtension && !aodJets[0]){
529e2916 631 aodJets[0] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[0].Data()));
632 aodJets[1] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[1].Data())); }
5bd732d4 633 if(fAODIn&&!aodJets[0]){
634 aodJets[0] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[0].Data()));
635 aodJets[1] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[1].Data())); }
636
ea693273 637
529e2916 638 //Double_t ptsub[aodJets[0]->GetEntriesFast()];
639 //Int_t inord[aodJets[0]->GetEntriesFast()];
640 //for(Int_t n=0;n<aodJets[0]->GetEntriesFast();n++){
641 // ptsub[n]=0;
642 // inord[n]=0;}
a9e585a7 643
ea693273 644 TList ParticleList;
645 Int_t nT = GetListOfTracks(&ParticleList);
646 for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
75bf77e3 647 fListJets[iJetType]->Clear();
648 if (!aodJets[iJetType]) continue;
75bf77e3 649 if(fDebug) Printf("%s: %d jets",fJetBranchName[iJetType].Data(),aodJets[iJetType]->GetEntriesFast());
75bf77e3 650 for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
651 AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
529e2916 652 if (jet) fListJets[iJetType]->Add(jet);
653 // if(iJetType==0){
654 // ptsub[iJet]=jet->Pt()-rho*jet->EffectiveAreaCharged();}
655 }}
75bf77e3 656
657 Double_t etabig=0;
658 Double_t ptbig=0;
659 Double_t areabig=0;
660 Double_t phibig=0.;
661 Double_t etasmall=0;
662 Double_t ptsmall=0;
663 Double_t areasmall=0;
75bf77e3 664 Double_t phismall=0.;
a9e585a7 665
da93bb11 666
d5e7475c 667
529e2916 668 Int_t iCount=0;
669 Int_t trigJet=-1;
670 Int_t trigBBTrack=-1;
671 Int_t trigInTrack=-1;
d5e7475c 672 fRPAngle = aod->GetHeader()->GetEventplane();
85b5b73e 673
41bcb7e7 674 AliVParticle *partback = (AliVParticle*)ParticleList.At(nT);
675 if(!partback){
676 PostData(1, fOutputList);
677 return;}
e22029ea 678
679
41bcb7e7 680 //for(Int_t tt=0;tt<ParticleList.GetEntries();tt++){
681 //if(fFlagOnlyHardest!=0){if(tt!=nT) continue;}
682 //AliVParticle *partback = (AliVParticle*)ParticleList.At(tt);
e22029ea 683
8205e054 684 fh2Ntriggers->Fill(centValue,partback->Pt());
d90d5d75 685 Double_t phiBinT = RelativePhi(partback->Phi(),fRPAngle);
686 if(centValue<20.) fh2RPTC20->Fill(TMath::Abs(phiBinT),partback->Pt());
687 if(centValue<10.) fh2RPTC10->Fill(TMath::Abs(phiBinT),partback->Pt());
8205e054 688 Double_t accep=2.*TMath::Pi()*1.8;
4e90d948 689 Int_t injet4=0;
690 Int_t injet=0;
41bcb7e7 691
75bf77e3 692 for(Int_t i=0; i<fListJets[0]->GetEntries(); ++i){
693 AliAODJet* jetbig = (AliAODJet*)(fListJets[0]->At(i));
694 etabig = jetbig->Eta();
695 phibig = jetbig->Phi();
696 ptbig = jetbig->Pt();
697 if(ptbig==0) continue;
d90d5d75 698 Double_t phiBin = RelativePhi(phibig,fRPAngle);
75bf77e3 699 areabig = jetbig->EffectiveAreaCharged();
ea693273 700 Double_t ptcorr=ptbig-rho*areabig;
75bf77e3 701 if((etabig<fJetEtaMin)||(etabig>fJetEtaMax)) continue;
ae1e07c1 702 if(areabig>=0.07) injet=injet+1;
8205e054 703 if(areabig>=0.4) injet4=injet4+1;
ca1009f1 704 Double_t dphi=RelativePhi(partback->Phi(),phibig);
705
ae1e07c1 706 if(fFlagEtaBkg==1){
ca1009f1 707 Double_t etadif= partback->Eta()-etabig;
708 if(TMath::Abs(etadif)<=0.5){
85b5b73e 709
ca1009f1 710 if(centValue<20.) fh3JetTrackC20->Fill(partback->Pt(),ptcorr,TMath::Abs(dphi));
711 if(centValue>30. && centValue<60.) fh3JetTrackC3060->Fill(partback->Pt(),ptcorr,TMath::Abs(dphi));}}
ca1009f1 712 if(fFlagEtaBkg==0){
e03c221a 713 if(centValue<20.) fh3JetTrackC20->Fill(partback->Pt(),ptcorr,TMath::Abs(dphi));
ca1009f1 714 if(centValue>30. && centValue<60.) fh3JetTrackC3060->Fill(partback->Pt(),ptcorr,TMath::Abs(dphi));}
715
716
16a4c8e8 717 if(fFlagJetHadron==0){
447d2a5b 718 if(fFlagPhiBkg!=0) if((TMath::Abs(dphi)<TMath::Pi()/2.-0.1)||(TMath::Abs(dphi)>TMath::Pi()/2.+0.1)) continue;
16a4c8e8 719 if(fFlagPhiBkg==0) if(TMath::Abs(dphi)<TMath::Pi()-0.6) continue;}
720
721 if(fFlagJetHadron!=0) if(TMath::Abs(dphi)>0.4) continue;
722
723
d90d5d75 724 if(centValue<10.) fh2RPJetsC10->Fill(TMath::Abs(phiBin), ptcorr);
725 if(centValue<20.) fh2RPJetsC20->Fill(TMath::Abs(phiBin), ptcorr);
75bf77e3 726 Double_t dismin=100.;
727 Double_t ptmax=-10.;
728 Int_t index1=-1;
729 Int_t index2=-1;
e22029ea 730
c15cc27b 731 if(centValue<10.) fh3spectriggeredC10->Fill(jetbig->EffectiveAreaCharged(),ptcorr,partback->Pt());
3353f803 732 if(centValue<20.) fh3spectriggeredC20->Fill(jetbig->EffectiveAreaCharged(),ptcorr,partback->Pt());
733 if(centValue>30. && centValue<60.) fh3spectriggeredC3060->Fill(jetbig->EffectiveAreaCharged(),ptcorr,partback->Pt());
734
735 if(ptcorr<=0) continue;
16a4c8e8 736
16a4c8e8 737 AliAODTrack* leadtrack=0;
da93bb11 738 Int_t ippt=0;
16a4c8e8 739 Double_t ppt=-10;
740 if(fFlagJetHadron==0){
529e2916 741 TRefArray *genTrackList = jetbig->GetRefTracks();
742 Int_t nTracksGenJet = genTrackList->GetEntriesFast();
743 AliAODTrack* genTrack;
da93bb11 744 for(Int_t ir=0; ir<nTracksGenJet; ++ir){
745 genTrack = (AliAODTrack*)(genTrackList->At(ir));
746 if(genTrack->Pt()>ppt){ppt=genTrack->Pt();
16a4c8e8 747 ippt=ir;}}
748 leadtrack=(AliAODTrack*)(genTrackList->At(ippt));
749 if(!leadtrack) continue;}
750
751 AliVParticle* leadtrackb=0;
752 if(fFlagJetHadron!=0){
753 Int_t nTb = GetHardestTrackBackToJet(jetbig);
754 leadtrackb = (AliVParticle*)ParticleList.At(nTb);
755 if(!leadtrackb) continue;
756 }
757
758
759
760
d5e7475c 761
529e2916 762 //store one trigger info
447d2a5b 763 if(iCount==0){
529e2916 764 trigJet=i;
8205e054 765 trigBBTrack=nT;
529e2916 766 trigInTrack=ippt;
767 iCount=iCount+1;}
768
d5e7475c 769
529e2916 770 if(fCheckMethods){
771 for(Int_t j=0; j<fListJets[1]->GetEntries(); ++j){
75bf77e3 772 AliAODJet* jetsmall = (AliAODJet*)(fListJets[1]->At(j));
773 etasmall = jetsmall->Eta();
774 phismall = jetsmall->Phi();
775 ptsmall = jetsmall->Pt();
776 areasmall = jetsmall->EffectiveAreaCharged();
ea693273 777 Double_t tmpDeltaR=(phismall-phibig)*(phismall-phibig)+(etasmall-etabig)*(etasmall-etabig);
778 tmpDeltaR=TMath::Sqrt(tmpDeltaR);
779 //Fraction in the jet core
780 if((ptsmall>ptmax)&&(tmpDeltaR<=fRadioFrac)){ptmax=ptsmall;
75bf77e3 781 index2=j;}
ea693273 782 if(tmpDeltaR<=dismin){ dismin=tmpDeltaR;
783 index1=j;}} //en of loop over R=0.2 jets
75bf77e3 784 //method1:most concentric jet=core
785 if(dismin<fMinDist){ AliAODJet* jetmethod1 = (AliAODJet*)(fListJets[1]->At(index1));
ea693273 786 if(centValue<10) fh2JetCoreMethod1C10->Fill(ptcorr,jetmethod1->Pt()/ptbig);
787 if((centValue>20)&&(centValue<40)) fh2JetCoreMethod1C20->Fill(ptcorr,jetmethod1->Pt()/ptbig);
788 if((centValue>30)&&(centValue<60)) fh2JetCoreMethod1C30->Fill(ptcorr,jetmethod1->Pt()/ptbig);
789 if(centValue>60) fh2JetCoreMethod1C60->Fill(ptcorr,jetmethod1->Pt()/ptbig); }
75bf77e3 790 //method2:hardest contained jet=core
791 if(index2!=-1){
792 AliAODJet* jetmethod2 = (AliAODJet*)(fListJets[1]->At(index2));
ea693273 793 if(centValue<10) fh2JetCoreMethod2C10->Fill(ptcorr,jetmethod2->Pt()/ptbig);
794 if((centValue>20)&&(centValue<40)) fh2JetCoreMethod2C20->Fill(ptcorr,jetmethod2->Pt()/ptbig);
795 if((centValue>30)&&(centValue<60)) fh2JetCoreMethod2C30->Fill(ptcorr,jetmethod2->Pt()/ptbig);
da93bb11 796 if(centValue>60) fh2JetCoreMethod2C60->Fill(ptcorr,jetmethod2->Pt()/ptbig); }}
d90d5d75 797 if(centValue<10) fh2Ntriggers2C10->Fill(leadtrack->Pt(),partback->Pt());
798 if(centValue<20) fh2Ntriggers2C20->Fill(leadtrack->Pt(),partback->Pt());
799 if(fDoEventMixing==0 && fFlagOnlyRecoil==0){
41bcb7e7 800 for(int it = 0;it<ParticleList.GetEntries();++it){
801 AliVParticle *part = (AliVParticle*)ParticleList.At(it);
802 Double_t deltaR = jetbig->DeltaR(part);
803 Double_t deltaEta = etabig-part->Eta();
d90d5d75 804
41bcb7e7 805 Double_t deltaPhi=phibig-part->Phi();
806 if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
807 if(deltaPhi>3./2.*TMath::Pi()) deltaPhi-=2.*TMath::Pi();
808 Double_t pTcont=0;
809 if(fFlagJetHadron==0) pTcont=leadtrack->Pt();
810 if(fFlagJetHadron!=0) pTcont=leadtrackb->Pt();
811 Double_t jetEntries[8] = {centValue,ptcorr,part->Pt(),deltaR,deltaEta,deltaPhi,pTcont,partback->Pt()};
812 fhnDeltaR->Fill(jetEntries);}
813
814
d90d5d75 815 }
85b5b73e 816
d5e7475c 817
d5e7475c 818 //end of track loop, we only do it if EM is switched off
819
ba143e7f 820
821
822
823
824
825
826
20dcc500 827
828 }
ae1e07c1 829 if(injet>0) fh3JetDensity->Fill(ParticleList.GetEntries(),injet/accep,partback->Pt());
830 if(injet4>0)fh3JetDensityA4->Fill(ParticleList.GetEntries(),injet4/accep,partback->Pt());
529e2916 831 //end of jet loop
832
41bcb7e7 833 //}
529e2916 834
835
447d2a5b 836 if(fDoEventMixing>0){
529e2916 837 //check before if the trigger exists
838 // fTrigBuffer[i][0] = zvtx
839 // fTrigBuffer[i][1] = phi
840 // fTrigBuffer[i][2] = eta
841 // fTrigBuffer[i][3] = pt_jet
842 // fTrigBuffer[i][4] = pt_trig
447d2a5b 843 // fTrigBuffer[i][5]= centrality
844 if(fTindex==10) fTindex=0;
529e2916 845 if(fTrigBuffer[fTindex][3]>0){
846 if (TMath::Abs(fTrigBuffer[fTindex][0]-primVtx->GetZ()<2.)){
447d2a5b 847 if (TMath::Abs(fTrigBuffer[fTindex][5]-centValue<5)){
529e2916 848
849 for(int it = 0;it<nT;++it){
20dcc500 850 AliVParticle *part = (AliVParticle*)ParticleList.At(it);
851 Double_t DPhi = fTrigBuffer[fTindex][1] - part->Phi();
529e2916 852 Double_t DEta = fTrigBuffer[fTindex][2] - part->Eta();
853 Double_t DR=TMath::Sqrt(DPhi*DPhi+DEta*DEta);
854 if(DPhi<-0.5*TMath::Pi()) DPhi+=2.*TMath::Pi();
855 if(DPhi>3./2.*TMath::Pi()) DPhi-=2.*TMath::Pi();
447d2a5b 856 Double_t triggerEntries[7] = {centValue,fTrigBuffer[fTindex][3],part->Pt(),DR,DEta,DPhi,fTrigBuffer[fTindex][4]};
529e2916 857 fhnMixedEvents->Fill(triggerEntries);
858 }
859 fNevents=fNevents+1;
447d2a5b 860 if(fNevents==10) fTindex=fTindex+1;
529e2916 861 }}}
ae1e07c1 862
447d2a5b 863 if(fTindex==10&&fNevents==10) fCountAgain=0;
529e2916 864
865 // Copy the triggers from the current event into the buffer.
866 //again, only if the trigger exists:
447d2a5b 867 if(fCountAgain==0){
529e2916 868 if(trigJet>-1){
447d2a5b 869 AliAODJet* jetT = (AliAODJet*)(fListJets[0]->At(trigJet)); AliVParticle *partT = (AliVParticle*)ParticleList.At(trigBBTrack);
529e2916 870 fTrigBuffer[fTrigBufferIndex][0] = primVtx->GetZ();
871 fTrigBuffer[fTrigBufferIndex][1] = jetT->Phi();
872 fTrigBuffer[fTrigBufferIndex][2] = jetT->Eta();
873 fTrigBuffer[fTrigBufferIndex][3] = jetT->Pt()-rho*jetT->EffectiveAreaCharged();
874 fTrigBuffer[fTrigBufferIndex][4] = partT->Pt();
447d2a5b 875 fTrigBuffer[fTrigBufferIndex][5] = centValue;
529e2916 876 fTrigBufferIndex++;
447d2a5b 877 if(fTrigBufferIndex==9) {fTrigBufferIndex=0;
878 fCountAgain=1;}
529e2916 879 }
447d2a5b 880 }
529e2916 881
447d2a5b 882 }
529e2916 883
884
da93bb11 885
ea693273 886 //////////////////ANGULAR STRUCTURE//////////////////////////////////////
887
888 //tracks up to R=0.8 distant from the jet axis
529e2916 889 // if(fAngStructCloseTracks==1){
890 // TList CloseTrackList;
891 // Int_t nn=GetListOfTracksCloseToJet(&CloseTrackList,jetbig);
892 // Double_t difR=0.04;
893 // for(Int_t l=0;l<15;l++){
894 // Double_t rr=l*0.1+0.1;
895 // for(int it = 0;it<nn;++it){
896 // AliVParticle *part1 = (AliVParticle*)CloseTrackList.At(it);
897 // for(int itu=it+1;itu<CloseTrackList.GetEntries();itu++){
898 // AliVParticle *part2 = (AliVParticle*)CloseTrackList.At(itu);
899 // Double_t ptm=part1->Pt();
900 // Double_t ptn=part2->Pt();
901 // Double_t Rnm = (part1->Eta()-part2->Eta())*(part1->Eta()-part2->Eta())+(part1->Phi()-part2->Phi())*(part1->Phi()-part2->Phi());
902 // Rnm=TMath::Sqrt(Rnm);
903 // Double_t deltag=(1./(TMath::Sqrt(2*TMath::Pi())*difR))*TMath::Exp(-1.*(rr-Rnm)*(rr-Rnm)/(2.*difR*difR));
904 // Double_t stepf=0.5*(1.+TMath::Erf((rr-Rnm)/(TMath::Sqrt(2.)*difR)));
905 // if((ptcorr<85.) && (ptcorr>=70.)){up1[l]=up1[l]+ptm*ptn*Rnm*Rnm*deltag;
906 // down1[l]=down1[l]+ptm*ptn*Rnm*Rnm*stepf;}
907 // if((ptcorr<100.) && (ptcorr>=85.)){up2[l]=up2[l]+ptm*ptn*Rnm*Rnm*deltag;
908 // down2[l]=down2[l]+ptm*ptn*Rnm*Rnm*stepf;}
909 // if((ptcorr<120.) && (ptcorr>=100.)){up3[l]=up3[l]+ptm*ptn*Rnm*Rnm*deltag;
910 // down3[l]=down3[l]+ptm*ptn*Rnm*Rnm*stepf;}
911 // if((ptcorr<140.) && (ptcorr>=120.)){up4[l]=up4[l]+ptm*ptn*Rnm*Rnm*deltag;
912 // down4[l]=down4[l]+ptm*ptn*Rnm*Rnm*stepf;}}}}
913 // }
ea693273 914
529e2916 915 // //only jet constituents
916 // if(fAngStructCloseTracks==2){
ea693273 917
529e2916 918 // Double_t difR=0.04;
919 // for(Int_t l=0;l<15;l++){
920 // Double_t rr=l*0.1+0.1;
ea693273 921
922
529e2916 923 // AliAODTrack* part1;
924 // AliAODTrack* part2;
da93bb11 925
529e2916 926 // TRefArray *genTrackListb = jetbig->GetRefTracks();
927 // Int_t nTracksGenJetb = genTrackListb->GetEntriesFast();
da93bb11 928
929
930
529e2916 931 // for(Int_t it=0; it<nTracksGenJetb; ++it){
932 // part1 = (AliAODTrack*)(genTrackListb->At(it));
933 // for(Int_t itu=0; itu<nTracksGenJetb; ++itu){
934 // part2 = (AliAODTrack*)(genTrackListb->At(itu));
935 // Double_t ptm=part1->Pt();
936 // Double_t ptn=part2->Pt();
937 // Double_t Rnm = (part1->Eta()-part2->Eta())*(part1->Eta()-part2->Eta())+(part1->Phi()-part2->Phi())*(part1->Phi()-part2->Phi());
938 // Rnm=TMath::Sqrt(Rnm);
939 // Double_t deltag=(1./(TMath::Sqrt(2*TMath::Pi())*difR))*TMath::Exp(-1.*(rr-Rnm)*(rr-Rnm)/(2.*difR*difR));
940 // Double_t stepf=0.5*(1.+TMath::Erf((rr-Rnm)/(TMath::Sqrt(2.)*difR)));
941 // if((ptcorr<85.) && (ptcorr>=70.)){up1[l]=up1[l]+ptm*ptn*Rnm*Rnm*deltag;
942 // down1[l]=down1[l]+ptm*ptn*Rnm*Rnm*stepf;}
943 // if((ptcorr<100.) && (ptcorr>=85.)){up2[l]=up2[l]+ptm*ptn*Rnm*Rnm*deltag;
944 // down2[l]=down2[l]+ptm*ptn*Rnm*Rnm*stepf;}
945 // if((ptcorr<120.) && (ptcorr>=100.)){up3[l]=up3[l]+ptm*ptn*Rnm*Rnm*deltag;
946 // down3[l]=down3[l]+ptm*ptn*Rnm*Rnm*stepf;}
947 // if((ptcorr<140.) && (ptcorr>=120.)){up4[l]=up4[l]+ptm*ptn*Rnm*Rnm*deltag;
948 // down4[l]=down4[l]+ptm*ptn*Rnm*Rnm*stepf;}}}}}
949 // }
950 // //end loop over R=0.4 jets
951 // if(fAngStructCloseTracks>0){
952 // for(Int_t l=0;l<15;l++){
953 // Double_t rr=l*0.1+0.1;
954 // if(down1[l]!=0){
955 // if(centValue<10.)fh2AngStructpt1C10->Fill(rr,rr*up1[l]/down1[l]);
956 // if(centValue>20. && centValue<40.) fh2AngStructpt1C20->Fill(rr,rr*up1[l]/down1[l]);
957 // if(centValue>30. && centValue<60.) fh2AngStructpt1C30->Fill(rr,rr*up1[l]/down1[l]);
958 // if(centValue>60.) fh2AngStructpt1C60->Fill(rr,rr*up1[l]/down1[l]);}
959 // if(down2[l]!=0){
960 // if(centValue<10.) fh2AngStructpt2C10->Fill(rr,rr*up2[l]/down2[l]);
961 // if(centValue>20. && centValue<40.) fh2AngStructpt2C20->Fill(rr,rr*up2[l]/down2[l]);
962 // if(centValue>30. && centValue<60.) fh2AngStructpt2C30->Fill(rr,rr*up2[l]/down2[l]);
963 // if(centValue>60.) fh2AngStructpt2C60->Fill(rr,rr*up2[l]/down2[l]);}
964 // if(down3[l]!=0){
965 // if(centValue<10.) fh2AngStructpt3C10->Fill(rr,rr*up3[l]/down3[l]);
966 // if(centValue>20. && centValue<40.) fh2AngStructpt3C20->Fill(rr,rr*up3[l]/down3[l]);
967 // if(centValue>30. && centValue<60.) fh2AngStructpt3C30->Fill(rr,rr*up3[l]/down3[l]);
968 // if(centValue>60.) fh2AngStructpt3C60->Fill(rr,rr*up3[l]/down3[l]);}
969 // if(down4[l]!=0){
970 // if(centValue<10.) fh2AngStructpt4C10->Fill(rr,rr*up4[l]/down4[l]);
971 // if(centValue>20. && centValue<40.) fh2AngStructpt4C20->Fill(rr,rr*up4[l]/down4[l]);
972 // if(centValue>30. && centValue<60.) fh2AngStructpt4C30->Fill(rr,rr*up4[l]/down4[l]);
973 // if(centValue>60.) fh2AngStructpt4C60->Fill(rr,rr*up4[l]/down4[l]);}}}
ea693273 974
975
976
977
978
75bf77e3 979
980
981 PostData(1, fOutputList);
da93bb11 982}
75bf77e3 983
984void AliAnalysisTaskJetCore::Terminate(const Option_t *)
985{
986 // Draw result to the screen
987 // Called once at the end of the query
988
989 if (!GetOutputData(1))
990 return;
991}
992
ea693273 993
994
5bd732d4 995
ea693273 996
997
998Int_t AliAnalysisTaskJetCore::GetListOfTracks(TList *list){
999
85b5b73e 1000 Int_t iCount = 0;
1001 AliAODEvent *aod = 0;
1002
ba143e7f 1003 if(!fESD)aod = fAODIn;
5bd732d4 1004 else aod = fAODOut;
85b5b73e 1005 Int_t index=-1;
8205e054 1006 Double_t ptmax=-10;
85b5b73e 1007
1008
1009
1010 for(int it = 0;it < aod->GetNumberOfTracks();++it){
5bd732d4 1011 AliAODTrack *tr = aod->GetTrack(it);
ae73f4c9 1012 Bool_t bGood = false;
1013 if(fFilterType == 0)bGood = true;
1014 else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
d90d5d75 1015 else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
1016 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
ae73f4c9 1017 if(bGood==false) continue;
ea693273 1018 if(TMath::Abs(tr->Eta())>0.9)continue;
1019 if(tr->Pt()<0.15)continue;
1020 list->Add(tr);
ea693273 1021 iCount++;
d90d5d75 1022 if(fFilterType==2 && fFilterMaskBestPt>0){// only set the trigger track index for good quality tracks
873306cf 1023 if(tr->TestFilterBit(fFilterMaskBestPt)){
1024 if(tr->Pt()>ptmax){
1025 ptmax=tr->Pt();
1026 index=iCount-1;
1027 }
1028 }
1029 }
1030 else{
1031 if(tr->Pt()>ptmax){
1032 ptmax=tr->Pt();
1033 index=iCount-1;
1034 }
1035 }
85b5b73e 1036 }
ea693273 1037
529e2916 1038
85b5b73e 1039 // else if (type == kTrackAODMCCharged) {
1040 // TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1041 // if(!tca)return iCount;
1042 // for(int it = 0;it < tca->GetEntriesFast();++it){
1043 // AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1044 // if(!part)continue;
1045 // if(part->Pt()<0.15)continue;
1046 // if(!part->IsPhysicalPrimary())continue;
1047 // if(part->Charge()==0)continue;
1048 // if(TMath::Abs(part->Eta())>0.9)continue;
1049 // list->Add(part);
1050 // iCount++;
1051 // if(part->Pt()>ptmax){ ptmax=part->Pt();
1052 // index=iCount-1;}}}
1053 return index;
ea693273 1054
1055}
1056
da93bb11 1057 Int_t AliAnalysisTaskJetCore::GetHardestTrackBackToJet(AliAODJet *jetbig){
5bd732d4 1058
1059 AliAODEvent *aod = 0;
1060 if(!fESD)aod = fAODIn;
1061 else aod = fAODOut;
da93bb11 1062 Int_t index=-1;
1063 Double_t ptmax=-10;
1064 Double_t dphi=0;
1065 Double_t dif=0;
1066 Int_t iCount=0;
5bd732d4 1067 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1068 AliAODTrack *tr = aod->GetTrack(it);
da93bb11 1069 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1070 if(TMath::Abs(tr->Eta())>0.9)continue;
1071 if(tr->Pt()<0.15)continue;
1072 iCount=iCount+1;
529e2916 1073 dphi=RelativePhi(tr->Phi(),jetbig->Phi());
16a4c8e8 1074 if(TMath::Abs(dphi)<TMath::Pi()-0.6) continue;
da93bb11 1075 if(tr->Pt()>ptmax){ ptmax=tr->Pt();
b6844e9e 1076 index=iCount-1;
1077 dif=dphi; }}
da93bb11 1078
1079 return index;
1080
1081 }
1082
1083
1084
1085
1086
1087
1088
1089
1090
ea693273 1091 Int_t AliAnalysisTaskJetCore::GetListOfTracksCloseToJet(TList *list,AliAODJet *jetbig){
1092
1093 Int_t iCount = 0;
5bd732d4 1094 AliAODEvent *aod = 0;
1095 if(!fESD)aod = fAODIn;
1096 else aod = fAODOut;
8b47ec90 1097
16a4c8e8 1098 for(int it = 0;it < aod->GetNumberOfTracks();++it){
5bd732d4 1099 AliAODTrack *tr = aod->GetTrack(it);
ea693273 1100 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1101 if(TMath::Abs(tr->Eta())>0.9)continue;
1102 if(tr->Pt()<0.15)continue;
1103 Double_t disR=jetbig->DeltaR(tr);
1104 if(disR>0.8) continue;
1105 list->Add(tr);
1106 //cout<<fAOD->GetNumberOfTracks()<<" "<<tr->Pt()<<endl;
1107 iCount++;
1108 }
1109
1110 list->Sort();
1111 return iCount;
1112
1113}
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
75bf77e3 1125Int_t AliAnalysisTaskJetCore::GetNInputTracks()
1126{
1127
1128 Int_t nInputTracks = 0;
5bd732d4 1129 AliAODEvent *aod = 0;
1130 if(!fESD)aod = fAODIn;
1131 else aod = fAODOut;
75bf77e3 1132 TString jbname(fJetBranchName[1]);
1133 //needs complete event, use jets without background subtraction
1134 for(Int_t i=1; i<=3; ++i){
1135 if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0");
1136 }
1137 // use only HI event
1138 if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
1139 if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
1140
1141 if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
5bd732d4 1142 TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(aod->FindListObject(jbname.Data()));
75bf77e3 1143 if(!tmpAODjets){
1144 Printf("Jet branch %s not found", jbname.Data());
1145 Printf("AliAnalysisTaskJetCore::GetNInputTracks FAILED");
1146 return -1;
1147 }
1148
1149 for (Int_t iJet=0; iJet<tmpAODjets->GetEntriesFast(); iJet++){
1150 AliAODJet *jet = dynamic_cast<AliAODJet*>((*tmpAODjets)[iJet]);
1151 if(!jet) continue;
1152 TRefArray *trackList = jet->GetRefTracks();
1153 Int_t nTracks = trackList->GetEntriesFast();
1154 nInputTracks += nTracks;
1155 if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks);
1156 }
1157 if(fDebug) Printf("---> input tracks: %d", nInputTracks);
1158
1159 return nInputTracks;
1160}
1161
1162
1163
ea693273 1164Double_t AliAnalysisTaskJetCore::RelativePhi(Double_t mphi,Double_t vphi){
1165
1166 if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
1167 else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
1168 if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
1169 else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
1170 double dphi = mphi-vphi;
1171 if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
1172 else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
1173 return dphi;//dphi in [-Pi, Pi]
1174}
1175
d5e7475c 1176Int_t AliAnalysisTaskJetCore::GetPhiBin(Double_t phi)
1177{
1178 Int_t phibin=-1;
1179 if(!(TMath::Abs(phi)<=2*TMath::Pi())){AliError("phi w.r.t. RP out of defined range");return -1;}
1180 Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi)));
1181 phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi()));
1182 if(phibin<0||phibin>=fNRPBins){AliError("Phi Bin not defined");}
1183 return phibin;
1184}
1185
1186
75bf77e3 1187
1188
da93bb11 1189THnSparse* AliAnalysisTaskJetCore::NewTHnSparseF(const char* name, UInt_t entries)
1190{
1191 // generate new THnSparseF, axes are defined in GetDimParams()
1192
1193 Int_t count = 0;
1194 UInt_t tmp = entries;
1195 while(tmp!=0){
1196 count++;
1197 tmp = tmp &~ -tmp; // clear lowest bit
1198 }
1199
1200 TString hnTitle(name);
1201 const Int_t dim = count;
1202 Int_t nbins[dim];
1203 Double_t xmin[dim];
1204 Double_t xmax[dim];
1205
1206 Int_t i=0;
1207 Int_t c=0;
1208 while(c<dim && i<32){
1209 if(entries&(1<<i)){
1210
1211 TString label("");
1212 GetDimParams(i, label, nbins[c], xmin[c], xmax[c]);
1213 hnTitle += Form(";%s",label.Data());
1214 c++;
1215 }
1216
1217 i++;
1218 }
1219 hnTitle += ";";
1220
1221 return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
1222}
1223
1224void AliAnalysisTaskJetCore::GetDimParams(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
1225{
1226 // stores label and binning of axis for THnSparse
1227
1228 const Double_t pi = TMath::Pi();
1229
1230 switch(iEntry){
1231
1232 case 0:
1233 label = "V0 centrality (%)";
1234
1235 nbins = 10;
1236 xmin = 0.;
1237 xmax = 100.;
1238 break;
1239
1240
1241 case 1:
1242 label = "corrected jet pt";
c42c6838 1243 nbins = 20;
da93bb11 1244 xmin = 0.;
1245 xmax = 200.;
1246 break;
1247
1248
1249 case 2:
1250 label = "track pT";
1251
ba143e7f 1252 nbins = 9;
da93bb11 1253 xmin = 0.;
5bd732d4 1254 xmax = 150;
da93bb11 1255 break;
1256
1257
529e2916 1258 case 3:
da93bb11 1259 label = "deltaR";
1260 nbins = 15;
1261 xmin = 0.;
1262 xmax = 1.5;
1263 break;
529e2916 1264
1265
1266
da93bb11 1267 case 4:
1268 label = "deltaEta";
ba143e7f 1269 nbins = 8;
1270 xmin = -1.6;
1271 xmax = 1.6;
da93bb11 1272 break;
1273
1274
5bd732d4 1275 case 5:
da93bb11 1276 label = "deltaPhi";
1277 nbins = 90;
1278 xmin = -0.5*pi;
1279 xmax = 1.5*pi;
1280 break;
1281
1282
529e2916 1283
1284 case 6:
da93bb11 1285 label = "leading track";
ba143e7f 1286 nbins = 13;
da93bb11 1287 xmin = 0;
20dcc500 1288 xmax = 50;
da93bb11 1289 break;
1290
529e2916 1291 case 7:
da93bb11 1292
1293 label = "trigger track";
ba143e7f 1294 nbins =10;
da93bb11 1295 xmin = 0;
1296 xmax = 50;
1297 break;
529e2916 1298
1299
1300
1301
1302
1303
1304
1305
da93bb11 1306 }
1307
1308}
1309