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