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