Trigger Associated with tracks 15<pT<100
[u/mrichter/AliRoot.git] / PWGJE / AliAnalysisTaskJetCorePP.cxx
... / ...
CommitLineData
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"
28#include "TClonesArray.h"
29#include "TTree.h"
30#include "TList.h"
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"
38#include "TArrayI.h"
39#include "TProfile.h"
40#include "TFile.h"
41#include "TKey.h"
42#include "TRandom3.h"
43
44#include "AliLog.h"
45
46#include "AliAnalysisTask.h"
47#include "AliAnalysisManager.h"
48
49#include "AliVEvent.h"
50#include "AliESDEvent.h"
51#include "AliMCEvent.h"
52#include "AliESDInputHandler.h"
53#include "AliCentrality.h"
54#include "AliAnalysisHelperJetTasks.h"
55#include "AliInputEventHandler.h"
56#include "AliAODJetEventBackground.h"
57#include "AliGenPythiaEventHeader.h"
58#include "AliAODMCParticle.h"
59#include "AliMCParticle.h"
60#include "AliAODEvent.h"
61#include "AliAODHandler.h"
62#include "AliAODJet.h"
63#include "AliVVertex.h"
64#include "AliVTrack.h"
65#include "AliAnalysisTaskJetCorePP.h"
66#include "AliHeader.h" //KF//
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),
82fMcEvent(0x0),
83fMcHandler(0x0),
84fJetBranchName(""),
85fJetBranchNameChargMC(""),
86fJetBranchNameKine(""),
87fJetBranchNameFullMC(""),
88fJetBranchNameBg(""),
89fJetBranchNameBgChargMC(""),
90fJetBranchNameBgKine(""),
91fListJets(0x0),
92fListJetsGen(0x0),
93fListJetsGenFull(0x0),
94fListJetsBg(0x0),
95fListJetsBgGen(0x0),
96fNonStdFile(""),
97fSystem(0), //pp=0 pPb=1
98fJetParamR(0.4),
99fBgJetParamR(0.3),
100fBgMaxJetPt(8.0),
101fBgConeR(0.4),
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),
114fUseExchContainer(0),
115fOutputList(0x0),
116fHistEvtSelection(0x0),
117fh2Ntriggers(0x0),
118fHJetSpec(0x0),
119fHJetSpecSubUeMedian(0x0),
120fHJetSpecSubUeCone(0x0),
121fHJetPhiCorr(0x0),
122fHJetPhiCorrSubUeMedian(0x0),
123fHJetPhiCorrSubUeCone(0x0),
124fHJetUeMedian(0x0),
125fHJetUeCone(0x0),
126fHRhoUeMedianVsCone(0x0),
127//fHJetDensity(0x0),
128//fHJetDensityA4(0x0),
129fhJetPhi(0x0),
130fhTriggerPhi(0x0),
131fhJetEta(0x0),
132fhTriggerEta(0x0),
133fhVertexZ(0x0),
134fhVertexZAccept(0x0),
135fhContribVtx(0x0),
136fhContribVtxAccept(0x0),
137fhDphiTriggerJet(0x0),
138fhDphiTriggerJetAccept(0x0),
139fhCentrality(0x0),
140fhCentralityAccept(0x0),
141fhNofMultipleTriggers(0x0),
142fhNofMultipleTriggersCone(0x0),
143fhDeltaRMultTriggers(0x0),
144//fHJetPtRaw(0x0),
145//fHLeadingJetPtRaw(0x0),
146//fHDphiVsJetPtAll(0x0),
147fhJetPtGenVsJetPtRec(0x0),
148fhJetPtGenVsJetPtRecSubUeMedian(0x0),
149fhJetPtGenVsJetPtRecSubUeCone(0x0),
150fhJetPtGen(0x0),
151fhJetPtSubUeMedianGen(0x0),
152fhJetPtSubUeConeGen(0x0),
153fhJetPtGenChargVsJetPtGenFull(0x0),
154fhJetPtGenFull(0x0),
155fh2NtriggersGen(0x0),
156fHJetSpecGen(0x0),
157fHJetSpecSubUeMedianGen(0x0),
158fHJetSpecSubUeConeGen(0x0),
159fHJetPhiCorrGen(0x0),
160fHJetPhiCorrSubUeMedianGen(0x0),
161fHJetPhiCorrSubUeConeGen(0x0),
162fHJetUeMedianGen(0x0),
163fHJetUeConeGen(0x0),
164fhPtTrkTruePrimRec(0x0),
165fhPtTrkTruePrimGen(0x0),
166fhPtTrkSecOrFakeRec(0x0),
167fHRhoUeMedianVsConeGen(0x0),
168fhEntriesToMedian(0x0),
169fhEntriesToMedianGen(0x0),
170fhCellAreaToMedian(0x0),
171fhCellAreaToMedianGen(0x0),
172fhNofMultipleTriggersGen(0x0),
173fhNofMultipleTriggersConeGen(0x0),
174fhDeltaRMultTriggersGen(0x0),
175fhNofMultipleTriggersConeGenA(0x0),
176fhNofMultipleTriggersConeGenB(0x0),
177fIsChargedMC(0),
178fIsKine(0),
179fIsFullMC(0),
180faGenIndex(0),
181faRecIndex(0),
182fkAcceptance(2.0*TMath::Pi()*1.8),
183fkDeltaPhiCut(TMath::Pi()-0.8),
184fh1Xsec(0x0),
185fh1Trials(0x0),
186fh1AvgTrials(0x0),
187fh1PtHard(0x0),
188fh1PtHardNoW(0x0),
189fh1PtHardTrials(0x0),
190fAvgTrials(1),
191fHardest(0),
192fEventNumberRangeLow(0),
193fEventNumberRangeHigh(99),
194fTriggerPtRangeLow(0.0),
195fTriggerPtRangeHigh(10000.0),
196fFillRespMx(0),
197fRandom(0x0),
198fnTrials(1000),
199fJetFreeAreaFrac(0.5),
200fnPhi(9),
201fnEta(2),
202fEtaSize(0.7),
203fPhiSize(2*TMath::Pi()/fnPhi),
204fCellArea(fPhiSize*fEtaSize),
205fSafetyMargin(1.1),
206fDoubleBinning(kFALSE)
207{
208 // default Constructor
209}
210
211//---------------------------------------------------------------------
212
213AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP(const char *name) :
214AliAnalysisTaskSE(name),
215fESD(0x0),
216fAODIn(0x0),
217fAODOut(0x0),
218fAODExtension(0x0),
219fMcEvent(0x0),
220fMcHandler(0x0),
221fJetBranchName(""),
222fJetBranchNameChargMC(""),
223fJetBranchNameKine(""),
224fJetBranchNameFullMC(""),
225fJetBranchNameBg(""),
226fJetBranchNameBgChargMC(""),
227fJetBranchNameBgKine(""),
228fListJets(0x0),
229fListJetsGen(0x0),
230fListJetsGenFull(0x0),
231fListJetsBg(0x0),
232fListJetsBgGen(0x0),
233fNonStdFile(""),
234fSystem(0), //pp=0 pPb=1
235fJetParamR(0.4),
236fBgJetParamR(0.3),
237fBgMaxJetPt(8.0),
238fBgConeR(0.4),
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),
251fUseExchContainer(0),
252fOutputList(0x0),
253fHistEvtSelection(0x0),
254fh2Ntriggers(0x0),
255fHJetSpec(0x0),
256fHJetSpecSubUeMedian(0x0),
257fHJetSpecSubUeCone(0x0),
258fHJetPhiCorr(0x0),
259fHJetPhiCorrSubUeMedian(0x0),
260fHJetPhiCorrSubUeCone(0x0),
261fHJetUeMedian(0x0),
262fHJetUeCone(0x0),
263fHRhoUeMedianVsCone(0x0),
264//fHJetDensity(0x0),
265//fHJetDensityA4(0x0),
266fhJetPhi(0x0),
267fhTriggerPhi(0x0),
268fhJetEta(0x0),
269fhTriggerEta(0x0),
270fhVertexZ(0x0),
271fhVertexZAccept(0x0),
272fhContribVtx(0x0),
273fhContribVtxAccept(0x0),
274fhDphiTriggerJet(0x0),
275fhDphiTriggerJetAccept(0x0),
276fhCentrality(0x0),
277fhCentralityAccept(0x0),
278fhNofMultipleTriggers(0x0),
279fhNofMultipleTriggersCone(0x0),
280fhDeltaRMultTriggers(0x0),
281//fHJetPtRaw(0x0),
282//fHLeadingJetPtRaw(0x0),
283//fHDphiVsJetPtAll(0x0),
284fhJetPtGenVsJetPtRec(0x0),
285fhJetPtGenVsJetPtRecSubUeMedian(0x0),
286fhJetPtGenVsJetPtRecSubUeCone(0x0),
287fhJetPtGen(0x0),
288fhJetPtSubUeMedianGen(0x0),
289fhJetPtSubUeConeGen(0x0),
290fhJetPtGenChargVsJetPtGenFull(0x0),
291fhJetPtGenFull(0x0),
292fh2NtriggersGen(0x0),
293fHJetSpecGen(0x0),
294fHJetSpecSubUeMedianGen(0x0),
295fHJetSpecSubUeConeGen(0x0),
296fHJetPhiCorrGen(0x0),
297fHJetPhiCorrSubUeMedianGen(0x0),
298fHJetPhiCorrSubUeConeGen(0x0),
299fHJetUeMedianGen(0x0),
300fHJetUeConeGen(0x0),
301fhPtTrkTruePrimRec(0x0),
302fhPtTrkTruePrimGen(0x0),
303fhPtTrkSecOrFakeRec(0x0),
304fHRhoUeMedianVsConeGen(0x0),
305fhEntriesToMedian(0x0),
306fhEntriesToMedianGen(0x0),
307fhCellAreaToMedian(0x0),
308fhCellAreaToMedianGen(0x0),
309fhNofMultipleTriggersGen(0x0),
310fhNofMultipleTriggersConeGen(0x0),
311fhDeltaRMultTriggersGen(0x0),
312fhNofMultipleTriggersConeGenA(0x0),
313fhNofMultipleTriggersConeGenB(0x0),
314fIsChargedMC(0),
315fIsKine(0),
316fIsFullMC(0),
317faGenIndex(0),
318faRecIndex(0),
319fkAcceptance(2.0*TMath::Pi()*1.8),
320fkDeltaPhiCut(TMath::Pi()-0.8),
321fh1Xsec(0x0),
322fh1Trials(0x0),
323fh1AvgTrials(0x0),
324fh1PtHard(0x0),
325fh1PtHardNoW(0x0),
326fh1PtHardTrials(0x0),
327fAvgTrials(1),
328fHardest(0),
329fEventNumberRangeLow(0),
330fEventNumberRangeHigh(99),
331fTriggerPtRangeLow(0.0),
332fTriggerPtRangeHigh(10000.0),
333fFillRespMx(0),
334fRandom(0x0),
335fnTrials(1000),
336fJetFreeAreaFrac(0.5),
337fnPhi(9),
338fnEta(2),
339fEtaSize(0.7),
340fPhiSize(2*TMath::Pi()/fnPhi),
341fCellArea(fPhiSize*fEtaSize),
342fSafetyMargin(1.1),
343fDoubleBinning(kFALSE)
344{
345// Constructor
346 DefineOutput(1, TList::Class());
347
348 TString dummy(name);
349 if(dummy.Contains("KINE")){
350 DefineInput(1, TClonesArray::Class());
351 DefineInput(2, TClonesArray::Class());
352 }
353}
354
355//--------------------------------------------------------------
356AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP(const AliAnalysisTaskJetCorePP& a):
357AliAnalysisTaskSE(a.GetName()),
358fESD(a.fESD),
359fAODIn(a.fAODIn),
360fAODOut(a.fAODOut),
361fAODExtension(a.fAODExtension),
362fMcEvent(a.fMcEvent),
363fMcHandler(a.fMcHandler),
364fJetBranchName(a.fJetBranchName),
365fJetBranchNameChargMC(a.fJetBranchNameChargMC),
366fJetBranchNameKine(a.fJetBranchNameKine),
367fJetBranchNameFullMC(a.fJetBranchNameFullMC),
368fJetBranchNameBg(a.fJetBranchNameBg),
369fJetBranchNameBgChargMC(a.fJetBranchNameBgChargMC),
370fJetBranchNameBgKine(a.fJetBranchNameBgKine),
371fListJets(a.fListJets),
372fListJetsGen(a.fListJetsGen),
373fListJetsGenFull(a.fListJetsGenFull),
374fListJetsBg(a.fListJetsBg),
375fListJetsBgGen(a.fListJetsBgGen),
376fNonStdFile(a.fNonStdFile),
377fSystem(a.fSystem),
378fJetParamR(a.fJetParamR),
379fBgJetParamR(a.fBgJetParamR),
380fBgMaxJetPt(a.fBgMaxJetPt),
381fBgConeR(a.fBgConeR),
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),
394fUseExchContainer(a.fUseExchContainer),
395fOutputList(a.fOutputList),
396fHistEvtSelection(a.fHistEvtSelection),
397fh2Ntriggers(a.fh2Ntriggers),
398fHJetSpec(a.fHJetSpec),
399fHJetSpecSubUeMedian(a.fHJetSpecSubUeMedian),
400fHJetSpecSubUeCone(a.fHJetSpecSubUeCone),
401fHJetPhiCorr(a.fHJetPhiCorr),
402fHJetPhiCorrSubUeMedian(a.fHJetPhiCorrSubUeMedian),
403fHJetPhiCorrSubUeCone(a.fHJetPhiCorrSubUeCone),
404fHJetUeMedian(a.fHJetUeMedian),
405fHJetUeCone(a.fHJetUeCone),
406fHRhoUeMedianVsCone(a.fHRhoUeMedianVsCone),
407//fHJetDensity(a.fHJetDensity),
408//fHJetDensityA4(a.fHJetDensityA4),
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),
421fhNofMultipleTriggers(a.fhNofMultipleTriggers),
422fhNofMultipleTriggersCone(a.fhNofMultipleTriggersCone),
423fhDeltaRMultTriggers(a.fhDeltaRMultTriggers),
424//fHJetPtRaw(a.fHJetPtRaw),
425//fHLeadingJetPtRaw(a.fHLeadingJetPtRaw),
426//fHDphiVsJetPtAll(a.fHDphiVsJetPtAll),
427fhJetPtGenVsJetPtRec(a.fhJetPtGenVsJetPtRec),
428fhJetPtGenVsJetPtRecSubUeMedian(a.fhJetPtGenVsJetPtRecSubUeMedian),
429fhJetPtGenVsJetPtRecSubUeCone(a.fhJetPtGenVsJetPtRecSubUeCone),
430fhJetPtGen(a.fhJetPtGen),
431fhJetPtSubUeMedianGen(a.fhJetPtSubUeMedianGen),
432fhJetPtSubUeConeGen(a.fhJetPtSubUeConeGen),
433fhJetPtGenChargVsJetPtGenFull(a.fhJetPtGenChargVsJetPtGenFull),
434fhJetPtGenFull(a.fhJetPtGenFull),
435fh2NtriggersGen(a.fh2NtriggersGen),
436fHJetSpecGen(a.fHJetSpecGen),
437fHJetSpecSubUeMedianGen(a.fHJetSpecSubUeMedianGen),
438fHJetSpecSubUeConeGen(a.fHJetSpecSubUeConeGen),
439fHJetPhiCorrGen(a.fHJetPhiCorrGen),
440fHJetPhiCorrSubUeMedianGen(a.fHJetPhiCorrSubUeMedianGen),
441fHJetPhiCorrSubUeConeGen(a.fHJetPhiCorrSubUeConeGen),
442fHJetUeMedianGen(a.fHJetUeMedianGen),
443fHJetUeConeGen(a.fHJetUeConeGen),
444fhPtTrkTruePrimRec(a.fhPtTrkTruePrimRec),
445fhPtTrkTruePrimGen(a.fhPtTrkTruePrimGen),
446fhPtTrkSecOrFakeRec(a.fhPtTrkSecOrFakeRec),
447fHRhoUeMedianVsConeGen(a.fHRhoUeMedianVsConeGen),
448fhEntriesToMedian(a.fhEntriesToMedian),
449fhEntriesToMedianGen(a.fhEntriesToMedianGen),
450fhCellAreaToMedian(a.fhCellAreaToMedian),
451fhCellAreaToMedianGen(a.fhCellAreaToMedianGen),
452fhNofMultipleTriggersGen(a.fhNofMultipleTriggersGen),
453fhNofMultipleTriggersConeGen(a.fhNofMultipleTriggersConeGen),
454fhDeltaRMultTriggersGen(a.fhDeltaRMultTriggersGen),
455fhNofMultipleTriggersConeGenA(a.fhNofMultipleTriggersConeGenA),
456fhNofMultipleTriggersConeGenB(a.fhNofMultipleTriggersConeGenB),
457fIsChargedMC(a.fIsChargedMC),
458fIsKine(a.fIsKine),
459fIsFullMC(a.fIsFullMC),
460faGenIndex(a.faGenIndex),
461faRecIndex(a.faRecIndex),
462fkAcceptance(a.fkAcceptance),
463fkDeltaPhiCut(a.fkDeltaPhiCut),
464fh1Xsec(a.fh1Xsec),
465fh1Trials(a.fh1Trials),
466fh1AvgTrials(a.fh1AvgTrials),
467fh1PtHard(a.fh1PtHard),
468fh1PtHardNoW(a.fh1PtHardNoW),
469fh1PtHardTrials(a.fh1PtHardTrials),
470fAvgTrials(a.fAvgTrials),
471fHardest(a.fHardest),
472fEventNumberRangeLow(a.fEventNumberRangeLow),
473fEventNumberRangeHigh(a.fEventNumberRangeHigh),
474fTriggerPtRangeLow(a.fTriggerPtRangeLow),
475fTriggerPtRangeHigh(a.fTriggerPtRangeHigh),
476fFillRespMx(a.fFillRespMx),
477fRandom(a.fRandom),
478fnTrials(a.fnTrials),
479fJetFreeAreaFrac(a.fJetFreeAreaFrac),
480fnPhi(a.fnPhi),
481fnEta(a.fnEta),
482fEtaSize(a.fEtaSize),
483fPhiSize(a.fPhiSize),
484fCellArea(a.fCellArea),
485fSafetyMargin(a.fSafetyMargin),
486fDoubleBinning(a.fDoubleBinning)
487{
488 //Copy Constructor
489 fRandom->SetSeed(0);
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;
505 delete fListJetsGen;
506 delete fListJetsGenFull;
507 delete fListJetsBg;
508 delete fListJetsBgGen;
509 delete fOutputList; // ????
510 delete fRandom;
511}
512
513//--------------------------------------------------------------
514
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()
521 if(!(fIsChargedMC || fIsKine)) return kFALSE;
522 Float_t xsection = 0;
523 Float_t trials = 1;
524 fAvgTrials = 1;
525
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 }
532
533 fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
534
535
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 }
544
545 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
546
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 }
567
568
569 return kTRUE;
570}
571//--------------------------------------------------------------
572
573void AliAnalysisTaskJetCorePP::Init()
574{
575 // check for jet branches
576 if(fJetBranchNameKine.Length()==0){
577 if(!strlen(fJetBranchName.Data())){
578 AliError("Jet branch name not set.");
579 }
580 }
581 fMcHandler = dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
582
583
584}
585
586//--------------------------------------------------------------
587
588void AliAnalysisTaskJetCorePP::UserCreateOutputObjects()
589{
590 // Create histograms and initilize variables
591 fSafetyMargin = fBgConeR*fBgConeR /(fBgJetParamR*fBgJetParamR);
592
593 // Called once
594 fListJets = new TList(); //reconstructed level antikt jets
595 fListJetsBg = new TList(); //reconstructed jets to be removed from bg
596
597 fIsChargedMC = (fJetBranchNameChargMC.Length()>0) ? kTRUE : kFALSE;
598 fIsKine = (fJetBranchNameKine.Length()>0) ? kTRUE : kFALSE;
599 fIsFullMC = (fJetBranchNameFullMC.Length()>0) ? kTRUE : kFALSE;
600
601 fRandom = new TRandom3(0);
602
603 if(fIsChargedMC || fIsKine){ //full MC or pure kine
604 fListJetsGen = new TList(); //generator level charged antikt jets
605 fListJetsBgGen = new TList(); //generator level jets to be removed from bg
606
607 if(fIsFullMC)
608 fListJetsGenFull = new TList(); //generator level full jets
609 }
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)");
623 fHistEvtSelection->GetXaxis()->SetBinLabel(6,"event number (rejected)");
624
625 fOutputList->Add(fHistEvtSelection);
626
627 Int_t nBinsCentrality = (fSystem==0) ? 1 : 10; // pp=1 else 10
628 //trigger pt spectrum (reconstructed)
629 fh2Ntriggers = new TH2F("fh2Ntriggers","# of triggers",
630 nBinsCentrality,0.0,100.0,50,0.0,50.0);
631 if(!fIsKine) fOutputList->Add(fh2Ntriggers);
632
633 Int_t bw = (fDoubleBinning==0) ? 1 : 2; //make larger bin width
634
635 //Centrality, A, pTjet, pTtrigg, dphi
636 const Int_t dimSpec = 5;
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()};
640 fHJetSpec = new THnSparseF("fHJetSpec",
641 "Recoil jet spectrum [cent,A,pTjet,pTtrig,dphi]",
642 dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
643 if(!fIsKine) fOutputList->Add(fHJetSpec);
644
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]");
648 if(!fIsKine) fOutputList->Add(fHJetSpecSubUeMedian);
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]");
652 if(!fIsKine) fOutputList->Add(fHJetSpecSubUeCone);
653
654
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);
674
675
676 //------------------- HISTOS FOR DIAGNOSTIC ----------------------
677 //A, pTjet, pTjet-pTUe, pTUe, rhoUe bg estimate from kT median
678 const Int_t dimSpecMed = 5;
679 const Int_t nBinsSpecMed[dimSpecMed] = {25, 50, 60, 20, 20};
680 const Double_t lowBinSpecMed[dimSpecMed] = {0.0, 0.0, -20.0, 0.0, 0.0};
681 const Double_t hiBinSpecMed[dimSpecMed] = {1.0, 100.0, 100.0, 10.0, 20.0};
682 fHJetUeMedian = new THnSparseF("fHJetUeMedian",
683 "Recoil jet spectrum [A,pTjet,pTjet-pTUe, pTUe, rhoUe]",
684 dimSpecMed, nBinsSpecMed, lowBinSpecMed, hiBinSpecMed);
685 if(!fIsKine) fOutputList->Add(fHJetUeMedian);
686
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]");
690 if(!fIsKine) fOutputList->Add(fHJetUeCone);
691
692 //rho bacground reconstructed data
693 const Int_t dimRho = 2;
694 const Int_t nBinsRho[dimRho] = {50 , 50};
695 const Double_t lowBinRho[dimRho] = {0.0 , 0.0};
696 const Double_t hiBinRho[dimRho] = {20.0 , 20.0};
697
698 fHRhoUeMedianVsCone = new THnSparseF("hRhoUeMedianVsCone","[Rho Cone, Rho Median]",
699 dimRho, nBinsRho, lowBinRho, hiBinRho);
700 if(!fIsKine) fOutputList->Add(fHRhoUeMedianVsCone);
701
702 //Jet number density histos [Trk Mult, jet density, pT trigger]
703 /*const Int_t dimJetDens = 3;
704 const Int_t nBinsJetDens[dimJetDens] = {100, 100, 10};
705 const Double_t lowBinJetDens[dimJetDens] = {0.0, 0.0, 0.0};
706 const Double_t hiBinJetDens[dimJetDens] = {500.0, 5.0, 50.0};
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
714 fOutputList->Add(fHJetDensity);
715 fOutputList->Add(fHJetDensityA4);
716 */
717
718 //inclusive azimuthal and pseudorapidity histograms
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",
722 50, 0, 50, 50,-TMath::Pi(),TMath::Pi());
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",
726 50, 0, 50, 40,-0.9,0.9);
727
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);
736 fhEntriesToMedian = new TH1D("fhEntriesToMedian","fhEntriesToMedian",30,0,30);
737 fhCellAreaToMedian = new TH1D("fhCellAreaToMedian", "fhCellAreaToMedian", 75,0,1.5);
738
739 fhNofMultipleTriggers = new TH1D("fhNofMultipleTriggers","fhNofMultipleTriggers",100,0,100);
740 fhNofMultipleTriggersCone = new TH1D("fhNofMultipleTriggersCone","fhNofMultipleTriggersCone R<0.4",100,0,100);
741 fhDeltaRMultTriggers = new TH1D("fhDeltaRMultTriggers","fhDeltaRMultTriggers", 100,0,4);
742
743 if(!fIsKine){
744 fOutputList->Add(fhJetPhi);
745 fOutputList->Add(fhTriggerPhi);
746 fOutputList->Add(fhJetEta);
747 fOutputList->Add(fhTriggerEta);
748 }
749 fOutputList->Add(fhVertexZ);
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);
760 fOutputList->Add(fhNofMultipleTriggers);
761 fOutputList->Add(fhNofMultipleTriggersCone);
762 fOutputList->Add(fhDeltaRMultTriggers);
763 }
764 // raw spectra of INCLUSIVE jets
765 //Centrality, pTjet, A
766 /*const Int_t dimRaw = 3;
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};
770 fHJetPtRaw = new THnSparseF("fHJetPtRaw",
771 "Incl. jet spectrum [cent,pTjet,A]",
772 dimRaw,nBinsRaw,lowBinRaw,hiBinRaw);
773 fOutputList->Add(fHJetPtRaw);
774
775 // raw spectra of LEADING jets
776 //Centrality, pTjet, A
777 fHLeadingJetPtRaw = new THnSparseF("fHLeadingJetPtRaw",
778 "Leading jet spectrum [cent,pTjet,A]",
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);
792 */
793
794 //analyze MC generator level
795 if(fIsChargedMC || fIsKine){
796 if(fFillRespMx){
797 //Fill response matrix only once
798 fhJetPtGenVsJetPtRec = new TH2D("fhJetPtGenVsJetPtRec","JetPtGenVsJetPtRec", bw*100,0,200, bw*100,0,200);
799 fOutputList->Add(fhJetPtGenVsJetPtRec); //gen MC charg jet pt spectrum versus rec charged jet pt spectrum
800 //....
801 fhJetPtGenVsJetPtRecSubUeMedian = new TH2D("fhJetPtGenVsJetPtRecSubUeMedian","fhJetPtGenVsJetPtRecSubUeMedian", bw*110,-20,200, bw*110,-20,200);
802 fOutputList->Add(fhJetPtGenVsJetPtRecSubUeMedian); // with kT median bg subtr
803 //....
804 fhJetPtGenVsJetPtRecSubUeCone=(TH2D*)fhJetPtGenVsJetPtRecSubUeMedian->Clone("fhJetPtGenVsJetPtRecSubUeCone");
805 fhJetPtGenVsJetPtRecSubUeCone->SetTitle("fhJetPtGenVsJetPtRecSubUeCone");
806 fOutputList->Add(fhJetPtGenVsJetPtRecSubUeCone); // with weighted kT median bg subtr
807 //....
808 fhJetPtGen = new TH1D("fhJetPtGen","Jet Pt (MC Gen)",bw*100,0,200); //MC generator charged jet pt spectrum
809 fOutputList->Add(fhJetPtGen);
810 //....
811 fhJetPtSubUeMedianGen = new TH1D("fhJetPtSubUeMedianGen","Jet Pt - UE pT (MC Gen)",bw*110,-20,200);
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){
819 fhJetPtGenChargVsJetPtGenFull = new TH2D("fhJetPtGenChargVsJetPtGenFull","fhJetPtGenChargVsJetPtGenFull", 100,0,200, 100,0,200);
820 fOutputList->Add(fhJetPtGenChargVsJetPtGenFull); //gen full MC jet pt versus gen charged jet MC pt
821 //....
822 fhJetPtGenFull = new TH1D("fhJetPtGenFull","Jet Pt (MC Full jets Gen)",100,0,200); //MC generator full jet pt spectrum
823 fOutputList->Add(fhJetPtGenFull);
824 }
825 }
826 //....
827 fh2NtriggersGen = (TH2F*) fh2Ntriggers->Clone("fh2NtriggersGen");
828 fh2NtriggersGen->SetTitle(Form("%s Gen MC",fh2Ntriggers->GetTitle()));
829 fOutputList->Add(fh2NtriggersGen);
830
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
836 fHJetSpecSubUeMedianGen = (THnSparseF*) fHJetSpecSubUeMedian->Clone("fHJetSpecSubUeMedianGen");
837 fHJetSpecSubUeMedianGen->SetTitle(Form("%s Gen MC",fHJetSpecSubUeMedian->GetTitle()));
838 fOutputList->Add(fHJetSpecSubUeMedianGen);
839
840 fHJetSpecSubUeConeGen = (THnSparseF*) fHJetSpecSubUeCone->Clone("fHJetSpecSubUeConeGen");
841 fHJetSpecSubUeConeGen->SetTitle(Form("%s Gen MC",fHJetSpecSubUeCone->GetTitle()));
842 fOutputList->Add(fHJetSpecSubUeConeGen);
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
856 //---
857 fHJetUeMedianGen = (THnSparseF*) fHJetUeMedian->Clone("fHJetUeMedianGen");
858 fHJetUeMedianGen->SetTitle(Form("%s Gen MC", fHJetUeMedian->GetTitle()));
859 fOutputList->Add(fHJetUeMedianGen);
860
861 fHJetUeConeGen = (THnSparseF*) fHJetUeCone->Clone("fHJetUeConeGen");
862 fHJetUeConeGen->SetTitle(Form("%s Gen MC", fHJetUeCone->GetTitle()));
863 fOutputList->Add(fHJetUeConeGen);
864
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);
869
870 fhPtTrkTruePrimGen = (TH2D*) fhPtTrkTruePrimRec->Clone("fhPtTrkTruePrimGen");
871 fhPtTrkTruePrimGen->SetTitle("PtTrkTruePrimGen");
872 fOutputList->Add(fhPtTrkTruePrimGen);
873
874 fhPtTrkSecOrFakeRec = (TH2D*) fhPtTrkTruePrimRec->Clone("fhPtTrkSecOrFakeRec");
875 fhPtTrkSecOrFakeRec->SetTitle("PtTrkSecOrFakeRec");
876 fOutputList->Add(fhPtTrkSecOrFakeRec);
877 }
878
879 fHRhoUeMedianVsConeGen = (THnSparseF*) fHRhoUeMedianVsCone->Clone("hRhoUeMedianVsConeGen");
880 fHRhoUeMedianVsConeGen->SetTitle(Form("%s Gen MC", fHRhoUeMedianVsCone->GetTitle()));
881 fOutputList->Add(fHRhoUeMedianVsConeGen);
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);
890
891 fhNofMultipleTriggersGen = (TH1D*) fhNofMultipleTriggers->Clone("fhNofMultipleTriggersGen");
892 fOutputList->Add(fhNofMultipleTriggersGen);
893
894 fhNofMultipleTriggersConeGen = (TH1D*) fhNofMultipleTriggersCone->Clone("fhNofMultipleTriggersConeGen");
895 fOutputList->Add(fhNofMultipleTriggersConeGen);
896
897 fhDeltaRMultTriggersGen = (TH1D*) fhDeltaRMultTriggers->Clone("fhDeltaRMultTriggersGen");
898 fOutputList->Add(fhDeltaRMultTriggersGen);
899
900 fhNofMultipleTriggersConeGenA = (TH1D*) fhNofMultipleTriggersConeGen->Clone("fhNofMultipleTriggersConeGen10");
901 fOutputList->Add(fhNofMultipleTriggersConeGenA);
902
903 fhNofMultipleTriggersConeGenB = (TH1D*) fhNofMultipleTriggersConeGen->Clone("fhNofMultipleTriggersConeGen5");
904 fOutputList->Add(fhNofMultipleTriggersConeGenB);
905
906 }
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
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{
957 //User Exec
958
959
960 //Event loop
961 Double_t eventW = 1.0;
962 Double_t ptHard = 0.0;
963 Double_t nTrials = 1.0; // Trials for MC trigger
964 if(fIsChargedMC || fIsKine) fh1AvgTrials->Fill("#sum{avg ntrials}",fAvgTrials);
965
966 if(TMath::Abs((Float_t) fJetParamR)<0.00001){
967 AliError("ANTIKT Cone radius is set to zero.");
968 return;
969 }
970
971 if(TMath::Abs((Float_t) fBgJetParamR)<0.00001){
972 AliError("ANTIKT Cone radius is set to zero.");
973 return;
974 }
975
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
986 if(!strlen(fJetBranchNameBgKine.Data())){
987 AliError("Name of jet bg branch for kine not set.");
988 return;
989 }
990
991 Init();
992 if(fMcHandler){
993 fMcEvent = fMcHandler->MCEvent();
994 }else{
995 if(fDebug > 1) printf("AliAnalysisTaskJetCorePP::Exec() fMcHandler=NULL\n");
996 PostData(1, fOutputList);
997 return;
998 }
999 if(!fMcEvent){
1000 if(fDebug > 1) printf("AliAnalysisTaskJetCorePP::Exec() fMcEvent=NULL \n");
1001 PostData(1, fOutputList);
1002 return;
1003 }
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);
1020 }
1021
1022
1023 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
1024 if(!fESD){
1025 if(fDebug>1) AliError("ESD not available");
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
1032 if(!aod && !fIsKine){
1033 if(!fESD) aod = fAODIn; //not ESD and not kine => input AOD
1034 else aod = fAODOut;// ESD or kine
1035 }
1036
1037
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
1047 // ----------------- EVENT SELECTION --------------------------
1048 fHistEvtSelection->Fill(1); // number of events before event selection
1049
1050 if(!fIsKine){
1051 // physics selection
1052 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
1053 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
1054
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
1062 //check AOD pointer
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 }
1069
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 }
1079
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);
1085
1086 if((nTracksPrim < fMinContribVtx) ||
1087 (vtxz < fVtxZMin) ||
1088 (vtxz > fVtxZMax)){
1089 if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",
1090 (char*)__FILE__,__LINE__,vtxz);
1091 fHistEvtSelection->Fill(3);
1092 PostData(1, fOutputList);
1093 return;
1094 }else{
1095 //Accepted events
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 }
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
1121 //------------------ CENTRALITY SELECTION ---------------
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{
1129 //centValue = aod->GetHeader()->GetCentrality();
1130 centValue = ((AliVAODHeader*)aod->GetHeader())->GetCentrality();
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
1146 //-----------------select disjunct event subsamples ----------------
1147 if(!fIsKine){ //reconstructed data
1148 //Int_t eventnum = aod->GetHeader()->GetEventNumberESDFile();
1149 AliAODHeader * header = dynamic_cast<AliAODHeader*>(aod->GetHeader());
1150 if(!header) AliFatal("Not a standard AOD");
1151
1152 Int_t eventnum = header->GetEventNumberESDFile();
1153 Int_t lastdigit = eventnum % 10;
1154 if(!(fEventNumberRangeLow<=lastdigit && lastdigit<=fEventNumberRangeHigh)){
1155 fHistEvtSelection->Fill(5);
1156 PostData(1, fOutputList);
1157 return;
1158 }
1159 }
1160
1161 if(fDebug) std::cout<<" ACCEPTED EVENT "<<endl;
1162 fHistEvtSelection->Fill(0); // accepted events
1163 // ==================== end event selection ============================
1164
1165 Double_t tmpArrayFive[5];
1166 Double_t tmpArrayFour[4];
1167
1168
1169 TList particleList; //list of tracks
1170 Int_t indexTrigg = -1;
1171 Double_t rhoFromCellMedian=0.0, rhoCone=0.0;
1172
1173 if(!fIsKine){
1174 //=============== Reconstructed antikt jets ===============
1175 ReadTClonesArray(fJetBranchName.Data() , fListJets);
1176 ReadTClonesArray(fJetBranchNameBg.Data(), fListJetsBg);
1177
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 }
1187 //============== analyze generator level MC ================
1188 TList particleListGen; //list of tracks in MC
1189
1190 if(fIsChargedMC || fIsKine){
1191
1192 if(fIsKine){ //pure kine
1193
1194 //================= generated charged antikt jets ================
1195 ReadTClonesArray(fJetBranchNameKine.Data(), fListJetsGen);
1196 ReadTClonesArray(fJetBranchNameBgKine.Data(), fListJetsBgGen);
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 }
1206 //========================================================
1207 //serarch for charged trigger at the MC generator level
1208
1209 Int_t indexTriggGen = -1;
1210 Double_t ptTriggGen = -1;
1211 Int_t iCounterGen = 0; //number of entries in particleListGen array
1212 Int_t triggersMC[200];//list of trigger candidates
1213 Int_t ntriggersMC = 0; //index in triggers array
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
1219
1220 if(!fIsKine){
1221 if(fESD){//ESD input
1222
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 }
1237
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 }
1280 }
1281 }else{ //analyze kine
1282
1283 for(Int_t it = 0; it < fMcEvent->GetNumberOfTracks(); it++){
1284 if(!fMcEvent->IsPhysicalPrimary(it)) continue;
1285 AliMCParticle* part = (AliMCParticle*) fMcEvent->GetTrack(it);
1286 if(!part) continue;
1287 if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){
1288
1289 if(fHardest==0 && ntriggersMC<200){//single inclusive trigger
1290 if(indexTriggGen > -1){//trigger candidate was found
1291 triggersMC[ntriggersMC] = indexTriggGen;
1292 ntriggersMC++;
1293 }
1294 }
1295
1296 iCounterGen++;//index in particleListGen array
1297 }
1298 }
1299 }
1300
1301 if(fHardest==0){
1302 Int_t npar = particleListGen.GetEntries();
1303 for(Int_t ip=0; ip < npar; ip++){
1304 AliVParticle *part = (AliVParticle*) particleListGen.At(ip);
1305 if(!part) continue;
1306
1307 Double_t pta = 0.9 * part->Pt(); //10% energy loss
1308 Double_t ptb = 0.95 * part->Pt(); //5% energy loss
1309 if(fTriggerPtRangeLow <= pta && pta < fTriggerPtRangeHigh && ntriggersMCa<200){
1310 triggersMCa[ntriggersMCa] = ip;
1311 ntriggersMCa++;
1312 }
1313
1314 if(fTriggerPtRangeLow <= ptb && ptb < fTriggerPtRangeHigh && ntriggersMCb<200){
1315 triggersMCb[ntriggersMCb] = ip;
1316 ntriggersMCb++;
1317 }
1318 }
1319
1320 if(ntriggersMCa>0){
1321 Int_t rnda = fRandom->Integer(ntriggersMCa); //0 to ntriggers-1
1322 Int_t indexTriggGena = triggersMCa[rnda];
1323
1324 Double_t deltaPhia, deltaEtaa, deltaRa;
1325 Int_t aa = 0;
1326
1327 //Correlation with single inclusive TRIGGER
1328 AliVParticle* tGenTa = (AliVParticle*) particleListGen.At(indexTriggGena);
1329 if(tGenTa){
1330 //for(Int_t ia=0; ia<ntriggersMCa; ia++)
1331 for(Int_t ia=0; ia< npar; ia++){
1332 //if(indexTriggGena == triggersMCa[ia]) continue;
1333 if(indexTriggGena == ia) continue;
1334
1335 //AliVParticle* tGenTz = (AliVParticle*) particleListGen.At(triggersMCa[ia]);
1336 AliVParticle* tGenTz = (AliVParticle*) particleListGen.At(ia);
1337 if(!tGenTz) continue;
1338 if(tGenTz->Pt()*0.9<15.0) continue;
1339 if(tGenTz->Pt()*0.9>100.0) continue;
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){
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;
1365
1366 //AliVParticle* tGenTz = (AliVParticle*) particleListGen.At(triggersMCb[ib]);
1367 AliVParticle* tGenTz = (AliVParticle*) particleListGen.At(ib);
1368 if(!tGenTz) continue;
1369 if(tGenTz->Pt()*0.95<15.0) continue;
1370 if(tGenTz->Pt()*0.95>100.0) continue;
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
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
1392 EstimateBgRhoMedian(fListJetsBgGen, &particleListGen, rhoFromCellMedianGen,1);//mc data
1393
1394 //============ Generator trigger+jet ==================
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];
1399
1400 fhNofMultipleTriggersGen->Fill(ntriggersMC-1);
1401
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){
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;
1412
1413 //AliVParticle* tGenT2 = (AliVParticle*) particleListGen.At(triggersMC[ia]);
1414 AliVParticle* tGenT2 = (AliVParticle*) particleListGen.At(ia);
1415 if(!tGenT2) continue;
1416 if(tGenT2->Pt()<15.0) continue;
1417 if(tGenT2->Pt()>100.0) continue;
1418 deltaPhi = RelativePhi(tGenT1->Phi(),tGenT2->Phi());
1419 deltaEta = tGenT1->Eta()-tGenT2->Eta();
1420 deltaR = sqrt(deltaPhi*deltaPhi + deltaEta*deltaEta);
1421
1422 if(deltaR<0.4) k++;
1423 }
1424 }
1425 fhNofMultipleTriggersConeGen->Fill(k);
1426
1427 if(ntriggersMC>1){
1428 //Correlation of each trigger with any other trigger
1429 for(Int_t ia=0; ia<ntriggersMC-1; ia++){
1430 AliVParticle* tGenI = (AliVParticle*) particleListGen.At(triggersMC[ia]);
1431 if(!tGenI) continue;
1432 for(Int_t ib=ia+1; ib<ntriggersMC; ib++){
1433 AliVParticle* tGenII = (AliVParticle*) particleListGen.At(triggersMC[ib]);
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
1444 }else{
1445 indexTriggGen = -1; //trigger not found
1446 }
1447 }
1448
1449 //-----------
1450 Int_t ilowGen = (fHardest==0 || fHardest==1) ? indexTriggGen : 0;
1451 Int_t ihighGen = (fHardest==0 || fHardest==1) ? indexTriggGen+1 : particleListGen.GetEntries();
1452 Bool_t fillOnceGen = kTRUE;
1453 //-----------
1454
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){
1463 if(triggerGen->Pt() < 6.0) continue; //all hadrons pt>6
1464 }
1465 if(TMath::Abs((Float_t) triggerGen->Eta()) > fTriggerEtaCut) continue;
1466
1467 ptTriggGen = triggerGen->Pt(); //count triggers
1468 fh2NtriggersGen->Fill(centValue, ptTriggGen);
1469
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();
1475 Double_t areaJetGen = jet->EffectiveAreaCharged();
1476
1477 if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){
1478
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
1506 Double_t dphi = RelativePhi(triggerGen->Phi(), jet->Phi());
1507 if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue;
1508
1509 //Centrality, A, pT, pTtrigg
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
1517 //Centrality, A, pTjet-pTbgCellMedian, pTtrigg, dphi
1518 tmpArrayFive[0] = centValue;
1519 tmpArrayFive[1] = areaJetGen;
1520 tmpArrayFive[2] = jet->Pt() - ptUeFromCellMedianGen;
1521 tmpArrayFive[3] = ptTriggGen;
1522 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1523 fHJetSpecSubUeMedianGen->Fill(tmpArrayFive);
1524
1525 //Centrality, A, pTjet-pTbgCone, pTtrigg, dphi
1526 tmpArrayFive[0] = centValue;
1527 tmpArrayFive[1] = areaJetGen;
1528 tmpArrayFive[2] = jet->Pt() - ptUeConeGen;
1529 tmpArrayFive[3] = ptTriggGen;
1530 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1531 fHJetSpecSubUeConeGen->Fill(tmpArrayFive);
1532
1533 //Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", kT median
1534 tmpArrayFive[0] = areaJetGen;
1535 tmpArrayFive[1] = jet->Pt();
1536 tmpArrayFive[2] = jet->Pt() - ptUeFromCellMedianGen;
1537 tmpArrayFive[3] = ptUeFromCellMedianGen;
1538 tmpArrayFive[4] = rhoFromCellMedianGen;
1539 fHJetUeMedianGen->Fill(tmpArrayFive);
1540
1541 //Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", perpendicular Cone
1542 tmpArrayFive[0] = areaJetGen;
1543 tmpArrayFive[1] = jet->Pt();
1544 tmpArrayFive[2] = jet->Pt() - ptUeConeGen;
1545 tmpArrayFive[3] = ptUeConeGen;
1546 tmpArrayFive[4] = rhoConeGen;
1547 fHJetUeConeGen->Fill(tmpArrayFive);
1548
1549 if(fillOnceGen){
1550 Double_t fillRhoGen[] = { rhoConeGen, rhoFromCellMedianGen};
1551 fHRhoUeMedianVsConeGen->Fill(fillRhoGen);
1552 fillOnceGen = kFALSE;
1553 }
1554 }//back to back jet-trigger pair
1555 }//jet loop
1556 }//trigger loop
1557
1558 if(fFillRespMx && !fIsKine){
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();
1566
1567 if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){
1568 fhJetPtGen->Fill(ptJetGen); // gen level pt spectum of jets response mx normalization
1569
1570 Double_t areaJetGen = jet->EffectiveAreaCharged();
1571 Double_t ptUeFromCellMedianGen = rhoFromCellMedianGen*areaJetGen;
1572 Double_t ptUeConeGen = rhoConeGen*areaJetGen;
1573
1574 fhJetPtSubUeMedianGen->Fill(ptJetGen - ptUeFromCellMedianGen);
1575 fhJetPtSubUeConeGen->Fill(ptJetGen - ptUeConeGen);
1576 }
1577 }
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();
1582
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
1586
1587 if(fDebug){
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());
1590 }
1591 //matching of MC genrator level and reconstructed jets
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();
1600 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1601
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();
1609
1610 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
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);
1622 }
1623 }//ig>=0
1624 }//rec jet in eta acceptance
1625 }//loop over reconstructed jets
1626 }// # of rec jets >0
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
1684 if(fIsKine){ //skip reconstructed data analysis in case of kine
1685 PostData(1, fOutputList);
1686 return;
1687 }
1688 //============= RECONSTRUCTED INCLUSIVE JETS ===============
1689
1690 Double_t etaJet = 0.0;
1691 Double_t pTJet = 0.0;
1692 Double_t areaJet = 0.0;
1693 Double_t phiJet = 0.0;
1694 //Int_t indexLeadingJet = -1;
1695 //Double_t pTLeadingJet = -10.0;
1696 //Double_t areaLeadingJet = -10.0;
1697
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;
1707 /*areaJet = jet->EffectiveAreaCharged();*/
1708
1709 //Jet Diagnostics---------------------------------
1710 fhJetPhi->Fill(pTJet, RelativePhi(phiJet,0.0)); //phi -pi,pi
1711 fhJetEta->Fill(pTJet, etaJet);
1712 //search for leading jet
1713 /*if(pTJet > pTLeadingJet){
1714 indexLeadingJet = ij;
1715 pTLeadingJet = pTJet;
1716 areaLeadingJet = areaJet;
1717 }
1718
1719 // raw spectra of INCLUSIVE jets
1720 //Centrality, pTjet, A
1721 Double_t fillraw[] = { centValue,
1722 pTJet,
1723 areaJet
1724 };
1725 fHJetPtRaw->Fill(fillraw);*/
1726 }
1727 /*
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 }
1737 */
1738
1739 // =============== RECONSTRUCTED TRIGGER-JET PAIRS ================
1740 if(fIsChargedMC && fFillRespMx){
1741 FillEffHistos(&particleList, &particleListGen); //Fill efficiency histos
1742 }
1743 Bool_t filledOnce = kTRUE; //fill rho histogram only once per event
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();
1747
1748 for(Int_t itrk= ilow; itrk< ihigh; itrk++){ //loop over possible trigger
1749 indexTrigg = itrk; //trigger hadron with pT >10 GeV
1750
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){
1759 if(triggerHadron->Pt() < 6.0) continue; //all hadrons pt>6
1760 }
1761 if(TMath::Abs((Float_t) triggerHadron->Eta())> fTriggerEtaCut) continue;
1762
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());
1767
1768
1769 //---------- make trigger-jet pairs ---------
1770 //Int_t injet4 = 0;
1771 //Int_t injet = 0;
1772
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();
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
1788 //if(areaJet >= 0.07) injet++;
1789 //if(areaJet >= 0.4) injet4++;
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
1796 /*Double_t filldp[] = { centValue,
1797 dphi,
1798 pTJet,
1799 triggerHadron->Pt()
1800 };
1801 fHDphiVsJetPtAll->Fill(filldp);
1802 */
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
1829 // Select back to back trigger - jet pairs
1830 if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue;
1831 fhDphiTriggerJetAccept->Fill(dphi); //Accepted
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
1842 //Centrality, A, pTjet-pTbgCellMedian, pTtrigg, dphi
1843 tmpArrayFive[0] = centValue;
1844 tmpArrayFive[1] = areaJet;
1845 tmpArrayFive[2] = pTJet-ptUeFromCellMedian;
1846 tmpArrayFive[3] = triggerHadron->Pt();
1847 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1848 fHJetSpecSubUeMedian->Fill(tmpArrayFive);
1849
1850 //Centrality, A, pTjet-pTbgCone, pTtrigg, dphi
1851 tmpArrayFive[0] = centValue;
1852 tmpArrayFive[1] = areaJet;
1853 tmpArrayFive[2] = pTJet - ptUeCone;
1854 tmpArrayFive[3] = triggerHadron->Pt();
1855 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1856 fHJetSpecSubUeCone->Fill(tmpArrayFive);
1857
1858 //Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", kT median
1859 tmpArrayFive[0] = areaJet;
1860 tmpArrayFive[1] = pTJet;
1861 tmpArrayFive[2] = pTJet - ptUeFromCellMedian;
1862 tmpArrayFive[3] = ptUeFromCellMedian;
1863 tmpArrayFive[4] = rhoFromCellMedian;
1864 fHJetUeMedian->Fill(tmpArrayFive);
1865
1866 //Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", Cone median
1867 tmpArrayFive[0] = areaJet;
1868 tmpArrayFive[1] = pTJet;
1869 tmpArrayFive[2] = pTJet - ptUeCone;
1870 tmpArrayFive[3] = ptUeCone;
1871 tmpArrayFive[4] = rhoCone;
1872 fHJetUeCone->Fill(tmpArrayFive);
1873
1874 if(filledOnce){ //fill for each event only once
1875 Double_t fillRho[] = { rhoCone,rhoFromCellMedian};
1876 fHRhoUeMedianVsCone->Fill(fillRho);
1877 filledOnce = kFALSE;
1878 }
1879 }//jet loop
1880
1881 //Fill Jet Density In the Event A>0.07
1882 /*if(injet>0){
1883 Double_t filldens[]={ (Double_t) particleList.GetEntries(),
1884 injet/fkAcceptance,
1885 triggerHadron->Pt()
1886 };
1887 fHJetDensity->Fill(filldens);
1888 }
1889
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);
1897 }*/
1898 }
1899
1900
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;
1929 Int_t triggers[200];
1930 Int_t ntriggers = 0; //index in triggers array
1931
1932
1933 for(Int_t it = 0; it < aodevt->GetNumberOfTracks(); it++){
1934 //AliAODTrack *tr = aodevt->GetTrack(it);
1935 AliAODTrack *tr = dynamic_cast<AliAODTrack*>(aodevt->GetTrack(it));
1936 if(!tr) AliFatal("Not a standard AOD");
1937
1938 if((fFilterMask > 0) && !(tr->TestFilterBit(fFilterMask))) continue;
1939 //if((fFilterMask > 0) && !(tr->IsHybridGlobalConstrainedGlobal())) continue;
1940 if(TMath::Abs((Float_t) tr->Eta()) > fTrackEtaCut) continue;
1941 if(tr->Pt() < fTrackLowPtCut) continue;
1942 list->Add(tr);
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 }
1958 }
1959
1960 iCount++;
1961 }
1962
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];
1966
1967 fhNofMultipleTriggers->Fill(ntriggers-1);
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){
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);
1980 if(!tGent2) continue;
1981 if(tGent2->Pt()<15.0) continue;
1982 if(tGent2->Pt()>100.0) continue;
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);
1990
1991 if(ntriggers>1){
1992 //Correlation with any other trigger
1993 for(Int_t ia=0; ia<ntriggers-1; ia++){
1994 AliVParticle* tGeni = (AliVParticle*) list->At(triggers[ia]);
1995 if(!tGeni) continue;
1996 for(Int_t ib=ia+1; ib<ntriggers; ib++){
1997 AliVParticle* tGenii = (AliVParticle*) list->At(triggers[ib]);
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 }
2006 }
2007
2008 return index;
2009}
2010
2011//----------------------------------------------------------------------------
2012/*
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}
2032*/
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
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);
2058 if(fFillRespMx) fhPtTrkTruePrimGen->Fill(trk->Pt(),trk->Eta()); //Efficiency denominator
2059
2060 //Search trigger:
2061 if(fHardest>0){ //leading particle
2062 if(ptLeading < trk->Pt()){
2063 index = counter;
2064 ptLeading = trk->Pt();
2065 }
2066 }
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
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}
2107//________________________________________________________________
2108void AliAnalysisTaskJetCorePP::EstimateBgRhoMedian(TList *listJet, TList* listPart, Double_t &rhoMedian, Int_t mode){
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++;
2144 }
2145 }
2146
2147
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;
2155 }
2156 }
2157
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
2187
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 }
2206 }
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;
2217 Double_t bufferArea=0.0, bufferPt=0.0; //sum cells where A< fJetFreeAreaFrac
2218 for(Int_t ip=0; ip<fnPhi; ip++){
2219 for(Int_t ie=0; ie<fnEta; ie++){
2220 relativeArea = nOutCone[ip][ie]/fnTrials;
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 }
2235 }
2236 }
2237
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 }
2247
2248}
2249
2250//__________________________________________________________________
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//__________________________________________________________________
2298
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;
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 }
2312 if(fJetBranchNameBgKine.Length()>0 && bname.CompareTo(fJetBranchNameBgKine.Data())==0){
2313 array = dynamic_cast<TClonesArray*>(GetInputData(2)); //connect exchange container slot 2
2314 }
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 }
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 }
2333
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}