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