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