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