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