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