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