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 // *******************************************
11 /**************************************************************************
12 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
14 * Author: The ALICE Off-line Project. *
15 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
35 #include "THnSparse.h"
45 #include "AliAnalysisTask.h"
46 #include "AliAnalysisManager.h"
48 #include "AliVEvent.h"
49 #include "AliESDEvent.h"
50 #include "AliMCEvent.h"
51 #include "AliESDInputHandler.h"
52 #include "AliCentrality.h"
53 #include "AliAnalysisHelperJetTasks.h"
54 #include "AliInputEventHandler.h"
55 #include "AliAODJetEventBackground.h"
56 #include "AliGenPythiaEventHeader.h"
57 #include "AliAODMCParticle.h"
58 #include "AliMCParticle.h"
59 #include "AliAODEvent.h"
60 #include "AliAODHandler.h"
61 #include "AliAODJet.h"
63 #include "AliAnalysisTaskJetCorePP.h"
68 ClassImp(AliAnalysisTaskJetCorePP)
70 //Filip Krizek 1st March 2013
72 //---------------------------------------------------------------------
73 AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP() :
84 fSystem(0), //pp=0 pPb=1
86 fOfflineTrgMask(AliVEvent::kAny),
99 fHistEvtSelection(0x0),
109 fhVertexZAccept(0x0),
111 fhContribVtxAccept(0x0),
112 fhDphiTriggerJet(0x0),
113 fhDphiTriggerJetAccept(0x0),
115 fhCentralityAccept(0x0),
117 fHLeadingJetPtRaw(0x0),
118 fHDphiVsJetPtAll(0x0),
119 fhJetPtGenVsJetPtRec(0x0),
121 fh2NtriggersGen(0x0),
123 fhPtTrkTruePrimRec(0x0),
124 fhPtTrkTruePrimGen(0x0),
125 fhPtTrkSecOrFakeRec(0x0),
129 fkAcceptance(2.0*TMath::Pi()*1.8),
130 fkDeltaPhiCut(TMath::Pi()-0.6),
136 fh1PtHardTrials(0x0),
139 fEventNumberRangeLow(0),
140 fEventNumberRangeHigh(99),
141 fTriggerPtRangeLow(0.0),
142 fTriggerPtRangeHigh(10000.0),
145 // default Constructor
148 //---------------------------------------------------------------------
150 AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP(const char *name) :
151 AliAnalysisTaskSE(name),
157 fJetBranchNameMC(""),
161 fSystem(0), //pp=0 pPb=1
163 fOfflineTrgMask(AliVEvent::kAny),
174 fTrackLowPtCut(0.15),
176 fHistEvtSelection(0x0),
186 fhVertexZAccept(0x0),
188 fhContribVtxAccept(0x0),
189 fhDphiTriggerJet(0x0),
190 fhDphiTriggerJetAccept(0x0),
192 fhCentralityAccept(0x0),
194 fHLeadingJetPtRaw(0x0),
195 fHDphiVsJetPtAll(0x0),
196 fhJetPtGenVsJetPtRec(0x0),
198 fh2NtriggersGen(0x0),
200 fhPtTrkTruePrimRec(0x0),
201 fhPtTrkTruePrimGen(0x0),
202 fhPtTrkSecOrFakeRec(0x0),
206 fkAcceptance(2.0*TMath::Pi()*1.8),
207 fkDeltaPhiCut(TMath::Pi()-0.6),
213 fh1PtHardTrials(0x0),
216 fEventNumberRangeLow(0),
217 fEventNumberRangeHigh(99),
218 fTriggerPtRangeLow(0.0),
219 fTriggerPtRangeHigh(10000.0),
224 DefineOutput(1, TList::Class());
227 //--------------------------------------------------------------
228 AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP(const AliAnalysisTaskJetCorePP& a):
229 AliAnalysisTaskSE(a.GetName()),
233 fAODExtension(a.fAODExtension),
234 fJetBranchName(a.fJetBranchName),
235 fJetBranchNameMC(a.fJetBranchNameMC),
236 fListJets(a.fListJets),
237 fListJetsGen(a.fListJetsGen),
238 fNonStdFile(a.fNonStdFile),
240 fJetParamR(a.fJetParamR),
241 fOfflineTrgMask(a.fOfflineTrgMask),
242 fMinContribVtx(a.fMinContribVtx),
243 fVtxZMin(a.fVtxZMin),
244 fVtxZMax(a.fVtxZMax),
245 fFilterMask(a.fFilterMask),
246 fCentMin(a.fCentMin),
247 fCentMax(a.fCentMax),
248 fJetEtaMin(a.fJetEtaMin),
249 fJetEtaMax(a.fJetEtaMax),
250 fTriggerEtaCut(a.fTriggerEtaCut),
251 fTrackEtaCut(a.fTrackEtaCut),
252 fTrackLowPtCut(a.fTrackLowPtCut),
253 fOutputList(a.fOutputList),
254 fHistEvtSelection(a.fHistEvtSelection),
255 fh2Ntriggers(a.fh2Ntriggers),
256 fHJetSpec(a.fHJetSpec),
257 fHJetDensity(a.fHJetDensity),
258 fHJetDensityA4(a.fHJetDensityA4),
259 fhJetPhi(a.fhJetPhi),
260 fhTriggerPhi(a.fhTriggerPhi),
261 fhJetEta(a.fhJetEta),
262 fhTriggerEta(a.fhTriggerEta),
263 fhVertexZ(a.fhVertexZ),
264 fhVertexZAccept(a.fhVertexZAccept),
265 fhContribVtx(a.fhContribVtx),
266 fhContribVtxAccept(a.fhContribVtxAccept),
267 fhDphiTriggerJet(a.fhDphiTriggerJet),
268 fhDphiTriggerJetAccept(a.fhDphiTriggerJetAccept),
269 fhCentrality(a.fhCentrality),
270 fhCentralityAccept(a.fhCentralityAccept),
271 fHJetPtRaw(a.fHJetPtRaw),
272 fHLeadingJetPtRaw(a.fHLeadingJetPtRaw),
273 fHDphiVsJetPtAll(a.fHDphiVsJetPtAll),
274 fhJetPtGenVsJetPtRec(a.fhJetPtGenVsJetPtRec),
275 fhJetPtGen(a.fhJetPtGen),
276 fh2NtriggersGen(a.fh2NtriggersGen),
277 fHJetSpecGen(a.fHJetSpecGen),
278 fhPtTrkTruePrimRec(a.fhPtTrkTruePrimRec),
279 fhPtTrkTruePrimGen(a.fhPtTrkTruePrimGen),
280 fhPtTrkSecOrFakeRec(a.fhPtTrkSecOrFakeRec),
282 faGenIndex(a.faGenIndex),
283 faRecIndex(a.faRecIndex),
284 fkAcceptance(a.fkAcceptance),
285 fkDeltaPhiCut(a.fkDeltaPhiCut),
287 fh1Trials(a.fh1Trials),
288 fh1AvgTrials(a.fh1AvgTrials),
289 fh1PtHard(a.fh1PtHard),
290 fh1PtHardNoW(a.fh1PtHardNoW),
291 fh1PtHardTrials(a.fh1PtHardTrials),
292 fAvgTrials(a.fAvgTrials),
293 fHardest(a.fHardest),
294 fEventNumberRangeLow(a.fEventNumberRangeLow),
295 fEventNumberRangeHigh(a.fEventNumberRangeHigh),
296 fTriggerPtRangeLow(a.fTriggerPtRangeLow),
297 fTriggerPtRangeHigh(a.fTriggerPtRangeHigh)
301 //--------------------------------------------------------------
303 AliAnalysisTaskJetCorePP& AliAnalysisTaskJetCorePP::operator = (const AliAnalysisTaskJetCorePP& a){
304 // assignment operator
305 this->~AliAnalysisTaskJetCorePP();
306 new(this) AliAnalysisTaskJetCorePP(a);
309 //--------------------------------------------------------------
311 AliAnalysisTaskJetCorePP::~AliAnalysisTaskJetCorePP()
316 delete fOutputList; // ????
320 //--------------------------------------------------------------
323 Bool_t AliAnalysisTaskJetCorePP::Notify()
325 //Implemented Notify() to read the cross sections
326 //and number of trials from pyxsec.root
327 //inspired by AliAnalysisTaskJetSpectrum2::Notify()
328 if(!fIsMC) return kFALSE;
330 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
332 if(fDebug>1) AliError("ESD not available");
333 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
336 fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
339 if(fNonStdFile.Length()!=0){
340 // case that we have an AOD extension we can fetch the jets from the extended output
341 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
342 fAODExtension = aodH ? aodH->GetExtension(fNonStdFile.Data()) : 0;
344 if(fDebug>1) Printf("AODExtension found for %s",fNonStdFile.Data());
348 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
349 Float_t xsection = 0;
354 TFile *curfile = tree->GetCurrentFile();
356 Error("Notify","No current file");
359 if(!fh1Xsec || !fh1Trials){
360 Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
363 AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
364 fh1Xsec->Fill("<#sigma>",xsection);
365 // construct a poor man average trials
366 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
367 if(ftrials>=nEntries && nEntries>0.) fAvgTrials = ftrials/nEntries;
368 fh1Trials->Fill("#sum{ntrials}",ftrials);
371 if(fDebug)Printf("Reading File %s",fInputHandler->GetTree()->GetCurrentFile()->GetName());
375 //--------------------------------------------------------------
377 void AliAnalysisTaskJetCorePP::Init()
379 // check for jet branches
380 if(!strlen(fJetBranchName.Data())){
381 AliError("Jet branch name not set.");
386 //--------------------------------------------------------------
388 void AliAnalysisTaskJetCorePP::UserCreateOutputObjects()
393 fListJets = new TList(); //reconstructed level
395 fIsMC = (fJetBranchNameMC.Length()>0) ? kTRUE : kFALSE;
397 fRandom = new TRandom3(0);
399 if(fIsMC) fListJetsGen = new TList(); //generator level
401 if(!fOutputList) fOutputList = new TList();
402 fOutputList->SetOwner(kTRUE);
404 Bool_t oldStatus = TH1::AddDirectoryStatus();
405 TH1::AddDirectory(kFALSE);
407 fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
408 fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
409 fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
410 fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
411 fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
412 fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
413 fHistEvtSelection->GetXaxis()->SetBinLabel(6,"event number (rejected)");
415 fOutputList->Add(fHistEvtSelection);
417 Int_t nBinsCentrality = (fSystem==0) ? 1 : 10; // pp=1 else 10
418 //trigger pt spectrum (reconstructed)
419 fh2Ntriggers = new TH2F("fh2Ntriggers","# of triggers",
420 nBinsCentrality,0.0,100.0,50,0.0,50.0);
421 fOutputList->Add(fh2Ntriggers);
423 //Centrality, A, pTjet, pTtrigg
424 const Int_t dimSpec = 4;
425 const Int_t nBinsSpec[dimSpec] = {nBinsCentrality, 100, 100, 50};
426 const Double_t lowBinSpec[dimSpec] = {0.0, 0.0, 0, 0.0};
427 const Double_t hiBinSpec[dimSpec] = {100.0, 1.0, 200.0,50.0};
428 fHJetSpec = new THnSparseF("fHJetSpec",
429 "Recoil jet spectrum [cent,A,pTjet,pTtrig]",
430 dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
431 fOutputList->Add(fHJetSpec);
433 //------------------- HISTOS FOR DIAGNOSTIC ----------------------
434 //Jet number density histos [Trk Mult, jet density, pT trigger]
435 const Int_t dimJetDens = 3;
436 const Int_t nBinsJetDens[dimJetDens] = {100, 100, 10};
437 const Double_t lowBinJetDens[dimJetDens] = {0.0, 0.0, 0.0};
438 const Double_t hiBinJetDens[dimJetDens] = {500.0, 5.0, 50.0};
440 fHJetDensity = new THnSparseF("fHJetDensity","Jet dens vs trk mult A>0.07",
441 dimJetDens,nBinsJetDens,lowBinJetDens,hiBinJetDens);
443 fHJetDensityA4 =new THnSparseF("fHJetDensityA4","Jet dens vs trk mult A>0.4",
444 dimJetDens,nBinsJetDens,lowBinJetDens,hiBinJetDens);
446 fOutputList->Add(fHJetDensity);
447 fOutputList->Add(fHJetDensityA4);
450 //inclusive azimuthal and pseudorapidity histograms
451 fhJetPhi = new TH2D("fhJetPhi","Azim dist jets vs pTjet",
452 50, 0, 100, 50,-TMath::Pi(),TMath::Pi());
453 fhTriggerPhi= new TH2D("fhTriggerPhi","azim dist trig had vs pTtrigg",
454 25, 0, 50, 50,-TMath::Pi(),TMath::Pi());
455 fhJetEta = new TH2D("fhJetEta","Eta dist jets vs pTjet",
456 50,0, 100, 40,-0.9,0.9);
457 fhTriggerEta = new TH2D("fhTriggerEta","Eta dist trig had vs pTtrigg",
458 25, 0, 50, 40,-0.9,0.9);
460 fhVertexZ = new TH1D("fhVertexZ","z vertex",40,-20,20);
461 fhVertexZAccept = new TH1D("fhVertexZAccept","z vertex after cut",40,-20,20);
462 fhContribVtx = new TH1D("fhContribVtx","contrib to vtx",200,0,200);
463 fhContribVtxAccept = new TH1D("fhContribVtxAccept","contrib to vtx after cut",200,0,200);
464 fhDphiTriggerJet = new TH1D("fhDphiTriggerJet","Deltaphi trig-jet",50, -TMath::Pi(),TMath::Pi());
465 fhDphiTriggerJetAccept = new TH1D("fhDphiTriggerJetAccept","Deltaphi trig-jet after cut",50, -TMath::Pi(),TMath::Pi());
466 fhCentrality = new TH1D("fhCentrality","Centrality",20,0,100);
467 fhCentralityAccept = new TH1D("fhCentralityAccept","Centrality after cut",20,0,100);
469 fOutputList->Add(fhJetPhi);
470 fOutputList->Add(fhTriggerPhi);
471 fOutputList->Add(fhJetEta);
472 fOutputList->Add(fhTriggerEta);
473 fOutputList->Add(fhVertexZ);
474 fOutputList->Add(fhVertexZAccept);
475 fOutputList->Add(fhContribVtx);
476 fOutputList->Add(fhContribVtxAccept);
477 fOutputList->Add(fhDphiTriggerJet);
478 fOutputList->Add(fhDphiTriggerJetAccept);
479 fOutputList->Add(fhCentrality);
480 fOutputList->Add(fhCentralityAccept);
482 // raw spectra of INCLUSIVE jets
483 //Centrality, pTjet, A
484 const Int_t dimRaw = 3;
485 const Int_t nBinsRaw[dimRaw] = {nBinsCentrality, 50, 100};
486 const Double_t lowBinRaw[dimRaw] = {0.0, 0.0, 0.0};
487 const Double_t hiBinRaw[dimRaw] = {100.0, 100, 1.0};
488 fHJetPtRaw = new THnSparseF("fHJetPtRaw",
489 "Incl. jet spectrum [cent,pTjet,A]",
490 dimRaw,nBinsRaw,lowBinRaw,hiBinRaw);
491 fOutputList->Add(fHJetPtRaw);
493 // raw spectra of LEADING jets
494 //Centrality, pTjet, A
495 fHLeadingJetPtRaw = new THnSparseF("fHLeadingJetPtRaw",
496 "Leading jet spectrum [cent,pTjet,A]",
497 dimRaw,nBinsRaw,lowBinRaw,hiBinRaw);
498 fOutputList->Add(fHLeadingJetPtRaw);
500 // Dphi versus pT jet
501 //Centrality, Dphi=phiTrig-phiJet, pTjet, pTtrigg
502 const Int_t dimDp = 4;
503 const Int_t nBinsDp[dimDp] = {nBinsCentrality, 50, 50, 50};
504 const Double_t lowBinDp[dimDp] = {0.0, -TMath::Pi(), 0.0, 0.0};
505 const Double_t hiBinDp[dimDp] = {100.0, TMath::Pi(), 100.0, 100.0};
506 fHDphiVsJetPtAll = new THnSparseF("fHDphiVsJetPtAll",
507 "Dphi vs jet pT [cent,Dphi,pTjet,pTtrigg]",
508 dimDp,nBinsDp,lowBinDp,hiBinDp);
509 fOutputList->Add(fHDphiVsJetPtAll);
512 //analyze MC generator level
514 fhJetPtGenVsJetPtRec = new TH2D("fhJetPtGenVsJetPtRec","JetPtGenVsJetPtRec", 100,0,200, 100,0,200);
515 fOutputList->Add(fhJetPtGenVsJetPtRec); //gen MC jet pt spectrum versus reconstructed pt spectrum
517 fhJetPtGen = new TH1D("fhJetPtGen","Jet Pt (MC Gen)",100,0,200); //MC generator jet pt spectrum
518 fOutputList->Add(fhJetPtGen);
520 fh2NtriggersGen = (TH2F*) fh2Ntriggers->Clone("fh2NtriggersGen");
521 fh2NtriggersGen->SetTitle(Form("%s Gen MC",fh2Ntriggers->GetTitle()));
522 fOutputList->Add(fh2NtriggersGen);
524 //Centrality, A, pT, pTtrigg
525 fHJetSpecGen = (THnSparseF*) fHJetSpec->Clone("fHJetSpecGen");
526 fHJetSpecGen->SetTitle(Form("%s Gen MC",fHJetSpec->GetTitle()));
527 fOutputList->Add(fHJetSpecGen);
529 //track efficiency/contamination histograms eta versus pT
530 fhPtTrkTruePrimRec = new TH2D("fhPtTrkTruePrimRec","PtTrkTruePrimRec",100,0,20,18,-0.9,0.9);
531 fOutputList->Add(fhPtTrkTruePrimRec);
533 fhPtTrkTruePrimGen = (TH2D*) fhPtTrkTruePrimRec->Clone("fhPtTrkTruePrimGen");
534 fhPtTrkTruePrimGen->SetTitle("PtTrkTruePrimGen");
535 fOutputList->Add(fhPtTrkTruePrimGen);
537 fhPtTrkSecOrFakeRec = (TH2D*) fhPtTrkTruePrimRec->Clone("fhPtTrkSecOrFakeRec");
538 fhPtTrkSecOrFakeRec->SetTitle("PtTrkSecOrFakeRec");
539 fOutputList->Add(fhPtTrkSecOrFakeRec);
541 //-------------------------------------
543 const Int_t nBinPt = 150;
544 Double_t binLimitsPt[nBinPt+1];
545 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
547 binLimitsPt[iPt] = -50.0;
549 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.;
553 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
554 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
555 fOutputList->Add(fh1Xsec);
556 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
557 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
558 fOutputList->Add(fh1Trials);
559 fh1AvgTrials = new TH1F("fh1AvgTrials","trials root file",1,0,1);
560 fh1AvgTrials->GetXaxis()->SetBinLabel(1,"#sum{avg ntrials}");
561 fOutputList->Add(fh1AvgTrials);
562 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
563 fOutputList->Add(fh1PtHard);
564 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
565 fOutputList->Add(fh1PtHardNoW);
566 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
567 fOutputList->Add(fh1PtHardTrials);
570 // =========== Switch on Sumw2 for all histos ===========
571 for(Int_t i=0; i<fOutputList->GetEntries(); i++){
572 TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
577 THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
582 TH1::AddDirectory(oldStatus);
584 PostData(1, fOutputList);
587 //--------------------------------------------------------------------
589 void AliAnalysisTaskJetCorePP::UserExec(Option_t *)
592 Double_t eventW = 1.0;
593 Double_t ptHard = 0.0;
594 Double_t nTrials = 1.0; // Trials for MC trigger
595 if(fIsMC) fh1AvgTrials->Fill("#sum{avg ntrials}",fAvgTrials);
597 if(TMath::Abs((Float_t) fJetParamR)<0.00001){
598 AliError("Cone radius is set to zero.");
601 if(!strlen(fJetBranchName.Data())){
602 AliError("Jet branch name not set.");
606 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
608 if(fDebug>1) AliError("ESD not available");
609 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
612 fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
613 AliAODEvent* aod = NULL;
614 // take all other information from the aod we take the tracks from
616 if(!fESD) aod = fAODIn;
620 if(fNonStdFile.Length()!=0){
621 // case that we have an AOD extension we can fetch the jets from the extended output
622 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
623 fAODExtension = aodH ? aodH->GetExtension(fNonStdFile.Data()) : 0;
625 if(fDebug>1) Printf("AODExtension found for %s",fNonStdFile.Data());
629 // ----------------- EVENT SELECTION --------------------------
630 fHistEvtSelection->Fill(1); // number of events before event selection
633 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
634 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
636 if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
637 if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
638 fHistEvtSelection->Fill(2);
639 PostData(1, fOutputList);
645 if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
646 fHistEvtSelection->Fill(3);
647 PostData(1, fOutputList);
652 AliAODVertex* primVtx = aod->GetPrimaryVertex();
655 if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
656 fHistEvtSelection->Fill(3);
657 PostData(1, fOutputList);
661 Int_t nTracksPrim = primVtx->GetNContributors();
662 Float_t vtxz = primVtx->GetZ();
664 fhContribVtx->Fill(nTracksPrim);
665 if( nTracksPrim > 0 ) fhVertexZ->Fill(vtxz);
667 if((nTracksPrim < fMinContribVtx) ||
670 if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",
671 (char*)__FILE__,__LINE__,vtxz);
672 fHistEvtSelection->Fill(3);
673 PostData(1, fOutputList);
677 fhContribVtxAccept->Fill(nTracksPrim);
678 fhVertexZAccept->Fill(vtxz);
680 //FK// No event class selection imposed
681 // event class selection (from jet helper task)
682 //Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
683 //if(fDebug) Printf("Event class %d", eventClass);
684 //if(eventClass < fEvtClassMin || eventClass > fEvtClassMax){
685 // fHistEvtSelection->Fill(4);
686 // PostData(1, fOutputList);
690 //------------------ CENTRALITY SELECTION ---------------
691 AliCentrality *cent = 0x0;
692 Double_t centValue = 0.0;
693 if(fSystem){ //fSystem=0 for pp, fSystem=1 for pPb
695 cent = fESD->GetCentrality();
696 if(cent) centValue = cent->GetCentralityPercentile("V0M");
698 centValue = aod->GetHeader()->GetCentrality();
700 if(fDebug) printf("centrality: %f\n", centValue);
702 fhCentrality->Fill(centValue);
704 if(centValue < fCentMin || centValue > fCentMax){
705 fHistEvtSelection->Fill(4);
706 PostData(1, fOutputList);
710 fhCentralityAccept->Fill( centValue );
714 //-----------------select disjunct event subsamples ----------------
715 Int_t eventnum = aod->GetHeader()->GetEventNumberESDFile();
716 Int_t lastdigit = eventnum % 10;
717 if(!(fEventNumberRangeLow<=lastdigit && lastdigit<=fEventNumberRangeHigh)){
718 fHistEvtSelection->Fill(5);
719 PostData(1, fOutputList);
723 if(fDebug) std::cout<<" ACCEPTED EVENT "<<endl;
724 fHistEvtSelection->Fill(0); // accepted events
725 // ------------------- end event selection --------------------
728 // fetch RECONSTRUCTED jets
729 TClonesArray *aodJets = 0x0;
731 if(fAODOut && !aodJets){
732 aodJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName.Data()));
734 if(fAODExtension && !aodJets){
735 aodJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName.Data()));
737 if(fAODIn && !aodJets){
738 aodJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName.Data()));
741 //--------- Fill list of RECONSTRUCTED jets --------------
744 if(fDebug) Printf("########## %s: %d jets",fJetBranchName.Data(), aodJets->GetEntriesFast());
745 for(Int_t iJet = 0; iJet < aodJets->GetEntriesFast(); iJet++) {
746 AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets)[iJet]);
747 if (jet) fListJets->Add(jet);
751 //--------- Fill list of MC GENERATOR LEVEL jets --------------
752 TList particleListGen; //list of tracks in MC
754 if(fIsMC){ //analyze generator level MC
755 // fetch MC generator level jets
756 TClonesArray *aodGenJets = NULL;
758 if(fAODOut&&!aodGenJets){
759 aodGenJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchNameMC.Data()));
761 if(fAODExtension&&!aodGenJets){
762 aodGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchNameMC.Data()));
764 if(fAODIn&&!aodGenJets){
765 aodGenJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchNameMC.Data()));
769 Printf("%s:%d no generated Jet array with name %s in AOD",
770 (char*)__FILE__,__LINE__, fJetBranchNameMC.Data());
771 PostData(1, fOutputList);
775 fListJetsGen->Clear();
777 //serarch for charged trigger at the MC generator level
778 Int_t indexTriggGen = -1;
779 Double_t ptTriggGen = -1;
780 Int_t iCounterGen = 0;
781 Int_t triggersMC[200];//list of trigger candidates
782 Int_t ntriggersMC = 0; //index in triggers array
786 AliMCEvent* mcEvent = MCEvent();
788 PostData(1, fOutputList);
792 AliGenPythiaEventHeader *pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
794 nTrials = pythiaGenHeader->Trials();
795 ptHard = pythiaGenHeader->GetPtHard();
796 fh1PtHard->Fill(ptHard,eventW);
797 fh1PtHardNoW->Fill(ptHard,1);
798 fh1PtHardTrials->Fill(ptHard,nTrials);
801 for(Int_t it = 0; it < mcEvent->GetNumberOfTracks(); it++){
802 if(!mcEvent->IsPhysicalPrimary(it)) continue;
803 AliMCParticle* part = (AliMCParticle*) mcEvent->GetTrack(it);
805 if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){
807 if(fHardest==0 && ntriggersMC<200){//single inclusive trigger
808 if(indexTriggGen > -1){//trigger candidater was found
809 triggersMC[ntriggersMC] = indexTriggGen;
819 PostData(1, fOutputList);
822 TClonesArray *tca = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(AliAODMCParticle::StdBranchName()));
824 PostData(1, fOutputList);
827 for(Int_t it = 0; it < tca->GetEntriesFast(); it++){
828 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
830 if(!part->IsPhysicalPrimary()) continue;
831 if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){
833 if(fHardest==0 && ntriggersMC<200){//single inclusive trigger
834 if(indexTriggGen > -1){ //trigger candidater was found
835 triggersMC[ntriggersMC] = indexTriggGen;
845 if(fHardest==0){ //single inclusive trigger
846 if(ntriggersMC>0){ //there is at least one trigger
847 Int_t rnd = fRandom->Integer(ntriggersMC); //0 to ntriggers-1
848 indexTriggGen = triggersMC[rnd];
850 indexTriggGen = -1; //trigger not found
854 //================== Fill jet list ===================
856 if(fDebug) Printf("########## %s: %d jets",fJetBranchNameMC.Data(), aodGenJets->GetEntriesFast());
858 for(Int_t igJet = 0; igJet < aodGenJets->GetEntriesFast(); igJet++) {
859 AliAODJet *jetGen = dynamic_cast<AliAODJet*>((*aodGenJets)[igJet]);
860 if(jetGen) fListJetsGen->Add(jetGen);
864 //============ Generator trigger+jet ==================
865 Int_t ilowGen = (fHardest==0 || fHardest==1) ? indexTriggGen : 0;
866 Int_t ihighGen = (fHardest==0 || fHardest==1) ? indexTriggGen+1 : particleListGen.GetEntries();
868 for(Int_t igen= ilowGen; igen< ihighGen; igen++){ //loop over possible trigger
869 indexTriggGen = igen; //trigger hadron
871 if(indexTriggGen == -1) continue;
872 AliVParticle* triggerGen=(AliVParticle*) particleListGen.At(indexTriggGen);
873 if(!triggerGen) continue;
876 if(triggerGen->Pt() < 10.0) continue; //all hadrons pt>10
878 if(TMath::Abs((Float_t) triggerGen->Eta()) > fTriggerEtaCut) continue;
880 ptTriggGen = triggerGen->Pt(); //count triggers
881 fh2NtriggersGen->Fill(centValue, ptTriggGen);
883 //Count jets and trigger-jet pairs at MC generator level
884 if(!aodGenJets) continue;
885 for(Int_t ij=0; ij<fListJetsGen->GetEntries(); ij++){
886 AliAODJet* jet = (AliAODJet*)(fListJetsGen->At(ij));
888 Double_t etaJetGen = jet->Eta();
890 if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){
892 Double_t dphi = RelativePhi(triggerGen->Phi(), jet->Phi());
893 if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue;
895 //Centrality, A, pT, pTtrigg
896 Double_t fillspecgen[] = { centValue,
897 jet->EffectiveAreaCharged(),
901 fHJetSpecGen->Fill(fillspecgen);
902 }//back to back jet-trigger pair
907 //================ RESPONSE MATRIX ===============
909 //Count jets and trigger-jet pairs at MC generator level
910 for(Int_t ij=0; ij<fListJetsGen->GetEntries(); ij++){
911 AliAODJet* jet = (AliAODJet*)(fListJetsGen->At(ij));
913 Double_t etaJetGen = jet->Eta();
914 Double_t ptJetGen = jet->Pt();
916 if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){
917 fhJetPtGen->Fill(ptJetGen); // generator level pt spectum of jets response mx normalization
920 if(fListJets->GetEntries()>0 && fListJetsGen->GetEntries()>0){ //at least some reconstructed jets
921 Int_t ng = (Int_t) fListJetsGen->GetEntries();
922 Int_t nr = (Int_t) fListJets->GetEntries();
924 //Find closest MC generator - reconstructed jets
925 if(faGenIndex.GetSize()<nr) faGenIndex.Set(nr); //idx of gen jet assoc to rec jet
926 if(faRecIndex.GetSize()<ng) faRecIndex.Set(ng); //idx of rec jet assoc to gen jet
929 Printf("New Rec List %d gen index Array %d",nr,faGenIndex.GetSize());
930 Printf("New Gen List %d rec index Array %d",ng,faRecIndex.GetSize());
932 //matching of MC genrator level and reconstructed jets
933 AliAnalysisHelperJetTasks::GetClosestJets(fListJetsGen,ng,fListJets,nr,faGenIndex,faRecIndex,fDebug);
935 // Fill response matrix
936 for(Int_t ir = 0; ir < nr; ir++){
937 AliAODJet *recJet = (AliAODJet*) fListJets->At(ir);
938 Double_t etaJetRec = recJet->Eta();
939 Double_t ptJetRec = recJet->Pt();
940 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
942 if((fJetEtaMin <= etaJetRec) && (etaJetRec <= fJetEtaMax)){
943 Int_t ig = faGenIndex[ir]; //associated generator level jet
944 if(ig >= 0 && ig < ng){
945 if(fDebug > 10) Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
946 AliAODJet *genJet = (AliAODJet*) fListJetsGen->At(ig);
947 Double_t ptJetGen = genJet->Pt();
948 Double_t etaJetGen = genJet->Eta();
950 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
951 if((fJetEtaMin <= etaJetGen) && (etaJetGen <= fJetEtaMax)){
952 fhJetPtGenVsJetPtRec->Fill(ptJetRec, ptJetGen);
955 }//rec jet in eta acceptance
956 }//loop over reconstructed jets
958 }//pointer MC generator jets
959 } //analyze generator level MC
961 //============= RECONSTRUCTED INCLUSIVE JETS ===============
963 Double_t etaJet = 0.0;
964 Double_t pTJet = 0.0;
965 Double_t areaJet = 0.0;
966 Double_t phiJet = 0.0;
967 Int_t indexLeadingJet = -1;
968 Double_t pTLeadingJet = -10.0;
969 Double_t areaLeadingJet = -10.0;
971 for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
972 AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
977 if(pTJet==0) continue;
979 if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
980 areaJet = jet->EffectiveAreaCharged();
982 //Jet Diagnostics---------------------------------
983 fhJetPhi->Fill(pTJet, RelativePhi(phiJet,0.0)); //phi -pi,pi
984 fhJetEta->Fill(pTJet, etaJet);
985 //search for leading jet
986 if(pTJet > pTLeadingJet){
987 indexLeadingJet = ij;
988 pTLeadingJet = pTJet;
989 areaLeadingJet = areaJet;
992 // raw spectra of INCLUSIVE jets
993 //Centrality, pTjet, A
994 Double_t fillraw[] = { centValue,
998 fHJetPtRaw->Fill(fillraw);
1001 if(indexLeadingJet > -1){
1002 // raw spectra of LEADING jets
1003 //Centrality, pTjet, A
1004 Double_t fillleading[] = { centValue,
1008 fHLeadingJetPtRaw->Fill(fillleading);
1012 // =============== RECONSTRUCTED TRIGGER-JET PAIRS ================
1013 //Find Hadron trigger
1014 TList particleList; //list of tracks
1015 Int_t indexTrigg = GetListOfTracks(&particleList); //index of trigger hadron in Particle list
1016 if(fIsMC) FillEffHistos(&particleList, &particleListGen); //Fill efficiency histos
1018 //set ranges of the trigger loop
1019 Int_t ilow = (fHardest==0 || fHardest==1) ? indexTrigg : 0;
1020 Int_t ihigh = (fHardest==0 || fHardest==1) ? indexTrigg+1 : particleList.GetEntries();
1022 for(Int_t itrk= ilow; itrk< ihigh; itrk++){ //loop over possible trigger
1023 indexTrigg = itrk; //trigger hadron with pT >10 GeV
1025 if(indexTrigg < 0) continue;
1027 AliVParticle *triggerHadron = (AliVParticle*) particleList.At(indexTrigg);
1029 PostData(1, fOutputList);
1033 if(triggerHadron->Pt() < 10.0) continue; //all hadrons pt>10
1035 if(TMath::Abs((Float_t) triggerHadron->Eta())> fTriggerEtaCut) continue;
1037 //Fill trigger histograms
1038 fh2Ntriggers->Fill(centValue,triggerHadron->Pt()); //trigger pT
1039 fhTriggerPhi->Fill(triggerHadron->Pt(),RelativePhi(triggerHadron->Phi(),0.0)); //phi -pi,pi
1040 fhTriggerEta->Fill(triggerHadron->Pt(),triggerHadron->Eta());
1043 //---------- make trigger-jet pairs ---------
1047 for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
1048 AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
1050 etaJet = jet->Eta();
1051 phiJet = jet->Phi();
1053 if(pTJet==0) continue;
1055 if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
1056 areaJet = jet->EffectiveAreaCharged();
1057 if(areaJet >= 0.07) injet++;
1058 if(areaJet >= 0.4) injet4++;
1060 Double_t dphi = RelativePhi(triggerHadron->Phi(), phiJet);
1061 fhDphiTriggerJet->Fill(dphi); //Input
1063 //Dphi versus jet pT
1064 //Centrality, Dphi=phiTrig-phiJet, pTjet, pTtrigg
1065 Double_t filldp[] = { centValue,
1070 fHDphiVsJetPtAll->Fill(filldp);
1072 // Select back to back trigger - jet pairs
1073 if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue;
1074 fhDphiTriggerJetAccept->Fill(dphi); //Accepted
1077 //Centrality, A, pTjet, pTtrigg
1078 Double_t fillspec[] = { centValue,
1083 fHJetSpec->Fill(fillspec);
1086 //Fill Jet Density In the Event A>0.07
1088 Double_t filldens[]={ (Double_t) particleList.GetEntries(),
1092 fHJetDensity->Fill(filldens);
1095 //Fill Jet Density In the Event A>0.4
1097 Double_t filldens4[]={ (Double_t) particleList.GetEntries(),
1098 injet4/fkAcceptance,
1101 fHJetDensityA4->Fill(filldens4);
1106 PostData(1, fOutputList);
1109 //----------------------------------------------------------------------------
1110 void AliAnalysisTaskJetCorePP::Terminate(const Option_t *)
1112 // Draw result to the screen
1113 // Called once at the end of the query
1115 if(fDebug) printf("AliAnalysisTaskJetCorePP DONE\n");
1116 if(!GetOutputData(1)) return;
1120 //----------------------------------------------------------------------------
1121 Int_t AliAnalysisTaskJetCorePP::GetListOfTracks(TList *list){
1122 //Fill the list of accepted tracks (passed track cut)
1123 //return consecutive index of the hardest ch hadron in the list
1125 AliAODEvent *aodevt = NULL;
1127 if(!fESD) aodevt = fAODIn;
1128 else aodevt = fAODOut;
1130 if(!aodevt) return -1;
1132 Int_t index = -1; //index of the highest particle in the list
1133 Double_t ptmax = -10;
1134 Int_t triggers[200];
1135 Int_t ntriggers = 0; //index in triggers array
1137 for(Int_t it = 0; it < aodevt->GetNumberOfTracks(); it++){
1138 AliAODTrack *tr = aodevt->GetTrack(it);
1140 //if((fFilterMask > 0) && !(tr->TestFilterBit(fFilterMask))) continue;
1141 if((fFilterMask > 0) && !(tr->IsHybridGlobalConstrainedGlobal())) continue;
1142 if(TMath::Abs((Float_t) tr->Eta()) > fTrackEtaCut) continue;
1143 if(tr->Pt() < fTrackLowPtCut) continue;
1147 if(fHardest>0){ //leading particle
1154 if(fHardest==0 && ntriggers<200){ //single inclusive
1155 if(fTriggerPtRangeLow <= tr->Pt() &&
1156 tr->Pt() < fTriggerPtRangeHigh){
1157 triggers[ntriggers] = iCount;
1165 if(fHardest==0 && ntriggers>0){ //select random inclusive trigger
1166 Int_t rnd = fRandom->Integer(ntriggers); //0 to ntriggers-1
1167 index = triggers[rnd];
1173 //----------------------------------------------------------------------------
1175 Double_t AliAnalysisTaskJetCorePP::GetBackgroundInPerpCone(Float_t jetR, Double_t jetPhi, Double_t jetEta, TList* trkList){
1176 //calculate sum of track pT in the cone perpendicular in phi to the jet
1177 //jetR = cone radius
1178 // jetPhi, jetEta = direction of the jet
1179 Int_t numberOfTrks = trkList->GetEntries();
1180 Double_t pTsum = 0.0;
1181 Double_t perpConePhi = jetPhi + TMath::Pi()/2;//perp cone w.r.t. jet in phi
1182 for(Int_t it=0; it<numberOfTrks; it++){
1183 AliVParticle *trk = (AliVParticle*) trkList->At(it);
1184 Double_t dphi = RelativePhi(perpConePhi,trk->Phi());
1185 Double_t deta = trk->Eta()-jetEta;
1187 if( (dphi*dphi + deta*deta)< (jetR*jetR)){
1195 //----------------------------------------------------------------------------
1197 Double_t AliAnalysisTaskJetCorePP::RelativePhi(Double_t mphi,Double_t vphi){
1198 //Get relative azimuthal angle of two particles -pi to pi
1199 if (vphi < -TMath::Pi()) vphi += TMath::TwoPi();
1200 else if (vphi > TMath::Pi()) vphi -= TMath::TwoPi();
1202 if (mphi < -TMath::Pi()) mphi += TMath::TwoPi();
1203 else if (mphi > TMath::Pi()) mphi -= TMath::TwoPi();
1205 Double_t dphi = mphi - vphi;
1206 if (dphi < -TMath::Pi()) dphi += TMath::TwoPi();
1207 else if (dphi > TMath::Pi()) dphi -= TMath::TwoPi();
1209 return dphi;//dphi in [-Pi, Pi]
1213 //----------------------------------------------------------------------------
1214 Bool_t AliAnalysisTaskJetCorePP::SelectMCGenTracks(AliVParticle *trk, TList *trkList, Double_t &ptLeading, Int_t &index, Int_t counter){
1215 //fill track efficiency denominator
1216 if(trk->Pt() < fTrackLowPtCut) return kFALSE;
1217 if(trk->Charge()==0) return kFALSE;
1218 if(TMath::Abs((Float_t) trk->Eta()) > fTrackEtaCut) return kFALSE;
1220 fhPtTrkTruePrimGen->Fill(trk->Pt(),trk->Eta()); //Efficiency denominator
1223 if(fHardest>0){ //leading particle
1224 if(ptLeading < trk->Pt()){
1226 ptLeading = trk->Pt();
1230 if(fHardest==0){ //single inclusive
1232 if(fTriggerPtRangeLow <= trk->Pt() &&
1233 trk->Pt() < fTriggerPtRangeHigh){
1241 //----------------------------------------------------------------------------
1242 void AliAnalysisTaskJetCorePP::FillEffHistos(TList *recList, TList *genList){
1243 //fill track efficiency numerator
1244 if(!recList) return;
1245 if(!genList) return;
1246 Int_t nRec = recList->GetEntries();
1247 Int_t nGen = genList->GetEntries();
1248 for(Int_t ir=0; ir<nRec; ir++){
1249 AliAODTrack *trkRec = (AliAODTrack*) recList->At(ir);
1250 if(!trkRec) continue;
1251 Int_t recLabel = TMath::Abs(trkRec->GetLabel());
1252 Bool_t matched = kFALSE;
1254 for(Int_t ig=0; ig<nGen; ig++){
1255 AliVParticle *trkGen = (AliVParticle*) genList->At(ig);
1256 if(!trkGen) continue;
1257 Int_t genLabel = TMath::Abs(trkGen->GetLabel());
1258 if(recLabel==genLabel){
1259 fhPtTrkTruePrimRec->Fill(trkGen->Pt(),trkGen->Eta());
1264 if(!matched) fhPtTrkSecOrFakeRec->Fill(trkRec->Pt(),trkRec->Eta());