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