]> 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
Updated TPC preprocessor to include gas composition entries
[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);
49dc2855 837 //---
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
255cb71c 850 //---
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
ac4f24f7 1302 Double_t deltaPhi, deltaEta, deltaR;
1303 Int_t k = 0;
1304
1305 //Correlation with single inclusive TRIGGER
1306 AliVParticle* tGenT1 = (AliVParticle*) particleListGen.At(indexTriggGen);
1307 if(tGenT1){
1308 for(Int_t ia=0; ia<ntriggersMC; ia++){
1309 if(indexTriggGen == triggersMC[ia]) continue;
1310
1311 AliVParticle* tGenT2 = (AliVParticle*) particleListGen.At(triggersMC[ia]);
1312 if(!tGenT2) continue;
1313
1314 deltaPhi = RelativePhi(tGenT1->Phi(),tGenT2->Phi());
1315 deltaEta = tGenT1->Eta()-tGenT2->Eta();
1316 deltaR = sqrt(deltaPhi*deltaPhi + deltaEta*deltaEta);
0fe881ea 1317
ac4f24f7 1318 if(deltaR<0.4) k++;
0fe881ea 1319 }
ac4f24f7 1320 }
1321 fhNofMultipleTriggersConeGen->Fill(k);
1322
1323 if(ntriggersMC>1){
0fe881ea 1324 //Correlation of each trigger with any other trigger
1e3ad34d 1325 for(Int_t ia=0; ia<ntriggersMC-1; ia++){
0fe881ea 1326 AliVParticle* tGenI = (AliVParticle*) particleListGen.At(triggersMC[ia]);
1e3ad34d 1327 if(!tGenI) continue;
1328 for(Int_t ib=ia+1; ib<ntriggersMC; ib++){
0fe881ea 1329 AliVParticle* tGenII = (AliVParticle*) particleListGen.At(triggersMC[ib]);
1e3ad34d 1330 if(!tGenII) continue;
1331
1332 deltaPhi = RelativePhi(tGenI->Phi(),tGenII->Phi());
1333 deltaEta = tGenI->Eta()-tGenII->Eta();
1334 deltaR = sqrt(deltaPhi*deltaPhi + deltaEta*deltaEta);
1335 fhDeltaRMultTriggersGen->Fill(deltaR);
1336 }
1337 }
1338 }
1339
80ac66f6 1340 }else{
8e103a56 1341 indexTriggGen = -1; //trigger not found
80ac66f6 1342 }
1343 }
8e103a56 1344
255cb71c 1345 //-----------
8e103a56 1346 Int_t ilowGen = (fHardest==0 || fHardest==1) ? indexTriggGen : 0;
1347 Int_t ihighGen = (fHardest==0 || fHardest==1) ? indexTriggGen+1 : particleListGen.GetEntries();
255cb71c 1348 Bool_t fillOnceGen = kTRUE;
1349 //-----------
1350
8e103a56 1351 for(Int_t igen= ilowGen; igen< ihighGen; igen++){ //loop over possible trigger
1352 indexTriggGen = igen; //trigger hadron
1353
1354 if(indexTriggGen == -1) continue;
1355 AliVParticle* triggerGen=(AliVParticle*) particleListGen.At(indexTriggGen);
1356 if(!triggerGen) continue;
1357
1358 if(fHardest >= 2){
6168afc3 1359 if(triggerGen->Pt() < 6.0) continue; //all hadrons pt>6
8e103a56 1360 }
1361 if(TMath::Abs((Float_t) triggerGen->Eta()) > fTriggerEtaCut) continue;
1362
1363 ptTriggGen = triggerGen->Pt(); //count triggers
1364 fh2NtriggersGen->Fill(centValue, ptTriggGen);
1365
80ac66f6 1366 //Count jets and trigger-jet pairs at MC generator level
1367 for(Int_t ij=0; ij<fListJetsGen->GetEntries(); ij++){
1368 AliAODJet* jet = (AliAODJet*)(fListJetsGen->At(ij));
1369 if(!jet) continue;
1370 Double_t etaJetGen = jet->Eta();
255cb71c 1371 Double_t areaJetGen = jet->EffectiveAreaCharged();
80ac66f6 1372
1373 if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){
80ac66f6 1374
49dc2855 1375 Double_t ptUeFromCellMedianGen = rhoFromCellMedianGen*areaJetGen;
1376 Double_t ptUeConeGen = rhoConeGen*areaJetGen;
1377
1378 //Rongrong's analysis
1379 Double_t dPhiGen = jet->Phi() - triggerGen->Phi();
1380 if(dPhiGen>2*TMath::Pi()) dPhiGen -= 2*TMath::Pi();
1381 if(dPhiGen<-2*TMath::Pi()) dPhiGen += 2*TMath::Pi();
1382 if(dPhiGen<-0.5*TMath::Pi()) dPhiGen += 2*TMath::Pi();
1383 if(dPhiGen>1.5*TMath::Pi()) dPhiGen -= 2*TMath::Pi();
1384
1385 //[A,pTjet,pTtrig,dphi]
1386 tmpArrayFour[0] = areaJetGen;
1387 tmpArrayFour[1] = jet->Pt();
1388 tmpArrayFour[2] = ptTriggGen;
1389 tmpArrayFour[3] = dPhiGen;
1390
1391 fHJetPhiCorrGen->Fill(tmpArrayFour); // Dphi distribution jet-triger
1392
1393 //[A,pTjet-pTUe,pTtrig,dphi]
1394 tmpArrayFour[1] = jet->Pt() - ptUeFromCellMedianGen;
1395 fHJetPhiCorrSubUeMedianGen->Fill(tmpArrayFour);
1396
1397 //[A,pTjet-pTUe,pTtrig,dphi]
1398 tmpArrayFour[1] = jet->Pt() - ptUeConeGen;
1399 fHJetPhiCorrSubUeConeGen->Fill(tmpArrayFour);
1400
1401 //Leticia's analysis
8e103a56 1402 Double_t dphi = RelativePhi(triggerGen->Phi(), jet->Phi());
1403 if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue;
80ac66f6 1404
8e103a56 1405 //Centrality, A, pT, pTtrigg
255cb71c 1406 tmpArrayFive[0] = centValue;
1407 tmpArrayFive[1] = areaJetGen;
1408 tmpArrayFive[2] = jet->Pt();
1409 tmpArrayFive[3] = ptTriggGen;
1410 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1411 fHJetSpecGen->Fill(tmpArrayFive);
1412
8aa19603 1413 //Centrality, A, pTjet-pTbgCellMedian, pTtrigg, dphi
255cb71c 1414 tmpArrayFive[0] = centValue;
1415 tmpArrayFive[1] = areaJetGen;
8aa19603 1416 tmpArrayFive[2] = jet->Pt() - ptUeFromCellMedianGen;
255cb71c 1417 tmpArrayFive[3] = ptTriggGen;
1418 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1419 fHJetSpecSubUeMedianGen->Fill(tmpArrayFive);
1420
8aa19603 1421 //Centrality, A, pTjet-pTbgCone, pTtrigg, dphi
255cb71c 1422 tmpArrayFive[0] = centValue;
1423 tmpArrayFive[1] = areaJetGen;
8aa19603 1424 tmpArrayFive[2] = jet->Pt() - ptUeConeGen;
255cb71c 1425 tmpArrayFive[3] = ptTriggGen;
1426 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
8aa19603 1427 fHJetSpecSubUeConeGen->Fill(tmpArrayFive);
255cb71c 1428
1429 //Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", kT median
1430 tmpArrayFive[0] = areaJetGen;
1431 tmpArrayFive[1] = jet->Pt();
8aa19603 1432 tmpArrayFive[2] = jet->Pt() - ptUeFromCellMedianGen;
1433 tmpArrayFive[3] = ptUeFromCellMedianGen;
1434 tmpArrayFive[4] = rhoFromCellMedianGen;
255cb71c 1435 fHJetUeMedianGen->Fill(tmpArrayFive);
1436
8aa19603 1437 //Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", perpendicular Cone
255cb71c 1438 tmpArrayFive[0] = areaJetGen;
1439 tmpArrayFive[1] = jet->Pt();
8aa19603 1440 tmpArrayFive[2] = jet->Pt() - ptUeConeGen;
1441 tmpArrayFive[3] = ptUeConeGen;
1442 tmpArrayFive[4] = rhoConeGen;
1443 fHJetUeConeGen->Fill(tmpArrayFive);
255cb71c 1444
1445 if(fillOnceGen){
8aa19603 1446 Double_t fillRhoGen[] = { rhoConeGen, rhoFromCellMedianGen};
1447 fHRhoUeMedianVsConeGen->Fill(fillRhoGen);
255cb71c 1448 fillOnceGen = kFALSE;
1449 }
8e103a56 1450 }//back to back jet-trigger pair
1451 }//jet loop
1452 }//trigger loop
1453
49dc2855 1454 if(fFillRespMx && !fIsKine){
8aa19603 1455 //================ RESPONSE MATRIX ===============
1456 //Count jets and trigger-jet pairs at MC generator level
1457 for(Int_t ij=0; ij<fListJetsGen->GetEntries(); ij++){
1458 AliAODJet* jet = (AliAODJet*)(fListJetsGen->At(ij));
1459 if(!jet) continue;
1460 Double_t etaJetGen = jet->Eta();
1461 Double_t ptJetGen = jet->Pt();
255cb71c 1462
8aa19603 1463 if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){
1464 fhJetPtGen->Fill(ptJetGen); // gen level pt spectum of jets response mx normalization
80ac66f6 1465
8aa19603 1466 Double_t areaJetGen = jet->EffectiveAreaCharged();
1467 Double_t ptUeFromCellMedianGen = rhoFromCellMedianGen*areaJetGen;
1468 Double_t ptUeConeGen = rhoConeGen*areaJetGen;
80ac66f6 1469
8aa19603 1470 fhJetPtSubUeMedianGen->Fill(ptJetGen - ptUeFromCellMedianGen);
1471 fhJetPtSubUeConeGen->Fill(ptJetGen - ptUeConeGen);
160dbdb8 1472 }
1473 }
8aa19603 1474 if(fListJets->GetEntries()>0 && fListJetsGen->GetEntries()>0){ //at least some reconstructed jets
1475
1476 Int_t ng = (Int_t) fListJetsGen->GetEntries();
1477 Int_t nr = (Int_t) fListJets->GetEntries();
160dbdb8 1478
8aa19603 1479 //Find closest MC generator - reconstructed jets
1480 if(faGenIndex.GetSize()<nr) faGenIndex.Set(nr); //idx of gen jet assoc to rec jet
1481 if(faRecIndex.GetSize()<ng) faRecIndex.Set(ng); //idx of rec jet assoc to gen jet
160dbdb8 1482
1483 if(fDebug){
8aa19603 1484 Printf("New Rec List %d gen index Array %d",nr,faGenIndex.GetSize());
1485 Printf("New Gen List %d rec index Array %d",ng,faRecIndex.GetSize());
160dbdb8 1486 }
1487 //matching of MC genrator level and reconstructed jets
8aa19603 1488 AliAnalysisHelperJetTasks::GetClosestJets(fListJetsGen,ng,fListJets,nr,faGenIndex,faRecIndex,fDebug);
1489
1490 // Fill response matrix
1491 for(Int_t ir = 0; ir < nr; ir++){
1492 AliAODJet *recJet = (AliAODJet*) fListJets->At(ir);
1493 Double_t etaJetRec = recJet->Eta();
1494 Double_t ptJetRec = recJet->Pt();
1495 Double_t areaJetRec = recJet->EffectiveAreaCharged();
160dbdb8 1496 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1497
8aa19603 1498 if((fJetEtaMin <= etaJetRec) && (etaJetRec <= fJetEtaMax)){
1499 Int_t ig = faGenIndex[ir]; //associated generator level jet
1500 if(ig >= 0 && ig < ng){
1501 if(fDebug > 10) Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1502 AliAODJet *genJet = (AliAODJet*) fListJetsGen->At(ig);
1503 Double_t ptJetGen = genJet->Pt();
1504 Double_t etaJetGen = genJet->Eta();
160dbdb8 1505
255cb71c 1506 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
8aa19603 1507 if((fJetEtaMin <= etaJetGen) && (etaJetGen <= fJetEtaMax)){
1508 fhJetPtGenVsJetPtRec->Fill(ptJetRec, ptJetGen);
1509
1510 Double_t areaJetGen = genJet->EffectiveAreaCharged();
1511 Double_t ptUeFromCellMedianGen = rhoFromCellMedianGen*areaJetGen;
1512 Double_t ptUeConeGen = rhoConeGen*areaJetGen;
1513 Double_t ptUeFromCellMedianRec = rhoFromCellMedian*areaJetRec;
1514 Double_t ptUeConeRec = rhoCone*areaJetRec;
1515 fhJetPtGenVsJetPtRecSubUeMedian->Fill(ptJetRec-ptUeFromCellMedianRec,
1516 ptJetGen-ptUeFromCellMedianGen);
1517 fhJetPtGenVsJetPtRecSubUeCone->Fill(ptJetRec-ptUeConeRec, ptJetGen-ptUeConeGen);
160dbdb8 1518 }
8aa19603 1519 }//ig>=0
160dbdb8 1520 }//rec jet in eta acceptance
1521 }//loop over reconstructed jets
1522 }// # of rec jets >0
8aa19603 1523
1524 //=========================== Full jet vs charged jet matrix ==========================
1525 if(fIsFullMC){
1526 //Count full jets and charged-jet pairs at MC generator level
1527 for(Int_t ij=0; ij<fListJetsGenFull->GetEntries(); ij++){
1528 AliAODJet* jetFull = (AliAODJet*)(fListJetsGenFull->At(ij));
1529 if(!jetFull) continue;
1530 Double_t etaJetFull = jetFull->Eta();
1531 Double_t ptJetFull = jetFull->Pt();
1532
1533 if((fJetEtaMin<=etaJetFull) && (etaJetFull<=fJetEtaMax)){
1534 fhJetPtGenFull->Fill(ptJetFull); // generator level pt spectum of full jets
1535 }
1536 }
1537 if(fListJetsGen->GetEntries()>0 && fListJetsGenFull->GetEntries()>0){ //at least some reconstructed jets
1538 Int_t nful = (Int_t) fListJetsGenFull->GetEntries();
1539 Int_t nchr = (Int_t) fListJetsGen->GetEntries();
1540
1541 //Find closest MC generator full - charged jet
1542 if(faGenIndex.GetSize()<nchr) faGenIndex.Set(nchr); //idx of gen FULL jet assoc to gen CHARGED jet
1543 if(faRecIndex.GetSize()<nful) faRecIndex.Set(nful); //idx of gen CHARGED jet assoc to gen FULL jet
1544
1545 if(fDebug){
1546 Printf("New Charg List %d Full index Array %d",nchr,faGenIndex.GetSize());
1547 Printf("New Full List %d Charg index Array %d",nful,faRecIndex.GetSize());
1548 }
1549 //matching of MC genrator level and reconstructed jets
1550 AliAnalysisHelperJetTasks::GetClosestJets(fListJetsGenFull,nful,fListJetsGen,nchr,faGenIndex,faRecIndex,fDebug);
1551
1552 // Fill response matrix
1553 for(Int_t ichr = 0; ichr < nchr; ichr++){ //charged jet loop
1554 AliAODJet *chJet = (AliAODJet*) fListJetsGen->At(ichr);
1555 Double_t etaJetCh = chJet->Eta();
1556 Double_t ptJetCh = chJet->Pt();
1557 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1558
1559 if((fJetEtaMin <= etaJetCh) && (etaJetCh <= fJetEtaMax)){
1560 Int_t iful = faGenIndex[ichr]; //associated generator level jet
1561 if(iful >= 0 && iful < nful){
1562 if(fDebug > 10) Printf("%s:%d iful = %d ichr = %d",(char*)__FILE__,__LINE__,iful,ichr);
1563 AliAODJet *genJetFull = (AliAODJet*) fListJetsGenFull->At(iful);
1564 Double_t ptJetFull = genJetFull->Pt();
1565 Double_t etaJetFull = genJetFull->Eta();
1566
1567 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1568 if((fJetEtaMin <= etaJetFull) && (etaJetFull <= fJetEtaMax)){
1569 fhJetPtGenChargVsJetPtGenFull->Fill(ptJetFull,ptJetCh);
1570 }
1571 }//iful>=0
1572 }//rec jet in eta acceptance
1573 }//loop over reconstructed jets
1574 }// # of rec jets >0
1575 }//pointer MC generator jets
1576 } //fill resp mx only for bin
1577 }//analyze generator level MC
1578
1579
49dc2855 1580 if(fIsKine){ //skip reconstructed data analysis in case of kine
1581 PostData(1, fOutputList);
1582 return;
1583 }
8e103a56 1584 //============= RECONSTRUCTED INCLUSIVE JETS ===============
ad869500 1585
ad869500 1586 Double_t etaJet = 0.0;
1587 Double_t pTJet = 0.0;
1588 Double_t areaJet = 0.0;
1589 Double_t phiJet = 0.0;
9f10e189 1590 //Int_t indexLeadingJet = -1;
1591 //Double_t pTLeadingJet = -10.0;
1592 //Double_t areaLeadingJet = -10.0;
8e103a56 1593
ad869500 1594 for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
1595 AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
1596 if(!jet) continue;
1597 etaJet = jet->Eta();
1598 phiJet = jet->Phi();
1599 pTJet = jet->Pt();
1600 if(pTJet==0) continue;
1601
1602 if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
255cb71c 1603 /*areaJet = jet->EffectiveAreaCharged();*/
7fc3d134 1604
ad869500 1605 //Jet Diagnostics---------------------------------
8e103a56 1606 fhJetPhi->Fill(pTJet, RelativePhi(phiJet,0.0)); //phi -pi,pi
1607 fhJetEta->Fill(pTJet, etaJet);
7fc3d134 1608 //search for leading jet
255cb71c 1609 /*if(pTJet > pTLeadingJet){
7fc3d134 1610 indexLeadingJet = ij;
1611 pTLeadingJet = pTJet;
7fc3d134 1612 areaLeadingJet = areaJet;
1613 }
1614
1615 // raw spectra of INCLUSIVE jets
80ac66f6 1616 //Centrality, pTjet, A
7fc3d134 1617 Double_t fillraw[] = { centValue,
1618 pTJet,
7fc3d134 1619 areaJet
1620 };
255cb71c 1621 fHJetPtRaw->Fill(fillraw);*/
8e103a56 1622 }
255cb71c 1623 /*
8e103a56 1624 if(indexLeadingJet > -1){
1625 // raw spectra of LEADING jets
1626 //Centrality, pTjet, A
1627 Double_t fillleading[] = { centValue,
1628 pTLeadingJet,
1629 areaLeadingJet
1630 };
1631 fHLeadingJetPtRaw->Fill(fillleading);
1632 }
255cb71c 1633 */
8e103a56 1634
1635 // =============== RECONSTRUCTED TRIGGER-JET PAIRS ================
8aa19603 1636 if(fIsChargedMC && fFillRespMx){
1637 FillEffHistos(&particleList, &particleListGen); //Fill efficiency histos
1638 }
255cb71c 1639 Bool_t filledOnce = kTRUE; //fill rho histogram only once per event
8e103a56 1640 //set ranges of the trigger loop
1641 Int_t ilow = (fHardest==0 || fHardest==1) ? indexTrigg : 0;
1642 Int_t ihigh = (fHardest==0 || fHardest==1) ? indexTrigg+1 : particleList.GetEntries();
7fc3d134 1643
8e103a56 1644 for(Int_t itrk= ilow; itrk< ihigh; itrk++){ //loop over possible trigger
1645 indexTrigg = itrk; //trigger hadron with pT >10 GeV
7fc3d134 1646
8e103a56 1647 if(indexTrigg < 0) continue;
1648
1649 AliVParticle *triggerHadron = (AliVParticle*) particleList.At(indexTrigg);
1650 if(!triggerHadron){
1651 PostData(1, fOutputList);
1652 return;
1653 }
1654 if(fHardest >= 2){
7d9fd386 1655 if(triggerHadron->Pt() < 6.0) continue; //all hadrons pt>6
8e103a56 1656 }
1657 if(TMath::Abs((Float_t) triggerHadron->Eta())> fTriggerEtaCut) continue;
7fc3d134 1658
8e103a56 1659 //Fill trigger histograms
1660 fh2Ntriggers->Fill(centValue,triggerHadron->Pt()); //trigger pT
1661 fhTriggerPhi->Fill(triggerHadron->Pt(),RelativePhi(triggerHadron->Phi(),0.0)); //phi -pi,pi
1662 fhTriggerEta->Fill(triggerHadron->Pt(),triggerHadron->Eta());
7fc3d134 1663
1664
8e103a56 1665 //---------- make trigger-jet pairs ---------
255cb71c 1666 //Int_t injet4 = 0;
1667 //Int_t injet = 0;
255cb71c 1668
8e103a56 1669 for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
1670 AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
1671 if(!jet) continue;
1672 etaJet = jet->Eta();
1673 phiJet = jet->Phi();
1674 pTJet = jet->Pt();
1675 if(pTJet==0) continue;
1676
1677 if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
1678 areaJet = jet->EffectiveAreaCharged();
49dc2855 1679
1680 //subtract bg using estinates base on median of kT jets
1681 Double_t ptUeFromCellMedian = rhoFromCellMedian*areaJet;
1682 Double_t ptUeCone = rhoCone*areaJet;
1683
255cb71c 1684 //if(areaJet >= 0.07) injet++;
1685 //if(areaJet >= 0.4) injet4++;
8e103a56 1686
1687 Double_t dphi = RelativePhi(triggerHadron->Phi(), phiJet);
1688 fhDphiTriggerJet->Fill(dphi); //Input
1689
1690 //Dphi versus jet pT
1691 //Centrality, Dphi=phiTrig-phiJet, pTjet, pTtrigg
255cb71c 1692 /*Double_t filldp[] = { centValue,
8e103a56 1693 dphi,
1694 pTJet,
1695 triggerHadron->Pt()
1696 };
1697 fHDphiVsJetPtAll->Fill(filldp);
255cb71c 1698 */
49dc2855 1699 //Rongrong's analysis
1700
1701 Double_t dPhi = phiJet - triggerHadron->Phi();
1702 if(dPhi>2*TMath::Pi()) dPhi -= 2*TMath::Pi();
1703 if(dPhi<-2*TMath::Pi()) dPhi += 2*TMath::Pi();
1704 if(dPhi<-0.5*TMath::Pi()) dPhi += 2*TMath::Pi();
1705 if(dPhi>1.5*TMath::Pi()) dPhi -= 2*TMath::Pi();
1706
1707 //[A,pTjet,pTtrig,dphi]
1708 tmpArrayFour[0] = areaJet;
1709 tmpArrayFour[1] = pTJet;
1710 tmpArrayFour[2] = triggerHadron->Pt();
1711 tmpArrayFour[3] = dPhi;
1712
1713 fHJetPhiCorr->Fill(tmpArrayFour); // Dphi distribution jet-triger
1714
1715 //[A,pTjet-pTUe,pTtrig,dphi]
1716 tmpArrayFour[1] = pTJet - ptUeFromCellMedian;
1717 fHJetPhiCorrSubUeMedian->Fill(tmpArrayFour);
1718
1719 //[A,pTjet-pTUe,pTtrig,dphi]
1720 tmpArrayFour[1] = pTJet - ptUeCone;
1721 fHJetPhiCorrSubUeCone->Fill(tmpArrayFour);
1722
1723 //Leticia's analysis
1724
8e103a56 1725 // Select back to back trigger - jet pairs
1726 if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue;
1727 fhDphiTriggerJetAccept->Fill(dphi); //Accepted
255cb71c 1728
1729 //NO bg subtraction
1730 //Centrality, A, pTjet, pTtrigg, dphi
1731 tmpArrayFive[0] = centValue;
1732 tmpArrayFive[1] = areaJet;
1733 tmpArrayFive[2] = pTJet;
1734 tmpArrayFive[3] = triggerHadron->Pt();
1735 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1736 fHJetSpec->Fill(tmpArrayFive);
1737
49dc2855 1738 //Centrality, A, pTjet-pTbgCellMedian, pTtrigg, dphi
255cb71c 1739 tmpArrayFive[0] = centValue;
1740 tmpArrayFive[1] = areaJet;
8aa19603 1741 tmpArrayFive[2] = pTJet-ptUeFromCellMedian;
255cb71c 1742 tmpArrayFive[3] = triggerHadron->Pt();
1743 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1744 fHJetSpecSubUeMedian->Fill(tmpArrayFive);
1745
8aa19603 1746 //Centrality, A, pTjet-pTbgCone, pTtrigg, dphi
255cb71c 1747 tmpArrayFive[0] = centValue;
1748 tmpArrayFive[1] = areaJet;
8aa19603 1749 tmpArrayFive[2] = pTJet - ptUeCone;
255cb71c 1750 tmpArrayFive[3] = triggerHadron->Pt();
1751 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
8aa19603 1752 fHJetSpecSubUeCone->Fill(tmpArrayFive);
255cb71c 1753
1754 //Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", kT median
1755 tmpArrayFive[0] = areaJet;
1756 tmpArrayFive[1] = pTJet;
8aa19603 1757 tmpArrayFive[2] = pTJet - ptUeFromCellMedian;
1758 tmpArrayFive[3] = ptUeFromCellMedian;
1759 tmpArrayFive[4] = rhoFromCellMedian;
255cb71c 1760 fHJetUeMedian->Fill(tmpArrayFive);
8e103a56 1761
8aa19603 1762 //Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", Cone median
255cb71c 1763 tmpArrayFive[0] = areaJet;
1764 tmpArrayFive[1] = pTJet;
8aa19603 1765 tmpArrayFive[2] = pTJet - ptUeCone;
1766 tmpArrayFive[3] = ptUeCone;
1767 tmpArrayFive[4] = rhoCone;
1768 fHJetUeCone->Fill(tmpArrayFive);
255cb71c 1769
1770 if(filledOnce){ //fill for each event only once
8aa19603 1771 Double_t fillRho[] = { rhoCone,rhoFromCellMedian};
1772 fHRhoUeMedianVsCone->Fill(fillRho);
255cb71c 1773 filledOnce = kFALSE;
1774 }
8e103a56 1775 }//jet loop
1776
1777 //Fill Jet Density In the Event A>0.07
255cb71c 1778 /*if(injet>0){
8e103a56 1779 Double_t filldens[]={ (Double_t) particleList.GetEntries(),
ad869500 1780 injet/fkAcceptance,
1781 triggerHadron->Pt()
1782 };
8e103a56 1783 fHJetDensity->Fill(filldens);
1784 }
ad869500 1785
8e103a56 1786 //Fill Jet Density In the Event A>0.4
1787 if(injet4>0){
1788 Double_t filldens4[]={ (Double_t) particleList.GetEntries(),
1789 injet4/fkAcceptance,
1790 triggerHadron->Pt()
1791 };
1792 fHJetDensityA4->Fill(filldens4);
255cb71c 1793 }*/
ad869500 1794 }
1795
8e103a56 1796
ad869500 1797 PostData(1, fOutputList);
1798}
1799
1800//----------------------------------------------------------------------------
1801void AliAnalysisTaskJetCorePP::Terminate(const Option_t *)
1802{
1803 // Draw result to the screen
1804 // Called once at the end of the query
1805
1806 if(fDebug) printf("AliAnalysisTaskJetCorePP DONE\n");
1807 if(!GetOutputData(1)) return;
1808}
1809
1810
1811//----------------------------------------------------------------------------
1812Int_t AliAnalysisTaskJetCorePP::GetListOfTracks(TList *list){
1813 //Fill the list of accepted tracks (passed track cut)
1814 //return consecutive index of the hardest ch hadron in the list
1815 Int_t iCount = 0;
1816 AliAODEvent *aodevt = NULL;
1817
1818 if(!fESD) aodevt = fAODIn;
1819 else aodevt = fAODOut;
1820
1821 if(!aodevt) return -1;
1822
1823 Int_t index = -1; //index of the highest particle in the list
1824 Double_t ptmax = -10;
8e103a56 1825 Int_t triggers[200];
1826 Int_t ntriggers = 0; //index in triggers array
1e3ad34d 1827
ad869500 1828
80ac66f6 1829 for(Int_t it = 0; it < aodevt->GetNumberOfTracks(); it++){
0fe881ea 1830 //AliAODTrack *tr = aodevt->GetTrack(it);
f15c1f69 1831 AliAODTrack *tr = dynamic_cast<AliAODTrack*>(aodevt->GetTrack(it));
1832 if(!tr) AliFatal("Not a standard AOD");
0fe881ea 1833
d1405a52
ML		 1834 if((fFilterMask > 0) && !(tr->TestFilterBit(fFilterMask))) continue;
1835 //if((fFilterMask > 0) && !(tr->IsHybridGlobalConstrainedGlobal())) continue;
ad869500 1836 if(TMath::Abs((Float_t) tr->Eta()) > fTrackEtaCut) continue;
1837 if(tr->Pt() < fTrackLowPtCut) continue;
1838 list->Add(tr);
8e103a56 1839
1840 //Search trigger:
1841 if(fHardest>0){ //leading particle
1842 if(tr->Pt()>ptmax){
1843 ptmax = tr->Pt();
1844 index = iCount;
1845 }
1846 }
1847
1848 if(fHardest==0 && ntriggers<200){ //single inclusive
1849 if(fTriggerPtRangeLow <= tr->Pt() &&
1850 tr->Pt() < fTriggerPtRangeHigh){
1851 triggers[ntriggers] = iCount;
1852 ntriggers++;
1853 }
ad869500 1854 }
8e103a56 1855
ad869500 1856 iCount++;
1857 }
1858
8e103a56 1859 if(fHardest==0 && ntriggers>0){ //select random inclusive trigger
1860 Int_t rnd = fRandom->Integer(ntriggers); //0 to ntriggers-1
1861 index = triggers[rnd];
1e3ad34d 1862
1863 fhNofMultipleTriggers->Fill(ntriggers-1);
ac4f24f7 1864
1865 Double_t deltaPhi, deltaEta, deltaR;
1866 Int_t k=0;
1867 //Correlation with single inclusive trigger
1868 AliVParticle* tGent1 = (AliVParticle*) list->At(index);
1869 if(tGent1){
1870 for(Int_t ia=0; ia<ntriggers; ia++){
1871 if(triggers[ia]==index) continue;
1872 AliVParticle* tGent2 = (AliVParticle*) list->At(triggers[ia]);
1873 if(!tGent2) continue;
1874 deltaPhi = RelativePhi(tGent1->Phi(),tGent2->Phi());
1875 deltaEta = tGent1->Eta()-tGent2->Eta();
1876 deltaR = sqrt(deltaPhi*deltaPhi + deltaEta*deltaEta);
1877 if(deltaR<0.4) k++;
1878 }
1879 }
1880 fhNofMultipleTriggersCone->Fill(k);
0fe881ea 1881
ac4f24f7 1882 if(ntriggers>1){
0fe881ea 1883 //Correlation with any other trigger
1e3ad34d 1884 for(Int_t ia=0; ia<ntriggers-1; ia++){
0fe881ea 1885 AliVParticle* tGeni = (AliVParticle*) list->At(triggers[ia]);
1e3ad34d 1886 if(!tGeni) continue;
1887 for(Int_t ib=ia+1; ib<ntriggers; ib++){
0fe881ea 1888 AliVParticle* tGenii = (AliVParticle*) list->At(triggers[ib]);
1e3ad34d 1889 if(!tGenii) continue;
1890 deltaPhi = RelativePhi(tGeni->Phi(),tGenii->Phi());
1891 deltaEta = tGeni->Eta()-tGenii->Eta();
1892 deltaR = sqrt(deltaPhi*deltaPhi + deltaEta*deltaEta);
1893 fhDeltaRMultTriggers->Fill(deltaR);
1894 }
1895 }
1896 }
ad869500 1897 }
8e103a56 1898
1899 return index;
ad869500 1900}
1901
1902//----------------------------------------------------------------------------
80ac66f6 1903/*
ad869500 1904Double_t AliAnalysisTaskJetCorePP::GetBackgroundInPerpCone(Float_t jetR, Double_t jetPhi, Double_t jetEta, TList* trkList){
1905 //calculate sum of track pT in the cone perpendicular in phi to the jet
1906 //jetR = cone radius
1907 // jetPhi, jetEta = direction of the jet
1908 Int_t numberOfTrks = trkList->GetEntries();
1909 Double_t pTsum = 0.0;
1910 Double_t perpConePhi = jetPhi + TMath::Pi()/2;//perp cone w.r.t. jet in phi
1911 for(Int_t it=0; it<numberOfTrks; it++){
1912 AliVParticle *trk = (AliVParticle*) trkList->At(it);
1913 Double_t dphi = RelativePhi(perpConePhi,trk->Phi());
1914 Double_t deta = trk->Eta()-jetEta;
1915
1916 if( (dphi*dphi + deta*deta)< (jetR*jetR)){
1917 pTsum += trk->Pt();
1918 }
1919 }
1920
1921 return pTsum;
1922}
80ac66f6 1923*/
ad869500 1924//----------------------------------------------------------------------------
1925
1926Double_t AliAnalysisTaskJetCorePP::RelativePhi(Double_t mphi,Double_t vphi){
1927 //Get relative azimuthal angle of two particles -pi to pi
1928 if (vphi < -TMath::Pi()) vphi += TMath::TwoPi();
1929 else if (vphi > TMath::Pi()) vphi -= TMath::TwoPi();
1930
1931 if (mphi < -TMath::Pi()) mphi += TMath::TwoPi();
1932 else if (mphi > TMath::Pi()) mphi -= TMath::TwoPi();
1933
1934 Double_t dphi = mphi - vphi;
1935 if (dphi < -TMath::Pi()) dphi += TMath::TwoPi();
1936 else if (dphi > TMath::Pi()) dphi -= TMath::TwoPi();
1937
1938 return dphi;//dphi in [-Pi, Pi]
1939}
1940
1941
80ac66f6 1942//----------------------------------------------------------------------------
1943Bool_t AliAnalysisTaskJetCorePP::SelectMCGenTracks(AliVParticle *trk, TList *trkList, Double_t &ptLeading, Int_t &index, Int_t counter){
1944 //fill track efficiency denominator
1945 if(trk->Pt() < fTrackLowPtCut) return kFALSE;
1946 if(trk->Charge()==0) return kFALSE;
1947 if(TMath::Abs((Float_t) trk->Eta()) > fTrackEtaCut) return kFALSE;
1948 trkList->Add(trk);
8aa19603 1949 if(fFillRespMx) fhPtTrkTruePrimGen->Fill(trk->Pt(),trk->Eta()); //Efficiency denominator
8e103a56 1950
1951 //Search trigger:
1952 if(fHardest>0){ //leading particle
1953 if(ptLeading < trk->Pt()){
1954 index = counter;
1955 ptLeading = trk->Pt();
1956 }
80ac66f6 1957 }
8e103a56 1958
1959 if(fHardest==0){ //single inclusive
1960 index = -1;
1961 if(fTriggerPtRangeLow <= trk->Pt() &&
1962 trk->Pt() < fTriggerPtRangeHigh){
1963 index = counter;
1964 }
1965 }
1966
80ac66f6 1967 return kTRUE;
1968}
1969
1970//----------------------------------------------------------------------------
1971void AliAnalysisTaskJetCorePP::FillEffHistos(TList *recList, TList *genList){
1972 //fill track efficiency numerator
1973 if(!recList) return;
1974 if(!genList) return;
1975 Int_t nRec = recList->GetEntries();
1976 Int_t nGen = genList->GetEntries();
1977 for(Int_t ir=0; ir<nRec; ir++){
1978 AliAODTrack *trkRec = (AliAODTrack*) recList->At(ir);
1979 if(!trkRec) continue;
1980 Int_t recLabel = TMath::Abs(trkRec->GetLabel());
1981 Bool_t matched = kFALSE;
1982
1983 for(Int_t ig=0; ig<nGen; ig++){
1984 AliVParticle *trkGen = (AliVParticle*) genList->At(ig);
1985 if(!trkGen) continue;
1986 Int_t genLabel = TMath::Abs(trkGen->GetLabel());
1987 if(recLabel==genLabel){
1988 fhPtTrkTruePrimRec->Fill(trkGen->Pt(),trkGen->Eta());
1989 matched = kTRUE;
1990 break;
1991 }
1992 }
1993 if(!matched) fhPtTrkSecOrFakeRec->Fill(trkRec->Pt(),trkRec->Eta());
1994 }
1995
1996 return;
1997}
255cb71c 1998//________________________________________________________________
d1405a52 1999void AliAnalysisTaskJetCorePP::EstimateBgRhoMedian(TList *listJet, TList* listPart, Double_t &rhoMedian, Int_t mode){
8aa19603 2000 //Estimate background rho by means of integrating track pT outside identified jet cones
2001
2002 rhoMedian = 0;
2003
2004 //identify leading jet in the full track acceptance
2005 Int_t idxLeadingJet = -1;
2006 Double_t pTleading = 0.0;
2007 AliAODJet* jet = NULL;
2008
2009 for(Int_t ij = 0; ij < listJet->GetEntries(); ij++){ //loop over bg kT jets
2010 jet = (AliAODJet*)(listJet->At(ij));
2011 if(!jet) continue;
2012 if(TMath::Abs(jet->Eta()) > fTrackEtaCut) continue; //track acceptance cut
2013 if(jet->Pt() > pTleading){
2014 idxLeadingJet = ij;
2015 pTleading = jet->Pt();
2016 }
2017 }
2018
2019 Int_t njets = 0;
2020 static Double_t jphi[100];
2021 static Double_t jeta[100];
2022 static Double_t jRsquared[100];
2023
2024 for(Int_t ij = 0; ij< listJet->GetEntries(); ij++){ //loop over bg kT jets
2025
2026 jet = (AliAODJet*)(listJet->At(ij));
2027 if(!jet) continue;
2028 if(TMath::Abs(jet->Eta()) > fTrackEtaCut) continue;
2029
2030 if((jet->Pt() > fBgMaxJetPt || ij== idxLeadingJet) && njets<100){
2031 jphi[njets] = RelativePhi(jet->Phi(),0.0); //-pi,pi
2032 jeta[njets] = jet->Eta();
2033 jRsquared[njets] = fSafetyMargin*jet->EffectiveAreaCharged()/TMath::Pi(); //scale jet radius
2034 njets++;
255cb71c 2035 }
2036 }
8aa19603 2037
255cb71c 2038
8aa19603 2039 static Double_t nOutCone[10][4];
2040 static Double_t sumPtOutOfCone[10][4];
2041
2042 for(Int_t ie=0; ie<fnEta; ie++){
2043 for(Int_t ip=0; ip<fnPhi; ip++){
2044 nOutCone[ip][ie] = 0.0; //initialize counter
2045 sumPtOutOfCone[ip][ie] = 0.0;
255cb71c 2046 }
8aa19603 2047 }
255cb71c 2048
8aa19603 2049 Double_t rndphi, rndeta;
2050 Double_t rndphishift, rndetashift;
2051 Double_t dphi, deta;
2052 Bool_t bIsInCone;
2053
2054 for(Int_t it=0; it<fnTrials; it++){
2055
2056 rndphi = fRandom->Uniform(0, fPhiSize);
2057 rndeta = fRandom->Uniform(0, fEtaSize);
2058
2059 for(Int_t ip=0; ip<fnPhi; ip++){ //move radom position to each cell
2060 rndphishift = rndphi + ip*fPhiSize - TMath::Pi();
2061 for(Int_t ie=0; ie<fnEta; ie++){
2062 rndetashift = rndeta + ie*fEtaSize - fEtaSize;
2063
2064 bIsInCone = 0; //tag if trial is in the jet cone
2065 for(Int_t ij=0; ij<njets; ij++){
2066 deta = jeta[ij] - rndetashift;
2067 dphi = RelativePhi(rndphishift,jphi[ij]);
2068 if((dphi*dphi + deta*deta) < jRsquared[ij]){
2069 bIsInCone = 1;
2070 break;
2071 }
2072 }
2073 if(!bIsInCone) nOutCone[ip][ie]++;
2074 }
2075 }
2076 }
2077
255cb71c 2078
8aa19603 2079 //Sum particle pT outside jets in cells
2080 Int_t npart = listPart->GetEntries();
2081 Int_t phicell,etacell;
2082 AliVParticle *part = NULL;
2083 for(Int_t ip=0; ip < npart; ip++){
2084
2085 part = (AliVParticle*) listPart->At(ip);
2086 if(!part) continue;
2087 if( TMath::Abs(part->Eta()) > fEtaSize) continue; //skip part outside +-0.7 in eta
2088
2089 bIsInCone = 0; //init
2090 for(Int_t ij=0; ij<njets; ij++){
2091 dphi = RelativePhi(jphi[ij], part->Phi());
2092 deta = jeta[ij] - part->Eta();
2093 if((dphi*dphi + deta*deta) < jRsquared[ij]){
2094 bIsInCone = 1;
2095 break;
2096 }
255cb71c 2097 }
8aa19603 2098 if(!bIsInCone){
2099 phicell = TMath::Nint(TMath::Floor((RelativePhi(part->Phi(),0.0) + TMath::Pi())/fPhiSize));
2100 etacell = TMath::Nint(TMath::Floor((part->Eta()+fEtaSize)/fEtaSize));
2101 sumPtOutOfCone[phicell][etacell]+= part->Pt();
2102 }
2103 }
2104
2105 static Double_t rhoInCells[20];
2106 Double_t relativeArea;
2107 Int_t nCells=0;
d1405a52 2108 Double_t bufferArea=0.0, bufferPt=0.0; //sum cells where A< fJetFreeAreaFrac
8aa19603 2109 for(Int_t ip=0; ip<fnPhi; ip++){
2110 for(Int_t ie=0; ie<fnEta; ie++){
2111 relativeArea = nOutCone[ip][ie]/fnTrials;
d1405a52
ML		 2112 //cout<<"RA= "<< relativeArea<<" BA="<<bufferArea<<" BPT="<< bufferPt <<endl;
2113
2114 bufferArea += relativeArea;
2115 bufferPt += sumPtOutOfCone[ip][ie];
2116 if(bufferArea > fJetFreeAreaFrac){
2117 rhoInCells[nCells] = bufferPt/(bufferArea*fCellArea);
2118
2119 if(mode==1) fhCellAreaToMedianGen->Fill(bufferArea*fCellArea);
2120 else fhCellAreaToMedian->Fill(bufferArea*fCellArea);
2121
2122 bufferArea = 0.0;
2123 bufferPt = 0.0;
2124 nCells++;
2125 }
8aa19603 2126 }
255cb71c 2127 }
8aa19603 2128
d1405a52
ML		 2129
2130 if(nCells>0){
2131 rhoMedian = TMath::Median(nCells, rhoInCells);
2132 if(mode==1){ //mc data
2133 fhEntriesToMedianGen->Fill(nCells);
2134 }else{ //real data
2135 fhEntriesToMedian->Fill(nCells);
2136 }
2137 }
8aa19603 2138
255cb71c 2139}
80ac66f6 2140
255cb71c 2141//__________________________________________________________________
8aa19603 2142void AliAnalysisTaskJetCorePP::EstimateBgCone(TList *listJet, TList* listPart, Double_t &rhoPerpCone){
2143
2144 rhoPerpCone = 0.0; //init
2145
2146 if(!listJet) return;
2147 Int_t njet = listJet->GetEntries();
2148 Int_t npart = listPart->GetEntries();
2149 Double_t pTleading = 0.0;
2150 Double_t phiLeading = 1000.;
2151 Double_t etaLeading = 1000.;
2152
2153 for(Int_t jsig=0; jsig < njet; jsig++){
2154 AliAODJet* signaljet = (AliAODJet*)(listJet->At(jsig));
2155 if(!signaljet) continue;
2156 if((signaljet->Eta()<fJetEtaMin) && (fJetEtaMax<signaljet->Eta())) continue; //acceptance cut
2157 if(signaljet->Pt() >= pTleading){ //replace leading and subleading jet
2158 pTleading = signaljet->Pt();
2159 phiLeading = signaljet->Phi();
2160 etaLeading = signaljet->Eta();
2161 }
2162 }
2163
2164 Double_t phileftcone = phiLeading + TMath::Pi()/2;
2165 Double_t phirightcone = phiLeading - TMath::Pi()/2;
2166 Double_t dp, de;
2167
2168 for(Int_t ip=0; ip < npart; ip++){
2169
2170 AliVParticle *part = (AliVParticle*) listPart->At(ip);
2171 if(!part){
2172 continue;
2173 }
2174
2175 dp = RelativePhi(phileftcone, part->Phi());
2176 de = etaLeading - part->Eta();
2177 if( sqrt(dp*dp + de*de)< fBgConeR) rhoPerpCone += part->Pt();
2178
2179 dp = RelativePhi(phirightcone, part->Phi());
2180 if( sqrt(dp*dp + de*de)< fBgConeR) rhoPerpCone += part->Pt();
2181
2182 }
2183
2184 //normalize total pT by two times cone are
2185 rhoPerpCone = rhoPerpCone/(2*TMath::Pi()*fBgConeR*fBgConeR);
2186
2187}
2188//__________________________________________________________________
80ac66f6 2189
255cb71c 2190void AliAnalysisTaskJetCorePP::ReadTClonesArray(TString bname, TList *list){
2191 //Convert TClones array of jets to TList
2192
2193 if(!strlen(bname.Data())){
2194 AliError(Form("Jet branch %s not set.", bname.Data()));
2195 return;
2196 }
2197
2198 TClonesArray *array=0x0;
49dc2855 2199 if(fUseExchContainer){ //take input from exchange containers
2200 if(fJetBranchNameKine.Length()>0 && bname.CompareTo(fJetBranchNameKine.Data())==0){
2201 array = dynamic_cast<TClonesArray*>(GetInputData(1)); //connect exchange container slot 1
2202 }
3d1fac8a 2203 if(fJetBranchNameBgKine.Length()>0 && bname.CompareTo(fJetBranchNameBgKine.Data())==0){
2204 array = dynamic_cast<TClonesArray*>(GetInputData(2)); //connect exchange container slot 2
2205 }
49dc2855 2206 }else{ //take input from AOD
2207 if(fAODOut&&!array){
2208 array = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(bname.Data()));
2209 }
2210 if(fAODExtension&&!array){
2211 array = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(bname.Data()));
2212 }
2213 if(fAODIn&&!array){
2214 array = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(bname.Data()));
2215 }
255cb71c 2216 }
2217
2218 list->Clear(); //clear list beforehand
2219
2220 if(array){
2221 if(fDebug){
2222 Printf("########## %s: %d jets",bname.Data(), array->GetEntriesFast());
2223 }
49dc2855 2224
255cb71c 2225 for(Int_t iJet = 0; iJet < array->GetEntriesFast(); iJet++) {
2226 AliAODJet *jet = dynamic_cast<AliAODJet*>((*array)[iJet]);
2227 if (jet) list->Add(jet);
2228 }
2229 }
2230
2231 return;
2232}