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