simplify the logic to determine the origin of a decay photon or merged photon, fix...
[u/mrichter/AliRoot.git] / PWGJE / AliAnalysisTaskJetCorePP.cxx
CommitLineData
ad869500 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"
80ac66f6 29#include "TList.h"
ad869500 30#include "TMath.h"
31#include "TH1F.h"
32#include "TH1D.h"
33#include "TH2F.h"
34#include "TH3F.h"
35#include "THnSparse.h"
36#include "TCanvas.h"
80ac66f6 37#include "TArrayI.h"
ea603b64 38#include "TProfile.h"
39#include "TFile.h"
40#include "TKey.h"
8e103a56 41#include "TRandom3.h"
ad869500 42
43#include "AliLog.h"
44
45#include "AliAnalysisTask.h"
46#include "AliAnalysisManager.h"
47
48#include "AliVEvent.h"
49#include "AliESDEvent.h"
80ac66f6 50#include "AliMCEvent.h"
ad869500 51#include "AliESDInputHandler.h"
52#include "AliCentrality.h"
53#include "AliAnalysisHelperJetTasks.h"
54#include "AliInputEventHandler.h"
55#include "AliAODJetEventBackground.h"
ea603b64 56#include "AliGenPythiaEventHeader.h"
ad869500 57#include "AliAODMCParticle.h"
80ac66f6 58#include "AliMCParticle.h"
ad869500 59#include "AliAODEvent.h"
60#include "AliAODHandler.h"
61#include "AliAODJet.h"
62
63#include "AliAnalysisTaskJetCorePP.h"
64
65using std::cout;
66using std::endl;
67
68ClassImp(AliAnalysisTaskJetCorePP)
69
70//Filip Krizek 1st March 2013
71
72//---------------------------------------------------------------------
73AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP() :
74AliAnalysisTaskSE(),
75fESD(0x0),
76fAODIn(0x0),
77fAODOut(0x0),
78fAODExtension(0x0),
79fJetBranchName(""),
80ac66f6 80fJetBranchNameMC(""),
7fc3d134 81fListJets(0x0),
80ac66f6 82fListJetsGen(0x0),
ad869500 83fNonStdFile(""),
84fSystem(0), //pp=0 pPb=1
85fJetParamR(0.4),
86fOfflineTrgMask(AliVEvent::kAny),
87fMinContribVtx(1),
88fVtxZMin(-10.0),
89fVtxZMax(10.0),
90fFilterMask(0),
91fCentMin(0.0),
92fCentMax(100.0),
93fJetEtaMin(-0.5),
94fJetEtaMax(0.5),
95fTriggerEtaCut(0.9),
96fTrackEtaCut(0.9),
97fTrackLowPtCut(0.15),
98fOutputList(0x0),
99fHistEvtSelection(0x0),
100fh2Ntriggers(0x0),
101fHJetSpec(0x0),
102fHJetDensity(0x0),
103fHJetDensityA4(0x0),
104fhJetPhi(0x0),
105fhTriggerPhi(0x0),
106fhJetEta(0x0),
107fhTriggerEta(0x0),
108fhVertexZ(0x0),
109fhVertexZAccept(0x0),
110fhContribVtx(0x0),
111fhContribVtxAccept(0x0),
112fhDphiTriggerJet(0x0),
113fhDphiTriggerJetAccept(0x0),
114fhCentrality(0x0),
115fhCentralityAccept(0x0),
7fc3d134 116fHJetPtRaw(0x0),
117fHLeadingJetPtRaw(0x0),
118fHDphiVsJetPtAll(0x0),
80ac66f6 119fhJetPtGenVsJetPtRec(0x0),
120fhJetPtGen(0x0),
121fh2NtriggersGen(0x0),
122fHJetSpecGen(0x0),
123fhPtTrkTruePrimRec(0x0),
124fhPtTrkTruePrimGen(0x0),
125fhPtTrkSecOrFakeRec(0x0),
126fIsMC(0),
127faGenIndex(0),
128faRecIndex(0),
7fc3d134 129fkAcceptance(2.0*TMath::Pi()*1.8),
ea603b64 130fkDeltaPhiCut(TMath::Pi()-0.6),
131fh1Xsec(0x0),
132fh1Trials(0x0),
133fh1AvgTrials(0x0),
134fh1PtHard(0x0),
135fh1PtHardNoW(0x0),
136fh1PtHardTrials(0x0),
8e103a56 137fAvgTrials(1),
138fHardest(0),
139fEventNumberRangeLow(0),
140fEventNumberRangeHigh(99),
141fTriggerPtRangeLow(0.0),
142fTriggerPtRangeHigh(10000.0),
143fRandom(0x0)
ad869500 144{
145 // default Constructor
ad869500 146}
147
148//---------------------------------------------------------------------
149
150AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP(const char *name) :
151AliAnalysisTaskSE(name),
152fESD(0x0),
153fAODIn(0x0),
154fAODOut(0x0),
155fAODExtension(0x0),
156fJetBranchName(""),
80ac66f6 157fJetBranchNameMC(""),
7fc3d134 158fListJets(0x0),
80ac66f6 159fListJetsGen(0x0),
ad869500 160fNonStdFile(""),
161fSystem(0), //pp=0 pPb=1
162fJetParamR(0.4),
163fOfflineTrgMask(AliVEvent::kAny),
164fMinContribVtx(1),
165fVtxZMin(-10.0),
166fVtxZMax(10.0),
167fFilterMask(0),
168fCentMin(0.0),
169fCentMax(100.0),
170fJetEtaMin(-0.5),
171fJetEtaMax(0.5),
172fTriggerEtaCut(0.9),
173fTrackEtaCut(0.9),
174fTrackLowPtCut(0.15),
175fOutputList(0x0),
176fHistEvtSelection(0x0),
177fh2Ntriggers(0x0),
178fHJetSpec(0x0),
179fHJetDensity(0x0),
180fHJetDensityA4(0x0),
181fhJetPhi(0x0),
182fhTriggerPhi(0x0),
183fhJetEta(0x0),
184fhTriggerEta(0x0),
185fhVertexZ(0x0),
186fhVertexZAccept(0x0),
187fhContribVtx(0x0),
188fhContribVtxAccept(0x0),
189fhDphiTriggerJet(0x0),
190fhDphiTriggerJetAccept(0x0),
191fhCentrality(0x0),
192fhCentralityAccept(0x0),
7fc3d134 193fHJetPtRaw(0x0),
194fHLeadingJetPtRaw(0x0),
195fHDphiVsJetPtAll(0x0),
80ac66f6 196fhJetPtGenVsJetPtRec(0x0),
197fhJetPtGen(0x0),
198fh2NtriggersGen(0x0),
199fHJetSpecGen(0x0),
200fhPtTrkTruePrimRec(0x0),
201fhPtTrkTruePrimGen(0x0),
202fhPtTrkSecOrFakeRec(0x0),
203fIsMC(0),
204faGenIndex(0),
205faRecIndex(0),
7fc3d134 206fkAcceptance(2.0*TMath::Pi()*1.8),
ea603b64 207fkDeltaPhiCut(TMath::Pi()-0.6),
208fh1Xsec(0x0),
209fh1Trials(0x0),
210fh1AvgTrials(0x0),
211fh1PtHard(0x0),
212fh1PtHardNoW(0x0),
213fh1PtHardTrials(0x0),
8e103a56 214fAvgTrials(1),
215fHardest(0),
216fEventNumberRangeLow(0),
217fEventNumberRangeHigh(99),
218fTriggerPtRangeLow(0.0),
219fTriggerPtRangeHigh(10000.0),
220fRandom(0x0)
ad869500 221{
7fc3d134 222// Constructor
ad869500 223
224 DefineOutput(1, TList::Class());
225}
226
227//--------------------------------------------------------------
228AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP(const AliAnalysisTaskJetCorePP& a):
229AliAnalysisTaskSE(a.GetName()),
230fESD(a.fESD),
231fAODIn(a.fAODIn),
232fAODOut(a.fAODOut),
233fAODExtension(a.fAODExtension),
234fJetBranchName(a.fJetBranchName),
80ac66f6 235fJetBranchNameMC(a.fJetBranchNameMC),
ad869500 236fListJets(a.fListJets),
80ac66f6 237fListJetsGen(a.fListJetsGen),
ad869500 238fNonStdFile(a.fNonStdFile),
239fSystem(a.fSystem),
240fJetParamR(a.fJetParamR),
241fOfflineTrgMask(a.fOfflineTrgMask),
242fMinContribVtx(a.fMinContribVtx),
243fVtxZMin(a.fVtxZMin),
244fVtxZMax(a.fVtxZMax),
245fFilterMask(a.fFilterMask),
246fCentMin(a.fCentMin),
247fCentMax(a.fCentMax),
248fJetEtaMin(a.fJetEtaMin),
249fJetEtaMax(a.fJetEtaMax),
250fTriggerEtaCut(a.fTriggerEtaCut),
251fTrackEtaCut(a.fTrackEtaCut),
252fTrackLowPtCut(a.fTrackLowPtCut),
253fOutputList(a.fOutputList),
254fHistEvtSelection(a.fHistEvtSelection),
255fh2Ntriggers(a.fh2Ntriggers),
256fHJetSpec(a.fHJetSpec),
257fHJetDensity(a.fHJetDensity),
258fHJetDensityA4(a.fHJetDensityA4),
259fhJetPhi(a.fhJetPhi),
260fhTriggerPhi(a.fhTriggerPhi),
261fhJetEta(a.fhJetEta),
262fhTriggerEta(a.fhTriggerEta),
263fhVertexZ(a.fhVertexZ),
264fhVertexZAccept(a.fhVertexZAccept),
265fhContribVtx(a.fhContribVtx),
266fhContribVtxAccept(a.fhContribVtxAccept),
267fhDphiTriggerJet(a.fhDphiTriggerJet),
268fhDphiTriggerJetAccept(a.fhDphiTriggerJetAccept),
269fhCentrality(a.fhCentrality),
270fhCentralityAccept(a.fhCentralityAccept),
7fc3d134 271fHJetPtRaw(a.fHJetPtRaw),
272fHLeadingJetPtRaw(a.fHLeadingJetPtRaw),
273fHDphiVsJetPtAll(a.fHDphiVsJetPtAll),
80ac66f6 274fhJetPtGenVsJetPtRec(a.fhJetPtGenVsJetPtRec),
275fhJetPtGen(a.fhJetPtGen),
276fh2NtriggersGen(a.fh2NtriggersGen),
277fHJetSpecGen(a.fHJetSpecGen),
278fhPtTrkTruePrimRec(a.fhPtTrkTruePrimRec),
279fhPtTrkTruePrimGen(a.fhPtTrkTruePrimGen),
280fhPtTrkSecOrFakeRec(a.fhPtTrkSecOrFakeRec),
281fIsMC(a.fIsMC),
282faGenIndex(a.faGenIndex),
283faRecIndex(a.faRecIndex),
7fc3d134 284fkAcceptance(a.fkAcceptance),
ea603b64 285fkDeltaPhiCut(a.fkDeltaPhiCut),
286fh1Xsec(a.fh1Xsec),
287fh1Trials(a.fh1Trials),
288fh1AvgTrials(a.fh1AvgTrials),
289fh1PtHard(a.fh1PtHard),
290fh1PtHardNoW(a.fh1PtHardNoW),
291fh1PtHardTrials(a.fh1PtHardTrials),
8e103a56 292fAvgTrials(a.fAvgTrials),
293fHardest(a.fHardest),
294fEventNumberRangeLow(a.fEventNumberRangeLow),
295fEventNumberRangeHigh(a.fEventNumberRangeHigh),
296fTriggerPtRangeLow(a.fTriggerPtRangeLow),
43457761 297fTriggerPtRangeHigh(a.fTriggerPtRangeHigh),
298fRandom(a.fRandom)
ad869500 299{
300 //Copy Constructor
43457761 301 fRandom->SetSeed(0);
ad869500 302}
303//--------------------------------------------------------------
304
305AliAnalysisTaskJetCorePP& AliAnalysisTaskJetCorePP::operator = (const AliAnalysisTaskJetCorePP& a){
306 // assignment operator
307 this->~AliAnalysisTaskJetCorePP();
308 new(this) AliAnalysisTaskJetCorePP(a);
309 return *this;
310}
311//--------------------------------------------------------------
312
313AliAnalysisTaskJetCorePP::~AliAnalysisTaskJetCorePP()
314{
315 //Destructor
316 delete fListJets;
80ac66f6 317 delete fListJetsGen;
ad869500 318 delete fOutputList; // ????
8e103a56 319 delete fRandom;
ad869500 320}
321
322//--------------------------------------------------------------
323
ea603b64 324
325Bool_t AliAnalysisTaskJetCorePP::Notify()
326{
327 //Implemented Notify() to read the cross sections
328 //and number of trials from pyxsec.root
329 //inspired by AliAnalysisTaskJetSpectrum2::Notify()
330 if(!fIsMC) return kFALSE;
331
332 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
333 if(!fESD){
334 if(fDebug>1) AliError("ESD not available");
335 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
336 }
337
338 fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
339
340
341 if(fNonStdFile.Length()!=0){
342 // case that we have an AOD extension we can fetch the jets from the extended output
343 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
344 fAODExtension = aodH ? aodH->GetExtension(fNonStdFile.Data()) : 0;
345 if(!fAODExtension){
346 if(fDebug>1) Printf("AODExtension found for %s",fNonStdFile.Data());
347 }
348 }
349
350 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
351 Float_t xsection = 0;
352 Float_t ftrials = 1;
353
354 fAvgTrials = 1;
355 if(tree){
356 TFile *curfile = tree->GetCurrentFile();
357 if(!curfile) {
358 Error("Notify","No current file");
359 return kFALSE;
360 }
361 if(!fh1Xsec || !fh1Trials){
362 Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
363 return kFALSE;
364 }
365 AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
366 fh1Xsec->Fill("<#sigma>",xsection);
367 // construct a poor man average trials
368 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
369 if(ftrials>=nEntries && nEntries>0.) fAvgTrials = ftrials/nEntries;
370 fh1Trials->Fill("#sum{ntrials}",ftrials);
371 }
372
373 if(fDebug)Printf("Reading File %s",fInputHandler->GetTree()->GetCurrentFile()->GetName());
374
375 return kTRUE;
376}
377//--------------------------------------------------------------
378
ad869500 379void AliAnalysisTaskJetCorePP::Init()
380{
381 // check for jet branches
382 if(!strlen(fJetBranchName.Data())){
383 AliError("Jet branch name not set.");
384 }
385
386}
387
388//--------------------------------------------------------------
389
390void AliAnalysisTaskJetCorePP::UserCreateOutputObjects()
391{
392
ad869500 393 // Create histograms
394 // Called once
80ac66f6 395 fListJets = new TList(); //reconstructed level
396
397 fIsMC = (fJetBranchNameMC.Length()>0) ? kTRUE : kFALSE;
398
8e103a56 399 fRandom = new TRandom3(0);
400
80ac66f6 401 if(fIsMC) fListJetsGen = new TList(); //generator level
ad869500 402 OpenFile(1);
403 if(!fOutputList) fOutputList = new TList();
404 fOutputList->SetOwner(kTRUE);
405
406 Bool_t oldStatus = TH1::AddDirectoryStatus();
407 TH1::AddDirectory(kFALSE);
408
409 fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
410 fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
411 fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
412 fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
413 fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
414 fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
8e103a56 415 fHistEvtSelection->GetXaxis()->SetBinLabel(6,"event number (rejected)");
ad869500 416
417 fOutputList->Add(fHistEvtSelection);
418
419 Int_t nBinsCentrality = (fSystem==0) ? 1 : 10; // pp=1 else 10
80ac66f6 420 //trigger pt spectrum (reconstructed)
ad869500 421 fh2Ntriggers = new TH2F("fh2Ntriggers","# of triggers",
422 nBinsCentrality,0.0,100.0,50,0.0,50.0);
7fc3d134 423 fOutputList->Add(fh2Ntriggers);
ad869500 424
80ac66f6 425 //Centrality, A, pTjet, pTtrigg
426 const Int_t dimSpec = 4;
427 const Int_t nBinsSpec[dimSpec] = {nBinsCentrality, 100, 100, 50};
428 const Double_t lowBinSpec[dimSpec] = {0.0, 0.0, 0, 0.0};
429 const Double_t hiBinSpec[dimSpec] = {100.0, 1.0, 200.0,50.0};
7fc3d134 430 fHJetSpec = new THnSparseF("fHJetSpec",
80ac66f6 431 "Recoil jet spectrum [cent,A,pTjet,pTtrig]",
7fc3d134 432 dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
ad869500 433 fOutputList->Add(fHJetSpec);
434
435 //------------------- HISTOS FOR DIAGNOSTIC ----------------------
436 //Jet number density histos [Trk Mult, jet density, pT trigger]
437 const Int_t dimJetDens = 3;
80ac66f6 438 const Int_t nBinsJetDens[dimJetDens] = {100, 100, 10};
ad869500 439 const Double_t lowBinJetDens[dimJetDens] = {0.0, 0.0, 0.0};
80ac66f6 440 const Double_t hiBinJetDens[dimJetDens] = {500.0, 5.0, 50.0};
ad869500 441
442 fHJetDensity = new THnSparseF("fHJetDensity","Jet dens vs trk mult A>0.07",
443 dimJetDens,nBinsJetDens,lowBinJetDens,hiBinJetDens);
444
445 fHJetDensityA4 =new THnSparseF("fHJetDensityA4","Jet dens vs trk mult A>0.4",
446 dimJetDens,nBinsJetDens,lowBinJetDens,hiBinJetDens);
447
ad869500 448 fOutputList->Add(fHJetDensity);
449 fOutputList->Add(fHJetDensityA4);
450
451
452 //inclusive azimuthal and pseudorapidity histograms
7fc3d134 453 fhJetPhi = new TH2D("fhJetPhi","Azim dist jets vs pTjet",
454 50, 0, 100, 50,-TMath::Pi(),TMath::Pi());
455 fhTriggerPhi= new TH2D("fhTriggerPhi","azim dist trig had vs pTtrigg",
456 25, 0, 50, 50,-TMath::Pi(),TMath::Pi());
457 fhJetEta = new TH2D("fhJetEta","Eta dist jets vs pTjet",
458 50,0, 100, 40,-0.9,0.9);
459 fhTriggerEta = new TH2D("fhTriggerEta","Eta dist trig had vs pTtrigg",
460 25, 0, 50, 40,-0.9,0.9);
461
ad869500 462 fhVertexZ = new TH1D("fhVertexZ","z vertex",40,-20,20);
463 fhVertexZAccept = new TH1D("fhVertexZAccept","z vertex after cut",40,-20,20);
464 fhContribVtx = new TH1D("fhContribVtx","contrib to vtx",200,0,200);
465 fhContribVtxAccept = new TH1D("fhContribVtxAccept","contrib to vtx after cut",200,0,200);
466 fhDphiTriggerJet = new TH1D("fhDphiTriggerJet","Deltaphi trig-jet",50, -TMath::Pi(),TMath::Pi());
467 fhDphiTriggerJetAccept = new TH1D("fhDphiTriggerJetAccept","Deltaphi trig-jet after cut",50, -TMath::Pi(),TMath::Pi());
468 fhCentrality = new TH1D("fhCentrality","Centrality",20,0,100);
469 fhCentralityAccept = new TH1D("fhCentralityAccept","Centrality after cut",20,0,100);
470
471 fOutputList->Add(fhJetPhi);
472 fOutputList->Add(fhTriggerPhi);
473 fOutputList->Add(fhJetEta);
474 fOutputList->Add(fhTriggerEta);
475 fOutputList->Add(fhVertexZ);
476 fOutputList->Add(fhVertexZAccept);
477 fOutputList->Add(fhContribVtx);
478 fOutputList->Add(fhContribVtxAccept);
479 fOutputList->Add(fhDphiTriggerJet);
480 fOutputList->Add(fhDphiTriggerJetAccept);
481 fOutputList->Add(fhCentrality);
482 fOutputList->Add(fhCentralityAccept);
483
7fc3d134 484 // raw spectra of INCLUSIVE jets
80ac66f6 485 //Centrality, pTjet, A
486 const Int_t dimRaw = 3;
487 const Int_t nBinsRaw[dimRaw] = {nBinsCentrality, 50, 100};
488 const Double_t lowBinRaw[dimRaw] = {0.0, 0.0, 0.0};
489 const Double_t hiBinRaw[dimRaw] = {100.0, 100, 1.0};
7fc3d134 490 fHJetPtRaw = new THnSparseF("fHJetPtRaw",
80ac66f6 491 "Incl. jet spectrum [cent,pTjet,A]",
7fc3d134 492 dimRaw,nBinsRaw,lowBinRaw,hiBinRaw);
493 fOutputList->Add(fHJetPtRaw);
494
495 // raw spectra of LEADING jets
80ac66f6 496 //Centrality, pTjet, A
7fc3d134 497 fHLeadingJetPtRaw = new THnSparseF("fHLeadingJetPtRaw",
80ac66f6 498 "Leading jet spectrum [cent,pTjet,A]",
7fc3d134 499 dimRaw,nBinsRaw,lowBinRaw,hiBinRaw);
500 fOutputList->Add(fHLeadingJetPtRaw);
501
502 // Dphi versus pT jet
503 //Centrality, Dphi=phiTrig-phiJet, pTjet, pTtrigg
504 const Int_t dimDp = 4;
505 const Int_t nBinsDp[dimDp] = {nBinsCentrality, 50, 50, 50};
506 const Double_t lowBinDp[dimDp] = {0.0, -TMath::Pi(), 0.0, 0.0};
507 const Double_t hiBinDp[dimDp] = {100.0, TMath::Pi(), 100.0, 100.0};
508 fHDphiVsJetPtAll = new THnSparseF("fHDphiVsJetPtAll",
509 "Dphi vs jet pT [cent,Dphi,pTjet,pTtrigg]",
510 dimDp,nBinsDp,lowBinDp,hiBinDp);
511 fOutputList->Add(fHDphiVsJetPtAll);
512
7fc3d134 513
80ac66f6 514 //analyze MC generator level
515 if(fIsMC){
516 fhJetPtGenVsJetPtRec = new TH2D("fhJetPtGenVsJetPtRec","JetPtGenVsJetPtRec", 100,0,200, 100,0,200);
517 fOutputList->Add(fhJetPtGenVsJetPtRec); //gen MC jet pt spectrum versus reconstructed pt spectrum
518
519 fhJetPtGen = new TH1D("fhJetPtGen","Jet Pt (MC Gen)",100,0,200); //MC generator jet pt spectrum
520 fOutputList->Add(fhJetPtGen);
7fc3d134 521
80ac66f6 522 fh2NtriggersGen = (TH2F*) fh2Ntriggers->Clone("fh2NtriggersGen");
523 fh2NtriggersGen->SetTitle(Form("%s Gen MC",fh2Ntriggers->GetTitle()));
524 fOutputList->Add(fh2NtriggersGen);
7fc3d134 525
80ac66f6 526 //Centrality, A, pT, pTtrigg
527 fHJetSpecGen = (THnSparseF*) fHJetSpec->Clone("fHJetSpecGen");
528 fHJetSpecGen->SetTitle(Form("%s Gen MC",fHJetSpec->GetTitle()));
529 fOutputList->Add(fHJetSpecGen);
530
531 //track efficiency/contamination histograms eta versus pT
532 fhPtTrkTruePrimRec = new TH2D("fhPtTrkTruePrimRec","PtTrkTruePrimRec",100,0,20,18,-0.9,0.9);
533 fOutputList->Add(fhPtTrkTruePrimRec);
534
535 fhPtTrkTruePrimGen = (TH2D*) fhPtTrkTruePrimRec->Clone("fhPtTrkTruePrimGen");
536 fhPtTrkTruePrimGen->SetTitle("PtTrkTruePrimGen");
537 fOutputList->Add(fhPtTrkTruePrimGen);
538
539 fhPtTrkSecOrFakeRec = (TH2D*) fhPtTrkTruePrimRec->Clone("fhPtTrkSecOrFakeRec");
540 fhPtTrkSecOrFakeRec->SetTitle("PtTrkSecOrFakeRec");
541 fOutputList->Add(fhPtTrkSecOrFakeRec);
542 }
ea603b64 543 //-------------------------------------
544 // pythia histograms
545 const Int_t nBinPt = 150;
546 Double_t binLimitsPt[nBinPt+1];
547 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
548 if(iPt == 0){
549 binLimitsPt[iPt] = -50.0;
550 }else{// 1.0
551 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.;
552 }
553 }
554
555 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
556 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
557 fOutputList->Add(fh1Xsec);
558 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
559 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
560 fOutputList->Add(fh1Trials);
561 fh1AvgTrials = new TH1F("fh1AvgTrials","trials root file",1,0,1);
562 fh1AvgTrials->GetXaxis()->SetBinLabel(1,"#sum{avg ntrials}");
563 fOutputList->Add(fh1AvgTrials);
564 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
565 fOutputList->Add(fh1PtHard);
566 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
567 fOutputList->Add(fh1PtHardNoW);
568 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
569 fOutputList->Add(fh1PtHardTrials);
570
ad869500 571
572 // =========== Switch on Sumw2 for all histos ===========
573 for(Int_t i=0; i<fOutputList->GetEntries(); i++){
574 TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
575 if(h1){
576 h1->Sumw2();
577 continue;
578 }
579 THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
580 if(hn){
581 hn->Sumw2();
582 }
583 }
584 TH1::AddDirectory(oldStatus);
585
586 PostData(1, fOutputList);
587}
588
589//--------------------------------------------------------------------
590
591void AliAnalysisTaskJetCorePP::UserExec(Option_t *)
592{
ad869500 593 //Event loop
ea603b64 594 Double_t eventW = 1.0;
595 Double_t ptHard = 0.0;
596 Double_t nTrials = 1.0; // Trials for MC trigger
597 if(fIsMC) fh1AvgTrials->Fill("#sum{avg ntrials}",fAvgTrials);
ad869500 598
599 if(TMath::Abs((Float_t) fJetParamR)<0.00001){
80ac66f6 600 AliError("Cone radius is set to zero.");
ad869500 601 return;
602 }
603 if(!strlen(fJetBranchName.Data())){
604 AliError("Jet branch name not set.");
605 return;
606 }
607
608 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
609 if(!fESD){
80ac66f6 610 if(fDebug>1) AliError("ESD not available");
ad869500 611 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
612 }
613
614 fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
615 AliAODEvent* aod = NULL;
616 // take all other information from the aod we take the tracks from
617 if(!aod){
618 if(!fESD) aod = fAODIn;
619 else aod = fAODOut;
620 }
621
622 if(fNonStdFile.Length()!=0){
623 // case that we have an AOD extension we can fetch the jets from the extended output
624 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
625 fAODExtension = aodH ? aodH->GetExtension(fNonStdFile.Data()) : 0;
626 if(!fAODExtension){
627 if(fDebug>1) Printf("AODExtension found for %s",fNonStdFile.Data());
628 }
629 }
630
80ac66f6 631 // ----------------- EVENT SELECTION --------------------------
ad869500 632 fHistEvtSelection->Fill(1); // number of events before event selection
633
634 // physics selection
635 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
636 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
637
638 if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
639 if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
640 fHistEvtSelection->Fill(2);
641 PostData(1, fOutputList);
642 return;
643 }
644
645 //check AOD pointer
646 if(!aod){
647 if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
648 fHistEvtSelection->Fill(3);
649 PostData(1, fOutputList);
650 return;
651 }
652
653 // vertex selection
654 AliAODVertex* primVtx = aod->GetPrimaryVertex();
655
656 if(!primVtx){
657 if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
658 fHistEvtSelection->Fill(3);
659 PostData(1, fOutputList);
660 return;
661 }
662
663 Int_t nTracksPrim = primVtx->GetNContributors();
664 Float_t vtxz = primVtx->GetZ();
665 //Input events
666 fhContribVtx->Fill(nTracksPrim);
667 if( nTracksPrim > 0 ) fhVertexZ->Fill(vtxz);
668
669 if((nTracksPrim < fMinContribVtx) ||
670 (vtxz < fVtxZMin) ||
671 (vtxz > fVtxZMax)){
672 if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",
673 (char*)__FILE__,__LINE__,vtxz);
674 fHistEvtSelection->Fill(3);
675 PostData(1, fOutputList);
676 return;
677 }else{
678 //Accepted events
679 fhContribVtxAccept->Fill(nTracksPrim);
680 fhVertexZAccept->Fill(vtxz);
681 }
682 //FK// No event class selection imposed
683 // event class selection (from jet helper task)
684 //Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
685 //if(fDebug) Printf("Event class %d", eventClass);
686 //if(eventClass < fEvtClassMin || eventClass > fEvtClassMax){
687 // fHistEvtSelection->Fill(4);
688 // PostData(1, fOutputList);
689 // return;
690 //}
691
80ac66f6 692 //------------------ CENTRALITY SELECTION ---------------
ad869500 693 AliCentrality *cent = 0x0;
694 Double_t centValue = 0.0;
695 if(fSystem){ //fSystem=0 for pp, fSystem=1 for pPb
696 if(fESD){
697 cent = fESD->GetCentrality();
698 if(cent) centValue = cent->GetCentralityPercentile("V0M");
699 }else{
700 centValue = aod->GetHeader()->GetCentrality();
701 }
702 if(fDebug) printf("centrality: %f\n", centValue);
703 //Input events
704 fhCentrality->Fill(centValue);
705
706 if(centValue < fCentMin || centValue > fCentMax){
707 fHistEvtSelection->Fill(4);
708 PostData(1, fOutputList);
709 return;
710 }else{
711 //Accepted events
712 fhCentralityAccept->Fill( centValue );
713 }
714 }
715
8e103a56 716 //-----------------select disjunct event subsamples ----------------
717 Int_t eventnum = aod->GetHeader()->GetEventNumberESDFile();
718 Int_t lastdigit = eventnum % 10;
719 if(!(fEventNumberRangeLow<=lastdigit && lastdigit<=fEventNumberRangeHigh)){
720 fHistEvtSelection->Fill(5);
721 PostData(1, fOutputList);
722 return;
723 }
724
ad869500 725 if(fDebug) std::cout<<" ACCEPTED EVENT "<<endl;
ad869500 726 fHistEvtSelection->Fill(0); // accepted events
ad869500 727 // ------------------- end event selection --------------------
80ac66f6 728
ad869500 729
80ac66f6 730 // fetch RECONSTRUCTED jets
ad869500 731 TClonesArray *aodJets = 0x0;
732
733 if(fAODOut && !aodJets){
734 aodJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName.Data()));
735 }
736 if(fAODExtension && !aodJets){
737 aodJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName.Data()));
738 }
739 if(fAODIn && !aodJets){
740 aodJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName.Data()));
741 }
742
80ac66f6 743 //--------- Fill list of RECONSTRUCTED jets --------------
744 fListJets->Clear();
745 if(aodJets){
746 if(fDebug) Printf("########## %s: %d jets",fJetBranchName.Data(), aodJets->GetEntriesFast());
747 for(Int_t iJet = 0; iJet < aodJets->GetEntriesFast(); iJet++) {
748 AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets)[iJet]);
749 if (jet) fListJets->Add(jet);
750 }
751 }
752
753 //--------- Fill list of MC GENERATOR LEVEL jets --------------
754 TList particleListGen; //list of tracks in MC
755
756 if(fIsMC){ //analyze generator level MC
757 // fetch MC generator level jets
758 TClonesArray *aodGenJets = NULL;
759
760 if(fAODOut&&!aodGenJets){
761 aodGenJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchNameMC.Data()));
762 }
763 if(fAODExtension&&!aodGenJets){
764 aodGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchNameMC.Data()));
765 }
766 if(fAODIn&&!aodGenJets){
767 aodGenJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchNameMC.Data()));
768 }
769
770 if(!aodGenJets){
771 Printf("%s:%d no generated Jet array with name %s in AOD",
772 (char*)__FILE__,__LINE__, fJetBranchNameMC.Data());
773 PostData(1, fOutputList);
774 return;
775 }
776
777 fListJetsGen->Clear();
778
779 //serarch for charged trigger at the MC generator level
8e103a56 780 Int_t indexTriggGen = -1;
781 Double_t ptTriggGen = -1;
782 Int_t iCounterGen = 0;
783 Int_t triggersMC[200];//list of trigger candidates
784 Int_t ntriggersMC = 0; //index in triggers array
785
80ac66f6 786 if(fESD){//ESD input
787
788 AliMCEvent* mcEvent = MCEvent();
789 if(!mcEvent){
790 PostData(1, fOutputList);
791 return;
ea603b64 792 }
793
794 AliGenPythiaEventHeader *pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
795 if(pythiaGenHeader){
796 nTrials = pythiaGenHeader->Trials();
797 ptHard = pythiaGenHeader->GetPtHard();
798 fh1PtHard->Fill(ptHard,eventW);
799 fh1PtHardNoW->Fill(ptHard,1);
800 fh1PtHardTrials->Fill(ptHard,nTrials);
801 }
802
80ac66f6 803 for(Int_t it = 0; it < mcEvent->GetNumberOfTracks(); it++){
804 if(!mcEvent->IsPhysicalPrimary(it)) continue;
805 AliMCParticle* part = (AliMCParticle*) mcEvent->GetTrack(it);
806 if(!part) continue;
807 if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){
8e103a56 808
809 if(fHardest==0 && ntriggersMC<200){//single inclusive trigger
810 if(indexTriggGen > -1){//trigger candidater was found
811 triggersMC[ntriggersMC] = indexTriggGen;
812 ntriggersMC++;
813 }
814 }
815
80ac66f6 816 iCounterGen++;
817 }
818 }
819 }else{ //AOD input
820 if(!fAODIn){
821 PostData(1, fOutputList);
822 return;
823 }
824 TClonesArray *tca = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(AliAODMCParticle::StdBranchName()));
825 if(!tca){
826 PostData(1, fOutputList);
827 return;
828 }
829 for(Int_t it = 0; it < tca->GetEntriesFast(); it++){
830 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
831 if(!part) continue;
832 if(!part->IsPhysicalPrimary()) continue;
833 if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){
834
8e103a56 835 if(fHardest==0 && ntriggersMC<200){//single inclusive trigger
836 if(indexTriggGen > -1){ //trigger candidater was found
837 triggersMC[ntriggersMC] = indexTriggGen;
838 ntriggersMC++;
839 }
840 }
841
80ac66f6 842 iCounterGen++;
843 }
844 }
845 }
8e103a56 846
847 if(fHardest==0){ //single inclusive trigger
848 if(ntriggersMC>0){ //there is at least one trigger
849 Int_t rnd = fRandom->Integer(ntriggersMC); //0 to ntriggers-1
850 indexTriggGen = triggersMC[rnd];
80ac66f6 851 }else{
8e103a56 852 indexTriggGen = -1; //trigger not found
80ac66f6 853 }
854 }
8e103a56 855
856 //================== Fill jet list ===================
80ac66f6 857 if(aodGenJets){
858 if(fDebug) Printf("########## %s: %d jets",fJetBranchNameMC.Data(), aodGenJets->GetEntriesFast());
859
860 for(Int_t igJet = 0; igJet < aodGenJets->GetEntriesFast(); igJet++) {
861 AliAODJet *jetGen = dynamic_cast<AliAODJet*>((*aodGenJets)[igJet]);
862 if(jetGen) fListJetsGen->Add(jetGen);
863 }
8e103a56 864 }
865
866 //============ Generator trigger+jet ==================
867 Int_t ilowGen = (fHardest==0 || fHardest==1) ? indexTriggGen : 0;
868 Int_t ihighGen = (fHardest==0 || fHardest==1) ? indexTriggGen+1 : particleListGen.GetEntries();
869
870 for(Int_t igen= ilowGen; igen< ihighGen; igen++){ //loop over possible trigger
871 indexTriggGen = igen; //trigger hadron
872
873 if(indexTriggGen == -1) continue;
874 AliVParticle* triggerGen=(AliVParticle*) particleListGen.At(indexTriggGen);
875 if(!triggerGen) continue;
876
877 if(fHardest >= 2){
8f9e503d 878 if(triggerGen->Pt() < 6.0) continue; //all hadrons pt>10
8e103a56 879 }
880 if(TMath::Abs((Float_t) triggerGen->Eta()) > fTriggerEtaCut) continue;
881
882 ptTriggGen = triggerGen->Pt(); //count triggers
883 fh2NtriggersGen->Fill(centValue, ptTriggGen);
884
80ac66f6 885 //Count jets and trigger-jet pairs at MC generator level
8e103a56 886 if(!aodGenJets) continue;
80ac66f6 887 for(Int_t ij=0; ij<fListJetsGen->GetEntries(); ij++){
888 AliAODJet* jet = (AliAODJet*)(fListJetsGen->At(ij));
889 if(!jet) continue;
890 Double_t etaJetGen = jet->Eta();
80ac66f6 891
892 if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){
80ac66f6 893
8e103a56 894 Double_t dphi = RelativePhi(triggerGen->Phi(), jet->Phi());
895 if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue;
80ac66f6 896
8e103a56 897 //Centrality, A, pT, pTtrigg
898 Double_t fillspecgen[] = { centValue,
80ac66f6 899 jet->EffectiveAreaCharged(),
900 jet->Pt(),
901 ptTriggGen
902 };
8e103a56 903 fHJetSpecGen->Fill(fillspecgen);
904 }//back to back jet-trigger pair
905 }//jet loop
906 }//trigger loop
907
908
909 //================ RESPONSE MATRIX ===============
910 if(aodGenJets){
911 //Count jets and trigger-jet pairs at MC generator level
912 for(Int_t ij=0; ij<fListJetsGen->GetEntries(); ij++){
913 AliAODJet* jet = (AliAODJet*)(fListJetsGen->At(ij));
914 if(!jet) continue;
915 Double_t etaJetGen = jet->Eta();
916 Double_t ptJetGen = jet->Pt();
917
918 if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){
919 fhJetPtGen->Fill(ptJetGen); // generator level pt spectum of jets response mx normalization
80ac66f6 920 }
921 }
80ac66f6 922 if(fListJets->GetEntries()>0 && fListJetsGen->GetEntries()>0){ //at least some reconstructed jets
923 Int_t ng = (Int_t) fListJetsGen->GetEntries();
924 Int_t nr = (Int_t) fListJets->GetEntries();
925
926 //Find closest MC generator - reconstructed jets
927 if(faGenIndex.GetSize()<nr) faGenIndex.Set(nr); //idx of gen jet assoc to rec jet
928 if(faRecIndex.GetSize()<ng) faRecIndex.Set(ng); //idx of rec jet assoc to gen jet
929
930 if(fDebug){
931 Printf("New Rec List %d gen index Array %d",nr,faGenIndex.GetSize());
932 Printf("New Gen List %d rec index Array %d",ng,faRecIndex.GetSize());
933 }
934 //matching of MC genrator level and reconstructed jets
935 AliAnalysisHelperJetTasks::GetClosestJets(fListJetsGen,ng,fListJets,nr,faGenIndex,faRecIndex,fDebug);
936
937 // Fill response matrix
938 for(Int_t ir = 0; ir < nr; ir++){
939 AliAODJet *recJet = (AliAODJet*) fListJets->At(ir);
940 Double_t etaJetRec = recJet->Eta();
941 Double_t ptJetRec = recJet->Pt();
942 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
943
944 if((fJetEtaMin <= etaJetRec) && (etaJetRec <= fJetEtaMax)){
945 Int_t ig = faGenIndex[ir]; //associated generator level jet
946 if(ig >= 0 && ig < ng){
947 if(fDebug > 10) Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
948 AliAODJet *genJet = (AliAODJet*) fListJetsGen->At(ig);
949 Double_t ptJetGen = genJet->Pt();
950 Double_t etaJetGen = genJet->Eta();
951
952 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
953 if((fJetEtaMin <= etaJetGen) && (etaJetGen <= fJetEtaMax)){
954 fhJetPtGenVsJetPtRec->Fill(ptJetRec, ptJetGen);
955 }
956 }//ig>=0
957 }//rec jet in eta acceptance
958 }//loop over reconstructed jets
959 }// # of rec jets >0
960 }//pointer MC generator jets
961 } //analyze generator level MC
962
8e103a56 963 //============= RECONSTRUCTED INCLUSIVE JETS ===============
ad869500 964
ad869500 965 Double_t etaJet = 0.0;
966 Double_t pTJet = 0.0;
967 Double_t areaJet = 0.0;
968 Double_t phiJet = 0.0;
7fc3d134 969 Int_t indexLeadingJet = -1;
970 Double_t pTLeadingJet = -10.0;
7fc3d134 971 Double_t areaLeadingJet = -10.0;
8e103a56 972
ad869500 973 for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
974 AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
975 if(!jet) continue;
976 etaJet = jet->Eta();
977 phiJet = jet->Phi();
978 pTJet = jet->Pt();
979 if(pTJet==0) continue;
980
981 if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
982 areaJet = jet->EffectiveAreaCharged();
7fc3d134 983
ad869500 984 //Jet Diagnostics---------------------------------
8e103a56 985 fhJetPhi->Fill(pTJet, RelativePhi(phiJet,0.0)); //phi -pi,pi
986 fhJetEta->Fill(pTJet, etaJet);
7fc3d134 987 //search for leading jet
988 if(pTJet > pTLeadingJet){
989 indexLeadingJet = ij;
990 pTLeadingJet = pTJet;
7fc3d134 991 areaLeadingJet = areaJet;
992 }
993
994 // raw spectra of INCLUSIVE jets
80ac66f6 995 //Centrality, pTjet, A
7fc3d134 996 Double_t fillraw[] = { centValue,
997 pTJet,
7fc3d134 998 areaJet
999 };
1000 fHJetPtRaw->Fill(fillraw);
8e103a56 1001 }
7fc3d134 1002
8e103a56 1003 if(indexLeadingJet > -1){
1004 // raw spectra of LEADING jets
1005 //Centrality, pTjet, A
1006 Double_t fillleading[] = { centValue,
1007 pTLeadingJet,
1008 areaLeadingJet
1009 };
1010 fHLeadingJetPtRaw->Fill(fillleading);
1011 }
1012
1013
1014 // =============== RECONSTRUCTED TRIGGER-JET PAIRS ================
1015 //Find Hadron trigger
1016 TList particleList; //list of tracks
1017 Int_t indexTrigg = GetListOfTracks(&particleList); //index of trigger hadron in Particle list
1018 if(fIsMC) FillEffHistos(&particleList, &particleListGen); //Fill efficiency histos
7fc3d134 1019
8e103a56 1020 //set ranges of the trigger loop
1021 Int_t ilow = (fHardest==0 || fHardest==1) ? indexTrigg : 0;
1022 Int_t ihigh = (fHardest==0 || fHardest==1) ? indexTrigg+1 : particleList.GetEntries();
7fc3d134 1023
8e103a56 1024 for(Int_t itrk= ilow; itrk< ihigh; itrk++){ //loop over possible trigger
1025 indexTrigg = itrk; //trigger hadron with pT >10 GeV
7fc3d134 1026
8e103a56 1027 if(indexTrigg < 0) continue;
1028
1029 AliVParticle *triggerHadron = (AliVParticle*) particleList.At(indexTrigg);
1030 if(!triggerHadron){
1031 PostData(1, fOutputList);
1032 return;
1033 }
1034 if(fHardest >= 2){
7d9fd386 1035 if(triggerHadron->Pt() < 6.0) continue; //all hadrons pt>6
8e103a56 1036 }
1037 if(TMath::Abs((Float_t) triggerHadron->Eta())> fTriggerEtaCut) continue;
7fc3d134 1038
8e103a56 1039 //Fill trigger histograms
1040 fh2Ntriggers->Fill(centValue,triggerHadron->Pt()); //trigger pT
1041 fhTriggerPhi->Fill(triggerHadron->Pt(),RelativePhi(triggerHadron->Phi(),0.0)); //phi -pi,pi
1042 fhTriggerEta->Fill(triggerHadron->Pt(),triggerHadron->Eta());
7fc3d134 1043
1044
8e103a56 1045 //---------- make trigger-jet pairs ---------
1046 Int_t injet4 = 0;
1047 Int_t injet = 0;
1048
1049 for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
1050 AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
1051 if(!jet) continue;
1052 etaJet = jet->Eta();
1053 phiJet = jet->Phi();
1054 pTJet = jet->Pt();
1055 if(pTJet==0) continue;
1056
1057 if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
1058 areaJet = jet->EffectiveAreaCharged();
1059 if(areaJet >= 0.07) injet++;
1060 if(areaJet >= 0.4) injet4++;
1061
1062 Double_t dphi = RelativePhi(triggerHadron->Phi(), phiJet);
1063 fhDphiTriggerJet->Fill(dphi); //Input
1064
1065 //Dphi versus jet pT
1066 //Centrality, Dphi=phiTrig-phiJet, pTjet, pTtrigg
1067 Double_t filldp[] = { centValue,
1068 dphi,
1069 pTJet,
1070 triggerHadron->Pt()
1071 };
1072 fHDphiVsJetPtAll->Fill(filldp);
1073
1074 // Select back to back trigger - jet pairs
1075 if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue;
1076 fhDphiTriggerJetAccept->Fill(dphi); //Accepted
1077
1078
1079 //Centrality, A, pTjet, pTtrigg
1080 Double_t fillspec[] = { centValue,
1081 areaJet,
1082 pTJet,
1083 triggerHadron->Pt()
1084 };
1085 fHJetSpec->Fill(fillspec);
1086 }//jet loop
1087
1088 //Fill Jet Density In the Event A>0.07
1089 if(injet>0){
1090 Double_t filldens[]={ (Double_t) particleList.GetEntries(),
ad869500 1091 injet/fkAcceptance,
1092 triggerHadron->Pt()
1093 };
8e103a56 1094 fHJetDensity->Fill(filldens);
1095 }
ad869500 1096
8e103a56 1097 //Fill Jet Density In the Event A>0.4
1098 if(injet4>0){
1099 Double_t filldens4[]={ (Double_t) particleList.GetEntries(),
1100 injet4/fkAcceptance,
1101 triggerHadron->Pt()
1102 };
1103 fHJetDensityA4->Fill(filldens4);
1104 }
ad869500 1105 }
1106
8e103a56 1107
ad869500 1108 PostData(1, fOutputList);
1109}
1110
1111//----------------------------------------------------------------------------
1112void AliAnalysisTaskJetCorePP::Terminate(const Option_t *)
1113{
1114 // Draw result to the screen
1115 // Called once at the end of the query
1116
1117 if(fDebug) printf("AliAnalysisTaskJetCorePP DONE\n");
1118 if(!GetOutputData(1)) return;
1119}
1120
1121
1122//----------------------------------------------------------------------------
1123Int_t AliAnalysisTaskJetCorePP::GetListOfTracks(TList *list){
1124 //Fill the list of accepted tracks (passed track cut)
1125 //return consecutive index of the hardest ch hadron in the list
1126 Int_t iCount = 0;
1127 AliAODEvent *aodevt = NULL;
1128
1129 if(!fESD) aodevt = fAODIn;
1130 else aodevt = fAODOut;
1131
1132 if(!aodevt) return -1;
1133
1134 Int_t index = -1; //index of the highest particle in the list
1135 Double_t ptmax = -10;
8e103a56 1136 Int_t triggers[200];
1137 Int_t ntriggers = 0; //index in triggers array
ad869500 1138
80ac66f6 1139 for(Int_t it = 0; it < aodevt->GetNumberOfTracks(); it++){
ad869500 1140 AliAODTrack *tr = aodevt->GetTrack(it);
1141
1142 //if((fFilterMask > 0) && !(tr->TestFilterBit(fFilterMask))) continue;
1143 if((fFilterMask > 0) && !(tr->IsHybridGlobalConstrainedGlobal())) continue;
1144 if(TMath::Abs((Float_t) tr->Eta()) > fTrackEtaCut) continue;
1145 if(tr->Pt() < fTrackLowPtCut) continue;
1146 list->Add(tr);
8e103a56 1147
1148 //Search trigger:
1149 if(fHardest>0){ //leading particle
1150 if(tr->Pt()>ptmax){
1151 ptmax = tr->Pt();
1152 index = iCount;
1153 }
1154 }
1155
1156 if(fHardest==0 && ntriggers<200){ //single inclusive
1157 if(fTriggerPtRangeLow <= tr->Pt() &&
1158 tr->Pt() < fTriggerPtRangeHigh){
1159 triggers[ntriggers] = iCount;
1160 ntriggers++;
1161 }
ad869500 1162 }
8e103a56 1163
ad869500 1164 iCount++;
1165 }
1166
8e103a56 1167 if(fHardest==0 && ntriggers>0){ //select random inclusive trigger
1168 Int_t rnd = fRandom->Integer(ntriggers); //0 to ntriggers-1
1169 index = triggers[rnd];
ad869500 1170 }
8e103a56 1171
1172 return index;
ad869500 1173}
1174
1175//----------------------------------------------------------------------------
80ac66f6 1176/*
ad869500 1177Double_t AliAnalysisTaskJetCorePP::GetBackgroundInPerpCone(Float_t jetR, Double_t jetPhi, Double_t jetEta, TList* trkList){
1178 //calculate sum of track pT in the cone perpendicular in phi to the jet
1179 //jetR = cone radius
1180 // jetPhi, jetEta = direction of the jet
1181 Int_t numberOfTrks = trkList->GetEntries();
1182 Double_t pTsum = 0.0;
1183 Double_t perpConePhi = jetPhi + TMath::Pi()/2;//perp cone w.r.t. jet in phi
1184 for(Int_t it=0; it<numberOfTrks; it++){
1185 AliVParticle *trk = (AliVParticle*) trkList->At(it);
1186 Double_t dphi = RelativePhi(perpConePhi,trk->Phi());
1187 Double_t deta = trk->Eta()-jetEta;
1188
1189 if( (dphi*dphi + deta*deta)< (jetR*jetR)){
1190 pTsum += trk->Pt();
1191 }
1192 }
1193
1194 return pTsum;
1195}
80ac66f6 1196*/
ad869500 1197//----------------------------------------------------------------------------
1198
1199Double_t AliAnalysisTaskJetCorePP::RelativePhi(Double_t mphi,Double_t vphi){
1200 //Get relative azimuthal angle of two particles -pi to pi
1201 if (vphi < -TMath::Pi()) vphi += TMath::TwoPi();
1202 else if (vphi > TMath::Pi()) vphi -= TMath::TwoPi();
1203
1204 if (mphi < -TMath::Pi()) mphi += TMath::TwoPi();
1205 else if (mphi > TMath::Pi()) mphi -= TMath::TwoPi();
1206
1207 Double_t dphi = mphi - vphi;
1208 if (dphi < -TMath::Pi()) dphi += TMath::TwoPi();
1209 else if (dphi > TMath::Pi()) dphi -= TMath::TwoPi();
1210
1211 return dphi;//dphi in [-Pi, Pi]
1212}
1213
1214
80ac66f6 1215//----------------------------------------------------------------------------
1216Bool_t AliAnalysisTaskJetCorePP::SelectMCGenTracks(AliVParticle *trk, TList *trkList, Double_t &ptLeading, Int_t &index, Int_t counter){
1217 //fill track efficiency denominator
1218 if(trk->Pt() < fTrackLowPtCut) return kFALSE;
1219 if(trk->Charge()==0) return kFALSE;
1220 if(TMath::Abs((Float_t) trk->Eta()) > fTrackEtaCut) return kFALSE;
1221 trkList->Add(trk);
1222 fhPtTrkTruePrimGen->Fill(trk->Pt(),trk->Eta()); //Efficiency denominator
8e103a56 1223
1224 //Search trigger:
1225 if(fHardest>0){ //leading particle
1226 if(ptLeading < trk->Pt()){
1227 index = counter;
1228 ptLeading = trk->Pt();
1229 }
80ac66f6 1230 }
8e103a56 1231
1232 if(fHardest==0){ //single inclusive
1233 index = -1;
1234 if(fTriggerPtRangeLow <= trk->Pt() &&
1235 trk->Pt() < fTriggerPtRangeHigh){
1236 index = counter;
1237 }
1238 }
1239
80ac66f6 1240 return kTRUE;
1241}
1242
1243//----------------------------------------------------------------------------
1244void AliAnalysisTaskJetCorePP::FillEffHistos(TList *recList, TList *genList){
1245 //fill track efficiency numerator
1246 if(!recList) return;
1247 if(!genList) return;
1248 Int_t nRec = recList->GetEntries();
1249 Int_t nGen = genList->GetEntries();
1250 for(Int_t ir=0; ir<nRec; ir++){
1251 AliAODTrack *trkRec = (AliAODTrack*) recList->At(ir);
1252 if(!trkRec) continue;
1253 Int_t recLabel = TMath::Abs(trkRec->GetLabel());
1254 Bool_t matched = kFALSE;
1255
1256 for(Int_t ig=0; ig<nGen; ig++){
1257 AliVParticle *trkGen = (AliVParticle*) genList->At(ig);
1258 if(!trkGen) continue;
1259 Int_t genLabel = TMath::Abs(trkGen->GetLabel());
1260 if(recLabel==genLabel){
1261 fhPtTrkTruePrimRec->Fill(trkGen->Pt(),trkGen->Eta());
1262 matched = kTRUE;
1263 break;
1264 }
1265 }
1266 if(!matched) fhPtTrkSecOrFakeRec->Fill(trkRec->Pt(),trkRec->Eta());
1267 }
1268
1269 return;
1270}
1271
1272