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