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"
44 #include "AliAnalysisTask.h"
45 #include "AliAnalysisManager.h"
47 #include "AliVEvent.h"
48 #include "AliESDEvent.h"
49 #include "AliMCEvent.h"
50 #include "AliESDInputHandler.h"
51 #include "AliCentrality.h"
52 #include "AliAnalysisHelperJetTasks.h"
53 #include "AliInputEventHandler.h"
54 #include "AliAODJetEventBackground.h"
55 #include "AliGenPythiaEventHeader.h"
56 #include "AliAODMCParticle.h"
57 #include "AliMCParticle.h"
58 #include "AliAODEvent.h"
59 #include "AliAODHandler.h"
60 #include "AliAODJet.h"
62 #include "AliAnalysisTaskJetCorePP.h"
67 ClassImp(AliAnalysisTaskJetCorePP)
69 //Filip Krizek 1st March 2013
71 //---------------------------------------------------------------------
72 AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP() :
83 fSystem(0), //pp=0 pPb=1
85 fOfflineTrgMask(AliVEvent::kAny),
98 fHistEvtSelection(0x0),
108 fhVertexZAccept(0x0),
110 fhContribVtxAccept(0x0),
111 fhDphiTriggerJet(0x0),
112 fhDphiTriggerJetAccept(0x0),
114 fhCentralityAccept(0x0),
116 fHLeadingJetPtRaw(0x0),
117 fHDphiVsJetPtAll(0x0),
118 fhJetPtGenVsJetPtRec(0x0),
120 fh2NtriggersGen(0x0),
122 fhPtTrkTruePrimRec(0x0),
123 fhPtTrkTruePrimGen(0x0),
124 fhPtTrkSecOrFakeRec(0x0),
128 fkAcceptance(2.0*TMath::Pi()*1.8),
129 fkDeltaPhiCut(TMath::Pi()-0.6),
135 fh1PtHardTrials(0x0),
138 // default Constructor
141 //---------------------------------------------------------------------
143 AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP(const char *name) :
144 AliAnalysisTaskSE(name),
150 fJetBranchNameMC(""),
154 fSystem(0), //pp=0 pPb=1
156 fOfflineTrgMask(AliVEvent::kAny),
167 fTrackLowPtCut(0.15),
169 fHistEvtSelection(0x0),
179 fhVertexZAccept(0x0),
181 fhContribVtxAccept(0x0),
182 fhDphiTriggerJet(0x0),
183 fhDphiTriggerJetAccept(0x0),
185 fhCentralityAccept(0x0),
187 fHLeadingJetPtRaw(0x0),
188 fHDphiVsJetPtAll(0x0),
189 fhJetPtGenVsJetPtRec(0x0),
191 fh2NtriggersGen(0x0),
193 fhPtTrkTruePrimRec(0x0),
194 fhPtTrkTruePrimGen(0x0),
195 fhPtTrkSecOrFakeRec(0x0),
199 fkAcceptance(2.0*TMath::Pi()*1.8),
200 fkDeltaPhiCut(TMath::Pi()-0.6),
206 fh1PtHardTrials(0x0),
211 DefineOutput(1, TList::Class());
214 //--------------------------------------------------------------
215 AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP(const AliAnalysisTaskJetCorePP& a):
216 AliAnalysisTaskSE(a.GetName()),
220 fAODExtension(a.fAODExtension),
221 fJetBranchName(a.fJetBranchName),
222 fJetBranchNameMC(a.fJetBranchNameMC),
223 fListJets(a.fListJets),
224 fListJetsGen(a.fListJetsGen),
225 fNonStdFile(a.fNonStdFile),
227 fJetParamR(a.fJetParamR),
228 fOfflineTrgMask(a.fOfflineTrgMask),
229 fMinContribVtx(a.fMinContribVtx),
230 fVtxZMin(a.fVtxZMin),
231 fVtxZMax(a.fVtxZMax),
232 fFilterMask(a.fFilterMask),
233 fCentMin(a.fCentMin),
234 fCentMax(a.fCentMax),
235 fJetEtaMin(a.fJetEtaMin),
236 fJetEtaMax(a.fJetEtaMax),
237 fTriggerEtaCut(a.fTriggerEtaCut),
238 fTrackEtaCut(a.fTrackEtaCut),
239 fTrackLowPtCut(a.fTrackLowPtCut),
240 fOutputList(a.fOutputList),
241 fHistEvtSelection(a.fHistEvtSelection),
242 fh2Ntriggers(a.fh2Ntriggers),
243 fHJetSpec(a.fHJetSpec),
244 fHJetDensity(a.fHJetDensity),
245 fHJetDensityA4(a.fHJetDensityA4),
246 fhJetPhi(a.fhJetPhi),
247 fhTriggerPhi(a.fhTriggerPhi),
248 fhJetEta(a.fhJetEta),
249 fhTriggerEta(a.fhTriggerEta),
250 fhVertexZ(a.fhVertexZ),
251 fhVertexZAccept(a.fhVertexZAccept),
252 fhContribVtx(a.fhContribVtx),
253 fhContribVtxAccept(a.fhContribVtxAccept),
254 fhDphiTriggerJet(a.fhDphiTriggerJet),
255 fhDphiTriggerJetAccept(a.fhDphiTriggerJetAccept),
256 fhCentrality(a.fhCentrality),
257 fhCentralityAccept(a.fhCentralityAccept),
258 fHJetPtRaw(a.fHJetPtRaw),
259 fHLeadingJetPtRaw(a.fHLeadingJetPtRaw),
260 fHDphiVsJetPtAll(a.fHDphiVsJetPtAll),
261 fhJetPtGenVsJetPtRec(a.fhJetPtGenVsJetPtRec),
262 fhJetPtGen(a.fhJetPtGen),
263 fh2NtriggersGen(a.fh2NtriggersGen),
264 fHJetSpecGen(a.fHJetSpecGen),
265 fhPtTrkTruePrimRec(a.fhPtTrkTruePrimRec),
266 fhPtTrkTruePrimGen(a.fhPtTrkTruePrimGen),
267 fhPtTrkSecOrFakeRec(a.fhPtTrkSecOrFakeRec),
269 faGenIndex(a.faGenIndex),
270 faRecIndex(a.faRecIndex),
271 fkAcceptance(a.fkAcceptance),
272 fkDeltaPhiCut(a.fkDeltaPhiCut),
274 fh1Trials(a.fh1Trials),
275 fh1AvgTrials(a.fh1AvgTrials),
276 fh1PtHard(a.fh1PtHard),
277 fh1PtHardNoW(a.fh1PtHardNoW),
278 fh1PtHardTrials(a.fh1PtHardTrials),
279 fAvgTrials(a.fAvgTrials)
283 //--------------------------------------------------------------
285 AliAnalysisTaskJetCorePP& AliAnalysisTaskJetCorePP::operator = (const AliAnalysisTaskJetCorePP& a){
286 // assignment operator
287 this->~AliAnalysisTaskJetCorePP();
288 new(this) AliAnalysisTaskJetCorePP(a);
291 //--------------------------------------------------------------
293 AliAnalysisTaskJetCorePP::~AliAnalysisTaskJetCorePP()
298 delete fOutputList; // ????
301 //--------------------------------------------------------------
304 Bool_t AliAnalysisTaskJetCorePP::Notify()
306 //Implemented Notify() to read the cross sections
307 //and number of trials from pyxsec.root
308 //inspired by AliAnalysisTaskJetSpectrum2::Notify()
309 if(!fIsMC) return kFALSE;
311 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
313 if(fDebug>1) AliError("ESD not available");
314 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
317 fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
320 if(fNonStdFile.Length()!=0){
321 // case that we have an AOD extension we can fetch the jets from the extended output
322 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
323 fAODExtension = aodH ? aodH->GetExtension(fNonStdFile.Data()) : 0;
325 if(fDebug>1) Printf("AODExtension found for %s",fNonStdFile.Data());
329 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
330 Float_t xsection = 0;
335 TFile *curfile = tree->GetCurrentFile();
337 Error("Notify","No current file");
340 if(!fh1Xsec || !fh1Trials){
341 Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
344 AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
345 fh1Xsec->Fill("<#sigma>",xsection);
346 // construct a poor man average trials
347 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
348 if(ftrials>=nEntries && nEntries>0.) fAvgTrials = ftrials/nEntries;
349 fh1Trials->Fill("#sum{ntrials}",ftrials);
352 if(fDebug)Printf("Reading File %s",fInputHandler->GetTree()->GetCurrentFile()->GetName());
356 //--------------------------------------------------------------
358 void AliAnalysisTaskJetCorePP::Init()
360 // check for jet branches
361 if(!strlen(fJetBranchName.Data())){
362 AliError("Jet branch name not set.");
367 //--------------------------------------------------------------
369 void AliAnalysisTaskJetCorePP::UserCreateOutputObjects()
374 fListJets = new TList(); //reconstructed level
376 fIsMC = (fJetBranchNameMC.Length()>0) ? kTRUE : kFALSE;
378 if(fIsMC) fListJetsGen = new TList(); //generator level
380 if(!fOutputList) fOutputList = new TList();
381 fOutputList->SetOwner(kTRUE);
383 Bool_t oldStatus = TH1::AddDirectoryStatus();
384 TH1::AddDirectory(kFALSE);
386 fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
387 fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
388 fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
389 fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
390 fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
391 fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
392 fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
394 fOutputList->Add(fHistEvtSelection);
396 Int_t nBinsCentrality = (fSystem==0) ? 1 : 10; // pp=1 else 10
397 //trigger pt spectrum (reconstructed)
398 fh2Ntriggers = new TH2F("fh2Ntriggers","# of triggers",
399 nBinsCentrality,0.0,100.0,50,0.0,50.0);
400 fOutputList->Add(fh2Ntriggers);
402 //Centrality, A, pTjet, pTtrigg
403 const Int_t dimSpec = 4;
404 const Int_t nBinsSpec[dimSpec] = {nBinsCentrality, 100, 100, 50};
405 const Double_t lowBinSpec[dimSpec] = {0.0, 0.0, 0, 0.0};
406 const Double_t hiBinSpec[dimSpec] = {100.0, 1.0, 200.0,50.0};
407 fHJetSpec = new THnSparseF("fHJetSpec",
408 "Recoil jet spectrum [cent,A,pTjet,pTtrig]",
409 dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
410 fOutputList->Add(fHJetSpec);
412 //------------------- HISTOS FOR DIAGNOSTIC ----------------------
413 //Jet number density histos [Trk Mult, jet density, pT trigger]
414 const Int_t dimJetDens = 3;
415 const Int_t nBinsJetDens[dimJetDens] = {100, 100, 10};
416 const Double_t lowBinJetDens[dimJetDens] = {0.0, 0.0, 0.0};
417 const Double_t hiBinJetDens[dimJetDens] = {500.0, 5.0, 50.0};
419 fHJetDensity = new THnSparseF("fHJetDensity","Jet dens vs trk mult A>0.07",
420 dimJetDens,nBinsJetDens,lowBinJetDens,hiBinJetDens);
422 fHJetDensityA4 =new THnSparseF("fHJetDensityA4","Jet dens vs trk mult A>0.4",
423 dimJetDens,nBinsJetDens,lowBinJetDens,hiBinJetDens);
425 fOutputList->Add(fHJetDensity);
426 fOutputList->Add(fHJetDensityA4);
429 //inclusive azimuthal and pseudorapidity histograms
430 fhJetPhi = new TH2D("fhJetPhi","Azim dist jets vs pTjet",
431 50, 0, 100, 50,-TMath::Pi(),TMath::Pi());
432 fhTriggerPhi= new TH2D("fhTriggerPhi","azim dist trig had vs pTtrigg",
433 25, 0, 50, 50,-TMath::Pi(),TMath::Pi());
434 fhJetEta = new TH2D("fhJetEta","Eta dist jets vs pTjet",
435 50,0, 100, 40,-0.9,0.9);
436 fhTriggerEta = new TH2D("fhTriggerEta","Eta dist trig had vs pTtrigg",
437 25, 0, 50, 40,-0.9,0.9);
439 fhVertexZ = new TH1D("fhVertexZ","z vertex",40,-20,20);
440 fhVertexZAccept = new TH1D("fhVertexZAccept","z vertex after cut",40,-20,20);
441 fhContribVtx = new TH1D("fhContribVtx","contrib to vtx",200,0,200);
442 fhContribVtxAccept = new TH1D("fhContribVtxAccept","contrib to vtx after cut",200,0,200);
443 fhDphiTriggerJet = new TH1D("fhDphiTriggerJet","Deltaphi trig-jet",50, -TMath::Pi(),TMath::Pi());
444 fhDphiTriggerJetAccept = new TH1D("fhDphiTriggerJetAccept","Deltaphi trig-jet after cut",50, -TMath::Pi(),TMath::Pi());
445 fhCentrality = new TH1D("fhCentrality","Centrality",20,0,100);
446 fhCentralityAccept = new TH1D("fhCentralityAccept","Centrality after cut",20,0,100);
448 fOutputList->Add(fhJetPhi);
449 fOutputList->Add(fhTriggerPhi);
450 fOutputList->Add(fhJetEta);
451 fOutputList->Add(fhTriggerEta);
452 fOutputList->Add(fhVertexZ);
453 fOutputList->Add(fhVertexZAccept);
454 fOutputList->Add(fhContribVtx);
455 fOutputList->Add(fhContribVtxAccept);
456 fOutputList->Add(fhDphiTriggerJet);
457 fOutputList->Add(fhDphiTriggerJetAccept);
458 fOutputList->Add(fhCentrality);
459 fOutputList->Add(fhCentralityAccept);
461 // raw spectra of INCLUSIVE jets
462 //Centrality, pTjet, A
463 const Int_t dimRaw = 3;
464 const Int_t nBinsRaw[dimRaw] = {nBinsCentrality, 50, 100};
465 const Double_t lowBinRaw[dimRaw] = {0.0, 0.0, 0.0};
466 const Double_t hiBinRaw[dimRaw] = {100.0, 100, 1.0};
467 fHJetPtRaw = new THnSparseF("fHJetPtRaw",
468 "Incl. jet spectrum [cent,pTjet,A]",
469 dimRaw,nBinsRaw,lowBinRaw,hiBinRaw);
470 fOutputList->Add(fHJetPtRaw);
472 // raw spectra of LEADING jets
473 //Centrality, pTjet, A
474 fHLeadingJetPtRaw = new THnSparseF("fHLeadingJetPtRaw",
475 "Leading jet spectrum [cent,pTjet,A]",
476 dimRaw,nBinsRaw,lowBinRaw,hiBinRaw);
477 fOutputList->Add(fHLeadingJetPtRaw);
479 // Dphi versus pT jet
480 //Centrality, Dphi=phiTrig-phiJet, pTjet, pTtrigg
481 const Int_t dimDp = 4;
482 const Int_t nBinsDp[dimDp] = {nBinsCentrality, 50, 50, 50};
483 const Double_t lowBinDp[dimDp] = {0.0, -TMath::Pi(), 0.0, 0.0};
484 const Double_t hiBinDp[dimDp] = {100.0, TMath::Pi(), 100.0, 100.0};
485 fHDphiVsJetPtAll = new THnSparseF("fHDphiVsJetPtAll",
486 "Dphi vs jet pT [cent,Dphi,pTjet,pTtrigg]",
487 dimDp,nBinsDp,lowBinDp,hiBinDp);
488 fOutputList->Add(fHDphiVsJetPtAll);
491 //analyze MC generator level
493 fhJetPtGenVsJetPtRec = new TH2D("fhJetPtGenVsJetPtRec","JetPtGenVsJetPtRec", 100,0,200, 100,0,200);
494 fOutputList->Add(fhJetPtGenVsJetPtRec); //gen MC jet pt spectrum versus reconstructed pt spectrum
496 fhJetPtGen = new TH1D("fhJetPtGen","Jet Pt (MC Gen)",100,0,200); //MC generator jet pt spectrum
497 fOutputList->Add(fhJetPtGen);
499 fh2NtriggersGen = (TH2F*) fh2Ntriggers->Clone("fh2NtriggersGen");
500 fh2NtriggersGen->SetTitle(Form("%s Gen MC",fh2Ntriggers->GetTitle()));
501 fOutputList->Add(fh2NtriggersGen);
503 //Centrality, A, pT, pTtrigg
504 fHJetSpecGen = (THnSparseF*) fHJetSpec->Clone("fHJetSpecGen");
505 fHJetSpecGen->SetTitle(Form("%s Gen MC",fHJetSpec->GetTitle()));
506 fOutputList->Add(fHJetSpecGen);
508 //track efficiency/contamination histograms eta versus pT
509 fhPtTrkTruePrimRec = new TH2D("fhPtTrkTruePrimRec","PtTrkTruePrimRec",100,0,20,18,-0.9,0.9);
510 fOutputList->Add(fhPtTrkTruePrimRec);
512 fhPtTrkTruePrimGen = (TH2D*) fhPtTrkTruePrimRec->Clone("fhPtTrkTruePrimGen");
513 fhPtTrkTruePrimGen->SetTitle("PtTrkTruePrimGen");
514 fOutputList->Add(fhPtTrkTruePrimGen);
516 fhPtTrkSecOrFakeRec = (TH2D*) fhPtTrkTruePrimRec->Clone("fhPtTrkSecOrFakeRec");
517 fhPtTrkSecOrFakeRec->SetTitle("PtTrkSecOrFakeRec");
518 fOutputList->Add(fhPtTrkSecOrFakeRec);
520 //-------------------------------------
522 const Int_t nBinPt = 150;
523 Double_t binLimitsPt[nBinPt+1];
524 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
526 binLimitsPt[iPt] = -50.0;
528 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.;
532 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
533 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
534 fOutputList->Add(fh1Xsec);
535 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
536 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
537 fOutputList->Add(fh1Trials);
538 fh1AvgTrials = new TH1F("fh1AvgTrials","trials root file",1,0,1);
539 fh1AvgTrials->GetXaxis()->SetBinLabel(1,"#sum{avg ntrials}");
540 fOutputList->Add(fh1AvgTrials);
541 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
542 fOutputList->Add(fh1PtHard);
543 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
544 fOutputList->Add(fh1PtHardNoW);
545 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
546 fOutputList->Add(fh1PtHardTrials);
549 // =========== Switch on Sumw2 for all histos ===========
550 for(Int_t i=0; i<fOutputList->GetEntries(); i++){
551 TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
556 THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
561 TH1::AddDirectory(oldStatus);
563 PostData(1, fOutputList);
566 //--------------------------------------------------------------------
568 void AliAnalysisTaskJetCorePP::UserExec(Option_t *)
571 Double_t eventW = 1.0;
572 Double_t ptHard = 0.0;
573 Double_t nTrials = 1.0; // Trials for MC trigger
574 if(fIsMC) fh1AvgTrials->Fill("#sum{avg ntrials}",fAvgTrials);
576 if(TMath::Abs((Float_t) fJetParamR)<0.00001){
577 AliError("Cone radius is set to zero.");
580 if(!strlen(fJetBranchName.Data())){
581 AliError("Jet branch name not set.");
585 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
587 if(fDebug>1) AliError("ESD not available");
588 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
591 fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
592 AliAODEvent* aod = NULL;
593 // take all other information from the aod we take the tracks from
595 if(!fESD) aod = fAODIn;
599 if(fNonStdFile.Length()!=0){
600 // case that we have an AOD extension we can fetch the jets from the extended output
601 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
602 fAODExtension = aodH ? aodH->GetExtension(fNonStdFile.Data()) : 0;
604 if(fDebug>1) Printf("AODExtension found for %s",fNonStdFile.Data());
608 // ----------------- EVENT SELECTION --------------------------
609 fHistEvtSelection->Fill(1); // number of events before event selection
612 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
613 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
615 if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
616 if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
617 fHistEvtSelection->Fill(2);
618 PostData(1, fOutputList);
624 if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
625 fHistEvtSelection->Fill(3);
626 PostData(1, fOutputList);
631 AliAODVertex* primVtx = aod->GetPrimaryVertex();
634 if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
635 fHistEvtSelection->Fill(3);
636 PostData(1, fOutputList);
640 Int_t nTracksPrim = primVtx->GetNContributors();
641 Float_t vtxz = primVtx->GetZ();
643 fhContribVtx->Fill(nTracksPrim);
644 if( nTracksPrim > 0 ) fhVertexZ->Fill(vtxz);
646 if((nTracksPrim < fMinContribVtx) ||
649 if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",
650 (char*)__FILE__,__LINE__,vtxz);
651 fHistEvtSelection->Fill(3);
652 PostData(1, fOutputList);
656 fhContribVtxAccept->Fill(nTracksPrim);
657 fhVertexZAccept->Fill(vtxz);
659 //FK// No event class selection imposed
660 // event class selection (from jet helper task)
661 //Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
662 //if(fDebug) Printf("Event class %d", eventClass);
663 //if(eventClass < fEvtClassMin || eventClass > fEvtClassMax){
664 // fHistEvtSelection->Fill(4);
665 // PostData(1, fOutputList);
669 //------------------ CENTRALITY SELECTION ---------------
670 AliCentrality *cent = 0x0;
671 Double_t centValue = 0.0;
672 if(fSystem){ //fSystem=0 for pp, fSystem=1 for pPb
674 cent = fESD->GetCentrality();
675 if(cent) centValue = cent->GetCentralityPercentile("V0M");
677 centValue = aod->GetHeader()->GetCentrality();
679 if(fDebug) printf("centrality: %f\n", centValue);
681 fhCentrality->Fill(centValue);
683 if(centValue < fCentMin || centValue > fCentMax){
684 fHistEvtSelection->Fill(4);
685 PostData(1, fOutputList);
689 fhCentralityAccept->Fill( centValue );
693 if(fDebug) std::cout<<" ACCEPTED EVENT "<<endl;
695 fHistEvtSelection->Fill(0); // accepted events
701 // ------------------- end event selection --------------------
704 // fetch RECONSTRUCTED jets
705 TClonesArray *aodJets = 0x0;
707 if(fAODOut && !aodJets){
708 aodJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName.Data()));
710 if(fAODExtension && !aodJets){
711 aodJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName.Data()));
713 if(fAODIn && !aodJets){
714 aodJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName.Data()));
717 //--------- Fill list of RECONSTRUCTED jets --------------
720 if(fDebug) Printf("########## %s: %d jets",fJetBranchName.Data(), aodJets->GetEntriesFast());
721 for(Int_t iJet = 0; iJet < aodJets->GetEntriesFast(); iJet++) {
722 AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets)[iJet]);
723 if (jet) fListJets->Add(jet);
727 //--------- Fill list of MC GENERATOR LEVEL jets --------------
728 TList particleListGen; //list of tracks in MC
730 if(fIsMC){ //analyze generator level MC
731 // fetch MC generator level jets
732 TClonesArray *aodGenJets = NULL;
734 if(fAODOut&&!aodGenJets){
735 aodGenJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchNameMC.Data()));
737 if(fAODExtension&&!aodGenJets){
738 aodGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchNameMC.Data()));
740 if(fAODIn&&!aodGenJets){
741 aodGenJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchNameMC.Data()));
745 Printf("%s:%d no generated Jet array with name %s in AOD",
746 (char*)__FILE__,__LINE__, fJetBranchNameMC.Data());
747 PostData(1, fOutputList);
751 fListJetsGen->Clear();
753 //serarch for charged trigger at the MC generator level
754 Int_t indexTriggGen = -1;
755 Double_t ptTriggGen = -1;
756 Int_t iCounterGen = 0;
759 AliMCEvent* mcEvent = MCEvent();
761 PostData(1, fOutputList);
765 AliGenPythiaEventHeader *pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
767 nTrials = pythiaGenHeader->Trials();
768 ptHard = pythiaGenHeader->GetPtHard();
769 fh1PtHard->Fill(ptHard,eventW);
770 fh1PtHardNoW->Fill(ptHard,1);
771 fh1PtHardTrials->Fill(ptHard,nTrials);
774 for(Int_t it = 0; it < mcEvent->GetNumberOfTracks(); it++){
775 if(!mcEvent->IsPhysicalPrimary(it)) continue;
776 AliMCParticle* part = (AliMCParticle*) mcEvent->GetTrack(it);
778 if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){
784 PostData(1, fOutputList);
787 TClonesArray *tca = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(AliAODMCParticle::StdBranchName()));
789 PostData(1, fOutputList);
792 for(Int_t it = 0; it < tca->GetEntriesFast(); it++){
793 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
795 if(!part->IsPhysicalPrimary()) continue;
796 if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){
803 if(indexTriggGen > -1){ //Fill trigger histogram
804 AliVParticle* triggerGen = (AliVParticle*) particleListGen.At(indexTriggGen);
805 if(triggerGen && TMath::Abs((Float_t) triggerGen->Eta()) < fTriggerEtaCut){
806 fh2NtriggersGen->Fill(centValue, ptTriggGen);
814 if(fDebug) Printf("########## %s: %d jets",fJetBranchNameMC.Data(), aodGenJets->GetEntriesFast());
816 for(Int_t igJet = 0; igJet < aodGenJets->GetEntriesFast(); igJet++) {
817 AliAODJet *jetGen = dynamic_cast<AliAODJet*>((*aodGenJets)[igJet]);
818 if(jetGen) fListJetsGen->Add(jetGen);
821 //Count jets and trigger-jet pairs at MC generator level
822 for(Int_t ij=0; ij<fListJetsGen->GetEntries(); ij++){
823 AliAODJet* jet = (AliAODJet*)(fListJetsGen->At(ij));
825 Double_t etaJetGen = jet->Eta();
826 Double_t ptJetGen = jet->Pt();
828 if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){
829 fhJetPtGen->Fill(ptJetGen); // generator level pt spectum of jets response mx normalization
831 if(indexTriggGen > -1){
832 AliVParticle* triggerGen = (AliVParticle*) particleListGen.At(indexTriggGen);
834 Double_t dphi = RelativePhi(triggerGen->Phi(), jet->Phi());
835 if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue;
837 //Centrality, A, pT, pTtrigg
838 Double_t fillspecgen[] = { centValue,
839 jet->EffectiveAreaCharged(),
843 fHJetSpecGen->Fill(fillspecgen);
848 if(fListJets->GetEntries()>0 && fListJetsGen->GetEntries()>0){ //at least some reconstructed jets
849 Int_t ng = (Int_t) fListJetsGen->GetEntries();
850 Int_t nr = (Int_t) fListJets->GetEntries();
852 //Find closest MC generator - reconstructed jets
853 if(faGenIndex.GetSize()<nr) faGenIndex.Set(nr); //idx of gen jet assoc to rec jet
854 if(faRecIndex.GetSize()<ng) faRecIndex.Set(ng); //idx of rec jet assoc to gen jet
857 Printf("New Rec List %d gen index Array %d",nr,faGenIndex.GetSize());
858 Printf("New Gen List %d rec index Array %d",ng,faRecIndex.GetSize());
860 //matching of MC genrator level and reconstructed jets
861 AliAnalysisHelperJetTasks::GetClosestJets(fListJetsGen,ng,fListJets,nr,faGenIndex,faRecIndex,fDebug);
863 // Fill response matrix
864 for(Int_t ir = 0; ir < nr; ir++){
865 AliAODJet *recJet = (AliAODJet*) fListJets->At(ir);
866 Double_t etaJetRec = recJet->Eta();
867 Double_t ptJetRec = recJet->Pt();
868 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
870 if((fJetEtaMin <= etaJetRec) && (etaJetRec <= fJetEtaMax)){
871 Int_t ig = faGenIndex[ir]; //associated generator level jet
872 if(ig >= 0 && ig < ng){
873 if(fDebug > 10) Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
874 AliAODJet *genJet = (AliAODJet*) fListJetsGen->At(ig);
875 Double_t ptJetGen = genJet->Pt();
876 Double_t etaJetGen = genJet->Eta();
878 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
879 if((fJetEtaMin <= etaJetGen) && (etaJetGen <= fJetEtaMax)){
880 fhJetPtGenVsJetPtRec->Fill(ptJetRec, ptJetGen);
883 }//rec jet in eta acceptance
884 }//loop over reconstructed jets
886 }//pointer MC generator jets
887 } //analyze generator level MC
889 // =============== RECONSTRUCTED TRACKS AND JETS ================
891 //Find Hadron trigger
892 TList particleList; //list of tracks
893 Int_t indexTrigg = GetListOfTracks(&particleList); //index of trigger hadron in Particle list
895 if(fIsMC) FillEffHistos(&particleList, &particleListGen); //Fill efficiency histos
898 PostData(1, fOutputList);
899 return; // no trigger track found above 150 MeV/c
902 AliVParticle *triggerHadron = (AliVParticle*) particleList.At(indexTrigg);
904 PostData(1, fOutputList);
909 fh2Ntriggers->Fill(centValue,triggerHadron->Pt()); //trigger pT
911 //Trigger Diagnostics---------------------------------
912 fhTriggerPhi->Fill(triggerHadron->Pt(),RelativePhi(triggerHadron->Phi(),0.0)); //phi -pi,pi
913 fhTriggerEta->Fill(triggerHadron->Pt(),triggerHadron->Eta());
916 Double_t etaJet = 0.0;
917 Double_t pTJet = 0.0;
918 Double_t areaJet = 0.0;
919 Double_t phiJet = 0.0;
922 Int_t indexLeadingJet = -1;
923 Double_t pTLeadingJet = -10.0;
924 Double_t areaLeadingJet = -10.0;
925 //---------- jet loop ---------
926 for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
927 AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
932 if(pTJet==0) continue;
934 if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
935 areaJet = jet->EffectiveAreaCharged();
937 //Jet Diagnostics---------------------------------
938 fhJetPhi->Fill(pTJet, RelativePhi(phiJet,0.0)); //phi -pi,pi
939 fhJetEta->Fill(pTJet, etaJet);
940 if(areaJet >= 0.07) injet++;
941 if(areaJet >= 0.4) injet4++;
942 //--------------------------------------------------
944 Double_t dphi = RelativePhi(triggerHadron->Phi(), phiJet);
946 fhDphiTriggerJet->Fill(dphi); //Input
948 //search for leading jet
949 if(pTJet > pTLeadingJet){
950 indexLeadingJet = ij;
951 pTLeadingJet = pTJet;
952 areaLeadingJet = areaJet;
955 // raw spectra of INCLUSIVE jets
956 //Centrality, pTjet, A
957 Double_t fillraw[] = { centValue,
961 fHJetPtRaw->Fill(fillraw);
964 //Centrality, Dphi=phiTrig-phiJet, pTjet, pTtrigg
965 Double_t filldp[] = { centValue,
970 fHDphiVsJetPtAll->Fill(filldp);
972 // Select back to back trigger - jet pairs
973 if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue;
974 fhDphiTriggerJetAccept->Fill(dphi); //Accepted
977 //Centrality, A, pTjet, pTtrigg
978 Double_t fillspec[] = { centValue,
983 fHJetSpec->Fill(fillspec);
987 if(indexLeadingJet > -1){
988 // raw spectra of LEADING jets
989 //Centrality, pTjet, A
990 Double_t fillleading[] = { centValue,
994 fHLeadingJetPtRaw->Fill(fillleading);
998 //Fill Jet Density In the Event A>0.07
1000 Double_t filldens[]={ (Double_t) particleList.GetEntries(),
1004 fHJetDensity->Fill(filldens);
1007 //Fill Jet Density In the Event A>0.4
1009 Double_t filldens4[]={ (Double_t) particleList.GetEntries(),
1010 injet4/fkAcceptance,
1013 fHJetDensityA4->Fill(filldens4);
1016 PostData(1, fOutputList);
1019 //----------------------------------------------------------------------------
1020 void AliAnalysisTaskJetCorePP::Terminate(const Option_t *)
1022 // Draw result to the screen
1023 // Called once at the end of the query
1025 if(fDebug) printf("AliAnalysisTaskJetCorePP DONE\n");
1026 if(!GetOutputData(1)) return;
1030 //----------------------------------------------------------------------------
1031 Int_t AliAnalysisTaskJetCorePP::GetListOfTracks(TList *list){
1032 //Fill the list of accepted tracks (passed track cut)
1033 //return consecutive index of the hardest ch hadron in the list
1035 AliAODEvent *aodevt = NULL;
1037 if(!fESD) aodevt = fAODIn;
1038 else aodevt = fAODOut;
1040 if(!aodevt) return -1;
1042 Int_t index = -1; //index of the highest particle in the list
1043 Double_t ptmax = -10;
1045 for(Int_t it = 0; it < aodevt->GetNumberOfTracks(); it++){
1046 AliAODTrack *tr = aodevt->GetTrack(it);
1048 //if((fFilterMask > 0) && !(tr->TestFilterBit(fFilterMask))) continue;
1049 if((fFilterMask > 0) && !(tr->IsHybridGlobalConstrainedGlobal())) continue;
1050 if(TMath::Abs((Float_t) tr->Eta()) > fTrackEtaCut) continue;
1051 if(tr->Pt() < fTrackLowPtCut) continue;
1060 if(index>-1){ //check pseudorapidity cut on trigger
1061 AliAODTrack *trigger = (AliAODTrack*) list->At(index);
1062 if(trigger && TMath::Abs((Float_t) trigger->Eta())< fTriggerEtaCut){ return index;}
1069 //----------------------------------------------------------------------------
1071 Double_t AliAnalysisTaskJetCorePP::GetBackgroundInPerpCone(Float_t jetR, Double_t jetPhi, Double_t jetEta, TList* trkList){
1072 //calculate sum of track pT in the cone perpendicular in phi to the jet
1073 //jetR = cone radius
1074 // jetPhi, jetEta = direction of the jet
1075 Int_t numberOfTrks = trkList->GetEntries();
1076 Double_t pTsum = 0.0;
1077 Double_t perpConePhi = jetPhi + TMath::Pi()/2;//perp cone w.r.t. jet in phi
1078 for(Int_t it=0; it<numberOfTrks; it++){
1079 AliVParticle *trk = (AliVParticle*) trkList->At(it);
1080 Double_t dphi = RelativePhi(perpConePhi,trk->Phi());
1081 Double_t deta = trk->Eta()-jetEta;
1083 if( (dphi*dphi + deta*deta)< (jetR*jetR)){
1091 //----------------------------------------------------------------------------
1093 Double_t AliAnalysisTaskJetCorePP::RelativePhi(Double_t mphi,Double_t vphi){
1094 //Get relative azimuthal angle of two particles -pi to pi
1095 if (vphi < -TMath::Pi()) vphi += TMath::TwoPi();
1096 else if (vphi > TMath::Pi()) vphi -= TMath::TwoPi();
1098 if (mphi < -TMath::Pi()) mphi += TMath::TwoPi();
1099 else if (mphi > TMath::Pi()) mphi -= TMath::TwoPi();
1101 Double_t dphi = mphi - vphi;
1102 if (dphi < -TMath::Pi()) dphi += TMath::TwoPi();
1103 else if (dphi > TMath::Pi()) dphi -= TMath::TwoPi();
1105 return dphi;//dphi in [-Pi, Pi]
1109 //----------------------------------------------------------------------------
1110 Bool_t AliAnalysisTaskJetCorePP::SelectMCGenTracks(AliVParticle *trk, TList *trkList, Double_t &ptLeading, Int_t &index, Int_t counter){
1111 //fill track efficiency denominator
1112 if(trk->Pt() < fTrackLowPtCut) return kFALSE;
1113 if(trk->Charge()==0) return kFALSE;
1114 if(TMath::Abs((Float_t) trk->Eta()) > fTrackEtaCut) return kFALSE;
1116 fhPtTrkTruePrimGen->Fill(trk->Pt(),trk->Eta()); //Efficiency denominator
1118 if(ptLeading < trk->Pt()){
1120 ptLeading = trk->Pt();
1125 //----------------------------------------------------------------------------
1126 void AliAnalysisTaskJetCorePP::FillEffHistos(TList *recList, TList *genList){
1127 //fill track efficiency numerator
1128 if(!recList) return;
1129 if(!genList) return;
1130 Int_t nRec = recList->GetEntries();
1131 Int_t nGen = genList->GetEntries();
1132 for(Int_t ir=0; ir<nRec; ir++){
1133 AliAODTrack *trkRec = (AliAODTrack*) recList->At(ir);
1134 if(!trkRec) continue;
1135 Int_t recLabel = TMath::Abs(trkRec->GetLabel());
1136 Bool_t matched = kFALSE;
1138 for(Int_t ig=0; ig<nGen; ig++){
1139 AliVParticle *trkGen = (AliVParticle*) genList->At(ig);
1140 if(!trkGen) continue;
1141 Int_t genLabel = TMath::Abs(trkGen->GetLabel());
1142 if(recLabel==genLabel){
1143 fhPtTrkTruePrimRec->Fill(trkGen->Pt(),trkGen->Eta());
1148 if(!matched) fhPtTrkSecOrFakeRec->Fill(trkRec->Pt(),trkRec->Eta());