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