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