]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/AliAnalysisTaskJetCorePP.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Use Vzero multiplicity from online trigger
[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"
49dc2855 28#include "TClonesArray.h"
ad869500 29#include "TTree.h"
80ac66f6 30#include "TList.h"
ad869500 31#include "TMath.h"
32#include "TH1F.h"
33#include "TH1D.h"
34#include "TH2F.h"
35#include "TH3F.h"
36#include "THnSparse.h"
37#include "TCanvas.h"
80ac66f6 38#include "TArrayI.h"
ea603b64 39#include "TProfile.h"
40#include "TFile.h"
41#include "TKey.h"
8e103a56 42#include "TRandom3.h"
ad869500 43
44#include "AliLog.h"
45
46#include "AliAnalysisTask.h"
47#include "AliAnalysisManager.h"
48
49#include "AliVEvent.h"
50#include "AliESDEvent.h"
80ac66f6 51#include "AliMCEvent.h"
ad869500 52#include "AliESDInputHandler.h"
53#include "AliCentrality.h"
54#include "AliAnalysisHelperJetTasks.h"
55#include "AliInputEventHandler.h"
56#include "AliAODJetEventBackground.h"
ea603b64 57#include "AliGenPythiaEventHeader.h"
ad869500 58#include "AliAODMCParticle.h"
80ac66f6 59#include "AliMCParticle.h"
ad869500 60#include "AliAODEvent.h"
61#include "AliAODHandler.h"
62#include "AliAODJet.h"
49dc2855 63#include "AliVVertex.h"
0fe881ea 64#include "AliVTrack.h"
ad869500 65#include "AliAnalysisTaskJetCorePP.h"
29ac2c9b 66#include "AliHeader.h" //KF//
ad869500 67
68using std::cout;
69using std::endl;
70
71ClassImp(AliAnalysisTaskJetCorePP)
72
73//Filip Krizek 1st March 2013
74
75//---------------------------------------------------------------------
76AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP() :
77AliAnalysisTaskSE(),
78fESD(0x0),
79fAODIn(0x0),
80fAODOut(0x0),
81fAODExtension(0x0),
49dc2855 82fMcEvent(0x0),
83fMcHandler(0x0),
ad869500 84fJetBranchName(""),
255cb71c 85fJetBranchNameChargMC(""),
49dc2855 86fJetBranchNameKine(""),
255cb71c 87fJetBranchNameFullMC(""),
88fJetBranchNameBg(""),
89fJetBranchNameBgChargMC(""),
49dc2855 90fJetBranchNameBgKine(""),
7fc3d134 91fListJets(0x0),
80ac66f6 92fListJetsGen(0x0),
160dbdb8 93fListJetsGenFull(0x0),
255cb71c 94fListJetsBg(0x0),
95fListJetsBgGen(0x0),
ad869500 96fNonStdFile(""),
97fSystem(0), //pp=0 pPb=1
98fJetParamR(0.4),
8aa19603 99fBgJetParamR(0.3),
100fBgMaxJetPt(8.0),
101fBgConeR(0.4),
ad869500 102fOfflineTrgMask(AliVEvent::kAny),
103fMinContribVtx(1),
104fVtxZMin(-10.0),
105fVtxZMax(10.0),
106fFilterMask(0),
107fCentMin(0.0),
108fCentMax(100.0),
109fJetEtaMin(-0.5),
110fJetEtaMax(0.5),
111fTriggerEtaCut(0.9),
112fTrackEtaCut(0.9),
113fTrackLowPtCut(0.15),
49dc2855 114fUseExchContainer(0),
ad869500 115fOutputList(0x0),
116fHistEvtSelection(0x0),
117fh2Ntriggers(0x0),
118fHJetSpec(0x0),
255cb71c 119fHJetSpecSubUeMedian(0x0),
8aa19603 120fHJetSpecSubUeCone(0x0),
49dc2855 121fHJetPhiCorr(0x0),
122fHJetPhiCorrSubUeMedian(0x0),
123fHJetPhiCorrSubUeCone(0x0),
255cb71c 124fHJetUeMedian(0x0),
8aa19603 125fHJetUeCone(0x0),
126fHRhoUeMedianVsCone(0x0),
255cb71c 127//fHJetDensity(0x0),
128//fHJetDensityA4(0x0),
ad869500 129fhJetPhi(0x0),
130fhTriggerPhi(0x0),
131fhJetEta(0x0),
132fhTriggerEta(0x0),
133fhVertexZ(0x0),
134fhVertexZAccept(0x0),
135fhContribVtx(0x0),
136fhContribVtxAccept(0x0),
137fhDphiTriggerJet(0x0),
138fhDphiTriggerJetAccept(0x0),
139fhCentrality(0x0),
140fhCentralityAccept(0x0),
1e3ad34d 141fhNofMultipleTriggers(0x0),
0fe881ea 142fhNofMultipleTriggersCone(0x0),
30557bf5 143fhNofMultipleTriggersConeLow(0x0),
144fhNofMultipleTriggersConeHigh(0x0),
1e3ad34d 145fhDeltaRMultTriggers(0x0),
30557bf5 146fhDeltaPhiMultTriggers(0x0),
255cb71c 147//fHJetPtRaw(0x0),
148//fHLeadingJetPtRaw(0x0),
149//fHDphiVsJetPtAll(0x0),
80ac66f6 150fhJetPtGenVsJetPtRec(0x0),
255cb71c 151fhJetPtGenVsJetPtRecSubUeMedian(0x0),
8aa19603 152fhJetPtGenVsJetPtRecSubUeCone(0x0),
255cb71c 153fhJetPtGen(0x0),
154fhJetPtSubUeMedianGen(0x0),
8aa19603 155fhJetPtSubUeConeGen(0x0),
160dbdb8 156fhJetPtGenChargVsJetPtGenFull(0x0),
157fhJetPtGenFull(0x0),
80ac66f6 158fh2NtriggersGen(0x0),
159fHJetSpecGen(0x0),
255cb71c 160fHJetSpecSubUeMedianGen(0x0),
8aa19603 161fHJetSpecSubUeConeGen(0x0),
49dc2855 162fHJetPhiCorrGen(0x0),
163fHJetPhiCorrSubUeMedianGen(0x0),
164fHJetPhiCorrSubUeConeGen(0x0),
255cb71c 165fHJetUeMedianGen(0x0),
8aa19603 166fHJetUeConeGen(0x0),
80ac66f6 167fhPtTrkTruePrimRec(0x0),
168fhPtTrkTruePrimGen(0x0),
169fhPtTrkSecOrFakeRec(0x0),
8aa19603 170fHRhoUeMedianVsConeGen(0x0),
d1405a52
ML		 171fhEntriesToMedian(0x0),
172fhEntriesToMedianGen(0x0),
173fhCellAreaToMedian(0x0),
174fhCellAreaToMedianGen(0x0),
1e3ad34d 175fhNofMultipleTriggersGen(0x0),
0fe881ea 176fhNofMultipleTriggersConeGen(0x0),
30557bf5 177fhNofMultipleTriggersConeGenLow(0x0),
178fhNofMultipleTriggersConeGenHigh(0x0),
1e3ad34d 179fhDeltaRMultTriggersGen(0x0),
30557bf5 180fhDeltaPhiMultTriggersGen(0x0),
6c1489e1 181fhNofMultipleTriggersConeGenA(0x0),
30557bf5 182fhNofMultipleTriggersConeGenALow(0x0),
183fhNofMultipleTriggersConeGenAHigh(0x0),
6c1489e1 184fhNofMultipleTriggersConeGenB(0x0),
30557bf5 185fhNofMultipleTriggersConeGenBLow(0x0),
186fhNofMultipleTriggersConeGenBHigh(0x0),
255cb71c 187fIsChargedMC(0),
49dc2855 188fIsKine(0),
255cb71c 189fIsFullMC(0),
80ac66f6 190faGenIndex(0),
191faRecIndex(0),
7fc3d134 192fkAcceptance(2.0*TMath::Pi()*1.8),
6168afc3 193fkDeltaPhiCut(TMath::Pi()-0.8),
ea603b64 194fh1Xsec(0x0),
195fh1Trials(0x0),
196fh1AvgTrials(0x0),
197fh1PtHard(0x0),
198fh1PtHardNoW(0x0),
199fh1PtHardTrials(0x0),
8e103a56 200fAvgTrials(1),
201fHardest(0),
202fEventNumberRangeLow(0),
203fEventNumberRangeHigh(99),
204fTriggerPtRangeLow(0.0),
205fTriggerPtRangeHigh(10000.0),
8aa19603 206fFillRespMx(0),
207fRandom(0x0),
208fnTrials(1000),
d1405a52 209fJetFreeAreaFrac(0.5),
8aa19603 210fnPhi(9),
211fnEta(2),
212fEtaSize(0.7),
213fPhiSize(2*TMath::Pi()/fnPhi),
214fCellArea(fPhiSize*fEtaSize),
49dc2855 215fSafetyMargin(1.1),
216fDoubleBinning(kFALSE)
ad869500 217{
218 // default Constructor
ad869500 219}
220
221//---------------------------------------------------------------------
222
223AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP(const char *name) :
224AliAnalysisTaskSE(name),
225fESD(0x0),
226fAODIn(0x0),
227fAODOut(0x0),
228fAODExtension(0x0),
49dc2855 229fMcEvent(0x0),
230fMcHandler(0x0),
ad869500 231fJetBranchName(""),
255cb71c 232fJetBranchNameChargMC(""),
49dc2855 233fJetBranchNameKine(""),
255cb71c 234fJetBranchNameFullMC(""),
235fJetBranchNameBg(""),
236fJetBranchNameBgChargMC(""),
49dc2855 237fJetBranchNameBgKine(""),
7fc3d134 238fListJets(0x0),
80ac66f6 239fListJetsGen(0x0),
160dbdb8 240fListJetsGenFull(0x0),
255cb71c 241fListJetsBg(0x0),
242fListJetsBgGen(0x0),
ad869500 243fNonStdFile(""),
244fSystem(0), //pp=0 pPb=1
245fJetParamR(0.4),
8aa19603 246fBgJetParamR(0.3),
247fBgMaxJetPt(8.0),
248fBgConeR(0.4),
ad869500 249fOfflineTrgMask(AliVEvent::kAny),
250fMinContribVtx(1),
251fVtxZMin(-10.0),
252fVtxZMax(10.0),
253fFilterMask(0),
254fCentMin(0.0),
255fCentMax(100.0),
256fJetEtaMin(-0.5),
257fJetEtaMax(0.5),
258fTriggerEtaCut(0.9),
259fTrackEtaCut(0.9),
260fTrackLowPtCut(0.15),
49dc2855 261fUseExchContainer(0),
ad869500 262fOutputList(0x0),
263fHistEvtSelection(0x0),
264fh2Ntriggers(0x0),
265fHJetSpec(0x0),
255cb71c 266fHJetSpecSubUeMedian(0x0),
8aa19603 267fHJetSpecSubUeCone(0x0),
49dc2855 268fHJetPhiCorr(0x0),
269fHJetPhiCorrSubUeMedian(0x0),
270fHJetPhiCorrSubUeCone(0x0),
255cb71c 271fHJetUeMedian(0x0),
8aa19603 272fHJetUeCone(0x0),
273fHRhoUeMedianVsCone(0x0),
255cb71c 274//fHJetDensity(0x0),
275//fHJetDensityA4(0x0),
ad869500 276fhJetPhi(0x0),
277fhTriggerPhi(0x0),
278fhJetEta(0x0),
279fhTriggerEta(0x0),
280fhVertexZ(0x0),
281fhVertexZAccept(0x0),
282fhContribVtx(0x0),
283fhContribVtxAccept(0x0),
284fhDphiTriggerJet(0x0),
285fhDphiTriggerJetAccept(0x0),
286fhCentrality(0x0),
287fhCentralityAccept(0x0),
1e3ad34d 288fhNofMultipleTriggers(0x0),
0fe881ea 289fhNofMultipleTriggersCone(0x0),
30557bf5 290fhNofMultipleTriggersConeLow(0x0),
291fhNofMultipleTriggersConeHigh(0x0),
1e3ad34d 292fhDeltaRMultTriggers(0x0),
30557bf5 293fhDeltaPhiMultTriggers(0x0),
255cb71c 294//fHJetPtRaw(0x0),
295//fHLeadingJetPtRaw(0x0),
296//fHDphiVsJetPtAll(0x0),
80ac66f6 297fhJetPtGenVsJetPtRec(0x0),
255cb71c 298fhJetPtGenVsJetPtRecSubUeMedian(0x0),
8aa19603 299fhJetPtGenVsJetPtRecSubUeCone(0x0),
80ac66f6 300fhJetPtGen(0x0),
255cb71c 301fhJetPtSubUeMedianGen(0x0),
8aa19603 302fhJetPtSubUeConeGen(0x0),
160dbdb8 303fhJetPtGenChargVsJetPtGenFull(0x0),
304fhJetPtGenFull(0x0),
80ac66f6 305fh2NtriggersGen(0x0),
306fHJetSpecGen(0x0),
255cb71c 307fHJetSpecSubUeMedianGen(0x0),
8aa19603 308fHJetSpecSubUeConeGen(0x0),
49dc2855 309fHJetPhiCorrGen(0x0),
310fHJetPhiCorrSubUeMedianGen(0x0),
311fHJetPhiCorrSubUeConeGen(0x0),
255cb71c 312fHJetUeMedianGen(0x0),
8aa19603 313fHJetUeConeGen(0x0),
80ac66f6 314fhPtTrkTruePrimRec(0x0),
315fhPtTrkTruePrimGen(0x0),
316fhPtTrkSecOrFakeRec(0x0),
8aa19603 317fHRhoUeMedianVsConeGen(0x0),
d1405a52
ML		 318fhEntriesToMedian(0x0),
319fhEntriesToMedianGen(0x0),
320fhCellAreaToMedian(0x0),
321fhCellAreaToMedianGen(0x0),
1e3ad34d 322fhNofMultipleTriggersGen(0x0),
0fe881ea 323fhNofMultipleTriggersConeGen(0x0),
30557bf5 324fhNofMultipleTriggersConeGenLow(0x0),
325fhNofMultipleTriggersConeGenHigh(0x0),
1e3ad34d 326fhDeltaRMultTriggersGen(0x0),
30557bf5 327fhDeltaPhiMultTriggersGen(0x0),
6c1489e1 328fhNofMultipleTriggersConeGenA(0x0),
30557bf5 329fhNofMultipleTriggersConeGenALow(0x0),
330fhNofMultipleTriggersConeGenAHigh(0x0),
6c1489e1 331fhNofMultipleTriggersConeGenB(0x0),
30557bf5 332fhNofMultipleTriggersConeGenBLow(0x0),
333fhNofMultipleTriggersConeGenBHigh(0x0),
255cb71c 334fIsChargedMC(0),
49dc2855 335fIsKine(0),
255cb71c 336fIsFullMC(0),
80ac66f6 337faGenIndex(0),
338faRecIndex(0),
7fc3d134 339fkAcceptance(2.0*TMath::Pi()*1.8),
6168afc3 340fkDeltaPhiCut(TMath::Pi()-0.8),
ea603b64 341fh1Xsec(0x0),
342fh1Trials(0x0),
343fh1AvgTrials(0x0),
344fh1PtHard(0x0),
345fh1PtHardNoW(0x0),
346fh1PtHardTrials(0x0),
8e103a56 347fAvgTrials(1),
348fHardest(0),
349fEventNumberRangeLow(0),
350fEventNumberRangeHigh(99),
351fTriggerPtRangeLow(0.0),
352fTriggerPtRangeHigh(10000.0),
8aa19603 353fFillRespMx(0),
354fRandom(0x0),
355fnTrials(1000),
d1405a52 356fJetFreeAreaFrac(0.5),
8aa19603 357fnPhi(9),
358fnEta(2),
359fEtaSize(0.7),
360fPhiSize(2*TMath::Pi()/fnPhi),
361fCellArea(fPhiSize*fEtaSize),
49dc2855 362fSafetyMargin(1.1),
363fDoubleBinning(kFALSE)
ad869500 364{
7fc3d134 365// Constructor
ad869500 366 DefineOutput(1, TList::Class());
49dc2855 367
368 TString dummy(name);
369 if(dummy.Contains("KINE")){
370 DefineInput(1, TClonesArray::Class());
3d1fac8a 371 DefineInput(2, TClonesArray::Class());
49dc2855 372 }
ad869500 373}
374
375//--------------------------------------------------------------
376AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP(const AliAnalysisTaskJetCorePP& a):
377AliAnalysisTaskSE(a.GetName()),
378fESD(a.fESD),
379fAODIn(a.fAODIn),
380fAODOut(a.fAODOut),
381fAODExtension(a.fAODExtension),
49dc2855 382fMcEvent(a.fMcEvent),
383fMcHandler(a.fMcHandler),
ad869500 384fJetBranchName(a.fJetBranchName),
255cb71c 385fJetBranchNameChargMC(a.fJetBranchNameChargMC),
49dc2855 386fJetBranchNameKine(a.fJetBranchNameKine),
255cb71c 387fJetBranchNameFullMC(a.fJetBranchNameFullMC),
388fJetBranchNameBg(a.fJetBranchNameBg),
389fJetBranchNameBgChargMC(a.fJetBranchNameBgChargMC),
49dc2855 390fJetBranchNameBgKine(a.fJetBranchNameBgKine),
ad869500 391fListJets(a.fListJets),
80ac66f6 392fListJetsGen(a.fListJetsGen),
160dbdb8 393fListJetsGenFull(a.fListJetsGenFull),
255cb71c 394fListJetsBg(a.fListJetsBg),
395fListJetsBgGen(a.fListJetsBgGen),
ad869500 396fNonStdFile(a.fNonStdFile),
397fSystem(a.fSystem),
398fJetParamR(a.fJetParamR),
8aa19603 399fBgJetParamR(a.fBgJetParamR),
400fBgMaxJetPt(a.fBgMaxJetPt),
401fBgConeR(a.fBgConeR),
ad869500 402fOfflineTrgMask(a.fOfflineTrgMask),
403fMinContribVtx(a.fMinContribVtx),
404fVtxZMin(a.fVtxZMin),
405fVtxZMax(a.fVtxZMax),
406fFilterMask(a.fFilterMask),
407fCentMin(a.fCentMin),
408fCentMax(a.fCentMax),
409fJetEtaMin(a.fJetEtaMin),
410fJetEtaMax(a.fJetEtaMax),
411fTriggerEtaCut(a.fTriggerEtaCut),
412fTrackEtaCut(a.fTrackEtaCut),
413fTrackLowPtCut(a.fTrackLowPtCut),
49dc2855 414fUseExchContainer(a.fUseExchContainer),
ad869500 415fOutputList(a.fOutputList),
416fHistEvtSelection(a.fHistEvtSelection),
417fh2Ntriggers(a.fh2Ntriggers),
418fHJetSpec(a.fHJetSpec),
255cb71c 419fHJetSpecSubUeMedian(a.fHJetSpecSubUeMedian),
8aa19603 420fHJetSpecSubUeCone(a.fHJetSpecSubUeCone),
49dc2855 421fHJetPhiCorr(a.fHJetPhiCorr),
422fHJetPhiCorrSubUeMedian(a.fHJetPhiCorrSubUeMedian),
423fHJetPhiCorrSubUeCone(a.fHJetPhiCorrSubUeCone),
255cb71c 424fHJetUeMedian(a.fHJetUeMedian),
8aa19603 425fHJetUeCone(a.fHJetUeCone),
426fHRhoUeMedianVsCone(a.fHRhoUeMedianVsCone),
255cb71c 427//fHJetDensity(a.fHJetDensity),
428//fHJetDensityA4(a.fHJetDensityA4),
ad869500 429fhJetPhi(a.fhJetPhi),
430fhTriggerPhi(a.fhTriggerPhi),
431fhJetEta(a.fhJetEta),
432fhTriggerEta(a.fhTriggerEta),
433fhVertexZ(a.fhVertexZ),
434fhVertexZAccept(a.fhVertexZAccept),
435fhContribVtx(a.fhContribVtx),
436fhContribVtxAccept(a.fhContribVtxAccept),
437fhDphiTriggerJet(a.fhDphiTriggerJet),
438fhDphiTriggerJetAccept(a.fhDphiTriggerJetAccept),
439fhCentrality(a.fhCentrality),
440fhCentralityAccept(a.fhCentralityAccept),
1e3ad34d 441fhNofMultipleTriggers(a.fhNofMultipleTriggers),
0fe881ea 442fhNofMultipleTriggersCone(a.fhNofMultipleTriggersCone),
30557bf5 443fhNofMultipleTriggersConeLow(a.fhNofMultipleTriggersConeLow),
444fhNofMultipleTriggersConeHigh(a.fhNofMultipleTriggersConeHigh),
1e3ad34d 445fhDeltaRMultTriggers(a.fhDeltaRMultTriggers),
30557bf5 446fhDeltaPhiMultTriggers(a.fhDeltaPhiMultTriggers),
255cb71c 447//fHJetPtRaw(a.fHJetPtRaw),
448//fHLeadingJetPtRaw(a.fHLeadingJetPtRaw),
449//fHDphiVsJetPtAll(a.fHDphiVsJetPtAll),
80ac66f6 450fhJetPtGenVsJetPtRec(a.fhJetPtGenVsJetPtRec),
255cb71c 451fhJetPtGenVsJetPtRecSubUeMedian(a.fhJetPtGenVsJetPtRecSubUeMedian),
8aa19603 452fhJetPtGenVsJetPtRecSubUeCone(a.fhJetPtGenVsJetPtRecSubUeCone),
80ac66f6 453fhJetPtGen(a.fhJetPtGen),
255cb71c 454fhJetPtSubUeMedianGen(a.fhJetPtSubUeMedianGen),
8aa19603 455fhJetPtSubUeConeGen(a.fhJetPtSubUeConeGen),
160dbdb8 456fhJetPtGenChargVsJetPtGenFull(a.fhJetPtGenChargVsJetPtGenFull),
457fhJetPtGenFull(a.fhJetPtGenFull),
80ac66f6 458fh2NtriggersGen(a.fh2NtriggersGen),
459fHJetSpecGen(a.fHJetSpecGen),
255cb71c 460fHJetSpecSubUeMedianGen(a.fHJetSpecSubUeMedianGen),
8aa19603 461fHJetSpecSubUeConeGen(a.fHJetSpecSubUeConeGen),
49dc2855 462fHJetPhiCorrGen(a.fHJetPhiCorrGen),
463fHJetPhiCorrSubUeMedianGen(a.fHJetPhiCorrSubUeMedianGen),
464fHJetPhiCorrSubUeConeGen(a.fHJetPhiCorrSubUeConeGen),
255cb71c 465fHJetUeMedianGen(a.fHJetUeMedianGen),
8aa19603 466fHJetUeConeGen(a.fHJetUeConeGen),
80ac66f6 467fhPtTrkTruePrimRec(a.fhPtTrkTruePrimRec),
468fhPtTrkTruePrimGen(a.fhPtTrkTruePrimGen),
469fhPtTrkSecOrFakeRec(a.fhPtTrkSecOrFakeRec),
8aa19603 470fHRhoUeMedianVsConeGen(a.fHRhoUeMedianVsConeGen),
d1405a52
ML		 471fhEntriesToMedian(a.fhEntriesToMedian),
472fhEntriesToMedianGen(a.fhEntriesToMedianGen),
473fhCellAreaToMedian(a.fhCellAreaToMedian),
474fhCellAreaToMedianGen(a.fhCellAreaToMedianGen),
1e3ad34d 475fhNofMultipleTriggersGen(a.fhNofMultipleTriggersGen),
0fe881ea 476fhNofMultipleTriggersConeGen(a.fhNofMultipleTriggersConeGen),
30557bf5 477fhNofMultipleTriggersConeGenLow(a.fhNofMultipleTriggersConeGenLow),
478fhNofMultipleTriggersConeGenHigh(a.fhNofMultipleTriggersConeGenHigh),
1e3ad34d 479fhDeltaRMultTriggersGen(a.fhDeltaRMultTriggersGen),
30557bf5 480fhDeltaPhiMultTriggersGen(a.fhDeltaPhiMultTriggersGen),
6c1489e1 481fhNofMultipleTriggersConeGenA(a.fhNofMultipleTriggersConeGenA),
30557bf5 482fhNofMultipleTriggersConeGenALow(a.fhNofMultipleTriggersConeGenALow),
483fhNofMultipleTriggersConeGenAHigh(a.fhNofMultipleTriggersConeGenAHigh),
6c1489e1 484fhNofMultipleTriggersConeGenB(a.fhNofMultipleTriggersConeGenB),
30557bf5 485fhNofMultipleTriggersConeGenBLow(a.fhNofMultipleTriggersConeGenBLow),
486fhNofMultipleTriggersConeGenBHigh(a.fhNofMultipleTriggersConeGenBHigh),
255cb71c 487fIsChargedMC(a.fIsChargedMC),
49dc2855 488fIsKine(a.fIsKine),
255cb71c 489fIsFullMC(a.fIsFullMC),
80ac66f6 490faGenIndex(a.faGenIndex),
491faRecIndex(a.faRecIndex),
7fc3d134 492fkAcceptance(a.fkAcceptance),
ea603b64 493fkDeltaPhiCut(a.fkDeltaPhiCut),
494fh1Xsec(a.fh1Xsec),
495fh1Trials(a.fh1Trials),
496fh1AvgTrials(a.fh1AvgTrials),
497fh1PtHard(a.fh1PtHard),
498fh1PtHardNoW(a.fh1PtHardNoW),
499fh1PtHardTrials(a.fh1PtHardTrials),
8e103a56 500fAvgTrials(a.fAvgTrials),
501fHardest(a.fHardest),
502fEventNumberRangeLow(a.fEventNumberRangeLow),
503fEventNumberRangeHigh(a.fEventNumberRangeHigh),
504fTriggerPtRangeLow(a.fTriggerPtRangeLow),
43457761 505fTriggerPtRangeHigh(a.fTriggerPtRangeHigh),
8aa19603 506fFillRespMx(a.fFillRespMx),
507fRandom(a.fRandom),
508fnTrials(a.fnTrials),
d1405a52 509fJetFreeAreaFrac(a.fJetFreeAreaFrac),
8aa19603 510fnPhi(a.fnPhi),
511fnEta(a.fnEta),
512fEtaSize(a.fEtaSize),
513fPhiSize(a.fPhiSize),
514fCellArea(a.fCellArea),
49dc2855 515fSafetyMargin(a.fSafetyMargin),
516fDoubleBinning(a.fDoubleBinning)
ad869500 517{
518 //Copy Constructor
43457761 519 fRandom->SetSeed(0);
ad869500 520}
521//--------------------------------------------------------------
522
523AliAnalysisTaskJetCorePP& AliAnalysisTaskJetCorePP::operator = (const AliAnalysisTaskJetCorePP& a){
524 // assignment operator
525 this->~AliAnalysisTaskJetCorePP();
526 new(this) AliAnalysisTaskJetCorePP(a);
527 return *this;
528}
529//--------------------------------------------------------------
530
531AliAnalysisTaskJetCorePP::~AliAnalysisTaskJetCorePP()
532{
533 //Destructor
534 delete fListJets;
80ac66f6 535 delete fListJetsGen;
160dbdb8 536 delete fListJetsGenFull;
255cb71c 537 delete fListJetsBg;
538 delete fListJetsBgGen;
ad869500 539 delete fOutputList; // ????
8e103a56 540 delete fRandom;
ad869500 541}
542
543//--------------------------------------------------------------
544
ea603b64 545
546Bool_t AliAnalysisTaskJetCorePP::Notify()
547{
548 //Implemented Notify() to read the cross sections
549 //and number of trials from pyxsec.root
550 //inspired by AliAnalysisTaskJetSpectrum2::Notify()
49dc2855 551 if(!(fIsChargedMC || fIsKine)) return kFALSE;
3d1fac8a 552 Float_t xsection = 0;
553 Float_t trials = 1;
554 fAvgTrials = 1;
ea603b64 555
3d1fac8a 556 if(fIsChargedMC){
557 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
558 if(!fESD){
559 if(fDebug>1) AliError("ESD not available");
560 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
561 }
ea603b64 562
3d1fac8a 563 fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
ea603b64 564
565
3d1fac8a 566 if(fNonStdFile.Length()!=0){
567 // case that we have an AOD extension we can fetch the jets from the extended output
568 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
569 fAODExtension = aodH ? aodH->GetExtension(fNonStdFile.Data()) : 0;
570 if(!fAODExtension){
571 if(fDebug>1) Printf("AODExtension found for %s",fNonStdFile.Data());
572 }
573 }
ea603b64 574
3d1fac8a 575 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
ea603b64 576
3d1fac8a 577 if(tree){
578 TFile *curfile = tree->GetCurrentFile();
579 if(!curfile) {
580 Error("Notify","No current file");
581 return kFALSE;
582 }
583 if(!fh1Xsec || !fh1Trials){
584 Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
585 return kFALSE;
586 }
587 AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,trials);
588 fh1Xsec->Fill("<#sigma>",xsection);
589 // construct a poor man average trials
590 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
591 if(trials>=nEntries && nEntries>0.) fAvgTrials = trials/nEntries;
592 fh1Trials->Fill("#sum{ntrials}",trials);
593 }
594
595 if(fDebug)Printf("Reading File %s",fInputHandler->GetTree()->GetCurrentFile()->GetName());
596 }
ea603b64 597
ea603b64 598
599 return kTRUE;
600}
601//--------------------------------------------------------------
602
ad869500 603void AliAnalysisTaskJetCorePP::Init()
604{
605 // check for jet branches
49dc2855 606 if(fJetBranchNameKine.Length()==0){
607 if(!strlen(fJetBranchName.Data())){
608 AliError("Jet branch name not set.");
609 }
ad869500 610 }
49dc2855 611 fMcHandler = dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
612
613
ad869500 614}
615
616//--------------------------------------------------------------
617
618void AliAnalysisTaskJetCorePP::UserCreateOutputObjects()
619{
8aa19603 620 // Create histograms and initilize variables
621 fSafetyMargin = fBgConeR*fBgConeR /(fBgJetParamR*fBgJetParamR);
49dc2855 622
ad869500 623 // Called once
255cb71c 624 fListJets = new TList(); //reconstructed level antikt jets
8aa19603 625 fListJetsBg = new TList(); //reconstructed jets to be removed from bg
80ac66f6 626
255cb71c 627 fIsChargedMC = (fJetBranchNameChargMC.Length()>0) ? kTRUE : kFALSE;
49dc2855 628 fIsKine = (fJetBranchNameKine.Length()>0) ? kTRUE : kFALSE;
255cb71c 629 fIsFullMC = (fJetBranchNameFullMC.Length()>0) ? kTRUE : kFALSE;
80ac66f6 630
8e103a56 631 fRandom = new TRandom3(0);
632
49dc2855 633 if(fIsChargedMC || fIsKine){ //full MC or pure kine
255cb71c 634 fListJetsGen = new TList(); //generator level charged antikt jets
8aa19603 635 fListJetsBgGen = new TList(); //generator level jets to be removed from bg
255cb71c 636
637 if(fIsFullMC)
160dbdb8 638 fListJetsGenFull = new TList(); //generator level full jets
639 }
ad869500 640 OpenFile(1);
641 if(!fOutputList) fOutputList = new TList();
642 fOutputList->SetOwner(kTRUE);
643
644 Bool_t oldStatus = TH1::AddDirectoryStatus();
645 TH1::AddDirectory(kFALSE);
646
647 fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
648 fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
649 fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
650 fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
651 fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
652 fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
8e103a56 653 fHistEvtSelection->GetXaxis()->SetBinLabel(6,"event number (rejected)");
ad869500 654
655 fOutputList->Add(fHistEvtSelection);
656
657 Int_t nBinsCentrality = (fSystem==0) ? 1 : 10; // pp=1 else 10
80ac66f6 658 //trigger pt spectrum (reconstructed)
ad869500 659 fh2Ntriggers = new TH2F("fh2Ntriggers","# of triggers",
660 nBinsCentrality,0.0,100.0,50,0.0,50.0);
49dc2855 661 if(!fIsKine) fOutputList->Add(fh2Ntriggers);
662
663 Int_t bw = (fDoubleBinning==0) ? 1 : 2; //make larger bin width
ad869500 664
6168afc3 665 //Centrality, A, pTjet, pTtrigg, dphi
666 const Int_t dimSpec = 5;
49dc2855 667 const Int_t nBinsSpec[dimSpec] = {nBinsCentrality, 50, bw*110, 50, TMath::Nint(10*(TMath::Pi()-fkDeltaPhiCut))};
668 const Double_t lowBinSpec[dimSpec] = {0.0, 0.0, -20, 0.0, fkDeltaPhiCut};
669 const Double_t hiBinSpec[dimSpec] = {100.0, 1.0, 200.0,50.0, TMath::Pi()};
7fc3d134 670 fHJetSpec = new THnSparseF("fHJetSpec",
6168afc3 671 "Recoil jet spectrum [cent,A,pTjet,pTtrig,dphi]",
7fc3d134 672 dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
49dc2855 673 if(!fIsKine) fOutputList->Add(fHJetSpec);
ad869500 674
255cb71c 675 //background estimated as median of kT jets
676 fHJetSpecSubUeMedian = (THnSparseF*) fHJetSpec->Clone("fHJetSpecSubUeMedian");
677 fHJetSpecSubUeMedian->SetTitle("Recoil jet spectrum [cent,A,pTjet-pTUe,pTtrig,dphi]");
49dc2855 678 if(!fIsKine) fOutputList->Add(fHJetSpecSubUeMedian);
8aa19603 679 //background estimated as weighted median of kT jets ala Cone
680 fHJetSpecSubUeCone = (THnSparseF*) fHJetSpec->Clone("fHJetSpecSubUeCone");
681 fHJetSpecSubUeCone->SetTitle("Recoil jet spectrum [cent,A,pTjet-pTUe,pTtrig,dphi]");
49dc2855 682 if(!fIsKine) fOutputList->Add(fHJetSpecSubUeCone);
683
255cb71c 684
49dc2855 685 //A, pTjet, pTtrigg, dphi
686 const Int_t dimCor = 4;
687 const Int_t nBinsCor[dimCor] = {50, 110, 50, 100};
688 const Double_t lowBinCor[dimCor] = {0.0, -20, 0.0, -0.5*TMath::Pi()};
689 const Double_t hiBinCor[dimCor] = {1.0, 200, 50.0, 1.5*TMath::Pi()};
690 fHJetPhiCorr = new THnSparseF("fHJetPhiCorr",
691 "Recoil jet spectrum [A,pTjet,pTtrig,dphi]",
692 dimCor,nBinsCor,lowBinCor,hiBinCor);
693
694 if(!fIsKine) fOutputList->Add(fHJetPhiCorr); // Dphi distribution jet-triger
695
696 //background estimated as median of kT jets
697 fHJetPhiCorrSubUeMedian = (THnSparseF*) fHJetPhiCorr->Clone("fHJetPhiCorrSubUeMedian");
698 fHJetPhiCorrSubUeMedian->SetTitle("Recoil jet spectrum [A,pTjet-pTUe,pTtrig,dphi]");
699 if(!fIsKine) fOutputList->Add(fHJetPhiCorrSubUeMedian);
700 //background estimated as weighted median of kT jets ala Cone
701 fHJetPhiCorrSubUeCone = (THnSparseF*) fHJetPhiCorr->Clone("fHJetPhiCorrSubUeCone");
702 fHJetPhiCorrSubUeCone->SetTitle("Recoil jet spectrum [A,pTjet-pTUe,pTtrig,dphi]");
703 if(!fIsKine) fOutputList->Add(fHJetPhiCorrSubUeCone);
255cb71c 704
705
ad869500 706 //------------------- HISTOS FOR DIAGNOSTIC ----------------------
255cb71c 707 //A, pTjet, pTjet-pTUe, pTUe, rhoUe bg estimate from kT median
708 const Int_t dimSpecMed = 5;
1afc79e9 709 const Int_t nBinsSpecMed[dimSpecMed] = {25, 50, 60, 20, 20};
8aa19603 710 const Double_t lowBinSpecMed[dimSpecMed] = {0.0, 0.0, -20.0, 0.0, 0.0};
1afc79e9 711 const Double_t hiBinSpecMed[dimSpecMed] = {1.0, 100.0, 100.0, 10.0, 20.0};
255cb71c 712 fHJetUeMedian = new THnSparseF("fHJetUeMedian",
713 "Recoil jet spectrum [A,pTjet,pTjet-pTUe, pTUe, rhoUe]",
714 dimSpecMed, nBinsSpecMed, lowBinSpecMed, hiBinSpecMed);
49dc2855 715 if(!fIsKine) fOutputList->Add(fHJetUeMedian);
255cb71c 716
8aa19603 717 //A, pTjet, pTjet-pTUe, pTUe, rhoUe bg estimate from kT median Cone
718 fHJetUeCone = (THnSparseF*) fHJetUeMedian->Clone("fHJetUeCone");
719 fHJetUeCone->SetTitle("Recoil jet spectrum [A,pTjet,pTjet-pTUe, pTUe, rhoUe]");
49dc2855 720 if(!fIsKine) fOutputList->Add(fHJetUeCone);
255cb71c 721
722 //rho bacground reconstructed data
723 const Int_t dimRho = 2;
1afc79e9 724 const Int_t nBinsRho[dimRho] = {50 , 50};
255cb71c 725 const Double_t lowBinRho[dimRho] = {0.0 , 0.0};
726 const Double_t hiBinRho[dimRho] = {20.0 , 20.0};
727
8aa19603 728 fHRhoUeMedianVsCone = new THnSparseF("hRhoUeMedianVsCone","[Rho Cone, Rho Median]",
255cb71c 729 dimRho, nBinsRho, lowBinRho, hiBinRho);
49dc2855 730 if(!fIsKine) fOutputList->Add(fHRhoUeMedianVsCone);
255cb71c 731
ad869500 732 //Jet number density histos [Trk Mult, jet density, pT trigger]
255cb71c 733 /*const Int_t dimJetDens = 3;
80ac66f6 734 const Int_t nBinsJetDens[dimJetDens] = {100, 100, 10};
ad869500 735 const Double_t lowBinJetDens[dimJetDens] = {0.0, 0.0, 0.0};
80ac66f6 736 const Double_t hiBinJetDens[dimJetDens] = {500.0, 5.0, 50.0};
ad869500 737
738 fHJetDensity = new THnSparseF("fHJetDensity","Jet dens vs trk mult A>0.07",
739 dimJetDens,nBinsJetDens,lowBinJetDens,hiBinJetDens);
740
741 fHJetDensityA4 =new THnSparseF("fHJetDensityA4","Jet dens vs trk mult A>0.4",
742 dimJetDens,nBinsJetDens,lowBinJetDens,hiBinJetDens);
743
ad869500 744 fOutputList->Add(fHJetDensity);
745 fOutputList->Add(fHJetDensityA4);
255cb71c 746 */
ad869500 747
748 //inclusive azimuthal and pseudorapidity histograms
7fc3d134 749 fhJetPhi = new TH2D("fhJetPhi","Azim dist jets vs pTjet",
750 50, 0, 100, 50,-TMath::Pi(),TMath::Pi());
751 fhTriggerPhi= new TH2D("fhTriggerPhi","azim dist trig had vs pTtrigg",
6168afc3 752 50, 0, 50, 50,-TMath::Pi(),TMath::Pi());
7fc3d134 753 fhJetEta = new TH2D("fhJetEta","Eta dist jets vs pTjet",
754 50,0, 100, 40,-0.9,0.9);
755 fhTriggerEta = new TH2D("fhTriggerEta","Eta dist trig had vs pTtrigg",
6168afc3 756 50, 0, 50, 40,-0.9,0.9);
7fc3d134 757
ad869500 758 fhVertexZ = new TH1D("fhVertexZ","z vertex",40,-20,20);
759 fhVertexZAccept = new TH1D("fhVertexZAccept","z vertex after cut",40,-20,20);
760 fhContribVtx = new TH1D("fhContribVtx","contrib to vtx",200,0,200);
761 fhContribVtxAccept = new TH1D("fhContribVtxAccept","contrib to vtx after cut",200,0,200);
762 fhDphiTriggerJet = new TH1D("fhDphiTriggerJet","Deltaphi trig-jet",50, -TMath::Pi(),TMath::Pi());
763 fhDphiTriggerJetAccept = new TH1D("fhDphiTriggerJetAccept","Deltaphi trig-jet after cut",50, -TMath::Pi(),TMath::Pi());
764 fhCentrality = new TH1D("fhCentrality","Centrality",20,0,100);
765 fhCentralityAccept = new TH1D("fhCentralityAccept","Centrality after cut",20,0,100);
d1405a52 766 fhEntriesToMedian = new TH1D("fhEntriesToMedian","fhEntriesToMedian",30,0,30);
49dc2855 767 fhCellAreaToMedian = new TH1D("fhCellAreaToMedian", "fhCellAreaToMedian", 75,0,1.5);
ad869500 768
30557bf5 769 fhNofMultipleTriggers = new TH1D("fhNofMultipleTriggers","fhNofMultipleTriggers",20,0,20);
770 fhNofMultipleTriggersCone = new TH1D("fhNofMultipleTriggersCone","fhNofMultipleTriggersCone R<0.4",20,0,20);
771 fhNofMultipleTriggersConeLow = (TH1D*) fhNofMultipleTriggersCone->Clone("fhNofMultipleTriggersConepTa15to20");
772 fhNofMultipleTriggersConeHigh = (TH1D*) fhNofMultipleTriggersCone->Clone("fhNofMultipleTriggersConepTa20to50");
1e3ad34d 773 fhDeltaRMultTriggers = new TH1D("fhDeltaRMultTriggers","fhDeltaRMultTriggers", 100,0,4);
30557bf5 774 fhDeltaPhiMultTriggers = new TH1D("fhDeltaPhiMultTriggers","fhDeltaPhiRultTriggers", 100,-TMath::Pi(),TMath::Pi());
1e3ad34d 775
49dc2855 776 if(!fIsKine){
777 fOutputList->Add(fhJetPhi);
778 fOutputList->Add(fhTriggerPhi);
779 fOutputList->Add(fhJetEta);
780 fOutputList->Add(fhTriggerEta);
781 }
ad869500 782 fOutputList->Add(fhVertexZ);
49dc2855 783 fOutputList->Add(fhVertexZAccept);
784 if(!fIsKine){
785 fOutputList->Add(fhContribVtx);
786 fOutputList->Add(fhContribVtxAccept);
787 fOutputList->Add(fhDphiTriggerJet);
788 fOutputList->Add(fhDphiTriggerJetAccept);
789 fOutputList->Add(fhCentrality);
790 fOutputList->Add(fhCentralityAccept);
791 fOutputList->Add(fhEntriesToMedian);
792 fOutputList->Add(fhCellAreaToMedian);
1e3ad34d 793 fOutputList->Add(fhNofMultipleTriggers);
0fe881ea 794 fOutputList->Add(fhNofMultipleTriggersCone);
30557bf5 795 fOutputList->Add(fhNofMultipleTriggersConeLow);
796 fOutputList->Add(fhNofMultipleTriggersConeHigh);
1e3ad34d 797 fOutputList->Add(fhDeltaRMultTriggers);
30557bf5 798 fOutputList->Add(fhDeltaPhiMultTriggers);
49dc2855 799 }
7fc3d134 800 // raw spectra of INCLUSIVE jets
80ac66f6 801 //Centrality, pTjet, A
255cb71c 802 /*const Int_t dimRaw = 3;
80ac66f6 803 const Int_t nBinsRaw[dimRaw] = {nBinsCentrality, 50, 100};
804 const Double_t lowBinRaw[dimRaw] = {0.0, 0.0, 0.0};
805 const Double_t hiBinRaw[dimRaw] = {100.0, 100, 1.0};
7fc3d134 806 fHJetPtRaw = new THnSparseF("fHJetPtRaw",
80ac66f6 807 "Incl. jet spectrum [cent,pTjet,A]",
7fc3d134 808 dimRaw,nBinsRaw,lowBinRaw,hiBinRaw);
809 fOutputList->Add(fHJetPtRaw);
810
811 // raw spectra of LEADING jets
80ac66f6 812 //Centrality, pTjet, A
7fc3d134 813 fHLeadingJetPtRaw = new THnSparseF("fHLeadingJetPtRaw",
80ac66f6 814 "Leading jet spectrum [cent,pTjet,A]",
7fc3d134 815 dimRaw,nBinsRaw,lowBinRaw,hiBinRaw);
816 fOutputList->Add(fHLeadingJetPtRaw);
817
818 // Dphi versus pT jet
819 //Centrality, Dphi=phiTrig-phiJet, pTjet, pTtrigg
820 const Int_t dimDp = 4;
821 const Int_t nBinsDp[dimDp] = {nBinsCentrality, 50, 50, 50};
822 const Double_t lowBinDp[dimDp] = {0.0, -TMath::Pi(), 0.0, 0.0};
823 const Double_t hiBinDp[dimDp] = {100.0, TMath::Pi(), 100.0, 100.0};
824 fHDphiVsJetPtAll = new THnSparseF("fHDphiVsJetPtAll",
825 "Dphi vs jet pT [cent,Dphi,pTjet,pTtrigg]",
826 dimDp,nBinsDp,lowBinDp,hiBinDp);
827 fOutputList->Add(fHDphiVsJetPtAll);
255cb71c 828 */
7fc3d134 829
80ac66f6 830 //analyze MC generator level
49dc2855 831 if(fIsChargedMC || fIsKine){
8aa19603 832 if(fFillRespMx){
833 //Fill response matrix only once
49dc2855 834 fhJetPtGenVsJetPtRec = new TH2D("fhJetPtGenVsJetPtRec","JetPtGenVsJetPtRec", bw*100,0,200, bw*100,0,200);
8aa19603 835 fOutputList->Add(fhJetPtGenVsJetPtRec); //gen MC charg jet pt spectrum versus rec charged jet pt spectrum
836 //....
49dc2855 837 fhJetPtGenVsJetPtRecSubUeMedian = new TH2D("fhJetPtGenVsJetPtRecSubUeMedian","fhJetPtGenVsJetPtRecSubUeMedian", bw*110,-20,200, bw*110,-20,200);
8aa19603 838 fOutputList->Add(fhJetPtGenVsJetPtRecSubUeMedian); // with kT median bg subtr
255cb71c 839 //....
8aa19603 840 fhJetPtGenVsJetPtRecSubUeCone=(TH2D*)fhJetPtGenVsJetPtRecSubUeMedian->Clone("fhJetPtGenVsJetPtRecSubUeCone");
841 fhJetPtGenVsJetPtRecSubUeCone->SetTitle("fhJetPtGenVsJetPtRecSubUeCone");
842 fOutputList->Add(fhJetPtGenVsJetPtRecSubUeCone); // with weighted kT median bg subtr
843 //....
49dc2855 844 fhJetPtGen = new TH1D("fhJetPtGen","Jet Pt (MC Gen)",bw*100,0,200); //MC generator charged jet pt spectrum
8aa19603 845 fOutputList->Add(fhJetPtGen);
846 //....
49dc2855 847 fhJetPtSubUeMedianGen = new TH1D("fhJetPtSubUeMedianGen","Jet Pt - UE pT (MC Gen)",bw*110,-20,200);
8aa19603 848 fOutputList->Add(fhJetPtSubUeMedianGen); // with kT median bg subtr
849 //....
850 fhJetPtSubUeConeGen = (TH1D*) fhJetPtSubUeMedianGen->Clone("fhJetPtSubUeConeGen");
851 fOutputList->Add(fhJetPtSubUeConeGen); // with weighted kT median bg subtr
852
853 //....
854 if(fIsFullMC){
1afc79e9 855 fhJetPtGenChargVsJetPtGenFull = new TH2D("fhJetPtGenChargVsJetPtGenFull","fhJetPtGenChargVsJetPtGenFull", 100,0,200, 100,0,200);
8aa19603 856 fOutputList->Add(fhJetPtGenChargVsJetPtGenFull); //gen full MC jet pt versus gen charged jet MC pt
857 //....
1afc79e9 858 fhJetPtGenFull = new TH1D("fhJetPtGenFull","Jet Pt (MC Full jets Gen)",100,0,200); //MC generator full jet pt spectrum
8aa19603 859 fOutputList->Add(fhJetPtGenFull);
860 }
255cb71c 861 }
862 //....
80ac66f6 863 fh2NtriggersGen = (TH2F*) fh2Ntriggers->Clone("fh2NtriggersGen");
864 fh2NtriggersGen->SetTitle(Form("%s Gen MC",fh2Ntriggers->GetTitle()));
865 fOutputList->Add(fh2NtriggersGen);
7fc3d134 866
80ac66f6 867 //Centrality, A, pT, pTtrigg
868 fHJetSpecGen = (THnSparseF*) fHJetSpec->Clone("fHJetSpecGen");
869 fHJetSpecGen->SetTitle(Form("%s Gen MC",fHJetSpec->GetTitle()));
870 fOutputList->Add(fHJetSpecGen);
871
255cb71c 872 fHJetSpecSubUeMedianGen = (THnSparseF*) fHJetSpecSubUeMedian->Clone("fHJetSpecSubUeMedianGen");
873 fHJetSpecSubUeMedianGen->SetTitle(Form("%s Gen MC",fHJetSpecSubUeMedian->GetTitle()));
874 fOutputList->Add(fHJetSpecSubUeMedianGen);
875
8aa19603 876 fHJetSpecSubUeConeGen = (THnSparseF*) fHJetSpecSubUeCone->Clone("fHJetSpecSubUeConeGen");
877 fHJetSpecSubUeConeGen->SetTitle(Form("%s Gen MC",fHJetSpecSubUeCone->GetTitle()));
878 fOutputList->Add(fHJetSpecSubUeConeGen);
49dc2855 879 //---
880 fHJetPhiCorrGen = (THnSparseF*) fHJetPhiCorr->Clone("fHJetPhiCorrGen");
881 fHJetPhiCorrGen->SetTitle(Form("%s Gen MC",fHJetPhiCorr->GetTitle()));
882 fOutputList->Add(fHJetPhiCorrGen);
883
884 fHJetPhiCorrSubUeMedianGen = (THnSparseF*) fHJetPhiCorrSubUeMedian->Clone("fHJetPhiCorrSubUeMedianGen");
885 fHJetPhiCorrSubUeMedianGen->SetTitle(Form("%s Gen MC",fHJetPhiCorrSubUeMedian->GetTitle()));
886 fOutputList->Add(fHJetPhiCorrSubUeMedianGen);
887
888 fHJetPhiCorrSubUeConeGen = (THnSparseF*) fHJetPhiCorrSubUeCone->Clone("fHJetPhiCorrSubUeConeGen");
889 fHJetPhiCorrSubUeConeGen->SetTitle(Form("%s Gen MC", fHJetPhiCorrSubUeCone->GetTitle()));
890 fOutputList->Add(fHJetPhiCorrSubUeConeGen);
891
255cb71c 892 //---
893 fHJetUeMedianGen = (THnSparseF*) fHJetUeMedian->Clone("fHJetUeMedianGen");
894 fHJetUeMedianGen->SetTitle(Form("%s Gen MC", fHJetUeMedian->GetTitle()));
895 fOutputList->Add(fHJetUeMedianGen);
896
8aa19603 897 fHJetUeConeGen = (THnSparseF*) fHJetUeCone->Clone("fHJetUeConeGen");
898 fHJetUeConeGen->SetTitle(Form("%s Gen MC", fHJetUeCone->GetTitle()));
899 fOutputList->Add(fHJetUeConeGen);
255cb71c 900
8aa19603 901 if(fFillRespMx){
902 //track efficiency/contamination histograms eta versus pT
903 fhPtTrkTruePrimRec = new TH2D("fhPtTrkTruePrimRec","PtTrkTruePrimRec",100,0,20,18,-0.9,0.9);
904 fOutputList->Add(fhPtTrkTruePrimRec);
80ac66f6 905
8aa19603 906 fhPtTrkTruePrimGen = (TH2D*) fhPtTrkTruePrimRec->Clone("fhPtTrkTruePrimGen");
907 fhPtTrkTruePrimGen->SetTitle("PtTrkTruePrimGen");
908 fOutputList->Add(fhPtTrkTruePrimGen);
80ac66f6 909
8aa19603 910 fhPtTrkSecOrFakeRec = (TH2D*) fhPtTrkTruePrimRec->Clone("fhPtTrkSecOrFakeRec");
911 fhPtTrkSecOrFakeRec->SetTitle("PtTrkSecOrFakeRec");
912 fOutputList->Add(fhPtTrkSecOrFakeRec);
913 }
255cb71c 914
8aa19603 915 fHRhoUeMedianVsConeGen = (THnSparseF*) fHRhoUeMedianVsCone->Clone("hRhoUeMedianVsConeGen");
916 fHRhoUeMedianVsConeGen->SetTitle(Form("%s Gen MC", fHRhoUeMedianVsCone->GetTitle()));
917 fOutputList->Add(fHRhoUeMedianVsConeGen);
d1405a52
ML		 918
919 fhEntriesToMedianGen = (TH1D*) fhEntriesToMedian->Clone("fhEntriesToMedianGen");
920 fhEntriesToMedianGen->SetTitle(Form("%s Gen MC", fhEntriesToMedian->GetTitle()));
921 fOutputList->Add(fhEntriesToMedianGen);
922
923 fhCellAreaToMedianGen = (TH1D*) fhCellAreaToMedian->Clone("fhCellAreaToMedianGen");
924 fhCellAreaToMedianGen->SetTitle(Form("%s Gen MC", fhCellAreaToMedian->GetTitle()));
925 fOutputList->Add(fhCellAreaToMedianGen);
1e3ad34d 926
927 fhNofMultipleTriggersGen = (TH1D*) fhNofMultipleTriggers->Clone("fhNofMultipleTriggersGen");
928 fOutputList->Add(fhNofMultipleTriggersGen);
929
0fe881ea 930 fhNofMultipleTriggersConeGen = (TH1D*) fhNofMultipleTriggersCone->Clone("fhNofMultipleTriggersConeGen");
931 fOutputList->Add(fhNofMultipleTriggersConeGen);
932
30557bf5 933 fhNofMultipleTriggersConeGenLow = (TH1D*) fhNofMultipleTriggersCone->Clone("fhNofMultipleTriggersConeGenpta15to20");
934 fOutputList->Add(fhNofMultipleTriggersConeGenLow);
935
936 fhNofMultipleTriggersConeGenHigh = (TH1D*) fhNofMultipleTriggersCone->Clone("fhNofMultipleTriggersConeGenpta20to50");
937 fOutputList->Add(fhNofMultipleTriggersConeGenHigh);
938
939
940
1e3ad34d 941 fhDeltaRMultTriggersGen = (TH1D*) fhDeltaRMultTriggers->Clone("fhDeltaRMultTriggersGen");
942 fOutputList->Add(fhDeltaRMultTriggersGen);
943
30557bf5 944 fhDeltaPhiMultTriggersGen = (TH1D*) fhDeltaPhiMultTriggers->Clone("fhDeltaPhiMultTriggersGen");
945 fOutputList->Add(fhDeltaPhiMultTriggersGen);
946
947
6c1489e1 948 fhNofMultipleTriggersConeGenA = (TH1D*) fhNofMultipleTriggersConeGen->Clone("fhNofMultipleTriggersConeGen10");
949 fOutputList->Add(fhNofMultipleTriggersConeGenA);
950
30557bf5 951 fhNofMultipleTriggersConeGenALow = (TH1D*) fhNofMultipleTriggersConeGen->Clone("fhNofMultipleTriggersConeGen10pTa15to20");
952 fOutputList->Add(fhNofMultipleTriggersConeGenALow);
953
954 fhNofMultipleTriggersConeGenAHigh = (TH1D*) fhNofMultipleTriggersConeGen->Clone("fhNofMultipleTriggersConeGen10pTa20to50");
955 fOutputList->Add(fhNofMultipleTriggersConeGenAHigh);
956
6c1489e1 957 fhNofMultipleTriggersConeGenB = (TH1D*) fhNofMultipleTriggersConeGen->Clone("fhNofMultipleTriggersConeGen5");
958 fOutputList->Add(fhNofMultipleTriggersConeGenB);
959
30557bf5 960 fhNofMultipleTriggersConeGenBLow = (TH1D*) fhNofMultipleTriggersConeGen->Clone("fhNofMultipleTriggersConeGen5pTa15to20");
961 fOutputList->Add(fhNofMultipleTriggersConeGenBLow);
962
963 fhNofMultipleTriggersConeGenBHigh = (TH1D*) fhNofMultipleTriggersConeGen->Clone("fhNofMultipleTriggersConeGen5pTa20to50");
964 fOutputList->Add(fhNofMultipleTriggersConeGenBHigh);
965
966
967
80ac66f6 968 }
ea603b64 969 //-------------------------------------
970 // pythia histograms
971 const Int_t nBinPt = 150;
972 Double_t binLimitsPt[nBinPt+1];
973 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
974 if(iPt == 0){
975 binLimitsPt[iPt] = -50.0;
976 }else{// 1.0
977 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.;
978 }
979 }
980
981 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
982 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
983 fOutputList->Add(fh1Xsec);
984 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
985 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
986 fOutputList->Add(fh1Trials);
987 fh1AvgTrials = new TH1F("fh1AvgTrials","trials root file",1,0,1);
988 fh1AvgTrials->GetXaxis()->SetBinLabel(1,"#sum{avg ntrials}");
989 fOutputList->Add(fh1AvgTrials);
990 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
991 fOutputList->Add(fh1PtHard);
992 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
993 fOutputList->Add(fh1PtHardNoW);
994 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
995 fOutputList->Add(fh1PtHardTrials);
996
ad869500 997
998 // =========== Switch on Sumw2 for all histos ===========
999 for(Int_t i=0; i<fOutputList->GetEntries(); i++){
1000 TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
1001 if(h1){
1002 h1->Sumw2();
1003 continue;
1004 }
1005 THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
1006 if(hn){
1007 hn->Sumw2();
1008 }
1009 }
1010 TH1::AddDirectory(oldStatus);
1011
1012 PostData(1, fOutputList);
1013}
1014
1015//--------------------------------------------------------------------
1016
1017void AliAnalysisTaskJetCorePP::UserExec(Option_t *)
1018{
3d1fac8a 1019 //User Exec
1020
3d1fac8a 1021
ad869500 1022 //Event loop
ea603b64 1023 Double_t eventW = 1.0;
1024 Double_t ptHard = 0.0;
1025 Double_t nTrials = 1.0; // Trials for MC trigger
49dc2855 1026 if(fIsChargedMC || fIsKine) fh1AvgTrials->Fill("#sum{avg ntrials}",fAvgTrials);
ad869500 1027
1028 if(TMath::Abs((Float_t) fJetParamR)<0.00001){
8aa19603 1029 AliError("ANTIKT Cone radius is set to zero.");
1030 return;
1031 }
49dc2855 1032
8aa19603 1033 if(TMath::Abs((Float_t) fBgJetParamR)<0.00001){
49dc2855 1034 AliError("ANTIKT Cone radius is set to zero.");
255cb71c 1035 return;
1036 }
1037
49dc2855 1038 if(!fIsKine){ //real data or full MC
1039 if(!strlen(fJetBranchName.Data())){
1040 AliError("Name of jet branch not set.");
1041 return;
1042 }
1043 if(!strlen(fJetBranchNameBg.Data())){
1044 AliError("Name of jet bg branch not set.");
1045 return;
1046 }
1047 }else{ //Kine
3d1fac8a 1048 if(!strlen(fJetBranchNameBgKine.Data())){
1049 AliError("Name of jet bg branch for kine not set.");
1050 return;
1051 }
49dc2855 1052
1053 Init();
1054 if(fMcHandler){
1055 fMcEvent = fMcHandler->MCEvent();
1056 }else{
3d1fac8a 1057 if(fDebug > 1) printf("AliAnalysisTaskJetCorePP::Exec() fMcHandler=NULL\n");
49dc2855 1058 PostData(1, fOutputList);
1059 return;
1060 }
1061 if(!fMcEvent){
3d1fac8a 1062 if(fDebug > 1) printf("AliAnalysisTaskJetCorePP::Exec() fMcEvent=NULL \n");
49dc2855 1063 PostData(1, fOutputList);
1064 return;
1065 }
29ac2c9b 1066
1067 Float_t xsection = 0;
1068 Float_t trials = 0;
1069
1070 AliGenPythiaEventHeader *genPH =
1071 dynamic_cast<AliGenPythiaEventHeader*> (fMcEvent->GenEventHeader());
1072 if(genPH){
1073 xsection = genPH->GetXsection();
1074 trials = genPH->Trials();
1075 ptHard = genPH->GetPtHard();
1076 }
1077 fh1Xsec->Fill("<#sigma>",xsection);
1078 fh1Trials->Fill("#sum{ntrials}",trials);
1079 fh1PtHard->Fill(ptHard,eventW);
1080 fh1PtHardNoW->Fill(ptHard,1);
1081 fh1PtHardTrials->Fill(ptHard,trials);
49dc2855 1082 }
ad869500 1083
29ac2c9b 1084
ad869500 1085 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
1086 if(!fESD){
80ac66f6 1087 if(fDebug>1) AliError("ESD not available");
ad869500 1088 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
1089 }
1090
1091 fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
1092 AliAODEvent* aod = NULL;
1093 // take all other information from the aod we take the tracks from
49dc2855 1094 if(!aod && !fIsKine){
1095 if(!fESD) aod = fAODIn; //not ESD and not kine => input AOD
1096 else aod = fAODOut;// ESD or kine
ad869500 1097 }
1098
49dc2855 1099
ad869500 1100 if(fNonStdFile.Length()!=0){
1101 // case that we have an AOD extension we can fetch the jets from the extended output
1102 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1103 fAODExtension = aodH ? aodH->GetExtension(fNonStdFile.Data()) : 0;
1104 if(!fAODExtension){
1105 if(fDebug>1) Printf("AODExtension found for %s",fNonStdFile.Data());
1106 }
1107 }
1108
80ac66f6 1109 // ----------------- EVENT SELECTION --------------------------
ad869500 1110 fHistEvtSelection->Fill(1); // number of events before event selection
1111
49dc2855 1112 if(!fIsKine){
1113 // physics selection
1114 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
ad869500 1115 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
1116
49dc2855 1117 if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
1118 if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
1119 fHistEvtSelection->Fill(2);
1120 PostData(1, fOutputList);
1121 return;
1122 }
1123
ad869500 1124 //check AOD pointer
49dc2855 1125 if(!aod){
1126 if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
1127 fHistEvtSelection->Fill(3);
1128 PostData(1, fOutputList);
1129 return;
1130 }
ad869500 1131
49dc2855 1132 // vertex selection for reconstructed data
1133 AliAODVertex* primVtx = aod->GetPrimaryVertex();
1134
1135 if(!primVtx){
1136 if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
1137 fHistEvtSelection->Fill(3);
1138 PostData(1, fOutputList);
1139 return;
1140 }
ad869500 1141
49dc2855 1142 Int_t nTracksPrim = primVtx->GetNContributors();
1143 Float_t vtxz = primVtx->GetZ();
1144 //Input events
1145 fhContribVtx->Fill(nTracksPrim);
1146 if( nTracksPrim > 0 ) fhVertexZ->Fill(vtxz);
ad869500 1147
49dc2855 1148 if((nTracksPrim < fMinContribVtx) ||
1149 (vtxz < fVtxZMin) ||
1150 (vtxz > fVtxZMax)){
1151 if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",
ad869500 1152 (char*)__FILE__,__LINE__,vtxz);
49dc2855 1153 fHistEvtSelection->Fill(3);
1154 PostData(1, fOutputList);
1155 return;
1156 }else{
ad869500 1157 //Accepted events
49dc2855 1158 fhContribVtxAccept->Fill(nTracksPrim);
1159 fhVertexZAccept->Fill(vtxz);
1160 }
1161 }else{ //KINE cut on vertex
1162 const AliVVertex *vtxMC = fMcEvent->GetPrimaryVertex();
1163 Float_t zVtx = vtxMC->GetZ();
1164 fhVertexZ->Fill(zVtx);
1165 if(zVtx < fVtxZMin || zVtx > fVtxZMax){ //vertex cut
1166 fHistEvtSelection->Fill(3);
1167 PostData(1, fOutputList);
1168 return;
1169 }else{
1170 fhVertexZAccept->Fill(zVtx);
1171 }
ad869500 1172 }
1173 //FK// No event class selection imposed
1174 // event class selection (from jet helper task)
1175 //Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
1176 //if(fDebug) Printf("Event class %d", eventClass);
1177 //if(eventClass < fEvtClassMin || eventClass > fEvtClassMax){
1178 // fHistEvtSelection->Fill(4);
1179 // PostData(1, fOutputList);
1180 // return;
1181 //}
1182
80ac66f6 1183 //------------------ CENTRALITY SELECTION ---------------
ad869500 1184 AliCentrality *cent = 0x0;
1185 Double_t centValue = 0.0;
1186 if(fSystem){ //fSystem=0 for pp, fSystem=1 for pPb
1187 if(fESD){
1188 cent = fESD->GetCentrality();
1189 if(cent) centValue = cent->GetCentralityPercentile("V0M");
1190 }else{
0fe881ea 1191 //centValue = aod->GetHeader()->GetCentrality();
0a918d8d 1192 centValue = ((AliVAODHeader*)aod->GetHeader())->GetCentrality();
ad869500 1193 }
1194 if(fDebug) printf("centrality: %f\n", centValue);
1195 //Input events
1196 fhCentrality->Fill(centValue);
1197
1198 if(centValue < fCentMin || centValue > fCentMax){
1199 fHistEvtSelection->Fill(4);
1200 PostData(1, fOutputList);
1201 return;
1202 }else{
1203 //Accepted events
1204 fhCentralityAccept->Fill( centValue );
1205 }
1206 }
1207
8e103a56 1208 //-----------------select disjunct event subsamples ----------------
49dc2855 1209 if(!fIsKine){ //reconstructed data
0fe881ea 1210 //Int_t eventnum = aod->GetHeader()->GetEventNumberESDFile();
0a918d8d 1211 AliAODHeader * header = dynamic_cast<AliAODHeader*>(aod->GetHeader());
1212 if(!header) AliFatal("Not a standard AOD");
1213
1214 Int_t eventnum = header->GetEventNumberESDFile();
49dc2855 1215 Int_t lastdigit = eventnum % 10;
1216 if(!(fEventNumberRangeLow<=lastdigit && lastdigit<=fEventNumberRangeHigh)){
1217 fHistEvtSelection->Fill(5);
1218 PostData(1, fOutputList);
1219 return;
0fe881ea 1220 }
49dc2855 1221 }
8e103a56 1222
ad869500 1223 if(fDebug) std::cout<<" ACCEPTED EVENT "<<endl;
ad869500 1224 fHistEvtSelection->Fill(0); // accepted events
255cb71c 1225 // ==================== end event selection ============================
1226
1227 Double_t tmpArrayFive[5];
49dc2855 1228 Double_t tmpArrayFour[4];
ad869500 1229
ad869500 1230
8aa19603 1231 TList particleList; //list of tracks
49dc2855 1232 Int_t indexTrigg = -1;
1233 Double_t rhoFromCellMedian=0.0, rhoCone=0.0;
80ac66f6 1234
49dc2855 1235 if(!fIsKine){
1236 //=============== Reconstructed antikt jets ===============
1237 ReadTClonesArray(fJetBranchName.Data() , fListJets);
1238 ReadTClonesArray(fJetBranchNameBg.Data(), fListJetsBg);
80ac66f6 1239
49dc2855 1240 //============ Estimate background in reconstructed events ===========
1241
1242 //Find Hadron trigger and Estimate rho from cone
1243 indexTrigg = GetListOfTracks(&particleList); //index of trigger hadron in Particle list
1244 EstimateBgCone(fListJets, &particleList, rhoCone);
1245
1246 //Estimate rho from cell median minus jets
1247 EstimateBgRhoMedian(fListJetsBg, &particleList, rhoFromCellMedian,0); //real data
1248 }
255cb71c 1249 //============== analyze generator level MC ================
1250 TList particleListGen; //list of tracks in MC
80ac66f6 1251
49dc2855 1252 if(fIsChargedMC || fIsKine){
1253
1254 if(fIsKine){ //pure kine
160dbdb8 1255
49dc2855 1256 //================= generated charged antikt jets ================
1257 ReadTClonesArray(fJetBranchNameKine.Data(), fListJetsGen);
3d1fac8a 1258 ReadTClonesArray(fJetBranchNameBgKine.Data(), fListJetsBgGen);
49dc2855 1259 }else{
1260 //================= generated charged antikt jets ================
1261 ReadTClonesArray(fJetBranchNameChargMC.Data(), fListJetsGen);
1262 ReadTClonesArray(fJetBranchNameBgChargMC.Data(), fListJetsBgGen);
1263
1264 if(fIsFullMC){ //generated full jets
1265 ReadTClonesArray(fJetBranchNameFullMC.Data(), fListJetsGenFull);
1266 }
1267 }
255cb71c 1268 //========================================================
80ac66f6 1269 //serarch for charged trigger at the MC generator level
255cb71c 1270
8e103a56 1271 Int_t indexTriggGen = -1;
1272 Double_t ptTriggGen = -1;
255cb71c 1273 Int_t iCounterGen = 0; //number of entries in particleListGen array
8e103a56 1274 Int_t triggersMC[200];//list of trigger candidates
1275 Int_t ntriggersMC = 0; //index in triggers array
6c1489e1 1276 Int_t triggersMCa[200]; //list of trigger candidates 10%eloss
1277 Int_t ntriggersMCa = 0; //index in triggers array 10%eloss
1278 Int_t triggersMCb[200]; //list of trigger candidates 5%eloss
1279 Int_t ntriggersMCb = 0; //index in triggers array 5%eloss
1280
8e103a56 1281
49dc2855 1282 if(!fIsKine){
1283 if(fESD){//ESD input
80ac66f6 1284
49dc2855 1285 AliMCEvent* mcEvent = MCEvent();
1286 if(!mcEvent){
1287 PostData(1, fOutputList);
1288 return;
1289 }
1290
1291 AliGenPythiaEventHeader *pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
1292 if(pythiaGenHeader){
1293 nTrials = pythiaGenHeader->Trials();
1294 ptHard = pythiaGenHeader->GetPtHard();
1295 fh1PtHard->Fill(ptHard,eventW);
1296 fh1PtHardNoW->Fill(ptHard,1);
1297 fh1PtHardTrials->Fill(ptHard,nTrials);
1298 }
ea603b64 1299
49dc2855 1300 for(Int_t it = 0; it < mcEvent->GetNumberOfTracks(); it++){
1301 if(!mcEvent->IsPhysicalPrimary(it)) continue;
1302 AliMCParticle* part = (AliMCParticle*) mcEvent->GetTrack(it);
1303 if(!part) continue;
1304 if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){
1305
1306 if(fHardest==0 && ntriggersMC<200){//single inclusive trigger
1307 if(indexTriggGen > -1){//trigger candidate was found
1308 triggersMC[ntriggersMC] = indexTriggGen;
1309 ntriggersMC++;
1310 }
1311 }
1312
1313 iCounterGen++;//index in particleListGen array
1314 }
1315 }
1316 }else{ //AOD input
1317 if(!fAODIn){
1318 PostData(1, fOutputList);
1319 return;
1320 }
1321 TClonesArray *tca = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(AliAODMCParticle::StdBranchName()));
1322 if(!tca){
1323 PostData(1, fOutputList);
1324 return;
1325 }
1326 for(Int_t it = 0; it < tca->GetEntriesFast(); it++){
1327 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1328 if(!part) continue;
1329 if(!part->IsPhysicalPrimary()) continue;
1330 if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){
1331
1332 if(fHardest==0 && ntriggersMC<200){//single inclusive trigger
1333 if(indexTriggGen > -1){ //trigger candidater was found
1334 triggersMC[ntriggersMC] = indexTriggGen;
1335 ntriggersMC++;
1336 }
1337 }
1338
1339 iCounterGen++;//number of entries in particleListGen array
1340 }
1341 }
ea603b64 1342 }
49dc2855 1343 }else{ //analyze kine
ea603b64 1344
49dc2855 1345 for(Int_t it = 0; it < fMcEvent->GetNumberOfTracks(); it++){
1346 if(!fMcEvent->IsPhysicalPrimary(it)) continue;
1347 AliMCParticle* part = (AliMCParticle*) fMcEvent->GetTrack(it);
80ac66f6 1348 if(!part) continue;
1349 if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){
49dc2855 1350
8e103a56 1351 if(fHardest==0 && ntriggersMC<200){//single inclusive trigger
255cb71c 1352 if(indexTriggGen > -1){//trigger candidate was found
8e103a56 1353 triggersMC[ntriggersMC] = indexTriggGen;
1354 ntriggersMC++;
1355 }
1356 }
49dc2855 1357
255cb71c 1358 iCounterGen++;//index in particleListGen array
80ac66f6 1359 }
1360 }
80ac66f6 1361 }
6c1489e1 1362
1363 if(fHardest==0){
1364 Int_t npar = particleListGen.GetEntries();
1365 for(Int_t ip=0; ip < npar; ip++){
1366 AliVParticle *part = (AliVParticle*) particleListGen.At(ip);
1367 if(!part) continue;
1368
1369 Double_t pta = 0.9 * part->Pt(); //10% energy loss
1370 Double_t ptb = 0.95 * part->Pt(); //5% energy loss
1371 if(fTriggerPtRangeLow <= pta && pta < fTriggerPtRangeHigh && ntriggersMCa<200){
1372 triggersMCa[ntriggersMCa] = ip;
1373 ntriggersMCa++;
1374 }
1375
1376 if(fTriggerPtRangeLow <= ptb && ptb < fTriggerPtRangeHigh && ntriggersMCb<200){
1377 triggersMCb[ntriggersMCb] = ip;
1378 ntriggersMCb++;
1379 }
1380 }
1381
1382 if(ntriggersMCa>0){
1383 Int_t rnda = fRandom->Integer(ntriggersMCa); //0 to ntriggers-1
1384 Int_t indexTriggGena = triggersMCa[rnda];
1385
1386 Double_t deltaPhia, deltaEtaa, deltaRa;
30557bf5 1387 Int_t aaLow = 0;
1388 Int_t aaHigh = 0;
6c1489e1 1389
1390 //Correlation with single inclusive TRIGGER
1391 AliVParticle* tGenTa = (AliVParticle*) particleListGen.At(indexTriggGena);
1392 if(tGenTa){
58003b11 1393 //for(Int_t ia=0; ia<ntriggersMCa; ia++)
1394 for(Int_t ia=0; ia< npar; ia++){
1395 //if(indexTriggGena == triggersMCa[ia]) continue;
1396 if(indexTriggGena == ia) continue;
6c1489e1 1397
58003b11 1398 //AliVParticle* tGenTz = (AliVParticle*) particleListGen.At(triggersMCa[ia]);
1399 AliVParticle* tGenTz = (AliVParticle*) particleListGen.At(ia);
6c1489e1 1400 if(!tGenTz) continue;
58003b11 1401 if(tGenTz->Pt()*0.9<15.0) continue;
30557bf5 1402 if(tGenTz->Pt()*0.9>50.0) continue;
6c1489e1 1403
1404 deltaPhia = RelativePhi(tGenTa->Phi(),tGenTz->Phi());
1405 deltaEtaa = tGenTa->Eta()-tGenTz->Eta();
1406 deltaRa = sqrt(deltaPhia*deltaPhia + deltaEtaa*deltaEtaa);
1407
30557bf5 1408 if(deltaRa<0.4){
1409 if(tGenTz->Pt()*0.9<20.0) aaLow++; //15-20
1410 else aaHigh++; //20-50
1411 }
6c1489e1 1412 }
1413 }
30557bf5 1414 fhNofMultipleTriggersConeGenA->Fill(aaLow+aaHigh); // 15-50
1415 fhNofMultipleTriggersConeGenALow->Fill(aaLow); //15-20
1416 fhNofMultipleTriggersConeGenAHigh->Fill(aaHigh);//20-50
6c1489e1 1417 }
1418
1419 if(ntriggersMCb>0){
1420 Int_t rndb = fRandom->Integer(ntriggersMCb); //0 to ntriggers-1
1421 Int_t indexTriggGenb = triggersMCb[rndb];
1422
1423 Double_t deltaPhib, deltaEtab, deltaRb;
30557bf5 1424 Int_t bbLow = 0;
1425 Int_t bbHigh = 0;
6c1489e1 1426
1427 //Correlation with single inclusive TRIGGER
1428 AliVParticle* tGenTb = (AliVParticle*) particleListGen.At(indexTriggGenb);
1429 if(tGenTb){
58003b11 1430 //for(Int_t ib=0; ib<ntriggersMCb; ib++)
1431 for(Int_t ib=0; ib<npar; ib++){
1432 //if(indexTriggGenb == triggersMCb[ib]) continue;
1433 if(indexTriggGenb == ib) continue;
6c1489e1 1434
58003b11 1435 //AliVParticle* tGenTz = (AliVParticle*) particleListGen.At(triggersMCb[ib]);
1436 AliVParticle* tGenTz = (AliVParticle*) particleListGen.At(ib);
6c1489e1 1437 if(!tGenTz) continue;
58003b11 1438 if(tGenTz->Pt()*0.95<15.0) continue;
30557bf5 1439 if(tGenTz->Pt()*0.95>50.0) continue;
6c1489e1 1440
1441 deltaPhib = RelativePhi(tGenTb->Phi(),tGenTz->Phi());
1442 deltaEtab = tGenTb->Eta()-tGenTz->Eta();
1443 deltaRb = sqrt(deltaPhib*deltaPhib + deltaEtab*deltaEtab);
1444
30557bf5 1445 if(deltaRb<0.4){
1446 if(tGenTz->Pt()*0.95<20.0) bbLow++; //15-20
1447 else bbHigh++;
1448 }
6c1489e1 1449 }
1450 }
30557bf5 1451 fhNofMultipleTriggersConeGenB->Fill(bbLow+bbHigh); //15-50
1452 fhNofMultipleTriggersConeGenBLow->Fill(bbLow);//15-20
1453 fhNofMultipleTriggersConeGenBHigh->Fill(bbHigh);//20-50
6c1489e1 1454 }
1455 }
1456
1457
1458
8aa19603 1459 //============== Estimate bg in generated events ==============
1460 Double_t rhoFromCellMedianGen=0.0, rhoConeGen=0.0;
1461
1462 //Estimate rho from cone
1463 EstimateBgCone(fListJetsGen, &particleListGen, rhoConeGen);
1464
1465 //Estimate rho from cell median minus jets
d1405a52 1466 EstimateBgRhoMedian(fListJetsBgGen, &particleListGen, rhoFromCellMedianGen,1);//mc data
8aa19603 1467
255cb71c 1468 //============ Generator trigger+jet ==================
8e103a56 1469 if(fHardest==0){ //single inclusive trigger
1470 if(ntriggersMC>0){ //there is at least one trigger
1471 Int_t rnd = fRandom->Integer(ntriggersMC); //0 to ntriggers-1
1472 indexTriggGen = triggersMC[rnd];
1e3ad34d 1473
1474 fhNofMultipleTriggersGen->Fill(ntriggersMC-1);
0fe881ea 1475
ac4f24f7 1476 Double_t deltaPhi, deltaEta, deltaR;
30557bf5 1477 Int_t iLow = 0;
1478 Int_t iHigh = 0;
ac4f24f7 1479
1480 //Correlation with single inclusive TRIGGER
1481 AliVParticle* tGenT1 = (AliVParticle*) particleListGen.At(indexTriggGen);
1482 if(tGenT1){
58003b11 1483 //for(Int_t ia=0; ia<ntriggersMC; ia++)
1484 for(Int_t ia=0; ia< particleListGen.GetEntries(); ia++){
1485 //if(indexTriggGen == triggersMC[ia]) continue;
1486 if(indexTriggGen == ia) continue;
ac4f24f7 1487
58003b11 1488 //AliVParticle* tGenT2 = (AliVParticle*) particleListGen.At(triggersMC[ia]);
1489 AliVParticle* tGenT2 = (AliVParticle*) particleListGen.At(ia);
ac4f24f7 1490 if(!tGenT2) continue;
58003b11 1491 if(tGenT2->Pt()<15.0) continue;
30557bf5 1492 if(tGenT2->Pt()>50.0) continue;
ac4f24f7 1493 deltaPhi = RelativePhi(tGenT1->Phi(),tGenT2->Phi());
1494 deltaEta = tGenT1->Eta()-tGenT2->Eta();
1495 deltaR = sqrt(deltaPhi*deltaPhi + deltaEta*deltaEta);
0fe881ea 1496
30557bf5 1497 fhDeltaRMultTriggersGen->Fill(deltaR);
1498 fhDeltaPhiMultTriggersGen->Fill(deltaPhi);
1499
1500 if(deltaR<0.4){
1501 if(tGenT2->Pt()<20.0) iLow++;
1502 else iHigh++;
1e3ad34d 1503 }
1504 }
30557bf5 1505 }
1506 fhNofMultipleTriggersConeGen->Fill(iLow+iHigh);
1507 fhNofMultipleTriggersConeGenLow->Fill(iLow);
1508 fhNofMultipleTriggersConeGenHigh->Fill(iHigh);
1e3ad34d 1509
80ac66f6 1510 }else{
8e103a56 1511 indexTriggGen = -1; //trigger not found
80ac66f6 1512 }
1513 }
8e103a56 1514
255cb71c 1515 //-----------
8e103a56 1516 Int_t ilowGen = (fHardest==0 || fHardest==1) ? indexTriggGen : 0;
1517 Int_t ihighGen = (fHardest==0 || fHardest==1) ? indexTriggGen+1 : particleListGen.GetEntries();
255cb71c 1518 Bool_t fillOnceGen = kTRUE;
1519 //-----------
1520
8e103a56 1521 for(Int_t igen= ilowGen; igen< ihighGen; igen++){ //loop over possible trigger
1522 indexTriggGen = igen; //trigger hadron
1523
1524 if(indexTriggGen == -1) continue;
1525 AliVParticle* triggerGen=(AliVParticle*) particleListGen.At(indexTriggGen);
1526 if(!triggerGen) continue;
1527
1528 if(fHardest >= 2){
6168afc3 1529 if(triggerGen->Pt() < 6.0) continue; //all hadrons pt>6
8e103a56 1530 }
1531 if(TMath::Abs((Float_t) triggerGen->Eta()) > fTriggerEtaCut) continue;
1532
1533 ptTriggGen = triggerGen->Pt(); //count triggers
1534 fh2NtriggersGen->Fill(centValue, ptTriggGen);
1535
80ac66f6 1536 //Count jets and trigger-jet pairs at MC generator level
1537 for(Int_t ij=0; ij<fListJetsGen->GetEntries(); ij++){
1538 AliAODJet* jet = (AliAODJet*)(fListJetsGen->At(ij));
1539 if(!jet) continue;
1540 Double_t etaJetGen = jet->Eta();
255cb71c 1541 Double_t areaJetGen = jet->EffectiveAreaCharged();
80ac66f6 1542
1543 if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){
80ac66f6 1544
49dc2855 1545 Double_t ptUeFromCellMedianGen = rhoFromCellMedianGen*areaJetGen;
1546 Double_t ptUeConeGen = rhoConeGen*areaJetGen;
1547
1548 //Rongrong's analysis
1549 Double_t dPhiGen = jet->Phi() - triggerGen->Phi();
1550 if(dPhiGen>2*TMath::Pi()) dPhiGen -= 2*TMath::Pi();
1551 if(dPhiGen<-2*TMath::Pi()) dPhiGen += 2*TMath::Pi();
1552 if(dPhiGen<-0.5*TMath::Pi()) dPhiGen += 2*TMath::Pi();
1553 if(dPhiGen>1.5*TMath::Pi()) dPhiGen -= 2*TMath::Pi();
1554
1555 //[A,pTjet,pTtrig,dphi]
1556 tmpArrayFour[0] = areaJetGen;
1557 tmpArrayFour[1] = jet->Pt();
1558 tmpArrayFour[2] = ptTriggGen;
1559 tmpArrayFour[3] = dPhiGen;
1560
1561 fHJetPhiCorrGen->Fill(tmpArrayFour); // Dphi distribution jet-triger
1562
1563 //[A,pTjet-pTUe,pTtrig,dphi]
1564 tmpArrayFour[1] = jet->Pt() - ptUeFromCellMedianGen;
1565 fHJetPhiCorrSubUeMedianGen->Fill(tmpArrayFour);
1566
1567 //[A,pTjet-pTUe,pTtrig,dphi]
1568 tmpArrayFour[1] = jet->Pt() - ptUeConeGen;
1569 fHJetPhiCorrSubUeConeGen->Fill(tmpArrayFour);
1570
1571 //Leticia's analysis
8e103a56 1572 Double_t dphi = RelativePhi(triggerGen->Phi(), jet->Phi());
1573 if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue;
80ac66f6 1574
8e103a56 1575 //Centrality, A, pT, pTtrigg
255cb71c 1576 tmpArrayFive[0] = centValue;
1577 tmpArrayFive[1] = areaJetGen;
1578 tmpArrayFive[2] = jet->Pt();
1579 tmpArrayFive[3] = ptTriggGen;
1580 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1581 fHJetSpecGen->Fill(tmpArrayFive);
1582
8aa19603 1583 //Centrality, A, pTjet-pTbgCellMedian, pTtrigg, dphi
255cb71c 1584 tmpArrayFive[0] = centValue;
1585 tmpArrayFive[1] = areaJetGen;
8aa19603 1586 tmpArrayFive[2] = jet->Pt() - ptUeFromCellMedianGen;
255cb71c 1587 tmpArrayFive[3] = ptTriggGen;
1588 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1589 fHJetSpecSubUeMedianGen->Fill(tmpArrayFive);
1590
8aa19603 1591 //Centrality, A, pTjet-pTbgCone, pTtrigg, dphi
255cb71c 1592 tmpArrayFive[0] = centValue;
1593 tmpArrayFive[1] = areaJetGen;
8aa19603 1594 tmpArrayFive[2] = jet->Pt() - ptUeConeGen;
255cb71c 1595 tmpArrayFive[3] = ptTriggGen;
1596 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
8aa19603 1597 fHJetSpecSubUeConeGen->Fill(tmpArrayFive);
255cb71c 1598
1599 //Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", kT median
1600 tmpArrayFive[0] = areaJetGen;
1601 tmpArrayFive[1] = jet->Pt();
8aa19603 1602 tmpArrayFive[2] = jet->Pt() - ptUeFromCellMedianGen;
1603 tmpArrayFive[3] = ptUeFromCellMedianGen;
1604 tmpArrayFive[4] = rhoFromCellMedianGen;
255cb71c 1605 fHJetUeMedianGen->Fill(tmpArrayFive);
1606
8aa19603 1607 //Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", perpendicular Cone
255cb71c 1608 tmpArrayFive[0] = areaJetGen;
1609 tmpArrayFive[1] = jet->Pt();
8aa19603 1610 tmpArrayFive[2] = jet->Pt() - ptUeConeGen;
1611 tmpArrayFive[3] = ptUeConeGen;
1612 tmpArrayFive[4] = rhoConeGen;
1613 fHJetUeConeGen->Fill(tmpArrayFive);
255cb71c 1614
1615 if(fillOnceGen){
8aa19603 1616 Double_t fillRhoGen[] = { rhoConeGen, rhoFromCellMedianGen};
1617 fHRhoUeMedianVsConeGen->Fill(fillRhoGen);
255cb71c 1618 fillOnceGen = kFALSE;
1619 }
8e103a56 1620 }//back to back jet-trigger pair
1621 }//jet loop
1622 }//trigger loop
1623
49dc2855 1624 if(fFillRespMx && !fIsKine){
8aa19603 1625 //================ RESPONSE MATRIX ===============
1626 //Count jets and trigger-jet pairs at MC generator level
1627 for(Int_t ij=0; ij<fListJetsGen->GetEntries(); ij++){
1628 AliAODJet* jet = (AliAODJet*)(fListJetsGen->At(ij));
1629 if(!jet) continue;
1630 Double_t etaJetGen = jet->Eta();
1631 Double_t ptJetGen = jet->Pt();
255cb71c 1632
8aa19603 1633 if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){
1634 fhJetPtGen->Fill(ptJetGen); // gen level pt spectum of jets response mx normalization
80ac66f6 1635
8aa19603 1636 Double_t areaJetGen = jet->EffectiveAreaCharged();
1637 Double_t ptUeFromCellMedianGen = rhoFromCellMedianGen*areaJetGen;
1638 Double_t ptUeConeGen = rhoConeGen*areaJetGen;
80ac66f6 1639
8aa19603 1640 fhJetPtSubUeMedianGen->Fill(ptJetGen - ptUeFromCellMedianGen);
1641 fhJetPtSubUeConeGen->Fill(ptJetGen - ptUeConeGen);
160dbdb8 1642 }
1643 }
8aa19603 1644 if(fListJets->GetEntries()>0 && fListJetsGen->GetEntries()>0){ //at least some reconstructed jets
1645
1646 Int_t ng = (Int_t) fListJetsGen->GetEntries();
1647 Int_t nr = (Int_t) fListJets->GetEntries();
160dbdb8 1648
8aa19603 1649 //Find closest MC generator - reconstructed jets
1650 if(faGenIndex.GetSize()<nr) faGenIndex.Set(nr); //idx of gen jet assoc to rec jet
1651 if(faRecIndex.GetSize()<ng) faRecIndex.Set(ng); //idx of rec jet assoc to gen jet
160dbdb8 1652
1653 if(fDebug){
8aa19603 1654 Printf("New Rec List %d gen index Array %d",nr,faGenIndex.GetSize());
1655 Printf("New Gen List %d rec index Array %d",ng,faRecIndex.GetSize());
160dbdb8 1656 }
1657 //matching of MC genrator level and reconstructed jets
8aa19603 1658 AliAnalysisHelperJetTasks::GetClosestJets(fListJetsGen,ng,fListJets,nr,faGenIndex,faRecIndex,fDebug);
1659
1660 // Fill response matrix
1661 for(Int_t ir = 0; ir < nr; ir++){
1662 AliAODJet *recJet = (AliAODJet*) fListJets->At(ir);
1663 Double_t etaJetRec = recJet->Eta();
1664 Double_t ptJetRec = recJet->Pt();
1665 Double_t areaJetRec = recJet->EffectiveAreaCharged();
160dbdb8 1666 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1667
8aa19603 1668 if((fJetEtaMin <= etaJetRec) && (etaJetRec <= fJetEtaMax)){
1669 Int_t ig = faGenIndex[ir]; //associated generator level jet
1670 if(ig >= 0 && ig < ng){
1671 if(fDebug > 10) Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1672 AliAODJet *genJet = (AliAODJet*) fListJetsGen->At(ig);
1673 Double_t ptJetGen = genJet->Pt();
1674 Double_t etaJetGen = genJet->Eta();
160dbdb8 1675
255cb71c 1676 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
8aa19603 1677 if((fJetEtaMin <= etaJetGen) && (etaJetGen <= fJetEtaMax)){
1678 fhJetPtGenVsJetPtRec->Fill(ptJetRec, ptJetGen);
1679
1680 Double_t areaJetGen = genJet->EffectiveAreaCharged();
1681 Double_t ptUeFromCellMedianGen = rhoFromCellMedianGen*areaJetGen;
1682 Double_t ptUeConeGen = rhoConeGen*areaJetGen;
1683 Double_t ptUeFromCellMedianRec = rhoFromCellMedian*areaJetRec;
1684 Double_t ptUeConeRec = rhoCone*areaJetRec;
1685 fhJetPtGenVsJetPtRecSubUeMedian->Fill(ptJetRec-ptUeFromCellMedianRec,
1686 ptJetGen-ptUeFromCellMedianGen);
1687 fhJetPtGenVsJetPtRecSubUeCone->Fill(ptJetRec-ptUeConeRec, ptJetGen-ptUeConeGen);
160dbdb8 1688 }
8aa19603 1689 }//ig>=0
160dbdb8 1690 }//rec jet in eta acceptance
1691 }//loop over reconstructed jets
1692 }// # of rec jets >0
8aa19603 1693
1694 //=========================== Full jet vs charged jet matrix ==========================
1695 if(fIsFullMC){
1696 //Count full jets and charged-jet pairs at MC generator level
1697 for(Int_t ij=0; ij<fListJetsGenFull->GetEntries(); ij++){
1698 AliAODJet* jetFull = (AliAODJet*)(fListJetsGenFull->At(ij));
1699 if(!jetFull) continue;
1700 Double_t etaJetFull = jetFull->Eta();
1701 Double_t ptJetFull = jetFull->Pt();
1702
1703 if((fJetEtaMin<=etaJetFull) && (etaJetFull<=fJetEtaMax)){
1704 fhJetPtGenFull->Fill(ptJetFull); // generator level pt spectum of full jets
1705 }
1706 }
1707 if(fListJetsGen->GetEntries()>0 && fListJetsGenFull->GetEntries()>0){ //at least some reconstructed jets
1708 Int_t nful = (Int_t) fListJetsGenFull->GetEntries();
1709 Int_t nchr = (Int_t) fListJetsGen->GetEntries();
1710
1711 //Find closest MC generator full - charged jet
1712 if(faGenIndex.GetSize()<nchr) faGenIndex.Set(nchr); //idx of gen FULL jet assoc to gen CHARGED jet
1713 if(faRecIndex.GetSize()<nful) faRecIndex.Set(nful); //idx of gen CHARGED jet assoc to gen FULL jet
1714
1715 if(fDebug){
1716 Printf("New Charg List %d Full index Array %d",nchr,faGenIndex.GetSize());
1717 Printf("New Full List %d Charg index Array %d",nful,faRecIndex.GetSize());
1718 }
1719 //matching of MC genrator level and reconstructed jets
1720 AliAnalysisHelperJetTasks::GetClosestJets(fListJetsGenFull,nful,fListJetsGen,nchr,faGenIndex,faRecIndex,fDebug);
1721
1722 // Fill response matrix
1723 for(Int_t ichr = 0; ichr < nchr; ichr++){ //charged jet loop
1724 AliAODJet *chJet = (AliAODJet*) fListJetsGen->At(ichr);
1725 Double_t etaJetCh = chJet->Eta();
1726 Double_t ptJetCh = chJet->Pt();
1727 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1728
1729 if((fJetEtaMin <= etaJetCh) && (etaJetCh <= fJetEtaMax)){
1730 Int_t iful = faGenIndex[ichr]; //associated generator level jet
1731 if(iful >= 0 && iful < nful){
1732 if(fDebug > 10) Printf("%s:%d iful = %d ichr = %d",(char*)__FILE__,__LINE__,iful,ichr);
1733 AliAODJet *genJetFull = (AliAODJet*) fListJetsGenFull->At(iful);
1734 Double_t ptJetFull = genJetFull->Pt();
1735 Double_t etaJetFull = genJetFull->Eta();
1736
1737 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1738 if((fJetEtaMin <= etaJetFull) && (etaJetFull <= fJetEtaMax)){
1739 fhJetPtGenChargVsJetPtGenFull->Fill(ptJetFull,ptJetCh);
1740 }
1741 }//iful>=0
1742 }//rec jet in eta acceptance
1743 }//loop over reconstructed jets
1744 }// # of rec jets >0
1745 }//pointer MC generator jets
1746 } //fill resp mx only for bin
1747 }//analyze generator level MC
1748
1749
49dc2855 1750 if(fIsKine){ //skip reconstructed data analysis in case of kine
1751 PostData(1, fOutputList);
1752 return;
1753 }
8e103a56 1754 //============= RECONSTRUCTED INCLUSIVE JETS ===============
ad869500 1755
ad869500 1756 Double_t etaJet = 0.0;
1757 Double_t pTJet = 0.0;
1758 Double_t areaJet = 0.0;
1759 Double_t phiJet = 0.0;
9f10e189 1760 //Int_t indexLeadingJet = -1;
1761 //Double_t pTLeadingJet = -10.0;
1762 //Double_t areaLeadingJet = -10.0;
8e103a56 1763
ad869500 1764 for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
1765 AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
1766 if(!jet) continue;
1767 etaJet = jet->Eta();
1768 phiJet = jet->Phi();
1769 pTJet = jet->Pt();
1770 if(pTJet==0) continue;
1771
1772 if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
255cb71c 1773 /*areaJet = jet->EffectiveAreaCharged();*/
7fc3d134 1774
ad869500 1775 //Jet Diagnostics---------------------------------
8e103a56 1776 fhJetPhi->Fill(pTJet, RelativePhi(phiJet,0.0)); //phi -pi,pi
1777 fhJetEta->Fill(pTJet, etaJet);
7fc3d134 1778 //search for leading jet
255cb71c 1779 /*if(pTJet > pTLeadingJet){
7fc3d134 1780 indexLeadingJet = ij;
1781 pTLeadingJet = pTJet;
7fc3d134 1782 areaLeadingJet = areaJet;
1783 }
1784
1785 // raw spectra of INCLUSIVE jets
80ac66f6 1786 //Centrality, pTjet, A
7fc3d134 1787 Double_t fillraw[] = { centValue,
1788 pTJet,
7fc3d134 1789 areaJet
1790 };
255cb71c 1791 fHJetPtRaw->Fill(fillraw);*/
8e103a56 1792 }
255cb71c 1793 /*
8e103a56 1794 if(indexLeadingJet > -1){
1795 // raw spectra of LEADING jets
1796 //Centrality, pTjet, A
1797 Double_t fillleading[] = { centValue,
1798 pTLeadingJet,
1799 areaLeadingJet
1800 };
1801 fHLeadingJetPtRaw->Fill(fillleading);
1802 }
255cb71c 1803 */
8e103a56 1804
1805 // =============== RECONSTRUCTED TRIGGER-JET PAIRS ================
8aa19603 1806 if(fIsChargedMC && fFillRespMx){
1807 FillEffHistos(&particleList, &particleListGen); //Fill efficiency histos
1808 }
255cb71c 1809 Bool_t filledOnce = kTRUE; //fill rho histogram only once per event
8e103a56 1810 //set ranges of the trigger loop
1811 Int_t ilow = (fHardest==0 || fHardest==1) ? indexTrigg : 0;
1812 Int_t ihigh = (fHardest==0 || fHardest==1) ? indexTrigg+1 : particleList.GetEntries();
7fc3d134 1813
8e103a56 1814 for(Int_t itrk= ilow; itrk< ihigh; itrk++){ //loop over possible trigger
1815 indexTrigg = itrk; //trigger hadron with pT >10 GeV
7fc3d134 1816
8e103a56 1817 if(indexTrigg < 0) continue;
1818
1819 AliVParticle *triggerHadron = (AliVParticle*) particleList.At(indexTrigg);
1820 if(!triggerHadron){
1821 PostData(1, fOutputList);
1822 return;
1823 }
1824 if(fHardest >= 2){
7d9fd386 1825 if(triggerHadron->Pt() < 6.0) continue; //all hadrons pt>6
8e103a56 1826 }
1827 if(TMath::Abs((Float_t) triggerHadron->Eta())> fTriggerEtaCut) continue;
7fc3d134 1828
8e103a56 1829 //Fill trigger histograms
1830 fh2Ntriggers->Fill(centValue,triggerHadron->Pt()); //trigger pT
1831 fhTriggerPhi->Fill(triggerHadron->Pt(),RelativePhi(triggerHadron->Phi(),0.0)); //phi -pi,pi
1832 fhTriggerEta->Fill(triggerHadron->Pt(),triggerHadron->Eta());
7fc3d134 1833
1834
8e103a56 1835 //---------- make trigger-jet pairs ---------
255cb71c 1836 //Int_t injet4 = 0;
1837 //Int_t injet = 0;
255cb71c 1838
8e103a56 1839 for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
1840 AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
1841 if(!jet) continue;
1842 etaJet = jet->Eta();
1843 phiJet = jet->Phi();
1844 pTJet = jet->Pt();
1845 if(pTJet==0) continue;
1846
1847 if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
1848 areaJet = jet->EffectiveAreaCharged();
49dc2855 1849
1850 //subtract bg using estinates base on median of kT jets
1851 Double_t ptUeFromCellMedian = rhoFromCellMedian*areaJet;
1852 Double_t ptUeCone = rhoCone*areaJet;
1853
255cb71c 1854 //if(areaJet >= 0.07) injet++;
1855 //if(areaJet >= 0.4) injet4++;
8e103a56 1856
1857 Double_t dphi = RelativePhi(triggerHadron->Phi(), phiJet);
1858 fhDphiTriggerJet->Fill(dphi); //Input
1859
1860 //Dphi versus jet pT
1861 //Centrality, Dphi=phiTrig-phiJet, pTjet, pTtrigg
255cb71c 1862 /*Double_t filldp[] = { centValue,
8e103a56 1863 dphi,
1864 pTJet,
1865 triggerHadron->Pt()
1866 };
1867 fHDphiVsJetPtAll->Fill(filldp);
255cb71c 1868 */
49dc2855 1869 //Rongrong's analysis
1870
1871 Double_t dPhi = phiJet - triggerHadron->Phi();
1872 if(dPhi>2*TMath::Pi()) dPhi -= 2*TMath::Pi();
1873 if(dPhi<-2*TMath::Pi()) dPhi += 2*TMath::Pi();
1874 if(dPhi<-0.5*TMath::Pi()) dPhi += 2*TMath::Pi();
1875 if(dPhi>1.5*TMath::Pi()) dPhi -= 2*TMath::Pi();
1876
1877 //[A,pTjet,pTtrig,dphi]
1878 tmpArrayFour[0] = areaJet;
1879 tmpArrayFour[1] = pTJet;
1880 tmpArrayFour[2] = triggerHadron->Pt();
1881 tmpArrayFour[3] = dPhi;
1882
1883 fHJetPhiCorr->Fill(tmpArrayFour); // Dphi distribution jet-triger
1884
1885 //[A,pTjet-pTUe,pTtrig,dphi]
1886 tmpArrayFour[1] = pTJet - ptUeFromCellMedian;
1887 fHJetPhiCorrSubUeMedian->Fill(tmpArrayFour);
1888
1889 //[A,pTjet-pTUe,pTtrig,dphi]
1890 tmpArrayFour[1] = pTJet - ptUeCone;
1891 fHJetPhiCorrSubUeCone->Fill(tmpArrayFour);
1892
1893 //Leticia's analysis
1894
8e103a56 1895 // Select back to back trigger - jet pairs
1896 if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue;
1897 fhDphiTriggerJetAccept->Fill(dphi); //Accepted
255cb71c 1898
1899 //NO bg subtraction
1900 //Centrality, A, pTjet, pTtrigg, dphi
1901 tmpArrayFive[0] = centValue;
1902 tmpArrayFive[1] = areaJet;
1903 tmpArrayFive[2] = pTJet;
1904 tmpArrayFive[3] = triggerHadron->Pt();
1905 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1906 fHJetSpec->Fill(tmpArrayFive);
1907
49dc2855 1908 //Centrality, A, pTjet-pTbgCellMedian, pTtrigg, dphi
255cb71c 1909 tmpArrayFive[0] = centValue;
1910 tmpArrayFive[1] = areaJet;
8aa19603 1911 tmpArrayFive[2] = pTJet-ptUeFromCellMedian;
255cb71c 1912 tmpArrayFive[3] = triggerHadron->Pt();
1913 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1914 fHJetSpecSubUeMedian->Fill(tmpArrayFive);
1915
8aa19603 1916 //Centrality, A, pTjet-pTbgCone, pTtrigg, dphi
255cb71c 1917 tmpArrayFive[0] = centValue;
1918 tmpArrayFive[1] = areaJet;
8aa19603 1919 tmpArrayFive[2] = pTJet - ptUeCone;
255cb71c 1920 tmpArrayFive[3] = triggerHadron->Pt();
1921 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
8aa19603 1922 fHJetSpecSubUeCone->Fill(tmpArrayFive);
255cb71c 1923
1924 //Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", kT median
1925 tmpArrayFive[0] = areaJet;
1926 tmpArrayFive[1] = pTJet;
8aa19603 1927 tmpArrayFive[2] = pTJet - ptUeFromCellMedian;
1928 tmpArrayFive[3] = ptUeFromCellMedian;
1929 tmpArrayFive[4] = rhoFromCellMedian;
255cb71c 1930 fHJetUeMedian->Fill(tmpArrayFive);
8e103a56 1931
8aa19603 1932 //Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", Cone median
255cb71c 1933 tmpArrayFive[0] = areaJet;
1934 tmpArrayFive[1] = pTJet;
8aa19603 1935 tmpArrayFive[2] = pTJet - ptUeCone;
1936 tmpArrayFive[3] = ptUeCone;
1937 tmpArrayFive[4] = rhoCone;
1938 fHJetUeCone->Fill(tmpArrayFive);
255cb71c 1939
1940 if(filledOnce){ //fill for each event only once
8aa19603 1941 Double_t fillRho[] = { rhoCone,rhoFromCellMedian};
1942 fHRhoUeMedianVsCone->Fill(fillRho);
255cb71c 1943 filledOnce = kFALSE;
1944 }
8e103a56 1945 }//jet loop
1946
1947 //Fill Jet Density In the Event A>0.07
255cb71c 1948 /*if(injet>0){
8e103a56 1949 Double_t filldens[]={ (Double_t) particleList.GetEntries(),
ad869500 1950 injet/fkAcceptance,
1951 triggerHadron->Pt()
1952 };
8e103a56 1953 fHJetDensity->Fill(filldens);
1954 }
ad869500 1955
8e103a56 1956 //Fill Jet Density In the Event A>0.4
1957 if(injet4>0){
1958 Double_t filldens4[]={ (Double_t) particleList.GetEntries(),
1959 injet4/fkAcceptance,
1960 triggerHadron->Pt()
1961 };
1962 fHJetDensityA4->Fill(filldens4);
255cb71c 1963 }*/
ad869500 1964 }
1965
8e103a56 1966
ad869500 1967 PostData(1, fOutputList);
1968}
1969
1970//----------------------------------------------------------------------------
1971void AliAnalysisTaskJetCorePP::Terminate(const Option_t *)
1972{
1973 // Draw result to the screen
1974 // Called once at the end of the query
1975
1976 if(fDebug) printf("AliAnalysisTaskJetCorePP DONE\n");
1977 if(!GetOutputData(1)) return;
1978}
1979
1980
1981//----------------------------------------------------------------------------
1982Int_t AliAnalysisTaskJetCorePP::GetListOfTracks(TList *list){
1983 //Fill the list of accepted tracks (passed track cut)
1984 //return consecutive index of the hardest ch hadron in the list
1985 Int_t iCount = 0;
1986 AliAODEvent *aodevt = NULL;
1987
1988 if(!fESD) aodevt = fAODIn;
1989 else aodevt = fAODOut;
1990
1991 if(!aodevt) return -1;
1992
1993 Int_t index = -1; //index of the highest particle in the list
1994 Double_t ptmax = -10;
8e103a56 1995 Int_t triggers[200];
1996 Int_t ntriggers = 0; //index in triggers array
1e3ad34d 1997
ad869500 1998
80ac66f6 1999 for(Int_t it = 0; it < aodevt->GetNumberOfTracks(); it++){
0fe881ea 2000 //AliAODTrack *tr = aodevt->GetTrack(it);
f15c1f69 2001 AliAODTrack *tr = dynamic_cast<AliAODTrack*>(aodevt->GetTrack(it));
2002 if(!tr) AliFatal("Not a standard AOD");
0fe881ea 2003
d1405a52
ML		 2004 if((fFilterMask > 0) && !(tr->TestFilterBit(fFilterMask))) continue;
2005 //if((fFilterMask > 0) && !(tr->IsHybridGlobalConstrainedGlobal())) continue;
ad869500 2006 if(TMath::Abs((Float_t) tr->Eta()) > fTrackEtaCut) continue;
2007 if(tr->Pt() < fTrackLowPtCut) continue;
2008 list->Add(tr);
8e103a56 2009
2010 //Search trigger:
2011 if(fHardest>0){ //leading particle
2012 if(tr->Pt()>ptmax){
2013 ptmax = tr->Pt();
2014 index = iCount;
2015 }
2016 }
2017
2018 if(fHardest==0 && ntriggers<200){ //single inclusive
2019 if(fTriggerPtRangeLow <= tr->Pt() &&
2020 tr->Pt() < fTriggerPtRangeHigh){
2021 triggers[ntriggers] = iCount;
2022 ntriggers++;
2023 }
ad869500 2024 }
8e103a56 2025
ad869500 2026 iCount++;
2027 }
2028
8e103a56 2029 if(fHardest==0 && ntriggers>0){ //select random inclusive trigger
2030 Int_t rnd = fRandom->Integer(ntriggers); //0 to ntriggers-1
2031 index = triggers[rnd];
1e3ad34d 2032
2033 fhNofMultipleTriggers->Fill(ntriggers-1);
ac4f24f7 2034
2035 Double_t deltaPhi, deltaEta, deltaR;
30557bf5 2036 Int_t iLow=0;
2037 Int_t iHigh=0;
ac4f24f7 2038 //Correlation with single inclusive trigger
2039 AliVParticle* tGent1 = (AliVParticle*) list->At(index);
2040 if(tGent1){
58003b11 2041 //for(Int_t ia=0; ia<ntriggers; ia++)
2042 for(Int_t ia=0; ia<list->GetEntries(); ia++){
2043 //if(triggers[ia]==index) continue;
2044 if(ia==index) continue;
2045 //AliVParticle* tGent2 = (AliVParticle*) list->At(triggers[ia]);
2046 AliVParticle* tGent2 = (AliVParticle*) list->At(ia);
ac4f24f7 2047 if(!tGent2) continue;
58003b11 2048 if(tGent2->Pt()<15.0) continue;
30557bf5 2049 if(tGent2->Pt()>50.0) continue;
ac4f24f7 2050 deltaPhi = RelativePhi(tGent1->Phi(),tGent2->Phi());
2051 deltaEta = tGent1->Eta()-tGent2->Eta();
2052 deltaR = sqrt(deltaPhi*deltaPhi + deltaEta*deltaEta);
30557bf5 2053
2054 fhDeltaRMultTriggers->Fill(deltaR);
2055 fhDeltaPhiMultTriggers->Fill(deltaPhi);
2056
2057 if(deltaR<0.4){
2058 if(tGent2->Pt()<20.0) iLow++;
2059 else iHigh++;
1e3ad34d 2060 }
30557bf5 2061 }
1e3ad34d 2062 }
30557bf5 2063 fhNofMultipleTriggersCone->Fill(iLow+iHigh);
2064 fhNofMultipleTriggersConeLow->Fill(iLow);
2065 fhNofMultipleTriggersConeHigh->Fill(iHigh);
2066
ad869500 2067 }
8e103a56 2068
2069 return index;
ad869500 2070}
2071
2072//----------------------------------------------------------------------------
80ac66f6 2073/*
ad869500 2074Double_t AliAnalysisTaskJetCorePP::GetBackgroundInPerpCone(Float_t jetR, Double_t jetPhi, Double_t jetEta, TList* trkList){
2075 //calculate sum of track pT in the cone perpendicular in phi to the jet
2076 //jetR = cone radius
2077 // jetPhi, jetEta = direction of the jet
2078 Int_t numberOfTrks = trkList->GetEntries();
2079 Double_t pTsum = 0.0;
2080 Double_t perpConePhi = jetPhi + TMath::Pi()/2;//perp cone w.r.t. jet in phi
2081 for(Int_t it=0; it<numberOfTrks; it++){
2082 AliVParticle *trk = (AliVParticle*) trkList->At(it);
2083 Double_t dphi = RelativePhi(perpConePhi,trk->Phi());
2084 Double_t deta = trk->Eta()-jetEta;
2085
2086 if( (dphi*dphi + deta*deta)< (jetR*jetR)){
2087 pTsum += trk->Pt();
2088 }
2089 }
2090
2091 return pTsum;
2092}
80ac66f6 2093*/
ad869500 2094//----------------------------------------------------------------------------
2095
2096Double_t AliAnalysisTaskJetCorePP::RelativePhi(Double_t mphi,Double_t vphi){
2097 //Get relative azimuthal angle of two particles -pi to pi
2098 if (vphi < -TMath::Pi()) vphi += TMath::TwoPi();
2099 else if (vphi > TMath::Pi()) vphi -= TMath::TwoPi();
2100
2101 if (mphi < -TMath::Pi()) mphi += TMath::TwoPi();
2102 else if (mphi > TMath::Pi()) mphi -= TMath::TwoPi();
2103
2104 Double_t dphi = mphi - vphi;
2105 if (dphi < -TMath::Pi()) dphi += TMath::TwoPi();
2106 else if (dphi > TMath::Pi()) dphi -= TMath::TwoPi();
2107
2108 return dphi;//dphi in [-Pi, Pi]
2109}
2110
2111
80ac66f6 2112//----------------------------------------------------------------------------
2113Bool_t AliAnalysisTaskJetCorePP::SelectMCGenTracks(AliVParticle *trk, TList *trkList, Double_t &ptLeading, Int_t &index, Int_t counter){
2114 //fill track efficiency denominator
2115 if(trk->Pt() < fTrackLowPtCut) return kFALSE;
2116 if(trk->Charge()==0) return kFALSE;
2117 if(TMath::Abs((Float_t) trk->Eta()) > fTrackEtaCut) return kFALSE;
2118 trkList->Add(trk);
8aa19603 2119 if(fFillRespMx) fhPtTrkTruePrimGen->Fill(trk->Pt(),trk->Eta()); //Efficiency denominator
8e103a56 2120
2121 //Search trigger:
2122 if(fHardest>0){ //leading particle
2123 if(ptLeading < trk->Pt()){
2124 index = counter;
2125 ptLeading = trk->Pt();
2126 }
80ac66f6 2127 }
8e103a56 2128
2129 if(fHardest==0){ //single inclusive
2130 index = -1;
2131 if(fTriggerPtRangeLow <= trk->Pt() &&
2132 trk->Pt() < fTriggerPtRangeHigh){
2133 index = counter;
2134 }
2135 }
2136
80ac66f6 2137 return kTRUE;
2138}
2139
2140//----------------------------------------------------------------------------
2141void AliAnalysisTaskJetCorePP::FillEffHistos(TList *recList, TList *genList){
2142 //fill track efficiency numerator
2143 if(!recList) return;
2144 if(!genList) return;
2145 Int_t nRec = recList->GetEntries();
2146 Int_t nGen = genList->GetEntries();
2147 for(Int_t ir=0; ir<nRec; ir++){
2148 AliAODTrack *trkRec = (AliAODTrack*) recList->At(ir);
2149 if(!trkRec) continue;
2150 Int_t recLabel = TMath::Abs(trkRec->GetLabel());
2151 Bool_t matched = kFALSE;
2152
2153 for(Int_t ig=0; ig<nGen; ig++){
2154 AliVParticle *trkGen = (AliVParticle*) genList->At(ig);
2155 if(!trkGen) continue;
2156 Int_t genLabel = TMath::Abs(trkGen->GetLabel());
2157 if(recLabel==genLabel){
2158 fhPtTrkTruePrimRec->Fill(trkGen->Pt(),trkGen->Eta());
2159 matched = kTRUE;
2160 break;
2161 }
2162 }
2163 if(!matched) fhPtTrkSecOrFakeRec->Fill(trkRec->Pt(),trkRec->Eta());
2164 }
2165
2166 return;
2167}
255cb71c 2168//________________________________________________________________
d1405a52 2169void AliAnalysisTaskJetCorePP::EstimateBgRhoMedian(TList *listJet, TList* listPart, Double_t &rhoMedian, Int_t mode){
8aa19603 2170 //Estimate background rho by means of integrating track pT outside identified jet cones
2171
2172 rhoMedian = 0;
2173
2174 //identify leading jet in the full track acceptance
2175 Int_t idxLeadingJet = -1;
2176 Double_t pTleading = 0.0;
2177 AliAODJet* jet = NULL;
2178
2179 for(Int_t ij = 0; ij < listJet->GetEntries(); ij++){ //loop over bg kT jets
2180 jet = (AliAODJet*)(listJet->At(ij));
2181 if(!jet) continue;
2182 if(TMath::Abs(jet->Eta()) > fTrackEtaCut) continue; //track acceptance cut
2183 if(jet->Pt() > pTleading){
2184 idxLeadingJet = ij;
2185 pTleading = jet->Pt();
2186 }
2187 }
2188
2189 Int_t njets = 0;
2190 static Double_t jphi[100];
2191 static Double_t jeta[100];
2192 static Double_t jRsquared[100];
2193
2194 for(Int_t ij = 0; ij< listJet->GetEntries(); ij++){ //loop over bg kT jets
2195
2196 jet = (AliAODJet*)(listJet->At(ij));
2197 if(!jet) continue;
2198 if(TMath::Abs(jet->Eta()) > fTrackEtaCut) continue;
2199
2200 if((jet->Pt() > fBgMaxJetPt || ij== idxLeadingJet) && njets<100){
2201 jphi[njets] = RelativePhi(jet->Phi(),0.0); //-pi,pi
2202 jeta[njets] = jet->Eta();
2203 jRsquared[njets] = fSafetyMargin*jet->EffectiveAreaCharged()/TMath::Pi(); //scale jet radius
2204 njets++;
255cb71c 2205 }
2206 }
8aa19603 2207
255cb71c 2208
8aa19603 2209 static Double_t nOutCone[10][4];
2210 static Double_t sumPtOutOfCone[10][4];
2211
2212 for(Int_t ie=0; ie<fnEta; ie++){
2213 for(Int_t ip=0; ip<fnPhi; ip++){
2214 nOutCone[ip][ie] = 0.0; //initialize counter
2215 sumPtOutOfCone[ip][ie] = 0.0;
255cb71c 2216 }
8aa19603 2217 }
255cb71c 2218
8aa19603 2219 Double_t rndphi, rndeta;
2220 Double_t rndphishift, rndetashift;
2221 Double_t dphi, deta;
2222 Bool_t bIsInCone;
2223
2224 for(Int_t it=0; it<fnTrials; it++){
2225
2226 rndphi = fRandom->Uniform(0, fPhiSize);
2227 rndeta = fRandom->Uniform(0, fEtaSize);
2228
2229 for(Int_t ip=0; ip<fnPhi; ip++){ //move radom position to each cell
2230 rndphishift = rndphi + ip*fPhiSize - TMath::Pi();
2231 for(Int_t ie=0; ie<fnEta; ie++){
2232 rndetashift = rndeta + ie*fEtaSize - fEtaSize;
2233
2234 bIsInCone = 0; //tag if trial is in the jet cone
2235 for(Int_t ij=0; ij<njets; ij++){
2236 deta = jeta[ij] - rndetashift;
2237 dphi = RelativePhi(rndphishift,jphi[ij]);
2238 if((dphi*dphi + deta*deta) < jRsquared[ij]){
2239 bIsInCone = 1;
2240 break;
2241 }
2242 }
2243 if(!bIsInCone) nOutCone[ip][ie]++;
2244 }
2245 }
2246 }
2247
255cb71c 2248
8aa19603 2249 //Sum particle pT outside jets in cells
2250 Int_t npart = listPart->GetEntries();
2251 Int_t phicell,etacell;
2252 AliVParticle *part = NULL;
2253 for(Int_t ip=0; ip < npart; ip++){
2254
2255 part = (AliVParticle*) listPart->At(ip);
2256 if(!part) continue;
2257 if( TMath::Abs(part->Eta()) > fEtaSize) continue; //skip part outside +-0.7 in eta
2258
2259 bIsInCone = 0; //init
2260 for(Int_t ij=0; ij<njets; ij++){
2261 dphi = RelativePhi(jphi[ij], part->Phi());
2262 deta = jeta[ij] - part->Eta();
2263 if((dphi*dphi + deta*deta) < jRsquared[ij]){
2264 bIsInCone = 1;
2265 break;
2266 }
255cb71c 2267 }
8aa19603 2268 if(!bIsInCone){
2269 phicell = TMath::Nint(TMath::Floor((RelativePhi(part->Phi(),0.0) + TMath::Pi())/fPhiSize));
2270 etacell = TMath::Nint(TMath::Floor((part->Eta()+fEtaSize)/fEtaSize));
2271 sumPtOutOfCone[phicell][etacell]+= part->Pt();
2272 }
2273 }
2274
2275 static Double_t rhoInCells[20];
2276 Double_t relativeArea;
2277 Int_t nCells=0;
d1405a52 2278 Double_t bufferArea=0.0, bufferPt=0.0; //sum cells where A< fJetFreeAreaFrac
8aa19603 2279 for(Int_t ip=0; ip<fnPhi; ip++){
2280 for(Int_t ie=0; ie<fnEta; ie++){
2281 relativeArea = nOutCone[ip][ie]/fnTrials;
d1405a52
ML		 2282 //cout<<"RA= "<< relativeArea<<" BA="<<bufferArea<<" BPT="<< bufferPt <<endl;
2283
2284 bufferArea += relativeArea;
2285 bufferPt += sumPtOutOfCone[ip][ie];
2286 if(bufferArea > fJetFreeAreaFrac){
2287 rhoInCells[nCells] = bufferPt/(bufferArea*fCellArea);
2288
2289 if(mode==1) fhCellAreaToMedianGen->Fill(bufferArea*fCellArea);
2290 else fhCellAreaToMedian->Fill(bufferArea*fCellArea);
2291
2292 bufferArea = 0.0;
2293 bufferPt = 0.0;
2294 nCells++;
2295 }
8aa19603 2296 }
255cb71c 2297 }
8aa19603 2298
d1405a52
ML		 2299
2300 if(nCells>0){
2301 rhoMedian = TMath::Median(nCells, rhoInCells);
2302 if(mode==1){ //mc data
2303 fhEntriesToMedianGen->Fill(nCells);
2304 }else{ //real data
2305 fhEntriesToMedian->Fill(nCells);
2306 }
2307 }
8aa19603 2308
255cb71c 2309}
80ac66f6 2310
255cb71c 2311//__________________________________________________________________
8aa19603 2312void AliAnalysisTaskJetCorePP::EstimateBgCone(TList *listJet, TList* listPart, Double_t &rhoPerpCone){
2313
2314 rhoPerpCone = 0.0; //init
2315
2316 if(!listJet) return;
2317 Int_t njet = listJet->GetEntries();
2318 Int_t npart = listPart->GetEntries();
2319 Double_t pTleading = 0.0;
2320 Double_t phiLeading = 1000.;
2321 Double_t etaLeading = 1000.;
2322
2323 for(Int_t jsig=0; jsig < njet; jsig++){
2324 AliAODJet* signaljet = (AliAODJet*)(listJet->At(jsig));
2325 if(!signaljet) continue;
2326 if((signaljet->Eta()<fJetEtaMin) && (fJetEtaMax<signaljet->Eta())) continue; //acceptance cut
2327 if(signaljet->Pt() >= pTleading){ //replace leading and subleading jet
2328 pTleading = signaljet->Pt();
2329 phiLeading = signaljet->Phi();
2330 etaLeading = signaljet->Eta();
2331 }
2332 }
2333
2334 Double_t phileftcone = phiLeading + TMath::Pi()/2;
2335 Double_t phirightcone = phiLeading - TMath::Pi()/2;
2336 Double_t dp, de;
2337
2338 for(Int_t ip=0; ip < npart; ip++){
2339
2340 AliVParticle *part = (AliVParticle*) listPart->At(ip);
2341 if(!part){
2342 continue;
2343 }
2344
2345 dp = RelativePhi(phileftcone, part->Phi());
2346 de = etaLeading - part->Eta();
2347 if( sqrt(dp*dp + de*de)< fBgConeR) rhoPerpCone += part->Pt();
2348
2349 dp = RelativePhi(phirightcone, part->Phi());
2350 if( sqrt(dp*dp + de*de)< fBgConeR) rhoPerpCone += part->Pt();
2351
2352 }
2353
2354 //normalize total pT by two times cone are
2355 rhoPerpCone = rhoPerpCone/(2*TMath::Pi()*fBgConeR*fBgConeR);
2356
2357}
2358//__________________________________________________________________
80ac66f6 2359
255cb71c 2360void AliAnalysisTaskJetCorePP::ReadTClonesArray(TString bname, TList *list){
2361 //Convert TClones array of jets to TList
2362
2363 if(!strlen(bname.Data())){
2364 AliError(Form("Jet branch %s not set.", bname.Data()));
2365 return;
2366 }
2367
2368 TClonesArray *array=0x0;
49dc2855 2369 if(fUseExchContainer){ //take input from exchange containers
2370 if(fJetBranchNameKine.Length()>0 && bname.CompareTo(fJetBranchNameKine.Data())==0){
2371 array = dynamic_cast<TClonesArray*>(GetInputData(1)); //connect exchange container slot 1
2372 }
3d1fac8a 2373 if(fJetBranchNameBgKine.Length()>0 && bname.CompareTo(fJetBranchNameBgKine.Data())==0){
2374 array = dynamic_cast<TClonesArray*>(GetInputData(2)); //connect exchange container slot 2
2375 }
49dc2855 2376 }else{ //take input from AOD
2377 if(fAODOut&&!array){
2378 array = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(bname.Data()));
2379 }
2380 if(fAODExtension&&!array){
2381 array = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(bname.Data()));
2382 }
2383 if(fAODIn&&!array){
2384 array = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(bname.Data()));
2385 }
255cb71c 2386 }
2387
2388 list->Clear(); //clear list beforehand
2389
2390 if(array){
2391 if(fDebug){
2392 Printf("########## %s: %d jets",bname.Data(), array->GetEntriesFast());
2393 }
49dc2855 2394
255cb71c 2395 for(Int_t iJet = 0; iJet < array->GetEntriesFast(); iJet++) {
2396 AliAODJet *jet = dynamic_cast<AliAODJet*>((*array)[iJet]);
2397 if (jet) list->Add(jet);
2398 }
2399 }
2400
2401 return;
2402}