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