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() :
80 fJetBranchNameChargMC(""),
81 fJetBranchNameFullMC(""),
83 fJetBranchNameBgChargMC(""),
86 fListJetsGenFull(0x0),
90 fSystem(0), //pp=0 pPb=1
95 fOfflineTrgMask(AliVEvent::kAny),
106 fTrackLowPtCut(0.15),
108 fHistEvtSelection(0x0),
111 fHJetSpecSubUeMedian(0x0),
112 fHJetSpecSubUeCone(0x0),
115 fHRhoUeMedianVsCone(0x0),
117 //fHJetDensityA4(0x0),
123 fhVertexZAccept(0x0),
125 fhContribVtxAccept(0x0),
126 fhDphiTriggerJet(0x0),
127 fhDphiTriggerJetAccept(0x0),
129 fhCentralityAccept(0x0),
131 //fHLeadingJetPtRaw(0x0),
132 //fHDphiVsJetPtAll(0x0),
133 fhJetPtGenVsJetPtRec(0x0),
134 fhJetPtGenVsJetPtRecSubUeMedian(0x0),
135 fhJetPtGenVsJetPtRecSubUeCone(0x0),
137 fhJetPtSubUeMedianGen(0x0),
138 fhJetPtSubUeConeGen(0x0),
139 fhJetPtGenChargVsJetPtGenFull(0x0),
141 fh2NtriggersGen(0x0),
143 fHJetSpecSubUeMedianGen(0x0),
144 fHJetSpecSubUeConeGen(0x0),
145 fHJetUeMedianGen(0x0),
147 fhPtTrkTruePrimRec(0x0),
148 fhPtTrkTruePrimGen(0x0),
149 fhPtTrkSecOrFakeRec(0x0),
150 fHRhoUeMedianVsConeGen(0x0),
151 fhEntriesToMedian(0x0),
152 fhEntriesToMedianGen(0x0),
153 fhCellAreaToMedian(0x0),
154 fhCellAreaToMedianGen(0x0),
159 fkAcceptance(2.0*TMath::Pi()*1.8),
160 fkDeltaPhiCut(TMath::Pi()-0.8),
166 fh1PtHardTrials(0x0),
169 fEventNumberRangeLow(0),
170 fEventNumberRangeHigh(99),
171 fTriggerPtRangeLow(0.0),
172 fTriggerPtRangeHigh(10000.0),
176 fJetFreeAreaFrac(0.5),
180 fPhiSize(2*TMath::Pi()/fnPhi),
181 fCellArea(fPhiSize*fEtaSize),
184 // default Constructor
187 //---------------------------------------------------------------------
189 AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP(const char *name) :
190 AliAnalysisTaskSE(name),
196 fJetBranchNameChargMC(""),
197 fJetBranchNameFullMC(""),
198 fJetBranchNameBg(""),
199 fJetBranchNameBgChargMC(""),
202 fListJetsGenFull(0x0),
206 fSystem(0), //pp=0 pPb=1
211 fOfflineTrgMask(AliVEvent::kAny),
222 fTrackLowPtCut(0.15),
224 fHistEvtSelection(0x0),
227 fHJetSpecSubUeMedian(0x0),
228 fHJetSpecSubUeCone(0x0),
231 fHRhoUeMedianVsCone(0x0),
233 //fHJetDensityA4(0x0),
239 fhVertexZAccept(0x0),
241 fhContribVtxAccept(0x0),
242 fhDphiTriggerJet(0x0),
243 fhDphiTriggerJetAccept(0x0),
245 fhCentralityAccept(0x0),
247 //fHLeadingJetPtRaw(0x0),
248 //fHDphiVsJetPtAll(0x0),
249 fhJetPtGenVsJetPtRec(0x0),
250 fhJetPtGenVsJetPtRecSubUeMedian(0x0),
251 fhJetPtGenVsJetPtRecSubUeCone(0x0),
253 fhJetPtSubUeMedianGen(0x0),
254 fhJetPtSubUeConeGen(0x0),
255 fhJetPtGenChargVsJetPtGenFull(0x0),
257 fh2NtriggersGen(0x0),
259 fHJetSpecSubUeMedianGen(0x0),
260 fHJetSpecSubUeConeGen(0x0),
261 fHJetUeMedianGen(0x0),
263 fhPtTrkTruePrimRec(0x0),
264 fhPtTrkTruePrimGen(0x0),
265 fhPtTrkSecOrFakeRec(0x0),
266 fHRhoUeMedianVsConeGen(0x0),
267 fhEntriesToMedian(0x0),
268 fhEntriesToMedianGen(0x0),
269 fhCellAreaToMedian(0x0),
270 fhCellAreaToMedianGen(0x0),
275 fkAcceptance(2.0*TMath::Pi()*1.8),
276 fkDeltaPhiCut(TMath::Pi()-0.8),
282 fh1PtHardTrials(0x0),
285 fEventNumberRangeLow(0),
286 fEventNumberRangeHigh(99),
287 fTriggerPtRangeLow(0.0),
288 fTriggerPtRangeHigh(10000.0),
292 fJetFreeAreaFrac(0.5),
296 fPhiSize(2*TMath::Pi()/fnPhi),
297 fCellArea(fPhiSize*fEtaSize),
302 DefineOutput(1, TList::Class());
305 //--------------------------------------------------------------
306 AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP(const AliAnalysisTaskJetCorePP& a):
307 AliAnalysisTaskSE(a.GetName()),
311 fAODExtension(a.fAODExtension),
312 fJetBranchName(a.fJetBranchName),
313 fJetBranchNameChargMC(a.fJetBranchNameChargMC),
314 fJetBranchNameFullMC(a.fJetBranchNameFullMC),
315 fJetBranchNameBg(a.fJetBranchNameBg),
316 fJetBranchNameBgChargMC(a.fJetBranchNameBgChargMC),
317 fListJets(a.fListJets),
318 fListJetsGen(a.fListJetsGen),
319 fListJetsGenFull(a.fListJetsGenFull),
320 fListJetsBg(a.fListJetsBg),
321 fListJetsBgGen(a.fListJetsBgGen),
322 fNonStdFile(a.fNonStdFile),
324 fJetParamR(a.fJetParamR),
325 fBgJetParamR(a.fBgJetParamR),
326 fBgMaxJetPt(a.fBgMaxJetPt),
327 fBgConeR(a.fBgConeR),
328 fOfflineTrgMask(a.fOfflineTrgMask),
329 fMinContribVtx(a.fMinContribVtx),
330 fVtxZMin(a.fVtxZMin),
331 fVtxZMax(a.fVtxZMax),
332 fFilterMask(a.fFilterMask),
333 fCentMin(a.fCentMin),
334 fCentMax(a.fCentMax),
335 fJetEtaMin(a.fJetEtaMin),
336 fJetEtaMax(a.fJetEtaMax),
337 fTriggerEtaCut(a.fTriggerEtaCut),
338 fTrackEtaCut(a.fTrackEtaCut),
339 fTrackLowPtCut(a.fTrackLowPtCut),
340 fOutputList(a.fOutputList),
341 fHistEvtSelection(a.fHistEvtSelection),
342 fh2Ntriggers(a.fh2Ntriggers),
343 fHJetSpec(a.fHJetSpec),
344 fHJetSpecSubUeMedian(a.fHJetSpecSubUeMedian),
345 fHJetSpecSubUeCone(a.fHJetSpecSubUeCone),
346 fHJetUeMedian(a.fHJetUeMedian),
347 fHJetUeCone(a.fHJetUeCone),
348 fHRhoUeMedianVsCone(a.fHRhoUeMedianVsCone),
349 //fHJetDensity(a.fHJetDensity),
350 //fHJetDensityA4(a.fHJetDensityA4),
351 fhJetPhi(a.fhJetPhi),
352 fhTriggerPhi(a.fhTriggerPhi),
353 fhJetEta(a.fhJetEta),
354 fhTriggerEta(a.fhTriggerEta),
355 fhVertexZ(a.fhVertexZ),
356 fhVertexZAccept(a.fhVertexZAccept),
357 fhContribVtx(a.fhContribVtx),
358 fhContribVtxAccept(a.fhContribVtxAccept),
359 fhDphiTriggerJet(a.fhDphiTriggerJet),
360 fhDphiTriggerJetAccept(a.fhDphiTriggerJetAccept),
361 fhCentrality(a.fhCentrality),
362 fhCentralityAccept(a.fhCentralityAccept),
363 //fHJetPtRaw(a.fHJetPtRaw),
364 //fHLeadingJetPtRaw(a.fHLeadingJetPtRaw),
365 //fHDphiVsJetPtAll(a.fHDphiVsJetPtAll),
366 fhJetPtGenVsJetPtRec(a.fhJetPtGenVsJetPtRec),
367 fhJetPtGenVsJetPtRecSubUeMedian(a.fhJetPtGenVsJetPtRecSubUeMedian),
368 fhJetPtGenVsJetPtRecSubUeCone(a.fhJetPtGenVsJetPtRecSubUeCone),
369 fhJetPtGen(a.fhJetPtGen),
370 fhJetPtSubUeMedianGen(a.fhJetPtSubUeMedianGen),
371 fhJetPtSubUeConeGen(a.fhJetPtSubUeConeGen),
372 fhJetPtGenChargVsJetPtGenFull(a.fhJetPtGenChargVsJetPtGenFull),
373 fhJetPtGenFull(a.fhJetPtGenFull),
374 fh2NtriggersGen(a.fh2NtriggersGen),
375 fHJetSpecGen(a.fHJetSpecGen),
376 fHJetSpecSubUeMedianGen(a.fHJetSpecSubUeMedianGen),
377 fHJetSpecSubUeConeGen(a.fHJetSpecSubUeConeGen),
378 fHJetUeMedianGen(a.fHJetUeMedianGen),
379 fHJetUeConeGen(a.fHJetUeConeGen),
380 fhPtTrkTruePrimRec(a.fhPtTrkTruePrimRec),
381 fhPtTrkTruePrimGen(a.fhPtTrkTruePrimGen),
382 fhPtTrkSecOrFakeRec(a.fhPtTrkSecOrFakeRec),
383 fHRhoUeMedianVsConeGen(a.fHRhoUeMedianVsConeGen),
384 fhEntriesToMedian(a.fhEntriesToMedian),
385 fhEntriesToMedianGen(a.fhEntriesToMedianGen),
386 fhCellAreaToMedian(a.fhCellAreaToMedian),
387 fhCellAreaToMedianGen(a.fhCellAreaToMedianGen),
388 fIsChargedMC(a.fIsChargedMC),
389 fIsFullMC(a.fIsFullMC),
390 faGenIndex(a.faGenIndex),
391 faRecIndex(a.faRecIndex),
392 fkAcceptance(a.fkAcceptance),
393 fkDeltaPhiCut(a.fkDeltaPhiCut),
395 fh1Trials(a.fh1Trials),
396 fh1AvgTrials(a.fh1AvgTrials),
397 fh1PtHard(a.fh1PtHard),
398 fh1PtHardNoW(a.fh1PtHardNoW),
399 fh1PtHardTrials(a.fh1PtHardTrials),
400 fAvgTrials(a.fAvgTrials),
401 fHardest(a.fHardest),
402 fEventNumberRangeLow(a.fEventNumberRangeLow),
403 fEventNumberRangeHigh(a.fEventNumberRangeHigh),
404 fTriggerPtRangeLow(a.fTriggerPtRangeLow),
405 fTriggerPtRangeHigh(a.fTriggerPtRangeHigh),
406 fFillRespMx(a.fFillRespMx),
408 fnTrials(a.fnTrials),
409 fJetFreeAreaFrac(a.fJetFreeAreaFrac),
412 fEtaSize(a.fEtaSize),
413 fPhiSize(a.fPhiSize),
414 fCellArea(a.fCellArea),
415 fSafetyMargin(a.fSafetyMargin)
420 //--------------------------------------------------------------
422 AliAnalysisTaskJetCorePP& AliAnalysisTaskJetCorePP::operator = (const AliAnalysisTaskJetCorePP& a){
423 // assignment operator
424 this->~AliAnalysisTaskJetCorePP();
425 new(this) AliAnalysisTaskJetCorePP(a);
428 //--------------------------------------------------------------
430 AliAnalysisTaskJetCorePP::~AliAnalysisTaskJetCorePP()
435 delete fListJetsGenFull;
437 delete fListJetsBgGen;
438 delete fOutputList; // ????
442 //--------------------------------------------------------------
445 Bool_t AliAnalysisTaskJetCorePP::Notify()
447 //Implemented Notify() to read the cross sections
448 //and number of trials from pyxsec.root
449 //inspired by AliAnalysisTaskJetSpectrum2::Notify()
450 if(!fIsChargedMC) return kFALSE;
452 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
454 if(fDebug>1) AliError("ESD not available");
455 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
458 fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
461 if(fNonStdFile.Length()!=0){
462 // case that we have an AOD extension we can fetch the jets from the extended output
463 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
464 fAODExtension = aodH ? aodH->GetExtension(fNonStdFile.Data()) : 0;
466 if(fDebug>1) Printf("AODExtension found for %s",fNonStdFile.Data());
470 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
471 Float_t xsection = 0;
476 TFile *curfile = tree->GetCurrentFile();
478 Error("Notify","No current file");
481 if(!fh1Xsec || !fh1Trials){
482 Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
485 AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
486 fh1Xsec->Fill("<#sigma>",xsection);
487 // construct a poor man average trials
488 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
489 if(ftrials>=nEntries && nEntries>0.) fAvgTrials = ftrials/nEntries;
490 fh1Trials->Fill("#sum{ntrials}",ftrials);
493 if(fDebug)Printf("Reading File %s",fInputHandler->GetTree()->GetCurrentFile()->GetName());
497 //--------------------------------------------------------------
499 void AliAnalysisTaskJetCorePP::Init()
501 // check for jet branches
502 if(!strlen(fJetBranchName.Data())){
503 AliError("Jet branch name not set.");
507 //--------------------------------------------------------------
509 void AliAnalysisTaskJetCorePP::UserCreateOutputObjects()
511 // Create histograms and initilize variables
512 fSafetyMargin = fBgConeR*fBgConeR /(fBgJetParamR*fBgJetParamR);
515 fListJets = new TList(); //reconstructed level antikt jets
516 fListJetsBg = new TList(); //reconstructed jets to be removed from bg
518 fIsChargedMC = (fJetBranchNameChargMC.Length()>0) ? kTRUE : kFALSE;
519 fIsFullMC = (fJetBranchNameFullMC.Length()>0) ? kTRUE : kFALSE;
521 fRandom = new TRandom3(0);
524 fListJetsGen = new TList(); //generator level charged antikt jets
525 fListJetsBgGen = new TList(); //generator level jets to be removed from bg
528 fListJetsGenFull = new TList(); //generator level full jets
531 if(!fOutputList) fOutputList = new TList();
532 fOutputList->SetOwner(kTRUE);
534 Bool_t oldStatus = TH1::AddDirectoryStatus();
535 TH1::AddDirectory(kFALSE);
537 fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
538 fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
539 fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
540 fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
541 fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
542 fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
543 fHistEvtSelection->GetXaxis()->SetBinLabel(6,"event number (rejected)");
545 fOutputList->Add(fHistEvtSelection);
547 Int_t nBinsCentrality = (fSystem==0) ? 1 : 10; // pp=1 else 10
548 //trigger pt spectrum (reconstructed)
549 fh2Ntriggers = new TH2F("fh2Ntriggers","# of triggers",
550 nBinsCentrality,0.0,100.0,50,0.0,50.0);
551 fOutputList->Add(fh2Ntriggers);
553 //Centrality, A, pTjet, pTtrigg, dphi
554 const Int_t dimSpec = 5;
555 const Int_t nBinsSpec[dimSpec] = {nBinsCentrality, 50, 110, 50, TMath::Nint(10*(TMath::Pi()-fkDeltaPhiCut))};
556 const Double_t lowBinSpec[dimSpec] = {0.0, 0.0, -20, 0.0, fkDeltaPhiCut};
557 const Double_t hiBinSpec[dimSpec] = {100.0, 1.0, 200.0,50.0, TMath::Pi()};
558 fHJetSpec = new THnSparseF("fHJetSpec",
559 "Recoil jet spectrum [cent,A,pTjet,pTtrig,dphi]",
560 dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
561 fOutputList->Add(fHJetSpec);
563 //background estimated as median of kT jets
564 fHJetSpecSubUeMedian = (THnSparseF*) fHJetSpec->Clone("fHJetSpecSubUeMedian");
565 fHJetSpecSubUeMedian->SetTitle("Recoil jet spectrum [cent,A,pTjet-pTUe,pTtrig,dphi]");
566 fOutputList->Add(fHJetSpecSubUeMedian);
567 //background estimated as weighted median of kT jets ala Cone
568 fHJetSpecSubUeCone = (THnSparseF*) fHJetSpec->Clone("fHJetSpecSubUeCone");
569 fHJetSpecSubUeCone->SetTitle("Recoil jet spectrum [cent,A,pTjet-pTUe,pTtrig,dphi]");
570 fOutputList->Add(fHJetSpecSubUeCone);
574 //------------------- HISTOS FOR DIAGNOSTIC ----------------------
575 //A, pTjet, pTjet-pTUe, pTUe, rhoUe bg estimate from kT median
576 const Int_t dimSpecMed = 5;
577 const Int_t nBinsSpecMed[dimSpecMed] = {25, 50, 60, 20, 20};
578 const Double_t lowBinSpecMed[dimSpecMed] = {0.0, 0.0, -20.0, 0.0, 0.0};
579 const Double_t hiBinSpecMed[dimSpecMed] = {1.0, 100.0, 100.0, 10.0, 20.0};
580 fHJetUeMedian = new THnSparseF("fHJetUeMedian",
581 "Recoil jet spectrum [A,pTjet,pTjet-pTUe, pTUe, rhoUe]",
582 dimSpecMed, nBinsSpecMed, lowBinSpecMed, hiBinSpecMed);
583 fOutputList->Add(fHJetUeMedian);
585 //A, pTjet, pTjet-pTUe, pTUe, rhoUe bg estimate from kT median Cone
586 fHJetUeCone = (THnSparseF*) fHJetUeMedian->Clone("fHJetUeCone");
587 fHJetUeCone->SetTitle("Recoil jet spectrum [A,pTjet,pTjet-pTUe, pTUe, rhoUe]");
588 fOutputList->Add(fHJetUeCone);
590 //rho bacground reconstructed data
591 const Int_t dimRho = 2;
592 const Int_t nBinsRho[dimRho] = {50 , 50};
593 const Double_t lowBinRho[dimRho] = {0.0 , 0.0};
594 const Double_t hiBinRho[dimRho] = {20.0 , 20.0};
596 fHRhoUeMedianVsCone = new THnSparseF("hRhoUeMedianVsCone","[Rho Cone, Rho Median]",
597 dimRho, nBinsRho, lowBinRho, hiBinRho);
598 fOutputList->Add(fHRhoUeMedianVsCone);
600 //Jet number density histos [Trk Mult, jet density, pT trigger]
601 /*const Int_t dimJetDens = 3;
602 const Int_t nBinsJetDens[dimJetDens] = {100, 100, 10};
603 const Double_t lowBinJetDens[dimJetDens] = {0.0, 0.0, 0.0};
604 const Double_t hiBinJetDens[dimJetDens] = {500.0, 5.0, 50.0};
606 fHJetDensity = new THnSparseF("fHJetDensity","Jet dens vs trk mult A>0.07",
607 dimJetDens,nBinsJetDens,lowBinJetDens,hiBinJetDens);
609 fHJetDensityA4 =new THnSparseF("fHJetDensityA4","Jet dens vs trk mult A>0.4",
610 dimJetDens,nBinsJetDens,lowBinJetDens,hiBinJetDens);
612 fOutputList->Add(fHJetDensity);
613 fOutputList->Add(fHJetDensityA4);
616 //inclusive azimuthal and pseudorapidity histograms
617 fhJetPhi = new TH2D("fhJetPhi","Azim dist jets vs pTjet",
618 50, 0, 100, 50,-TMath::Pi(),TMath::Pi());
619 fhTriggerPhi= new TH2D("fhTriggerPhi","azim dist trig had vs pTtrigg",
620 50, 0, 50, 50,-TMath::Pi(),TMath::Pi());
621 fhJetEta = new TH2D("fhJetEta","Eta dist jets vs pTjet",
622 50,0, 100, 40,-0.9,0.9);
623 fhTriggerEta = new TH2D("fhTriggerEta","Eta dist trig had vs pTtrigg",
624 50, 0, 50, 40,-0.9,0.9);
626 fhVertexZ = new TH1D("fhVertexZ","z vertex",40,-20,20);
627 fhVertexZAccept = new TH1D("fhVertexZAccept","z vertex after cut",40,-20,20);
628 fhContribVtx = new TH1D("fhContribVtx","contrib to vtx",200,0,200);
629 fhContribVtxAccept = new TH1D("fhContribVtxAccept","contrib to vtx after cut",200,0,200);
630 fhDphiTriggerJet = new TH1D("fhDphiTriggerJet","Deltaphi trig-jet",50, -TMath::Pi(),TMath::Pi());
631 fhDphiTriggerJetAccept = new TH1D("fhDphiTriggerJetAccept","Deltaphi trig-jet after cut",50, -TMath::Pi(),TMath::Pi());
632 fhCentrality = new TH1D("fhCentrality","Centrality",20,0,100);
633 fhCentralityAccept = new TH1D("fhCentralityAccept","Centrality after cut",20,0,100);
634 fhEntriesToMedian = new TH1D("fhEntriesToMedian","fhEntriesToMedian",30,0,30);
635 fhCellAreaToMedian = new TH1D("fhCellAreaToMedian", "fhCellAreaToMedian", 75,0,1.5),
637 fOutputList->Add(fhJetPhi);
638 fOutputList->Add(fhTriggerPhi);
639 fOutputList->Add(fhJetEta);
640 fOutputList->Add(fhTriggerEta);
641 fOutputList->Add(fhVertexZ);
642 fOutputList->Add(fhVertexZAccept);
643 fOutputList->Add(fhContribVtx);
644 fOutputList->Add(fhContribVtxAccept);
645 fOutputList->Add(fhDphiTriggerJet);
646 fOutputList->Add(fhDphiTriggerJetAccept);
647 fOutputList->Add(fhCentrality);
648 fOutputList->Add(fhCentralityAccept);
649 fOutputList->Add(fhEntriesToMedian);
650 fOutputList->Add(fhCellAreaToMedian);
651 // raw spectra of INCLUSIVE jets
652 //Centrality, pTjet, A
653 /*const Int_t dimRaw = 3;
654 const Int_t nBinsRaw[dimRaw] = {nBinsCentrality, 50, 100};
655 const Double_t lowBinRaw[dimRaw] = {0.0, 0.0, 0.0};
656 const Double_t hiBinRaw[dimRaw] = {100.0, 100, 1.0};
657 fHJetPtRaw = new THnSparseF("fHJetPtRaw",
658 "Incl. jet spectrum [cent,pTjet,A]",
659 dimRaw,nBinsRaw,lowBinRaw,hiBinRaw);
660 fOutputList->Add(fHJetPtRaw);
662 // raw spectra of LEADING jets
663 //Centrality, pTjet, A
664 fHLeadingJetPtRaw = new THnSparseF("fHLeadingJetPtRaw",
665 "Leading jet spectrum [cent,pTjet,A]",
666 dimRaw,nBinsRaw,lowBinRaw,hiBinRaw);
667 fOutputList->Add(fHLeadingJetPtRaw);
669 // Dphi versus pT jet
670 //Centrality, Dphi=phiTrig-phiJet, pTjet, pTtrigg
671 const Int_t dimDp = 4;
672 const Int_t nBinsDp[dimDp] = {nBinsCentrality, 50, 50, 50};
673 const Double_t lowBinDp[dimDp] = {0.0, -TMath::Pi(), 0.0, 0.0};
674 const Double_t hiBinDp[dimDp] = {100.0, TMath::Pi(), 100.0, 100.0};
675 fHDphiVsJetPtAll = new THnSparseF("fHDphiVsJetPtAll",
676 "Dphi vs jet pT [cent,Dphi,pTjet,pTtrigg]",
677 dimDp,nBinsDp,lowBinDp,hiBinDp);
678 fOutputList->Add(fHDphiVsJetPtAll);
681 //analyze MC generator level
684 //Fill response matrix only once
685 fhJetPtGenVsJetPtRec = new TH2D("fhJetPtGenVsJetPtRec","JetPtGenVsJetPtRec", 100,0,200, 100,0,200);
686 fOutputList->Add(fhJetPtGenVsJetPtRec); //gen MC charg jet pt spectrum versus rec charged jet pt spectrum
688 fhJetPtGenVsJetPtRecSubUeMedian = new TH2D("fhJetPtGenVsJetPtRecSubUeMedian","fhJetPtGenVsJetPtRecSubUeMedian", 110,-20,200, 110,-20,200);
689 fOutputList->Add(fhJetPtGenVsJetPtRecSubUeMedian); // with kT median bg subtr
691 fhJetPtGenVsJetPtRecSubUeCone=(TH2D*)fhJetPtGenVsJetPtRecSubUeMedian->Clone("fhJetPtGenVsJetPtRecSubUeCone");
692 fhJetPtGenVsJetPtRecSubUeCone->SetTitle("fhJetPtGenVsJetPtRecSubUeCone");
693 fOutputList->Add(fhJetPtGenVsJetPtRecSubUeCone); // with weighted kT median bg subtr
695 fhJetPtGen = new TH1D("fhJetPtGen","Jet Pt (MC Gen)",100,0,200); //MC generator charged jet pt spectrum
696 fOutputList->Add(fhJetPtGen);
698 fhJetPtSubUeMedianGen = new TH1D("fhJetPtSubUeMedianGen","Jet Pt - UE pT (MC Gen)",110,-20,200);
699 fOutputList->Add(fhJetPtSubUeMedianGen); // with kT median bg subtr
701 fhJetPtSubUeConeGen = (TH1D*) fhJetPtSubUeMedianGen->Clone("fhJetPtSubUeConeGen");
702 fOutputList->Add(fhJetPtSubUeConeGen); // with weighted kT median bg subtr
706 fhJetPtGenChargVsJetPtGenFull = new TH2D("fhJetPtGenChargVsJetPtGenFull","fhJetPtGenChargVsJetPtGenFull", 100,0,200, 100,0,200);
707 fOutputList->Add(fhJetPtGenChargVsJetPtGenFull); //gen full MC jet pt versus gen charged jet MC pt
709 fhJetPtGenFull = new TH1D("fhJetPtGenFull","Jet Pt (MC Full jets Gen)",100,0,200); //MC generator full jet pt spectrum
710 fOutputList->Add(fhJetPtGenFull);
714 fh2NtriggersGen = (TH2F*) fh2Ntriggers->Clone("fh2NtriggersGen");
715 fh2NtriggersGen->SetTitle(Form("%s Gen MC",fh2Ntriggers->GetTitle()));
716 fOutputList->Add(fh2NtriggersGen);
718 //Centrality, A, pT, pTtrigg
719 fHJetSpecGen = (THnSparseF*) fHJetSpec->Clone("fHJetSpecGen");
720 fHJetSpecGen->SetTitle(Form("%s Gen MC",fHJetSpec->GetTitle()));
721 fOutputList->Add(fHJetSpecGen);
723 fHJetSpecSubUeMedianGen = (THnSparseF*) fHJetSpecSubUeMedian->Clone("fHJetSpecSubUeMedianGen");
724 fHJetSpecSubUeMedianGen->SetTitle(Form("%s Gen MC",fHJetSpecSubUeMedian->GetTitle()));
725 fOutputList->Add(fHJetSpecSubUeMedianGen);
727 fHJetSpecSubUeConeGen = (THnSparseF*) fHJetSpecSubUeCone->Clone("fHJetSpecSubUeConeGen");
728 fHJetSpecSubUeConeGen->SetTitle(Form("%s Gen MC",fHJetSpecSubUeCone->GetTitle()));
729 fOutputList->Add(fHJetSpecSubUeConeGen);
731 fHJetUeMedianGen = (THnSparseF*) fHJetUeMedian->Clone("fHJetUeMedianGen");
732 fHJetUeMedianGen->SetTitle(Form("%s Gen MC", fHJetUeMedian->GetTitle()));
733 fOutputList->Add(fHJetUeMedianGen);
735 fHJetUeConeGen = (THnSparseF*) fHJetUeCone->Clone("fHJetUeConeGen");
736 fHJetUeConeGen->SetTitle(Form("%s Gen MC", fHJetUeCone->GetTitle()));
737 fOutputList->Add(fHJetUeConeGen);
740 //track efficiency/contamination histograms eta versus pT
741 fhPtTrkTruePrimRec = new TH2D("fhPtTrkTruePrimRec","PtTrkTruePrimRec",100,0,20,18,-0.9,0.9);
742 fOutputList->Add(fhPtTrkTruePrimRec);
744 fhPtTrkTruePrimGen = (TH2D*) fhPtTrkTruePrimRec->Clone("fhPtTrkTruePrimGen");
745 fhPtTrkTruePrimGen->SetTitle("PtTrkTruePrimGen");
746 fOutputList->Add(fhPtTrkTruePrimGen);
748 fhPtTrkSecOrFakeRec = (TH2D*) fhPtTrkTruePrimRec->Clone("fhPtTrkSecOrFakeRec");
749 fhPtTrkSecOrFakeRec->SetTitle("PtTrkSecOrFakeRec");
750 fOutputList->Add(fhPtTrkSecOrFakeRec);
753 fHRhoUeMedianVsConeGen = (THnSparseF*) fHRhoUeMedianVsCone->Clone("hRhoUeMedianVsConeGen");
754 fHRhoUeMedianVsConeGen->SetTitle(Form("%s Gen MC", fHRhoUeMedianVsCone->GetTitle()));
755 fOutputList->Add(fHRhoUeMedianVsConeGen);
757 fhEntriesToMedianGen = (TH1D*) fhEntriesToMedian->Clone("fhEntriesToMedianGen");
758 fhEntriesToMedianGen->SetTitle(Form("%s Gen MC", fhEntriesToMedian->GetTitle()));
759 fOutputList->Add(fhEntriesToMedianGen);
761 fhCellAreaToMedianGen = (TH1D*) fhCellAreaToMedian->Clone("fhCellAreaToMedianGen");
762 fhCellAreaToMedianGen->SetTitle(Form("%s Gen MC", fhCellAreaToMedian->GetTitle()));
763 fOutputList->Add(fhCellAreaToMedianGen);
765 //-------------------------------------
767 const Int_t nBinPt = 150;
768 Double_t binLimitsPt[nBinPt+1];
769 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
771 binLimitsPt[iPt] = -50.0;
773 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.;
777 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
778 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
779 fOutputList->Add(fh1Xsec);
780 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
781 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
782 fOutputList->Add(fh1Trials);
783 fh1AvgTrials = new TH1F("fh1AvgTrials","trials root file",1,0,1);
784 fh1AvgTrials->GetXaxis()->SetBinLabel(1,"#sum{avg ntrials}");
785 fOutputList->Add(fh1AvgTrials);
786 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
787 fOutputList->Add(fh1PtHard);
788 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
789 fOutputList->Add(fh1PtHardNoW);
790 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
791 fOutputList->Add(fh1PtHardTrials);
794 // =========== Switch on Sumw2 for all histos ===========
795 for(Int_t i=0; i<fOutputList->GetEntries(); i++){
796 TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
801 THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
806 TH1::AddDirectory(oldStatus);
808 PostData(1, fOutputList);
811 //--------------------------------------------------------------------
813 void AliAnalysisTaskJetCorePP::UserExec(Option_t *)
816 Double_t eventW = 1.0;
817 Double_t ptHard = 0.0;
818 Double_t nTrials = 1.0; // Trials for MC trigger
819 if(fIsChargedMC) fh1AvgTrials->Fill("#sum{avg ntrials}",fAvgTrials);
821 if(TMath::Abs((Float_t) fJetParamR)<0.00001){
822 AliError("ANTIKT Cone radius is set to zero.");
825 if(TMath::Abs((Float_t) fBgJetParamR)<0.00001){
826 AliError("KT Cone radius is set to zero.");
829 if(!strlen(fJetBranchName.Data())){
830 AliError("Jet branch name not set.");
833 if(!strlen(fJetBranchNameBg.Data())){
834 AliError("Jet bg branch name not set.");
839 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
841 if(fDebug>1) AliError("ESD not available");
842 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
845 fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
846 AliAODEvent* aod = NULL;
847 // take all other information from the aod we take the tracks from
849 if(!fESD) aod = fAODIn;
853 if(fNonStdFile.Length()!=0){
854 // case that we have an AOD extension we can fetch the jets from the extended output
855 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
856 fAODExtension = aodH ? aodH->GetExtension(fNonStdFile.Data()) : 0;
858 if(fDebug>1) Printf("AODExtension found for %s",fNonStdFile.Data());
862 // ----------------- EVENT SELECTION --------------------------
863 fHistEvtSelection->Fill(1); // number of events before event selection
866 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
867 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
869 if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
870 if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
871 fHistEvtSelection->Fill(2);
872 PostData(1, fOutputList);
878 if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
879 fHistEvtSelection->Fill(3);
880 PostData(1, fOutputList);
885 AliAODVertex* primVtx = aod->GetPrimaryVertex();
888 if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
889 fHistEvtSelection->Fill(3);
890 PostData(1, fOutputList);
894 Int_t nTracksPrim = primVtx->GetNContributors();
895 Float_t vtxz = primVtx->GetZ();
897 fhContribVtx->Fill(nTracksPrim);
898 if( nTracksPrim > 0 ) fhVertexZ->Fill(vtxz);
900 if((nTracksPrim < fMinContribVtx) ||
903 if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",
904 (char*)__FILE__,__LINE__,vtxz);
905 fHistEvtSelection->Fill(3);
906 PostData(1, fOutputList);
910 fhContribVtxAccept->Fill(nTracksPrim);
911 fhVertexZAccept->Fill(vtxz);
913 //FK// No event class selection imposed
914 // event class selection (from jet helper task)
915 //Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
916 //if(fDebug) Printf("Event class %d", eventClass);
917 //if(eventClass < fEvtClassMin || eventClass > fEvtClassMax){
918 // fHistEvtSelection->Fill(4);
919 // PostData(1, fOutputList);
923 //------------------ CENTRALITY SELECTION ---------------
924 AliCentrality *cent = 0x0;
925 Double_t centValue = 0.0;
926 if(fSystem){ //fSystem=0 for pp, fSystem=1 for pPb
928 cent = fESD->GetCentrality();
929 if(cent) centValue = cent->GetCentralityPercentile("V0M");
931 centValue = aod->GetHeader()->GetCentrality();
933 if(fDebug) printf("centrality: %f\n", centValue);
935 fhCentrality->Fill(centValue);
937 if(centValue < fCentMin || centValue > fCentMax){
938 fHistEvtSelection->Fill(4);
939 PostData(1, fOutputList);
943 fhCentralityAccept->Fill( centValue );
947 //-----------------select disjunct event subsamples ----------------
948 Int_t eventnum = aod->GetHeader()->GetEventNumberESDFile();
949 Int_t lastdigit = eventnum % 10;
950 if(!(fEventNumberRangeLow<=lastdigit && lastdigit<=fEventNumberRangeHigh)){
951 fHistEvtSelection->Fill(5);
952 PostData(1, fOutputList);
956 if(fDebug) std::cout<<" ACCEPTED EVENT "<<endl;
957 fHistEvtSelection->Fill(0); // accepted events
958 // ==================== end event selection ============================
960 Double_t tmpArrayFive[5];
962 //=============== Reconstructed antikt jets ===============
963 ReadTClonesArray(fJetBranchName.Data() , fListJets);
964 ReadTClonesArray(fJetBranchNameBg.Data(), fListJetsBg);
966 //============ Estimate background in reconstructed events ===========
967 Double_t rhoFromCellMedian=0.0, rhoCone=0.0;
969 //Find Hadron trigger and Estimate rho from cone
970 TList particleList; //list of tracks
971 Int_t indexTrigg = GetListOfTracks(&particleList); //index of trigger hadron in Particle list
972 EstimateBgCone(fListJets, &particleList, rhoCone);
974 //Estimate rho from cell median minus jets
975 EstimateBgRhoMedian(fListJetsBg, &particleList, rhoFromCellMedian,0); //real data
976 //fListJetsBg->Clear(); //this list is further not needed
978 //============== analyze generator level MC ================
979 TList particleListGen; //list of tracks in MC
982 //================= generated charged antikt jets ================
983 ReadTClonesArray(fJetBranchNameChargMC.Data(), fListJetsGen);
984 ReadTClonesArray(fJetBranchNameBgChargMC.Data(), fListJetsBgGen);
986 if(fIsFullMC){ //generated full jets
987 ReadTClonesArray(fJetBranchNameFullMC.Data(), fListJetsGenFull);
990 //========================================================
991 //serarch for charged trigger at the MC generator level
993 Int_t indexTriggGen = -1;
994 Double_t ptTriggGen = -1;
995 Int_t iCounterGen = 0; //number of entries in particleListGen array
996 Int_t triggersMC[200];//list of trigger candidates
997 Int_t ntriggersMC = 0; //index in triggers array
1001 AliMCEvent* mcEvent = MCEvent();
1003 PostData(1, fOutputList);
1007 AliGenPythiaEventHeader *pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
1008 if(pythiaGenHeader){
1009 nTrials = pythiaGenHeader->Trials();
1010 ptHard = pythiaGenHeader->GetPtHard();
1011 fh1PtHard->Fill(ptHard,eventW);
1012 fh1PtHardNoW->Fill(ptHard,1);
1013 fh1PtHardTrials->Fill(ptHard,nTrials);
1016 for(Int_t it = 0; it < mcEvent->GetNumberOfTracks(); it++){
1017 if(!mcEvent->IsPhysicalPrimary(it)) continue;
1018 AliMCParticle* part = (AliMCParticle*) mcEvent->GetTrack(it);
1020 if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){
1022 if(fHardest==0 && ntriggersMC<200){//single inclusive trigger
1023 if(indexTriggGen > -1){//trigger candidate was found
1024 triggersMC[ntriggersMC] = indexTriggGen;
1029 iCounterGen++;//index in particleListGen array
1034 PostData(1, fOutputList);
1037 TClonesArray *tca = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(AliAODMCParticle::StdBranchName()));
1039 PostData(1, fOutputList);
1042 for(Int_t it = 0; it < tca->GetEntriesFast(); it++){
1043 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1045 if(!part->IsPhysicalPrimary()) continue;
1046 if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){
1048 if(fHardest==0 && ntriggersMC<200){//single inclusive trigger
1049 if(indexTriggGen > -1){ //trigger candidater was found
1050 triggersMC[ntriggersMC] = indexTriggGen;
1055 iCounterGen++;//number of entries in particleListGen array
1060 //============== Estimate bg in generated events ==============
1061 Double_t rhoFromCellMedianGen=0.0, rhoConeGen=0.0;
1063 //Estimate rho from cone
1064 EstimateBgCone(fListJetsGen, &particleListGen, rhoConeGen);
1066 //Estimate rho from cell median minus jets
1067 EstimateBgRhoMedian(fListJetsBgGen, &particleListGen, rhoFromCellMedianGen,1);//mc data
1068 //fListJetsBgGen->Clear();
1070 //============ Generator trigger+jet ==================
1071 if(fHardest==0){ //single inclusive trigger
1072 if(ntriggersMC>0){ //there is at least one trigger
1073 Int_t rnd = fRandom->Integer(ntriggersMC); //0 to ntriggers-1
1074 indexTriggGen = triggersMC[rnd];
1076 indexTriggGen = -1; //trigger not found
1081 Int_t ilowGen = (fHardest==0 || fHardest==1) ? indexTriggGen : 0;
1082 Int_t ihighGen = (fHardest==0 || fHardest==1) ? indexTriggGen+1 : particleListGen.GetEntries();
1083 Bool_t fillOnceGen = kTRUE;
1086 for(Int_t igen= ilowGen; igen< ihighGen; igen++){ //loop over possible trigger
1087 indexTriggGen = igen; //trigger hadron
1089 if(indexTriggGen == -1) continue;
1090 AliVParticle* triggerGen=(AliVParticle*) particleListGen.At(indexTriggGen);
1091 if(!triggerGen) continue;
1094 if(triggerGen->Pt() < 6.0) continue; //all hadrons pt>6
1096 if(TMath::Abs((Float_t) triggerGen->Eta()) > fTriggerEtaCut) continue;
1098 ptTriggGen = triggerGen->Pt(); //count triggers
1099 fh2NtriggersGen->Fill(centValue, ptTriggGen);
1101 //Count jets and trigger-jet pairs at MC generator level
1102 for(Int_t ij=0; ij<fListJetsGen->GetEntries(); ij++){
1103 AliAODJet* jet = (AliAODJet*)(fListJetsGen->At(ij));
1105 Double_t etaJetGen = jet->Eta();
1106 Double_t areaJetGen = jet->EffectiveAreaCharged();
1108 if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){
1110 Double_t dphi = RelativePhi(triggerGen->Phi(), jet->Phi());
1111 if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue;
1113 //Centrality, A, pT, pTtrigg
1114 tmpArrayFive[0] = centValue;
1115 tmpArrayFive[1] = areaJetGen;
1116 tmpArrayFive[2] = jet->Pt();
1117 tmpArrayFive[3] = ptTriggGen;
1118 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1119 fHJetSpecGen->Fill(tmpArrayFive);
1122 Double_t ptUeFromCellMedianGen = rhoFromCellMedianGen*areaJetGen;
1123 Double_t ptUeConeGen = rhoConeGen*areaJetGen;
1125 //Centrality, A, pTjet-pTbgCellMedian, pTtrigg, dphi
1126 tmpArrayFive[0] = centValue;
1127 tmpArrayFive[1] = areaJetGen;
1128 tmpArrayFive[2] = jet->Pt() - ptUeFromCellMedianGen;
1129 tmpArrayFive[3] = ptTriggGen;
1130 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1131 fHJetSpecSubUeMedianGen->Fill(tmpArrayFive);
1133 //Centrality, A, pTjet-pTbgCone, pTtrigg, dphi
1134 tmpArrayFive[0] = centValue;
1135 tmpArrayFive[1] = areaJetGen;
1136 tmpArrayFive[2] = jet->Pt() - ptUeConeGen;
1137 tmpArrayFive[3] = ptTriggGen;
1138 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1139 fHJetSpecSubUeConeGen->Fill(tmpArrayFive);
1141 //Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", kT median
1142 tmpArrayFive[0] = areaJetGen;
1143 tmpArrayFive[1] = jet->Pt();
1144 tmpArrayFive[2] = jet->Pt() - ptUeFromCellMedianGen;
1145 tmpArrayFive[3] = ptUeFromCellMedianGen;
1146 tmpArrayFive[4] = rhoFromCellMedianGen;
1147 fHJetUeMedianGen->Fill(tmpArrayFive);
1149 //Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", perpendicular Cone
1150 tmpArrayFive[0] = areaJetGen;
1151 tmpArrayFive[1] = jet->Pt();
1152 tmpArrayFive[2] = jet->Pt() - ptUeConeGen;
1153 tmpArrayFive[3] = ptUeConeGen;
1154 tmpArrayFive[4] = rhoConeGen;
1155 fHJetUeConeGen->Fill(tmpArrayFive);
1158 Double_t fillRhoGen[] = { rhoConeGen, rhoFromCellMedianGen};
1159 fHRhoUeMedianVsConeGen->Fill(fillRhoGen);
1160 fillOnceGen = kFALSE;
1162 }//back to back jet-trigger pair
1167 //================ RESPONSE MATRIX ===============
1168 //Count jets and trigger-jet pairs at MC generator level
1169 for(Int_t ij=0; ij<fListJetsGen->GetEntries(); ij++){
1170 AliAODJet* jet = (AliAODJet*)(fListJetsGen->At(ij));
1172 Double_t etaJetGen = jet->Eta();
1173 Double_t ptJetGen = jet->Pt();
1175 if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){
1176 fhJetPtGen->Fill(ptJetGen); // gen level pt spectum of jets response mx normalization
1178 Double_t areaJetGen = jet->EffectiveAreaCharged();
1179 Double_t ptUeFromCellMedianGen = rhoFromCellMedianGen*areaJetGen;
1180 Double_t ptUeConeGen = rhoConeGen*areaJetGen;
1182 fhJetPtSubUeMedianGen->Fill(ptJetGen - ptUeFromCellMedianGen);
1183 fhJetPtSubUeConeGen->Fill(ptJetGen - ptUeConeGen);
1186 if(fListJets->GetEntries()>0 && fListJetsGen->GetEntries()>0){ //at least some reconstructed jets
1188 Int_t ng = (Int_t) fListJetsGen->GetEntries();
1189 Int_t nr = (Int_t) fListJets->GetEntries();
1191 //Find closest MC generator - reconstructed jets
1192 if(faGenIndex.GetSize()<nr) faGenIndex.Set(nr); //idx of gen jet assoc to rec jet
1193 if(faRecIndex.GetSize()<ng) faRecIndex.Set(ng); //idx of rec jet assoc to gen jet
1196 Printf("New Rec List %d gen index Array %d",nr,faGenIndex.GetSize());
1197 Printf("New Gen List %d rec index Array %d",ng,faRecIndex.GetSize());
1199 //matching of MC genrator level and reconstructed jets
1200 AliAnalysisHelperJetTasks::GetClosestJets(fListJetsGen,ng,fListJets,nr,faGenIndex,faRecIndex,fDebug);
1202 // Fill response matrix
1203 for(Int_t ir = 0; ir < nr; ir++){
1204 AliAODJet *recJet = (AliAODJet*) fListJets->At(ir);
1205 Double_t etaJetRec = recJet->Eta();
1206 Double_t ptJetRec = recJet->Pt();
1207 Double_t areaJetRec = recJet->EffectiveAreaCharged();
1208 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1210 if((fJetEtaMin <= etaJetRec) && (etaJetRec <= fJetEtaMax)){
1211 Int_t ig = faGenIndex[ir]; //associated generator level jet
1212 if(ig >= 0 && ig < ng){
1213 if(fDebug > 10) Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1214 AliAODJet *genJet = (AliAODJet*) fListJetsGen->At(ig);
1215 Double_t ptJetGen = genJet->Pt();
1216 Double_t etaJetGen = genJet->Eta();
1218 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1219 if((fJetEtaMin <= etaJetGen) && (etaJetGen <= fJetEtaMax)){
1220 fhJetPtGenVsJetPtRec->Fill(ptJetRec, ptJetGen);
1222 Double_t areaJetGen = genJet->EffectiveAreaCharged();
1223 Double_t ptUeFromCellMedianGen = rhoFromCellMedianGen*areaJetGen;
1224 Double_t ptUeConeGen = rhoConeGen*areaJetGen;
1225 Double_t ptUeFromCellMedianRec = rhoFromCellMedian*areaJetRec;
1226 Double_t ptUeConeRec = rhoCone*areaJetRec;
1227 fhJetPtGenVsJetPtRecSubUeMedian->Fill(ptJetRec-ptUeFromCellMedianRec,
1228 ptJetGen-ptUeFromCellMedianGen);
1229 fhJetPtGenVsJetPtRecSubUeCone->Fill(ptJetRec-ptUeConeRec, ptJetGen-ptUeConeGen);
1232 }//rec jet in eta acceptance
1233 }//loop over reconstructed jets
1234 }// # of rec jets >0
1236 //=========================== Full jet vs charged jet matrix ==========================
1238 //Count full jets and charged-jet pairs at MC generator level
1239 for(Int_t ij=0; ij<fListJetsGenFull->GetEntries(); ij++){
1240 AliAODJet* jetFull = (AliAODJet*)(fListJetsGenFull->At(ij));
1241 if(!jetFull) continue;
1242 Double_t etaJetFull = jetFull->Eta();
1243 Double_t ptJetFull = jetFull->Pt();
1245 if((fJetEtaMin<=etaJetFull) && (etaJetFull<=fJetEtaMax)){
1246 fhJetPtGenFull->Fill(ptJetFull); // generator level pt spectum of full jets
1249 if(fListJetsGen->GetEntries()>0 && fListJetsGenFull->GetEntries()>0){ //at least some reconstructed jets
1250 Int_t nful = (Int_t) fListJetsGenFull->GetEntries();
1251 Int_t nchr = (Int_t) fListJetsGen->GetEntries();
1253 //Find closest MC generator full - charged jet
1254 if(faGenIndex.GetSize()<nchr) faGenIndex.Set(nchr); //idx of gen FULL jet assoc to gen CHARGED jet
1255 if(faRecIndex.GetSize()<nful) faRecIndex.Set(nful); //idx of gen CHARGED jet assoc to gen FULL jet
1258 Printf("New Charg List %d Full index Array %d",nchr,faGenIndex.GetSize());
1259 Printf("New Full List %d Charg index Array %d",nful,faRecIndex.GetSize());
1261 //matching of MC genrator level and reconstructed jets
1262 AliAnalysisHelperJetTasks::GetClosestJets(fListJetsGenFull,nful,fListJetsGen,nchr,faGenIndex,faRecIndex,fDebug);
1264 // Fill response matrix
1265 for(Int_t ichr = 0; ichr < nchr; ichr++){ //charged jet loop
1266 AliAODJet *chJet = (AliAODJet*) fListJetsGen->At(ichr);
1267 Double_t etaJetCh = chJet->Eta();
1268 Double_t ptJetCh = chJet->Pt();
1269 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1271 if((fJetEtaMin <= etaJetCh) && (etaJetCh <= fJetEtaMax)){
1272 Int_t iful = faGenIndex[ichr]; //associated generator level jet
1273 if(iful >= 0 && iful < nful){
1274 if(fDebug > 10) Printf("%s:%d iful = %d ichr = %d",(char*)__FILE__,__LINE__,iful,ichr);
1275 AliAODJet *genJetFull = (AliAODJet*) fListJetsGenFull->At(iful);
1276 Double_t ptJetFull = genJetFull->Pt();
1277 Double_t etaJetFull = genJetFull->Eta();
1279 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1280 if((fJetEtaMin <= etaJetFull) && (etaJetFull <= fJetEtaMax)){
1281 fhJetPtGenChargVsJetPtGenFull->Fill(ptJetFull,ptJetCh);
1284 }//rec jet in eta acceptance
1285 }//loop over reconstructed jets
1286 }// # of rec jets >0
1287 }//pointer MC generator jets
1288 } //fill resp mx only for bin
1289 }//analyze generator level MC
1293 //============= RECONSTRUCTED INCLUSIVE JETS ===============
1295 Double_t etaJet = 0.0;
1296 Double_t pTJet = 0.0;
1297 Double_t areaJet = 0.0;
1298 Double_t phiJet = 0.0;
1299 //Int_t indexLeadingJet = -1;
1300 //Double_t pTLeadingJet = -10.0;
1301 //Double_t areaLeadingJet = -10.0;
1303 for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
1304 AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
1306 etaJet = jet->Eta();
1307 phiJet = jet->Phi();
1309 if(pTJet==0) continue;
1311 if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
1312 /*areaJet = jet->EffectiveAreaCharged();*/
1314 //Jet Diagnostics---------------------------------
1315 fhJetPhi->Fill(pTJet, RelativePhi(phiJet,0.0)); //phi -pi,pi
1316 fhJetEta->Fill(pTJet, etaJet);
1317 //search for leading jet
1318 /*if(pTJet > pTLeadingJet){
1319 indexLeadingJet = ij;
1320 pTLeadingJet = pTJet;
1321 areaLeadingJet = areaJet;
1324 // raw spectra of INCLUSIVE jets
1325 //Centrality, pTjet, A
1326 Double_t fillraw[] = { centValue,
1330 fHJetPtRaw->Fill(fillraw);*/
1333 if(indexLeadingJet > -1){
1334 // raw spectra of LEADING jets
1335 //Centrality, pTjet, A
1336 Double_t fillleading[] = { centValue,
1340 fHLeadingJetPtRaw->Fill(fillleading);
1344 // =============== RECONSTRUCTED TRIGGER-JET PAIRS ================
1345 if(fIsChargedMC && fFillRespMx){
1346 FillEffHistos(&particleList, &particleListGen); //Fill efficiency histos
1348 Bool_t filledOnce = kTRUE; //fill rho histogram only once per event
1349 //set ranges of the trigger loop
1350 Int_t ilow = (fHardest==0 || fHardest==1) ? indexTrigg : 0;
1351 Int_t ihigh = (fHardest==0 || fHardest==1) ? indexTrigg+1 : particleList.GetEntries();
1353 for(Int_t itrk= ilow; itrk< ihigh; itrk++){ //loop over possible trigger
1354 indexTrigg = itrk; //trigger hadron with pT >10 GeV
1356 if(indexTrigg < 0) continue;
1358 AliVParticle *triggerHadron = (AliVParticle*) particleList.At(indexTrigg);
1360 PostData(1, fOutputList);
1364 if(triggerHadron->Pt() < 6.0) continue; //all hadrons pt>6
1366 if(TMath::Abs((Float_t) triggerHadron->Eta())> fTriggerEtaCut) continue;
1368 //Fill trigger histograms
1369 fh2Ntriggers->Fill(centValue,triggerHadron->Pt()); //trigger pT
1370 fhTriggerPhi->Fill(triggerHadron->Pt(),RelativePhi(triggerHadron->Phi(),0.0)); //phi -pi,pi
1371 fhTriggerEta->Fill(triggerHadron->Pt(),triggerHadron->Eta());
1374 //---------- make trigger-jet pairs ---------
1378 for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
1379 AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
1381 etaJet = jet->Eta();
1382 phiJet = jet->Phi();
1384 if(pTJet==0) continue;
1386 if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
1387 areaJet = jet->EffectiveAreaCharged();
1388 //if(areaJet >= 0.07) injet++;
1389 //if(areaJet >= 0.4) injet4++;
1391 Double_t dphi = RelativePhi(triggerHadron->Phi(), phiJet);
1392 fhDphiTriggerJet->Fill(dphi); //Input
1394 //Dphi versus jet pT
1395 //Centrality, Dphi=phiTrig-phiJet, pTjet, pTtrigg
1396 /*Double_t filldp[] = { centValue,
1401 fHDphiVsJetPtAll->Fill(filldp);
1403 // Select back to back trigger - jet pairs
1404 if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue;
1405 fhDphiTriggerJetAccept->Fill(dphi); //Accepted
1408 //Centrality, A, pTjet, pTtrigg, dphi
1409 tmpArrayFive[0] = centValue;
1410 tmpArrayFive[1] = areaJet;
1411 tmpArrayFive[2] = pTJet;
1412 tmpArrayFive[3] = triggerHadron->Pt();
1413 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1414 fHJetSpec->Fill(tmpArrayFive);
1416 //subtract bg using estinates base on median of kT jets
1417 Double_t ptUeFromCellMedian = rhoFromCellMedian*areaJet;
1418 Double_t ptUeCone = rhoCone*areaJet;
1420 //Centrality, A, pTjet-pTbgCellMedian, pTtrigg, dphi
1421 tmpArrayFive[0] = centValue;
1422 tmpArrayFive[1] = areaJet;
1423 tmpArrayFive[2] = pTJet-ptUeFromCellMedian;
1424 tmpArrayFive[3] = triggerHadron->Pt();
1425 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1426 fHJetSpecSubUeMedian->Fill(tmpArrayFive);
1428 //Centrality, A, pTjet-pTbgCone, pTtrigg, dphi
1429 tmpArrayFive[0] = centValue;
1430 tmpArrayFive[1] = areaJet;
1431 tmpArrayFive[2] = pTJet - ptUeCone;
1432 tmpArrayFive[3] = triggerHadron->Pt();
1433 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1434 fHJetSpecSubUeCone->Fill(tmpArrayFive);
1436 //Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", kT median
1437 tmpArrayFive[0] = areaJet;
1438 tmpArrayFive[1] = pTJet;
1439 tmpArrayFive[2] = pTJet - ptUeFromCellMedian;
1440 tmpArrayFive[3] = ptUeFromCellMedian;
1441 tmpArrayFive[4] = rhoFromCellMedian;
1442 fHJetUeMedian->Fill(tmpArrayFive);
1444 //Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", Cone median
1445 tmpArrayFive[0] = areaJet;
1446 tmpArrayFive[1] = pTJet;
1447 tmpArrayFive[2] = pTJet - ptUeCone;
1448 tmpArrayFive[3] = ptUeCone;
1449 tmpArrayFive[4] = rhoCone;
1450 fHJetUeCone->Fill(tmpArrayFive);
1452 if(filledOnce){ //fill for each event only once
1453 Double_t fillRho[] = { rhoCone,rhoFromCellMedian};
1454 fHRhoUeMedianVsCone->Fill(fillRho);
1455 filledOnce = kFALSE;
1459 //Fill Jet Density In the Event A>0.07
1461 Double_t filldens[]={ (Double_t) particleList.GetEntries(),
1465 fHJetDensity->Fill(filldens);
1468 //Fill Jet Density In the Event A>0.4
1470 Double_t filldens4[]={ (Double_t) particleList.GetEntries(),
1471 injet4/fkAcceptance,
1474 fHJetDensityA4->Fill(filldens4);
1479 PostData(1, fOutputList);
1482 //----------------------------------------------------------------------------
1483 void AliAnalysisTaskJetCorePP::Terminate(const Option_t *)
1485 // Draw result to the screen
1486 // Called once at the end of the query
1488 if(fDebug) printf("AliAnalysisTaskJetCorePP DONE\n");
1489 if(!GetOutputData(1)) return;
1493 //----------------------------------------------------------------------------
1494 Int_t AliAnalysisTaskJetCorePP::GetListOfTracks(TList *list){
1495 //Fill the list of accepted tracks (passed track cut)
1496 //return consecutive index of the hardest ch hadron in the list
1498 AliAODEvent *aodevt = NULL;
1500 if(!fESD) aodevt = fAODIn;
1501 else aodevt = fAODOut;
1503 if(!aodevt) return -1;
1505 Int_t index = -1; //index of the highest particle in the list
1506 Double_t ptmax = -10;
1507 Int_t triggers[200];
1508 Int_t ntriggers = 0; //index in triggers array
1510 for(Int_t it = 0; it < aodevt->GetNumberOfTracks(); it++){
1511 AliAODTrack *tr = aodevt->GetTrack(it);
1513 if((fFilterMask > 0) && !(tr->TestFilterBit(fFilterMask))) continue;
1514 //if((fFilterMask > 0) && !(tr->IsHybridGlobalConstrainedGlobal())) continue;
1515 if(TMath::Abs((Float_t) tr->Eta()) > fTrackEtaCut) continue;
1516 if(tr->Pt() < fTrackLowPtCut) continue;
1520 if(fHardest>0){ //leading particle
1527 if(fHardest==0 && ntriggers<200){ //single inclusive
1528 if(fTriggerPtRangeLow <= tr->Pt() &&
1529 tr->Pt() < fTriggerPtRangeHigh){
1530 triggers[ntriggers] = iCount;
1538 if(fHardest==0 && ntriggers>0){ //select random inclusive trigger
1539 Int_t rnd = fRandom->Integer(ntriggers); //0 to ntriggers-1
1540 index = triggers[rnd];
1546 //----------------------------------------------------------------------------
1548 Double_t AliAnalysisTaskJetCorePP::GetBackgroundInPerpCone(Float_t jetR, Double_t jetPhi, Double_t jetEta, TList* trkList){
1549 //calculate sum of track pT in the cone perpendicular in phi to the jet
1550 //jetR = cone radius
1551 // jetPhi, jetEta = direction of the jet
1552 Int_t numberOfTrks = trkList->GetEntries();
1553 Double_t pTsum = 0.0;
1554 Double_t perpConePhi = jetPhi + TMath::Pi()/2;//perp cone w.r.t. jet in phi
1555 for(Int_t it=0; it<numberOfTrks; it++){
1556 AliVParticle *trk = (AliVParticle*) trkList->At(it);
1557 Double_t dphi = RelativePhi(perpConePhi,trk->Phi());
1558 Double_t deta = trk->Eta()-jetEta;
1560 if( (dphi*dphi + deta*deta)< (jetR*jetR)){
1568 //----------------------------------------------------------------------------
1570 Double_t AliAnalysisTaskJetCorePP::RelativePhi(Double_t mphi,Double_t vphi){
1571 //Get relative azimuthal angle of two particles -pi to pi
1572 if (vphi < -TMath::Pi()) vphi += TMath::TwoPi();
1573 else if (vphi > TMath::Pi()) vphi -= TMath::TwoPi();
1575 if (mphi < -TMath::Pi()) mphi += TMath::TwoPi();
1576 else if (mphi > TMath::Pi()) mphi -= TMath::TwoPi();
1578 Double_t dphi = mphi - vphi;
1579 if (dphi < -TMath::Pi()) dphi += TMath::TwoPi();
1580 else if (dphi > TMath::Pi()) dphi -= TMath::TwoPi();
1582 return dphi;//dphi in [-Pi, Pi]
1586 //----------------------------------------------------------------------------
1587 Bool_t AliAnalysisTaskJetCorePP::SelectMCGenTracks(AliVParticle *trk, TList *trkList, Double_t &ptLeading, Int_t &index, Int_t counter){
1588 //fill track efficiency denominator
1589 if(trk->Pt() < fTrackLowPtCut) return kFALSE;
1590 if(trk->Charge()==0) return kFALSE;
1591 if(TMath::Abs((Float_t) trk->Eta()) > fTrackEtaCut) return kFALSE;
1593 if(fFillRespMx) fhPtTrkTruePrimGen->Fill(trk->Pt(),trk->Eta()); //Efficiency denominator
1596 if(fHardest>0){ //leading particle
1597 if(ptLeading < trk->Pt()){
1599 ptLeading = trk->Pt();
1603 if(fHardest==0){ //single inclusive
1605 if(fTriggerPtRangeLow <= trk->Pt() &&
1606 trk->Pt() < fTriggerPtRangeHigh){
1614 //----------------------------------------------------------------------------
1615 void AliAnalysisTaskJetCorePP::FillEffHistos(TList *recList, TList *genList){
1616 //fill track efficiency numerator
1617 if(!recList) return;
1618 if(!genList) return;
1619 Int_t nRec = recList->GetEntries();
1620 Int_t nGen = genList->GetEntries();
1621 for(Int_t ir=0; ir<nRec; ir++){
1622 AliAODTrack *trkRec = (AliAODTrack*) recList->At(ir);
1623 if(!trkRec) continue;
1624 Int_t recLabel = TMath::Abs(trkRec->GetLabel());
1625 Bool_t matched = kFALSE;
1627 for(Int_t ig=0; ig<nGen; ig++){
1628 AliVParticle *trkGen = (AliVParticle*) genList->At(ig);
1629 if(!trkGen) continue;
1630 Int_t genLabel = TMath::Abs(trkGen->GetLabel());
1631 if(recLabel==genLabel){
1632 fhPtTrkTruePrimRec->Fill(trkGen->Pt(),trkGen->Eta());
1637 if(!matched) fhPtTrkSecOrFakeRec->Fill(trkRec->Pt(),trkRec->Eta());
1642 //________________________________________________________________
1643 void AliAnalysisTaskJetCorePP::EstimateBgRhoMedian(TList *listJet, TList* listPart, Double_t &rhoMedian, Int_t mode){
1644 //Estimate background rho by means of integrating track pT outside identified jet cones
1648 //identify leading jet in the full track acceptance
1649 Int_t idxLeadingJet = -1;
1650 Double_t pTleading = 0.0;
1651 AliAODJet* jet = NULL;
1653 for(Int_t ij = 0; ij < listJet->GetEntries(); ij++){ //loop over bg kT jets
1654 jet = (AliAODJet*)(listJet->At(ij));
1656 if(TMath::Abs(jet->Eta()) > fTrackEtaCut) continue; //track acceptance cut
1657 if(jet->Pt() > pTleading){
1659 pTleading = jet->Pt();
1664 static Double_t jphi[100];
1665 static Double_t jeta[100];
1666 static Double_t jRsquared[100];
1668 for(Int_t ij = 0; ij< listJet->GetEntries(); ij++){ //loop over bg kT jets
1670 jet = (AliAODJet*)(listJet->At(ij));
1672 if(TMath::Abs(jet->Eta()) > fTrackEtaCut) continue;
1674 if((jet->Pt() > fBgMaxJetPt || ij== idxLeadingJet) && njets<100){
1675 jphi[njets] = RelativePhi(jet->Phi(),0.0); //-pi,pi
1676 jeta[njets] = jet->Eta();
1677 jRsquared[njets] = fSafetyMargin*jet->EffectiveAreaCharged()/TMath::Pi(); //scale jet radius
1683 static Double_t nOutCone[10][4];
1684 static Double_t sumPtOutOfCone[10][4];
1686 for(Int_t ie=0; ie<fnEta; ie++){
1687 for(Int_t ip=0; ip<fnPhi; ip++){
1688 nOutCone[ip][ie] = 0.0; //initialize counter
1689 sumPtOutOfCone[ip][ie] = 0.0;
1693 Double_t rndphi, rndeta;
1694 Double_t rndphishift, rndetashift;
1695 Double_t dphi, deta;
1698 for(Int_t it=0; it<fnTrials; it++){
1700 rndphi = fRandom->Uniform(0, fPhiSize);
1701 rndeta = fRandom->Uniform(0, fEtaSize);
1703 for(Int_t ip=0; ip<fnPhi; ip++){ //move radom position to each cell
1704 rndphishift = rndphi + ip*fPhiSize - TMath::Pi();
1705 for(Int_t ie=0; ie<fnEta; ie++){
1706 rndetashift = rndeta + ie*fEtaSize - fEtaSize;
1708 bIsInCone = 0; //tag if trial is in the jet cone
1709 for(Int_t ij=0; ij<njets; ij++){
1710 deta = jeta[ij] - rndetashift;
1711 dphi = RelativePhi(rndphishift,jphi[ij]);
1712 if((dphi*dphi + deta*deta) < jRsquared[ij]){
1717 if(!bIsInCone) nOutCone[ip][ie]++;
1723 //Sum particle pT outside jets in cells
1724 Int_t npart = listPart->GetEntries();
1725 Int_t phicell,etacell;
1726 AliVParticle *part = NULL;
1727 for(Int_t ip=0; ip < npart; ip++){
1729 part = (AliVParticle*) listPart->At(ip);
1731 if( TMath::Abs(part->Eta()) > fEtaSize) continue; //skip part outside +-0.7 in eta
1733 bIsInCone = 0; //init
1734 for(Int_t ij=0; ij<njets; ij++){
1735 dphi = RelativePhi(jphi[ij], part->Phi());
1736 deta = jeta[ij] - part->Eta();
1737 if((dphi*dphi + deta*deta) < jRsquared[ij]){
1743 phicell = TMath::Nint(TMath::Floor((RelativePhi(part->Phi(),0.0) + TMath::Pi())/fPhiSize));
1744 etacell = TMath::Nint(TMath::Floor((part->Eta()+fEtaSize)/fEtaSize));
1745 sumPtOutOfCone[phicell][etacell]+= part->Pt();
1749 static Double_t rhoInCells[20];
1750 Double_t relativeArea;
1752 Double_t bufferArea=0.0, bufferPt=0.0; //sum cells where A< fJetFreeAreaFrac
1753 for(Int_t ip=0; ip<fnPhi; ip++){
1754 for(Int_t ie=0; ie<fnEta; ie++){
1755 relativeArea = nOutCone[ip][ie]/fnTrials;
1756 //cout<<"RA= "<< relativeArea<<" BA="<<bufferArea<<" BPT="<< bufferPt <<endl;
1758 bufferArea += relativeArea;
1759 bufferPt += sumPtOutOfCone[ip][ie];
1760 if(bufferArea > fJetFreeAreaFrac){
1761 rhoInCells[nCells] = bufferPt/(bufferArea*fCellArea);
1763 if(mode==1) fhCellAreaToMedianGen->Fill(bufferArea*fCellArea);
1764 else fhCellAreaToMedian->Fill(bufferArea*fCellArea);
1775 rhoMedian = TMath::Median(nCells, rhoInCells);
1776 if(mode==1){ //mc data
1777 fhEntriesToMedianGen->Fill(nCells);
1779 fhEntriesToMedian->Fill(nCells);
1785 //__________________________________________________________________
1786 void AliAnalysisTaskJetCorePP::EstimateBgCone(TList *listJet, TList* listPart, Double_t &rhoPerpCone){
1788 rhoPerpCone = 0.0; //init
1790 if(!listJet) return;
1791 Int_t njet = listJet->GetEntries();
1792 Int_t npart = listPart->GetEntries();
1793 Double_t pTleading = 0.0;
1794 Double_t phiLeading = 1000.;
1795 Double_t etaLeading = 1000.;
1797 for(Int_t jsig=0; jsig < njet; jsig++){
1798 AliAODJet* signaljet = (AliAODJet*)(listJet->At(jsig));
1799 if(!signaljet) continue;
1800 if((signaljet->Eta()<fJetEtaMin) && (fJetEtaMax<signaljet->Eta())) continue; //acceptance cut
1801 if(signaljet->Pt() >= pTleading){ //replace leading and subleading jet
1802 pTleading = signaljet->Pt();
1803 phiLeading = signaljet->Phi();
1804 etaLeading = signaljet->Eta();
1808 Double_t phileftcone = phiLeading + TMath::Pi()/2;
1809 Double_t phirightcone = phiLeading - TMath::Pi()/2;
1812 for(Int_t ip=0; ip < npart; ip++){
1814 AliVParticle *part = (AliVParticle*) listPart->At(ip);
1819 dp = RelativePhi(phileftcone, part->Phi());
1820 de = etaLeading - part->Eta();
1821 if( sqrt(dp*dp + de*de)< fBgConeR) rhoPerpCone += part->Pt();
1823 dp = RelativePhi(phirightcone, part->Phi());
1824 if( sqrt(dp*dp + de*de)< fBgConeR) rhoPerpCone += part->Pt();
1828 //normalize total pT by two times cone are
1829 rhoPerpCone = rhoPerpCone/(2*TMath::Pi()*fBgConeR*fBgConeR);
1832 //__________________________________________________________________
1834 void AliAnalysisTaskJetCorePP::ReadTClonesArray(TString bname, TList *list){
1835 //Convert TClones array of jets to TList
1837 if(!strlen(bname.Data())){
1838 AliError(Form("Jet branch %s not set.", bname.Data()));
1842 TClonesArray *array=0x0;
1844 if(fAODOut&&!array){
1845 array = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(bname.Data()));
1847 if(fAODExtension&&!array){
1848 array = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(bname.Data()));
1851 array = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(bname.Data()));
1854 list->Clear(); //clear list beforehand
1858 Printf("########## %s: %d jets",bname.Data(), array->GetEntriesFast());
1860 for(Int_t iJet = 0; iJet < array->GetEntriesFast(); iJet++) {
1861 AliAODJet *jet = dynamic_cast<AliAODJet*>((*array)[iJet]);
1862 if (jet) list->Add(jet);