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