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