]> 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
New histograms for 2particle correl added
[u/mrichter/AliRoot.git] / PWGJE / AliAnalysisTaskJetCorePP.cxx
CommitLineData
ad869500 1
2// ******************************************
3// This task computes several jet observables like
4// the fraction of energy in inner and outer coronnas,
5// jet-track correlations,triggered jet shapes and
6// correlation strength distribution of particles inside jets.
7// Author: lcunquei@cern.ch
8// *******************************************
9
10
11/**************************************************************************
12 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
13 * *
14 * Author: The ALICE Off-line Project. *
15 * Contributors are mentioned in the code where appropriate. *
16 * *
17 * Permission to use, copy, modify and distribute this software and its *
18 * documentation strictly for non-commercial purposes is hereby granted *
19 * without fee, provided that the above copyright notice appears in all *
20 * copies and that both the copyright notice and this permission notice *
21 * appear in the supporting documentation. The authors make no claims *
22 * about the suitability of this software for any purpose. It is *
23 * provided "as is" without express or implied warranty. *
24 **************************************************************************/
25
26
27#include "TChain.h"
49dc2855 28#include "TClonesArray.h"
ad869500 29#include "TTree.h"
80ac66f6 30#include "TList.h"
ad869500 31#include "TMath.h"
32#include "TH1F.h"
33#include "TH1D.h"
34#include "TH2F.h"
35#include "TH3F.h"
36#include "THnSparse.h"
37#include "TCanvas.h"
80ac66f6 38#include "TArrayI.h"
ea603b64 39#include "TProfile.h"
40#include "TFile.h"
41#include "TKey.h"
8e103a56 42#include "TRandom3.h"
ad869500 43
44#include "AliLog.h"
45
46#include "AliAnalysisTask.h"
47#include "AliAnalysisManager.h"
48
49#include "AliVEvent.h"
50#include "AliESDEvent.h"
80ac66f6 51#include "AliMCEvent.h"
ad869500 52#include "AliESDInputHandler.h"
53#include "AliCentrality.h"
54#include "AliAnalysisHelperJetTasks.h"
55#include "AliInputEventHandler.h"
56#include "AliAODJetEventBackground.h"
ea603b64 57#include "AliGenPythiaEventHeader.h"
ad869500 58#include "AliAODMCParticle.h"
80ac66f6 59#include "AliMCParticle.h"
ad869500 60#include "AliAODEvent.h"
61#include "AliAODHandler.h"
62#include "AliAODJet.h"
49dc2855 63#include "AliVVertex.h"
0fe881ea 64#include "AliVTrack.h"
ad869500 65#include "AliAnalysisTaskJetCorePP.h"
29ac2c9b 66#include "AliHeader.h" //KF//
ad869500 67
68using std::cout;
69using std::endl;
70
71ClassImp(AliAnalysisTaskJetCorePP)
72
73//Filip Krizek 1st March 2013
74
75//---------------------------------------------------------------------
76AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP() :
77AliAnalysisTaskSE(),
78fESD(0x0),
79fAODIn(0x0),
80fAODOut(0x0),
81fAODExtension(0x0),
49dc2855 82fMcEvent(0x0),
83fMcHandler(0x0),
ad869500 84fJetBranchName(""),
255cb71c 85fJetBranchNameChargMC(""),
49dc2855 86fJetBranchNameKine(""),
255cb71c 87fJetBranchNameFullMC(""),
88fJetBranchNameBg(""),
89fJetBranchNameBgChargMC(""),
49dc2855 90fJetBranchNameBgKine(""),
7fc3d134 91fListJets(0x0),
80ac66f6 92fListJetsGen(0x0),
160dbdb8 93fListJetsGenFull(0x0),
255cb71c 94fListJetsBg(0x0),
95fListJetsBgGen(0x0),
ad869500 96fNonStdFile(""),
97fSystem(0), //pp=0 pPb=1
98fJetParamR(0.4),
8aa19603 99fBgJetParamR(0.3),
100fBgMaxJetPt(8.0),
101fBgConeR(0.4),
ad869500 102fOfflineTrgMask(AliVEvent::kAny),
103fMinContribVtx(1),
104fVtxZMin(-10.0),
105fVtxZMax(10.0),
106fFilterMask(0),
107fCentMin(0.0),
108fCentMax(100.0),
109fJetEtaMin(-0.5),
110fJetEtaMax(0.5),
111fTriggerEtaCut(0.9),
112fTrackEtaCut(0.9),
113fTrackLowPtCut(0.15),
49dc2855 114fUseExchContainer(0),
ad869500 115fOutputList(0x0),
116fHistEvtSelection(0x0),
117fh2Ntriggers(0x0),
118fHJetSpec(0x0),
255cb71c 119fHJetSpecSubUeMedian(0x0),
8aa19603 120fHJetSpecSubUeCone(0x0),
49dc2855 121fHJetPhiCorr(0x0),
122fHJetPhiCorrSubUeMedian(0x0),
123fHJetPhiCorrSubUeCone(0x0),
255cb71c 124fHJetUeMedian(0x0),
8aa19603 125fHJetUeCone(0x0),
126fHRhoUeMedianVsCone(0x0),
255cb71c 127//fHJetDensity(0x0),
128//fHJetDensityA4(0x0),
ad869500 129fhJetPhi(0x0),
130fhTriggerPhi(0x0),
131fhJetEta(0x0),
132fhTriggerEta(0x0),
133fhVertexZ(0x0),
134fhVertexZAccept(0x0),
135fhContribVtx(0x0),
136fhContribVtxAccept(0x0),
137fhDphiTriggerJet(0x0),
138fhDphiTriggerJetAccept(0x0),
139fhCentrality(0x0),
140fhCentralityAccept(0x0),
1e3ad34d 141fhNofMultipleTriggers(0x0),
0fe881ea 142fhNofMultipleTriggersCone(0x0),
30557bf5 143fhNofMultipleTriggersConeLow(0x0),
144fhNofMultipleTriggersConeHigh(0x0),
364e98aa 145fhDeltaRMultTriggersLow(0x0),
146fhDeltaRMultTriggersHigh(0x0),
f1e6baa6 147fhDeltaPhiMultTriggersLow(0x0),
148fhDeltaPhiMultTriggersHigh(0x0),
149fhDeltaPhiMultTriggersInclLow(0x0),
150fhDeltaPhiMultTriggersInclHigh(0x0),
151fhInclTrigCounter(0x0),
255cb71c 152//fHJetPtRaw(0x0),
153//fHLeadingJetPtRaw(0x0),
154//fHDphiVsJetPtAll(0x0),
80ac66f6 155fhJetPtGenVsJetPtRec(0x0),
255cb71c 156fhJetPtGenVsJetPtRecSubUeMedian(0x0),
8aa19603 157fhJetPtGenVsJetPtRecSubUeCone(0x0),
255cb71c 158fhJetPtGen(0x0),
159fhJetPtSubUeMedianGen(0x0),
f1e6baa6 160fhJetPtSubUeConeGen(0x0),
161fhJetPtResolutionVsPtGen(0x0),
162fhJetPtResolutionVsPtConeGen(0x0),
160dbdb8 163fhJetPtGenChargVsJetPtGenFull(0x0),
164fhJetPtGenFull(0x0),
80ac66f6 165fh2NtriggersGen(0x0),
166fHJetSpecGen(0x0),
255cb71c 167fHJetSpecSubUeMedianGen(0x0),
8aa19603 168fHJetSpecSubUeConeGen(0x0),
49dc2855 169fHJetPhiCorrGen(0x0),
170fHJetPhiCorrSubUeMedianGen(0x0),
171fHJetPhiCorrSubUeConeGen(0x0),
255cb71c 172fHJetUeMedianGen(0x0),
8aa19603 173fHJetUeConeGen(0x0),
80ac66f6 174fhPtTrkTruePrimRec(0x0),
175fhPtTrkTruePrimGen(0x0),
176fhPtTrkSecOrFakeRec(0x0),
8aa19603 177fHRhoUeMedianVsConeGen(0x0),
d1405a52
ML		 178fhEntriesToMedian(0x0),
179fhEntriesToMedianGen(0x0),
180fhCellAreaToMedian(0x0),
181fhCellAreaToMedianGen(0x0),
1e3ad34d 182fhNofMultipleTriggersGen(0x0),
0fe881ea 183fhNofMultipleTriggersConeGen(0x0),
30557bf5 184fhNofMultipleTriggersConeGenLow(0x0),
185fhNofMultipleTriggersConeGenHigh(0x0),
364e98aa 186fhDeltaRMultTriggersGenLow(0x0),
187fhDeltaRMultTriggersGenHigh(0x0),
f1e6baa6 188fhDeltaPhiMultTriggersGenLow(0x0),
189fhDeltaPhiMultTriggersGenHigh(0x0),
190fhDeltaPhiMultTriggersInclGenLow(0x0),
191fhDeltaPhiMultTriggersInclGenHigh(0x0),
192fhInclTrigCounterGen(0x0),
6c1489e1 193fhNofMultipleTriggersConeGenA(0x0),
30557bf5 194fhNofMultipleTriggersConeGenALow(0x0),
195fhNofMultipleTriggersConeGenAHigh(0x0),
364e98aa 196fhDeltaRMultTriggersGenALow(0x0),
197fhDeltaPhiMultTriggersGenALow(0x0),
198fhDeltaRMultTriggersGenAHigh(0x0),
199fhDeltaPhiMultTriggersGenAHigh(0x0),
200fhNofTriggersGenA(0x0),
6c1489e1 201fhNofMultipleTriggersConeGenB(0x0),
30557bf5 202fhNofMultipleTriggersConeGenBLow(0x0),
203fhNofMultipleTriggersConeGenBHigh(0x0),
364e98aa 204fhDeltaRMultTriggersGenBLow(0x0),
205fhDeltaPhiMultTriggersGenBLow(0x0),
206fhDeltaRMultTriggersGenBHigh(0x0),
207fhDeltaPhiMultTriggersGenBHigh(0x0),
208fhNofTriggersGenB(0x0),
209fhTriggerCounterGenLevel(0x0),
210fhDeltaRMultTriggersGenLevelLow(0x0),
211fhDeltaPhiMultTriggersGenLevelLow(0x0),
212fhDeltaRMultTriggersGenLevelHigh(0x0),
213fhDeltaPhiMultTriggersGenLevelHigh(0x0),
255cb71c 214fIsChargedMC(0),
49dc2855 215fIsKine(0),
255cb71c 216fIsFullMC(0),
80ac66f6 217faGenIndex(0),
218faRecIndex(0),
7fc3d134 219fkAcceptance(2.0*TMath::Pi()*1.8),
6168afc3 220fkDeltaPhiCut(TMath::Pi()-0.8),
ea603b64 221fh1Xsec(0x0),
222fh1Trials(0x0),
223fh1AvgTrials(0x0),
224fh1PtHard(0x0),
225fh1PtHardNoW(0x0),
226fh1PtHardTrials(0x0),
8e103a56 227fAvgTrials(1),
228fHardest(0),
229fEventNumberRangeLow(0),
230fEventNumberRangeHigh(99),
231fTriggerPtRangeLow(0.0),
232fTriggerPtRangeHigh(10000.0),
8aa19603 233fFillRespMx(0),
234fRandom(0x0),
235fnTrials(1000),
d1405a52 236fJetFreeAreaFrac(0.5),
8aa19603 237fnPhi(9),
238fnEta(2),
239fEtaSize(0.7),
240fPhiSize(2*TMath::Pi()/fnPhi),
241fCellArea(fPhiSize*fEtaSize),
49dc2855 242fSafetyMargin(1.1),
243fDoubleBinning(kFALSE)
ad869500 244{
245 // default Constructor
ad869500 246}
247
248//---------------------------------------------------------------------
249
250AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP(const char *name) :
251AliAnalysisTaskSE(name),
252fESD(0x0),
253fAODIn(0x0),
254fAODOut(0x0),
255fAODExtension(0x0),
49dc2855 256fMcEvent(0x0),
257fMcHandler(0x0),
ad869500 258fJetBranchName(""),
255cb71c 259fJetBranchNameChargMC(""),
49dc2855 260fJetBranchNameKine(""),
255cb71c 261fJetBranchNameFullMC(""),
262fJetBranchNameBg(""),
263fJetBranchNameBgChargMC(""),
49dc2855 264fJetBranchNameBgKine(""),
7fc3d134 265fListJets(0x0),
80ac66f6 266fListJetsGen(0x0),
160dbdb8 267fListJetsGenFull(0x0),
255cb71c 268fListJetsBg(0x0),
269fListJetsBgGen(0x0),
ad869500 270fNonStdFile(""),
271fSystem(0), //pp=0 pPb=1
272fJetParamR(0.4),
8aa19603 273fBgJetParamR(0.3),
274fBgMaxJetPt(8.0),
275fBgConeR(0.4),
ad869500 276fOfflineTrgMask(AliVEvent::kAny),
277fMinContribVtx(1),
278fVtxZMin(-10.0),
279fVtxZMax(10.0),
280fFilterMask(0),
281fCentMin(0.0),
282fCentMax(100.0),
283fJetEtaMin(-0.5),
284fJetEtaMax(0.5),
285fTriggerEtaCut(0.9),
286fTrackEtaCut(0.9),
287fTrackLowPtCut(0.15),
49dc2855 288fUseExchContainer(0),
ad869500 289fOutputList(0x0),
290fHistEvtSelection(0x0),
291fh2Ntriggers(0x0),
292fHJetSpec(0x0),
255cb71c 293fHJetSpecSubUeMedian(0x0),
8aa19603 294fHJetSpecSubUeCone(0x0),
49dc2855 295fHJetPhiCorr(0x0),
296fHJetPhiCorrSubUeMedian(0x0),
297fHJetPhiCorrSubUeCone(0x0),
255cb71c 298fHJetUeMedian(0x0),
8aa19603 299fHJetUeCone(0x0),
300fHRhoUeMedianVsCone(0x0),
255cb71c 301//fHJetDensity(0x0),
302//fHJetDensityA4(0x0),
ad869500 303fhJetPhi(0x0),
304fhTriggerPhi(0x0),
305fhJetEta(0x0),
306fhTriggerEta(0x0),
307fhVertexZ(0x0),
308fhVertexZAccept(0x0),
309fhContribVtx(0x0),
310fhContribVtxAccept(0x0),
311fhDphiTriggerJet(0x0),
312fhDphiTriggerJetAccept(0x0),
313fhCentrality(0x0),
314fhCentralityAccept(0x0),
1e3ad34d 315fhNofMultipleTriggers(0x0),
0fe881ea 316fhNofMultipleTriggersCone(0x0),
30557bf5 317fhNofMultipleTriggersConeLow(0x0),
318fhNofMultipleTriggersConeHigh(0x0),
364e98aa 319fhDeltaRMultTriggersLow(0x0),
320fhDeltaRMultTriggersHigh(0x0),
f1e6baa6 321fhDeltaPhiMultTriggersLow(0x0),
322fhDeltaPhiMultTriggersHigh(0x0),
323fhDeltaPhiMultTriggersInclLow(0x0),
324fhDeltaPhiMultTriggersInclHigh(0x0),
325fhInclTrigCounter(0x0),
255cb71c 326//fHJetPtRaw(0x0),
327//fHLeadingJetPtRaw(0x0),
328//fHDphiVsJetPtAll(0x0),
80ac66f6 329fhJetPtGenVsJetPtRec(0x0),
255cb71c 330fhJetPtGenVsJetPtRecSubUeMedian(0x0),
8aa19603 331fhJetPtGenVsJetPtRecSubUeCone(0x0),
80ac66f6 332fhJetPtGen(0x0),
255cb71c 333fhJetPtSubUeMedianGen(0x0),
8aa19603 334fhJetPtSubUeConeGen(0x0),
f1e6baa6 335fhJetPtResolutionVsPtGen(0x0),
336fhJetPtResolutionVsPtConeGen(0x0),
160dbdb8 337fhJetPtGenChargVsJetPtGenFull(0x0),
338fhJetPtGenFull(0x0),
80ac66f6 339fh2NtriggersGen(0x0),
340fHJetSpecGen(0x0),
255cb71c 341fHJetSpecSubUeMedianGen(0x0),
8aa19603 342fHJetSpecSubUeConeGen(0x0),
49dc2855 343fHJetPhiCorrGen(0x0),
344fHJetPhiCorrSubUeMedianGen(0x0),
345fHJetPhiCorrSubUeConeGen(0x0),
255cb71c 346fHJetUeMedianGen(0x0),
8aa19603 347fHJetUeConeGen(0x0),
80ac66f6 348fhPtTrkTruePrimRec(0x0),
349fhPtTrkTruePrimGen(0x0),
350fhPtTrkSecOrFakeRec(0x0),
8aa19603 351fHRhoUeMedianVsConeGen(0x0),
d1405a52
ML		 352fhEntriesToMedian(0x0),
353fhEntriesToMedianGen(0x0),
354fhCellAreaToMedian(0x0),
355fhCellAreaToMedianGen(0x0),
1e3ad34d 356fhNofMultipleTriggersGen(0x0),
0fe881ea 357fhNofMultipleTriggersConeGen(0x0),
30557bf5 358fhNofMultipleTriggersConeGenLow(0x0),
359fhNofMultipleTriggersConeGenHigh(0x0),
364e98aa 360fhDeltaRMultTriggersGenLow(0x0),
361fhDeltaRMultTriggersGenHigh(0x0),
f1e6baa6 362fhDeltaPhiMultTriggersGenLow(0x0),
363fhDeltaPhiMultTriggersGenHigh(0x0),
364fhDeltaPhiMultTriggersInclGenLow(0x0),
365fhDeltaPhiMultTriggersInclGenHigh(0x0),
366fhInclTrigCounterGen(0x0),
6c1489e1 367fhNofMultipleTriggersConeGenA(0x0),
30557bf5 368fhNofMultipleTriggersConeGenALow(0x0),
369fhNofMultipleTriggersConeGenAHigh(0x0),
364e98aa 370fhDeltaRMultTriggersGenALow(0x0),
371fhDeltaPhiMultTriggersGenALow(0x0),
372fhDeltaRMultTriggersGenAHigh(0x0),
373fhDeltaPhiMultTriggersGenAHigh(0x0),
374fhNofTriggersGenA(0x0),
6c1489e1 375fhNofMultipleTriggersConeGenB(0x0),
30557bf5 376fhNofMultipleTriggersConeGenBLow(0x0),
377fhNofMultipleTriggersConeGenBHigh(0x0),
364e98aa 378fhDeltaRMultTriggersGenBLow(0x0),
379fhDeltaPhiMultTriggersGenBLow(0x0),
380fhDeltaRMultTriggersGenBHigh(0x0),
381fhDeltaPhiMultTriggersGenBHigh(0x0),
382fhNofTriggersGenB(0x0),
383fhTriggerCounterGenLevel(0x0),
384fhDeltaRMultTriggersGenLevelLow(0x0),
385fhDeltaPhiMultTriggersGenLevelLow(0x0),
386fhDeltaRMultTriggersGenLevelHigh(0x0),
387fhDeltaPhiMultTriggersGenLevelHigh(0x0),
255cb71c 388fIsChargedMC(0),
49dc2855 389fIsKine(0),
255cb71c 390fIsFullMC(0),
80ac66f6 391faGenIndex(0),
392faRecIndex(0),
7fc3d134 393fkAcceptance(2.0*TMath::Pi()*1.8),
6168afc3 394fkDeltaPhiCut(TMath::Pi()-0.8),
ea603b64 395fh1Xsec(0x0),
396fh1Trials(0x0),
397fh1AvgTrials(0x0),
398fh1PtHard(0x0),
399fh1PtHardNoW(0x0),
400fh1PtHardTrials(0x0),
8e103a56 401fAvgTrials(1),
402fHardest(0),
403fEventNumberRangeLow(0),
404fEventNumberRangeHigh(99),
405fTriggerPtRangeLow(0.0),
406fTriggerPtRangeHigh(10000.0),
8aa19603 407fFillRespMx(0),
408fRandom(0x0),
409fnTrials(1000),
d1405a52 410fJetFreeAreaFrac(0.5),
8aa19603 411fnPhi(9),
412fnEta(2),
413fEtaSize(0.7),
414fPhiSize(2*TMath::Pi()/fnPhi),
415fCellArea(fPhiSize*fEtaSize),
49dc2855 416fSafetyMargin(1.1),
417fDoubleBinning(kFALSE)
ad869500 418{
7fc3d134 419// Constructor
ad869500 420 DefineOutput(1, TList::Class());
49dc2855 421
422 TString dummy(name);
423 if(dummy.Contains("KINE")){
424 DefineInput(1, TClonesArray::Class());
3d1fac8a 425 DefineInput(2, TClonesArray::Class());
49dc2855 426 }
ad869500 427}
428
429//--------------------------------------------------------------
430AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP(const AliAnalysisTaskJetCorePP& a):
431AliAnalysisTaskSE(a.GetName()),
432fESD(a.fESD),
433fAODIn(a.fAODIn),
434fAODOut(a.fAODOut),
435fAODExtension(a.fAODExtension),
49dc2855 436fMcEvent(a.fMcEvent),
437fMcHandler(a.fMcHandler),
ad869500 438fJetBranchName(a.fJetBranchName),
255cb71c 439fJetBranchNameChargMC(a.fJetBranchNameChargMC),
49dc2855 440fJetBranchNameKine(a.fJetBranchNameKine),
255cb71c 441fJetBranchNameFullMC(a.fJetBranchNameFullMC),
442fJetBranchNameBg(a.fJetBranchNameBg),
443fJetBranchNameBgChargMC(a.fJetBranchNameBgChargMC),
49dc2855 444fJetBranchNameBgKine(a.fJetBranchNameBgKine),
ad869500 445fListJets(a.fListJets),
80ac66f6 446fListJetsGen(a.fListJetsGen),
160dbdb8 447fListJetsGenFull(a.fListJetsGenFull),
255cb71c 448fListJetsBg(a.fListJetsBg),
449fListJetsBgGen(a.fListJetsBgGen),
ad869500 450fNonStdFile(a.fNonStdFile),
451fSystem(a.fSystem),
452fJetParamR(a.fJetParamR),
8aa19603 453fBgJetParamR(a.fBgJetParamR),
454fBgMaxJetPt(a.fBgMaxJetPt),
455fBgConeR(a.fBgConeR),
ad869500 456fOfflineTrgMask(a.fOfflineTrgMask),
457fMinContribVtx(a.fMinContribVtx),
458fVtxZMin(a.fVtxZMin),
459fVtxZMax(a.fVtxZMax),
460fFilterMask(a.fFilterMask),
461fCentMin(a.fCentMin),
462fCentMax(a.fCentMax),
463fJetEtaMin(a.fJetEtaMin),
464fJetEtaMax(a.fJetEtaMax),
465fTriggerEtaCut(a.fTriggerEtaCut),
466fTrackEtaCut(a.fTrackEtaCut),
467fTrackLowPtCut(a.fTrackLowPtCut),
49dc2855 468fUseExchContainer(a.fUseExchContainer),
ad869500 469fOutputList(a.fOutputList),
470fHistEvtSelection(a.fHistEvtSelection),
471fh2Ntriggers(a.fh2Ntriggers),
472fHJetSpec(a.fHJetSpec),
255cb71c 473fHJetSpecSubUeMedian(a.fHJetSpecSubUeMedian),
8aa19603 474fHJetSpecSubUeCone(a.fHJetSpecSubUeCone),
49dc2855 475fHJetPhiCorr(a.fHJetPhiCorr),
476fHJetPhiCorrSubUeMedian(a.fHJetPhiCorrSubUeMedian),
477fHJetPhiCorrSubUeCone(a.fHJetPhiCorrSubUeCone),
255cb71c 478fHJetUeMedian(a.fHJetUeMedian),
8aa19603 479fHJetUeCone(a.fHJetUeCone),
480fHRhoUeMedianVsCone(a.fHRhoUeMedianVsCone),
255cb71c 481//fHJetDensity(a.fHJetDensity),
482//fHJetDensityA4(a.fHJetDensityA4),
ad869500 483fhJetPhi(a.fhJetPhi),
484fhTriggerPhi(a.fhTriggerPhi),
485fhJetEta(a.fhJetEta),
486fhTriggerEta(a.fhTriggerEta),
487fhVertexZ(a.fhVertexZ),
488fhVertexZAccept(a.fhVertexZAccept),
489fhContribVtx(a.fhContribVtx),
490fhContribVtxAccept(a.fhContribVtxAccept),
491fhDphiTriggerJet(a.fhDphiTriggerJet),
492fhDphiTriggerJetAccept(a.fhDphiTriggerJetAccept),
493fhCentrality(a.fhCentrality),
494fhCentralityAccept(a.fhCentralityAccept),
1e3ad34d 495fhNofMultipleTriggers(a.fhNofMultipleTriggers),
0fe881ea 496fhNofMultipleTriggersCone(a.fhNofMultipleTriggersCone),
30557bf5 497fhNofMultipleTriggersConeLow(a.fhNofMultipleTriggersConeLow),
498fhNofMultipleTriggersConeHigh(a.fhNofMultipleTriggersConeHigh),
364e98aa 499fhDeltaRMultTriggersLow(a.fhDeltaRMultTriggersLow),
500fhDeltaRMultTriggersHigh(a.fhDeltaRMultTriggersHigh),
f1e6baa6 501fhDeltaPhiMultTriggersLow(a.fhDeltaPhiMultTriggersLow),
502fhDeltaPhiMultTriggersHigh(a.fhDeltaPhiMultTriggersHigh),
503fhDeltaPhiMultTriggersInclLow(a.fhDeltaPhiMultTriggersInclLow),
504fhDeltaPhiMultTriggersInclHigh(a.fhDeltaPhiMultTriggersInclHigh),
505fhInclTrigCounter(a.fhInclTrigCounter),
255cb71c 506//fHJetPtRaw(a.fHJetPtRaw),
507//fHLeadingJetPtRaw(a.fHLeadingJetPtRaw),
508//fHDphiVsJetPtAll(a.fHDphiVsJetPtAll),
80ac66f6 509fhJetPtGenVsJetPtRec(a.fhJetPtGenVsJetPtRec),
255cb71c 510fhJetPtGenVsJetPtRecSubUeMedian(a.fhJetPtGenVsJetPtRecSubUeMedian),
8aa19603 511fhJetPtGenVsJetPtRecSubUeCone(a.fhJetPtGenVsJetPtRecSubUeCone),
80ac66f6 512fhJetPtGen(a.fhJetPtGen),
255cb71c 513fhJetPtSubUeMedianGen(a.fhJetPtSubUeMedianGen),
8aa19603 514fhJetPtSubUeConeGen(a.fhJetPtSubUeConeGen),
f1e6baa6 515fhJetPtResolutionVsPtGen(a.fhJetPtResolutionVsPtGen),
516fhJetPtResolutionVsPtConeGen(a.fhJetPtResolutionVsPtConeGen),
160dbdb8 517fhJetPtGenChargVsJetPtGenFull(a.fhJetPtGenChargVsJetPtGenFull),
518fhJetPtGenFull(a.fhJetPtGenFull),
80ac66f6 519fh2NtriggersGen(a.fh2NtriggersGen),
520fHJetSpecGen(a.fHJetSpecGen),
255cb71c 521fHJetSpecSubUeMedianGen(a.fHJetSpecSubUeMedianGen),
8aa19603 522fHJetSpecSubUeConeGen(a.fHJetSpecSubUeConeGen),
49dc2855 523fHJetPhiCorrGen(a.fHJetPhiCorrGen),
524fHJetPhiCorrSubUeMedianGen(a.fHJetPhiCorrSubUeMedianGen),
525fHJetPhiCorrSubUeConeGen(a.fHJetPhiCorrSubUeConeGen),
255cb71c 526fHJetUeMedianGen(a.fHJetUeMedianGen),
8aa19603 527fHJetUeConeGen(a.fHJetUeConeGen),
80ac66f6 528fhPtTrkTruePrimRec(a.fhPtTrkTruePrimRec),
529fhPtTrkTruePrimGen(a.fhPtTrkTruePrimGen),
530fhPtTrkSecOrFakeRec(a.fhPtTrkSecOrFakeRec),
8aa19603 531fHRhoUeMedianVsConeGen(a.fHRhoUeMedianVsConeGen),
d1405a52
ML		 532fhEntriesToMedian(a.fhEntriesToMedian),
533fhEntriesToMedianGen(a.fhEntriesToMedianGen),
534fhCellAreaToMedian(a.fhCellAreaToMedian),
535fhCellAreaToMedianGen(a.fhCellAreaToMedianGen),
1e3ad34d 536fhNofMultipleTriggersGen(a.fhNofMultipleTriggersGen),
0fe881ea 537fhNofMultipleTriggersConeGen(a.fhNofMultipleTriggersConeGen),
30557bf5 538fhNofMultipleTriggersConeGenLow(a.fhNofMultipleTriggersConeGenLow),
539fhNofMultipleTriggersConeGenHigh(a.fhNofMultipleTriggersConeGenHigh),
364e98aa 540fhDeltaRMultTriggersGenLow(a.fhDeltaRMultTriggersGenLow),
541fhDeltaRMultTriggersGenHigh(a.fhDeltaRMultTriggersGenHigh),
f1e6baa6 542fhDeltaPhiMultTriggersGenLow(a.fhDeltaPhiMultTriggersGenLow),
543fhDeltaPhiMultTriggersGenHigh(a.fhDeltaPhiMultTriggersGenHigh),
544fhDeltaPhiMultTriggersInclGenLow(a.fhDeltaPhiMultTriggersInclGenLow),
545fhDeltaPhiMultTriggersInclGenHigh(a.fhDeltaPhiMultTriggersInclGenHigh),
546fhInclTrigCounterGen(a.fhInclTrigCounterGen),
6c1489e1 547fhNofMultipleTriggersConeGenA(a.fhNofMultipleTriggersConeGenA),
30557bf5 548fhNofMultipleTriggersConeGenALow(a.fhNofMultipleTriggersConeGenALow),
549fhNofMultipleTriggersConeGenAHigh(a.fhNofMultipleTriggersConeGenAHigh),
364e98aa 550fhDeltaRMultTriggersGenALow(a.fhDeltaRMultTriggersGenALow),
551fhDeltaPhiMultTriggersGenALow(a.fhDeltaPhiMultTriggersGenALow),
552fhDeltaRMultTriggersGenAHigh(a.fhDeltaRMultTriggersGenAHigh),
553fhDeltaPhiMultTriggersGenAHigh(a.fhDeltaPhiMultTriggersGenAHigh),
554fhNofTriggersGenA(a.fhNofTriggersGenA),
6c1489e1 555fhNofMultipleTriggersConeGenB(a.fhNofMultipleTriggersConeGenB),
30557bf5 556fhNofMultipleTriggersConeGenBLow(a.fhNofMultipleTriggersConeGenBLow),
557fhNofMultipleTriggersConeGenBHigh(a.fhNofMultipleTriggersConeGenBHigh),
364e98aa 558fhDeltaRMultTriggersGenBLow(a.fhDeltaRMultTriggersGenBLow),
559fhDeltaPhiMultTriggersGenBLow(a.fhDeltaPhiMultTriggersGenBLow),
560fhDeltaRMultTriggersGenBHigh(a.fhDeltaRMultTriggersGenBHigh),
561fhDeltaPhiMultTriggersGenBHigh(a.fhDeltaPhiMultTriggersGenBHigh),
562fhNofTriggersGenB(a.fhNofTriggersGenB),
563fhTriggerCounterGenLevel(a.fhTriggerCounterGenLevel),
564fhDeltaRMultTriggersGenLevelLow(a.fhDeltaRMultTriggersGenLevelLow),
565fhDeltaPhiMultTriggersGenLevelLow(a.fhDeltaPhiMultTriggersGenLevelLow),
566fhDeltaRMultTriggersGenLevelHigh(a.fhDeltaRMultTriggersGenLevelHigh),
567fhDeltaPhiMultTriggersGenLevelHigh(a.fhDeltaPhiMultTriggersGenLevelHigh),
255cb71c 568fIsChargedMC(a.fIsChargedMC),
49dc2855 569fIsKine(a.fIsKine),
255cb71c 570fIsFullMC(a.fIsFullMC),
80ac66f6 571faGenIndex(a.faGenIndex),
572faRecIndex(a.faRecIndex),
7fc3d134 573fkAcceptance(a.fkAcceptance),
ea603b64 574fkDeltaPhiCut(a.fkDeltaPhiCut),
575fh1Xsec(a.fh1Xsec),
576fh1Trials(a.fh1Trials),
577fh1AvgTrials(a.fh1AvgTrials),
578fh1PtHard(a.fh1PtHard),
579fh1PtHardNoW(a.fh1PtHardNoW),
580fh1PtHardTrials(a.fh1PtHardTrials),
8e103a56 581fAvgTrials(a.fAvgTrials),
582fHardest(a.fHardest),
583fEventNumberRangeLow(a.fEventNumberRangeLow),
584fEventNumberRangeHigh(a.fEventNumberRangeHigh),
585fTriggerPtRangeLow(a.fTriggerPtRangeLow),
43457761 586fTriggerPtRangeHigh(a.fTriggerPtRangeHigh),
8aa19603 587fFillRespMx(a.fFillRespMx),
588fRandom(a.fRandom),
589fnTrials(a.fnTrials),
d1405a52 590fJetFreeAreaFrac(a.fJetFreeAreaFrac),
8aa19603 591fnPhi(a.fnPhi),
592fnEta(a.fnEta),
593fEtaSize(a.fEtaSize),
594fPhiSize(a.fPhiSize),
595fCellArea(a.fCellArea),
49dc2855 596fSafetyMargin(a.fSafetyMargin),
597fDoubleBinning(a.fDoubleBinning)
ad869500 598{
599 //Copy Constructor
43457761 600 fRandom->SetSeed(0);
ad869500 601}
602//--------------------------------------------------------------
603
604AliAnalysisTaskJetCorePP& AliAnalysisTaskJetCorePP::operator = (const AliAnalysisTaskJetCorePP& a){
605 // assignment operator
606 this->~AliAnalysisTaskJetCorePP();
607 new(this) AliAnalysisTaskJetCorePP(a);
608 return *this;
609}
610//--------------------------------------------------------------
611
612AliAnalysisTaskJetCorePP::~AliAnalysisTaskJetCorePP()
613{
614 //Destructor
615 delete fListJets;
80ac66f6 616 delete fListJetsGen;
160dbdb8 617 delete fListJetsGenFull;
255cb71c 618 delete fListJetsBg;
619 delete fListJetsBgGen;
ad869500 620 delete fOutputList; // ????
8e103a56 621 delete fRandom;
ad869500 622}
623
624//--------------------------------------------------------------
625
ea603b64 626
627Bool_t AliAnalysisTaskJetCorePP::Notify()
628{
629 //Implemented Notify() to read the cross sections
630 //and number of trials from pyxsec.root
631 //inspired by AliAnalysisTaskJetSpectrum2::Notify()
49dc2855 632 if(!(fIsChargedMC || fIsKine)) return kFALSE;
3d1fac8a 633 Float_t xsection = 0;
634 Float_t trials = 1;
635 fAvgTrials = 1;
ea603b64 636
3d1fac8a 637 if(fIsChargedMC){
638 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
639 if(!fESD){
640 if(fDebug>1) AliError("ESD not available");
641 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
642 }
ea603b64 643
3d1fac8a 644 fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
ea603b64 645
646
3d1fac8a 647 if(fNonStdFile.Length()!=0){
648 // case that we have an AOD extension we can fetch the jets from the extended output
649 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
650 fAODExtension = aodH ? aodH->GetExtension(fNonStdFile.Data()) : 0;
651 if(!fAODExtension){
652 if(fDebug>1) Printf("AODExtension found for %s",fNonStdFile.Data());
653 }
654 }
ea603b64 655
3d1fac8a 656 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
ea603b64 657
3d1fac8a 658 if(tree){
659 TFile *curfile = tree->GetCurrentFile();
660 if(!curfile) {
661 Error("Notify","No current file");
662 return kFALSE;
663 }
664 if(!fh1Xsec || !fh1Trials){
665 Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
666 return kFALSE;
667 }
668 AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,trials);
669 fh1Xsec->Fill("<#sigma>",xsection);
670 // construct a poor man average trials
671 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
672 if(trials>=nEntries && nEntries>0.) fAvgTrials = trials/nEntries;
673 fh1Trials->Fill("#sum{ntrials}",trials);
674 }
675
676 if(fDebug)Printf("Reading File %s",fInputHandler->GetTree()->GetCurrentFile()->GetName());
677 }
ea603b64 678
ea603b64 679
680 return kTRUE;
681}
682//--------------------------------------------------------------
683
ad869500 684void AliAnalysisTaskJetCorePP::Init()
685{
686 // check for jet branches
49dc2855 687 if(fJetBranchNameKine.Length()==0){
688 if(!strlen(fJetBranchName.Data())){
689 AliError("Jet branch name not set.");
690 }
ad869500 691 }
49dc2855 692 fMcHandler = dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
693
694
ad869500 695}
696
697//--------------------------------------------------------------
698
699void AliAnalysisTaskJetCorePP::UserCreateOutputObjects()
700{
8aa19603 701 // Create histograms and initilize variables
702 fSafetyMargin = fBgConeR*fBgConeR /(fBgJetParamR*fBgJetParamR);
49dc2855 703
ad869500 704 // Called once
255cb71c 705 fListJets = new TList(); //reconstructed level antikt jets
8aa19603 706 fListJetsBg = new TList(); //reconstructed jets to be removed from bg
80ac66f6 707
255cb71c 708 fIsChargedMC = (fJetBranchNameChargMC.Length()>0) ? kTRUE : kFALSE;
49dc2855 709 fIsKine = (fJetBranchNameKine.Length()>0) ? kTRUE : kFALSE;
255cb71c 710 fIsFullMC = (fJetBranchNameFullMC.Length()>0) ? kTRUE : kFALSE;
80ac66f6 711
8e103a56 712 fRandom = new TRandom3(0);
713
49dc2855 714 if(fIsChargedMC || fIsKine){ //full MC or pure kine
255cb71c 715 fListJetsGen = new TList(); //generator level charged antikt jets
8aa19603 716 fListJetsBgGen = new TList(); //generator level jets to be removed from bg
255cb71c 717
718 if(fIsFullMC)
160dbdb8 719 fListJetsGenFull = new TList(); //generator level full jets
720 }
ad869500 721 OpenFile(1);
722 if(!fOutputList) fOutputList = new TList();
723 fOutputList->SetOwner(kTRUE);
724
725 Bool_t oldStatus = TH1::AddDirectoryStatus();
726 TH1::AddDirectory(kFALSE);
727
728 fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
729 fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
730 fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
731 fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
732 fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
733 fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
8e103a56 734 fHistEvtSelection->GetXaxis()->SetBinLabel(6,"event number (rejected)");
ad869500 735
736 fOutputList->Add(fHistEvtSelection);
737
738 Int_t nBinsCentrality = (fSystem==0) ? 1 : 10; // pp=1 else 10
80ac66f6 739 //trigger pt spectrum (reconstructed)
ad869500 740 fh2Ntriggers = new TH2F("fh2Ntriggers","# of triggers",
741 nBinsCentrality,0.0,100.0,50,0.0,50.0);
49dc2855 742 if(!fIsKine) fOutputList->Add(fh2Ntriggers);
743
744 Int_t bw = (fDoubleBinning==0) ? 1 : 2; //make larger bin width
ad869500 745
6168afc3 746 //Centrality, A, pTjet, pTtrigg, dphi
747 const Int_t dimSpec = 5;
49dc2855 748 const Int_t nBinsSpec[dimSpec] = {nBinsCentrality, 50, bw*110, 50, TMath::Nint(10*(TMath::Pi()-fkDeltaPhiCut))};
749 const Double_t lowBinSpec[dimSpec] = {0.0, 0.0, -20, 0.0, fkDeltaPhiCut};
750 const Double_t hiBinSpec[dimSpec] = {100.0, 1.0, 200.0,50.0, TMath::Pi()};
7fc3d134 751 fHJetSpec = new THnSparseF("fHJetSpec",
6168afc3 752 "Recoil jet spectrum [cent,A,pTjet,pTtrig,dphi]",
7fc3d134 753 dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
49dc2855 754 if(!fIsKine) fOutputList->Add(fHJetSpec);
ad869500 755
255cb71c 756 //background estimated as median of kT jets
757 fHJetSpecSubUeMedian = (THnSparseF*) fHJetSpec->Clone("fHJetSpecSubUeMedian");
758 fHJetSpecSubUeMedian->SetTitle("Recoil jet spectrum [cent,A,pTjet-pTUe,pTtrig,dphi]");
49dc2855 759 if(!fIsKine) fOutputList->Add(fHJetSpecSubUeMedian);
8aa19603 760 //background estimated as weighted median of kT jets ala Cone
761 fHJetSpecSubUeCone = (THnSparseF*) fHJetSpec->Clone("fHJetSpecSubUeCone");
762 fHJetSpecSubUeCone->SetTitle("Recoil jet spectrum [cent,A,pTjet-pTUe,pTtrig,dphi]");
49dc2855 763 if(!fIsKine) fOutputList->Add(fHJetSpecSubUeCone);
764
255cb71c 765
49dc2855 766 //A, pTjet, pTtrigg, dphi
767 const Int_t dimCor = 4;
768 const Int_t nBinsCor[dimCor] = {50, 110, 50, 100};
769 const Double_t lowBinCor[dimCor] = {0.0, -20, 0.0, -0.5*TMath::Pi()};
770 const Double_t hiBinCor[dimCor] = {1.0, 200, 50.0, 1.5*TMath::Pi()};
771 fHJetPhiCorr = new THnSparseF("fHJetPhiCorr",
772 "Recoil jet spectrum [A,pTjet,pTtrig,dphi]",
773 dimCor,nBinsCor,lowBinCor,hiBinCor);
774
775 if(!fIsKine) fOutputList->Add(fHJetPhiCorr); // Dphi distribution jet-triger
776
777 //background estimated as median of kT jets
778 fHJetPhiCorrSubUeMedian = (THnSparseF*) fHJetPhiCorr->Clone("fHJetPhiCorrSubUeMedian");
779 fHJetPhiCorrSubUeMedian->SetTitle("Recoil jet spectrum [A,pTjet-pTUe,pTtrig,dphi]");
780 if(!fIsKine) fOutputList->Add(fHJetPhiCorrSubUeMedian);
781 //background estimated as weighted median of kT jets ala Cone
782 fHJetPhiCorrSubUeCone = (THnSparseF*) fHJetPhiCorr->Clone("fHJetPhiCorrSubUeCone");
783 fHJetPhiCorrSubUeCone->SetTitle("Recoil jet spectrum [A,pTjet-pTUe,pTtrig,dphi]");
784 if(!fIsKine) fOutputList->Add(fHJetPhiCorrSubUeCone);
255cb71c 785
786
ad869500 787 //------------------- HISTOS FOR DIAGNOSTIC ----------------------
255cb71c 788 //A, pTjet, pTjet-pTUe, pTUe, rhoUe bg estimate from kT median
789 const Int_t dimSpecMed = 5;
1afc79e9 790 const Int_t nBinsSpecMed[dimSpecMed] = {25, 50, 60, 20, 20};
8aa19603 791 const Double_t lowBinSpecMed[dimSpecMed] = {0.0, 0.0, -20.0, 0.0, 0.0};
1afc79e9 792 const Double_t hiBinSpecMed[dimSpecMed] = {1.0, 100.0, 100.0, 10.0, 20.0};
255cb71c 793 fHJetUeMedian = new THnSparseF("fHJetUeMedian",
794 "Recoil jet spectrum [A,pTjet,pTjet-pTUe, pTUe, rhoUe]",
795 dimSpecMed, nBinsSpecMed, lowBinSpecMed, hiBinSpecMed);
49dc2855 796 if(!fIsKine) fOutputList->Add(fHJetUeMedian);
255cb71c 797
8aa19603 798 //A, pTjet, pTjet-pTUe, pTUe, rhoUe bg estimate from kT median Cone
799 fHJetUeCone = (THnSparseF*) fHJetUeMedian->Clone("fHJetUeCone");
800 fHJetUeCone->SetTitle("Recoil jet spectrum [A,pTjet,pTjet-pTUe, pTUe, rhoUe]");
49dc2855 801 if(!fIsKine) fOutputList->Add(fHJetUeCone);
255cb71c 802
803 //rho bacground reconstructed data
804 const Int_t dimRho = 2;
1afc79e9 805 const Int_t nBinsRho[dimRho] = {50 , 50};
255cb71c 806 const Double_t lowBinRho[dimRho] = {0.0 , 0.0};
807 const Double_t hiBinRho[dimRho] = {20.0 , 20.0};
808
8aa19603 809 fHRhoUeMedianVsCone = new THnSparseF("hRhoUeMedianVsCone","[Rho Cone, Rho Median]",
255cb71c 810 dimRho, nBinsRho, lowBinRho, hiBinRho);
49dc2855 811 if(!fIsKine) fOutputList->Add(fHRhoUeMedianVsCone);
255cb71c 812
ad869500 813 //Jet number density histos [Trk Mult, jet density, pT trigger]
255cb71c 814 /*const Int_t dimJetDens = 3;
80ac66f6 815 const Int_t nBinsJetDens[dimJetDens] = {100, 100, 10};
ad869500 816 const Double_t lowBinJetDens[dimJetDens] = {0.0, 0.0, 0.0};
80ac66f6 817 const Double_t hiBinJetDens[dimJetDens] = {500.0, 5.0, 50.0};
ad869500 818
819 fHJetDensity = new THnSparseF("fHJetDensity","Jet dens vs trk mult A>0.07",
820 dimJetDens,nBinsJetDens,lowBinJetDens,hiBinJetDens);
821
822 fHJetDensityA4 =new THnSparseF("fHJetDensityA4","Jet dens vs trk mult A>0.4",
823 dimJetDens,nBinsJetDens,lowBinJetDens,hiBinJetDens);
824
ad869500 825 fOutputList->Add(fHJetDensity);
826 fOutputList->Add(fHJetDensityA4);
255cb71c 827 */
ad869500 828
829 //inclusive azimuthal and pseudorapidity histograms
7fc3d134 830 fhJetPhi = new TH2D("fhJetPhi","Azim dist jets vs pTjet",
831 50, 0, 100, 50,-TMath::Pi(),TMath::Pi());
832 fhTriggerPhi= new TH2D("fhTriggerPhi","azim dist trig had vs pTtrigg",
6168afc3 833 50, 0, 50, 50,-TMath::Pi(),TMath::Pi());
7fc3d134 834 fhJetEta = new TH2D("fhJetEta","Eta dist jets vs pTjet",
835 50,0, 100, 40,-0.9,0.9);
836 fhTriggerEta = new TH2D("fhTriggerEta","Eta dist trig had vs pTtrigg",
6168afc3 837 50, 0, 50, 40,-0.9,0.9);
7fc3d134 838
ad869500 839 fhVertexZ = new TH1D("fhVertexZ","z vertex",40,-20,20);
840 fhVertexZAccept = new TH1D("fhVertexZAccept","z vertex after cut",40,-20,20);
841 fhContribVtx = new TH1D("fhContribVtx","contrib to vtx",200,0,200);
842 fhContribVtxAccept = new TH1D("fhContribVtxAccept","contrib to vtx after cut",200,0,200);
843 fhDphiTriggerJet = new TH1D("fhDphiTriggerJet","Deltaphi trig-jet",50, -TMath::Pi(),TMath::Pi());
844 fhDphiTriggerJetAccept = new TH1D("fhDphiTriggerJetAccept","Deltaphi trig-jet after cut",50, -TMath::Pi(),TMath::Pi());
845 fhCentrality = new TH1D("fhCentrality","Centrality",20,0,100);
846 fhCentralityAccept = new TH1D("fhCentralityAccept","Centrality after cut",20,0,100);
d1405a52 847 fhEntriesToMedian = new TH1D("fhEntriesToMedian","fhEntriesToMedian",30,0,30);
49dc2855 848 fhCellAreaToMedian = new TH1D("fhCellAreaToMedian", "fhCellAreaToMedian", 75,0,1.5);
ad869500 849
30557bf5 850 fhNofMultipleTriggers = new TH1D("fhNofMultipleTriggers","fhNofMultipleTriggers",20,0,20);
851 fhNofMultipleTriggersCone = new TH1D("fhNofMultipleTriggersCone","fhNofMultipleTriggersCone R<0.4",20,0,20);
852 fhNofMultipleTriggersConeLow = (TH1D*) fhNofMultipleTriggersCone->Clone("fhNofMultipleTriggersConepTa15to20");
853 fhNofMultipleTriggersConeHigh = (TH1D*) fhNofMultipleTriggersCone->Clone("fhNofMultipleTriggersConepTa20to50");
364e98aa 854
855 fhDeltaRMultTriggersLow = new TH1D("fhDeltaRMultTriggersLow","fhDeltaRMultTriggersLow", 200,0,4);
856 fhDeltaRMultTriggersHigh = (TH1D*) fhDeltaRMultTriggersLow->Clone("fhDeltaRMultTriggersHigh");
857
f1e6baa6 858 fhDeltaPhiMultTriggersLow = new TH1D("fhDeltaPhiMultTriggersLow","fhDeltaPhiRultTriggers", 320,-3.2,3.2); //single incl trigger
859 fhDeltaPhiMultTriggersHigh = (TH1D*) fhDeltaPhiMultTriggersLow->Clone("fhDeltaPhiMultTriggersHigh"); //single incl trigger
860
861 fhDeltaPhiMultTriggersInclLow = (TH1D*) fhDeltaPhiMultTriggersLow->Clone("fhDeltaPhiMultTriggersInclLow");
862 fhDeltaPhiMultTriggersInclHigh = (TH1D*) fhDeltaPhiMultTriggersLow->Clone("fhDeltaPhiMultTriggersInclHigh");
863 fhInclTrigCounter = new TH1D("fhInclTrigCounter","fhInclTrigCounter",1,0,1);
1e3ad34d 864
49dc2855 865 if(!fIsKine){
866 fOutputList->Add(fhJetPhi);
867 fOutputList->Add(fhTriggerPhi);
868 fOutputList->Add(fhJetEta);
869 fOutputList->Add(fhTriggerEta);
870 }
ad869500 871 fOutputList->Add(fhVertexZ);
49dc2855 872 fOutputList->Add(fhVertexZAccept);
873 if(!fIsKine){
874 fOutputList->Add(fhContribVtx);
875 fOutputList->Add(fhContribVtxAccept);
876 fOutputList->Add(fhDphiTriggerJet);
877 fOutputList->Add(fhDphiTriggerJetAccept);
878 fOutputList->Add(fhCentrality);
879 fOutputList->Add(fhCentralityAccept);
880 fOutputList->Add(fhEntriesToMedian);
881 fOutputList->Add(fhCellAreaToMedian);
1e3ad34d 882 fOutputList->Add(fhNofMultipleTriggers);
0fe881ea 883 fOutputList->Add(fhNofMultipleTriggersCone);
30557bf5 884 fOutputList->Add(fhNofMultipleTriggersConeLow);
885 fOutputList->Add(fhNofMultipleTriggersConeHigh);
364e98aa 886 fOutputList->Add(fhDeltaRMultTriggersLow);
887 fOutputList->Add(fhDeltaRMultTriggersHigh);
f1e6baa6 888 fOutputList->Add(fhDeltaPhiMultTriggersLow);
889 fOutputList->Add(fhDeltaPhiMultTriggersHigh);
890 fOutputList->Add(fhDeltaPhiMultTriggersInclLow);
891 fOutputList->Add(fhDeltaPhiMultTriggersInclHigh);
892 fOutputList->Add(fhInclTrigCounter);
49dc2855 893 }
7fc3d134 894 // raw spectra of INCLUSIVE jets
80ac66f6 895 //Centrality, pTjet, A
255cb71c 896 /*const Int_t dimRaw = 3;
80ac66f6 897 const Int_t nBinsRaw[dimRaw] = {nBinsCentrality, 50, 100};
898 const Double_t lowBinRaw[dimRaw] = {0.0, 0.0, 0.0};
899 const Double_t hiBinRaw[dimRaw] = {100.0, 100, 1.0};
7fc3d134 900 fHJetPtRaw = new THnSparseF("fHJetPtRaw",
80ac66f6 901 "Incl. jet spectrum [cent,pTjet,A]",
7fc3d134 902 dimRaw,nBinsRaw,lowBinRaw,hiBinRaw);
903 fOutputList->Add(fHJetPtRaw);
904
905 // raw spectra of LEADING jets
80ac66f6 906 //Centrality, pTjet, A
7fc3d134 907 fHLeadingJetPtRaw = new THnSparseF("fHLeadingJetPtRaw",
80ac66f6 908 "Leading jet spectrum [cent,pTjet,A]",
7fc3d134 909 dimRaw,nBinsRaw,lowBinRaw,hiBinRaw);
910 fOutputList->Add(fHLeadingJetPtRaw);
911
912 // Dphi versus pT jet
913 //Centrality, Dphi=phiTrig-phiJet, pTjet, pTtrigg
914 const Int_t dimDp = 4;
915 const Int_t nBinsDp[dimDp] = {nBinsCentrality, 50, 50, 50};
916 const Double_t lowBinDp[dimDp] = {0.0, -TMath::Pi(), 0.0, 0.0};
917 const Double_t hiBinDp[dimDp] = {100.0, TMath::Pi(), 100.0, 100.0};
918 fHDphiVsJetPtAll = new THnSparseF("fHDphiVsJetPtAll",
919 "Dphi vs jet pT [cent,Dphi,pTjet,pTtrigg]",
920 dimDp,nBinsDp,lowBinDp,hiBinDp);
921 fOutputList->Add(fHDphiVsJetPtAll);
255cb71c 922 */
7fc3d134 923
80ac66f6 924 //analyze MC generator level
49dc2855 925 if(fIsChargedMC || fIsKine){
8aa19603 926 if(fFillRespMx){
927 //Fill response matrix only once
49dc2855 928 fhJetPtGenVsJetPtRec = new TH2D("fhJetPtGenVsJetPtRec","JetPtGenVsJetPtRec", bw*100,0,200, bw*100,0,200);
8aa19603 929 fOutputList->Add(fhJetPtGenVsJetPtRec); //gen MC charg jet pt spectrum versus rec charged jet pt spectrum
930 //....
49dc2855 931 fhJetPtGenVsJetPtRecSubUeMedian = new TH2D("fhJetPtGenVsJetPtRecSubUeMedian","fhJetPtGenVsJetPtRecSubUeMedian", bw*110,-20,200, bw*110,-20,200);
8aa19603 932 fOutputList->Add(fhJetPtGenVsJetPtRecSubUeMedian); // with kT median bg subtr
255cb71c 933 //....
8aa19603 934 fhJetPtGenVsJetPtRecSubUeCone=(TH2D*)fhJetPtGenVsJetPtRecSubUeMedian->Clone("fhJetPtGenVsJetPtRecSubUeCone");
935 fhJetPtGenVsJetPtRecSubUeCone->SetTitle("fhJetPtGenVsJetPtRecSubUeCone");
936 fOutputList->Add(fhJetPtGenVsJetPtRecSubUeCone); // with weighted kT median bg subtr
937 //....
49dc2855 938 fhJetPtGen = new TH1D("fhJetPtGen","Jet Pt (MC Gen)",bw*100,0,200); //MC generator charged jet pt spectrum
8aa19603 939 fOutputList->Add(fhJetPtGen);
940 //....
49dc2855 941 fhJetPtSubUeMedianGen = new TH1D("fhJetPtSubUeMedianGen","Jet Pt - UE pT (MC Gen)",bw*110,-20,200);
8aa19603 942 fOutputList->Add(fhJetPtSubUeMedianGen); // with kT median bg subtr
943 //....
944 fhJetPtSubUeConeGen = (TH1D*) fhJetPtSubUeMedianGen->Clone("fhJetPtSubUeConeGen");
945 fOutputList->Add(fhJetPtSubUeConeGen); // with weighted kT median bg subtr
f1e6baa6 946 //....
947 fhJetPtResolutionVsPtGen = new TH2D("fhJetPtResolutionVsPtGen","fhJetPtResolutionVsPtGen",20,0,100, 35,-1.,0.4);
948 fOutputList->Add(fhJetPtResolutionVsPtGen); // with weighted kT median bg subtr
949 //....
950 fhJetPtResolutionVsPtConeGen = (TH2D*) fhJetPtResolutionVsPtGen->Clone("fhJetPtResolutionVsPtConeGen");
951 fOutputList->Add(fhJetPtResolutionVsPtConeGen); // with weighted kT median bg subtr
952
8aa19603 953 //....
954 if(fIsFullMC){
1afc79e9 955 fhJetPtGenChargVsJetPtGenFull = new TH2D("fhJetPtGenChargVsJetPtGenFull","fhJetPtGenChargVsJetPtGenFull", 100,0,200, 100,0,200);
8aa19603 956 fOutputList->Add(fhJetPtGenChargVsJetPtGenFull); //gen full MC jet pt versus gen charged jet MC pt
957 //....
1afc79e9 958 fhJetPtGenFull = new TH1D("fhJetPtGenFull","Jet Pt (MC Full jets Gen)",100,0,200); //MC generator full jet pt spectrum
8aa19603 959 fOutputList->Add(fhJetPtGenFull);
960 }
255cb71c 961 }
962 //....
80ac66f6 963 fh2NtriggersGen = (TH2F*) fh2Ntriggers->Clone("fh2NtriggersGen");
964 fh2NtriggersGen->SetTitle(Form("%s Gen MC",fh2Ntriggers->GetTitle()));
965 fOutputList->Add(fh2NtriggersGen);
7fc3d134 966
80ac66f6 967 //Centrality, A, pT, pTtrigg
968 fHJetSpecGen = (THnSparseF*) fHJetSpec->Clone("fHJetSpecGen");
969 fHJetSpecGen->SetTitle(Form("%s Gen MC",fHJetSpec->GetTitle()));
970 fOutputList->Add(fHJetSpecGen);
971
255cb71c 972 fHJetSpecSubUeMedianGen = (THnSparseF*) fHJetSpecSubUeMedian->Clone("fHJetSpecSubUeMedianGen");
973 fHJetSpecSubUeMedianGen->SetTitle(Form("%s Gen MC",fHJetSpecSubUeMedian->GetTitle()));
974 fOutputList->Add(fHJetSpecSubUeMedianGen);
975
8aa19603 976 fHJetSpecSubUeConeGen = (THnSparseF*) fHJetSpecSubUeCone->Clone("fHJetSpecSubUeConeGen");
977 fHJetSpecSubUeConeGen->SetTitle(Form("%s Gen MC",fHJetSpecSubUeCone->GetTitle()));
978 fOutputList->Add(fHJetSpecSubUeConeGen);
49dc2855 979 //---
980 fHJetPhiCorrGen = (THnSparseF*) fHJetPhiCorr->Clone("fHJetPhiCorrGen");
981 fHJetPhiCorrGen->SetTitle(Form("%s Gen MC",fHJetPhiCorr->GetTitle()));
982 fOutputList->Add(fHJetPhiCorrGen);
983
984 fHJetPhiCorrSubUeMedianGen = (THnSparseF*) fHJetPhiCorrSubUeMedian->Clone("fHJetPhiCorrSubUeMedianGen");
985 fHJetPhiCorrSubUeMedianGen->SetTitle(Form("%s Gen MC",fHJetPhiCorrSubUeMedian->GetTitle()));
986 fOutputList->Add(fHJetPhiCorrSubUeMedianGen);
987
988 fHJetPhiCorrSubUeConeGen = (THnSparseF*) fHJetPhiCorrSubUeCone->Clone("fHJetPhiCorrSubUeConeGen");
989 fHJetPhiCorrSubUeConeGen->SetTitle(Form("%s Gen MC", fHJetPhiCorrSubUeCone->GetTitle()));
990 fOutputList->Add(fHJetPhiCorrSubUeConeGen);
991
255cb71c 992 //---
993 fHJetUeMedianGen = (THnSparseF*) fHJetUeMedian->Clone("fHJetUeMedianGen");
994 fHJetUeMedianGen->SetTitle(Form("%s Gen MC", fHJetUeMedian->GetTitle()));
995 fOutputList->Add(fHJetUeMedianGen);
996
8aa19603 997 fHJetUeConeGen = (THnSparseF*) fHJetUeCone->Clone("fHJetUeConeGen");
998 fHJetUeConeGen->SetTitle(Form("%s Gen MC", fHJetUeCone->GetTitle()));
999 fOutputList->Add(fHJetUeConeGen);
255cb71c 1000
8aa19603 1001 if(fFillRespMx){
1002 //track efficiency/contamination histograms eta versus pT
f1e6baa6 1003 Double_t bins [] = {0, 0.2,0.4,0.6, 0.8, 1., 1.2, 1.4, 1.6, 1.8, 2., 2.5, 3., 3.5, 4., 5., 6., 8., 10., 20., 50.};
1004 Int_t nbins = sizeof(bins)/sizeof(Double_t)-1;
1005 fhPtTrkTruePrimRec = new TH2D("fhPtTrkTruePrimRec","PtTrkTruePrimRec",nbins, bins, 18,-0.9,0.9);
8aa19603 1006 fOutputList->Add(fhPtTrkTruePrimRec);
80ac66f6 1007
8aa19603 1008 fhPtTrkTruePrimGen = (TH2D*) fhPtTrkTruePrimRec->Clone("fhPtTrkTruePrimGen");
1009 fhPtTrkTruePrimGen->SetTitle("PtTrkTruePrimGen");
1010 fOutputList->Add(fhPtTrkTruePrimGen);
80ac66f6 1011
8aa19603 1012 fhPtTrkSecOrFakeRec = (TH2D*) fhPtTrkTruePrimRec->Clone("fhPtTrkSecOrFakeRec");
1013 fhPtTrkSecOrFakeRec->SetTitle("PtTrkSecOrFakeRec");
1014 fOutputList->Add(fhPtTrkSecOrFakeRec);
1015 }
255cb71c 1016
8aa19603 1017 fHRhoUeMedianVsConeGen = (THnSparseF*) fHRhoUeMedianVsCone->Clone("hRhoUeMedianVsConeGen");
1018 fHRhoUeMedianVsConeGen->SetTitle(Form("%s Gen MC", fHRhoUeMedianVsCone->GetTitle()));
1019 fOutputList->Add(fHRhoUeMedianVsConeGen);
d1405a52
ML		 1020
1021 fhEntriesToMedianGen = (TH1D*) fhEntriesToMedian->Clone("fhEntriesToMedianGen");
1022 fhEntriesToMedianGen->SetTitle(Form("%s Gen MC", fhEntriesToMedian->GetTitle()));
1023 fOutputList->Add(fhEntriesToMedianGen);
1024
1025 fhCellAreaToMedianGen = (TH1D*) fhCellAreaToMedian->Clone("fhCellAreaToMedianGen");
1026 fhCellAreaToMedianGen->SetTitle(Form("%s Gen MC", fhCellAreaToMedian->GetTitle()));
1027 fOutputList->Add(fhCellAreaToMedianGen);
1e3ad34d 1028
1029 fhNofMultipleTriggersGen = (TH1D*) fhNofMultipleTriggers->Clone("fhNofMultipleTriggersGen");
1030 fOutputList->Add(fhNofMultipleTriggersGen);
1031
0fe881ea 1032 fhNofMultipleTriggersConeGen = (TH1D*) fhNofMultipleTriggersCone->Clone("fhNofMultipleTriggersConeGen");
1033 fOutputList->Add(fhNofMultipleTriggersConeGen);
1034
30557bf5 1035 fhNofMultipleTriggersConeGenLow = (TH1D*) fhNofMultipleTriggersCone->Clone("fhNofMultipleTriggersConeGenpta15to20");
1036 fOutputList->Add(fhNofMultipleTriggersConeGenLow);
1037
1038 fhNofMultipleTriggersConeGenHigh = (TH1D*) fhNofMultipleTriggersCone->Clone("fhNofMultipleTriggersConeGenpta20to50");
1039 fOutputList->Add(fhNofMultipleTriggersConeGenHigh);
1040
364e98aa 1041 fhDeltaRMultTriggersGenLow = (TH1D*) fhDeltaRMultTriggersLow->Clone("fhDeltaRMultTriggersGenLow");
1042 fOutputList->Add(fhDeltaRMultTriggersGenLow);
30557bf5 1043
364e98aa 1044 fhDeltaRMultTriggersGenHigh = (TH1D*) fhDeltaRMultTriggersLow->Clone("fhDeltaRMultTriggersGenHigh");
1045 fOutputList->Add(fhDeltaRMultTriggersGenHigh);
1e3ad34d 1046
f1e6baa6 1047 fhDeltaPhiMultTriggersGenLow = (TH1D*) fhDeltaPhiMultTriggersLow->Clone("fhDeltaPhiMultTriggersGenLow");
1048 fOutputList->Add(fhDeltaPhiMultTriggersGenLow);
1049
1050 fhDeltaPhiMultTriggersGenHigh = (TH1D*) fhDeltaPhiMultTriggersLow->Clone("fhDeltaPhiMultTriggersGenHigh");
1051 fOutputList->Add(fhDeltaPhiMultTriggersGenHigh);
30557bf5 1052
f1e6baa6 1053 fhDeltaPhiMultTriggersInclGenLow = (TH1D*) fhDeltaPhiMultTriggersInclLow->Clone("fhDeltaPhiMultTriggersInclGenLow");
1054 fOutputList->Add(fhDeltaPhiMultTriggersInclGenLow);
1055
1056 fhDeltaPhiMultTriggersInclGenHigh = (TH1D*) fhDeltaPhiMultTriggersInclHigh->Clone("fhDeltaPhiMultTriggersInclGenHigh");
1057 fOutputList->Add(fhDeltaPhiMultTriggersInclGenHigh);
1058
1059 fhInclTrigCounterGen= (TH1D*) fhInclTrigCounter->Clone("fhInclTrigCounterGen");
1060 fOutputList->Add(fhInclTrigCounterGen);
30557bf5 1061
6c1489e1 1062 fhNofMultipleTriggersConeGenA = (TH1D*) fhNofMultipleTriggersConeGen->Clone("fhNofMultipleTriggersConeGen10");
1063 fOutputList->Add(fhNofMultipleTriggersConeGenA);
1064
30557bf5 1065 fhNofMultipleTriggersConeGenALow = (TH1D*) fhNofMultipleTriggersConeGen->Clone("fhNofMultipleTriggersConeGen10pTa15to20");
1066 fOutputList->Add(fhNofMultipleTriggersConeGenALow);
1067
1068 fhNofMultipleTriggersConeGenAHigh = (TH1D*) fhNofMultipleTriggersConeGen->Clone("fhNofMultipleTriggersConeGen10pTa20to50");
1069 fOutputList->Add(fhNofMultipleTriggersConeGenAHigh);
1070
364e98aa 1071 fhDeltaRMultTriggersGenALow = (TH1D*) fhDeltaRMultTriggersLow->Clone("fhDeltaRMultTriggersGen10pTa15to20");
1072 fOutputList->Add(fhDeltaRMultTriggersGenALow);
1073
1074 fhDeltaPhiMultTriggersGenALow = (TH1D*) fhDeltaPhiMultTriggersLow->Clone("fhDeltaPhiMultTriggersGen10pTa15to20");
1075 fOutputList->Add(fhDeltaPhiMultTriggersGenALow);
1076
1077 fhDeltaRMultTriggersGenAHigh = (TH1D*) fhDeltaRMultTriggersLow->Clone("fhDeltaRMultTriggersGen10pTa20to50");
1078 fOutputList->Add(fhDeltaRMultTriggersGenAHigh);
1079
1080 fhDeltaPhiMultTriggersGenAHigh = (TH1D*) fhDeltaPhiMultTriggersLow->Clone("fhDeltaPhiMultTriggersGen10pTa20to50");
1081 fOutputList->Add(fhDeltaPhiMultTriggersGenAHigh);
1082
1083 fhNofTriggersGenA = (TH1D*) fhInclTrigCounter->Clone("fhNofTriggersGen10");
1084 fOutputList->Add(fhNofTriggersGenA);
1085
1086
6c1489e1 1087 fhNofMultipleTriggersConeGenB = (TH1D*) fhNofMultipleTriggersConeGen->Clone("fhNofMultipleTriggersConeGen5");
1088 fOutputList->Add(fhNofMultipleTriggersConeGenB);
1089
30557bf5 1090 fhNofMultipleTriggersConeGenBLow = (TH1D*) fhNofMultipleTriggersConeGen->Clone("fhNofMultipleTriggersConeGen5pTa15to20");
1091 fOutputList->Add(fhNofMultipleTriggersConeGenBLow);
1092
1093 fhNofMultipleTriggersConeGenBHigh = (TH1D*) fhNofMultipleTriggersConeGen->Clone("fhNofMultipleTriggersConeGen5pTa20to50");
1094 fOutputList->Add(fhNofMultipleTriggersConeGenBHigh);
1095
364e98aa 1096 fhDeltaRMultTriggersGenBLow = (TH1D*) fhDeltaRMultTriggersLow->Clone("fhDeltaRMultTriggersGen5pTa15to20");
1097 fOutputList->Add(fhDeltaRMultTriggersGenBLow);
1098
1099 fhDeltaPhiMultTriggersGenBLow = (TH1D*) fhDeltaPhiMultTriggersLow->Clone("fhDeltaPhiMultTriggersGen5pTa15to20");
1100 fOutputList->Add(fhDeltaPhiMultTriggersGenBLow);
1101
1102 fhDeltaRMultTriggersGenBHigh = (TH1D*) fhDeltaRMultTriggersLow->Clone("fhDeltaRMultTriggersGen5pTa20to50");
1103 fOutputList->Add(fhDeltaRMultTriggersGenBHigh);
1104
1105 fhDeltaPhiMultTriggersGenBHigh = (TH1D*) fhDeltaPhiMultTriggersLow->Clone("fhDeltaPhiMultTriggersGen5pTa20to50");
1106 fOutputList->Add(fhDeltaPhiMultTriggersGenBHigh);
1107
1108 fhNofTriggersGenB = (TH1D*) fhInclTrigCounter->Clone("fhNofTriggersGen5");
1109 fOutputList->Add(fhNofTriggersGenB);
1110
1111
1112
1113
1114 if(fIsChargedMC){ // Filled from gen level if well reconstructed trigger exists
1115 fhTriggerCounterGenLevel = (TH1D*) fhInclTrigCounter->Clone("fhTriggerCounterGenLevel");
1116 fOutputList->Add(fhTriggerCounterGenLevel);
30557bf5 1117
364e98aa 1118 fhDeltaRMultTriggersGenLevelLow = (TH1D*) fhDeltaRMultTriggersLow->Clone("fhDeltaRMultTriggersGenLevelLow");
1119 fOutputList->Add(fhDeltaRMultTriggersGenLevelLow);
30557bf5 1120
364e98aa 1121 fhDeltaPhiMultTriggersGenLevelLow = (TH1D*) fhDeltaPhiMultTriggersLow->Clone("fhDeltaPhiMultTriggersGenLevelLow");
1122 fOutputList->Add(fhDeltaPhiMultTriggersGenLevelLow);
1123
1124 fhDeltaRMultTriggersGenLevelHigh = (TH1D*) fhDeltaRMultTriggersLow->Clone("fhDeltaRMultTriggersGenLevelHigh");
1125 fOutputList->Add(fhDeltaRMultTriggersGenLevelHigh);
1126
1127 fhDeltaPhiMultTriggersGenLevelHigh = (TH1D*) fhDeltaPhiMultTriggersLow->Clone("fhDeltaPhiMultTriggersGenLevelHigh");
1128 fOutputList->Add(fhDeltaPhiMultTriggersGenLevelHigh);
1129 }
80ac66f6 1130 }
ea603b64 1131 //-------------------------------------
1132 // pythia histograms
1133 const Int_t nBinPt = 150;
1134 Double_t binLimitsPt[nBinPt+1];
1135 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
1136 if(iPt == 0){
1137 binLimitsPt[iPt] = -50.0;
1138 }else{// 1.0
1139 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.;
1140 }
1141 }
1142
1143 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
1144 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
1145 fOutputList->Add(fh1Xsec);
1146 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
1147 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
1148 fOutputList->Add(fh1Trials);
1149 fh1AvgTrials = new TH1F("fh1AvgTrials","trials root file",1,0,1);
1150 fh1AvgTrials->GetXaxis()->SetBinLabel(1,"#sum{avg ntrials}");
1151 fOutputList->Add(fh1AvgTrials);
1152 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
1153 fOutputList->Add(fh1PtHard);
1154 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
1155 fOutputList->Add(fh1PtHardNoW);
1156 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
1157 fOutputList->Add(fh1PtHardTrials);
1158
ad869500 1159
1160 // =========== Switch on Sumw2 for all histos ===========
1161 for(Int_t i=0; i<fOutputList->GetEntries(); i++){
1162 TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
1163 if(h1){
1164 h1->Sumw2();
1165 continue;
1166 }
1167 THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
1168 if(hn){
1169 hn->Sumw2();
1170 }
1171 }
1172 TH1::AddDirectory(oldStatus);
1173
1174 PostData(1, fOutputList);
1175}
1176
1177//--------------------------------------------------------------------
1178
1179void AliAnalysisTaskJetCorePP::UserExec(Option_t *)
1180{
3d1fac8a 1181 //User Exec
1182
3d1fac8a 1183
ad869500 1184 //Event loop
ea603b64 1185 Double_t eventW = 1.0;
1186 Double_t ptHard = 0.0;
1187 Double_t nTrials = 1.0; // Trials for MC trigger
49dc2855 1188 if(fIsChargedMC || fIsKine) fh1AvgTrials->Fill("#sum{avg ntrials}",fAvgTrials);
ad869500 1189
1190 if(TMath::Abs((Float_t) fJetParamR)<0.00001){
8aa19603 1191 AliError("ANTIKT Cone radius is set to zero.");
1192 return;
1193 }
49dc2855 1194
8aa19603 1195 if(TMath::Abs((Float_t) fBgJetParamR)<0.00001){
49dc2855 1196 AliError("ANTIKT Cone radius is set to zero.");
255cb71c 1197 return;
1198 }
1199
49dc2855 1200 if(!fIsKine){ //real data or full MC
1201 if(!strlen(fJetBranchName.Data())){
1202 AliError("Name of jet branch not set.");
1203 return;
1204 }
1205 if(!strlen(fJetBranchNameBg.Data())){
1206 AliError("Name of jet bg branch not set.");
1207 return;
1208 }
1209 }else{ //Kine
3d1fac8a 1210 if(!strlen(fJetBranchNameBgKine.Data())){
1211 AliError("Name of jet bg branch for kine not set.");
1212 return;
1213 }
49dc2855 1214
1215 Init();
1216 if(fMcHandler){
1217 fMcEvent = fMcHandler->MCEvent();
1218 }else{
3d1fac8a 1219 if(fDebug > 1) printf("AliAnalysisTaskJetCorePP::Exec() fMcHandler=NULL\n");
49dc2855 1220 PostData(1, fOutputList);
1221 return;
1222 }
1223 if(!fMcEvent){
3d1fac8a 1224 if(fDebug > 1) printf("AliAnalysisTaskJetCorePP::Exec() fMcEvent=NULL \n");
49dc2855 1225 PostData(1, fOutputList);
1226 return;
1227 }
29ac2c9b 1228
1229 Float_t xsection = 0;
1230 Float_t trials = 0;
1231
1232 AliGenPythiaEventHeader *genPH =
1233 dynamic_cast<AliGenPythiaEventHeader*> (fMcEvent->GenEventHeader());
1234 if(genPH){
1235 xsection = genPH->GetXsection();
1236 trials = genPH->Trials();
1237 ptHard = genPH->GetPtHard();
1238 }
1239 fh1Xsec->Fill("<#sigma>",xsection);
1240 fh1Trials->Fill("#sum{ntrials}",trials);
1241 fh1PtHard->Fill(ptHard,eventW);
1242 fh1PtHardNoW->Fill(ptHard,1);
1243 fh1PtHardTrials->Fill(ptHard,trials);
49dc2855 1244 }
ad869500 1245
29ac2c9b 1246
ad869500 1247 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
1248 if(!fESD){
80ac66f6 1249 if(fDebug>1) AliError("ESD not available");
ad869500 1250 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
1251 }
1252
1253 fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
1254 AliAODEvent* aod = NULL;
1255 // take all other information from the aod we take the tracks from
49dc2855 1256 if(!aod && !fIsKine){
1257 if(!fESD) aod = fAODIn; //not ESD and not kine => input AOD
1258 else aod = fAODOut;// ESD or kine
ad869500 1259 }
1260
49dc2855 1261
ad869500 1262 if(fNonStdFile.Length()!=0){
1263 // case that we have an AOD extension we can fetch the jets from the extended output
1264 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1265 fAODExtension = aodH ? aodH->GetExtension(fNonStdFile.Data()) : 0;
1266 if(!fAODExtension){
1267 if(fDebug>1) Printf("AODExtension found for %s",fNonStdFile.Data());
1268 }
1269 }
1270
80ac66f6 1271 // ----------------- EVENT SELECTION --------------------------
ad869500 1272 fHistEvtSelection->Fill(1); // number of events before event selection
1273
49dc2855 1274 if(!fIsKine){
1275 // physics selection
1276 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
ad869500 1277 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
1278
49dc2855 1279 if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
1280 if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
1281 fHistEvtSelection->Fill(2);
1282 PostData(1, fOutputList);
1283 return;
1284 }
1285
ad869500 1286 //check AOD pointer
49dc2855 1287 if(!aod){
1288 if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
1289 fHistEvtSelection->Fill(3);
1290 PostData(1, fOutputList);
1291 return;
1292 }
ad869500 1293
49dc2855 1294 // vertex selection for reconstructed data
1295 AliAODVertex* primVtx = aod->GetPrimaryVertex();
1296
1297 if(!primVtx){
1298 if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
1299 fHistEvtSelection->Fill(3);
1300 PostData(1, fOutputList);
1301 return;
1302 }
ad869500 1303
49dc2855 1304 Int_t nTracksPrim = primVtx->GetNContributors();
1305 Float_t vtxz = primVtx->GetZ();
1306 //Input events
1307 fhContribVtx->Fill(nTracksPrim);
1308 if( nTracksPrim > 0 ) fhVertexZ->Fill(vtxz);
ad869500 1309
49dc2855 1310 if((nTracksPrim < fMinContribVtx) ||
1311 (vtxz < fVtxZMin) ||
1312 (vtxz > fVtxZMax)){
1313 if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",
ad869500 1314 (char*)__FILE__,__LINE__,vtxz);
49dc2855 1315 fHistEvtSelection->Fill(3);
1316 PostData(1, fOutputList);
1317 return;
1318 }else{
ad869500 1319 //Accepted events
49dc2855 1320 fhContribVtxAccept->Fill(nTracksPrim);
1321 fhVertexZAccept->Fill(vtxz);
1322 }
1323 }else{ //KINE cut on vertex
1324 const AliVVertex *vtxMC = fMcEvent->GetPrimaryVertex();
1325 Float_t zVtx = vtxMC->GetZ();
1326 fhVertexZ->Fill(zVtx);
1327 if(zVtx < fVtxZMin || zVtx > fVtxZMax){ //vertex cut
1328 fHistEvtSelection->Fill(3);
1329 PostData(1, fOutputList);
1330 return;
1331 }else{
1332 fhVertexZAccept->Fill(zVtx);
1333 }
ad869500 1334 }
1335 //FK// No event class selection imposed
1336 // event class selection (from jet helper task)
1337 //Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
1338 //if(fDebug) Printf("Event class %d", eventClass);
1339 //if(eventClass < fEvtClassMin || eventClass > fEvtClassMax){
1340 // fHistEvtSelection->Fill(4);
1341 // PostData(1, fOutputList);
1342 // return;
1343 //}
1344
80ac66f6 1345 //------------------ CENTRALITY SELECTION ---------------
ad869500 1346 AliCentrality *cent = 0x0;
1347 Double_t centValue = 0.0;
1348 if(fSystem){ //fSystem=0 for pp, fSystem=1 for pPb
1349 if(fESD){
1350 cent = fESD->GetCentrality();
1351 if(cent) centValue = cent->GetCentralityPercentile("V0M");
1352 }else{
0fe881ea 1353 //centValue = aod->GetHeader()->GetCentrality();
0a918d8d 1354 centValue = ((AliVAODHeader*)aod->GetHeader())->GetCentrality();
ad869500 1355 }
1356 if(fDebug) printf("centrality: %f\n", centValue);
1357 //Input events
1358 fhCentrality->Fill(centValue);
1359
1360 if(centValue < fCentMin || centValue > fCentMax){
1361 fHistEvtSelection->Fill(4);
1362 PostData(1, fOutputList);
1363 return;
1364 }else{
1365 //Accepted events
1366 fhCentralityAccept->Fill( centValue );
1367 }
1368 }
1369
8e103a56 1370 //-----------------select disjunct event subsamples ----------------
49dc2855 1371 if(!fIsKine){ //reconstructed data
0fe881ea 1372 //Int_t eventnum = aod->GetHeader()->GetEventNumberESDFile();
0a918d8d 1373 AliAODHeader * header = dynamic_cast<AliAODHeader*>(aod->GetHeader());
1374 if(!header) AliFatal("Not a standard AOD");
1375
1376 Int_t eventnum = header->GetEventNumberESDFile();
49dc2855 1377 Int_t lastdigit = eventnum % 10;
1378 if(!(fEventNumberRangeLow<=lastdigit && lastdigit<=fEventNumberRangeHigh)){
1379 fHistEvtSelection->Fill(5);
1380 PostData(1, fOutputList);
1381 return;
0fe881ea 1382 }
49dc2855 1383 }
8e103a56 1384
ad869500 1385 if(fDebug) std::cout<<" ACCEPTED EVENT "<<endl;
ad869500 1386 fHistEvtSelection->Fill(0); // accepted events
255cb71c 1387 // ==================== end event selection ============================
1388
1389 Double_t tmpArrayFive[5];
49dc2855 1390 Double_t tmpArrayFour[4];
ad869500 1391
ad869500 1392
8aa19603 1393 TList particleList; //list of tracks
49dc2855 1394 Int_t indexTrigg = -1;
1395 Double_t rhoFromCellMedian=0.0, rhoCone=0.0;
80ac66f6 1396
49dc2855 1397 if(!fIsKine){
1398 //=============== Reconstructed antikt jets ===============
1399 ReadTClonesArray(fJetBranchName.Data() , fListJets);
1400 ReadTClonesArray(fJetBranchNameBg.Data(), fListJetsBg);
80ac66f6 1401
49dc2855 1402 //============ Estimate background in reconstructed events ===========
1403
1404 //Find Hadron trigger and Estimate rho from cone
1405 indexTrigg = GetListOfTracks(&particleList); //index of trigger hadron in Particle list
1406 EstimateBgCone(fListJets, &particleList, rhoCone);
1407
1408 //Estimate rho from cell median minus jets
1409 EstimateBgRhoMedian(fListJetsBg, &particleList, rhoFromCellMedian,0); //real data
1410 }
255cb71c 1411 //============== analyze generator level MC ================
1412 TList particleListGen; //list of tracks in MC
80ac66f6 1413
49dc2855 1414 if(fIsChargedMC || fIsKine){
1415
1416 if(fIsKine){ //pure kine
160dbdb8 1417
49dc2855 1418 //================= generated charged antikt jets ================
1419 ReadTClonesArray(fJetBranchNameKine.Data(), fListJetsGen);
3d1fac8a 1420 ReadTClonesArray(fJetBranchNameBgKine.Data(), fListJetsBgGen);
49dc2855 1421 }else{
1422 //================= generated charged antikt jets ================
1423 ReadTClonesArray(fJetBranchNameChargMC.Data(), fListJetsGen);
1424 ReadTClonesArray(fJetBranchNameBgChargMC.Data(), fListJetsBgGen);
1425
1426 if(fIsFullMC){ //generated full jets
1427 ReadTClonesArray(fJetBranchNameFullMC.Data(), fListJetsGenFull);
1428 }
1429 }
255cb71c 1430 //========================================================
80ac66f6 1431 //serarch for charged trigger at the MC generator level
255cb71c 1432
8e103a56 1433 Int_t indexTriggGen = -1;
1434 Double_t ptTriggGen = -1;
255cb71c 1435 Int_t iCounterGen = 0; //number of entries in particleListGen array
8e103a56 1436 Int_t triggersMC[200];//list of trigger candidates
1437 Int_t ntriggersMC = 0; //index in triggers array
6c1489e1 1438 Int_t triggersMCa[200]; //list of trigger candidates 10%eloss
1439 Int_t ntriggersMCa = 0; //index in triggers array 10%eloss
1440 Int_t triggersMCb[200]; //list of trigger candidates 5%eloss
1441 Int_t ntriggersMCb = 0; //index in triggers array 5%eloss
1442
8e103a56 1443
49dc2855 1444 if(!fIsKine){
1445 if(fESD){//ESD input
80ac66f6 1446
49dc2855 1447 AliMCEvent* mcEvent = MCEvent();
1448 if(!mcEvent){
1449 PostData(1, fOutputList);
1450 return;
1451 }
1452
1453 AliGenPythiaEventHeader *pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
1454 if(pythiaGenHeader){
1455 nTrials = pythiaGenHeader->Trials();
1456 ptHard = pythiaGenHeader->GetPtHard();
1457 fh1PtHard->Fill(ptHard,eventW);
1458 fh1PtHardNoW->Fill(ptHard,1);
1459 fh1PtHardTrials->Fill(ptHard,nTrials);
1460 }
ea603b64 1461
49dc2855 1462 for(Int_t it = 0; it < mcEvent->GetNumberOfTracks(); it++){
1463 if(!mcEvent->IsPhysicalPrimary(it)) continue;
1464 AliMCParticle* part = (AliMCParticle*) mcEvent->GetTrack(it);
1465 if(!part) continue;
1466 if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){
1467
1468 if(fHardest==0 && ntriggersMC<200){//single inclusive trigger
1469 if(indexTriggGen > -1){//trigger candidate was found
1470 triggersMC[ntriggersMC] = indexTriggGen;
1471 ntriggersMC++;
1472 }
1473 }
1474
1475 iCounterGen++;//index in particleListGen array
1476 }
1477 }
1478 }else{ //AOD input
1479 if(!fAODIn){
1480 PostData(1, fOutputList);
1481 return;
1482 }
1483 TClonesArray *tca = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(AliAODMCParticle::StdBranchName()));
1484 if(!tca){
1485 PostData(1, fOutputList);
1486 return;
1487 }
1488 for(Int_t it = 0; it < tca->GetEntriesFast(); it++){
1489 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1490 if(!part) continue;
1491 if(!part->IsPhysicalPrimary()) continue;
1492 if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){
1493
1494 if(fHardest==0 && ntriggersMC<200){//single inclusive trigger
1495 if(indexTriggGen > -1){ //trigger candidater was found
1496 triggersMC[ntriggersMC] = indexTriggGen;
1497 ntriggersMC++;
1498 }
1499 }
1500
1501 iCounterGen++;//number of entries in particleListGen array
1502 }
1503 }
ea603b64 1504 }
49dc2855 1505 }else{ //analyze kine
ea603b64 1506
49dc2855 1507 for(Int_t it = 0; it < fMcEvent->GetNumberOfTracks(); it++){
1508 if(!fMcEvent->IsPhysicalPrimary(it)) continue;
1509 AliMCParticle* part = (AliMCParticle*) fMcEvent->GetTrack(it);
80ac66f6 1510 if(!part) continue;
1511 if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){
49dc2855 1512
8e103a56 1513 if(fHardest==0 && ntriggersMC<200){//single inclusive trigger
255cb71c 1514 if(indexTriggGen > -1){//trigger candidate was found
8e103a56 1515 triggersMC[ntriggersMC] = indexTriggGen;
1516 ntriggersMC++;
1517 }
1518 }
49dc2855 1519
255cb71c 1520 iCounterGen++;//index in particleListGen array
80ac66f6 1521 }
1522 }
80ac66f6 1523 }
6c1489e1 1524
1525 if(fHardest==0){
1526 Int_t npar = particleListGen.GetEntries();
1527 for(Int_t ip=0; ip < npar; ip++){
1528 AliVParticle *part = (AliVParticle*) particleListGen.At(ip);
1529 if(!part) continue;
1530
1531 Double_t pta = 0.9 * part->Pt(); //10% energy loss
1532 Double_t ptb = 0.95 * part->Pt(); //5% energy loss
f1e6baa6 1533 if(fTriggerPtRangeLow <= pta && pta < fTriggerPtRangeHigh &&
1534 TMath::Abs(part->Eta()) < fTriggerEtaCut && ntriggersMCa<200){
6c1489e1 1535 triggersMCa[ntriggersMCa] = ip;
1536 ntriggersMCa++;
1537 }
1538
f1e6baa6 1539 if(fTriggerPtRangeLow <= ptb && ptb < fTriggerPtRangeHigh &&
1540 TMath::Abs(part->Eta()) < fTriggerEtaCut && ntriggersMCb<200){
6c1489e1 1541 triggersMCb[ntriggersMCb] = ip;
1542 ntriggersMCb++;
1543 }
1544 }
1545
1546 if(ntriggersMCa>0){
1547 Int_t rnda = fRandom->Integer(ntriggersMCa); //0 to ntriggers-1
1548 Int_t indexTriggGena = triggersMCa[rnda];
1549
1550 Double_t deltaPhia, deltaEtaa, deltaRa;
30557bf5 1551 Int_t aaLow = 0;
1552 Int_t aaHigh = 0;
6c1489e1 1553
1554 //Correlation with single inclusive TRIGGER
1555 AliVParticle* tGenTa = (AliVParticle*) particleListGen.At(indexTriggGena);
1556 if(tGenTa){
364e98aa 1557 fhNofTriggersGenA->Fill(0.5); // 15-50
58003b11 1558 //for(Int_t ia=0; ia<ntriggersMCa; ia++)
1559 for(Int_t ia=0; ia< npar; ia++){
1560 //if(indexTriggGena == triggersMCa[ia]) continue;
1561 if(indexTriggGena == ia) continue;
6c1489e1 1562
58003b11 1563 //AliVParticle* tGenTz = (AliVParticle*) particleListGen.At(triggersMCa[ia]);
1564 AliVParticle* tGenTz = (AliVParticle*) particleListGen.At(ia);
6c1489e1 1565 if(!tGenTz) continue;
58003b11 1566 if(tGenTz->Pt()*0.9<15.0) continue;
30557bf5 1567 if(tGenTz->Pt()*0.9>50.0) continue;
6c1489e1 1568
1569 deltaPhia = RelativePhi(tGenTa->Phi(),tGenTz->Phi());
1570 deltaEtaa = tGenTa->Eta()-tGenTz->Eta();
1571 deltaRa = sqrt(deltaPhia*deltaPhia + deltaEtaa*deltaEtaa);
364e98aa 1572
1573 if(tGenTz->Pt()*0.9<20.0){
1574 fhDeltaRMultTriggersGenALow->Fill(deltaRa);
1575 fhDeltaPhiMultTriggersGenALow->Fill(deltaPhia);
1576 }else{
1577 fhDeltaRMultTriggersGenAHigh->Fill(deltaRa);
1578 fhDeltaPhiMultTriggersGenAHigh->Fill(deltaPhia);
1579 }
1580
30557bf5 1581 if(deltaRa<0.4){
1582 if(tGenTz->Pt()*0.9<20.0) aaLow++; //15-20
1583 else aaHigh++; //20-50
1584 }
6c1489e1 1585 }
1586 }
30557bf5 1587 fhNofMultipleTriggersConeGenA->Fill(aaLow+aaHigh); // 15-50
1588 fhNofMultipleTriggersConeGenALow->Fill(aaLow); //15-20
1589 fhNofMultipleTriggersConeGenAHigh->Fill(aaHigh);//20-50
6c1489e1 1590 }
1591
1592 if(ntriggersMCb>0){
1593 Int_t rndb = fRandom->Integer(ntriggersMCb); //0 to ntriggers-1
1594 Int_t indexTriggGenb = triggersMCb[rndb];
1595
1596 Double_t deltaPhib, deltaEtab, deltaRb;
30557bf5 1597 Int_t bbLow = 0;
1598 Int_t bbHigh = 0;
6c1489e1 1599
1600 //Correlation with single inclusive TRIGGER
1601 AliVParticle* tGenTb = (AliVParticle*) particleListGen.At(indexTriggGenb);
1602 if(tGenTb){
364e98aa 1603 fhNofTriggersGenB->Fill(0.5); // 20-50
58003b11 1604 //for(Int_t ib=0; ib<ntriggersMCb; ib++)
1605 for(Int_t ib=0; ib<npar; ib++){
1606 //if(indexTriggGenb == triggersMCb[ib]) continue;
1607 if(indexTriggGenb == ib) continue;
6c1489e1 1608
58003b11 1609 //AliVParticle* tGenTz = (AliVParticle*) particleListGen.At(triggersMCb[ib]);
1610 AliVParticle* tGenTz = (AliVParticle*) particleListGen.At(ib);
6c1489e1 1611 if(!tGenTz) continue;
58003b11 1612 if(tGenTz->Pt()*0.95<15.0) continue;
30557bf5 1613 if(tGenTz->Pt()*0.95>50.0) continue;
6c1489e1 1614
1615 deltaPhib = RelativePhi(tGenTb->Phi(),tGenTz->Phi());
1616 deltaEtab = tGenTb->Eta()-tGenTz->Eta();
1617 deltaRb = sqrt(deltaPhib*deltaPhib + deltaEtab*deltaEtab);
364e98aa 1618
1619 if(tGenTz->Pt()*0.95<20.0){
1620 fhDeltaRMultTriggersGenBLow->Fill(deltaRb);
1621 fhDeltaPhiMultTriggersGenBLow->Fill(deltaPhib);
1622 }else{
1623 fhDeltaRMultTriggersGenBHigh->Fill(deltaRb);
1624 fhDeltaPhiMultTriggersGenBHigh->Fill(deltaPhib);
1625 }
6c1489e1 1626
30557bf5 1627 if(deltaRb<0.4){
1628 if(tGenTz->Pt()*0.95<20.0) bbLow++; //15-20
1629 else bbHigh++;
1630 }
6c1489e1 1631 }
1632 }
30557bf5 1633 fhNofMultipleTriggersConeGenB->Fill(bbLow+bbHigh); //15-50
1634 fhNofMultipleTriggersConeGenBLow->Fill(bbLow);//15-20
1635 fhNofMultipleTriggersConeGenBHigh->Fill(bbHigh);//20-50
6c1489e1 1636 }
1637 }
1638
1639
1640
8aa19603 1641 //============== Estimate bg in generated events ==============
1642 Double_t rhoFromCellMedianGen=0.0, rhoConeGen=0.0;
1643
1644 //Estimate rho from cone
1645 EstimateBgCone(fListJetsGen, &particleListGen, rhoConeGen);
1646
1647 //Estimate rho from cell median minus jets
d1405a52 1648 EstimateBgRhoMedian(fListJetsBgGen, &particleListGen, rhoFromCellMedianGen,1);//mc data
8aa19603 1649
f1e6baa6 1650 if(ntriggersMC > 0){ //select random inclusive trigger
1651 fhInclTrigCounterGen->Fill(0.5,ntriggersMC); //count inclusive triggers
1652
1653 //Lower pTa bin 15-20
1654 for(Int_t it1=0; it1<ntriggersMC; it1++){
1655 AliVParticle* tGent1 = (AliVParticle*) particleListGen.At(triggersMC[it1]);
1656 if(!tGent1) continue;
1657 for(Int_t ia=0; ia<particleListGen.GetEntries(); ia++){
1658 if(ia == triggersMC[it1]) continue;
1659 AliVParticle* tGent2 = (AliVParticle*) particleListGen.At(ia);
1660 if(!tGent2) continue;
1661 if(tGent2->Pt()<15.0) continue;
1662 if(tGent2->Pt()>20.0) continue;
1663 fhDeltaPhiMultTriggersInclGenLow->Fill(RelativePhi(tGent1->Phi(),tGent2->Phi()));
1664 }
1665 }
1666 //Higher pTa bin 20-50
1667 for(Int_t it1=0; it1<ntriggersMC-1; it1++){
1668 AliVParticle* tGent1 = (AliVParticle*) particleListGen.At(triggersMC[it1]);
1669 if(!tGent1) continue;
1670
1671 for(Int_t it2=it1+1; it2<ntriggersMC; it2++){
1672 AliVParticle* tGent2 = (AliVParticle*) particleListGen.At(triggersMC[it2]);
1673 if(!tGent2) continue;
1674 fhDeltaPhiMultTriggersInclGenHigh->Fill(RelativePhi(tGent1->Phi(),tGent2->Phi()));
1675 }
1676 }
1677 }
1678
1679
255cb71c 1680 //============ Generator trigger+jet ==================
8e103a56 1681 if(fHardest==0){ //single inclusive trigger
1682 if(ntriggersMC>0){ //there is at least one trigger
1683 Int_t rnd = fRandom->Integer(ntriggersMC); //0 to ntriggers-1
1684 indexTriggGen = triggersMC[rnd];
1e3ad34d 1685
1686 fhNofMultipleTriggersGen->Fill(ntriggersMC-1);
0fe881ea 1687
ac4f24f7 1688 Double_t deltaPhi, deltaEta, deltaR;
30557bf5 1689 Int_t iLow = 0;
1690 Int_t iHigh = 0;
ac4f24f7 1691
1692 //Correlation with single inclusive TRIGGER
1693 AliVParticle* tGenT1 = (AliVParticle*) particleListGen.At(indexTriggGen);
1694 if(tGenT1){
58003b11 1695 //for(Int_t ia=0; ia<ntriggersMC; ia++)
1696 for(Int_t ia=0; ia< particleListGen.GetEntries(); ia++){
1697 //if(indexTriggGen == triggersMC[ia]) continue;
1698 if(indexTriggGen == ia) continue;
ac4f24f7 1699
58003b11 1700 //AliVParticle* tGenT2 = (AliVParticle*) particleListGen.At(triggersMC[ia]);
1701 AliVParticle* tGenT2 = (AliVParticle*) particleListGen.At(ia);
ac4f24f7 1702 if(!tGenT2) continue;
58003b11 1703 if(tGenT2->Pt()<15.0) continue;
30557bf5 1704 if(tGenT2->Pt()>50.0) continue;
ac4f24f7 1705 deltaPhi = RelativePhi(tGenT1->Phi(),tGenT2->Phi());
1706 deltaEta = tGenT1->Eta()-tGenT2->Eta();
1707 deltaR = sqrt(deltaPhi*deltaPhi + deltaEta*deltaEta);
0fe881ea 1708
f1e6baa6 1709 if(tGenT2->Pt()<20.0){ //single inclusive trigger + assoc 15-20
364e98aa 1710 fhDeltaRMultTriggersGenLow->Fill(deltaR);
f1e6baa6 1711 fhDeltaPhiMultTriggersGenLow->Fill(deltaPhi);
1712 }else{ //single inclusive trigger + assoc 20-50
364e98aa 1713 fhDeltaRMultTriggersGenHigh->Fill(deltaR);
f1e6baa6 1714 fhDeltaPhiMultTriggersGenHigh->Fill(deltaPhi);
1715 }
30557bf5 1716
1717 if(deltaR<0.4){
1718 if(tGenT2->Pt()<20.0) iLow++;
1719 else iHigh++;
1e3ad34d 1720 }
1721 }
30557bf5 1722 }
1723 fhNofMultipleTriggersConeGen->Fill(iLow+iHigh);
1724 fhNofMultipleTriggersConeGenLow->Fill(iLow);
1725 fhNofMultipleTriggersConeGenHigh->Fill(iHigh);
1e3ad34d 1726
f1e6baa6 1727 //Single inclusive trigger within |Delta eta
1728
1729
80ac66f6 1730 }else{
8e103a56 1731 indexTriggGen = -1; //trigger not found
80ac66f6 1732 }
1733 }
8e103a56 1734
255cb71c 1735 //-----------
8e103a56 1736 Int_t ilowGen = (fHardest==0 || fHardest==1) ? indexTriggGen : 0;
1737 Int_t ihighGen = (fHardest==0 || fHardest==1) ? indexTriggGen+1 : particleListGen.GetEntries();
255cb71c 1738 Bool_t fillOnceGen = kTRUE;
1739 //-----------
1740
8e103a56 1741 for(Int_t igen= ilowGen; igen< ihighGen; igen++){ //loop over possible trigger
1742 indexTriggGen = igen; //trigger hadron
1743
1744 if(indexTriggGen == -1) continue;
1745 AliVParticle* triggerGen=(AliVParticle*) particleListGen.At(indexTriggGen);
1746 if(!triggerGen) continue;
1747
1748 if(fHardest >= 2){
6168afc3 1749 if(triggerGen->Pt() < 6.0) continue; //all hadrons pt>6
8e103a56 1750 }
1751 if(TMath::Abs((Float_t) triggerGen->Eta()) > fTriggerEtaCut) continue;
1752
1753 ptTriggGen = triggerGen->Pt(); //count triggers
1754 fh2NtriggersGen->Fill(centValue, ptTriggGen);
1755
80ac66f6 1756 //Count jets and trigger-jet pairs at MC generator level
1757 for(Int_t ij=0; ij<fListJetsGen->GetEntries(); ij++){
1758 AliAODJet* jet = (AliAODJet*)(fListJetsGen->At(ij));
1759 if(!jet) continue;
1760 Double_t etaJetGen = jet->Eta();
255cb71c 1761 Double_t areaJetGen = jet->EffectiveAreaCharged();
80ac66f6 1762
1763 if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){
80ac66f6 1764
49dc2855 1765 Double_t ptUeFromCellMedianGen = rhoFromCellMedianGen*areaJetGen;
1766 Double_t ptUeConeGen = rhoConeGen*areaJetGen;
1767
1768 //Rongrong's analysis
1769 Double_t dPhiGen = jet->Phi() - triggerGen->Phi();
1770 if(dPhiGen>2*TMath::Pi()) dPhiGen -= 2*TMath::Pi();
1771 if(dPhiGen<-2*TMath::Pi()) dPhiGen += 2*TMath::Pi();
1772 if(dPhiGen<-0.5*TMath::Pi()) dPhiGen += 2*TMath::Pi();
1773 if(dPhiGen>1.5*TMath::Pi()) dPhiGen -= 2*TMath::Pi();
1774
1775 //[A,pTjet,pTtrig,dphi]
1776 tmpArrayFour[0] = areaJetGen;
1777 tmpArrayFour[1] = jet->Pt();
1778 tmpArrayFour[2] = ptTriggGen;
1779 tmpArrayFour[3] = dPhiGen;
1780
1781 fHJetPhiCorrGen->Fill(tmpArrayFour); // Dphi distribution jet-triger
1782
1783 //[A,pTjet-pTUe,pTtrig,dphi]
1784 tmpArrayFour[1] = jet->Pt() - ptUeFromCellMedianGen;
1785 fHJetPhiCorrSubUeMedianGen->Fill(tmpArrayFour);
1786
1787 //[A,pTjet-pTUe,pTtrig,dphi]
1788 tmpArrayFour[1] = jet->Pt() - ptUeConeGen;
1789 fHJetPhiCorrSubUeConeGen->Fill(tmpArrayFour);
1790
1791 //Leticia's analysis
8e103a56 1792 Double_t dphi = RelativePhi(triggerGen->Phi(), jet->Phi());
1793 if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue;
80ac66f6 1794
8e103a56 1795 //Centrality, A, pT, pTtrigg
255cb71c 1796 tmpArrayFive[0] = centValue;
1797 tmpArrayFive[1] = areaJetGen;
1798 tmpArrayFive[2] = jet->Pt();
1799 tmpArrayFive[3] = ptTriggGen;
1800 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1801 fHJetSpecGen->Fill(tmpArrayFive);
1802
8aa19603 1803 //Centrality, A, pTjet-pTbgCellMedian, pTtrigg, dphi
255cb71c 1804 tmpArrayFive[0] = centValue;
1805 tmpArrayFive[1] = areaJetGen;
8aa19603 1806 tmpArrayFive[2] = jet->Pt() - ptUeFromCellMedianGen;
255cb71c 1807 tmpArrayFive[3] = ptTriggGen;
1808 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1809 fHJetSpecSubUeMedianGen->Fill(tmpArrayFive);
1810
8aa19603 1811 //Centrality, A, pTjet-pTbgCone, pTtrigg, dphi
255cb71c 1812 tmpArrayFive[0] = centValue;
1813 tmpArrayFive[1] = areaJetGen;
8aa19603 1814 tmpArrayFive[2] = jet->Pt() - ptUeConeGen;
255cb71c 1815 tmpArrayFive[3] = ptTriggGen;
1816 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
8aa19603 1817 fHJetSpecSubUeConeGen->Fill(tmpArrayFive);
255cb71c 1818
1819 //Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", kT median
1820 tmpArrayFive[0] = areaJetGen;
1821 tmpArrayFive[1] = jet->Pt();
8aa19603 1822 tmpArrayFive[2] = jet->Pt() - ptUeFromCellMedianGen;
1823 tmpArrayFive[3] = ptUeFromCellMedianGen;
1824 tmpArrayFive[4] = rhoFromCellMedianGen;
255cb71c 1825 fHJetUeMedianGen->Fill(tmpArrayFive);
1826
8aa19603 1827 //Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", perpendicular Cone
255cb71c 1828 tmpArrayFive[0] = areaJetGen;
1829 tmpArrayFive[1] = jet->Pt();
8aa19603 1830 tmpArrayFive[2] = jet->Pt() - ptUeConeGen;
1831 tmpArrayFive[3] = ptUeConeGen;
1832 tmpArrayFive[4] = rhoConeGen;
1833 fHJetUeConeGen->Fill(tmpArrayFive);
255cb71c 1834
1835 if(fillOnceGen){
8aa19603 1836 Double_t fillRhoGen[] = { rhoConeGen, rhoFromCellMedianGen};
1837 fHRhoUeMedianVsConeGen->Fill(fillRhoGen);
255cb71c 1838 fillOnceGen = kFALSE;
1839 }
8e103a56 1840 }//back to back jet-trigger pair
1841 }//jet loop
1842 }//trigger loop
1843
49dc2855 1844 if(fFillRespMx && !fIsKine){
8aa19603 1845 //================ RESPONSE MATRIX ===============
1846 //Count jets and trigger-jet pairs at MC generator level
1847 for(Int_t ij=0; ij<fListJetsGen->GetEntries(); ij++){
1848 AliAODJet* jet = (AliAODJet*)(fListJetsGen->At(ij));
1849 if(!jet) continue;
1850 Double_t etaJetGen = jet->Eta();
1851 Double_t ptJetGen = jet->Pt();
255cb71c 1852
8aa19603 1853 if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){
1854 fhJetPtGen->Fill(ptJetGen); // gen level pt spectum of jets response mx normalization
80ac66f6 1855
8aa19603 1856 Double_t areaJetGen = jet->EffectiveAreaCharged();
1857 Double_t ptUeFromCellMedianGen = rhoFromCellMedianGen*areaJetGen;
1858 Double_t ptUeConeGen = rhoConeGen*areaJetGen;
80ac66f6 1859
8aa19603 1860 fhJetPtSubUeMedianGen->Fill(ptJetGen - ptUeFromCellMedianGen);
1861 fhJetPtSubUeConeGen->Fill(ptJetGen - ptUeConeGen);
160dbdb8 1862 }
1863 }
8aa19603 1864 if(fListJets->GetEntries()>0 && fListJetsGen->GetEntries()>0){ //at least some reconstructed jets
1865
1866 Int_t ng = (Int_t) fListJetsGen->GetEntries();
1867 Int_t nr = (Int_t) fListJets->GetEntries();
160dbdb8 1868
8aa19603 1869 //Find closest MC generator - reconstructed jets
1870 if(faGenIndex.GetSize()<nr) faGenIndex.Set(nr); //idx of gen jet assoc to rec jet
1871 if(faRecIndex.GetSize()<ng) faRecIndex.Set(ng); //idx of rec jet assoc to gen jet
160dbdb8 1872
1873 if(fDebug){
8aa19603 1874 Printf("New Rec List %d gen index Array %d",nr,faGenIndex.GetSize());
1875 Printf("New Gen List %d rec index Array %d",ng,faRecIndex.GetSize());
160dbdb8 1876 }
1877 //matching of MC genrator level and reconstructed jets
8aa19603 1878 AliAnalysisHelperJetTasks::GetClosestJets(fListJetsGen,ng,fListJets,nr,faGenIndex,faRecIndex,fDebug);
1879
1880 // Fill response matrix
1881 for(Int_t ir = 0; ir < nr; ir++){
1882 AliAODJet *recJet = (AliAODJet*) fListJets->At(ir);
1883 Double_t etaJetRec = recJet->Eta();
1884 Double_t ptJetRec = recJet->Pt();
1885 Double_t areaJetRec = recJet->EffectiveAreaCharged();
160dbdb8 1886 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1887
8aa19603 1888 if((fJetEtaMin <= etaJetRec) && (etaJetRec <= fJetEtaMax)){
1889 Int_t ig = faGenIndex[ir]; //associated generator level jet
1890 if(ig >= 0 && ig < ng){
1891 if(fDebug > 10) Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1892 AliAODJet *genJet = (AliAODJet*) fListJetsGen->At(ig);
1893 Double_t ptJetGen = genJet->Pt();
1894 Double_t etaJetGen = genJet->Eta();
160dbdb8 1895
255cb71c 1896 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
8aa19603 1897 if((fJetEtaMin <= etaJetGen) && (etaJetGen <= fJetEtaMax)){
1898 fhJetPtGenVsJetPtRec->Fill(ptJetRec, ptJetGen);
f1e6baa6 1899 if(ptJetGen>0){
1900 fhJetPtResolutionVsPtGen->Fill(ptJetGen,(ptJetRec-ptJetGen)/ptJetGen);
1901 }
8aa19603 1902 Double_t areaJetGen = genJet->EffectiveAreaCharged();
1903 Double_t ptUeFromCellMedianGen = rhoFromCellMedianGen*areaJetGen;
1904 Double_t ptUeConeGen = rhoConeGen*areaJetGen;
1905 Double_t ptUeFromCellMedianRec = rhoFromCellMedian*areaJetRec;
1906 Double_t ptUeConeRec = rhoCone*areaJetRec;
1907 fhJetPtGenVsJetPtRecSubUeMedian->Fill(ptJetRec-ptUeFromCellMedianRec,
1908 ptJetGen-ptUeFromCellMedianGen);
1909 fhJetPtGenVsJetPtRecSubUeCone->Fill(ptJetRec-ptUeConeRec, ptJetGen-ptUeConeGen);
f1e6baa6 1910
1911 if((ptJetGen-ptUeConeGen)>0){
1912 fhJetPtResolutionVsPtConeGen->Fill(ptJetGen-ptUeConeGen,((ptJetRec-ptUeConeRec)- (ptJetGen-ptUeConeGen))/ (ptJetGen-ptUeConeGen));
1913 }
160dbdb8 1914 }
8aa19603 1915 }//ig>=0
160dbdb8 1916 }//rec jet in eta acceptance
1917 }//loop over reconstructed jets
1918 }// # of rec jets >0
8aa19603 1919
1920 //=========================== Full jet vs charged jet matrix ==========================
1921 if(fIsFullMC){
1922 //Count full jets and charged-jet pairs at MC generator level
1923 for(Int_t ij=0; ij<fListJetsGenFull->GetEntries(); ij++){
1924 AliAODJet* jetFull = (AliAODJet*)(fListJetsGenFull->At(ij));
1925 if(!jetFull) continue;
1926 Double_t etaJetFull = jetFull->Eta();
1927 Double_t ptJetFull = jetFull->Pt();
1928
1929 if((fJetEtaMin<=etaJetFull) && (etaJetFull<=fJetEtaMax)){
1930 fhJetPtGenFull->Fill(ptJetFull); // generator level pt spectum of full jets
1931 }
1932 }
1933 if(fListJetsGen->GetEntries()>0 && fListJetsGenFull->GetEntries()>0){ //at least some reconstructed jets
1934 Int_t nful = (Int_t) fListJetsGenFull->GetEntries();
1935 Int_t nchr = (Int_t) fListJetsGen->GetEntries();
1936
1937 //Find closest MC generator full - charged jet
1938 if(faGenIndex.GetSize()<nchr) faGenIndex.Set(nchr); //idx of gen FULL jet assoc to gen CHARGED jet
1939 if(faRecIndex.GetSize()<nful) faRecIndex.Set(nful); //idx of gen CHARGED jet assoc to gen FULL jet
1940
1941 if(fDebug){
1942 Printf("New Charg List %d Full index Array %d",nchr,faGenIndex.GetSize());
1943 Printf("New Full List %d Charg index Array %d",nful,faRecIndex.GetSize());
1944 }
1945 //matching of MC genrator level and reconstructed jets
1946 AliAnalysisHelperJetTasks::GetClosestJets(fListJetsGenFull,nful,fListJetsGen,nchr,faGenIndex,faRecIndex,fDebug);
1947
1948 // Fill response matrix
1949 for(Int_t ichr = 0; ichr < nchr; ichr++){ //charged jet loop
1950 AliAODJet *chJet = (AliAODJet*) fListJetsGen->At(ichr);
1951 Double_t etaJetCh = chJet->Eta();
1952 Double_t ptJetCh = chJet->Pt();
1953 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1954
1955 if((fJetEtaMin <= etaJetCh) && (etaJetCh <= fJetEtaMax)){
1956 Int_t iful = faGenIndex[ichr]; //associated generator level jet
1957 if(iful >= 0 && iful < nful){
1958 if(fDebug > 10) Printf("%s:%d iful = %d ichr = %d",(char*)__FILE__,__LINE__,iful,ichr);
1959 AliAODJet *genJetFull = (AliAODJet*) fListJetsGenFull->At(iful);
1960 Double_t ptJetFull = genJetFull->Pt();
1961 Double_t etaJetFull = genJetFull->Eta();
1962
1963 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1964 if((fJetEtaMin <= etaJetFull) && (etaJetFull <= fJetEtaMax)){
1965 fhJetPtGenChargVsJetPtGenFull->Fill(ptJetFull,ptJetCh);
1966 }
1967 }//iful>=0
1968 }//rec jet in eta acceptance
1969 }//loop over reconstructed jets
1970 }// # of rec jets >0
1971 }//pointer MC generator jets
1972 } //fill resp mx only for bin
1973 }//analyze generator level MC
1974
1975
49dc2855 1976 if(fIsKine){ //skip reconstructed data analysis in case of kine
1977 PostData(1, fOutputList);
1978 return;
1979 }
8e103a56 1980 //============= RECONSTRUCTED INCLUSIVE JETS ===============
ad869500 1981
ad869500 1982 Double_t etaJet = 0.0;
1983 Double_t pTJet = 0.0;
1984 Double_t areaJet = 0.0;
1985 Double_t phiJet = 0.0;
9f10e189 1986 //Int_t indexLeadingJet = -1;
1987 //Double_t pTLeadingJet = -10.0;
1988 //Double_t areaLeadingJet = -10.0;
8e103a56 1989
ad869500 1990 for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
1991 AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
1992 if(!jet) continue;
1993 etaJet = jet->Eta();
1994 phiJet = jet->Phi();
1995 pTJet = jet->Pt();
1996 if(pTJet==0) continue;
1997
1998 if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
255cb71c 1999 /*areaJet = jet->EffectiveAreaCharged();*/
7fc3d134 2000
ad869500 2001 //Jet Diagnostics---------------------------------
8e103a56 2002 fhJetPhi->Fill(pTJet, RelativePhi(phiJet,0.0)); //phi -pi,pi
2003 fhJetEta->Fill(pTJet, etaJet);
7fc3d134 2004 //search for leading jet
255cb71c 2005 /*if(pTJet > pTLeadingJet){
7fc3d134 2006 indexLeadingJet = ij;
2007 pTLeadingJet = pTJet;
7fc3d134 2008 areaLeadingJet = areaJet;
2009 }
2010
2011 // raw spectra of INCLUSIVE jets
80ac66f6 2012 //Centrality, pTjet, A
7fc3d134 2013 Double_t fillraw[] = { centValue,
2014 pTJet,
7fc3d134 2015 areaJet
2016 };
255cb71c 2017 fHJetPtRaw->Fill(fillraw);*/
8e103a56 2018 }
255cb71c 2019 /*
8e103a56 2020 if(indexLeadingJet > -1){
2021 // raw spectra of LEADING jets
2022 //Centrality, pTjet, A
2023 Double_t fillleading[] = { centValue,
2024 pTLeadingJet,
2025 areaLeadingJet
2026 };
2027 fHLeadingJetPtRaw->Fill(fillleading);
2028 }
255cb71c 2029 */
364e98aa 2030 //================ TWO PARTICLE CORRELATION CORRECTION ===============
2031 if(fIsChargedMC && fHardest==0 ){
2032 if(indexTrigg>=0){ // Reconstructed trigger
2033 AliVParticle *triggerHadron = (AliVParticle*) particleList.At(indexTrigg);
2034 Int_t genTriggIndex = -1;
2035 if(triggerHadron){
2036 Int_t trigLabel = TMath::Abs(triggerHadron->GetLabel());
2037
2038 Int_t nGen = particleListGen.GetEntries();
2039 for(Int_t ig=0; ig<nGen; ig++){
2040 AliVParticle *trkGen = (AliVParticle*) particleListGen.At(ig);
2041 if(!trkGen) continue;
2042
2043 Int_t genLabel = TMath::Abs(trkGen->GetLabel());
2044 if(trigLabel==genLabel){
2045 if(fTriggerPtRangeLow <= trkGen->Pt() &&
2046 trkGen->Pt() < fTriggerPtRangeHigh &&
2047 TMath::Abs((Float_t) trkGen->Eta()) < fTriggerEtaCut){
2048 genTriggIndex = ig; //There is a gen level particle track matching to rec trigger
2049 break;
2050 }
2051 }
2052 }
2053 if(genTriggIndex>-1){ // make generator level TT x Assoc
2054 AliVParticle *trkGenTT = (AliVParticle*) particleListGen.At(genTriggIndex);
2055 if(trkGenTT){
2056 fhTriggerCounterGenLevel->Fill(0.5); // count gen triggers that have reconstructed analog
2057 Double_t deltaPhi,deltaEta,deltaR;
2058
2059 for(Int_t ig=0; ig<nGen; ig++){
2060 if(ig==genTriggIndex) continue; //skip trigger
2061 AliVParticle *trkGenA = (AliVParticle*) particleListGen.At(ig);
2062 if(!trkGenA) continue;
2063 if(trkGenA->Pt()<15.0) continue;
2064 if(trkGenA->Pt()>50.0) continue;
2065 deltaPhi = RelativePhi(trkGenTT->Phi(), trkGenA->Phi());
2066 deltaEta = trkGenTT->Eta() - trkGenA->Eta();
2067 deltaR = sqrt(deltaPhi*deltaPhi + deltaEta*deltaEta);
2068
2069 if(trkGenA->Pt()<20.0){ // gen level assoc particles
2070 fhDeltaRMultTriggersGenLevelLow->Fill(deltaR);
2071 fhDeltaPhiMultTriggersGenLevelLow->Fill(deltaPhi);
2072 }else{
2073 fhDeltaRMultTriggersGenLevelHigh->Fill(deltaR);
2074 fhDeltaPhiMultTriggersGenLevelHigh->Fill(deltaPhi);
2075 }
2076 }
2077 }
2078 }
2079 }
2080 }
2081 }
8e103a56 2082
2083 // =============== RECONSTRUCTED TRIGGER-JET PAIRS ================
8aa19603 2084 if(fIsChargedMC && fFillRespMx){
2085 FillEffHistos(&particleList, &particleListGen); //Fill efficiency histos
2086 }
255cb71c 2087 Bool_t filledOnce = kTRUE; //fill rho histogram only once per event
8e103a56 2088 //set ranges of the trigger loop
2089 Int_t ilow = (fHardest==0 || fHardest==1) ? indexTrigg : 0;
2090 Int_t ihigh = (fHardest==0 || fHardest==1) ? indexTrigg+1 : particleList.GetEntries();
7fc3d134 2091
8e103a56 2092 for(Int_t itrk= ilow; itrk< ihigh; itrk++){ //loop over possible trigger
2093 indexTrigg = itrk; //trigger hadron with pT >10 GeV
7fc3d134 2094
8e103a56 2095 if(indexTrigg < 0) continue;
2096
2097 AliVParticle *triggerHadron = (AliVParticle*) particleList.At(indexTrigg);
2098 if(!triggerHadron){
2099 PostData(1, fOutputList);
2100 return;
2101 }
2102 if(fHardest >= 2){
7d9fd386 2103 if(triggerHadron->Pt() < 6.0) continue; //all hadrons pt>6
8e103a56 2104 }
2105 if(TMath::Abs((Float_t) triggerHadron->Eta())> fTriggerEtaCut) continue;
7fc3d134 2106
8e103a56 2107 //Fill trigger histograms
2108 fh2Ntriggers->Fill(centValue,triggerHadron->Pt()); //trigger pT
2109 fhTriggerPhi->Fill(triggerHadron->Pt(),RelativePhi(triggerHadron->Phi(),0.0)); //phi -pi,pi
2110 fhTriggerEta->Fill(triggerHadron->Pt(),triggerHadron->Eta());
7fc3d134 2111
2112
8e103a56 2113 //---------- make trigger-jet pairs ---------
255cb71c 2114 //Int_t injet4 = 0;
2115 //Int_t injet = 0;
255cb71c 2116
8e103a56 2117 for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
2118 AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
2119 if(!jet) continue;
2120 etaJet = jet->Eta();
2121 phiJet = jet->Phi();
2122 pTJet = jet->Pt();
2123 if(pTJet==0) continue;
2124
2125 if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
2126 areaJet = jet->EffectiveAreaCharged();
49dc2855 2127
2128 //subtract bg using estinates base on median of kT jets
2129 Double_t ptUeFromCellMedian = rhoFromCellMedian*areaJet;
2130 Double_t ptUeCone = rhoCone*areaJet;
2131
255cb71c 2132 //if(areaJet >= 0.07) injet++;
2133 //if(areaJet >= 0.4) injet4++;
8e103a56 2134
2135 Double_t dphi = RelativePhi(triggerHadron->Phi(), phiJet);
2136 fhDphiTriggerJet->Fill(dphi); //Input
2137
2138 //Dphi versus jet pT
2139 //Centrality, Dphi=phiTrig-phiJet, pTjet, pTtrigg
255cb71c 2140 /*Double_t filldp[] = { centValue,
8e103a56 2141 dphi,
2142 pTJet,
2143 triggerHadron->Pt()
2144 };
2145 fHDphiVsJetPtAll->Fill(filldp);
255cb71c 2146 */
49dc2855 2147 //Rongrong's analysis
2148
2149 Double_t dPhi = phiJet - triggerHadron->Phi();
2150 if(dPhi>2*TMath::Pi()) dPhi -= 2*TMath::Pi();
2151 if(dPhi<-2*TMath::Pi()) dPhi += 2*TMath::Pi();
2152 if(dPhi<-0.5*TMath::Pi()) dPhi += 2*TMath::Pi();
2153 if(dPhi>1.5*TMath::Pi()) dPhi -= 2*TMath::Pi();
2154
2155 //[A,pTjet,pTtrig,dphi]
2156 tmpArrayFour[0] = areaJet;
2157 tmpArrayFour[1] = pTJet;
2158 tmpArrayFour[2] = triggerHadron->Pt();
2159 tmpArrayFour[3] = dPhi;
2160
2161 fHJetPhiCorr->Fill(tmpArrayFour); // Dphi distribution jet-triger
2162
2163 //[A,pTjet-pTUe,pTtrig,dphi]
2164 tmpArrayFour[1] = pTJet - ptUeFromCellMedian;
2165 fHJetPhiCorrSubUeMedian->Fill(tmpArrayFour);
2166
2167 //[A,pTjet-pTUe,pTtrig,dphi]
2168 tmpArrayFour[1] = pTJet - ptUeCone;
2169 fHJetPhiCorrSubUeCone->Fill(tmpArrayFour);
2170
2171 //Leticia's analysis
2172
8e103a56 2173 // Select back to back trigger - jet pairs
2174 if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue;
2175 fhDphiTriggerJetAccept->Fill(dphi); //Accepted
255cb71c 2176
2177 //NO bg subtraction
2178 //Centrality, A, pTjet, pTtrigg, dphi
2179 tmpArrayFive[0] = centValue;
2180 tmpArrayFive[1] = areaJet;
2181 tmpArrayFive[2] = pTJet;
2182 tmpArrayFive[3] = triggerHadron->Pt();
2183 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
2184 fHJetSpec->Fill(tmpArrayFive);
2185
49dc2855 2186 //Centrality, A, pTjet-pTbgCellMedian, pTtrigg, dphi
255cb71c 2187 tmpArrayFive[0] = centValue;
2188 tmpArrayFive[1] = areaJet;
8aa19603 2189 tmpArrayFive[2] = pTJet-ptUeFromCellMedian;
255cb71c 2190 tmpArrayFive[3] = triggerHadron->Pt();
2191 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
2192 fHJetSpecSubUeMedian->Fill(tmpArrayFive);
2193
8aa19603 2194 //Centrality, A, pTjet-pTbgCone, pTtrigg, dphi
255cb71c 2195 tmpArrayFive[0] = centValue;
2196 tmpArrayFive[1] = areaJet;
8aa19603 2197 tmpArrayFive[2] = pTJet - ptUeCone;
255cb71c 2198 tmpArrayFive[3] = triggerHadron->Pt();
2199 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
8aa19603 2200 fHJetSpecSubUeCone->Fill(tmpArrayFive);
255cb71c 2201
2202 //Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", kT median
2203 tmpArrayFive[0] = areaJet;
2204 tmpArrayFive[1] = pTJet;
8aa19603 2205 tmpArrayFive[2] = pTJet - ptUeFromCellMedian;
2206 tmpArrayFive[3] = ptUeFromCellMedian;
2207 tmpArrayFive[4] = rhoFromCellMedian;
255cb71c 2208 fHJetUeMedian->Fill(tmpArrayFive);
8e103a56 2209
8aa19603 2210 //Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", Cone median
255cb71c 2211 tmpArrayFive[0] = areaJet;
2212 tmpArrayFive[1] = pTJet;
8aa19603 2213 tmpArrayFive[2] = pTJet - ptUeCone;
2214 tmpArrayFive[3] = ptUeCone;
2215 tmpArrayFive[4] = rhoCone;
2216 fHJetUeCone->Fill(tmpArrayFive);
255cb71c 2217
2218 if(filledOnce){ //fill for each event only once
8aa19603 2219 Double_t fillRho[] = { rhoCone,rhoFromCellMedian};
2220 fHRhoUeMedianVsCone->Fill(fillRho);
255cb71c 2221 filledOnce = kFALSE;
2222 }
8e103a56 2223 }//jet loop
2224
2225 //Fill Jet Density In the Event A>0.07
255cb71c 2226 /*if(injet>0){
8e103a56 2227 Double_t filldens[]={ (Double_t) particleList.GetEntries(),
ad869500 2228 injet/fkAcceptance,
2229 triggerHadron->Pt()
2230 };
8e103a56 2231 fHJetDensity->Fill(filldens);
2232 }
ad869500 2233
8e103a56 2234 //Fill Jet Density In the Event A>0.4
2235 if(injet4>0){
2236 Double_t filldens4[]={ (Double_t) particleList.GetEntries(),
2237 injet4/fkAcceptance,
2238 triggerHadron->Pt()
2239 };
2240 fHJetDensityA4->Fill(filldens4);
255cb71c 2241 }*/
ad869500 2242 }
2243
8e103a56 2244
ad869500 2245 PostData(1, fOutputList);
2246}
2247
2248//----------------------------------------------------------------------------
2249void AliAnalysisTaskJetCorePP::Terminate(const Option_t *)
2250{
2251 // Draw result to the screen
2252 // Called once at the end of the query
2253
2254 if(fDebug) printf("AliAnalysisTaskJetCorePP DONE\n");
2255 if(!GetOutputData(1)) return;
2256}
2257
2258
2259//----------------------------------------------------------------------------
2260Int_t AliAnalysisTaskJetCorePP::GetListOfTracks(TList *list){
2261 //Fill the list of accepted tracks (passed track cut)
2262 //return consecutive index of the hardest ch hadron in the list
2263 Int_t iCount = 0;
2264 AliAODEvent *aodevt = NULL;
2265
2266 if(!fESD) aodevt = fAODIn;
2267 else aodevt = fAODOut;
2268
2269 if(!aodevt) return -1;
2270
2271 Int_t index = -1; //index of the highest particle in the list
2272 Double_t ptmax = -10;
8e103a56 2273 Int_t triggers[200];
2274 Int_t ntriggers = 0; //index in triggers array
1e3ad34d 2275
ad869500 2276
80ac66f6 2277 for(Int_t it = 0; it < aodevt->GetNumberOfTracks(); it++){
0fe881ea 2278 //AliAODTrack *tr = aodevt->GetTrack(it);
f15c1f69 2279 AliAODTrack *tr = dynamic_cast<AliAODTrack*>(aodevt->GetTrack(it));
2280 if(!tr) AliFatal("Not a standard AOD");
0fe881ea 2281
d1405a52
ML		 2282 if((fFilterMask > 0) && !(tr->TestFilterBit(fFilterMask))) continue;
2283 //if((fFilterMask > 0) && !(tr->IsHybridGlobalConstrainedGlobal())) continue;
ad869500 2284 if(TMath::Abs((Float_t) tr->Eta()) > fTrackEtaCut) continue;
2285 if(tr->Pt() < fTrackLowPtCut) continue;
2286 list->Add(tr);
8e103a56 2287
2288 //Search trigger:
2289 if(fHardest>0){ //leading particle
f1e6baa6 2290 if(tr->Pt()>ptmax &&
2291 TMath::Abs((Float_t) tr->Eta()) < fTriggerEtaCut){
8e103a56 2292 ptmax = tr->Pt();
2293 index = iCount;
2294 }
2295 }
2296
2297 if(fHardest==0 && ntriggers<200){ //single inclusive
2298 if(fTriggerPtRangeLow <= tr->Pt() &&
f1e6baa6 2299 tr->Pt() < fTriggerPtRangeHigh &&
2300 TMath::Abs((Float_t) tr->Eta()) < fTriggerEtaCut){
8e103a56 2301 triggers[ntriggers] = iCount;
2302 ntriggers++;
2303 }
ad869500 2304 }
8e103a56 2305
ad869500 2306 iCount++;
2307 }
f1e6baa6 2308 //inclusive trigger
2309 if(ntriggers>0){ //select random inclusive trigger
2310 fhInclTrigCounter->Fill(0.5,ntriggers); //count inclusive triggers
2311
2312 //Lower pTa bin
2313 for(Int_t it1=0; it1<ntriggers; it1++){
2314 AliVParticle* tGent1 = (AliVParticle*) list->At(triggers[it1]);
2315 if(!tGent1) continue;
2316 for(Int_t ia=0; ia<list->GetEntries(); ia++){
2317 if(ia == triggers[it1]) continue;
2318 AliVParticle* tGent2 = (AliVParticle*) list->At(ia);
2319 if(!tGent2) continue;
2320 if(tGent2->Pt()<15.0) continue;
2321 if(tGent2->Pt()>20.0) continue;
2322 fhDeltaPhiMultTriggersInclLow->Fill(RelativePhi(tGent1->Phi(),tGent2->Phi()));
2323 }
2324 }
2325 //Higher pTa bin
2326 for(Int_t it1=0; it1<ntriggers-1; it1++){
2327 AliVParticle* tGent1 = (AliVParticle*) list->At(triggers[it1]);
2328 if(!tGent1) continue;
2329
2330 for(Int_t it2=it1+1; it2<ntriggers; it2++){
2331 AliVParticle* tGent2 = (AliVParticle*) list->At(triggers[it2]);
2332 if(!tGent2) continue;
2333 fhDeltaPhiMultTriggersInclHigh->Fill(RelativePhi(tGent1->Phi(),tGent2->Phi()));
2334 }
2335 }
2336 }
ad869500 2337
f1e6baa6 2338 //single inclusive trigger
8e103a56 2339 if(fHardest==0 && ntriggers>0){ //select random inclusive trigger
2340 Int_t rnd = fRandom->Integer(ntriggers); //0 to ntriggers-1
2341 index = triggers[rnd];
1e3ad34d 2342
2343 fhNofMultipleTriggers->Fill(ntriggers-1);
ac4f24f7 2344
2345 Double_t deltaPhi, deltaEta, deltaR;
30557bf5 2346 Int_t iLow=0;
2347 Int_t iHigh=0;
ac4f24f7 2348 //Correlation with single inclusive trigger
2349 AliVParticle* tGent1 = (AliVParticle*) list->At(index);
2350 if(tGent1){
58003b11 2351 //for(Int_t ia=0; ia<ntriggers; ia++)
2352 for(Int_t ia=0; ia<list->GetEntries(); ia++){
2353 //if(triggers[ia]==index) continue;
2354 if(ia==index) continue;
2355 //AliVParticle* tGent2 = (AliVParticle*) list->At(triggers[ia]);
2356 AliVParticle* tGent2 = (AliVParticle*) list->At(ia);
ac4f24f7 2357 if(!tGent2) continue;
58003b11 2358 if(tGent2->Pt()<15.0) continue;
30557bf5 2359 if(tGent2->Pt()>50.0) continue;
ac4f24f7 2360 deltaPhi = RelativePhi(tGent1->Phi(),tGent2->Phi());
2361 deltaEta = tGent1->Eta()-tGent2->Eta();
2362 deltaR = sqrt(deltaPhi*deltaPhi + deltaEta*deltaEta);
30557bf5 2363
f1e6baa6 2364 if(tGent2->Pt()<20.0){
364e98aa 2365 fhDeltaRMultTriggersLow->Fill(deltaR);
f1e6baa6 2366 fhDeltaPhiMultTriggersLow->Fill(deltaPhi);
2367 }else{
364e98aa 2368 fhDeltaRMultTriggersHigh->Fill(deltaR);
f1e6baa6 2369 fhDeltaPhiMultTriggersHigh->Fill(deltaPhi);
2370 }
2371
30557bf5 2372 if(deltaR<0.4){
2373 if(tGent2->Pt()<20.0) iLow++;
2374 else iHigh++;
1e3ad34d 2375 }
30557bf5 2376 }
1e3ad34d 2377 }
30557bf5 2378 fhNofMultipleTriggersCone->Fill(iLow+iHigh);
2379 fhNofMultipleTriggersConeLow->Fill(iLow);
2380 fhNofMultipleTriggersConeHigh->Fill(iHigh);
2381
ad869500 2382 }
8e103a56 2383
2384 return index;
ad869500 2385}
2386
2387//----------------------------------------------------------------------------
80ac66f6 2388/*
ad869500 2389Double_t AliAnalysisTaskJetCorePP::GetBackgroundInPerpCone(Float_t jetR, Double_t jetPhi, Double_t jetEta, TList* trkList){
2390 //calculate sum of track pT in the cone perpendicular in phi to the jet
2391 //jetR = cone radius
2392 // jetPhi, jetEta = direction of the jet
2393 Int_t numberOfTrks = trkList->GetEntries();
2394 Double_t pTsum = 0.0;
2395 Double_t perpConePhi = jetPhi + TMath::Pi()/2;//perp cone w.r.t. jet in phi
2396 for(Int_t it=0; it<numberOfTrks; it++){
2397 AliVParticle *trk = (AliVParticle*) trkList->At(it);
2398 Double_t dphi = RelativePhi(perpConePhi,trk->Phi());
2399 Double_t deta = trk->Eta()-jetEta;
2400
2401 if( (dphi*dphi + deta*deta)< (jetR*jetR)){
2402 pTsum += trk->Pt();
2403 }
2404 }
2405
2406 return pTsum;
2407}
80ac66f6 2408*/
ad869500 2409//----------------------------------------------------------------------------
2410
2411Double_t AliAnalysisTaskJetCorePP::RelativePhi(Double_t mphi,Double_t vphi){
2412 //Get relative azimuthal angle of two particles -pi to pi
2413 if (vphi < -TMath::Pi()) vphi += TMath::TwoPi();
2414 else if (vphi > TMath::Pi()) vphi -= TMath::TwoPi();
2415
2416 if (mphi < -TMath::Pi()) mphi += TMath::TwoPi();
2417 else if (mphi > TMath::Pi()) mphi -= TMath::TwoPi();
2418
2419 Double_t dphi = mphi - vphi;
2420 if (dphi < -TMath::Pi()) dphi += TMath::TwoPi();
2421 else if (dphi > TMath::Pi()) dphi -= TMath::TwoPi();
2422
2423 return dphi;//dphi in [-Pi, Pi]
2424}
2425
2426
80ac66f6 2427//----------------------------------------------------------------------------
2428Bool_t AliAnalysisTaskJetCorePP::SelectMCGenTracks(AliVParticle *trk, TList *trkList, Double_t &ptLeading, Int_t &index, Int_t counter){
2429 //fill track efficiency denominator
2430 if(trk->Pt() < fTrackLowPtCut) return kFALSE;
2431 if(trk->Charge()==0) return kFALSE;
2432 if(TMath::Abs((Float_t) trk->Eta()) > fTrackEtaCut) return kFALSE;
2433 trkList->Add(trk);
8aa19603 2434 if(fFillRespMx) fhPtTrkTruePrimGen->Fill(trk->Pt(),trk->Eta()); //Efficiency denominator
8e103a56 2435
2436 //Search trigger:
2437 if(fHardest>0){ //leading particle
f1e6baa6 2438 if(ptLeading < trk->Pt() &&
2439 TMath::Abs((Float_t) trk->Eta()) < fTriggerEtaCut){
8e103a56 2440 index = counter;
2441 ptLeading = trk->Pt();
2442 }
80ac66f6 2443 }
8e103a56 2444
2445 if(fHardest==0){ //single inclusive
2446 index = -1;
2447 if(fTriggerPtRangeLow <= trk->Pt() &&
f1e6baa6 2448 trk->Pt() < fTriggerPtRangeHigh &&
2449 TMath::Abs((Float_t) trk->Eta()) < fTriggerEtaCut){
8e103a56 2450 index = counter;
2451 }
2452 }
2453
80ac66f6 2454 return kTRUE;
2455}
2456
2457//----------------------------------------------------------------------------
2458void AliAnalysisTaskJetCorePP::FillEffHistos(TList *recList, TList *genList){
2459 //fill track efficiency numerator
2460 if(!recList) return;
2461 if(!genList) return;
2462 Int_t nRec = recList->GetEntries();
2463 Int_t nGen = genList->GetEntries();
2464 for(Int_t ir=0; ir<nRec; ir++){
2465 AliAODTrack *trkRec = (AliAODTrack*) recList->At(ir);
2466 if(!trkRec) continue;
2467 Int_t recLabel = TMath::Abs(trkRec->GetLabel());
2468 Bool_t matched = kFALSE;
2469
2470 for(Int_t ig=0; ig<nGen; ig++){
2471 AliVParticle *trkGen = (AliVParticle*) genList->At(ig);
2472 if(!trkGen) continue;
2473 Int_t genLabel = TMath::Abs(trkGen->GetLabel());
2474 if(recLabel==genLabel){
2475 fhPtTrkTruePrimRec->Fill(trkGen->Pt(),trkGen->Eta());
2476 matched = kTRUE;
2477 break;
2478 }
2479 }
2480 if(!matched) fhPtTrkSecOrFakeRec->Fill(trkRec->Pt(),trkRec->Eta());
2481 }
2482
2483 return;
2484}
255cb71c 2485//________________________________________________________________
d1405a52 2486void AliAnalysisTaskJetCorePP::EstimateBgRhoMedian(TList *listJet, TList* listPart, Double_t &rhoMedian, Int_t mode){
8aa19603 2487 //Estimate background rho by means of integrating track pT outside identified jet cones
2488
2489 rhoMedian = 0;
2490
2491 //identify leading jet in the full track acceptance
2492 Int_t idxLeadingJet = -1;
2493 Double_t pTleading = 0.0;
2494 AliAODJet* jet = NULL;
2495
2496 for(Int_t ij = 0; ij < listJet->GetEntries(); ij++){ //loop over bg kT jets
2497 jet = (AliAODJet*)(listJet->At(ij));
2498 if(!jet) continue;
2499 if(TMath::Abs(jet->Eta()) > fTrackEtaCut) continue; //track acceptance cut
2500 if(jet->Pt() > pTleading){
2501 idxLeadingJet = ij;
2502 pTleading = jet->Pt();
2503 }
2504 }
2505
2506 Int_t njets = 0;
2507 static Double_t jphi[100];
2508 static Double_t jeta[100];
2509 static Double_t jRsquared[100];
2510
2511 for(Int_t ij = 0; ij< listJet->GetEntries(); ij++){ //loop over bg kT jets
2512
2513 jet = (AliAODJet*)(listJet->At(ij));
2514 if(!jet) continue;
2515 if(TMath::Abs(jet->Eta()) > fTrackEtaCut) continue;
2516
2517 if((jet->Pt() > fBgMaxJetPt || ij== idxLeadingJet) && njets<100){
2518 jphi[njets] = RelativePhi(jet->Phi(),0.0); //-pi,pi
2519 jeta[njets] = jet->Eta();
2520 jRsquared[njets] = fSafetyMargin*jet->EffectiveAreaCharged()/TMath::Pi(); //scale jet radius
2521 njets++;
255cb71c 2522 }
2523 }
8aa19603 2524
255cb71c 2525
8aa19603 2526 static Double_t nOutCone[10][4];
2527 static Double_t sumPtOutOfCone[10][4];
2528
2529 for(Int_t ie=0; ie<fnEta; ie++){
2530 for(Int_t ip=0; ip<fnPhi; ip++){
2531 nOutCone[ip][ie] = 0.0; //initialize counter
2532 sumPtOutOfCone[ip][ie] = 0.0;
255cb71c 2533 }
8aa19603 2534 }
255cb71c 2535
8aa19603 2536 Double_t rndphi, rndeta;
2537 Double_t rndphishift, rndetashift;
2538 Double_t dphi, deta;
2539 Bool_t bIsInCone;
2540
2541 for(Int_t it=0; it<fnTrials; it++){
2542
2543 rndphi = fRandom->Uniform(0, fPhiSize);
2544 rndeta = fRandom->Uniform(0, fEtaSize);
2545
2546 for(Int_t ip=0; ip<fnPhi; ip++){ //move radom position to each cell
2547 rndphishift = rndphi + ip*fPhiSize - TMath::Pi();
2548 for(Int_t ie=0; ie<fnEta; ie++){
2549 rndetashift = rndeta + ie*fEtaSize - fEtaSize;
2550
2551 bIsInCone = 0; //tag if trial is in the jet cone
2552 for(Int_t ij=0; ij<njets; ij++){
2553 deta = jeta[ij] - rndetashift;
2554 dphi = RelativePhi(rndphishift,jphi[ij]);
2555 if((dphi*dphi + deta*deta) < jRsquared[ij]){
2556 bIsInCone = 1;
2557 break;
2558 }
2559 }
2560 if(!bIsInCone) nOutCone[ip][ie]++;
2561 }
2562 }
2563 }
2564
255cb71c 2565
8aa19603 2566 //Sum particle pT outside jets in cells
2567 Int_t npart = listPart->GetEntries();
2568 Int_t phicell,etacell;
2569 AliVParticle *part = NULL;
2570 for(Int_t ip=0; ip < npart; ip++){
2571
2572 part = (AliVParticle*) listPart->At(ip);
2573 if(!part) continue;
2574 if( TMath::Abs(part->Eta()) > fEtaSize) continue; //skip part outside +-0.7 in eta
2575
2576 bIsInCone = 0; //init
2577 for(Int_t ij=0; ij<njets; ij++){
2578 dphi = RelativePhi(jphi[ij], part->Phi());
2579 deta = jeta[ij] - part->Eta();
2580 if((dphi*dphi + deta*deta) < jRsquared[ij]){
2581 bIsInCone = 1;
2582 break;
2583 }
255cb71c 2584 }
8aa19603 2585 if(!bIsInCone){
2586 phicell = TMath::Nint(TMath::Floor((RelativePhi(part->Phi(),0.0) + TMath::Pi())/fPhiSize));
2587 etacell = TMath::Nint(TMath::Floor((part->Eta()+fEtaSize)/fEtaSize));
2588 sumPtOutOfCone[phicell][etacell]+= part->Pt();
2589 }
2590 }
2591
2592 static Double_t rhoInCells[20];
2593 Double_t relativeArea;
2594 Int_t nCells=0;
d1405a52 2595 Double_t bufferArea=0.0, bufferPt=0.0; //sum cells where A< fJetFreeAreaFrac
8aa19603 2596 for(Int_t ip=0; ip<fnPhi; ip++){
2597 for(Int_t ie=0; ie<fnEta; ie++){
2598 relativeArea = nOutCone[ip][ie]/fnTrials;
d1405a52
ML		 2599 //cout<<"RA= "<< relativeArea<<" BA="<<bufferArea<<" BPT="<< bufferPt <<endl;
2600
2601 bufferArea += relativeArea;
2602 bufferPt += sumPtOutOfCone[ip][ie];
2603 if(bufferArea > fJetFreeAreaFrac){
2604 rhoInCells[nCells] = bufferPt/(bufferArea*fCellArea);
2605
2606 if(mode==1) fhCellAreaToMedianGen->Fill(bufferArea*fCellArea);
2607 else fhCellAreaToMedian->Fill(bufferArea*fCellArea);
2608
2609 bufferArea = 0.0;
2610 bufferPt = 0.0;
2611 nCells++;
2612 }
8aa19603 2613 }
255cb71c 2614 }
8aa19603 2615
d1405a52
ML		 2616
2617 if(nCells>0){
2618 rhoMedian = TMath::Median(nCells, rhoInCells);
2619 if(mode==1){ //mc data
2620 fhEntriesToMedianGen->Fill(nCells);
2621 }else{ //real data
2622 fhEntriesToMedian->Fill(nCells);
2623 }
2624 }
8aa19603 2625
255cb71c 2626}
80ac66f6 2627
255cb71c 2628//__________________________________________________________________
8aa19603 2629void AliAnalysisTaskJetCorePP::EstimateBgCone(TList *listJet, TList* listPart, Double_t &rhoPerpCone){
2630
2631 rhoPerpCone = 0.0; //init
2632
2633 if(!listJet) return;
2634 Int_t njet = listJet->GetEntries();
2635 Int_t npart = listPart->GetEntries();
2636 Double_t pTleading = 0.0;
2637 Double_t phiLeading = 1000.;
2638 Double_t etaLeading = 1000.;
2639
2640 for(Int_t jsig=0; jsig < njet; jsig++){
2641 AliAODJet* signaljet = (AliAODJet*)(listJet->At(jsig));
2642 if(!signaljet) continue;
2643 if((signaljet->Eta()<fJetEtaMin) && (fJetEtaMax<signaljet->Eta())) continue; //acceptance cut
2644 if(signaljet->Pt() >= pTleading){ //replace leading and subleading jet
2645 pTleading = signaljet->Pt();
2646 phiLeading = signaljet->Phi();
2647 etaLeading = signaljet->Eta();
2648 }
2649 }
2650
2651 Double_t phileftcone = phiLeading + TMath::Pi()/2;
2652 Double_t phirightcone = phiLeading - TMath::Pi()/2;
2653 Double_t dp, de;
2654
2655 for(Int_t ip=0; ip < npart; ip++){
2656
2657 AliVParticle *part = (AliVParticle*) listPart->At(ip);
2658 if(!part){
2659 continue;
2660 }
2661
2662 dp = RelativePhi(phileftcone, part->Phi());
2663 de = etaLeading - part->Eta();
2664 if( sqrt(dp*dp + de*de)< fBgConeR) rhoPerpCone += part->Pt();
2665
2666 dp = RelativePhi(phirightcone, part->Phi());
2667 if( sqrt(dp*dp + de*de)< fBgConeR) rhoPerpCone += part->Pt();
2668
2669 }
2670
2671 //normalize total pT by two times cone are
2672 rhoPerpCone = rhoPerpCone/(2*TMath::Pi()*fBgConeR*fBgConeR);
2673
2674}
2675//__________________________________________________________________
80ac66f6 2676
255cb71c 2677void AliAnalysisTaskJetCorePP::ReadTClonesArray(TString bname, TList *list){
2678 //Convert TClones array of jets to TList
2679
2680 if(!strlen(bname.Data())){
2681 AliError(Form("Jet branch %s not set.", bname.Data()));
2682 return;
2683 }
2684
2685 TClonesArray *array=0x0;
49dc2855 2686 if(fUseExchContainer){ //take input from exchange containers
2687 if(fJetBranchNameKine.Length()>0 && bname.CompareTo(fJetBranchNameKine.Data())==0){
2688 array = dynamic_cast<TClonesArray*>(GetInputData(1)); //connect exchange container slot 1
2689 }
3d1fac8a 2690 if(fJetBranchNameBgKine.Length()>0 && bname.CompareTo(fJetBranchNameBgKine.Data())==0){
2691 array = dynamic_cast<TClonesArray*>(GetInputData(2)); //connect exchange container slot 2
2692 }
49dc2855 2693 }else{ //take input from AOD
2694 if(fAODOut&&!array){
2695 array = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(bname.Data()));
2696 }
2697 if(fAODExtension&&!array){
2698 array = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(bname.Data()));
2699 }
2700 if(fAODIn&&!array){
2701 array = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(bname.Data()));
2702 }
255cb71c 2703 }
2704
2705 list->Clear(); //clear list beforehand
2706
2707 if(array){
2708 if(fDebug){
2709 Printf("########## %s: %d jets",bname.Data(), array->GetEntriesFast());
2710 }
49dc2855 2711
255cb71c 2712 for(Int_t iJet = 0; iJet < array->GetEntriesFast(); iJet++) {
2713 AliAODJet *jet = dynamic_cast<AliAODJet*>((*array)[iJet]);
2714 if (jet) list->Add(jet);
2715 }
2716 }
2717
2718 return;
2719}