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