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),
303 //--------------------------------------------------------------
305 AliAnalysisTaskJetCorePP& AliAnalysisTaskJetCorePP::operator = (const AliAnalysisTaskJetCorePP& a){
306 // assignment operator
307 this->~AliAnalysisTaskJetCorePP();
308 new(this) AliAnalysisTaskJetCorePP(a);
311 //--------------------------------------------------------------
313 AliAnalysisTaskJetCorePP::~AliAnalysisTaskJetCorePP()
318 delete fOutputList; // ????
322 //--------------------------------------------------------------
325 Bool_t AliAnalysisTaskJetCorePP::Notify()
327 //Implemented Notify() to read the cross sections
328 //and number of trials from pyxsec.root
329 //inspired by AliAnalysisTaskJetSpectrum2::Notify()
330 if(!fIsMC) return kFALSE;
332 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
334 if(fDebug>1) AliError("ESD not available");
335 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
338 fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
341 if(fNonStdFile.Length()!=0){
342 // case that we have an AOD extension we can fetch the jets from the extended output
343 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
344 fAODExtension = aodH ? aodH->GetExtension(fNonStdFile.Data()) : 0;
346 if(fDebug>1) Printf("AODExtension found for %s",fNonStdFile.Data());
350 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
351 Float_t xsection = 0;
356 TFile *curfile = tree->GetCurrentFile();
358 Error("Notify","No current file");
361 if(!fh1Xsec || !fh1Trials){
362 Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
365 AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
366 fh1Xsec->Fill("<#sigma>",xsection);
367 // construct a poor man average trials
368 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
369 if(ftrials>=nEntries && nEntries>0.) fAvgTrials = ftrials/nEntries;
370 fh1Trials->Fill("#sum{ntrials}",ftrials);
373 if(fDebug)Printf("Reading File %s",fInputHandler->GetTree()->GetCurrentFile()->GetName());
377 //--------------------------------------------------------------
379 void AliAnalysisTaskJetCorePP::Init()
381 // check for jet branches
382 if(!strlen(fJetBranchName.Data())){
383 AliError("Jet branch name not set.");
388 //--------------------------------------------------------------
390 void AliAnalysisTaskJetCorePP::UserCreateOutputObjects()
395 fListJets = new TList(); //reconstructed level
397 fIsMC = (fJetBranchNameMC.Length()>0) ? kTRUE : kFALSE;
399 fRandom = new TRandom3(0);
401 if(fIsMC) fListJetsGen = new TList(); //generator level
403 if(!fOutputList) fOutputList = new TList();
404 fOutputList->SetOwner(kTRUE);
406 Bool_t oldStatus = TH1::AddDirectoryStatus();
407 TH1::AddDirectory(kFALSE);
409 fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
410 fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
411 fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
412 fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
413 fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
414 fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
415 fHistEvtSelection->GetXaxis()->SetBinLabel(6,"event number (rejected)");
417 fOutputList->Add(fHistEvtSelection);
419 Int_t nBinsCentrality = (fSystem==0) ? 1 : 10; // pp=1 else 10
420 //trigger pt spectrum (reconstructed)
421 fh2Ntriggers = new TH2F("fh2Ntriggers","# of triggers",
422 nBinsCentrality,0.0,100.0,50,0.0,50.0);
423 fOutputList->Add(fh2Ntriggers);
425 //Centrality, A, pTjet, pTtrigg
426 const Int_t dimSpec = 4;
427 const Int_t nBinsSpec[dimSpec] = {nBinsCentrality, 100, 100, 50};
428 const Double_t lowBinSpec[dimSpec] = {0.0, 0.0, 0, 0.0};
429 const Double_t hiBinSpec[dimSpec] = {100.0, 1.0, 200.0,50.0};
430 fHJetSpec = new THnSparseF("fHJetSpec",
431 "Recoil jet spectrum [cent,A,pTjet,pTtrig]",
432 dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
433 fOutputList->Add(fHJetSpec);
435 //------------------- HISTOS FOR DIAGNOSTIC ----------------------
436 //Jet number density histos [Trk Mult, jet density, pT trigger]
437 const Int_t dimJetDens = 3;
438 const Int_t nBinsJetDens[dimJetDens] = {100, 100, 10};
439 const Double_t lowBinJetDens[dimJetDens] = {0.0, 0.0, 0.0};
440 const Double_t hiBinJetDens[dimJetDens] = {500.0, 5.0, 50.0};
442 fHJetDensity = new THnSparseF("fHJetDensity","Jet dens vs trk mult A>0.07",
443 dimJetDens,nBinsJetDens,lowBinJetDens,hiBinJetDens);
445 fHJetDensityA4 =new THnSparseF("fHJetDensityA4","Jet dens vs trk mult A>0.4",
446 dimJetDens,nBinsJetDens,lowBinJetDens,hiBinJetDens);
448 fOutputList->Add(fHJetDensity);
449 fOutputList->Add(fHJetDensityA4);
452 //inclusive azimuthal and pseudorapidity histograms
453 fhJetPhi = new TH2D("fhJetPhi","Azim dist jets vs pTjet",
454 50, 0, 100, 50,-TMath::Pi(),TMath::Pi());
455 fhTriggerPhi= new TH2D("fhTriggerPhi","azim dist trig had vs pTtrigg",
456 25, 0, 50, 50,-TMath::Pi(),TMath::Pi());
457 fhJetEta = new TH2D("fhJetEta","Eta dist jets vs pTjet",
458 50,0, 100, 40,-0.9,0.9);
459 fhTriggerEta = new TH2D("fhTriggerEta","Eta dist trig had vs pTtrigg",
460 25, 0, 50, 40,-0.9,0.9);
462 fhVertexZ = new TH1D("fhVertexZ","z vertex",40,-20,20);
463 fhVertexZAccept = new TH1D("fhVertexZAccept","z vertex after cut",40,-20,20);
464 fhContribVtx = new TH1D("fhContribVtx","contrib to vtx",200,0,200);
465 fhContribVtxAccept = new TH1D("fhContribVtxAccept","contrib to vtx after cut",200,0,200);
466 fhDphiTriggerJet = new TH1D("fhDphiTriggerJet","Deltaphi trig-jet",50, -TMath::Pi(),TMath::Pi());
467 fhDphiTriggerJetAccept = new TH1D("fhDphiTriggerJetAccept","Deltaphi trig-jet after cut",50, -TMath::Pi(),TMath::Pi());
468 fhCentrality = new TH1D("fhCentrality","Centrality",20,0,100);
469 fhCentralityAccept = new TH1D("fhCentralityAccept","Centrality after cut",20,0,100);
471 fOutputList->Add(fhJetPhi);
472 fOutputList->Add(fhTriggerPhi);
473 fOutputList->Add(fhJetEta);
474 fOutputList->Add(fhTriggerEta);
475 fOutputList->Add(fhVertexZ);
476 fOutputList->Add(fhVertexZAccept);
477 fOutputList->Add(fhContribVtx);
478 fOutputList->Add(fhContribVtxAccept);
479 fOutputList->Add(fhDphiTriggerJet);
480 fOutputList->Add(fhDphiTriggerJetAccept);
481 fOutputList->Add(fhCentrality);
482 fOutputList->Add(fhCentralityAccept);
484 // raw spectra of INCLUSIVE jets
485 //Centrality, pTjet, A
486 const Int_t dimRaw = 3;
487 const Int_t nBinsRaw[dimRaw] = {nBinsCentrality, 50, 100};
488 const Double_t lowBinRaw[dimRaw] = {0.0, 0.0, 0.0};
489 const Double_t hiBinRaw[dimRaw] = {100.0, 100, 1.0};
490 fHJetPtRaw = new THnSparseF("fHJetPtRaw",
491 "Incl. jet spectrum [cent,pTjet,A]",
492 dimRaw,nBinsRaw,lowBinRaw,hiBinRaw);
493 fOutputList->Add(fHJetPtRaw);
495 // raw spectra of LEADING jets
496 //Centrality, pTjet, A
497 fHLeadingJetPtRaw = new THnSparseF("fHLeadingJetPtRaw",
498 "Leading jet spectrum [cent,pTjet,A]",
499 dimRaw,nBinsRaw,lowBinRaw,hiBinRaw);
500 fOutputList->Add(fHLeadingJetPtRaw);
502 // Dphi versus pT jet
503 //Centrality, Dphi=phiTrig-phiJet, pTjet, pTtrigg
504 const Int_t dimDp = 4;
505 const Int_t nBinsDp[dimDp] = {nBinsCentrality, 50, 50, 50};
506 const Double_t lowBinDp[dimDp] = {0.0, -TMath::Pi(), 0.0, 0.0};
507 const Double_t hiBinDp[dimDp] = {100.0, TMath::Pi(), 100.0, 100.0};
508 fHDphiVsJetPtAll = new THnSparseF("fHDphiVsJetPtAll",
509 "Dphi vs jet pT [cent,Dphi,pTjet,pTtrigg]",
510 dimDp,nBinsDp,lowBinDp,hiBinDp);
511 fOutputList->Add(fHDphiVsJetPtAll);
514 //analyze MC generator level
516 fhJetPtGenVsJetPtRec = new TH2D("fhJetPtGenVsJetPtRec","JetPtGenVsJetPtRec", 100,0,200, 100,0,200);
517 fOutputList->Add(fhJetPtGenVsJetPtRec); //gen MC jet pt spectrum versus reconstructed pt spectrum
519 fhJetPtGen = new TH1D("fhJetPtGen","Jet Pt (MC Gen)",100,0,200); //MC generator jet pt spectrum
520 fOutputList->Add(fhJetPtGen);
522 fh2NtriggersGen = (TH2F*) fh2Ntriggers->Clone("fh2NtriggersGen");
523 fh2NtriggersGen->SetTitle(Form("%s Gen MC",fh2Ntriggers->GetTitle()));
524 fOutputList->Add(fh2NtriggersGen);
526 //Centrality, A, pT, pTtrigg
527 fHJetSpecGen = (THnSparseF*) fHJetSpec->Clone("fHJetSpecGen");
528 fHJetSpecGen->SetTitle(Form("%s Gen MC",fHJetSpec->GetTitle()));
529 fOutputList->Add(fHJetSpecGen);
531 //track efficiency/contamination histograms eta versus pT
532 fhPtTrkTruePrimRec = new TH2D("fhPtTrkTruePrimRec","PtTrkTruePrimRec",100,0,20,18,-0.9,0.9);
533 fOutputList->Add(fhPtTrkTruePrimRec);
535 fhPtTrkTruePrimGen = (TH2D*) fhPtTrkTruePrimRec->Clone("fhPtTrkTruePrimGen");
536 fhPtTrkTruePrimGen->SetTitle("PtTrkTruePrimGen");
537 fOutputList->Add(fhPtTrkTruePrimGen);
539 fhPtTrkSecOrFakeRec = (TH2D*) fhPtTrkTruePrimRec->Clone("fhPtTrkSecOrFakeRec");
540 fhPtTrkSecOrFakeRec->SetTitle("PtTrkSecOrFakeRec");
541 fOutputList->Add(fhPtTrkSecOrFakeRec);
543 //-------------------------------------
545 const Int_t nBinPt = 150;
546 Double_t binLimitsPt[nBinPt+1];
547 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
549 binLimitsPt[iPt] = -50.0;
551 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.;
555 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
556 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
557 fOutputList->Add(fh1Xsec);
558 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
559 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
560 fOutputList->Add(fh1Trials);
561 fh1AvgTrials = new TH1F("fh1AvgTrials","trials root file",1,0,1);
562 fh1AvgTrials->GetXaxis()->SetBinLabel(1,"#sum{avg ntrials}");
563 fOutputList->Add(fh1AvgTrials);
564 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
565 fOutputList->Add(fh1PtHard);
566 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
567 fOutputList->Add(fh1PtHardNoW);
568 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
569 fOutputList->Add(fh1PtHardTrials);
572 // =========== Switch on Sumw2 for all histos ===========
573 for(Int_t i=0; i<fOutputList->GetEntries(); i++){
574 TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
579 THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
584 TH1::AddDirectory(oldStatus);
586 PostData(1, fOutputList);
589 //--------------------------------------------------------------------
591 void AliAnalysisTaskJetCorePP::UserExec(Option_t *)
594 Double_t eventW = 1.0;
595 Double_t ptHard = 0.0;
596 Double_t nTrials = 1.0; // Trials for MC trigger
597 if(fIsMC) fh1AvgTrials->Fill("#sum{avg ntrials}",fAvgTrials);
599 if(TMath::Abs((Float_t) fJetParamR)<0.00001){
600 AliError("Cone radius is set to zero.");
603 if(!strlen(fJetBranchName.Data())){
604 AliError("Jet branch name not set.");
608 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
610 if(fDebug>1) AliError("ESD not available");
611 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
614 fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
615 AliAODEvent* aod = NULL;
616 // take all other information from the aod we take the tracks from
618 if(!fESD) aod = fAODIn;
622 if(fNonStdFile.Length()!=0){
623 // case that we have an AOD extension we can fetch the jets from the extended output
624 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
625 fAODExtension = aodH ? aodH->GetExtension(fNonStdFile.Data()) : 0;
627 if(fDebug>1) Printf("AODExtension found for %s",fNonStdFile.Data());
631 // ----------------- EVENT SELECTION --------------------------
632 fHistEvtSelection->Fill(1); // number of events before event selection
635 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
636 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
638 if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
639 if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
640 fHistEvtSelection->Fill(2);
641 PostData(1, fOutputList);
647 if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
648 fHistEvtSelection->Fill(3);
649 PostData(1, fOutputList);
654 AliAODVertex* primVtx = aod->GetPrimaryVertex();
657 if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
658 fHistEvtSelection->Fill(3);
659 PostData(1, fOutputList);
663 Int_t nTracksPrim = primVtx->GetNContributors();
664 Float_t vtxz = primVtx->GetZ();
666 fhContribVtx->Fill(nTracksPrim);
667 if( nTracksPrim > 0 ) fhVertexZ->Fill(vtxz);
669 if((nTracksPrim < fMinContribVtx) ||
672 if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",
673 (char*)__FILE__,__LINE__,vtxz);
674 fHistEvtSelection->Fill(3);
675 PostData(1, fOutputList);
679 fhContribVtxAccept->Fill(nTracksPrim);
680 fhVertexZAccept->Fill(vtxz);
682 //FK// No event class selection imposed
683 // event class selection (from jet helper task)
684 //Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
685 //if(fDebug) Printf("Event class %d", eventClass);
686 //if(eventClass < fEvtClassMin || eventClass > fEvtClassMax){
687 // fHistEvtSelection->Fill(4);
688 // PostData(1, fOutputList);
692 //------------------ CENTRALITY SELECTION ---------------
693 AliCentrality *cent = 0x0;
694 Double_t centValue = 0.0;
695 if(fSystem){ //fSystem=0 for pp, fSystem=1 for pPb
697 cent = fESD->GetCentrality();
698 if(cent) centValue = cent->GetCentralityPercentile("V0M");
700 centValue = aod->GetHeader()->GetCentrality();
702 if(fDebug) printf("centrality: %f\n", centValue);
704 fhCentrality->Fill(centValue);
706 if(centValue < fCentMin || centValue > fCentMax){
707 fHistEvtSelection->Fill(4);
708 PostData(1, fOutputList);
712 fhCentralityAccept->Fill( centValue );
716 //-----------------select disjunct event subsamples ----------------
717 Int_t eventnum = aod->GetHeader()->GetEventNumberESDFile();
718 Int_t lastdigit = eventnum % 10;
719 if(!(fEventNumberRangeLow<=lastdigit && lastdigit<=fEventNumberRangeHigh)){
720 fHistEvtSelection->Fill(5);
721 PostData(1, fOutputList);
725 if(fDebug) std::cout<<" ACCEPTED EVENT "<<endl;
726 fHistEvtSelection->Fill(0); // accepted events
727 // ------------------- end event selection --------------------
730 // fetch RECONSTRUCTED jets
731 TClonesArray *aodJets = 0x0;
733 if(fAODOut && !aodJets){
734 aodJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName.Data()));
736 if(fAODExtension && !aodJets){
737 aodJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName.Data()));
739 if(fAODIn && !aodJets){
740 aodJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName.Data()));
743 //--------- Fill list of RECONSTRUCTED jets --------------
746 if(fDebug) Printf("########## %s: %d jets",fJetBranchName.Data(), aodJets->GetEntriesFast());
747 for(Int_t iJet = 0; iJet < aodJets->GetEntriesFast(); iJet++) {
748 AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets)[iJet]);
749 if (jet) fListJets->Add(jet);
753 //--------- Fill list of MC GENERATOR LEVEL jets --------------
754 TList particleListGen; //list of tracks in MC
756 if(fIsMC){ //analyze generator level MC
757 // fetch MC generator level jets
758 TClonesArray *aodGenJets = NULL;
760 if(fAODOut&&!aodGenJets){
761 aodGenJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchNameMC.Data()));
763 if(fAODExtension&&!aodGenJets){
764 aodGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchNameMC.Data()));
766 if(fAODIn&&!aodGenJets){
767 aodGenJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchNameMC.Data()));
771 Printf("%s:%d no generated Jet array with name %s in AOD",
772 (char*)__FILE__,__LINE__, fJetBranchNameMC.Data());
773 PostData(1, fOutputList);
777 fListJetsGen->Clear();
779 //serarch for charged trigger at the MC generator level
780 Int_t indexTriggGen = -1;
781 Double_t ptTriggGen = -1;
782 Int_t iCounterGen = 0;
783 Int_t triggersMC[200];//list of trigger candidates
784 Int_t ntriggersMC = 0; //index in triggers array
788 AliMCEvent* mcEvent = MCEvent();
790 PostData(1, fOutputList);
794 AliGenPythiaEventHeader *pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
796 nTrials = pythiaGenHeader->Trials();
797 ptHard = pythiaGenHeader->GetPtHard();
798 fh1PtHard->Fill(ptHard,eventW);
799 fh1PtHardNoW->Fill(ptHard,1);
800 fh1PtHardTrials->Fill(ptHard,nTrials);
803 for(Int_t it = 0; it < mcEvent->GetNumberOfTracks(); it++){
804 if(!mcEvent->IsPhysicalPrimary(it)) continue;
805 AliMCParticle* part = (AliMCParticle*) mcEvent->GetTrack(it);
807 if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){
809 if(fHardest==0 && ntriggersMC<200){//single inclusive trigger
810 if(indexTriggGen > -1){//trigger candidater was found
811 triggersMC[ntriggersMC] = indexTriggGen;
821 PostData(1, fOutputList);
824 TClonesArray *tca = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(AliAODMCParticle::StdBranchName()));
826 PostData(1, fOutputList);
829 for(Int_t it = 0; it < tca->GetEntriesFast(); it++){
830 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
832 if(!part->IsPhysicalPrimary()) continue;
833 if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){
835 if(fHardest==0 && ntriggersMC<200){//single inclusive trigger
836 if(indexTriggGen > -1){ //trigger candidater was found
837 triggersMC[ntriggersMC] = indexTriggGen;
847 if(fHardest==0){ //single inclusive trigger
848 if(ntriggersMC>0){ //there is at least one trigger
849 Int_t rnd = fRandom->Integer(ntriggersMC); //0 to ntriggers-1
850 indexTriggGen = triggersMC[rnd];
852 indexTriggGen = -1; //trigger not found
856 //================== Fill jet list ===================
858 if(fDebug) Printf("########## %s: %d jets",fJetBranchNameMC.Data(), aodGenJets->GetEntriesFast());
860 for(Int_t igJet = 0; igJet < aodGenJets->GetEntriesFast(); igJet++) {
861 AliAODJet *jetGen = dynamic_cast<AliAODJet*>((*aodGenJets)[igJet]);
862 if(jetGen) fListJetsGen->Add(jetGen);
866 //============ Generator trigger+jet ==================
867 Int_t ilowGen = (fHardest==0 || fHardest==1) ? indexTriggGen : 0;
868 Int_t ihighGen = (fHardest==0 || fHardest==1) ? indexTriggGen+1 : particleListGen.GetEntries();
870 for(Int_t igen= ilowGen; igen< ihighGen; igen++){ //loop over possible trigger
871 indexTriggGen = igen; //trigger hadron
873 if(indexTriggGen == -1) continue;
874 AliVParticle* triggerGen=(AliVParticle*) particleListGen.At(indexTriggGen);
875 if(!triggerGen) continue;
878 if(triggerGen->Pt() < 6.0) continue; //all hadrons pt>10
880 if(TMath::Abs((Float_t) triggerGen->Eta()) > fTriggerEtaCut) continue;
882 ptTriggGen = triggerGen->Pt(); //count triggers
883 fh2NtriggersGen->Fill(centValue, ptTriggGen);
885 //Count jets and trigger-jet pairs at MC generator level
886 if(!aodGenJets) continue;
887 for(Int_t ij=0; ij<fListJetsGen->GetEntries(); ij++){
888 AliAODJet* jet = (AliAODJet*)(fListJetsGen->At(ij));
890 Double_t etaJetGen = jet->Eta();
892 if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){
894 Double_t dphi = RelativePhi(triggerGen->Phi(), jet->Phi());
895 if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue;
897 //Centrality, A, pT, pTtrigg
898 Double_t fillspecgen[] = { centValue,
899 jet->EffectiveAreaCharged(),
903 fHJetSpecGen->Fill(fillspecgen);
904 }//back to back jet-trigger pair
909 //================ RESPONSE MATRIX ===============
911 //Count jets and trigger-jet pairs at MC generator level
912 for(Int_t ij=0; ij<fListJetsGen->GetEntries(); ij++){
913 AliAODJet* jet = (AliAODJet*)(fListJetsGen->At(ij));
915 Double_t etaJetGen = jet->Eta();
916 Double_t ptJetGen = jet->Pt();
918 if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){
919 fhJetPtGen->Fill(ptJetGen); // generator level pt spectum of jets response mx normalization
922 if(fListJets->GetEntries()>0 && fListJetsGen->GetEntries()>0){ //at least some reconstructed jets
923 Int_t ng = (Int_t) fListJetsGen->GetEntries();
924 Int_t nr = (Int_t) fListJets->GetEntries();
926 //Find closest MC generator - reconstructed jets
927 if(faGenIndex.GetSize()<nr) faGenIndex.Set(nr); //idx of gen jet assoc to rec jet
928 if(faRecIndex.GetSize()<ng) faRecIndex.Set(ng); //idx of rec jet assoc to gen jet
931 Printf("New Rec List %d gen index Array %d",nr,faGenIndex.GetSize());
932 Printf("New Gen List %d rec index Array %d",ng,faRecIndex.GetSize());
934 //matching of MC genrator level and reconstructed jets
935 AliAnalysisHelperJetTasks::GetClosestJets(fListJetsGen,ng,fListJets,nr,faGenIndex,faRecIndex,fDebug);
937 // Fill response matrix
938 for(Int_t ir = 0; ir < nr; ir++){
939 AliAODJet *recJet = (AliAODJet*) fListJets->At(ir);
940 Double_t etaJetRec = recJet->Eta();
941 Double_t ptJetRec = recJet->Pt();
942 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
944 if((fJetEtaMin <= etaJetRec) && (etaJetRec <= fJetEtaMax)){
945 Int_t ig = faGenIndex[ir]; //associated generator level jet
946 if(ig >= 0 && ig < ng){
947 if(fDebug > 10) Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
948 AliAODJet *genJet = (AliAODJet*) fListJetsGen->At(ig);
949 Double_t ptJetGen = genJet->Pt();
950 Double_t etaJetGen = genJet->Eta();
952 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
953 if((fJetEtaMin <= etaJetGen) && (etaJetGen <= fJetEtaMax)){
954 fhJetPtGenVsJetPtRec->Fill(ptJetRec, ptJetGen);
957 }//rec jet in eta acceptance
958 }//loop over reconstructed jets
960 }//pointer MC generator jets
961 } //analyze generator level MC
963 //============= RECONSTRUCTED INCLUSIVE JETS ===============
965 Double_t etaJet = 0.0;
966 Double_t pTJet = 0.0;
967 Double_t areaJet = 0.0;
968 Double_t phiJet = 0.0;
969 Int_t indexLeadingJet = -1;
970 Double_t pTLeadingJet = -10.0;
971 Double_t areaLeadingJet = -10.0;
973 for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
974 AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
979 if(pTJet==0) continue;
981 if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
982 areaJet = jet->EffectiveAreaCharged();
984 //Jet Diagnostics---------------------------------
985 fhJetPhi->Fill(pTJet, RelativePhi(phiJet,0.0)); //phi -pi,pi
986 fhJetEta->Fill(pTJet, etaJet);
987 //search for leading jet
988 if(pTJet > pTLeadingJet){
989 indexLeadingJet = ij;
990 pTLeadingJet = pTJet;
991 areaLeadingJet = areaJet;
994 // raw spectra of INCLUSIVE jets
995 //Centrality, pTjet, A
996 Double_t fillraw[] = { centValue,
1000 fHJetPtRaw->Fill(fillraw);
1003 if(indexLeadingJet > -1){
1004 // raw spectra of LEADING jets
1005 //Centrality, pTjet, A
1006 Double_t fillleading[] = { centValue,
1010 fHLeadingJetPtRaw->Fill(fillleading);
1014 // =============== RECONSTRUCTED TRIGGER-JET PAIRS ================
1015 //Find Hadron trigger
1016 TList particleList; //list of tracks
1017 Int_t indexTrigg = GetListOfTracks(&particleList); //index of trigger hadron in Particle list
1018 if(fIsMC) FillEffHistos(&particleList, &particleListGen); //Fill efficiency histos
1020 //set ranges of the trigger loop
1021 Int_t ilow = (fHardest==0 || fHardest==1) ? indexTrigg : 0;
1022 Int_t ihigh = (fHardest==0 || fHardest==1) ? indexTrigg+1 : particleList.GetEntries();
1024 for(Int_t itrk= ilow; itrk< ihigh; itrk++){ //loop over possible trigger
1025 indexTrigg = itrk; //trigger hadron with pT >10 GeV
1027 if(indexTrigg < 0) continue;
1029 AliVParticle *triggerHadron = (AliVParticle*) particleList.At(indexTrigg);
1031 PostData(1, fOutputList);
1035 if(triggerHadron->Pt() < 6.0) continue; //all hadrons pt>6
1037 if(TMath::Abs((Float_t) triggerHadron->Eta())> fTriggerEtaCut) continue;
1039 //Fill trigger histograms
1040 fh2Ntriggers->Fill(centValue,triggerHadron->Pt()); //trigger pT
1041 fhTriggerPhi->Fill(triggerHadron->Pt(),RelativePhi(triggerHadron->Phi(),0.0)); //phi -pi,pi
1042 fhTriggerEta->Fill(triggerHadron->Pt(),triggerHadron->Eta());
1045 //---------- make trigger-jet pairs ---------
1049 for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
1050 AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
1052 etaJet = jet->Eta();
1053 phiJet = jet->Phi();
1055 if(pTJet==0) continue;
1057 if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
1058 areaJet = jet->EffectiveAreaCharged();
1059 if(areaJet >= 0.07) injet++;
1060 if(areaJet >= 0.4) injet4++;
1062 Double_t dphi = RelativePhi(triggerHadron->Phi(), phiJet);
1063 fhDphiTriggerJet->Fill(dphi); //Input
1065 //Dphi versus jet pT
1066 //Centrality, Dphi=phiTrig-phiJet, pTjet, pTtrigg
1067 Double_t filldp[] = { centValue,
1072 fHDphiVsJetPtAll->Fill(filldp);
1074 // Select back to back trigger - jet pairs
1075 if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue;
1076 fhDphiTriggerJetAccept->Fill(dphi); //Accepted
1079 //Centrality, A, pTjet, pTtrigg
1080 Double_t fillspec[] = { centValue,
1085 fHJetSpec->Fill(fillspec);
1088 //Fill Jet Density In the Event A>0.07
1090 Double_t filldens[]={ (Double_t) particleList.GetEntries(),
1094 fHJetDensity->Fill(filldens);
1097 //Fill Jet Density In the Event A>0.4
1099 Double_t filldens4[]={ (Double_t) particleList.GetEntries(),
1100 injet4/fkAcceptance,
1103 fHJetDensityA4->Fill(filldens4);
1108 PostData(1, fOutputList);
1111 //----------------------------------------------------------------------------
1112 void AliAnalysisTaskJetCorePP::Terminate(const Option_t *)
1114 // Draw result to the screen
1115 // Called once at the end of the query
1117 if(fDebug) printf("AliAnalysisTaskJetCorePP DONE\n");
1118 if(!GetOutputData(1)) return;
1122 //----------------------------------------------------------------------------
1123 Int_t AliAnalysisTaskJetCorePP::GetListOfTracks(TList *list){
1124 //Fill the list of accepted tracks (passed track cut)
1125 //return consecutive index of the hardest ch hadron in the list
1127 AliAODEvent *aodevt = NULL;
1129 if(!fESD) aodevt = fAODIn;
1130 else aodevt = fAODOut;
1132 if(!aodevt) return -1;
1134 Int_t index = -1; //index of the highest particle in the list
1135 Double_t ptmax = -10;
1136 Int_t triggers[200];
1137 Int_t ntriggers = 0; //index in triggers array
1139 for(Int_t it = 0; it < aodevt->GetNumberOfTracks(); it++){
1140 AliAODTrack *tr = aodevt->GetTrack(it);
1142 //if((fFilterMask > 0) && !(tr->TestFilterBit(fFilterMask))) continue;
1143 if((fFilterMask > 0) && !(tr->IsHybridGlobalConstrainedGlobal())) continue;
1144 if(TMath::Abs((Float_t) tr->Eta()) > fTrackEtaCut) continue;
1145 if(tr->Pt() < fTrackLowPtCut) continue;
1149 if(fHardest>0){ //leading particle
1156 if(fHardest==0 && ntriggers<200){ //single inclusive
1157 if(fTriggerPtRangeLow <= tr->Pt() &&
1158 tr->Pt() < fTriggerPtRangeHigh){
1159 triggers[ntriggers] = iCount;
1167 if(fHardest==0 && ntriggers>0){ //select random inclusive trigger
1168 Int_t rnd = fRandom->Integer(ntriggers); //0 to ntriggers-1
1169 index = triggers[rnd];
1175 //----------------------------------------------------------------------------
1177 Double_t AliAnalysisTaskJetCorePP::GetBackgroundInPerpCone(Float_t jetR, Double_t jetPhi, Double_t jetEta, TList* trkList){
1178 //calculate sum of track pT in the cone perpendicular in phi to the jet
1179 //jetR = cone radius
1180 // jetPhi, jetEta = direction of the jet
1181 Int_t numberOfTrks = trkList->GetEntries();
1182 Double_t pTsum = 0.0;
1183 Double_t perpConePhi = jetPhi + TMath::Pi()/2;//perp cone w.r.t. jet in phi
1184 for(Int_t it=0; it<numberOfTrks; it++){
1185 AliVParticle *trk = (AliVParticle*) trkList->At(it);
1186 Double_t dphi = RelativePhi(perpConePhi,trk->Phi());
1187 Double_t deta = trk->Eta()-jetEta;
1189 if( (dphi*dphi + deta*deta)< (jetR*jetR)){
1197 //----------------------------------------------------------------------------
1199 Double_t AliAnalysisTaskJetCorePP::RelativePhi(Double_t mphi,Double_t vphi){
1200 //Get relative azimuthal angle of two particles -pi to pi
1201 if (vphi < -TMath::Pi()) vphi += TMath::TwoPi();
1202 else if (vphi > TMath::Pi()) vphi -= TMath::TwoPi();
1204 if (mphi < -TMath::Pi()) mphi += TMath::TwoPi();
1205 else if (mphi > TMath::Pi()) mphi -= TMath::TwoPi();
1207 Double_t dphi = mphi - vphi;
1208 if (dphi < -TMath::Pi()) dphi += TMath::TwoPi();
1209 else if (dphi > TMath::Pi()) dphi -= TMath::TwoPi();
1211 return dphi;//dphi in [-Pi, Pi]
1215 //----------------------------------------------------------------------------
1216 Bool_t AliAnalysisTaskJetCorePP::SelectMCGenTracks(AliVParticle *trk, TList *trkList, Double_t &ptLeading, Int_t &index, Int_t counter){
1217 //fill track efficiency denominator
1218 if(trk->Pt() < fTrackLowPtCut) return kFALSE;
1219 if(trk->Charge()==0) return kFALSE;
1220 if(TMath::Abs((Float_t) trk->Eta()) > fTrackEtaCut) return kFALSE;
1222 fhPtTrkTruePrimGen->Fill(trk->Pt(),trk->Eta()); //Efficiency denominator
1225 if(fHardest>0){ //leading particle
1226 if(ptLeading < trk->Pt()){
1228 ptLeading = trk->Pt();
1232 if(fHardest==0){ //single inclusive
1234 if(fTriggerPtRangeLow <= trk->Pt() &&
1235 trk->Pt() < fTriggerPtRangeHigh){
1243 //----------------------------------------------------------------------------
1244 void AliAnalysisTaskJetCorePP::FillEffHistos(TList *recList, TList *genList){
1245 //fill track efficiency numerator
1246 if(!recList) return;
1247 if(!genList) return;
1248 Int_t nRec = recList->GetEntries();
1249 Int_t nGen = genList->GetEntries();
1250 for(Int_t ir=0; ir<nRec; ir++){
1251 AliAODTrack *trkRec = (AliAODTrack*) recList->At(ir);
1252 if(!trkRec) continue;
1253 Int_t recLabel = TMath::Abs(trkRec->GetLabel());
1254 Bool_t matched = kFALSE;
1256 for(Int_t ig=0; ig<nGen; ig++){
1257 AliVParticle *trkGen = (AliVParticle*) genList->At(ig);
1258 if(!trkGen) continue;
1259 Int_t genLabel = TMath::Abs(trkGen->GetLabel());
1260 if(recLabel==genLabel){
1261 fhPtTrkTruePrimRec->Fill(trkGen->Pt(),trkGen->Eta());
1266 if(!matched) fhPtTrkSecOrFakeRec->Fill(trkRec->Pt(),trkRec->Eta());