including extra debugging message
[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(""),
80ac66f6 80fJetBranchNameMC(""),
160dbdb8 81fJetBranchNameMCFull(""),
7fc3d134 82fListJets(0x0),
80ac66f6 83fListJetsGen(0x0),
160dbdb8 84fListJetsGenFull(0x0),
ad869500 85fNonStdFile(""),
86fSystem(0), //pp=0 pPb=1
87fJetParamR(0.4),
88fOfflineTrgMask(AliVEvent::kAny),
89fMinContribVtx(1),
90fVtxZMin(-10.0),
91fVtxZMax(10.0),
92fFilterMask(0),
93fCentMin(0.0),
94fCentMax(100.0),
95fJetEtaMin(-0.5),
96fJetEtaMax(0.5),
97fTriggerEtaCut(0.9),
98fTrackEtaCut(0.9),
99fTrackLowPtCut(0.15),
100fOutputList(0x0),
101fHistEvtSelection(0x0),
102fh2Ntriggers(0x0),
103fHJetSpec(0x0),
104fHJetDensity(0x0),
105fHJetDensityA4(0x0),
106fhJetPhi(0x0),
107fhTriggerPhi(0x0),
108fhJetEta(0x0),
109fhTriggerEta(0x0),
110fhVertexZ(0x0),
111fhVertexZAccept(0x0),
112fhContribVtx(0x0),
113fhContribVtxAccept(0x0),
114fhDphiTriggerJet(0x0),
115fhDphiTriggerJetAccept(0x0),
116fhCentrality(0x0),
117fhCentralityAccept(0x0),
7fc3d134 118fHJetPtRaw(0x0),
119fHLeadingJetPtRaw(0x0),
120fHDphiVsJetPtAll(0x0),
80ac66f6 121fhJetPtGenVsJetPtRec(0x0),
160dbdb8 122fhJetPtGen(0x0),
123fhJetPtGenChargVsJetPtGenFull(0x0),
124fhJetPtGenFull(0x0),
80ac66f6 125fh2NtriggersGen(0x0),
126fHJetSpecGen(0x0),
127fhPtTrkTruePrimRec(0x0),
128fhPtTrkTruePrimGen(0x0),
129fhPtTrkSecOrFakeRec(0x0),
130fIsMC(0),
131faGenIndex(0),
132faRecIndex(0),
7fc3d134 133fkAcceptance(2.0*TMath::Pi()*1.8),
6168afc3 134fkDeltaPhiCut(TMath::Pi()-0.8),
ea603b64 135fh1Xsec(0x0),
136fh1Trials(0x0),
137fh1AvgTrials(0x0),
138fh1PtHard(0x0),
139fh1PtHardNoW(0x0),
140fh1PtHardTrials(0x0),
8e103a56 141fAvgTrials(1),
142fHardest(0),
143fEventNumberRangeLow(0),
144fEventNumberRangeHigh(99),
145fTriggerPtRangeLow(0.0),
146fTriggerPtRangeHigh(10000.0),
147fRandom(0x0)
ad869500 148{
149 // default Constructor
ad869500 150}
151
152//---------------------------------------------------------------------
153
154AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP(const char *name) :
155AliAnalysisTaskSE(name),
156fESD(0x0),
157fAODIn(0x0),
158fAODOut(0x0),
159fAODExtension(0x0),
160fJetBranchName(""),
80ac66f6 161fJetBranchNameMC(""),
160dbdb8 162fJetBranchNameMCFull(""),
7fc3d134 163fListJets(0x0),
80ac66f6 164fListJetsGen(0x0),
160dbdb8 165fListJetsGenFull(0x0),
ad869500 166fNonStdFile(""),
167fSystem(0), //pp=0 pPb=1
168fJetParamR(0.4),
169fOfflineTrgMask(AliVEvent::kAny),
170fMinContribVtx(1),
171fVtxZMin(-10.0),
172fVtxZMax(10.0),
173fFilterMask(0),
174fCentMin(0.0),
175fCentMax(100.0),
176fJetEtaMin(-0.5),
177fJetEtaMax(0.5),
178fTriggerEtaCut(0.9),
179fTrackEtaCut(0.9),
180fTrackLowPtCut(0.15),
181fOutputList(0x0),
182fHistEvtSelection(0x0),
183fh2Ntriggers(0x0),
184fHJetSpec(0x0),
185fHJetDensity(0x0),
186fHJetDensityA4(0x0),
187fhJetPhi(0x0),
188fhTriggerPhi(0x0),
189fhJetEta(0x0),
190fhTriggerEta(0x0),
191fhVertexZ(0x0),
192fhVertexZAccept(0x0),
193fhContribVtx(0x0),
194fhContribVtxAccept(0x0),
195fhDphiTriggerJet(0x0),
196fhDphiTriggerJetAccept(0x0),
197fhCentrality(0x0),
198fhCentralityAccept(0x0),
7fc3d134 199fHJetPtRaw(0x0),
200fHLeadingJetPtRaw(0x0),
201fHDphiVsJetPtAll(0x0),
80ac66f6 202fhJetPtGenVsJetPtRec(0x0),
203fhJetPtGen(0x0),
160dbdb8 204fhJetPtGenChargVsJetPtGenFull(0x0),
205fhJetPtGenFull(0x0),
80ac66f6 206fh2NtriggersGen(0x0),
207fHJetSpecGen(0x0),
208fhPtTrkTruePrimRec(0x0),
209fhPtTrkTruePrimGen(0x0),
210fhPtTrkSecOrFakeRec(0x0),
211fIsMC(0),
212faGenIndex(0),
213faRecIndex(0),
7fc3d134 214fkAcceptance(2.0*TMath::Pi()*1.8),
6168afc3 215fkDeltaPhiCut(TMath::Pi()-0.8),
ea603b64 216fh1Xsec(0x0),
217fh1Trials(0x0),
218fh1AvgTrials(0x0),
219fh1PtHard(0x0),
220fh1PtHardNoW(0x0),
221fh1PtHardTrials(0x0),
8e103a56 222fAvgTrials(1),
223fHardest(0),
224fEventNumberRangeLow(0),
225fEventNumberRangeHigh(99),
226fTriggerPtRangeLow(0.0),
227fTriggerPtRangeHigh(10000.0),
228fRandom(0x0)
ad869500 229{
7fc3d134 230// Constructor
ad869500 231
232 DefineOutput(1, TList::Class());
233}
234
235//--------------------------------------------------------------
236AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP(const AliAnalysisTaskJetCorePP& a):
237AliAnalysisTaskSE(a.GetName()),
238fESD(a.fESD),
239fAODIn(a.fAODIn),
240fAODOut(a.fAODOut),
241fAODExtension(a.fAODExtension),
242fJetBranchName(a.fJetBranchName),
80ac66f6 243fJetBranchNameMC(a.fJetBranchNameMC),
160dbdb8 244fJetBranchNameMCFull(a.fJetBranchNameMCFull),
ad869500 245fListJets(a.fListJets),
80ac66f6 246fListJetsGen(a.fListJetsGen),
160dbdb8 247fListJetsGenFull(a.fListJetsGenFull),
ad869500 248fNonStdFile(a.fNonStdFile),
249fSystem(a.fSystem),
250fJetParamR(a.fJetParamR),
251fOfflineTrgMask(a.fOfflineTrgMask),
252fMinContribVtx(a.fMinContribVtx),
253fVtxZMin(a.fVtxZMin),
254fVtxZMax(a.fVtxZMax),
255fFilterMask(a.fFilterMask),
256fCentMin(a.fCentMin),
257fCentMax(a.fCentMax),
258fJetEtaMin(a.fJetEtaMin),
259fJetEtaMax(a.fJetEtaMax),
260fTriggerEtaCut(a.fTriggerEtaCut),
261fTrackEtaCut(a.fTrackEtaCut),
262fTrackLowPtCut(a.fTrackLowPtCut),
263fOutputList(a.fOutputList),
264fHistEvtSelection(a.fHistEvtSelection),
265fh2Ntriggers(a.fh2Ntriggers),
266fHJetSpec(a.fHJetSpec),
267fHJetDensity(a.fHJetDensity),
268fHJetDensityA4(a.fHJetDensityA4),
269fhJetPhi(a.fhJetPhi),
270fhTriggerPhi(a.fhTriggerPhi),
271fhJetEta(a.fhJetEta),
272fhTriggerEta(a.fhTriggerEta),
273fhVertexZ(a.fhVertexZ),
274fhVertexZAccept(a.fhVertexZAccept),
275fhContribVtx(a.fhContribVtx),
276fhContribVtxAccept(a.fhContribVtxAccept),
277fhDphiTriggerJet(a.fhDphiTriggerJet),
278fhDphiTriggerJetAccept(a.fhDphiTriggerJetAccept),
279fhCentrality(a.fhCentrality),
280fhCentralityAccept(a.fhCentralityAccept),
7fc3d134 281fHJetPtRaw(a.fHJetPtRaw),
282fHLeadingJetPtRaw(a.fHLeadingJetPtRaw),
283fHDphiVsJetPtAll(a.fHDphiVsJetPtAll),
80ac66f6 284fhJetPtGenVsJetPtRec(a.fhJetPtGenVsJetPtRec),
285fhJetPtGen(a.fhJetPtGen),
160dbdb8 286fhJetPtGenChargVsJetPtGenFull(a.fhJetPtGenChargVsJetPtGenFull),
287fhJetPtGenFull(a.fhJetPtGenFull),
80ac66f6 288fh2NtriggersGen(a.fh2NtriggersGen),
289fHJetSpecGen(a.fHJetSpecGen),
290fhPtTrkTruePrimRec(a.fhPtTrkTruePrimRec),
291fhPtTrkTruePrimGen(a.fhPtTrkTruePrimGen),
292fhPtTrkSecOrFakeRec(a.fhPtTrkSecOrFakeRec),
293fIsMC(a.fIsMC),
294faGenIndex(a.faGenIndex),
295faRecIndex(a.faRecIndex),
7fc3d134 296fkAcceptance(a.fkAcceptance),
ea603b64 297fkDeltaPhiCut(a.fkDeltaPhiCut),
298fh1Xsec(a.fh1Xsec),
299fh1Trials(a.fh1Trials),
300fh1AvgTrials(a.fh1AvgTrials),
301fh1PtHard(a.fh1PtHard),
302fh1PtHardNoW(a.fh1PtHardNoW),
303fh1PtHardTrials(a.fh1PtHardTrials),
8e103a56 304fAvgTrials(a.fAvgTrials),
305fHardest(a.fHardest),
306fEventNumberRangeLow(a.fEventNumberRangeLow),
307fEventNumberRangeHigh(a.fEventNumberRangeHigh),
308fTriggerPtRangeLow(a.fTriggerPtRangeLow),
43457761 309fTriggerPtRangeHigh(a.fTriggerPtRangeHigh),
310fRandom(a.fRandom)
ad869500 311{
312 //Copy Constructor
43457761 313 fRandom->SetSeed(0);
ad869500 314}
315//--------------------------------------------------------------
316
317AliAnalysisTaskJetCorePP& AliAnalysisTaskJetCorePP::operator = (const AliAnalysisTaskJetCorePP& a){
318 // assignment operator
319 this->~AliAnalysisTaskJetCorePP();
320 new(this) AliAnalysisTaskJetCorePP(a);
321 return *this;
322}
323//--------------------------------------------------------------
324
325AliAnalysisTaskJetCorePP::~AliAnalysisTaskJetCorePP()
326{
327 //Destructor
328 delete fListJets;
80ac66f6 329 delete fListJetsGen;
160dbdb8 330 delete fListJetsGenFull;
ad869500 331 delete fOutputList; // ????
8e103a56 332 delete fRandom;
ad869500 333}
334
335//--------------------------------------------------------------
336
ea603b64 337
338Bool_t AliAnalysisTaskJetCorePP::Notify()
339{
340 //Implemented Notify() to read the cross sections
341 //and number of trials from pyxsec.root
342 //inspired by AliAnalysisTaskJetSpectrum2::Notify()
343 if(!fIsMC) return kFALSE;
344
345 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
346 if(!fESD){
347 if(fDebug>1) AliError("ESD not available");
348 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
349 }
350
351 fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
352
353
354 if(fNonStdFile.Length()!=0){
355 // case that we have an AOD extension we can fetch the jets from the extended output
356 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
357 fAODExtension = aodH ? aodH->GetExtension(fNonStdFile.Data()) : 0;
358 if(!fAODExtension){
359 if(fDebug>1) Printf("AODExtension found for %s",fNonStdFile.Data());
360 }
361 }
362
363 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
364 Float_t xsection = 0;
365 Float_t ftrials = 1;
366
367 fAvgTrials = 1;
368 if(tree){
369 TFile *curfile = tree->GetCurrentFile();
370 if(!curfile) {
371 Error("Notify","No current file");
372 return kFALSE;
373 }
374 if(!fh1Xsec || !fh1Trials){
375 Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
376 return kFALSE;
377 }
378 AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
379 fh1Xsec->Fill("<#sigma>",xsection);
380 // construct a poor man average trials
381 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
382 if(ftrials>=nEntries && nEntries>0.) fAvgTrials = ftrials/nEntries;
383 fh1Trials->Fill("#sum{ntrials}",ftrials);
384 }
385
386 if(fDebug)Printf("Reading File %s",fInputHandler->GetTree()->GetCurrentFile()->GetName());
387
388 return kTRUE;
389}
390//--------------------------------------------------------------
391
ad869500 392void AliAnalysisTaskJetCorePP::Init()
393{
394 // check for jet branches
395 if(!strlen(fJetBranchName.Data())){
396 AliError("Jet branch name not set.");
397 }
ad869500 398}
399
400//--------------------------------------------------------------
401
402void AliAnalysisTaskJetCorePP::UserCreateOutputObjects()
403{
ad869500 404 // Create histograms
405 // Called once
80ac66f6 406 fListJets = new TList(); //reconstructed level
407
408 fIsMC = (fJetBranchNameMC.Length()>0) ? kTRUE : kFALSE;
409
8e103a56 410 fRandom = new TRandom3(0);
411
160dbdb8 412 if(fIsMC){
413 fListJetsGen = new TList(); //generator level charged jets
414 if(fJetBranchNameMCFull.Length()>0)
415 fListJetsGenFull = new TList(); //generator level full jets
416 }
ad869500 417 OpenFile(1);
418 if(!fOutputList) fOutputList = new TList();
419 fOutputList->SetOwner(kTRUE);
420
421 Bool_t oldStatus = TH1::AddDirectoryStatus();
422 TH1::AddDirectory(kFALSE);
423
424 fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
425 fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
426 fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
427 fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
428 fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
429 fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
8e103a56 430 fHistEvtSelection->GetXaxis()->SetBinLabel(6,"event number (rejected)");
ad869500 431
432 fOutputList->Add(fHistEvtSelection);
433
434 Int_t nBinsCentrality = (fSystem==0) ? 1 : 10; // pp=1 else 10
80ac66f6 435 //trigger pt spectrum (reconstructed)
ad869500 436 fh2Ntriggers = new TH2F("fh2Ntriggers","# of triggers",
437 nBinsCentrality,0.0,100.0,50,0.0,50.0);
7fc3d134 438 fOutputList->Add(fh2Ntriggers);
ad869500 439
6168afc3 440 //Centrality, A, pTjet, pTtrigg, dphi
441 const Int_t dimSpec = 5;
160dbdb8 442 const Int_t nBinsSpec[dimSpec] = {nBinsCentrality, 100, 200, 50, TMath::Nint(10*(TMath::Pi()-fkDeltaPhiCut))};
6168afc3 443 const Double_t lowBinSpec[dimSpec] = {0.0, 0.0, 0, 0.0, fkDeltaPhiCut};
444 const Double_t hiBinSpec[dimSpec] = {100.0, 1.0, 200.0,50.0, TMath::Pi()};
7fc3d134 445 fHJetSpec = new THnSparseF("fHJetSpec",
6168afc3 446 "Recoil jet spectrum [cent,A,pTjet,pTtrig,dphi]",
7fc3d134 447 dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
ad869500 448 fOutputList->Add(fHJetSpec);
449
450 //------------------- HISTOS FOR DIAGNOSTIC ----------------------
451 //Jet number density histos [Trk Mult, jet density, pT trigger]
452 const Int_t dimJetDens = 3;
80ac66f6 453 const Int_t nBinsJetDens[dimJetDens] = {100, 100, 10};
ad869500 454 const Double_t lowBinJetDens[dimJetDens] = {0.0, 0.0, 0.0};
80ac66f6 455 const Double_t hiBinJetDens[dimJetDens] = {500.0, 5.0, 50.0};
ad869500 456
457 fHJetDensity = new THnSparseF("fHJetDensity","Jet dens vs trk mult A>0.07",
458 dimJetDens,nBinsJetDens,lowBinJetDens,hiBinJetDens);
459
460 fHJetDensityA4 =new THnSparseF("fHJetDensityA4","Jet dens vs trk mult A>0.4",
461 dimJetDens,nBinsJetDens,lowBinJetDens,hiBinJetDens);
462
ad869500 463 fOutputList->Add(fHJetDensity);
464 fOutputList->Add(fHJetDensityA4);
465
466
467 //inclusive azimuthal and pseudorapidity histograms
7fc3d134 468 fhJetPhi = new TH2D("fhJetPhi","Azim dist jets vs pTjet",
469 50, 0, 100, 50,-TMath::Pi(),TMath::Pi());
470 fhTriggerPhi= new TH2D("fhTriggerPhi","azim dist trig had vs pTtrigg",
6168afc3 471 50, 0, 50, 50,-TMath::Pi(),TMath::Pi());
7fc3d134 472 fhJetEta = new TH2D("fhJetEta","Eta dist jets vs pTjet",
473 50,0, 100, 40,-0.9,0.9);
474 fhTriggerEta = new TH2D("fhTriggerEta","Eta dist trig had vs pTtrigg",
6168afc3 475 50, 0, 50, 40,-0.9,0.9);
7fc3d134 476
ad869500 477 fhVertexZ = new TH1D("fhVertexZ","z vertex",40,-20,20);
478 fhVertexZAccept = new TH1D("fhVertexZAccept","z vertex after cut",40,-20,20);
479 fhContribVtx = new TH1D("fhContribVtx","contrib to vtx",200,0,200);
480 fhContribVtxAccept = new TH1D("fhContribVtxAccept","contrib to vtx after cut",200,0,200);
481 fhDphiTriggerJet = new TH1D("fhDphiTriggerJet","Deltaphi trig-jet",50, -TMath::Pi(),TMath::Pi());
482 fhDphiTriggerJetAccept = new TH1D("fhDphiTriggerJetAccept","Deltaphi trig-jet after cut",50, -TMath::Pi(),TMath::Pi());
483 fhCentrality = new TH1D("fhCentrality","Centrality",20,0,100);
484 fhCentralityAccept = new TH1D("fhCentralityAccept","Centrality after cut",20,0,100);
485
486 fOutputList->Add(fhJetPhi);
487 fOutputList->Add(fhTriggerPhi);
488 fOutputList->Add(fhJetEta);
489 fOutputList->Add(fhTriggerEta);
490 fOutputList->Add(fhVertexZ);
491 fOutputList->Add(fhVertexZAccept);
492 fOutputList->Add(fhContribVtx);
493 fOutputList->Add(fhContribVtxAccept);
494 fOutputList->Add(fhDphiTriggerJet);
495 fOutputList->Add(fhDphiTriggerJetAccept);
496 fOutputList->Add(fhCentrality);
497 fOutputList->Add(fhCentralityAccept);
498
7fc3d134 499 // raw spectra of INCLUSIVE jets
80ac66f6 500 //Centrality, pTjet, A
501 const Int_t dimRaw = 3;
502 const Int_t nBinsRaw[dimRaw] = {nBinsCentrality, 50, 100};
503 const Double_t lowBinRaw[dimRaw] = {0.0, 0.0, 0.0};
504 const Double_t hiBinRaw[dimRaw] = {100.0, 100, 1.0};
7fc3d134 505 fHJetPtRaw = new THnSparseF("fHJetPtRaw",
80ac66f6 506 "Incl. jet spectrum [cent,pTjet,A]",
7fc3d134 507 dimRaw,nBinsRaw,lowBinRaw,hiBinRaw);
508 fOutputList->Add(fHJetPtRaw);
509
510 // raw spectra of LEADING jets
80ac66f6 511 //Centrality, pTjet, A
7fc3d134 512 fHLeadingJetPtRaw = new THnSparseF("fHLeadingJetPtRaw",
80ac66f6 513 "Leading jet spectrum [cent,pTjet,A]",
7fc3d134 514 dimRaw,nBinsRaw,lowBinRaw,hiBinRaw);
515 fOutputList->Add(fHLeadingJetPtRaw);
516
517 // Dphi versus pT jet
518 //Centrality, Dphi=phiTrig-phiJet, pTjet, pTtrigg
519 const Int_t dimDp = 4;
520 const Int_t nBinsDp[dimDp] = {nBinsCentrality, 50, 50, 50};
521 const Double_t lowBinDp[dimDp] = {0.0, -TMath::Pi(), 0.0, 0.0};
522 const Double_t hiBinDp[dimDp] = {100.0, TMath::Pi(), 100.0, 100.0};
523 fHDphiVsJetPtAll = new THnSparseF("fHDphiVsJetPtAll",
524 "Dphi vs jet pT [cent,Dphi,pTjet,pTtrigg]",
525 dimDp,nBinsDp,lowBinDp,hiBinDp);
526 fOutputList->Add(fHDphiVsJetPtAll);
527
7fc3d134 528
80ac66f6 529 //analyze MC generator level
530 if(fIsMC){
160dbdb8 531 fhJetPtGenVsJetPtRec = new TH2D("fhJetPtGenVsJetPtRec","JetPtGenVsJetPtRec", 200,0,200, 200,0,200);
80ac66f6 532 fOutputList->Add(fhJetPtGenVsJetPtRec); //gen MC jet pt spectrum versus reconstructed pt spectrum
533
160dbdb8 534 fhJetPtGen = new TH1D("fhJetPtGen","Jet Pt (MC charged jets Gen)",200,0,200); //MC generator jet pt spectrum
80ac66f6 535 fOutputList->Add(fhJetPtGen);
7fc3d134 536
160dbdb8 537 fhJetPtGenChargVsJetPtGenFull = new TH2D("fhJetPtGenChargVsJetPtGenFull","fhJetPtGenChargVsJetPtGenFull", 200,0,200, 200,0,200);
538 fOutputList->Add(fhJetPtGenChargVsJetPtGenFull); //gen MC jet pt spectrum versus reconstructed pt spectrum
539
540 fhJetPtGenFull = new TH1D("fhJetPtGenFull","Jet Pt (MC Full jets Gen)",200,0,200); //MC generator jet pt spectrum
541 fOutputList->Add(fhJetPtGenFull);
542
80ac66f6 543 fh2NtriggersGen = (TH2F*) fh2Ntriggers->Clone("fh2NtriggersGen");
544 fh2NtriggersGen->SetTitle(Form("%s Gen MC",fh2Ntriggers->GetTitle()));
545 fOutputList->Add(fh2NtriggersGen);
7fc3d134 546
80ac66f6 547 //Centrality, A, pT, pTtrigg
548 fHJetSpecGen = (THnSparseF*) fHJetSpec->Clone("fHJetSpecGen");
549 fHJetSpecGen->SetTitle(Form("%s Gen MC",fHJetSpec->GetTitle()));
550 fOutputList->Add(fHJetSpecGen);
551
552 //track efficiency/contamination histograms eta versus pT
553 fhPtTrkTruePrimRec = new TH2D("fhPtTrkTruePrimRec","PtTrkTruePrimRec",100,0,20,18,-0.9,0.9);
554 fOutputList->Add(fhPtTrkTruePrimRec);
555
556 fhPtTrkTruePrimGen = (TH2D*) fhPtTrkTruePrimRec->Clone("fhPtTrkTruePrimGen");
557 fhPtTrkTruePrimGen->SetTitle("PtTrkTruePrimGen");
558 fOutputList->Add(fhPtTrkTruePrimGen);
559
560 fhPtTrkSecOrFakeRec = (TH2D*) fhPtTrkTruePrimRec->Clone("fhPtTrkSecOrFakeRec");
561 fhPtTrkSecOrFakeRec->SetTitle("PtTrkSecOrFakeRec");
562 fOutputList->Add(fhPtTrkSecOrFakeRec);
563 }
ea603b64 564 //-------------------------------------
565 // pythia histograms
566 const Int_t nBinPt = 150;
567 Double_t binLimitsPt[nBinPt+1];
568 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
569 if(iPt == 0){
570 binLimitsPt[iPt] = -50.0;
571 }else{// 1.0
572 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.;
573 }
574 }
575
576 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
577 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
578 fOutputList->Add(fh1Xsec);
579 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
580 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
581 fOutputList->Add(fh1Trials);
582 fh1AvgTrials = new TH1F("fh1AvgTrials","trials root file",1,0,1);
583 fh1AvgTrials->GetXaxis()->SetBinLabel(1,"#sum{avg ntrials}");
584 fOutputList->Add(fh1AvgTrials);
585 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
586 fOutputList->Add(fh1PtHard);
587 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
588 fOutputList->Add(fh1PtHardNoW);
589 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
590 fOutputList->Add(fh1PtHardTrials);
591
ad869500 592
593 // =========== Switch on Sumw2 for all histos ===========
594 for(Int_t i=0; i<fOutputList->GetEntries(); i++){
595 TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
596 if(h1){
597 h1->Sumw2();
598 continue;
599 }
600 THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
601 if(hn){
602 hn->Sumw2();
603 }
604 }
605 TH1::AddDirectory(oldStatus);
606
607 PostData(1, fOutputList);
608}
609
610//--------------------------------------------------------------------
611
612void AliAnalysisTaskJetCorePP::UserExec(Option_t *)
613{
ad869500 614 //Event loop
ea603b64 615 Double_t eventW = 1.0;
616 Double_t ptHard = 0.0;
617 Double_t nTrials = 1.0; // Trials for MC trigger
618 if(fIsMC) fh1AvgTrials->Fill("#sum{avg ntrials}",fAvgTrials);
ad869500 619
620 if(TMath::Abs((Float_t) fJetParamR)<0.00001){
80ac66f6 621 AliError("Cone radius is set to zero.");
ad869500 622 return;
623 }
624 if(!strlen(fJetBranchName.Data())){
625 AliError("Jet branch name not set.");
626 return;
627 }
628
629 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
630 if(!fESD){
80ac66f6 631 if(fDebug>1) AliError("ESD not available");
ad869500 632 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
633 }
634
635 fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
636 AliAODEvent* aod = NULL;
637 // take all other information from the aod we take the tracks from
638 if(!aod){
639 if(!fESD) aod = fAODIn;
640 else aod = fAODOut;
641 }
642
643 if(fNonStdFile.Length()!=0){
644 // case that we have an AOD extension we can fetch the jets from the extended output
645 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
646 fAODExtension = aodH ? aodH->GetExtension(fNonStdFile.Data()) : 0;
647 if(!fAODExtension){
648 if(fDebug>1) Printf("AODExtension found for %s",fNonStdFile.Data());
649 }
650 }
651
80ac66f6 652 // ----------------- EVENT SELECTION --------------------------
ad869500 653 fHistEvtSelection->Fill(1); // number of events before event selection
654
655 // physics selection
656 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
657 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
658
659 if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
660 if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
661 fHistEvtSelection->Fill(2);
662 PostData(1, fOutputList);
663 return;
664 }
665
666 //check AOD pointer
667 if(!aod){
668 if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
669 fHistEvtSelection->Fill(3);
670 PostData(1, fOutputList);
671 return;
672 }
673
674 // vertex selection
675 AliAODVertex* primVtx = aod->GetPrimaryVertex();
676
677 if(!primVtx){
678 if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
679 fHistEvtSelection->Fill(3);
680 PostData(1, fOutputList);
681 return;
682 }
683
684 Int_t nTracksPrim = primVtx->GetNContributors();
685 Float_t vtxz = primVtx->GetZ();
686 //Input events
687 fhContribVtx->Fill(nTracksPrim);
688 if( nTracksPrim > 0 ) fhVertexZ->Fill(vtxz);
689
690 if((nTracksPrim < fMinContribVtx) ||
691 (vtxz < fVtxZMin) ||
692 (vtxz > fVtxZMax)){
693 if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",
694 (char*)__FILE__,__LINE__,vtxz);
695 fHistEvtSelection->Fill(3);
696 PostData(1, fOutputList);
697 return;
698 }else{
699 //Accepted events
700 fhContribVtxAccept->Fill(nTracksPrim);
701 fhVertexZAccept->Fill(vtxz);
702 }
703 //FK// No event class selection imposed
704 // event class selection (from jet helper task)
705 //Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
706 //if(fDebug) Printf("Event class %d", eventClass);
707 //if(eventClass < fEvtClassMin || eventClass > fEvtClassMax){
708 // fHistEvtSelection->Fill(4);
709 // PostData(1, fOutputList);
710 // return;
711 //}
712
80ac66f6 713 //------------------ CENTRALITY SELECTION ---------------
ad869500 714 AliCentrality *cent = 0x0;
715 Double_t centValue = 0.0;
716 if(fSystem){ //fSystem=0 for pp, fSystem=1 for pPb
717 if(fESD){
718 cent = fESD->GetCentrality();
719 if(cent) centValue = cent->GetCentralityPercentile("V0M");
720 }else{
721 centValue = aod->GetHeader()->GetCentrality();
722 }
723 if(fDebug) printf("centrality: %f\n", centValue);
724 //Input events
725 fhCentrality->Fill(centValue);
726
727 if(centValue < fCentMin || centValue > fCentMax){
728 fHistEvtSelection->Fill(4);
729 PostData(1, fOutputList);
730 return;
731 }else{
732 //Accepted events
733 fhCentralityAccept->Fill( centValue );
734 }
735 }
736
8e103a56 737 //-----------------select disjunct event subsamples ----------------
738 Int_t eventnum = aod->GetHeader()->GetEventNumberESDFile();
739 Int_t lastdigit = eventnum % 10;
740 if(!(fEventNumberRangeLow<=lastdigit && lastdigit<=fEventNumberRangeHigh)){
741 fHistEvtSelection->Fill(5);
742 PostData(1, fOutputList);
743 return;
744 }
745
ad869500 746 if(fDebug) std::cout<<" ACCEPTED EVENT "<<endl;
ad869500 747 fHistEvtSelection->Fill(0); // accepted events
ad869500 748 // ------------------- end event selection --------------------
80ac66f6 749
ad869500 750
80ac66f6 751 // fetch RECONSTRUCTED jets
ad869500 752 TClonesArray *aodJets = 0x0;
753
754 if(fAODOut && !aodJets){
755 aodJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName.Data()));
756 }
757 if(fAODExtension && !aodJets){
758 aodJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName.Data()));
759 }
760 if(fAODIn && !aodJets){
761 aodJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName.Data()));
762 }
763
80ac66f6 764 //--------- Fill list of RECONSTRUCTED jets --------------
765 fListJets->Clear();
766 if(aodJets){
767 if(fDebug) Printf("########## %s: %d jets",fJetBranchName.Data(), aodJets->GetEntriesFast());
768 for(Int_t iJet = 0; iJet < aodJets->GetEntriesFast(); iJet++) {
769 AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets)[iJet]);
770 if (jet) fListJets->Add(jet);
771 }
772 }
773
774 //--------- Fill list of MC GENERATOR LEVEL jets --------------
775 TList particleListGen; //list of tracks in MC
776
777 if(fIsMC){ //analyze generator level MC
778 // fetch MC generator level jets
779 TClonesArray *aodGenJets = NULL;
160dbdb8 780 TClonesArray *aodGenJetsFull = NULL;
781
782 //charged jets
80ac66f6 783 if(fAODOut&&!aodGenJets){
784 aodGenJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchNameMC.Data()));
785 }
786 if(fAODExtension&&!aodGenJets){
787 aodGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchNameMC.Data()));
788 }
789 if(fAODIn&&!aodGenJets){
790 aodGenJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchNameMC.Data()));
791 }
792
793 if(!aodGenJets){
794 Printf("%s:%d no generated Jet array with name %s in AOD",
795 (char*)__FILE__,__LINE__, fJetBranchNameMC.Data());
796 PostData(1, fOutputList);
797 return;
798 }
799
800 fListJetsGen->Clear();
801
160dbdb8 802 //full jets
803 if(fJetBranchNameMCFull.Length()>0){
804
805 if(fAODOut&&!aodGenJetsFull){
806 aodGenJetsFull = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchNameMCFull.Data()));
807 }
808 if(fAODExtension&&!aodGenJetsFull){
809 aodGenJetsFull = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchNameMCFull.Data()));
810 }
811 if(fAODIn&&!aodGenJetsFull){
812 aodGenJetsFull = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchNameMCFull.Data()));
813 }
814
815 if(!aodGenJetsFull){
816 Printf("%s:%d no generated Jet array with name %s in AOD",
817 (char*)__FILE__,__LINE__, fJetBranchNameMCFull.Data());
818 PostData(1, fOutputList);
819 return;
820 }
821
822 fListJetsGenFull->Clear();
823 }
824
825
80ac66f6 826 //serarch for charged trigger at the MC generator level
8e103a56 827 Int_t indexTriggGen = -1;
828 Double_t ptTriggGen = -1;
829 Int_t iCounterGen = 0;
830 Int_t triggersMC[200];//list of trigger candidates
831 Int_t ntriggersMC = 0; //index in triggers array
832
80ac66f6 833 if(fESD){//ESD input
834
835 AliMCEvent* mcEvent = MCEvent();
836 if(!mcEvent){
837 PostData(1, fOutputList);
838 return;
ea603b64 839 }
840
841 AliGenPythiaEventHeader *pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
842 if(pythiaGenHeader){
843 nTrials = pythiaGenHeader->Trials();
844 ptHard = pythiaGenHeader->GetPtHard();
845 fh1PtHard->Fill(ptHard,eventW);
846 fh1PtHardNoW->Fill(ptHard,1);
847 fh1PtHardTrials->Fill(ptHard,nTrials);
848 }
849
80ac66f6 850 for(Int_t it = 0; it < mcEvent->GetNumberOfTracks(); it++){
851 if(!mcEvent->IsPhysicalPrimary(it)) continue;
852 AliMCParticle* part = (AliMCParticle*) mcEvent->GetTrack(it);
853 if(!part) continue;
854 if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){
8e103a56 855
856 if(fHardest==0 && ntriggersMC<200){//single inclusive trigger
857 if(indexTriggGen > -1){//trigger candidater was found
858 triggersMC[ntriggersMC] = indexTriggGen;
859 ntriggersMC++;
860 }
861 }
862
80ac66f6 863 iCounterGen++;
864 }
865 }
866 }else{ //AOD input
867 if(!fAODIn){
868 PostData(1, fOutputList);
869 return;
870 }
871 TClonesArray *tca = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(AliAODMCParticle::StdBranchName()));
872 if(!tca){
873 PostData(1, fOutputList);
874 return;
875 }
876 for(Int_t it = 0; it < tca->GetEntriesFast(); it++){
877 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
878 if(!part) continue;
879 if(!part->IsPhysicalPrimary()) continue;
880 if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){
881
8e103a56 882 if(fHardest==0 && ntriggersMC<200){//single inclusive trigger
883 if(indexTriggGen > -1){ //trigger candidater was found
884 triggersMC[ntriggersMC] = indexTriggGen;
885 ntriggersMC++;
886 }
887 }
888
80ac66f6 889 iCounterGen++;
890 }
891 }
892 }
8e103a56 893
894 if(fHardest==0){ //single inclusive trigger
895 if(ntriggersMC>0){ //there is at least one trigger
896 Int_t rnd = fRandom->Integer(ntriggersMC); //0 to ntriggers-1
897 indexTriggGen = triggersMC[rnd];
80ac66f6 898 }else{
8e103a56 899 indexTriggGen = -1; //trigger not found
80ac66f6 900 }
901 }
8e103a56 902
903 //================== Fill jet list ===================
160dbdb8 904 if(aodGenJets){ //charged jets
80ac66f6 905 if(fDebug) Printf("########## %s: %d jets",fJetBranchNameMC.Data(), aodGenJets->GetEntriesFast());
906
907 for(Int_t igJet = 0; igJet < aodGenJets->GetEntriesFast(); igJet++) {
908 AliAODJet *jetGen = dynamic_cast<AliAODJet*>((*aodGenJets)[igJet]);
909 if(jetGen) fListJetsGen->Add(jetGen);
910 }
8e103a56 911 }
912
160dbdb8 913 if(aodGenJetsFull){ //full jets
914 if(fDebug) Printf("########## %s: %d jets",fJetBranchNameMCFull.Data(), aodGenJetsFull->GetEntriesFast());
915
916 for(Int_t igJet = 0; igJet < aodGenJetsFull->GetEntriesFast(); igJet++) {
917 AliAODJet *jetGen = dynamic_cast<AliAODJet*>((*aodGenJetsFull)[igJet]);
918 if(jetGen) fListJetsGenFull->Add(jetGen);
919 }
920 }
921
922
8e103a56 923 //============ Generator trigger+jet ==================
924 Int_t ilowGen = (fHardest==0 || fHardest==1) ? indexTriggGen : 0;
925 Int_t ihighGen = (fHardest==0 || fHardest==1) ? indexTriggGen+1 : particleListGen.GetEntries();
926
927 for(Int_t igen= ilowGen; igen< ihighGen; igen++){ //loop over possible trigger
928 indexTriggGen = igen; //trigger hadron
929
930 if(indexTriggGen == -1) continue;
931 AliVParticle* triggerGen=(AliVParticle*) particleListGen.At(indexTriggGen);
932 if(!triggerGen) continue;
933
934 if(fHardest >= 2){
6168afc3 935 if(triggerGen->Pt() < 6.0) continue; //all hadrons pt>6
8e103a56 936 }
937 if(TMath::Abs((Float_t) triggerGen->Eta()) > fTriggerEtaCut) continue;
938
939 ptTriggGen = triggerGen->Pt(); //count triggers
940 fh2NtriggersGen->Fill(centValue, ptTriggGen);
941
80ac66f6 942 //Count jets and trigger-jet pairs at MC generator level
8e103a56 943 if(!aodGenJets) continue;
80ac66f6 944 for(Int_t ij=0; ij<fListJetsGen->GetEntries(); ij++){
945 AliAODJet* jet = (AliAODJet*)(fListJetsGen->At(ij));
946 if(!jet) continue;
947 Double_t etaJetGen = jet->Eta();
80ac66f6 948
949 if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){
80ac66f6 950
8e103a56 951 Double_t dphi = RelativePhi(triggerGen->Phi(), jet->Phi());
952 if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue;
80ac66f6 953
8e103a56 954 //Centrality, A, pT, pTtrigg
955 Double_t fillspecgen[] = { centValue,
80ac66f6 956 jet->EffectiveAreaCharged(),
957 jet->Pt(),
6168afc3 958 ptTriggGen,
959 TMath::Abs((Double_t) dphi)
80ac66f6 960 };
8e103a56 961 fHJetSpecGen->Fill(fillspecgen);
962 }//back to back jet-trigger pair
963 }//jet loop
964 }//trigger loop
965
966
967 //================ RESPONSE MATRIX ===============
968 if(aodGenJets){
969 //Count jets and trigger-jet pairs at MC generator level
970 for(Int_t ij=0; ij<fListJetsGen->GetEntries(); ij++){
971 AliAODJet* jet = (AliAODJet*)(fListJetsGen->At(ij));
972 if(!jet) continue;
973 Double_t etaJetGen = jet->Eta();
974 Double_t ptJetGen = jet->Pt();
975
976 if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){
977 fhJetPtGen->Fill(ptJetGen); // generator level pt spectum of jets response mx normalization
80ac66f6 978 }
979 }
80ac66f6 980 if(fListJets->GetEntries()>0 && fListJetsGen->GetEntries()>0){ //at least some reconstructed jets
981 Int_t ng = (Int_t) fListJetsGen->GetEntries();
982 Int_t nr = (Int_t) fListJets->GetEntries();
983
984 //Find closest MC generator - reconstructed jets
985 if(faGenIndex.GetSize()<nr) faGenIndex.Set(nr); //idx of gen jet assoc to rec jet
986 if(faRecIndex.GetSize()<ng) faRecIndex.Set(ng); //idx of rec jet assoc to gen jet
987
988 if(fDebug){
989 Printf("New Rec List %d gen index Array %d",nr,faGenIndex.GetSize());
990 Printf("New Gen List %d rec index Array %d",ng,faRecIndex.GetSize());
991 }
992 //matching of MC genrator level and reconstructed jets
993 AliAnalysisHelperJetTasks::GetClosestJets(fListJetsGen,ng,fListJets,nr,faGenIndex,faRecIndex,fDebug);
994
995 // Fill response matrix
996 for(Int_t ir = 0; ir < nr; ir++){
997 AliAODJet *recJet = (AliAODJet*) fListJets->At(ir);
998 Double_t etaJetRec = recJet->Eta();
999 Double_t ptJetRec = recJet->Pt();
1000 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1001
1002 if((fJetEtaMin <= etaJetRec) && (etaJetRec <= fJetEtaMax)){
1003 Int_t ig = faGenIndex[ir]; //associated generator level jet
1004 if(ig >= 0 && ig < ng){
1005 if(fDebug > 10) Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1006 AliAODJet *genJet = (AliAODJet*) fListJetsGen->At(ig);
1007 Double_t ptJetGen = genJet->Pt();
1008 Double_t etaJetGen = genJet->Eta();
1009
1010 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1011 if((fJetEtaMin <= etaJetGen) && (etaJetGen <= fJetEtaMax)){
1012 fhJetPtGenVsJetPtRec->Fill(ptJetRec, ptJetGen);
1013 }
1014 }//ig>=0
1015 }//rec jet in eta acceptance
1016 }//loop over reconstructed jets
1017 }// # of rec jets >0
1018 }//pointer MC generator jets
160dbdb8 1019
1020 //=========================== Full jet charged jet matrix ==========================
1021 if(aodGenJetsFull){
1022 //Count full jets and charged-jet pairs at MC generator level
1023 for(Int_t ij=0; ij<fListJetsGenFull->GetEntries(); ij++){
1024 AliAODJet* jetFull = (AliAODJet*)(fListJetsGenFull->At(ij));
1025 if(!jetFull) continue;
1026 Double_t etaJetFull = jetFull->Eta();
1027 Double_t ptJetFull = jetFull->Pt();
1028
1029 if((fJetEtaMin<=etaJetFull) && (etaJetFull<=fJetEtaMax)){
1030 fhJetPtGenFull->Fill(ptJetFull); // generator level pt spectum of full jets
1031 }
1032 }
1033 if(fListJetsGen->GetEntries()>0 && fListJetsGenFull->GetEntries()>0){ //at least some reconstructed jets
1034 Int_t nful = (Int_t) fListJetsGenFull->GetEntries();
1035 Int_t nchr = (Int_t) fListJetsGen->GetEntries();
1036
1037 //Find closest MC generator full - charged jet
1038 if(faGenIndex.GetSize()<nchr) faGenIndex.Set(nchr); //idx of gen FULL jet assoc to gen CHARGED jet
1039 if(faRecIndex.GetSize()<nful) faRecIndex.Set(nful); //idx of gen CHARGED jet assoc to gen FULL jet
1040
1041 if(fDebug){
1042 Printf("New Rec List %d gen index Array %d",nchr,faGenIndex.GetSize());
1043 Printf("New Gen List %d rec index Array %d",nful,faRecIndex.GetSize());
1044 }
1045 //matching of MC genrator level and reconstructed jets
1046 AliAnalysisHelperJetTasks::GetClosestJets(fListJetsGenFull,nful,fListJetsGen,nchr,faGenIndex,faRecIndex,fDebug);
1047
1048 // Fill response matrix
1049 for(Int_t ichr = 0; ichr < nchr; ichr++){ //charged jet loop
1050 AliAODJet *chJet = (AliAODJet*) fListJetsGen->At(ichr);
1051 Double_t etaJetCh = chJet->Eta();
1052 Double_t ptJetCh = chJet->Pt();
1053 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1054
1055 if((fJetEtaMin <= etaJetCh) && (etaJetCh <= fJetEtaMax)){
1056 Int_t iful = faGenIndex[ichr]; //associated generator level jet
1057 if(iful >= 0 && iful < nful){
1058 if(fDebug > 10) Printf("%s:%d iful = %d ichr = %d",(char*)__FILE__,__LINE__,iful,ichr);
1059 AliAODJet *genJetFull = (AliAODJet*) fListJetsGenFull->At(iful);
1060 Double_t ptJetFull = genJetFull->Pt();
1061 Double_t etaJetFull = genJetFull->Eta();
1062
1063 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1064 if((fJetEtaMin <= etaJetFull) && (etaJetFull <= fJetEtaMax)){
1065 fhJetPtGenChargVsJetPtGenFull->Fill(ptJetFull,ptJetCh);
1066 }
1067 }//iful>=0
1068 }//rec jet in eta acceptance
1069 }//loop over reconstructed jets
1070 }// # of rec jets >0
1071 }//pointer MC generator jets
1072
1073
1074
80ac66f6 1075 } //analyze generator level MC
1076
8e103a56 1077 //============= RECONSTRUCTED INCLUSIVE JETS ===============
ad869500 1078
ad869500 1079 Double_t etaJet = 0.0;
1080 Double_t pTJet = 0.0;
1081 Double_t areaJet = 0.0;
1082 Double_t phiJet = 0.0;
7fc3d134 1083 Int_t indexLeadingJet = -1;
1084 Double_t pTLeadingJet = -10.0;
7fc3d134 1085 Double_t areaLeadingJet = -10.0;
8e103a56 1086
ad869500 1087 for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
1088 AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
1089 if(!jet) continue;
1090 etaJet = jet->Eta();
1091 phiJet = jet->Phi();
1092 pTJet = jet->Pt();
1093 if(pTJet==0) continue;
1094
1095 if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
1096 areaJet = jet->EffectiveAreaCharged();
7fc3d134 1097
ad869500 1098 //Jet Diagnostics---------------------------------
8e103a56 1099 fhJetPhi->Fill(pTJet, RelativePhi(phiJet,0.0)); //phi -pi,pi
1100 fhJetEta->Fill(pTJet, etaJet);
7fc3d134 1101 //search for leading jet
1102 if(pTJet > pTLeadingJet){
1103 indexLeadingJet = ij;
1104 pTLeadingJet = pTJet;
7fc3d134 1105 areaLeadingJet = areaJet;
1106 }
1107
1108 // raw spectra of INCLUSIVE jets
80ac66f6 1109 //Centrality, pTjet, A
7fc3d134 1110 Double_t fillraw[] = { centValue,
1111 pTJet,
7fc3d134 1112 areaJet
1113 };
1114 fHJetPtRaw->Fill(fillraw);
8e103a56 1115 }
7fc3d134 1116
8e103a56 1117 if(indexLeadingJet > -1){
1118 // raw spectra of LEADING jets
1119 //Centrality, pTjet, A
1120 Double_t fillleading[] = { centValue,
1121 pTLeadingJet,
1122 areaLeadingJet
1123 };
1124 fHLeadingJetPtRaw->Fill(fillleading);
1125 }
1126
1127
1128 // =============== RECONSTRUCTED TRIGGER-JET PAIRS ================
1129 //Find Hadron trigger
1130 TList particleList; //list of tracks
1131 Int_t indexTrigg = GetListOfTracks(&particleList); //index of trigger hadron in Particle list
1132 if(fIsMC) FillEffHistos(&particleList, &particleListGen); //Fill efficiency histos
7fc3d134 1133
8e103a56 1134 //set ranges of the trigger loop
1135 Int_t ilow = (fHardest==0 || fHardest==1) ? indexTrigg : 0;
1136 Int_t ihigh = (fHardest==0 || fHardest==1) ? indexTrigg+1 : particleList.GetEntries();
7fc3d134 1137
8e103a56 1138 for(Int_t itrk= ilow; itrk< ihigh; itrk++){ //loop over possible trigger
1139 indexTrigg = itrk; //trigger hadron with pT >10 GeV
7fc3d134 1140
8e103a56 1141 if(indexTrigg < 0) continue;
1142
1143 AliVParticle *triggerHadron = (AliVParticle*) particleList.At(indexTrigg);
1144 if(!triggerHadron){
1145 PostData(1, fOutputList);
1146 return;
1147 }
1148 if(fHardest >= 2){
7d9fd386 1149 if(triggerHadron->Pt() < 6.0) continue; //all hadrons pt>6
8e103a56 1150 }
1151 if(TMath::Abs((Float_t) triggerHadron->Eta())> fTriggerEtaCut) continue;
7fc3d134 1152
8e103a56 1153 //Fill trigger histograms
1154 fh2Ntriggers->Fill(centValue,triggerHadron->Pt()); //trigger pT
1155 fhTriggerPhi->Fill(triggerHadron->Pt(),RelativePhi(triggerHadron->Phi(),0.0)); //phi -pi,pi
1156 fhTriggerEta->Fill(triggerHadron->Pt(),triggerHadron->Eta());
7fc3d134 1157
1158
8e103a56 1159 //---------- make trigger-jet pairs ---------
1160 Int_t injet4 = 0;
1161 Int_t injet = 0;
1162
1163 for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
1164 AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
1165 if(!jet) continue;
1166 etaJet = jet->Eta();
1167 phiJet = jet->Phi();
1168 pTJet = jet->Pt();
1169 if(pTJet==0) continue;
1170
1171 if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
1172 areaJet = jet->EffectiveAreaCharged();
1173 if(areaJet >= 0.07) injet++;
1174 if(areaJet >= 0.4) injet4++;
1175
1176 Double_t dphi = RelativePhi(triggerHadron->Phi(), phiJet);
1177 fhDphiTriggerJet->Fill(dphi); //Input
1178
1179 //Dphi versus jet pT
1180 //Centrality, Dphi=phiTrig-phiJet, pTjet, pTtrigg
1181 Double_t filldp[] = { centValue,
1182 dphi,
1183 pTJet,
1184 triggerHadron->Pt()
1185 };
1186 fHDphiVsJetPtAll->Fill(filldp);
1187
1188 // Select back to back trigger - jet pairs
1189 if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue;
1190 fhDphiTriggerJetAccept->Fill(dphi); //Accepted
1191
1192
1193 //Centrality, A, pTjet, pTtrigg
1194 Double_t fillspec[] = { centValue,
1195 areaJet,
1196 pTJet,
6168afc3 1197 triggerHadron->Pt(),
1198 TMath::Abs((Double_t) dphi)
8e103a56 1199 };
1200 fHJetSpec->Fill(fillspec);
1201 }//jet loop
1202
1203 //Fill Jet Density In the Event A>0.07
1204 if(injet>0){
1205 Double_t filldens[]={ (Double_t) particleList.GetEntries(),
ad869500 1206 injet/fkAcceptance,
1207 triggerHadron->Pt()
1208 };
8e103a56 1209 fHJetDensity->Fill(filldens);
1210 }
ad869500 1211
8e103a56 1212 //Fill Jet Density In the Event A>0.4
1213 if(injet4>0){
1214 Double_t filldens4[]={ (Double_t) particleList.GetEntries(),
1215 injet4/fkAcceptance,
1216 triggerHadron->Pt()
1217 };
1218 fHJetDensityA4->Fill(filldens4);
1219 }
ad869500 1220 }
1221
8e103a56 1222
ad869500 1223 PostData(1, fOutputList);
1224}
1225
1226//----------------------------------------------------------------------------
1227void AliAnalysisTaskJetCorePP::Terminate(const Option_t *)
1228{
1229 // Draw result to the screen
1230 // Called once at the end of the query
1231
1232 if(fDebug) printf("AliAnalysisTaskJetCorePP DONE\n");
1233 if(!GetOutputData(1)) return;
1234}
1235
1236
1237//----------------------------------------------------------------------------
1238Int_t AliAnalysisTaskJetCorePP::GetListOfTracks(TList *list){
1239 //Fill the list of accepted tracks (passed track cut)
1240 //return consecutive index of the hardest ch hadron in the list
1241 Int_t iCount = 0;
1242 AliAODEvent *aodevt = NULL;
1243
1244 if(!fESD) aodevt = fAODIn;
1245 else aodevt = fAODOut;
1246
1247 if(!aodevt) return -1;
1248
1249 Int_t index = -1; //index of the highest particle in the list
1250 Double_t ptmax = -10;
8e103a56 1251 Int_t triggers[200];
1252 Int_t ntriggers = 0; //index in triggers array
ad869500 1253
80ac66f6 1254 for(Int_t it = 0; it < aodevt->GetNumberOfTracks(); it++){
ad869500 1255 AliAODTrack *tr = aodevt->GetTrack(it);
1256
1257 //if((fFilterMask > 0) && !(tr->TestFilterBit(fFilterMask))) continue;
1258 if((fFilterMask > 0) && !(tr->IsHybridGlobalConstrainedGlobal())) continue;
1259 if(TMath::Abs((Float_t) tr->Eta()) > fTrackEtaCut) continue;
1260 if(tr->Pt() < fTrackLowPtCut) continue;
1261 list->Add(tr);
8e103a56 1262
1263 //Search trigger:
1264 if(fHardest>0){ //leading particle
1265 if(tr->Pt()>ptmax){
1266 ptmax = tr->Pt();
1267 index = iCount;
1268 }
1269 }
1270
1271 if(fHardest==0 && ntriggers<200){ //single inclusive
1272 if(fTriggerPtRangeLow <= tr->Pt() &&
1273 tr->Pt() < fTriggerPtRangeHigh){
1274 triggers[ntriggers] = iCount;
1275 ntriggers++;
1276 }
ad869500 1277 }
8e103a56 1278
ad869500 1279 iCount++;
1280 }
1281
8e103a56 1282 if(fHardest==0 && ntriggers>0){ //select random inclusive trigger
1283 Int_t rnd = fRandom->Integer(ntriggers); //0 to ntriggers-1
1284 index = triggers[rnd];
ad869500 1285 }
8e103a56 1286
1287 return index;
ad869500 1288}
1289
1290//----------------------------------------------------------------------------
80ac66f6 1291/*
ad869500 1292Double_t AliAnalysisTaskJetCorePP::GetBackgroundInPerpCone(Float_t jetR, Double_t jetPhi, Double_t jetEta, TList* trkList){
1293 //calculate sum of track pT in the cone perpendicular in phi to the jet
1294 //jetR = cone radius
1295 // jetPhi, jetEta = direction of the jet
1296 Int_t numberOfTrks = trkList->GetEntries();
1297 Double_t pTsum = 0.0;
1298 Double_t perpConePhi = jetPhi + TMath::Pi()/2;//perp cone w.r.t. jet in phi
1299 for(Int_t it=0; it<numberOfTrks; it++){
1300 AliVParticle *trk = (AliVParticle*) trkList->At(it);
1301 Double_t dphi = RelativePhi(perpConePhi,trk->Phi());
1302 Double_t deta = trk->Eta()-jetEta;
1303
1304 if( (dphi*dphi + deta*deta)< (jetR*jetR)){
1305 pTsum += trk->Pt();
1306 }
1307 }
1308
1309 return pTsum;
1310}
80ac66f6 1311*/
ad869500 1312//----------------------------------------------------------------------------
1313
1314Double_t AliAnalysisTaskJetCorePP::RelativePhi(Double_t mphi,Double_t vphi){
1315 //Get relative azimuthal angle of two particles -pi to pi
1316 if (vphi < -TMath::Pi()) vphi += TMath::TwoPi();
1317 else if (vphi > TMath::Pi()) vphi -= TMath::TwoPi();
1318
1319 if (mphi < -TMath::Pi()) mphi += TMath::TwoPi();
1320 else if (mphi > TMath::Pi()) mphi -= TMath::TwoPi();
1321
1322 Double_t dphi = mphi - vphi;
1323 if (dphi < -TMath::Pi()) dphi += TMath::TwoPi();
1324 else if (dphi > TMath::Pi()) dphi -= TMath::TwoPi();
1325
1326 return dphi;//dphi in [-Pi, Pi]
1327}
1328
1329
80ac66f6 1330//----------------------------------------------------------------------------
1331Bool_t AliAnalysisTaskJetCorePP::SelectMCGenTracks(AliVParticle *trk, TList *trkList, Double_t &ptLeading, Int_t &index, Int_t counter){
1332 //fill track efficiency denominator
1333 if(trk->Pt() < fTrackLowPtCut) return kFALSE;
1334 if(trk->Charge()==0) return kFALSE;
1335 if(TMath::Abs((Float_t) trk->Eta()) > fTrackEtaCut) return kFALSE;
1336 trkList->Add(trk);
1337 fhPtTrkTruePrimGen->Fill(trk->Pt(),trk->Eta()); //Efficiency denominator
8e103a56 1338
1339 //Search trigger:
1340 if(fHardest>0){ //leading particle
1341 if(ptLeading < trk->Pt()){
1342 index = counter;
1343 ptLeading = trk->Pt();
1344 }
80ac66f6 1345 }
8e103a56 1346
1347 if(fHardest==0){ //single inclusive
1348 index = -1;
1349 if(fTriggerPtRangeLow <= trk->Pt() &&
1350 trk->Pt() < fTriggerPtRangeHigh){
1351 index = counter;
1352 }
1353 }
1354
80ac66f6 1355 return kTRUE;
1356}
1357
1358//----------------------------------------------------------------------------
1359void AliAnalysisTaskJetCorePP::FillEffHistos(TList *recList, TList *genList){
1360 //fill track efficiency numerator
1361 if(!recList) return;
1362 if(!genList) return;
1363 Int_t nRec = recList->GetEntries();
1364 Int_t nGen = genList->GetEntries();
1365 for(Int_t ir=0; ir<nRec; ir++){
1366 AliAODTrack *trkRec = (AliAODTrack*) recList->At(ir);
1367 if(!trkRec) continue;
1368 Int_t recLabel = TMath::Abs(trkRec->GetLabel());
1369 Bool_t matched = kFALSE;
1370
1371 for(Int_t ig=0; ig<nGen; ig++){
1372 AliVParticle *trkGen = (AliVParticle*) genList->At(ig);
1373 if(!trkGen) continue;
1374 Int_t genLabel = TMath::Abs(trkGen->GetLabel());
1375 if(recLabel==genLabel){
1376 fhPtTrkTruePrimRec->Fill(trkGen->Pt(),trkGen->Eta());
1377 matched = kTRUE;
1378 break;
1379 }
1380 }
1381 if(!matched) fhPtTrkSecOrFakeRec->Fill(trkRec->Pt(),trkRec->Eta());
1382 }
1383
1384 return;
1385}
1386
1387