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