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