]> 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
fix AliHLTGlobalEsdConverterComponent
[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){
58003b11 1330 //for(Int_t ia=0; ia<ntriggersMCa; ia++)
1331 for(Int_t ia=0; ia< npar; ia++){
1332 //if(indexTriggGena == triggersMCa[ia]) continue;
1333 if(indexTriggGena == ia) continue;
6c1489e1 1334
58003b11 1335 //AliVParticle* tGenTz = (AliVParticle*) particleListGen.At(triggersMCa[ia]);
1336 AliVParticle* tGenTz = (AliVParticle*) particleListGen.At(ia);
6c1489e1 1337 if(!tGenTz) continue;
58003b11 1338 if(tGenTz->Pt()*0.9<15.0) continue;
1339 if(tGenTz->Pt()*0.9>100.0) continue;
6c1489e1 1340
1341 deltaPhia = RelativePhi(tGenTa->Phi(),tGenTz->Phi());
1342 deltaEtaa = tGenTa->Eta()-tGenTz->Eta();
1343 deltaRa = sqrt(deltaPhia*deltaPhia + deltaEtaa*deltaEtaa);
1344
1345 if(deltaRa<0.4) aa++;
1346 }
1347 }
1348 fhNofMultipleTriggersConeGenA->Fill(aa);
1349 }
1350
1351 if(ntriggersMCb>0){
1352 Int_t rndb = fRandom->Integer(ntriggersMCb); //0 to ntriggers-1
1353 Int_t indexTriggGenb = triggersMCb[rndb];
1354
1355 Double_t deltaPhib, deltaEtab, deltaRb;
1356 Int_t bb = 0;
1357
1358 //Correlation with single inclusive TRIGGER
1359 AliVParticle* tGenTb = (AliVParticle*) particleListGen.At(indexTriggGenb);
1360 if(tGenTb){
58003b11 1361 //for(Int_t ib=0; ib<ntriggersMCb; ib++)
1362 for(Int_t ib=0; ib<npar; ib++){
1363 //if(indexTriggGenb == triggersMCb[ib]) continue;
1364 if(indexTriggGenb == ib) continue;
6c1489e1 1365
58003b11 1366 //AliVParticle* tGenTz = (AliVParticle*) particleListGen.At(triggersMCb[ib]);
1367 AliVParticle* tGenTz = (AliVParticle*) particleListGen.At(ib);
6c1489e1 1368 if(!tGenTz) continue;
58003b11 1369 if(tGenTz->Pt()*0.95<15.0) continue;
1370 if(tGenTz->Pt()*0.95>100.0) continue;
6c1489e1 1371
1372 deltaPhib = RelativePhi(tGenTb->Phi(),tGenTz->Phi());
1373 deltaEtab = tGenTb->Eta()-tGenTz->Eta();
1374 deltaRb = sqrt(deltaPhib*deltaPhib + deltaEtab*deltaEtab);
1375
1376 if(deltaRb<0.4) bb++;
1377 }
1378 }
1379 fhNofMultipleTriggersConeGenB->Fill(bb);
1380 }
1381 }
1382
1383
1384
8aa19603 1385 //============== Estimate bg in generated events ==============
1386 Double_t rhoFromCellMedianGen=0.0, rhoConeGen=0.0;
1387
1388 //Estimate rho from cone
1389 EstimateBgCone(fListJetsGen, &particleListGen, rhoConeGen);
1390
1391 //Estimate rho from cell median minus jets
d1405a52 1392 EstimateBgRhoMedian(fListJetsBgGen, &particleListGen, rhoFromCellMedianGen,1);//mc data
8aa19603 1393
255cb71c 1394 //============ Generator trigger+jet ==================
8e103a56 1395 if(fHardest==0){ //single inclusive trigger
1396 if(ntriggersMC>0){ //there is at least one trigger
1397 Int_t rnd = fRandom->Integer(ntriggersMC); //0 to ntriggers-1
1398 indexTriggGen = triggersMC[rnd];
1e3ad34d 1399
1400 fhNofMultipleTriggersGen->Fill(ntriggersMC-1);
0fe881ea 1401
ac4f24f7 1402 Double_t deltaPhi, deltaEta, deltaR;
1403 Int_t k = 0;
1404
1405 //Correlation with single inclusive TRIGGER
1406 AliVParticle* tGenT1 = (AliVParticle*) particleListGen.At(indexTriggGen);
1407 if(tGenT1){
58003b11 1408 //for(Int_t ia=0; ia<ntriggersMC; ia++)
1409 for(Int_t ia=0; ia< particleListGen.GetEntries(); ia++){
1410 //if(indexTriggGen == triggersMC[ia]) continue;
1411 if(indexTriggGen == ia) continue;
ac4f24f7 1412
58003b11 1413 //AliVParticle* tGenT2 = (AliVParticle*) particleListGen.At(triggersMC[ia]);
1414 AliVParticle* tGenT2 = (AliVParticle*) particleListGen.At(ia);
ac4f24f7 1415 if(!tGenT2) continue;
58003b11 1416 if(tGenT2->Pt()<15.0) continue;
1417 if(tGenT2->Pt()>100.0) continue;
ac4f24f7 1418 deltaPhi = RelativePhi(tGenT1->Phi(),tGenT2->Phi());
1419 deltaEta = tGenT1->Eta()-tGenT2->Eta();
1420 deltaR = sqrt(deltaPhi*deltaPhi + deltaEta*deltaEta);
0fe881ea 1421
ac4f24f7 1422 if(deltaR<0.4) k++;
0fe881ea 1423 }
ac4f24f7 1424 }
1425 fhNofMultipleTriggersConeGen->Fill(k);
1426
1427 if(ntriggersMC>1){
0fe881ea 1428 //Correlation of each trigger with any other trigger
1e3ad34d 1429 for(Int_t ia=0; ia<ntriggersMC-1; ia++){
0fe881ea 1430 AliVParticle* tGenI = (AliVParticle*) particleListGen.At(triggersMC[ia]);
1e3ad34d 1431 if(!tGenI) continue;
1432 for(Int_t ib=ia+1; ib<ntriggersMC; ib++){
0fe881ea 1433 AliVParticle* tGenII = (AliVParticle*) particleListGen.At(triggersMC[ib]);
1e3ad34d 1434 if(!tGenII) continue;
1435
1436 deltaPhi = RelativePhi(tGenI->Phi(),tGenII->Phi());
1437 deltaEta = tGenI->Eta()-tGenII->Eta();
1438 deltaR = sqrt(deltaPhi*deltaPhi + deltaEta*deltaEta);
1439 fhDeltaRMultTriggersGen->Fill(deltaR);
1440 }
1441 }
1442 }
1443
80ac66f6 1444 }else{
8e103a56 1445 indexTriggGen = -1; //trigger not found
80ac66f6 1446 }
1447 }
8e103a56 1448
255cb71c 1449 //-----------
8e103a56 1450 Int_t ilowGen = (fHardest==0 || fHardest==1) ? indexTriggGen : 0;
1451 Int_t ihighGen = (fHardest==0 || fHardest==1) ? indexTriggGen+1 : particleListGen.GetEntries();
255cb71c 1452 Bool_t fillOnceGen = kTRUE;
1453 //-----------
1454
8e103a56 1455 for(Int_t igen= ilowGen; igen< ihighGen; igen++){ //loop over possible trigger
1456 indexTriggGen = igen; //trigger hadron
1457
1458 if(indexTriggGen == -1) continue;
1459 AliVParticle* triggerGen=(AliVParticle*) particleListGen.At(indexTriggGen);
1460 if(!triggerGen) continue;
1461
1462 if(fHardest >= 2){
6168afc3 1463 if(triggerGen->Pt() < 6.0) continue; //all hadrons pt>6
8e103a56 1464 }
1465 if(TMath::Abs((Float_t) triggerGen->Eta()) > fTriggerEtaCut) continue;
1466
1467 ptTriggGen = triggerGen->Pt(); //count triggers
1468 fh2NtriggersGen->Fill(centValue, ptTriggGen);
1469
80ac66f6 1470 //Count jets and trigger-jet pairs at MC generator level
1471 for(Int_t ij=0; ij<fListJetsGen->GetEntries(); ij++){
1472 AliAODJet* jet = (AliAODJet*)(fListJetsGen->At(ij));
1473 if(!jet) continue;
1474 Double_t etaJetGen = jet->Eta();
255cb71c 1475 Double_t areaJetGen = jet->EffectiveAreaCharged();
80ac66f6 1476
1477 if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){
80ac66f6 1478
49dc2855 1479 Double_t ptUeFromCellMedianGen = rhoFromCellMedianGen*areaJetGen;
1480 Double_t ptUeConeGen = rhoConeGen*areaJetGen;
1481
1482 //Rongrong's analysis
1483 Double_t dPhiGen = jet->Phi() - triggerGen->Phi();
1484 if(dPhiGen>2*TMath::Pi()) dPhiGen -= 2*TMath::Pi();
1485 if(dPhiGen<-2*TMath::Pi()) dPhiGen += 2*TMath::Pi();
1486 if(dPhiGen<-0.5*TMath::Pi()) dPhiGen += 2*TMath::Pi();
1487 if(dPhiGen>1.5*TMath::Pi()) dPhiGen -= 2*TMath::Pi();
1488
1489 //[A,pTjet,pTtrig,dphi]
1490 tmpArrayFour[0] = areaJetGen;
1491 tmpArrayFour[1] = jet->Pt();
1492 tmpArrayFour[2] = ptTriggGen;
1493 tmpArrayFour[3] = dPhiGen;
1494
1495 fHJetPhiCorrGen->Fill(tmpArrayFour); // Dphi distribution jet-triger
1496
1497 //[A,pTjet-pTUe,pTtrig,dphi]
1498 tmpArrayFour[1] = jet->Pt() - ptUeFromCellMedianGen;
1499 fHJetPhiCorrSubUeMedianGen->Fill(tmpArrayFour);
1500
1501 //[A,pTjet-pTUe,pTtrig,dphi]
1502 tmpArrayFour[1] = jet->Pt() - ptUeConeGen;
1503 fHJetPhiCorrSubUeConeGen->Fill(tmpArrayFour);
1504
1505 //Leticia's analysis
8e103a56 1506 Double_t dphi = RelativePhi(triggerGen->Phi(), jet->Phi());
1507 if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue;
80ac66f6 1508
8e103a56 1509 //Centrality, A, pT, pTtrigg
255cb71c 1510 tmpArrayFive[0] = centValue;
1511 tmpArrayFive[1] = areaJetGen;
1512 tmpArrayFive[2] = jet->Pt();
1513 tmpArrayFive[3] = ptTriggGen;
1514 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1515 fHJetSpecGen->Fill(tmpArrayFive);
1516
8aa19603 1517 //Centrality, A, pTjet-pTbgCellMedian, pTtrigg, dphi
255cb71c 1518 tmpArrayFive[0] = centValue;
1519 tmpArrayFive[1] = areaJetGen;
8aa19603 1520 tmpArrayFive[2] = jet->Pt() - ptUeFromCellMedianGen;
255cb71c 1521 tmpArrayFive[3] = ptTriggGen;
1522 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1523 fHJetSpecSubUeMedianGen->Fill(tmpArrayFive);
1524
8aa19603 1525 //Centrality, A, pTjet-pTbgCone, pTtrigg, dphi
255cb71c 1526 tmpArrayFive[0] = centValue;
1527 tmpArrayFive[1] = areaJetGen;
8aa19603 1528 tmpArrayFive[2] = jet->Pt() - ptUeConeGen;
255cb71c 1529 tmpArrayFive[3] = ptTriggGen;
1530 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
8aa19603 1531 fHJetSpecSubUeConeGen->Fill(tmpArrayFive);
255cb71c 1532
1533 //Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", kT median
1534 tmpArrayFive[0] = areaJetGen;
1535 tmpArrayFive[1] = jet->Pt();
8aa19603 1536 tmpArrayFive[2] = jet->Pt() - ptUeFromCellMedianGen;
1537 tmpArrayFive[3] = ptUeFromCellMedianGen;
1538 tmpArrayFive[4] = rhoFromCellMedianGen;
255cb71c 1539 fHJetUeMedianGen->Fill(tmpArrayFive);
1540
8aa19603 1541 //Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", perpendicular Cone
255cb71c 1542 tmpArrayFive[0] = areaJetGen;
1543 tmpArrayFive[1] = jet->Pt();
8aa19603 1544 tmpArrayFive[2] = jet->Pt() - ptUeConeGen;
1545 tmpArrayFive[3] = ptUeConeGen;
1546 tmpArrayFive[4] = rhoConeGen;
1547 fHJetUeConeGen->Fill(tmpArrayFive);
255cb71c 1548
1549 if(fillOnceGen){
8aa19603 1550 Double_t fillRhoGen[] = { rhoConeGen, rhoFromCellMedianGen};
1551 fHRhoUeMedianVsConeGen->Fill(fillRhoGen);
255cb71c 1552 fillOnceGen = kFALSE;
1553 }
8e103a56 1554 }//back to back jet-trigger pair
1555 }//jet loop
1556 }//trigger loop
1557
49dc2855 1558 if(fFillRespMx && !fIsKine){
8aa19603 1559 //================ RESPONSE MATRIX ===============
1560 //Count jets and trigger-jet pairs at MC generator level
1561 for(Int_t ij=0; ij<fListJetsGen->GetEntries(); ij++){
1562 AliAODJet* jet = (AliAODJet*)(fListJetsGen->At(ij));
1563 if(!jet) continue;
1564 Double_t etaJetGen = jet->Eta();
1565 Double_t ptJetGen = jet->Pt();
255cb71c 1566
8aa19603 1567 if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){
1568 fhJetPtGen->Fill(ptJetGen); // gen level pt spectum of jets response mx normalization
80ac66f6 1569
8aa19603 1570 Double_t areaJetGen = jet->EffectiveAreaCharged();
1571 Double_t ptUeFromCellMedianGen = rhoFromCellMedianGen*areaJetGen;
1572 Double_t ptUeConeGen = rhoConeGen*areaJetGen;
80ac66f6 1573
8aa19603 1574 fhJetPtSubUeMedianGen->Fill(ptJetGen - ptUeFromCellMedianGen);
1575 fhJetPtSubUeConeGen->Fill(ptJetGen - ptUeConeGen);
160dbdb8 1576 }
1577 }
8aa19603 1578 if(fListJets->GetEntries()>0 && fListJetsGen->GetEntries()>0){ //at least some reconstructed jets
1579
1580 Int_t ng = (Int_t) fListJetsGen->GetEntries();
1581 Int_t nr = (Int_t) fListJets->GetEntries();
160dbdb8 1582
8aa19603 1583 //Find closest MC generator - reconstructed jets
1584 if(faGenIndex.GetSize()<nr) faGenIndex.Set(nr); //idx of gen jet assoc to rec jet
1585 if(faRecIndex.GetSize()<ng) faRecIndex.Set(ng); //idx of rec jet assoc to gen jet
160dbdb8 1586
1587 if(fDebug){
8aa19603 1588 Printf("New Rec List %d gen index Array %d",nr,faGenIndex.GetSize());
1589 Printf("New Gen List %d rec index Array %d",ng,faRecIndex.GetSize());
160dbdb8 1590 }
1591 //matching of MC genrator level and reconstructed jets
8aa19603 1592 AliAnalysisHelperJetTasks::GetClosestJets(fListJetsGen,ng,fListJets,nr,faGenIndex,faRecIndex,fDebug);
1593
1594 // Fill response matrix
1595 for(Int_t ir = 0; ir < nr; ir++){
1596 AliAODJet *recJet = (AliAODJet*) fListJets->At(ir);
1597 Double_t etaJetRec = recJet->Eta();
1598 Double_t ptJetRec = recJet->Pt();
1599 Double_t areaJetRec = recJet->EffectiveAreaCharged();
160dbdb8 1600 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1601
8aa19603 1602 if((fJetEtaMin <= etaJetRec) && (etaJetRec <= fJetEtaMax)){
1603 Int_t ig = faGenIndex[ir]; //associated generator level jet
1604 if(ig >= 0 && ig < ng){
1605 if(fDebug > 10) Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1606 AliAODJet *genJet = (AliAODJet*) fListJetsGen->At(ig);
1607 Double_t ptJetGen = genJet->Pt();
1608 Double_t etaJetGen = genJet->Eta();
160dbdb8 1609
255cb71c 1610 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
8aa19603 1611 if((fJetEtaMin <= etaJetGen) && (etaJetGen <= fJetEtaMax)){
1612 fhJetPtGenVsJetPtRec->Fill(ptJetRec, ptJetGen);
1613
1614 Double_t areaJetGen = genJet->EffectiveAreaCharged();
1615 Double_t ptUeFromCellMedianGen = rhoFromCellMedianGen*areaJetGen;
1616 Double_t ptUeConeGen = rhoConeGen*areaJetGen;
1617 Double_t ptUeFromCellMedianRec = rhoFromCellMedian*areaJetRec;
1618 Double_t ptUeConeRec = rhoCone*areaJetRec;
1619 fhJetPtGenVsJetPtRecSubUeMedian->Fill(ptJetRec-ptUeFromCellMedianRec,
1620 ptJetGen-ptUeFromCellMedianGen);
1621 fhJetPtGenVsJetPtRecSubUeCone->Fill(ptJetRec-ptUeConeRec, ptJetGen-ptUeConeGen);
160dbdb8 1622 }
8aa19603 1623 }//ig>=0
160dbdb8 1624 }//rec jet in eta acceptance
1625 }//loop over reconstructed jets
1626 }// # of rec jets >0
8aa19603 1627
1628 //=========================== Full jet vs charged jet matrix ==========================
1629 if(fIsFullMC){
1630 //Count full jets and charged-jet pairs at MC generator level
1631 for(Int_t ij=0; ij<fListJetsGenFull->GetEntries(); ij++){
1632 AliAODJet* jetFull = (AliAODJet*)(fListJetsGenFull->At(ij));
1633 if(!jetFull) continue;
1634 Double_t etaJetFull = jetFull->Eta();
1635 Double_t ptJetFull = jetFull->Pt();
1636
1637 if((fJetEtaMin<=etaJetFull) && (etaJetFull<=fJetEtaMax)){
1638 fhJetPtGenFull->Fill(ptJetFull); // generator level pt spectum of full jets
1639 }
1640 }
1641 if(fListJetsGen->GetEntries()>0 && fListJetsGenFull->GetEntries()>0){ //at least some reconstructed jets
1642 Int_t nful = (Int_t) fListJetsGenFull->GetEntries();
1643 Int_t nchr = (Int_t) fListJetsGen->GetEntries();
1644
1645 //Find closest MC generator full - charged jet
1646 if(faGenIndex.GetSize()<nchr) faGenIndex.Set(nchr); //idx of gen FULL jet assoc to gen CHARGED jet
1647 if(faRecIndex.GetSize()<nful) faRecIndex.Set(nful); //idx of gen CHARGED jet assoc to gen FULL jet
1648
1649 if(fDebug){
1650 Printf("New Charg List %d Full index Array %d",nchr,faGenIndex.GetSize());
1651 Printf("New Full List %d Charg index Array %d",nful,faRecIndex.GetSize());
1652 }
1653 //matching of MC genrator level and reconstructed jets
1654 AliAnalysisHelperJetTasks::GetClosestJets(fListJetsGenFull,nful,fListJetsGen,nchr,faGenIndex,faRecIndex,fDebug);
1655
1656 // Fill response matrix
1657 for(Int_t ichr = 0; ichr < nchr; ichr++){ //charged jet loop
1658 AliAODJet *chJet = (AliAODJet*) fListJetsGen->At(ichr);
1659 Double_t etaJetCh = chJet->Eta();
1660 Double_t ptJetCh = chJet->Pt();
1661 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1662
1663 if((fJetEtaMin <= etaJetCh) && (etaJetCh <= fJetEtaMax)){
1664 Int_t iful = faGenIndex[ichr]; //associated generator level jet
1665 if(iful >= 0 && iful < nful){
1666 if(fDebug > 10) Printf("%s:%d iful = %d ichr = %d",(char*)__FILE__,__LINE__,iful,ichr);
1667 AliAODJet *genJetFull = (AliAODJet*) fListJetsGenFull->At(iful);
1668 Double_t ptJetFull = genJetFull->Pt();
1669 Double_t etaJetFull = genJetFull->Eta();
1670
1671 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1672 if((fJetEtaMin <= etaJetFull) && (etaJetFull <= fJetEtaMax)){
1673 fhJetPtGenChargVsJetPtGenFull->Fill(ptJetFull,ptJetCh);
1674 }
1675 }//iful>=0
1676 }//rec jet in eta acceptance
1677 }//loop over reconstructed jets
1678 }// # of rec jets >0
1679 }//pointer MC generator jets
1680 } //fill resp mx only for bin
1681 }//analyze generator level MC
1682
1683
49dc2855 1684 if(fIsKine){ //skip reconstructed data analysis in case of kine
1685 PostData(1, fOutputList);
1686 return;
1687 }
8e103a56 1688 //============= RECONSTRUCTED INCLUSIVE JETS ===============
ad869500 1689
ad869500 1690 Double_t etaJet = 0.0;
1691 Double_t pTJet = 0.0;
1692 Double_t areaJet = 0.0;
1693 Double_t phiJet = 0.0;
9f10e189 1694 //Int_t indexLeadingJet = -1;
1695 //Double_t pTLeadingJet = -10.0;
1696 //Double_t areaLeadingJet = -10.0;
8e103a56 1697
ad869500 1698 for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
1699 AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
1700 if(!jet) continue;
1701 etaJet = jet->Eta();
1702 phiJet = jet->Phi();
1703 pTJet = jet->Pt();
1704 if(pTJet==0) continue;
1705
1706 if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
255cb71c 1707 /*areaJet = jet->EffectiveAreaCharged();*/
7fc3d134 1708
ad869500 1709 //Jet Diagnostics---------------------------------
8e103a56 1710 fhJetPhi->Fill(pTJet, RelativePhi(phiJet,0.0)); //phi -pi,pi
1711 fhJetEta->Fill(pTJet, etaJet);
7fc3d134 1712 //search for leading jet
255cb71c 1713 /*if(pTJet > pTLeadingJet){
7fc3d134 1714 indexLeadingJet = ij;
1715 pTLeadingJet = pTJet;
7fc3d134 1716 areaLeadingJet = areaJet;
1717 }
1718
1719 // raw spectra of INCLUSIVE jets
80ac66f6 1720 //Centrality, pTjet, A
7fc3d134 1721 Double_t fillraw[] = { centValue,
1722 pTJet,
7fc3d134 1723 areaJet
1724 };
255cb71c 1725 fHJetPtRaw->Fill(fillraw);*/
8e103a56 1726 }
255cb71c 1727 /*
8e103a56 1728 if(indexLeadingJet > -1){
1729 // raw spectra of LEADING jets
1730 //Centrality, pTjet, A
1731 Double_t fillleading[] = { centValue,
1732 pTLeadingJet,
1733 areaLeadingJet
1734 };
1735 fHLeadingJetPtRaw->Fill(fillleading);
1736 }
255cb71c 1737 */
8e103a56 1738
1739 // =============== RECONSTRUCTED TRIGGER-JET PAIRS ================
8aa19603 1740 if(fIsChargedMC && fFillRespMx){
1741 FillEffHistos(&particleList, &particleListGen); //Fill efficiency histos
1742 }
255cb71c 1743 Bool_t filledOnce = kTRUE; //fill rho histogram only once per event
8e103a56 1744 //set ranges of the trigger loop
1745 Int_t ilow = (fHardest==0 || fHardest==1) ? indexTrigg : 0;
1746 Int_t ihigh = (fHardest==0 || fHardest==1) ? indexTrigg+1 : particleList.GetEntries();
7fc3d134 1747
8e103a56 1748 for(Int_t itrk= ilow; itrk< ihigh; itrk++){ //loop over possible trigger
1749 indexTrigg = itrk; //trigger hadron with pT >10 GeV
7fc3d134 1750
8e103a56 1751 if(indexTrigg < 0) continue;
1752
1753 AliVParticle *triggerHadron = (AliVParticle*) particleList.At(indexTrigg);
1754 if(!triggerHadron){
1755 PostData(1, fOutputList);
1756 return;
1757 }
1758 if(fHardest >= 2){
7d9fd386 1759 if(triggerHadron->Pt() < 6.0) continue; //all hadrons pt>6
8e103a56 1760 }
1761 if(TMath::Abs((Float_t) triggerHadron->Eta())> fTriggerEtaCut) continue;
7fc3d134 1762
8e103a56 1763 //Fill trigger histograms
1764 fh2Ntriggers->Fill(centValue,triggerHadron->Pt()); //trigger pT
1765 fhTriggerPhi->Fill(triggerHadron->Pt(),RelativePhi(triggerHadron->Phi(),0.0)); //phi -pi,pi
1766 fhTriggerEta->Fill(triggerHadron->Pt(),triggerHadron->Eta());
7fc3d134 1767
1768
8e103a56 1769 //---------- make trigger-jet pairs ---------
255cb71c 1770 //Int_t injet4 = 0;
1771 //Int_t injet = 0;
255cb71c 1772
8e103a56 1773 for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
1774 AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
1775 if(!jet) continue;
1776 etaJet = jet->Eta();
1777 phiJet = jet->Phi();
1778 pTJet = jet->Pt();
1779 if(pTJet==0) continue;
1780
1781 if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
1782 areaJet = jet->EffectiveAreaCharged();
49dc2855 1783
1784 //subtract bg using estinates base on median of kT jets
1785 Double_t ptUeFromCellMedian = rhoFromCellMedian*areaJet;
1786 Double_t ptUeCone = rhoCone*areaJet;
1787
255cb71c 1788 //if(areaJet >= 0.07) injet++;
1789 //if(areaJet >= 0.4) injet4++;
8e103a56 1790
1791 Double_t dphi = RelativePhi(triggerHadron->Phi(), phiJet);
1792 fhDphiTriggerJet->Fill(dphi); //Input
1793
1794 //Dphi versus jet pT
1795 //Centrality, Dphi=phiTrig-phiJet, pTjet, pTtrigg
255cb71c 1796 /*Double_t filldp[] = { centValue,
8e103a56 1797 dphi,
1798 pTJet,
1799 triggerHadron->Pt()
1800 };
1801 fHDphiVsJetPtAll->Fill(filldp);
255cb71c 1802 */
49dc2855 1803 //Rongrong's analysis
1804
1805 Double_t dPhi = phiJet - triggerHadron->Phi();
1806 if(dPhi>2*TMath::Pi()) dPhi -= 2*TMath::Pi();
1807 if(dPhi<-2*TMath::Pi()) dPhi += 2*TMath::Pi();
1808 if(dPhi<-0.5*TMath::Pi()) dPhi += 2*TMath::Pi();
1809 if(dPhi>1.5*TMath::Pi()) dPhi -= 2*TMath::Pi();
1810
1811 //[A,pTjet,pTtrig,dphi]
1812 tmpArrayFour[0] = areaJet;
1813 tmpArrayFour[1] = pTJet;
1814 tmpArrayFour[2] = triggerHadron->Pt();
1815 tmpArrayFour[3] = dPhi;
1816
1817 fHJetPhiCorr->Fill(tmpArrayFour); // Dphi distribution jet-triger
1818
1819 //[A,pTjet-pTUe,pTtrig,dphi]
1820 tmpArrayFour[1] = pTJet - ptUeFromCellMedian;
1821 fHJetPhiCorrSubUeMedian->Fill(tmpArrayFour);
1822
1823 //[A,pTjet-pTUe,pTtrig,dphi]
1824 tmpArrayFour[1] = pTJet - ptUeCone;
1825 fHJetPhiCorrSubUeCone->Fill(tmpArrayFour);
1826
1827 //Leticia's analysis
1828
8e103a56 1829 // Select back to back trigger - jet pairs
1830 if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue;
1831 fhDphiTriggerJetAccept->Fill(dphi); //Accepted
255cb71c 1832
1833 //NO bg subtraction
1834 //Centrality, A, pTjet, pTtrigg, dphi
1835 tmpArrayFive[0] = centValue;
1836 tmpArrayFive[1] = areaJet;
1837 tmpArrayFive[2] = pTJet;
1838 tmpArrayFive[3] = triggerHadron->Pt();
1839 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1840 fHJetSpec->Fill(tmpArrayFive);
1841
49dc2855 1842 //Centrality, A, pTjet-pTbgCellMedian, pTtrigg, dphi
255cb71c 1843 tmpArrayFive[0] = centValue;
1844 tmpArrayFive[1] = areaJet;
8aa19603 1845 tmpArrayFive[2] = pTJet-ptUeFromCellMedian;
255cb71c 1846 tmpArrayFive[3] = triggerHadron->Pt();
1847 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1848 fHJetSpecSubUeMedian->Fill(tmpArrayFive);
1849
8aa19603 1850 //Centrality, A, pTjet-pTbgCone, pTtrigg, dphi
255cb71c 1851 tmpArrayFive[0] = centValue;
1852 tmpArrayFive[1] = areaJet;
8aa19603 1853 tmpArrayFive[2] = pTJet - ptUeCone;
255cb71c 1854 tmpArrayFive[3] = triggerHadron->Pt();
1855 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
8aa19603 1856 fHJetSpecSubUeCone->Fill(tmpArrayFive);
255cb71c 1857
1858 //Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", kT median
1859 tmpArrayFive[0] = areaJet;
1860 tmpArrayFive[1] = pTJet;
8aa19603 1861 tmpArrayFive[2] = pTJet - ptUeFromCellMedian;
1862 tmpArrayFive[3] = ptUeFromCellMedian;
1863 tmpArrayFive[4] = rhoFromCellMedian;
255cb71c 1864 fHJetUeMedian->Fill(tmpArrayFive);
8e103a56 1865
8aa19603 1866 //Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", Cone median
255cb71c 1867 tmpArrayFive[0] = areaJet;
1868 tmpArrayFive[1] = pTJet;
8aa19603 1869 tmpArrayFive[2] = pTJet - ptUeCone;
1870 tmpArrayFive[3] = ptUeCone;
1871 tmpArrayFive[4] = rhoCone;
1872 fHJetUeCone->Fill(tmpArrayFive);
255cb71c 1873
1874 if(filledOnce){ //fill for each event only once
8aa19603 1875 Double_t fillRho[] = { rhoCone,rhoFromCellMedian};
1876 fHRhoUeMedianVsCone->Fill(fillRho);
255cb71c 1877 filledOnce = kFALSE;
1878 }
8e103a56 1879 }//jet loop
1880
1881 //Fill Jet Density In the Event A>0.07
255cb71c 1882 /*if(injet>0){
8e103a56 1883 Double_t filldens[]={ (Double_t) particleList.GetEntries(),
ad869500 1884 injet/fkAcceptance,
1885 triggerHadron->Pt()
1886 };
8e103a56 1887 fHJetDensity->Fill(filldens);
1888 }
ad869500 1889
8e103a56 1890 //Fill Jet Density In the Event A>0.4
1891 if(injet4>0){
1892 Double_t filldens4[]={ (Double_t) particleList.GetEntries(),
1893 injet4/fkAcceptance,
1894 triggerHadron->Pt()
1895 };
1896 fHJetDensityA4->Fill(filldens4);
255cb71c 1897 }*/
ad869500 1898 }
1899
8e103a56 1900
ad869500 1901 PostData(1, fOutputList);
1902}
1903
1904//----------------------------------------------------------------------------
1905void AliAnalysisTaskJetCorePP::Terminate(const Option_t *)
1906{
1907 // Draw result to the screen
1908 // Called once at the end of the query
1909
1910 if(fDebug) printf("AliAnalysisTaskJetCorePP DONE\n");
1911 if(!GetOutputData(1)) return;
1912}
1913
1914
1915//----------------------------------------------------------------------------
1916Int_t AliAnalysisTaskJetCorePP::GetListOfTracks(TList *list){
1917 //Fill the list of accepted tracks (passed track cut)
1918 //return consecutive index of the hardest ch hadron in the list
1919 Int_t iCount = 0;
1920 AliAODEvent *aodevt = NULL;
1921
1922 if(!fESD) aodevt = fAODIn;
1923 else aodevt = fAODOut;
1924
1925 if(!aodevt) return -1;
1926
1927 Int_t index = -1; //index of the highest particle in the list
1928 Double_t ptmax = -10;
8e103a56 1929 Int_t triggers[200];
1930 Int_t ntriggers = 0; //index in triggers array
1e3ad34d 1931
ad869500 1932
80ac66f6 1933 for(Int_t it = 0; it < aodevt->GetNumberOfTracks(); it++){
0fe881ea 1934 //AliAODTrack *tr = aodevt->GetTrack(it);
f15c1f69 1935 AliAODTrack *tr = dynamic_cast<AliAODTrack*>(aodevt->GetTrack(it));
1936 if(!tr) AliFatal("Not a standard AOD");
0fe881ea 1937
d1405a52
ML		 1938 if((fFilterMask > 0) && !(tr->TestFilterBit(fFilterMask))) continue;
1939 //if((fFilterMask > 0) && !(tr->IsHybridGlobalConstrainedGlobal())) continue;
ad869500 1940 if(TMath::Abs((Float_t) tr->Eta()) > fTrackEtaCut) continue;
1941 if(tr->Pt() < fTrackLowPtCut) continue;
1942 list->Add(tr);
8e103a56 1943
1944 //Search trigger:
1945 if(fHardest>0){ //leading particle
1946 if(tr->Pt()>ptmax){
1947 ptmax = tr->Pt();
1948 index = iCount;
1949 }
1950 }
1951
1952 if(fHardest==0 && ntriggers<200){ //single inclusive
1953 if(fTriggerPtRangeLow <= tr->Pt() &&
1954 tr->Pt() < fTriggerPtRangeHigh){
1955 triggers[ntriggers] = iCount;
1956 ntriggers++;
1957 }
ad869500 1958 }
8e103a56 1959
ad869500 1960 iCount++;
1961 }
1962
8e103a56 1963 if(fHardest==0 && ntriggers>0){ //select random inclusive trigger
1964 Int_t rnd = fRandom->Integer(ntriggers); //0 to ntriggers-1
1965 index = triggers[rnd];
1e3ad34d 1966
1967 fhNofMultipleTriggers->Fill(ntriggers-1);
ac4f24f7 1968
1969 Double_t deltaPhi, deltaEta, deltaR;
1970 Int_t k=0;
1971 //Correlation with single inclusive trigger
1972 AliVParticle* tGent1 = (AliVParticle*) list->At(index);
1973 if(tGent1){
58003b11 1974 //for(Int_t ia=0; ia<ntriggers; ia++)
1975 for(Int_t ia=0; ia<list->GetEntries(); ia++){
1976 //if(triggers[ia]==index) continue;
1977 if(ia==index) continue;
1978 //AliVParticle* tGent2 = (AliVParticle*) list->At(triggers[ia]);
1979 AliVParticle* tGent2 = (AliVParticle*) list->At(ia);
ac4f24f7 1980 if(!tGent2) continue;
58003b11 1981 if(tGent2->Pt()<15.0) continue;
1982 if(tGent2->Pt()>100.0) continue;
ac4f24f7 1983 deltaPhi = RelativePhi(tGent1->Phi(),tGent2->Phi());
1984 deltaEta = tGent1->Eta()-tGent2->Eta();
1985 deltaR = sqrt(deltaPhi*deltaPhi + deltaEta*deltaEta);
1986 if(deltaR<0.4) k++;
1987 }
1988 }
1989 fhNofMultipleTriggersCone->Fill(k);
0fe881ea 1990
ac4f24f7 1991 if(ntriggers>1){
0fe881ea 1992 //Correlation with any other trigger
1e3ad34d 1993 for(Int_t ia=0; ia<ntriggers-1; ia++){
0fe881ea 1994 AliVParticle* tGeni = (AliVParticle*) list->At(triggers[ia]);
1e3ad34d 1995 if(!tGeni) continue;
1996 for(Int_t ib=ia+1; ib<ntriggers; ib++){
0fe881ea 1997 AliVParticle* tGenii = (AliVParticle*) list->At(triggers[ib]);
1e3ad34d 1998 if(!tGenii) continue;
1999 deltaPhi = RelativePhi(tGeni->Phi(),tGenii->Phi());
2000 deltaEta = tGeni->Eta()-tGenii->Eta();
2001 deltaR = sqrt(deltaPhi*deltaPhi + deltaEta*deltaEta);
2002 fhDeltaRMultTriggers->Fill(deltaR);
2003 }
2004 }
2005 }
ad869500 2006 }
8e103a56 2007
2008 return index;
ad869500 2009}
2010
2011//----------------------------------------------------------------------------
80ac66f6 2012/*
ad869500 2013Double_t AliAnalysisTaskJetCorePP::GetBackgroundInPerpCone(Float_t jetR, Double_t jetPhi, Double_t jetEta, TList* trkList){
2014 //calculate sum of track pT in the cone perpendicular in phi to the jet
2015 //jetR = cone radius
2016 // jetPhi, jetEta = direction of the jet
2017 Int_t numberOfTrks = trkList->GetEntries();
2018 Double_t pTsum = 0.0;
2019 Double_t perpConePhi = jetPhi + TMath::Pi()/2;//perp cone w.r.t. jet in phi
2020 for(Int_t it=0; it<numberOfTrks; it++){
2021 AliVParticle *trk = (AliVParticle*) trkList->At(it);
2022 Double_t dphi = RelativePhi(perpConePhi,trk->Phi());
2023 Double_t deta = trk->Eta()-jetEta;
2024
2025 if( (dphi*dphi + deta*deta)< (jetR*jetR)){
2026 pTsum += trk->Pt();
2027 }
2028 }
2029
2030 return pTsum;
2031}
80ac66f6 2032*/
ad869500 2033//----------------------------------------------------------------------------
2034
2035Double_t AliAnalysisTaskJetCorePP::RelativePhi(Double_t mphi,Double_t vphi){
2036 //Get relative azimuthal angle of two particles -pi to pi
2037 if (vphi < -TMath::Pi()) vphi += TMath::TwoPi();
2038 else if (vphi > TMath::Pi()) vphi -= TMath::TwoPi();
2039
2040 if (mphi < -TMath::Pi()) mphi += TMath::TwoPi();
2041 else if (mphi > TMath::Pi()) mphi -= TMath::TwoPi();
2042
2043 Double_t dphi = mphi - vphi;
2044 if (dphi < -TMath::Pi()) dphi += TMath::TwoPi();
2045 else if (dphi > TMath::Pi()) dphi -= TMath::TwoPi();
2046
2047 return dphi;//dphi in [-Pi, Pi]
2048}
2049
2050
80ac66f6 2051//----------------------------------------------------------------------------
2052Bool_t AliAnalysisTaskJetCorePP::SelectMCGenTracks(AliVParticle *trk, TList *trkList, Double_t &ptLeading, Int_t &index, Int_t counter){
2053 //fill track efficiency denominator
2054 if(trk->Pt() < fTrackLowPtCut) return kFALSE;
2055 if(trk->Charge()==0) return kFALSE;
2056 if(TMath::Abs((Float_t) trk->Eta()) > fTrackEtaCut) return kFALSE;
2057 trkList->Add(trk);
8aa19603 2058 if(fFillRespMx) fhPtTrkTruePrimGen->Fill(trk->Pt(),trk->Eta()); //Efficiency denominator
8e103a56 2059
2060 //Search trigger:
2061 if(fHardest>0){ //leading particle
2062 if(ptLeading < trk->Pt()){
2063 index = counter;
2064 ptLeading = trk->Pt();
2065 }
80ac66f6 2066 }
8e103a56 2067
2068 if(fHardest==0){ //single inclusive
2069 index = -1;
2070 if(fTriggerPtRangeLow <= trk->Pt() &&
2071 trk->Pt() < fTriggerPtRangeHigh){
2072 index = counter;
2073 }
2074 }
2075
80ac66f6 2076 return kTRUE;
2077}
2078
2079//----------------------------------------------------------------------------
2080void AliAnalysisTaskJetCorePP::FillEffHistos(TList *recList, TList *genList){
2081 //fill track efficiency numerator
2082 if(!recList) return;
2083 if(!genList) return;
2084 Int_t nRec = recList->GetEntries();
2085 Int_t nGen = genList->GetEntries();
2086 for(Int_t ir=0; ir<nRec; ir++){
2087 AliAODTrack *trkRec = (AliAODTrack*) recList->At(ir);
2088 if(!trkRec) continue;
2089 Int_t recLabel = TMath::Abs(trkRec->GetLabel());
2090 Bool_t matched = kFALSE;
2091
2092 for(Int_t ig=0; ig<nGen; ig++){
2093 AliVParticle *trkGen = (AliVParticle*) genList->At(ig);
2094 if(!trkGen) continue;
2095 Int_t genLabel = TMath::Abs(trkGen->GetLabel());
2096 if(recLabel==genLabel){
2097 fhPtTrkTruePrimRec->Fill(trkGen->Pt(),trkGen->Eta());
2098 matched = kTRUE;
2099 break;
2100 }
2101 }
2102 if(!matched) fhPtTrkSecOrFakeRec->Fill(trkRec->Pt(),trkRec->Eta());
2103 }
2104
2105 return;
2106}
255cb71c 2107//________________________________________________________________
d1405a52 2108void AliAnalysisTaskJetCorePP::EstimateBgRhoMedian(TList *listJet, TList* listPart, Double_t &rhoMedian, Int_t mode){
8aa19603 2109 //Estimate background rho by means of integrating track pT outside identified jet cones
2110
2111 rhoMedian = 0;
2112
2113 //identify leading jet in the full track acceptance
2114 Int_t idxLeadingJet = -1;
2115 Double_t pTleading = 0.0;
2116 AliAODJet* jet = NULL;
2117
2118 for(Int_t ij = 0; ij < listJet->GetEntries(); ij++){ //loop over bg kT jets
2119 jet = (AliAODJet*)(listJet->At(ij));
2120 if(!jet) continue;
2121 if(TMath::Abs(jet->Eta()) > fTrackEtaCut) continue; //track acceptance cut
2122 if(jet->Pt() > pTleading){
2123 idxLeadingJet = ij;
2124 pTleading = jet->Pt();
2125 }
2126 }
2127
2128 Int_t njets = 0;
2129 static Double_t jphi[100];
2130 static Double_t jeta[100];
2131 static Double_t jRsquared[100];
2132
2133 for(Int_t ij = 0; ij< listJet->GetEntries(); ij++){ //loop over bg kT jets
2134
2135 jet = (AliAODJet*)(listJet->At(ij));
2136 if(!jet) continue;
2137 if(TMath::Abs(jet->Eta()) > fTrackEtaCut) continue;
2138
2139 if((jet->Pt() > fBgMaxJetPt || ij== idxLeadingJet) && njets<100){
2140 jphi[njets] = RelativePhi(jet->Phi(),0.0); //-pi,pi
2141 jeta[njets] = jet->Eta();
2142 jRsquared[njets] = fSafetyMargin*jet->EffectiveAreaCharged()/TMath::Pi(); //scale jet radius
2143 njets++;
255cb71c 2144 }
2145 }
8aa19603 2146
255cb71c 2147
8aa19603 2148 static Double_t nOutCone[10][4];
2149 static Double_t sumPtOutOfCone[10][4];
2150
2151 for(Int_t ie=0; ie<fnEta; ie++){
2152 for(Int_t ip=0; ip<fnPhi; ip++){
2153 nOutCone[ip][ie] = 0.0; //initialize counter
2154 sumPtOutOfCone[ip][ie] = 0.0;
255cb71c 2155 }
8aa19603 2156 }
255cb71c 2157
8aa19603 2158 Double_t rndphi, rndeta;
2159 Double_t rndphishift, rndetashift;
2160 Double_t dphi, deta;
2161 Bool_t bIsInCone;
2162
2163 for(Int_t it=0; it<fnTrials; it++){
2164
2165 rndphi = fRandom->Uniform(0, fPhiSize);
2166 rndeta = fRandom->Uniform(0, fEtaSize);
2167
2168 for(Int_t ip=0; ip<fnPhi; ip++){ //move radom position to each cell
2169 rndphishift = rndphi + ip*fPhiSize - TMath::Pi();
2170 for(Int_t ie=0; ie<fnEta; ie++){
2171 rndetashift = rndeta + ie*fEtaSize - fEtaSize;
2172
2173 bIsInCone = 0; //tag if trial is in the jet cone
2174 for(Int_t ij=0; ij<njets; ij++){
2175 deta = jeta[ij] - rndetashift;
2176 dphi = RelativePhi(rndphishift,jphi[ij]);
2177 if((dphi*dphi + deta*deta) < jRsquared[ij]){
2178 bIsInCone = 1;
2179 break;
2180 }
2181 }
2182 if(!bIsInCone) nOutCone[ip][ie]++;
2183 }
2184 }
2185 }
2186
255cb71c 2187
8aa19603 2188 //Sum particle pT outside jets in cells
2189 Int_t npart = listPart->GetEntries();
2190 Int_t phicell,etacell;
2191 AliVParticle *part = NULL;
2192 for(Int_t ip=0; ip < npart; ip++){
2193
2194 part = (AliVParticle*) listPart->At(ip);
2195 if(!part) continue;
2196 if( TMath::Abs(part->Eta()) > fEtaSize) continue; //skip part outside +-0.7 in eta
2197
2198 bIsInCone = 0; //init
2199 for(Int_t ij=0; ij<njets; ij++){
2200 dphi = RelativePhi(jphi[ij], part->Phi());
2201 deta = jeta[ij] - part->Eta();
2202 if((dphi*dphi + deta*deta) < jRsquared[ij]){
2203 bIsInCone = 1;
2204 break;
2205 }
255cb71c 2206 }
8aa19603 2207 if(!bIsInCone){
2208 phicell = TMath::Nint(TMath::Floor((RelativePhi(part->Phi(),0.0) + TMath::Pi())/fPhiSize));
2209 etacell = TMath::Nint(TMath::Floor((part->Eta()+fEtaSize)/fEtaSize));
2210 sumPtOutOfCone[phicell][etacell]+= part->Pt();
2211 }
2212 }
2213
2214 static Double_t rhoInCells[20];
2215 Double_t relativeArea;
2216 Int_t nCells=0;
d1405a52 2217 Double_t bufferArea=0.0, bufferPt=0.0; //sum cells where A< fJetFreeAreaFrac
8aa19603 2218 for(Int_t ip=0; ip<fnPhi; ip++){
2219 for(Int_t ie=0; ie<fnEta; ie++){
2220 relativeArea = nOutCone[ip][ie]/fnTrials;
d1405a52
ML		 2221 //cout<<"RA= "<< relativeArea<<" BA="<<bufferArea<<" BPT="<< bufferPt <<endl;
2222
2223 bufferArea += relativeArea;
2224 bufferPt += sumPtOutOfCone[ip][ie];
2225 if(bufferArea > fJetFreeAreaFrac){
2226 rhoInCells[nCells] = bufferPt/(bufferArea*fCellArea);
2227
2228 if(mode==1) fhCellAreaToMedianGen->Fill(bufferArea*fCellArea);
2229 else fhCellAreaToMedian->Fill(bufferArea*fCellArea);
2230
2231 bufferArea = 0.0;
2232 bufferPt = 0.0;
2233 nCells++;
2234 }
8aa19603 2235 }
255cb71c 2236 }
8aa19603 2237
d1405a52
ML		 2238
2239 if(nCells>0){
2240 rhoMedian = TMath::Median(nCells, rhoInCells);
2241 if(mode==1){ //mc data
2242 fhEntriesToMedianGen->Fill(nCells);
2243 }else{ //real data
2244 fhEntriesToMedian->Fill(nCells);
2245 }
2246 }
8aa19603 2247
255cb71c 2248}
80ac66f6 2249
255cb71c 2250//__________________________________________________________________
8aa19603 2251void AliAnalysisTaskJetCorePP::EstimateBgCone(TList *listJet, TList* listPart, Double_t &rhoPerpCone){
2252
2253 rhoPerpCone = 0.0; //init
2254
2255 if(!listJet) return;
2256 Int_t njet = listJet->GetEntries();
2257 Int_t npart = listPart->GetEntries();
2258 Double_t pTleading = 0.0;
2259 Double_t phiLeading = 1000.;
2260 Double_t etaLeading = 1000.;
2261
2262 for(Int_t jsig=0; jsig < njet; jsig++){
2263 AliAODJet* signaljet = (AliAODJet*)(listJet->At(jsig));
2264 if(!signaljet) continue;
2265 if((signaljet->Eta()<fJetEtaMin) && (fJetEtaMax<signaljet->Eta())) continue; //acceptance cut
2266 if(signaljet->Pt() >= pTleading){ //replace leading and subleading jet
2267 pTleading = signaljet->Pt();
2268 phiLeading = signaljet->Phi();
2269 etaLeading = signaljet->Eta();
2270 }
2271 }
2272
2273 Double_t phileftcone = phiLeading + TMath::Pi()/2;
2274 Double_t phirightcone = phiLeading - TMath::Pi()/2;
2275 Double_t dp, de;
2276
2277 for(Int_t ip=0; ip < npart; ip++){
2278
2279 AliVParticle *part = (AliVParticle*) listPart->At(ip);
2280 if(!part){
2281 continue;
2282 }
2283
2284 dp = RelativePhi(phileftcone, part->Phi());
2285 de = etaLeading - part->Eta();
2286 if( sqrt(dp*dp + de*de)< fBgConeR) rhoPerpCone += part->Pt();
2287
2288 dp = RelativePhi(phirightcone, part->Phi());
2289 if( sqrt(dp*dp + de*de)< fBgConeR) rhoPerpCone += part->Pt();
2290
2291 }
2292
2293 //normalize total pT by two times cone are
2294 rhoPerpCone = rhoPerpCone/(2*TMath::Pi()*fBgConeR*fBgConeR);
2295
2296}
2297//__________________________________________________________________
80ac66f6 2298
255cb71c 2299void AliAnalysisTaskJetCorePP::ReadTClonesArray(TString bname, TList *list){
2300 //Convert TClones array of jets to TList
2301
2302 if(!strlen(bname.Data())){
2303 AliError(Form("Jet branch %s not set.", bname.Data()));
2304 return;
2305 }
2306
2307 TClonesArray *array=0x0;
49dc2855 2308 if(fUseExchContainer){ //take input from exchange containers
2309 if(fJetBranchNameKine.Length()>0 && bname.CompareTo(fJetBranchNameKine.Data())==0){
2310 array = dynamic_cast<TClonesArray*>(GetInputData(1)); //connect exchange container slot 1
2311 }
3d1fac8a 2312 if(fJetBranchNameBgKine.Length()>0 && bname.CompareTo(fJetBranchNameBgKine.Data())==0){
2313 array = dynamic_cast<TClonesArray*>(GetInputData(2)); //connect exchange container slot 2
2314 }
49dc2855 2315 }else{ //take input from AOD
2316 if(fAODOut&&!array){
2317 array = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(bname.Data()));
2318 }
2319 if(fAODExtension&&!array){
2320 array = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(bname.Data()));
2321 }
2322 if(fAODIn&&!array){
2323 array = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(bname.Data()));
2324 }
255cb71c 2325 }
2326
2327 list->Clear(); //clear list beforehand
2328
2329 if(array){
2330 if(fDebug){
2331 Printf("########## %s: %d jets",bname.Data(), array->GetEntriesFast());
2332 }
49dc2855 2333
255cb71c 2334 for(Int_t iJet = 0; iJet < array->GetEntriesFast(); iJet++) {
2335 AliAODJet *jet = dynamic_cast<AliAODJet*>((*array)[iJet]);
2336 if (jet) list->Add(jet);
2337 }
2338 }
2339
2340 return;
2341}