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