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