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