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 **************************************************************************/
28 #include "TClonesArray.h"
36 #include "THnSparse.h"
46 #include "AliAnalysisTask.h"
47 #include "AliAnalysisManager.h"
49 #include "AliVEvent.h"
50 #include "AliESDEvent.h"
51 #include "AliMCEvent.h"
52 #include "AliESDInputHandler.h"
53 #include "AliCentrality.h"
54 #include "AliAnalysisHelperJetTasks.h"
55 #include "AliInputEventHandler.h"
56 #include "AliAODJetEventBackground.h"
57 #include "AliGenPythiaEventHeader.h"
58 #include "AliAODMCParticle.h"
59 #include "AliMCParticle.h"
60 #include "AliAODEvent.h"
61 #include "AliAODHandler.h"
62 #include "AliAODJet.h"
63 #include "AliVVertex.h"
64 #include "AliAnalysisTaskJetCorePP.h"
65 #include "AliHeader.h" //KF//
70 ClassImp(AliAnalysisTaskJetCorePP)
72 //Filip Krizek 1st March 2013
74 //---------------------------------------------------------------------
75 AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP() :
84 fJetBranchNameChargMC(""),
85 fJetBranchNameKine(""),
86 fJetBranchNameFullMC(""),
88 fJetBranchNameBgChargMC(""),
89 fJetBranchNameBgKine(""),
92 fListJetsGenFull(0x0),
96 fSystem(0), //pp=0 pPb=1
101 fOfflineTrgMask(AliVEvent::kAny),
112 fTrackLowPtCut(0.15),
113 fUseExchContainer(0),
115 fHistEvtSelection(0x0),
118 fHJetSpecSubUeMedian(0x0),
119 fHJetSpecSubUeCone(0x0),
121 fHJetPhiCorrSubUeMedian(0x0),
122 fHJetPhiCorrSubUeCone(0x0),
125 fHRhoUeMedianVsCone(0x0),
127 //fHJetDensityA4(0x0),
133 fhVertexZAccept(0x0),
135 fhContribVtxAccept(0x0),
136 fhDphiTriggerJet(0x0),
137 fhDphiTriggerJetAccept(0x0),
139 fhCentralityAccept(0x0),
141 //fHLeadingJetPtRaw(0x0),
142 //fHDphiVsJetPtAll(0x0),
143 fhJetPtGenVsJetPtRec(0x0),
144 fhJetPtGenVsJetPtRecSubUeMedian(0x0),
145 fhJetPtGenVsJetPtRecSubUeCone(0x0),
147 fhJetPtSubUeMedianGen(0x0),
148 fhJetPtSubUeConeGen(0x0),
149 fhJetPtGenChargVsJetPtGenFull(0x0),
151 fh2NtriggersGen(0x0),
153 fHJetSpecSubUeMedianGen(0x0),
154 fHJetSpecSubUeConeGen(0x0),
155 fHJetPhiCorrGen(0x0),
156 fHJetPhiCorrSubUeMedianGen(0x0),
157 fHJetPhiCorrSubUeConeGen(0x0),
158 fHJetUeMedianGen(0x0),
160 fhPtTrkTruePrimRec(0x0),
161 fhPtTrkTruePrimGen(0x0),
162 fhPtTrkSecOrFakeRec(0x0),
163 fHRhoUeMedianVsConeGen(0x0),
164 fhEntriesToMedian(0x0),
165 fhEntriesToMedianGen(0x0),
166 fhCellAreaToMedian(0x0),
167 fhCellAreaToMedianGen(0x0),
173 fkAcceptance(2.0*TMath::Pi()*1.8),
174 fkDeltaPhiCut(TMath::Pi()-0.8),
180 fh1PtHardTrials(0x0),
183 fEventNumberRangeLow(0),
184 fEventNumberRangeHigh(99),
185 fTriggerPtRangeLow(0.0),
186 fTriggerPtRangeHigh(10000.0),
190 fJetFreeAreaFrac(0.5),
194 fPhiSize(2*TMath::Pi()/fnPhi),
195 fCellArea(fPhiSize*fEtaSize),
197 fDoubleBinning(kFALSE)
199 // default Constructor
202 //---------------------------------------------------------------------
204 AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP(const char *name) :
205 AliAnalysisTaskSE(name),
213 fJetBranchNameChargMC(""),
214 fJetBranchNameKine(""),
215 fJetBranchNameFullMC(""),
216 fJetBranchNameBg(""),
217 fJetBranchNameBgChargMC(""),
218 fJetBranchNameBgKine(""),
221 fListJetsGenFull(0x0),
225 fSystem(0), //pp=0 pPb=1
230 fOfflineTrgMask(AliVEvent::kAny),
241 fTrackLowPtCut(0.15),
242 fUseExchContainer(0),
244 fHistEvtSelection(0x0),
247 fHJetSpecSubUeMedian(0x0),
248 fHJetSpecSubUeCone(0x0),
250 fHJetPhiCorrSubUeMedian(0x0),
251 fHJetPhiCorrSubUeCone(0x0),
254 fHRhoUeMedianVsCone(0x0),
256 //fHJetDensityA4(0x0),
262 fhVertexZAccept(0x0),
264 fhContribVtxAccept(0x0),
265 fhDphiTriggerJet(0x0),
266 fhDphiTriggerJetAccept(0x0),
268 fhCentralityAccept(0x0),
270 //fHLeadingJetPtRaw(0x0),
271 //fHDphiVsJetPtAll(0x0),
272 fhJetPtGenVsJetPtRec(0x0),
273 fhJetPtGenVsJetPtRecSubUeMedian(0x0),
274 fhJetPtGenVsJetPtRecSubUeCone(0x0),
276 fhJetPtSubUeMedianGen(0x0),
277 fhJetPtSubUeConeGen(0x0),
278 fhJetPtGenChargVsJetPtGenFull(0x0),
280 fh2NtriggersGen(0x0),
282 fHJetSpecSubUeMedianGen(0x0),
283 fHJetSpecSubUeConeGen(0x0),
284 fHJetPhiCorrGen(0x0),
285 fHJetPhiCorrSubUeMedianGen(0x0),
286 fHJetPhiCorrSubUeConeGen(0x0),
287 fHJetUeMedianGen(0x0),
289 fhPtTrkTruePrimRec(0x0),
290 fhPtTrkTruePrimGen(0x0),
291 fhPtTrkSecOrFakeRec(0x0),
292 fHRhoUeMedianVsConeGen(0x0),
293 fhEntriesToMedian(0x0),
294 fhEntriesToMedianGen(0x0),
295 fhCellAreaToMedian(0x0),
296 fhCellAreaToMedianGen(0x0),
302 fkAcceptance(2.0*TMath::Pi()*1.8),
303 fkDeltaPhiCut(TMath::Pi()-0.8),
309 fh1PtHardTrials(0x0),
312 fEventNumberRangeLow(0),
313 fEventNumberRangeHigh(99),
314 fTriggerPtRangeLow(0.0),
315 fTriggerPtRangeHigh(10000.0),
319 fJetFreeAreaFrac(0.5),
323 fPhiSize(2*TMath::Pi()/fnPhi),
324 fCellArea(fPhiSize*fEtaSize),
326 fDoubleBinning(kFALSE)
329 DefineOutput(1, TList::Class());
332 if(dummy.Contains("KINE")){
333 DefineInput(1, TClonesArray::Class());
334 DefineInput(2, TClonesArray::Class());
338 //--------------------------------------------------------------
339 AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP(const AliAnalysisTaskJetCorePP& a):
340 AliAnalysisTaskSE(a.GetName()),
344 fAODExtension(a.fAODExtension),
345 fMcEvent(a.fMcEvent),
346 fMcHandler(a.fMcHandler),
347 fJetBranchName(a.fJetBranchName),
348 fJetBranchNameChargMC(a.fJetBranchNameChargMC),
349 fJetBranchNameKine(a.fJetBranchNameKine),
350 fJetBranchNameFullMC(a.fJetBranchNameFullMC),
351 fJetBranchNameBg(a.fJetBranchNameBg),
352 fJetBranchNameBgChargMC(a.fJetBranchNameBgChargMC),
353 fJetBranchNameBgKine(a.fJetBranchNameBgKine),
354 fListJets(a.fListJets),
355 fListJetsGen(a.fListJetsGen),
356 fListJetsGenFull(a.fListJetsGenFull),
357 fListJetsBg(a.fListJetsBg),
358 fListJetsBgGen(a.fListJetsBgGen),
359 fNonStdFile(a.fNonStdFile),
361 fJetParamR(a.fJetParamR),
362 fBgJetParamR(a.fBgJetParamR),
363 fBgMaxJetPt(a.fBgMaxJetPt),
364 fBgConeR(a.fBgConeR),
365 fOfflineTrgMask(a.fOfflineTrgMask),
366 fMinContribVtx(a.fMinContribVtx),
367 fVtxZMin(a.fVtxZMin),
368 fVtxZMax(a.fVtxZMax),
369 fFilterMask(a.fFilterMask),
370 fCentMin(a.fCentMin),
371 fCentMax(a.fCentMax),
372 fJetEtaMin(a.fJetEtaMin),
373 fJetEtaMax(a.fJetEtaMax),
374 fTriggerEtaCut(a.fTriggerEtaCut),
375 fTrackEtaCut(a.fTrackEtaCut),
376 fTrackLowPtCut(a.fTrackLowPtCut),
377 fUseExchContainer(a.fUseExchContainer),
378 fOutputList(a.fOutputList),
379 fHistEvtSelection(a.fHistEvtSelection),
380 fh2Ntriggers(a.fh2Ntriggers),
381 fHJetSpec(a.fHJetSpec),
382 fHJetSpecSubUeMedian(a.fHJetSpecSubUeMedian),
383 fHJetSpecSubUeCone(a.fHJetSpecSubUeCone),
384 fHJetPhiCorr(a.fHJetPhiCorr),
385 fHJetPhiCorrSubUeMedian(a.fHJetPhiCorrSubUeMedian),
386 fHJetPhiCorrSubUeCone(a.fHJetPhiCorrSubUeCone),
387 fHJetUeMedian(a.fHJetUeMedian),
388 fHJetUeCone(a.fHJetUeCone),
389 fHRhoUeMedianVsCone(a.fHRhoUeMedianVsCone),
390 //fHJetDensity(a.fHJetDensity),
391 //fHJetDensityA4(a.fHJetDensityA4),
392 fhJetPhi(a.fhJetPhi),
393 fhTriggerPhi(a.fhTriggerPhi),
394 fhJetEta(a.fhJetEta),
395 fhTriggerEta(a.fhTriggerEta),
396 fhVertexZ(a.fhVertexZ),
397 fhVertexZAccept(a.fhVertexZAccept),
398 fhContribVtx(a.fhContribVtx),
399 fhContribVtxAccept(a.fhContribVtxAccept),
400 fhDphiTriggerJet(a.fhDphiTriggerJet),
401 fhDphiTriggerJetAccept(a.fhDphiTriggerJetAccept),
402 fhCentrality(a.fhCentrality),
403 fhCentralityAccept(a.fhCentralityAccept),
404 //fHJetPtRaw(a.fHJetPtRaw),
405 //fHLeadingJetPtRaw(a.fHLeadingJetPtRaw),
406 //fHDphiVsJetPtAll(a.fHDphiVsJetPtAll),
407 fhJetPtGenVsJetPtRec(a.fhJetPtGenVsJetPtRec),
408 fhJetPtGenVsJetPtRecSubUeMedian(a.fhJetPtGenVsJetPtRecSubUeMedian),
409 fhJetPtGenVsJetPtRecSubUeCone(a.fhJetPtGenVsJetPtRecSubUeCone),
410 fhJetPtGen(a.fhJetPtGen),
411 fhJetPtSubUeMedianGen(a.fhJetPtSubUeMedianGen),
412 fhJetPtSubUeConeGen(a.fhJetPtSubUeConeGen),
413 fhJetPtGenChargVsJetPtGenFull(a.fhJetPtGenChargVsJetPtGenFull),
414 fhJetPtGenFull(a.fhJetPtGenFull),
415 fh2NtriggersGen(a.fh2NtriggersGen),
416 fHJetSpecGen(a.fHJetSpecGen),
417 fHJetSpecSubUeMedianGen(a.fHJetSpecSubUeMedianGen),
418 fHJetSpecSubUeConeGen(a.fHJetSpecSubUeConeGen),
419 fHJetPhiCorrGen(a.fHJetPhiCorrGen),
420 fHJetPhiCorrSubUeMedianGen(a.fHJetPhiCorrSubUeMedianGen),
421 fHJetPhiCorrSubUeConeGen(a.fHJetPhiCorrSubUeConeGen),
422 fHJetUeMedianGen(a.fHJetUeMedianGen),
423 fHJetUeConeGen(a.fHJetUeConeGen),
424 fhPtTrkTruePrimRec(a.fhPtTrkTruePrimRec),
425 fhPtTrkTruePrimGen(a.fhPtTrkTruePrimGen),
426 fhPtTrkSecOrFakeRec(a.fhPtTrkSecOrFakeRec),
427 fHRhoUeMedianVsConeGen(a.fHRhoUeMedianVsConeGen),
428 fhEntriesToMedian(a.fhEntriesToMedian),
429 fhEntriesToMedianGen(a.fhEntriesToMedianGen),
430 fhCellAreaToMedian(a.fhCellAreaToMedian),
431 fhCellAreaToMedianGen(a.fhCellAreaToMedianGen),
432 fIsChargedMC(a.fIsChargedMC),
434 fIsFullMC(a.fIsFullMC),
435 faGenIndex(a.faGenIndex),
436 faRecIndex(a.faRecIndex),
437 fkAcceptance(a.fkAcceptance),
438 fkDeltaPhiCut(a.fkDeltaPhiCut),
440 fh1Trials(a.fh1Trials),
441 fh1AvgTrials(a.fh1AvgTrials),
442 fh1PtHard(a.fh1PtHard),
443 fh1PtHardNoW(a.fh1PtHardNoW),
444 fh1PtHardTrials(a.fh1PtHardTrials),
445 fAvgTrials(a.fAvgTrials),
446 fHardest(a.fHardest),
447 fEventNumberRangeLow(a.fEventNumberRangeLow),
448 fEventNumberRangeHigh(a.fEventNumberRangeHigh),
449 fTriggerPtRangeLow(a.fTriggerPtRangeLow),
450 fTriggerPtRangeHigh(a.fTriggerPtRangeHigh),
451 fFillRespMx(a.fFillRespMx),
453 fnTrials(a.fnTrials),
454 fJetFreeAreaFrac(a.fJetFreeAreaFrac),
457 fEtaSize(a.fEtaSize),
458 fPhiSize(a.fPhiSize),
459 fCellArea(a.fCellArea),
460 fSafetyMargin(a.fSafetyMargin),
461 fDoubleBinning(a.fDoubleBinning)
466 //--------------------------------------------------------------
468 AliAnalysisTaskJetCorePP& AliAnalysisTaskJetCorePP::operator = (const AliAnalysisTaskJetCorePP& a){
469 // assignment operator
470 this->~AliAnalysisTaskJetCorePP();
471 new(this) AliAnalysisTaskJetCorePP(a);
474 //--------------------------------------------------------------
476 AliAnalysisTaskJetCorePP::~AliAnalysisTaskJetCorePP()
481 delete fListJetsGenFull;
483 delete fListJetsBgGen;
484 delete fOutputList; // ????
488 //--------------------------------------------------------------
491 Bool_t AliAnalysisTaskJetCorePP::Notify()
493 //Implemented Notify() to read the cross sections
494 //and number of trials from pyxsec.root
495 //inspired by AliAnalysisTaskJetSpectrum2::Notify()
496 if(!(fIsChargedMC || fIsKine)) return kFALSE;
497 Float_t xsection = 0;
502 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
504 if(fDebug>1) AliError("ESD not available");
505 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
508 fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
511 if(fNonStdFile.Length()!=0){
512 // case that we have an AOD extension we can fetch the jets from the extended output
513 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
514 fAODExtension = aodH ? aodH->GetExtension(fNonStdFile.Data()) : 0;
516 if(fDebug>1) Printf("AODExtension found for %s",fNonStdFile.Data());
520 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
523 TFile *curfile = tree->GetCurrentFile();
525 Error("Notify","No current file");
528 if(!fh1Xsec || !fh1Trials){
529 Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
532 AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,trials);
533 fh1Xsec->Fill("<#sigma>",xsection);
534 // construct a poor man average trials
535 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
536 if(trials>=nEntries && nEntries>0.) fAvgTrials = trials/nEntries;
537 fh1Trials->Fill("#sum{ntrials}",trials);
540 if(fDebug)Printf("Reading File %s",fInputHandler->GetTree()->GetCurrentFile()->GetName());
546 //--------------------------------------------------------------
548 void AliAnalysisTaskJetCorePP::Init()
550 // check for jet branches
551 if(fJetBranchNameKine.Length()==0){
552 if(!strlen(fJetBranchName.Data())){
553 AliError("Jet branch name not set.");
556 fMcHandler = dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
561 //--------------------------------------------------------------
563 void AliAnalysisTaskJetCorePP::UserCreateOutputObjects()
565 // Create histograms and initilize variables
566 fSafetyMargin = fBgConeR*fBgConeR /(fBgJetParamR*fBgJetParamR);
569 fListJets = new TList(); //reconstructed level antikt jets
570 fListJetsBg = new TList(); //reconstructed jets to be removed from bg
572 fIsChargedMC = (fJetBranchNameChargMC.Length()>0) ? kTRUE : kFALSE;
573 fIsKine = (fJetBranchNameKine.Length()>0) ? kTRUE : kFALSE;
574 fIsFullMC = (fJetBranchNameFullMC.Length()>0) ? kTRUE : kFALSE;
576 fRandom = new TRandom3(0);
578 if(fIsChargedMC || fIsKine){ //full MC or pure kine
579 fListJetsGen = new TList(); //generator level charged antikt jets
580 fListJetsBgGen = new TList(); //generator level jets to be removed from bg
583 fListJetsGenFull = new TList(); //generator level full jets
586 if(!fOutputList) fOutputList = new TList();
587 fOutputList->SetOwner(kTRUE);
589 Bool_t oldStatus = TH1::AddDirectoryStatus();
590 TH1::AddDirectory(kFALSE);
592 fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
593 fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
594 fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
595 fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
596 fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
597 fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
598 fHistEvtSelection->GetXaxis()->SetBinLabel(6,"event number (rejected)");
600 fOutputList->Add(fHistEvtSelection);
602 Int_t nBinsCentrality = (fSystem==0) ? 1 : 10; // pp=1 else 10
603 //trigger pt spectrum (reconstructed)
604 fh2Ntriggers = new TH2F("fh2Ntriggers","# of triggers",
605 nBinsCentrality,0.0,100.0,50,0.0,50.0);
606 if(!fIsKine) fOutputList->Add(fh2Ntriggers);
608 Int_t bw = (fDoubleBinning==0) ? 1 : 2; //make larger bin width
610 //Centrality, A, pTjet, pTtrigg, dphi
611 const Int_t dimSpec = 5;
612 const Int_t nBinsSpec[dimSpec] = {nBinsCentrality, 50, bw*110, 50, TMath::Nint(10*(TMath::Pi()-fkDeltaPhiCut))};
613 const Double_t lowBinSpec[dimSpec] = {0.0, 0.0, -20, 0.0, fkDeltaPhiCut};
614 const Double_t hiBinSpec[dimSpec] = {100.0, 1.0, 200.0,50.0, TMath::Pi()};
615 fHJetSpec = new THnSparseF("fHJetSpec",
616 "Recoil jet spectrum [cent,A,pTjet,pTtrig,dphi]",
617 dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
618 if(!fIsKine) fOutputList->Add(fHJetSpec);
620 //background estimated as median of kT jets
621 fHJetSpecSubUeMedian = (THnSparseF*) fHJetSpec->Clone("fHJetSpecSubUeMedian");
622 fHJetSpecSubUeMedian->SetTitle("Recoil jet spectrum [cent,A,pTjet-pTUe,pTtrig,dphi]");
623 if(!fIsKine) fOutputList->Add(fHJetSpecSubUeMedian);
624 //background estimated as weighted median of kT jets ala Cone
625 fHJetSpecSubUeCone = (THnSparseF*) fHJetSpec->Clone("fHJetSpecSubUeCone");
626 fHJetSpecSubUeCone->SetTitle("Recoil jet spectrum [cent,A,pTjet-pTUe,pTtrig,dphi]");
627 if(!fIsKine) fOutputList->Add(fHJetSpecSubUeCone);
630 //A, pTjet, pTtrigg, dphi
631 const Int_t dimCor = 4;
632 const Int_t nBinsCor[dimCor] = {50, 110, 50, 100};
633 const Double_t lowBinCor[dimCor] = {0.0, -20, 0.0, -0.5*TMath::Pi()};
634 const Double_t hiBinCor[dimCor] = {1.0, 200, 50.0, 1.5*TMath::Pi()};
635 fHJetPhiCorr = new THnSparseF("fHJetPhiCorr",
636 "Recoil jet spectrum [A,pTjet,pTtrig,dphi]",
637 dimCor,nBinsCor,lowBinCor,hiBinCor);
639 if(!fIsKine) fOutputList->Add(fHJetPhiCorr); // Dphi distribution jet-triger
641 //background estimated as median of kT jets
642 fHJetPhiCorrSubUeMedian = (THnSparseF*) fHJetPhiCorr->Clone("fHJetPhiCorrSubUeMedian");
643 fHJetPhiCorrSubUeMedian->SetTitle("Recoil jet spectrum [A,pTjet-pTUe,pTtrig,dphi]");
644 if(!fIsKine) fOutputList->Add(fHJetPhiCorrSubUeMedian);
645 //background estimated as weighted median of kT jets ala Cone
646 fHJetPhiCorrSubUeCone = (THnSparseF*) fHJetPhiCorr->Clone("fHJetPhiCorrSubUeCone");
647 fHJetPhiCorrSubUeCone->SetTitle("Recoil jet spectrum [A,pTjet-pTUe,pTtrig,dphi]");
648 if(!fIsKine) fOutputList->Add(fHJetPhiCorrSubUeCone);
651 //------------------- HISTOS FOR DIAGNOSTIC ----------------------
652 //A, pTjet, pTjet-pTUe, pTUe, rhoUe bg estimate from kT median
653 const Int_t dimSpecMed = 5;
654 const Int_t nBinsSpecMed[dimSpecMed] = {25, 50, 60, 20, 20};
655 const Double_t lowBinSpecMed[dimSpecMed] = {0.0, 0.0, -20.0, 0.0, 0.0};
656 const Double_t hiBinSpecMed[dimSpecMed] = {1.0, 100.0, 100.0, 10.0, 20.0};
657 fHJetUeMedian = new THnSparseF("fHJetUeMedian",
658 "Recoil jet spectrum [A,pTjet,pTjet-pTUe, pTUe, rhoUe]",
659 dimSpecMed, nBinsSpecMed, lowBinSpecMed, hiBinSpecMed);
660 if(!fIsKine) fOutputList->Add(fHJetUeMedian);
662 //A, pTjet, pTjet-pTUe, pTUe, rhoUe bg estimate from kT median Cone
663 fHJetUeCone = (THnSparseF*) fHJetUeMedian->Clone("fHJetUeCone");
664 fHJetUeCone->SetTitle("Recoil jet spectrum [A,pTjet,pTjet-pTUe, pTUe, rhoUe]");
665 if(!fIsKine) fOutputList->Add(fHJetUeCone);
667 //rho bacground reconstructed data
668 const Int_t dimRho = 2;
669 const Int_t nBinsRho[dimRho] = {50 , 50};
670 const Double_t lowBinRho[dimRho] = {0.0 , 0.0};
671 const Double_t hiBinRho[dimRho] = {20.0 , 20.0};
673 fHRhoUeMedianVsCone = new THnSparseF("hRhoUeMedianVsCone","[Rho Cone, Rho Median]",
674 dimRho, nBinsRho, lowBinRho, hiBinRho);
675 if(!fIsKine) fOutputList->Add(fHRhoUeMedianVsCone);
677 //Jet number density histos [Trk Mult, jet density, pT trigger]
678 /*const Int_t dimJetDens = 3;
679 const Int_t nBinsJetDens[dimJetDens] = {100, 100, 10};
680 const Double_t lowBinJetDens[dimJetDens] = {0.0, 0.0, 0.0};
681 const Double_t hiBinJetDens[dimJetDens] = {500.0, 5.0, 50.0};
683 fHJetDensity = new THnSparseF("fHJetDensity","Jet dens vs trk mult A>0.07",
684 dimJetDens,nBinsJetDens,lowBinJetDens,hiBinJetDens);
686 fHJetDensityA4 =new THnSparseF("fHJetDensityA4","Jet dens vs trk mult A>0.4",
687 dimJetDens,nBinsJetDens,lowBinJetDens,hiBinJetDens);
689 fOutputList->Add(fHJetDensity);
690 fOutputList->Add(fHJetDensityA4);
693 //inclusive azimuthal and pseudorapidity histograms
694 fhJetPhi = new TH2D("fhJetPhi","Azim dist jets vs pTjet",
695 50, 0, 100, 50,-TMath::Pi(),TMath::Pi());
696 fhTriggerPhi= new TH2D("fhTriggerPhi","azim dist trig had vs pTtrigg",
697 50, 0, 50, 50,-TMath::Pi(),TMath::Pi());
698 fhJetEta = new TH2D("fhJetEta","Eta dist jets vs pTjet",
699 50,0, 100, 40,-0.9,0.9);
700 fhTriggerEta = new TH2D("fhTriggerEta","Eta dist trig had vs pTtrigg",
701 50, 0, 50, 40,-0.9,0.9);
703 fhVertexZ = new TH1D("fhVertexZ","z vertex",40,-20,20);
704 fhVertexZAccept = new TH1D("fhVertexZAccept","z vertex after cut",40,-20,20);
705 fhContribVtx = new TH1D("fhContribVtx","contrib to vtx",200,0,200);
706 fhContribVtxAccept = new TH1D("fhContribVtxAccept","contrib to vtx after cut",200,0,200);
707 fhDphiTriggerJet = new TH1D("fhDphiTriggerJet","Deltaphi trig-jet",50, -TMath::Pi(),TMath::Pi());
708 fhDphiTriggerJetAccept = new TH1D("fhDphiTriggerJetAccept","Deltaphi trig-jet after cut",50, -TMath::Pi(),TMath::Pi());
709 fhCentrality = new TH1D("fhCentrality","Centrality",20,0,100);
710 fhCentralityAccept = new TH1D("fhCentralityAccept","Centrality after cut",20,0,100);
711 fhEntriesToMedian = new TH1D("fhEntriesToMedian","fhEntriesToMedian",30,0,30);
712 fhCellAreaToMedian = new TH1D("fhCellAreaToMedian", "fhCellAreaToMedian", 75,0,1.5);
715 fOutputList->Add(fhJetPhi);
716 fOutputList->Add(fhTriggerPhi);
717 fOutputList->Add(fhJetEta);
718 fOutputList->Add(fhTriggerEta);
720 fOutputList->Add(fhVertexZ);
721 fOutputList->Add(fhVertexZAccept);
723 fOutputList->Add(fhContribVtx);
724 fOutputList->Add(fhContribVtxAccept);
725 fOutputList->Add(fhDphiTriggerJet);
726 fOutputList->Add(fhDphiTriggerJetAccept);
727 fOutputList->Add(fhCentrality);
728 fOutputList->Add(fhCentralityAccept);
729 fOutputList->Add(fhEntriesToMedian);
730 fOutputList->Add(fhCellAreaToMedian);
732 // raw spectra of INCLUSIVE jets
733 //Centrality, pTjet, A
734 /*const Int_t dimRaw = 3;
735 const Int_t nBinsRaw[dimRaw] = {nBinsCentrality, 50, 100};
736 const Double_t lowBinRaw[dimRaw] = {0.0, 0.0, 0.0};
737 const Double_t hiBinRaw[dimRaw] = {100.0, 100, 1.0};
738 fHJetPtRaw = new THnSparseF("fHJetPtRaw",
739 "Incl. jet spectrum [cent,pTjet,A]",
740 dimRaw,nBinsRaw,lowBinRaw,hiBinRaw);
741 fOutputList->Add(fHJetPtRaw);
743 // raw spectra of LEADING jets
744 //Centrality, pTjet, A
745 fHLeadingJetPtRaw = new THnSparseF("fHLeadingJetPtRaw",
746 "Leading jet spectrum [cent,pTjet,A]",
747 dimRaw,nBinsRaw,lowBinRaw,hiBinRaw);
748 fOutputList->Add(fHLeadingJetPtRaw);
750 // Dphi versus pT jet
751 //Centrality, Dphi=phiTrig-phiJet, pTjet, pTtrigg
752 const Int_t dimDp = 4;
753 const Int_t nBinsDp[dimDp] = {nBinsCentrality, 50, 50, 50};
754 const Double_t lowBinDp[dimDp] = {0.0, -TMath::Pi(), 0.0, 0.0};
755 const Double_t hiBinDp[dimDp] = {100.0, TMath::Pi(), 100.0, 100.0};
756 fHDphiVsJetPtAll = new THnSparseF("fHDphiVsJetPtAll",
757 "Dphi vs jet pT [cent,Dphi,pTjet,pTtrigg]",
758 dimDp,nBinsDp,lowBinDp,hiBinDp);
759 fOutputList->Add(fHDphiVsJetPtAll);
762 //analyze MC generator level
763 if(fIsChargedMC || fIsKine){
765 //Fill response matrix only once
766 fhJetPtGenVsJetPtRec = new TH2D("fhJetPtGenVsJetPtRec","JetPtGenVsJetPtRec", bw*100,0,200, bw*100,0,200);
767 fOutputList->Add(fhJetPtGenVsJetPtRec); //gen MC charg jet pt spectrum versus rec charged jet pt spectrum
769 fhJetPtGenVsJetPtRecSubUeMedian = new TH2D("fhJetPtGenVsJetPtRecSubUeMedian","fhJetPtGenVsJetPtRecSubUeMedian", bw*110,-20,200, bw*110,-20,200);
770 fOutputList->Add(fhJetPtGenVsJetPtRecSubUeMedian); // with kT median bg subtr
772 fhJetPtGenVsJetPtRecSubUeCone=(TH2D*)fhJetPtGenVsJetPtRecSubUeMedian->Clone("fhJetPtGenVsJetPtRecSubUeCone");
773 fhJetPtGenVsJetPtRecSubUeCone->SetTitle("fhJetPtGenVsJetPtRecSubUeCone");
774 fOutputList->Add(fhJetPtGenVsJetPtRecSubUeCone); // with weighted kT median bg subtr
776 fhJetPtGen = new TH1D("fhJetPtGen","Jet Pt (MC Gen)",bw*100,0,200); //MC generator charged jet pt spectrum
777 fOutputList->Add(fhJetPtGen);
779 fhJetPtSubUeMedianGen = new TH1D("fhJetPtSubUeMedianGen","Jet Pt - UE pT (MC Gen)",bw*110,-20,200);
780 fOutputList->Add(fhJetPtSubUeMedianGen); // with kT median bg subtr
782 fhJetPtSubUeConeGen = (TH1D*) fhJetPtSubUeMedianGen->Clone("fhJetPtSubUeConeGen");
783 fOutputList->Add(fhJetPtSubUeConeGen); // with weighted kT median bg subtr
787 fhJetPtGenChargVsJetPtGenFull = new TH2D("fhJetPtGenChargVsJetPtGenFull","fhJetPtGenChargVsJetPtGenFull", 100,0,200, 100,0,200);
788 fOutputList->Add(fhJetPtGenChargVsJetPtGenFull); //gen full MC jet pt versus gen charged jet MC pt
790 fhJetPtGenFull = new TH1D("fhJetPtGenFull","Jet Pt (MC Full jets Gen)",100,0,200); //MC generator full jet pt spectrum
791 fOutputList->Add(fhJetPtGenFull);
795 fh2NtriggersGen = (TH2F*) fh2Ntriggers->Clone("fh2NtriggersGen");
796 fh2NtriggersGen->SetTitle(Form("%s Gen MC",fh2Ntriggers->GetTitle()));
797 fOutputList->Add(fh2NtriggersGen);
799 //Centrality, A, pT, pTtrigg
800 fHJetSpecGen = (THnSparseF*) fHJetSpec->Clone("fHJetSpecGen");
801 fHJetSpecGen->SetTitle(Form("%s Gen MC",fHJetSpec->GetTitle()));
802 fOutputList->Add(fHJetSpecGen);
804 fHJetSpecSubUeMedianGen = (THnSparseF*) fHJetSpecSubUeMedian->Clone("fHJetSpecSubUeMedianGen");
805 fHJetSpecSubUeMedianGen->SetTitle(Form("%s Gen MC",fHJetSpecSubUeMedian->GetTitle()));
806 fOutputList->Add(fHJetSpecSubUeMedianGen);
808 fHJetSpecSubUeConeGen = (THnSparseF*) fHJetSpecSubUeCone->Clone("fHJetSpecSubUeConeGen");
809 fHJetSpecSubUeConeGen->SetTitle(Form("%s Gen MC",fHJetSpecSubUeCone->GetTitle()));
810 fOutputList->Add(fHJetSpecSubUeConeGen);
812 fHJetPhiCorrGen = (THnSparseF*) fHJetPhiCorr->Clone("fHJetPhiCorrGen");
813 fHJetPhiCorrGen->SetTitle(Form("%s Gen MC",fHJetPhiCorr->GetTitle()));
814 fOutputList->Add(fHJetPhiCorrGen);
816 fHJetPhiCorrSubUeMedianGen = (THnSparseF*) fHJetPhiCorrSubUeMedian->Clone("fHJetPhiCorrSubUeMedianGen");
817 fHJetPhiCorrSubUeMedianGen->SetTitle(Form("%s Gen MC",fHJetPhiCorrSubUeMedian->GetTitle()));
818 fOutputList->Add(fHJetPhiCorrSubUeMedianGen);
820 fHJetPhiCorrSubUeConeGen = (THnSparseF*) fHJetPhiCorrSubUeCone->Clone("fHJetPhiCorrSubUeConeGen");
821 fHJetPhiCorrSubUeConeGen->SetTitle(Form("%s Gen MC", fHJetPhiCorrSubUeCone->GetTitle()));
822 fOutputList->Add(fHJetPhiCorrSubUeConeGen);
825 fHJetUeMedianGen = (THnSparseF*) fHJetUeMedian->Clone("fHJetUeMedianGen");
826 fHJetUeMedianGen->SetTitle(Form("%s Gen MC", fHJetUeMedian->GetTitle()));
827 fOutputList->Add(fHJetUeMedianGen);
829 fHJetUeConeGen = (THnSparseF*) fHJetUeCone->Clone("fHJetUeConeGen");
830 fHJetUeConeGen->SetTitle(Form("%s Gen MC", fHJetUeCone->GetTitle()));
831 fOutputList->Add(fHJetUeConeGen);
834 //track efficiency/contamination histograms eta versus pT
835 fhPtTrkTruePrimRec = new TH2D("fhPtTrkTruePrimRec","PtTrkTruePrimRec",100,0,20,18,-0.9,0.9);
836 fOutputList->Add(fhPtTrkTruePrimRec);
838 fhPtTrkTruePrimGen = (TH2D*) fhPtTrkTruePrimRec->Clone("fhPtTrkTruePrimGen");
839 fhPtTrkTruePrimGen->SetTitle("PtTrkTruePrimGen");
840 fOutputList->Add(fhPtTrkTruePrimGen);
842 fhPtTrkSecOrFakeRec = (TH2D*) fhPtTrkTruePrimRec->Clone("fhPtTrkSecOrFakeRec");
843 fhPtTrkSecOrFakeRec->SetTitle("PtTrkSecOrFakeRec");
844 fOutputList->Add(fhPtTrkSecOrFakeRec);
847 fHRhoUeMedianVsConeGen = (THnSparseF*) fHRhoUeMedianVsCone->Clone("hRhoUeMedianVsConeGen");
848 fHRhoUeMedianVsConeGen->SetTitle(Form("%s Gen MC", fHRhoUeMedianVsCone->GetTitle()));
849 fOutputList->Add(fHRhoUeMedianVsConeGen);
851 fhEntriesToMedianGen = (TH1D*) fhEntriesToMedian->Clone("fhEntriesToMedianGen");
852 fhEntriesToMedianGen->SetTitle(Form("%s Gen MC", fhEntriesToMedian->GetTitle()));
853 fOutputList->Add(fhEntriesToMedianGen);
855 fhCellAreaToMedianGen = (TH1D*) fhCellAreaToMedian->Clone("fhCellAreaToMedianGen");
856 fhCellAreaToMedianGen->SetTitle(Form("%s Gen MC", fhCellAreaToMedian->GetTitle()));
857 fOutputList->Add(fhCellAreaToMedianGen);
859 //-------------------------------------
861 const Int_t nBinPt = 150;
862 Double_t binLimitsPt[nBinPt+1];
863 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
865 binLimitsPt[iPt] = -50.0;
867 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.;
871 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
872 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
873 fOutputList->Add(fh1Xsec);
874 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
875 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
876 fOutputList->Add(fh1Trials);
877 fh1AvgTrials = new TH1F("fh1AvgTrials","trials root file",1,0,1);
878 fh1AvgTrials->GetXaxis()->SetBinLabel(1,"#sum{avg ntrials}");
879 fOutputList->Add(fh1AvgTrials);
880 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
881 fOutputList->Add(fh1PtHard);
882 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
883 fOutputList->Add(fh1PtHardNoW);
884 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
885 fOutputList->Add(fh1PtHardTrials);
888 // =========== Switch on Sumw2 for all histos ===========
889 for(Int_t i=0; i<fOutputList->GetEntries(); i++){
890 TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
895 THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
900 TH1::AddDirectory(oldStatus);
902 PostData(1, fOutputList);
905 //--------------------------------------------------------------------
907 void AliAnalysisTaskJetCorePP::UserExec(Option_t *)
913 Double_t eventW = 1.0;
914 Double_t ptHard = 0.0;
915 Double_t nTrials = 1.0; // Trials for MC trigger
916 if(fIsChargedMC || fIsKine) fh1AvgTrials->Fill("#sum{avg ntrials}",fAvgTrials);
918 if(TMath::Abs((Float_t) fJetParamR)<0.00001){
919 AliError("ANTIKT Cone radius is set to zero.");
923 if(TMath::Abs((Float_t) fBgJetParamR)<0.00001){
924 AliError("ANTIKT Cone radius is set to zero.");
928 if(!fIsKine){ //real data or full MC
929 if(!strlen(fJetBranchName.Data())){
930 AliError("Name of jet branch not set.");
933 if(!strlen(fJetBranchNameBg.Data())){
934 AliError("Name of jet bg branch not set.");
938 if(!strlen(fJetBranchNameBgKine.Data())){
939 AliError("Name of jet bg branch for kine not set.");
945 fMcEvent = fMcHandler->MCEvent();
947 if(fDebug > 1) printf("AliAnalysisTaskJetCorePP::Exec() fMcHandler=NULL\n");
948 PostData(1, fOutputList);
952 if(fDebug > 1) printf("AliAnalysisTaskJetCorePP::Exec() fMcEvent=NULL \n");
953 PostData(1, fOutputList);
957 Float_t xsection = 0;
960 AliGenPythiaEventHeader *genPH =
961 dynamic_cast<AliGenPythiaEventHeader*> (fMcEvent->GenEventHeader());
963 xsection = genPH->GetXsection();
964 trials = genPH->Trials();
965 ptHard = genPH->GetPtHard();
967 fh1Xsec->Fill("<#sigma>",xsection);
968 fh1Trials->Fill("#sum{ntrials}",trials);
969 fh1PtHard->Fill(ptHard,eventW);
970 fh1PtHardNoW->Fill(ptHard,1);
971 fh1PtHardTrials->Fill(ptHard,trials);
975 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
977 if(fDebug>1) AliError("ESD not available");
978 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
981 fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
982 AliAODEvent* aod = NULL;
983 // take all other information from the aod we take the tracks from
984 if(!aod && !fIsKine){
985 if(!fESD) aod = fAODIn; //not ESD and not kine => input AOD
986 else aod = fAODOut;// ESD or kine
990 if(fNonStdFile.Length()!=0){
991 // case that we have an AOD extension we can fetch the jets from the extended output
992 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
993 fAODExtension = aodH ? aodH->GetExtension(fNonStdFile.Data()) : 0;
995 if(fDebug>1) Printf("AODExtension found for %s",fNonStdFile.Data());
999 // ----------------- EVENT SELECTION --------------------------
1000 fHistEvtSelection->Fill(1); // number of events before event selection
1003 // physics selection
1004 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
1005 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
1007 if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
1008 if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
1009 fHistEvtSelection->Fill(2);
1010 PostData(1, fOutputList);
1016 if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
1017 fHistEvtSelection->Fill(3);
1018 PostData(1, fOutputList);
1022 // vertex selection for reconstructed data
1023 AliAODVertex* primVtx = aod->GetPrimaryVertex();
1026 if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
1027 fHistEvtSelection->Fill(3);
1028 PostData(1, fOutputList);
1032 Int_t nTracksPrim = primVtx->GetNContributors();
1033 Float_t vtxz = primVtx->GetZ();
1035 fhContribVtx->Fill(nTracksPrim);
1036 if( nTracksPrim > 0 ) fhVertexZ->Fill(vtxz);
1038 if((nTracksPrim < fMinContribVtx) ||
1039 (vtxz < fVtxZMin) ||
1041 if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",
1042 (char*)__FILE__,__LINE__,vtxz);
1043 fHistEvtSelection->Fill(3);
1044 PostData(1, fOutputList);
1048 fhContribVtxAccept->Fill(nTracksPrim);
1049 fhVertexZAccept->Fill(vtxz);
1051 }else{ //KINE cut on vertex
1052 const AliVVertex *vtxMC = fMcEvent->GetPrimaryVertex();
1053 Float_t zVtx = vtxMC->GetZ();
1054 fhVertexZ->Fill(zVtx);
1055 if(zVtx < fVtxZMin || zVtx > fVtxZMax){ //vertex cut
1056 fHistEvtSelection->Fill(3);
1057 PostData(1, fOutputList);
1060 fhVertexZAccept->Fill(zVtx);
1063 //FK// No event class selection imposed
1064 // event class selection (from jet helper task)
1065 //Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
1066 //if(fDebug) Printf("Event class %d", eventClass);
1067 //if(eventClass < fEvtClassMin || eventClass > fEvtClassMax){
1068 // fHistEvtSelection->Fill(4);
1069 // PostData(1, fOutputList);
1073 //------------------ CENTRALITY SELECTION ---------------
1074 AliCentrality *cent = 0x0;
1075 Double_t centValue = 0.0;
1076 if(fSystem){ //fSystem=0 for pp, fSystem=1 for pPb
1078 cent = fESD->GetCentrality();
1079 if(cent) centValue = cent->GetCentralityPercentile("V0M");
1081 centValue = aod->GetHeader()->GetCentrality();
1083 if(fDebug) printf("centrality: %f\n", centValue);
1085 fhCentrality->Fill(centValue);
1087 if(centValue < fCentMin || centValue > fCentMax){
1088 fHistEvtSelection->Fill(4);
1089 PostData(1, fOutputList);
1093 fhCentralityAccept->Fill( centValue );
1097 //-----------------select disjunct event subsamples ----------------
1098 if(!fIsKine){ //reconstructed data
1099 Int_t eventnum = aod->GetHeader()->GetEventNumberESDFile();
1100 Int_t lastdigit = eventnum % 10;
1101 if(!(fEventNumberRangeLow<=lastdigit && lastdigit<=fEventNumberRangeHigh)){
1102 fHistEvtSelection->Fill(5);
1103 PostData(1, fOutputList);
1108 if(fDebug) std::cout<<" ACCEPTED EVENT "<<endl;
1109 fHistEvtSelection->Fill(0); // accepted events
1110 // ==================== end event selection ============================
1112 Double_t tmpArrayFive[5];
1113 Double_t tmpArrayFour[4];
1116 TList particleList; //list of tracks
1117 Int_t indexTrigg = -1;
1118 Double_t rhoFromCellMedian=0.0, rhoCone=0.0;
1121 //=============== Reconstructed antikt jets ===============
1122 ReadTClonesArray(fJetBranchName.Data() , fListJets);
1123 ReadTClonesArray(fJetBranchNameBg.Data(), fListJetsBg);
1125 //============ Estimate background in reconstructed events ===========
1127 //Find Hadron trigger and Estimate rho from cone
1128 indexTrigg = GetListOfTracks(&particleList); //index of trigger hadron in Particle list
1129 EstimateBgCone(fListJets, &particleList, rhoCone);
1131 //Estimate rho from cell median minus jets
1132 EstimateBgRhoMedian(fListJetsBg, &particleList, rhoFromCellMedian,0); //real data
1134 //============== analyze generator level MC ================
1135 TList particleListGen; //list of tracks in MC
1138 if(fIsChargedMC || fIsKine){
1140 if(fIsKine){ //pure kine
1142 //================= generated charged antikt jets ================
1143 ReadTClonesArray(fJetBranchNameKine.Data(), fListJetsGen);
1144 ReadTClonesArray(fJetBranchNameBgKine.Data(), fListJetsBgGen);
1146 //================= generated charged antikt jets ================
1147 ReadTClonesArray(fJetBranchNameChargMC.Data(), fListJetsGen);
1148 ReadTClonesArray(fJetBranchNameBgChargMC.Data(), fListJetsBgGen);
1150 if(fIsFullMC){ //generated full jets
1151 ReadTClonesArray(fJetBranchNameFullMC.Data(), fListJetsGenFull);
1154 //========================================================
1155 //serarch for charged trigger at the MC generator level
1157 Int_t indexTriggGen = -1;
1158 Double_t ptTriggGen = -1;
1159 Int_t iCounterGen = 0; //number of entries in particleListGen array
1160 Int_t triggersMC[200];//list of trigger candidates
1161 Int_t ntriggersMC = 0; //index in triggers array
1164 if(fESD){//ESD input
1166 AliMCEvent* mcEvent = MCEvent();
1168 PostData(1, fOutputList);
1172 AliGenPythiaEventHeader *pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
1173 if(pythiaGenHeader){
1174 nTrials = pythiaGenHeader->Trials();
1175 ptHard = pythiaGenHeader->GetPtHard();
1176 fh1PtHard->Fill(ptHard,eventW);
1177 fh1PtHardNoW->Fill(ptHard,1);
1178 fh1PtHardTrials->Fill(ptHard,nTrials);
1181 for(Int_t it = 0; it < mcEvent->GetNumberOfTracks(); it++){
1182 if(!mcEvent->IsPhysicalPrimary(it)) continue;
1183 AliMCParticle* part = (AliMCParticle*) mcEvent->GetTrack(it);
1185 if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){
1187 if(fHardest==0 && ntriggersMC<200){//single inclusive trigger
1188 if(indexTriggGen > -1){//trigger candidate was found
1189 triggersMC[ntriggersMC] = indexTriggGen;
1194 iCounterGen++;//index in particleListGen array
1199 PostData(1, fOutputList);
1202 TClonesArray *tca = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(AliAODMCParticle::StdBranchName()));
1204 PostData(1, fOutputList);
1207 for(Int_t it = 0; it < tca->GetEntriesFast(); it++){
1208 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1210 if(!part->IsPhysicalPrimary()) continue;
1211 if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){
1213 if(fHardest==0 && ntriggersMC<200){//single inclusive trigger
1214 if(indexTriggGen > -1){ //trigger candidater was found
1215 triggersMC[ntriggersMC] = indexTriggGen;
1220 iCounterGen++;//number of entries in particleListGen array
1224 }else{ //analyze kine
1226 for(Int_t it = 0; it < fMcEvent->GetNumberOfTracks(); it++){
1227 if(!fMcEvent->IsPhysicalPrimary(it)) continue;
1228 AliMCParticle* part = (AliMCParticle*) fMcEvent->GetTrack(it);
1230 if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){
1232 if(fHardest==0 && ntriggersMC<200){//single inclusive trigger
1233 if(indexTriggGen > -1){//trigger candidate was found
1234 triggersMC[ntriggersMC] = indexTriggGen;
1239 iCounterGen++;//index in particleListGen array
1244 //============== Estimate bg in generated events ==============
1245 Double_t rhoFromCellMedianGen=0.0, rhoConeGen=0.0;
1247 //Estimate rho from cone
1248 EstimateBgCone(fListJetsGen, &particleListGen, rhoConeGen);
1250 //Estimate rho from cell median minus jets
1251 EstimateBgRhoMedian(fListJetsBgGen, &particleListGen, rhoFromCellMedianGen,1);//mc data
1253 //============ Generator trigger+jet ==================
1254 if(fHardest==0){ //single inclusive trigger
1255 if(ntriggersMC>0){ //there is at least one trigger
1256 Int_t rnd = fRandom->Integer(ntriggersMC); //0 to ntriggers-1
1257 indexTriggGen = triggersMC[rnd];
1259 indexTriggGen = -1; //trigger not found
1264 Int_t ilowGen = (fHardest==0 || fHardest==1) ? indexTriggGen : 0;
1265 Int_t ihighGen = (fHardest==0 || fHardest==1) ? indexTriggGen+1 : particleListGen.GetEntries();
1266 Bool_t fillOnceGen = kTRUE;
1269 for(Int_t igen= ilowGen; igen< ihighGen; igen++){ //loop over possible trigger
1270 indexTriggGen = igen; //trigger hadron
1272 if(indexTriggGen == -1) continue;
1273 AliVParticle* triggerGen=(AliVParticle*) particleListGen.At(indexTriggGen);
1274 if(!triggerGen) continue;
1277 if(triggerGen->Pt() < 6.0) continue; //all hadrons pt>6
1279 if(TMath::Abs((Float_t) triggerGen->Eta()) > fTriggerEtaCut) continue;
1281 ptTriggGen = triggerGen->Pt(); //count triggers
1282 fh2NtriggersGen->Fill(centValue, ptTriggGen);
1284 //Count jets and trigger-jet pairs at MC generator level
1285 for(Int_t ij=0; ij<fListJetsGen->GetEntries(); ij++){
1286 AliAODJet* jet = (AliAODJet*)(fListJetsGen->At(ij));
1288 Double_t etaJetGen = jet->Eta();
1289 Double_t areaJetGen = jet->EffectiveAreaCharged();
1291 if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){
1293 Double_t ptUeFromCellMedianGen = rhoFromCellMedianGen*areaJetGen;
1294 Double_t ptUeConeGen = rhoConeGen*areaJetGen;
1296 //Rongrong's analysis
1297 Double_t dPhiGen = jet->Phi() - triggerGen->Phi();
1298 if(dPhiGen>2*TMath::Pi()) dPhiGen -= 2*TMath::Pi();
1299 if(dPhiGen<-2*TMath::Pi()) dPhiGen += 2*TMath::Pi();
1300 if(dPhiGen<-0.5*TMath::Pi()) dPhiGen += 2*TMath::Pi();
1301 if(dPhiGen>1.5*TMath::Pi()) dPhiGen -= 2*TMath::Pi();
1303 //[A,pTjet,pTtrig,dphi]
1304 tmpArrayFour[0] = areaJetGen;
1305 tmpArrayFour[1] = jet->Pt();
1306 tmpArrayFour[2] = ptTriggGen;
1307 tmpArrayFour[3] = dPhiGen;
1309 fHJetPhiCorrGen->Fill(tmpArrayFour); // Dphi distribution jet-triger
1311 //[A,pTjet-pTUe,pTtrig,dphi]
1312 tmpArrayFour[1] = jet->Pt() - ptUeFromCellMedianGen;
1313 fHJetPhiCorrSubUeMedianGen->Fill(tmpArrayFour);
1315 //[A,pTjet-pTUe,pTtrig,dphi]
1316 tmpArrayFour[1] = jet->Pt() - ptUeConeGen;
1317 fHJetPhiCorrSubUeConeGen->Fill(tmpArrayFour);
1319 //Leticia's analysis
1320 Double_t dphi = RelativePhi(triggerGen->Phi(), jet->Phi());
1321 if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue;
1323 //Centrality, A, pT, pTtrigg
1324 tmpArrayFive[0] = centValue;
1325 tmpArrayFive[1] = areaJetGen;
1326 tmpArrayFive[2] = jet->Pt();
1327 tmpArrayFive[3] = ptTriggGen;
1328 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1329 fHJetSpecGen->Fill(tmpArrayFive);
1331 //Centrality, A, pTjet-pTbgCellMedian, pTtrigg, dphi
1332 tmpArrayFive[0] = centValue;
1333 tmpArrayFive[1] = areaJetGen;
1334 tmpArrayFive[2] = jet->Pt() - ptUeFromCellMedianGen;
1335 tmpArrayFive[3] = ptTriggGen;
1336 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1337 fHJetSpecSubUeMedianGen->Fill(tmpArrayFive);
1339 //Centrality, A, pTjet-pTbgCone, pTtrigg, dphi
1340 tmpArrayFive[0] = centValue;
1341 tmpArrayFive[1] = areaJetGen;
1342 tmpArrayFive[2] = jet->Pt() - ptUeConeGen;
1343 tmpArrayFive[3] = ptTriggGen;
1344 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1345 fHJetSpecSubUeConeGen->Fill(tmpArrayFive);
1347 //Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", kT median
1348 tmpArrayFive[0] = areaJetGen;
1349 tmpArrayFive[1] = jet->Pt();
1350 tmpArrayFive[2] = jet->Pt() - ptUeFromCellMedianGen;
1351 tmpArrayFive[3] = ptUeFromCellMedianGen;
1352 tmpArrayFive[4] = rhoFromCellMedianGen;
1353 fHJetUeMedianGen->Fill(tmpArrayFive);
1355 //Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", perpendicular Cone
1356 tmpArrayFive[0] = areaJetGen;
1357 tmpArrayFive[1] = jet->Pt();
1358 tmpArrayFive[2] = jet->Pt() - ptUeConeGen;
1359 tmpArrayFive[3] = ptUeConeGen;
1360 tmpArrayFive[4] = rhoConeGen;
1361 fHJetUeConeGen->Fill(tmpArrayFive);
1364 Double_t fillRhoGen[] = { rhoConeGen, rhoFromCellMedianGen};
1365 fHRhoUeMedianVsConeGen->Fill(fillRhoGen);
1366 fillOnceGen = kFALSE;
1368 }//back to back jet-trigger pair
1372 if(fFillRespMx && !fIsKine){
1373 //================ RESPONSE MATRIX ===============
1374 //Count jets and trigger-jet pairs at MC generator level
1375 for(Int_t ij=0; ij<fListJetsGen->GetEntries(); ij++){
1376 AliAODJet* jet = (AliAODJet*)(fListJetsGen->At(ij));
1378 Double_t etaJetGen = jet->Eta();
1379 Double_t ptJetGen = jet->Pt();
1381 if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){
1382 fhJetPtGen->Fill(ptJetGen); // gen level pt spectum of jets response mx normalization
1384 Double_t areaJetGen = jet->EffectiveAreaCharged();
1385 Double_t ptUeFromCellMedianGen = rhoFromCellMedianGen*areaJetGen;
1386 Double_t ptUeConeGen = rhoConeGen*areaJetGen;
1388 fhJetPtSubUeMedianGen->Fill(ptJetGen - ptUeFromCellMedianGen);
1389 fhJetPtSubUeConeGen->Fill(ptJetGen - ptUeConeGen);
1392 if(fListJets->GetEntries()>0 && fListJetsGen->GetEntries()>0){ //at least some reconstructed jets
1394 Int_t ng = (Int_t) fListJetsGen->GetEntries();
1395 Int_t nr = (Int_t) fListJets->GetEntries();
1397 //Find closest MC generator - reconstructed jets
1398 if(faGenIndex.GetSize()<nr) faGenIndex.Set(nr); //idx of gen jet assoc to rec jet
1399 if(faRecIndex.GetSize()<ng) faRecIndex.Set(ng); //idx of rec jet assoc to gen jet
1402 Printf("New Rec List %d gen index Array %d",nr,faGenIndex.GetSize());
1403 Printf("New Gen List %d rec index Array %d",ng,faRecIndex.GetSize());
1405 //matching of MC genrator level and reconstructed jets
1406 AliAnalysisHelperJetTasks::GetClosestJets(fListJetsGen,ng,fListJets,nr,faGenIndex,faRecIndex,fDebug);
1408 // Fill response matrix
1409 for(Int_t ir = 0; ir < nr; ir++){
1410 AliAODJet *recJet = (AliAODJet*) fListJets->At(ir);
1411 Double_t etaJetRec = recJet->Eta();
1412 Double_t ptJetRec = recJet->Pt();
1413 Double_t areaJetRec = recJet->EffectiveAreaCharged();
1414 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1416 if((fJetEtaMin <= etaJetRec) && (etaJetRec <= fJetEtaMax)){
1417 Int_t ig = faGenIndex[ir]; //associated generator level jet
1418 if(ig >= 0 && ig < ng){
1419 if(fDebug > 10) Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1420 AliAODJet *genJet = (AliAODJet*) fListJetsGen->At(ig);
1421 Double_t ptJetGen = genJet->Pt();
1422 Double_t etaJetGen = genJet->Eta();
1424 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1425 if((fJetEtaMin <= etaJetGen) && (etaJetGen <= fJetEtaMax)){
1426 fhJetPtGenVsJetPtRec->Fill(ptJetRec, ptJetGen);
1428 Double_t areaJetGen = genJet->EffectiveAreaCharged();
1429 Double_t ptUeFromCellMedianGen = rhoFromCellMedianGen*areaJetGen;
1430 Double_t ptUeConeGen = rhoConeGen*areaJetGen;
1431 Double_t ptUeFromCellMedianRec = rhoFromCellMedian*areaJetRec;
1432 Double_t ptUeConeRec = rhoCone*areaJetRec;
1433 fhJetPtGenVsJetPtRecSubUeMedian->Fill(ptJetRec-ptUeFromCellMedianRec,
1434 ptJetGen-ptUeFromCellMedianGen);
1435 fhJetPtGenVsJetPtRecSubUeCone->Fill(ptJetRec-ptUeConeRec, ptJetGen-ptUeConeGen);
1438 }//rec jet in eta acceptance
1439 }//loop over reconstructed jets
1440 }// # of rec jets >0
1442 //=========================== Full jet vs charged jet matrix ==========================
1444 //Count full jets and charged-jet pairs at MC generator level
1445 for(Int_t ij=0; ij<fListJetsGenFull->GetEntries(); ij++){
1446 AliAODJet* jetFull = (AliAODJet*)(fListJetsGenFull->At(ij));
1447 if(!jetFull) continue;
1448 Double_t etaJetFull = jetFull->Eta();
1449 Double_t ptJetFull = jetFull->Pt();
1451 if((fJetEtaMin<=etaJetFull) && (etaJetFull<=fJetEtaMax)){
1452 fhJetPtGenFull->Fill(ptJetFull); // generator level pt spectum of full jets
1455 if(fListJetsGen->GetEntries()>0 && fListJetsGenFull->GetEntries()>0){ //at least some reconstructed jets
1456 Int_t nful = (Int_t) fListJetsGenFull->GetEntries();
1457 Int_t nchr = (Int_t) fListJetsGen->GetEntries();
1459 //Find closest MC generator full - charged jet
1460 if(faGenIndex.GetSize()<nchr) faGenIndex.Set(nchr); //idx of gen FULL jet assoc to gen CHARGED jet
1461 if(faRecIndex.GetSize()<nful) faRecIndex.Set(nful); //idx of gen CHARGED jet assoc to gen FULL jet
1464 Printf("New Charg List %d Full index Array %d",nchr,faGenIndex.GetSize());
1465 Printf("New Full List %d Charg index Array %d",nful,faRecIndex.GetSize());
1467 //matching of MC genrator level and reconstructed jets
1468 AliAnalysisHelperJetTasks::GetClosestJets(fListJetsGenFull,nful,fListJetsGen,nchr,faGenIndex,faRecIndex,fDebug);
1470 // Fill response matrix
1471 for(Int_t ichr = 0; ichr < nchr; ichr++){ //charged jet loop
1472 AliAODJet *chJet = (AliAODJet*) fListJetsGen->At(ichr);
1473 Double_t etaJetCh = chJet->Eta();
1474 Double_t ptJetCh = chJet->Pt();
1475 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1477 if((fJetEtaMin <= etaJetCh) && (etaJetCh <= fJetEtaMax)){
1478 Int_t iful = faGenIndex[ichr]; //associated generator level jet
1479 if(iful >= 0 && iful < nful){
1480 if(fDebug > 10) Printf("%s:%d iful = %d ichr = %d",(char*)__FILE__,__LINE__,iful,ichr);
1481 AliAODJet *genJetFull = (AliAODJet*) fListJetsGenFull->At(iful);
1482 Double_t ptJetFull = genJetFull->Pt();
1483 Double_t etaJetFull = genJetFull->Eta();
1485 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1486 if((fJetEtaMin <= etaJetFull) && (etaJetFull <= fJetEtaMax)){
1487 fhJetPtGenChargVsJetPtGenFull->Fill(ptJetFull,ptJetCh);
1490 }//rec jet in eta acceptance
1491 }//loop over reconstructed jets
1492 }// # of rec jets >0
1493 }//pointer MC generator jets
1494 } //fill resp mx only for bin
1495 }//analyze generator level MC
1498 if(fIsKine){ //skip reconstructed data analysis in case of kine
1499 PostData(1, fOutputList);
1502 //============= RECONSTRUCTED INCLUSIVE JETS ===============
1504 Double_t etaJet = 0.0;
1505 Double_t pTJet = 0.0;
1506 Double_t areaJet = 0.0;
1507 Double_t phiJet = 0.0;
1508 //Int_t indexLeadingJet = -1;
1509 //Double_t pTLeadingJet = -10.0;
1510 //Double_t areaLeadingJet = -10.0;
1512 for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
1513 AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
1515 etaJet = jet->Eta();
1516 phiJet = jet->Phi();
1518 if(pTJet==0) continue;
1520 if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
1521 /*areaJet = jet->EffectiveAreaCharged();*/
1523 //Jet Diagnostics---------------------------------
1524 fhJetPhi->Fill(pTJet, RelativePhi(phiJet,0.0)); //phi -pi,pi
1525 fhJetEta->Fill(pTJet, etaJet);
1526 //search for leading jet
1527 /*if(pTJet > pTLeadingJet){
1528 indexLeadingJet = ij;
1529 pTLeadingJet = pTJet;
1530 areaLeadingJet = areaJet;
1533 // raw spectra of INCLUSIVE jets
1534 //Centrality, pTjet, A
1535 Double_t fillraw[] = { centValue,
1539 fHJetPtRaw->Fill(fillraw);*/
1542 if(indexLeadingJet > -1){
1543 // raw spectra of LEADING jets
1544 //Centrality, pTjet, A
1545 Double_t fillleading[] = { centValue,
1549 fHLeadingJetPtRaw->Fill(fillleading);
1553 // =============== RECONSTRUCTED TRIGGER-JET PAIRS ================
1554 if(fIsChargedMC && fFillRespMx){
1555 FillEffHistos(&particleList, &particleListGen); //Fill efficiency histos
1557 Bool_t filledOnce = kTRUE; //fill rho histogram only once per event
1558 //set ranges of the trigger loop
1559 Int_t ilow = (fHardest==0 || fHardest==1) ? indexTrigg : 0;
1560 Int_t ihigh = (fHardest==0 || fHardest==1) ? indexTrigg+1 : particleList.GetEntries();
1562 for(Int_t itrk= ilow; itrk< ihigh; itrk++){ //loop over possible trigger
1563 indexTrigg = itrk; //trigger hadron with pT >10 GeV
1565 if(indexTrigg < 0) continue;
1567 AliVParticle *triggerHadron = (AliVParticle*) particleList.At(indexTrigg);
1569 PostData(1, fOutputList);
1573 if(triggerHadron->Pt() < 6.0) continue; //all hadrons pt>6
1575 if(TMath::Abs((Float_t) triggerHadron->Eta())> fTriggerEtaCut) continue;
1577 //Fill trigger histograms
1578 fh2Ntriggers->Fill(centValue,triggerHadron->Pt()); //trigger pT
1579 fhTriggerPhi->Fill(triggerHadron->Pt(),RelativePhi(triggerHadron->Phi(),0.0)); //phi -pi,pi
1580 fhTriggerEta->Fill(triggerHadron->Pt(),triggerHadron->Eta());
1583 //---------- make trigger-jet pairs ---------
1587 for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
1588 AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
1590 etaJet = jet->Eta();
1591 phiJet = jet->Phi();
1593 if(pTJet==0) continue;
1595 if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
1596 areaJet = jet->EffectiveAreaCharged();
1598 //subtract bg using estinates base on median of kT jets
1599 Double_t ptUeFromCellMedian = rhoFromCellMedian*areaJet;
1600 Double_t ptUeCone = rhoCone*areaJet;
1602 //if(areaJet >= 0.07) injet++;
1603 //if(areaJet >= 0.4) injet4++;
1605 Double_t dphi = RelativePhi(triggerHadron->Phi(), phiJet);
1606 fhDphiTriggerJet->Fill(dphi); //Input
1608 //Dphi versus jet pT
1609 //Centrality, Dphi=phiTrig-phiJet, pTjet, pTtrigg
1610 /*Double_t filldp[] = { centValue,
1615 fHDphiVsJetPtAll->Fill(filldp);
1617 //Rongrong's analysis
1619 Double_t dPhi = phiJet - triggerHadron->Phi();
1620 if(dPhi>2*TMath::Pi()) dPhi -= 2*TMath::Pi();
1621 if(dPhi<-2*TMath::Pi()) dPhi += 2*TMath::Pi();
1622 if(dPhi<-0.5*TMath::Pi()) dPhi += 2*TMath::Pi();
1623 if(dPhi>1.5*TMath::Pi()) dPhi -= 2*TMath::Pi();
1625 //[A,pTjet,pTtrig,dphi]
1626 tmpArrayFour[0] = areaJet;
1627 tmpArrayFour[1] = pTJet;
1628 tmpArrayFour[2] = triggerHadron->Pt();
1629 tmpArrayFour[3] = dPhi;
1631 fHJetPhiCorr->Fill(tmpArrayFour); // Dphi distribution jet-triger
1633 //[A,pTjet-pTUe,pTtrig,dphi]
1634 tmpArrayFour[1] = pTJet - ptUeFromCellMedian;
1635 fHJetPhiCorrSubUeMedian->Fill(tmpArrayFour);
1637 //[A,pTjet-pTUe,pTtrig,dphi]
1638 tmpArrayFour[1] = pTJet - ptUeCone;
1639 fHJetPhiCorrSubUeCone->Fill(tmpArrayFour);
1641 //Leticia's analysis
1643 // Select back to back trigger - jet pairs
1644 if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue;
1645 fhDphiTriggerJetAccept->Fill(dphi); //Accepted
1648 //Centrality, A, pTjet, pTtrigg, dphi
1649 tmpArrayFive[0] = centValue;
1650 tmpArrayFive[1] = areaJet;
1651 tmpArrayFive[2] = pTJet;
1652 tmpArrayFive[3] = triggerHadron->Pt();
1653 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1654 fHJetSpec->Fill(tmpArrayFive);
1656 //Centrality, A, pTjet-pTbgCellMedian, pTtrigg, dphi
1657 tmpArrayFive[0] = centValue;
1658 tmpArrayFive[1] = areaJet;
1659 tmpArrayFive[2] = pTJet-ptUeFromCellMedian;
1660 tmpArrayFive[3] = triggerHadron->Pt();
1661 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1662 fHJetSpecSubUeMedian->Fill(tmpArrayFive);
1664 //Centrality, A, pTjet-pTbgCone, pTtrigg, dphi
1665 tmpArrayFive[0] = centValue;
1666 tmpArrayFive[1] = areaJet;
1667 tmpArrayFive[2] = pTJet - ptUeCone;
1668 tmpArrayFive[3] = triggerHadron->Pt();
1669 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1670 fHJetSpecSubUeCone->Fill(tmpArrayFive);
1672 //Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", kT median
1673 tmpArrayFive[0] = areaJet;
1674 tmpArrayFive[1] = pTJet;
1675 tmpArrayFive[2] = pTJet - ptUeFromCellMedian;
1676 tmpArrayFive[3] = ptUeFromCellMedian;
1677 tmpArrayFive[4] = rhoFromCellMedian;
1678 fHJetUeMedian->Fill(tmpArrayFive);
1680 //Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", Cone median
1681 tmpArrayFive[0] = areaJet;
1682 tmpArrayFive[1] = pTJet;
1683 tmpArrayFive[2] = pTJet - ptUeCone;
1684 tmpArrayFive[3] = ptUeCone;
1685 tmpArrayFive[4] = rhoCone;
1686 fHJetUeCone->Fill(tmpArrayFive);
1688 if(filledOnce){ //fill for each event only once
1689 Double_t fillRho[] = { rhoCone,rhoFromCellMedian};
1690 fHRhoUeMedianVsCone->Fill(fillRho);
1691 filledOnce = kFALSE;
1695 //Fill Jet Density In the Event A>0.07
1697 Double_t filldens[]={ (Double_t) particleList.GetEntries(),
1701 fHJetDensity->Fill(filldens);
1704 //Fill Jet Density In the Event A>0.4
1706 Double_t filldens4[]={ (Double_t) particleList.GetEntries(),
1707 injet4/fkAcceptance,
1710 fHJetDensityA4->Fill(filldens4);
1715 PostData(1, fOutputList);
1718 //----------------------------------------------------------------------------
1719 void AliAnalysisTaskJetCorePP::Terminate(const Option_t *)
1721 // Draw result to the screen
1722 // Called once at the end of the query
1724 if(fDebug) printf("AliAnalysisTaskJetCorePP DONE\n");
1725 if(!GetOutputData(1)) return;
1729 //----------------------------------------------------------------------------
1730 Int_t AliAnalysisTaskJetCorePP::GetListOfTracks(TList *list){
1731 //Fill the list of accepted tracks (passed track cut)
1732 //return consecutive index of the hardest ch hadron in the list
1734 AliAODEvent *aodevt = NULL;
1736 if(!fESD) aodevt = fAODIn;
1737 else aodevt = fAODOut;
1739 if(!aodevt) return -1;
1741 Int_t index = -1; //index of the highest particle in the list
1742 Double_t ptmax = -10;
1743 Int_t triggers[200];
1744 Int_t ntriggers = 0; //index in triggers array
1746 for(Int_t it = 0; it < aodevt->GetNumberOfTracks(); it++){
1747 AliAODTrack *tr = aodevt->GetTrack(it);
1749 if((fFilterMask > 0) && !(tr->TestFilterBit(fFilterMask))) continue;
1750 //if((fFilterMask > 0) && !(tr->IsHybridGlobalConstrainedGlobal())) continue;
1751 if(TMath::Abs((Float_t) tr->Eta()) > fTrackEtaCut) continue;
1752 if(tr->Pt() < fTrackLowPtCut) continue;
1756 if(fHardest>0){ //leading particle
1763 if(fHardest==0 && ntriggers<200){ //single inclusive
1764 if(fTriggerPtRangeLow <= tr->Pt() &&
1765 tr->Pt() < fTriggerPtRangeHigh){
1766 triggers[ntriggers] = iCount;
1774 if(fHardest==0 && ntriggers>0){ //select random inclusive trigger
1775 Int_t rnd = fRandom->Integer(ntriggers); //0 to ntriggers-1
1776 index = triggers[rnd];
1782 //----------------------------------------------------------------------------
1784 Double_t AliAnalysisTaskJetCorePP::GetBackgroundInPerpCone(Float_t jetR, Double_t jetPhi, Double_t jetEta, TList* trkList){
1785 //calculate sum of track pT in the cone perpendicular in phi to the jet
1786 //jetR = cone radius
1787 // jetPhi, jetEta = direction of the jet
1788 Int_t numberOfTrks = trkList->GetEntries();
1789 Double_t pTsum = 0.0;
1790 Double_t perpConePhi = jetPhi + TMath::Pi()/2;//perp cone w.r.t. jet in phi
1791 for(Int_t it=0; it<numberOfTrks; it++){
1792 AliVParticle *trk = (AliVParticle*) trkList->At(it);
1793 Double_t dphi = RelativePhi(perpConePhi,trk->Phi());
1794 Double_t deta = trk->Eta()-jetEta;
1796 if( (dphi*dphi + deta*deta)< (jetR*jetR)){
1804 //----------------------------------------------------------------------------
1806 Double_t AliAnalysisTaskJetCorePP::RelativePhi(Double_t mphi,Double_t vphi){
1807 //Get relative azimuthal angle of two particles -pi to pi
1808 if (vphi < -TMath::Pi()) vphi += TMath::TwoPi();
1809 else if (vphi > TMath::Pi()) vphi -= TMath::TwoPi();
1811 if (mphi < -TMath::Pi()) mphi += TMath::TwoPi();
1812 else if (mphi > TMath::Pi()) mphi -= TMath::TwoPi();
1814 Double_t dphi = mphi - vphi;
1815 if (dphi < -TMath::Pi()) dphi += TMath::TwoPi();
1816 else if (dphi > TMath::Pi()) dphi -= TMath::TwoPi();
1818 return dphi;//dphi in [-Pi, Pi]
1822 //----------------------------------------------------------------------------
1823 Bool_t AliAnalysisTaskJetCorePP::SelectMCGenTracks(AliVParticle *trk, TList *trkList, Double_t &ptLeading, Int_t &index, Int_t counter){
1824 //fill track efficiency denominator
1825 if(trk->Pt() < fTrackLowPtCut) return kFALSE;
1826 if(trk->Charge()==0) return kFALSE;
1827 if(TMath::Abs((Float_t) trk->Eta()) > fTrackEtaCut) return kFALSE;
1829 if(fFillRespMx) fhPtTrkTruePrimGen->Fill(trk->Pt(),trk->Eta()); //Efficiency denominator
1832 if(fHardest>0){ //leading particle
1833 if(ptLeading < trk->Pt()){
1835 ptLeading = trk->Pt();
1839 if(fHardest==0){ //single inclusive
1841 if(fTriggerPtRangeLow <= trk->Pt() &&
1842 trk->Pt() < fTriggerPtRangeHigh){
1850 //----------------------------------------------------------------------------
1851 void AliAnalysisTaskJetCorePP::FillEffHistos(TList *recList, TList *genList){
1852 //fill track efficiency numerator
1853 if(!recList) return;
1854 if(!genList) return;
1855 Int_t nRec = recList->GetEntries();
1856 Int_t nGen = genList->GetEntries();
1857 for(Int_t ir=0; ir<nRec; ir++){
1858 AliAODTrack *trkRec = (AliAODTrack*) recList->At(ir);
1859 if(!trkRec) continue;
1860 Int_t recLabel = TMath::Abs(trkRec->GetLabel());
1861 Bool_t matched = kFALSE;
1863 for(Int_t ig=0; ig<nGen; ig++){
1864 AliVParticle *trkGen = (AliVParticle*) genList->At(ig);
1865 if(!trkGen) continue;
1866 Int_t genLabel = TMath::Abs(trkGen->GetLabel());
1867 if(recLabel==genLabel){
1868 fhPtTrkTruePrimRec->Fill(trkGen->Pt(),trkGen->Eta());
1873 if(!matched) fhPtTrkSecOrFakeRec->Fill(trkRec->Pt(),trkRec->Eta());
1878 //________________________________________________________________
1879 void AliAnalysisTaskJetCorePP::EstimateBgRhoMedian(TList *listJet, TList* listPart, Double_t &rhoMedian, Int_t mode){
1880 //Estimate background rho by means of integrating track pT outside identified jet cones
1884 //identify leading jet in the full track acceptance
1885 Int_t idxLeadingJet = -1;
1886 Double_t pTleading = 0.0;
1887 AliAODJet* jet = NULL;
1889 for(Int_t ij = 0; ij < listJet->GetEntries(); ij++){ //loop over bg kT jets
1890 jet = (AliAODJet*)(listJet->At(ij));
1892 if(TMath::Abs(jet->Eta()) > fTrackEtaCut) continue; //track acceptance cut
1893 if(jet->Pt() > pTleading){
1895 pTleading = jet->Pt();
1900 static Double_t jphi[100];
1901 static Double_t jeta[100];
1902 static Double_t jRsquared[100];
1904 for(Int_t ij = 0; ij< listJet->GetEntries(); ij++){ //loop over bg kT jets
1906 jet = (AliAODJet*)(listJet->At(ij));
1908 if(TMath::Abs(jet->Eta()) > fTrackEtaCut) continue;
1910 if((jet->Pt() > fBgMaxJetPt || ij== idxLeadingJet) && njets<100){
1911 jphi[njets] = RelativePhi(jet->Phi(),0.0); //-pi,pi
1912 jeta[njets] = jet->Eta();
1913 jRsquared[njets] = fSafetyMargin*jet->EffectiveAreaCharged()/TMath::Pi(); //scale jet radius
1919 static Double_t nOutCone[10][4];
1920 static Double_t sumPtOutOfCone[10][4];
1922 for(Int_t ie=0; ie<fnEta; ie++){
1923 for(Int_t ip=0; ip<fnPhi; ip++){
1924 nOutCone[ip][ie] = 0.0; //initialize counter
1925 sumPtOutOfCone[ip][ie] = 0.0;
1929 Double_t rndphi, rndeta;
1930 Double_t rndphishift, rndetashift;
1931 Double_t dphi, deta;
1934 for(Int_t it=0; it<fnTrials; it++){
1936 rndphi = fRandom->Uniform(0, fPhiSize);
1937 rndeta = fRandom->Uniform(0, fEtaSize);
1939 for(Int_t ip=0; ip<fnPhi; ip++){ //move radom position to each cell
1940 rndphishift = rndphi + ip*fPhiSize - TMath::Pi();
1941 for(Int_t ie=0; ie<fnEta; ie++){
1942 rndetashift = rndeta + ie*fEtaSize - fEtaSize;
1944 bIsInCone = 0; //tag if trial is in the jet cone
1945 for(Int_t ij=0; ij<njets; ij++){
1946 deta = jeta[ij] - rndetashift;
1947 dphi = RelativePhi(rndphishift,jphi[ij]);
1948 if((dphi*dphi + deta*deta) < jRsquared[ij]){
1953 if(!bIsInCone) nOutCone[ip][ie]++;
1959 //Sum particle pT outside jets in cells
1960 Int_t npart = listPart->GetEntries();
1961 Int_t phicell,etacell;
1962 AliVParticle *part = NULL;
1963 for(Int_t ip=0; ip < npart; ip++){
1965 part = (AliVParticle*) listPart->At(ip);
1967 if( TMath::Abs(part->Eta()) > fEtaSize) continue; //skip part outside +-0.7 in eta
1969 bIsInCone = 0; //init
1970 for(Int_t ij=0; ij<njets; ij++){
1971 dphi = RelativePhi(jphi[ij], part->Phi());
1972 deta = jeta[ij] - part->Eta();
1973 if((dphi*dphi + deta*deta) < jRsquared[ij]){
1979 phicell = TMath::Nint(TMath::Floor((RelativePhi(part->Phi(),0.0) + TMath::Pi())/fPhiSize));
1980 etacell = TMath::Nint(TMath::Floor((part->Eta()+fEtaSize)/fEtaSize));
1981 sumPtOutOfCone[phicell][etacell]+= part->Pt();
1985 static Double_t rhoInCells[20];
1986 Double_t relativeArea;
1988 Double_t bufferArea=0.0, bufferPt=0.0; //sum cells where A< fJetFreeAreaFrac
1989 for(Int_t ip=0; ip<fnPhi; ip++){
1990 for(Int_t ie=0; ie<fnEta; ie++){
1991 relativeArea = nOutCone[ip][ie]/fnTrials;
1992 //cout<<"RA= "<< relativeArea<<" BA="<<bufferArea<<" BPT="<< bufferPt <<endl;
1994 bufferArea += relativeArea;
1995 bufferPt += sumPtOutOfCone[ip][ie];
1996 if(bufferArea > fJetFreeAreaFrac){
1997 rhoInCells[nCells] = bufferPt/(bufferArea*fCellArea);
1999 if(mode==1) fhCellAreaToMedianGen->Fill(bufferArea*fCellArea);
2000 else fhCellAreaToMedian->Fill(bufferArea*fCellArea);
2011 rhoMedian = TMath::Median(nCells, rhoInCells);
2012 if(mode==1){ //mc data
2013 fhEntriesToMedianGen->Fill(nCells);
2015 fhEntriesToMedian->Fill(nCells);
2021 //__________________________________________________________________
2022 void AliAnalysisTaskJetCorePP::EstimateBgCone(TList *listJet, TList* listPart, Double_t &rhoPerpCone){
2024 rhoPerpCone = 0.0; //init
2026 if(!listJet) return;
2027 Int_t njet = listJet->GetEntries();
2028 Int_t npart = listPart->GetEntries();
2029 Double_t pTleading = 0.0;
2030 Double_t phiLeading = 1000.;
2031 Double_t etaLeading = 1000.;
2033 for(Int_t jsig=0; jsig < njet; jsig++){
2034 AliAODJet* signaljet = (AliAODJet*)(listJet->At(jsig));
2035 if(!signaljet) continue;
2036 if((signaljet->Eta()<fJetEtaMin) && (fJetEtaMax<signaljet->Eta())) continue; //acceptance cut
2037 if(signaljet->Pt() >= pTleading){ //replace leading and subleading jet
2038 pTleading = signaljet->Pt();
2039 phiLeading = signaljet->Phi();
2040 etaLeading = signaljet->Eta();
2044 Double_t phileftcone = phiLeading + TMath::Pi()/2;
2045 Double_t phirightcone = phiLeading - TMath::Pi()/2;
2048 for(Int_t ip=0; ip < npart; ip++){
2050 AliVParticle *part = (AliVParticle*) listPart->At(ip);
2055 dp = RelativePhi(phileftcone, part->Phi());
2056 de = etaLeading - part->Eta();
2057 if( sqrt(dp*dp + de*de)< fBgConeR) rhoPerpCone += part->Pt();
2059 dp = RelativePhi(phirightcone, part->Phi());
2060 if( sqrt(dp*dp + de*de)< fBgConeR) rhoPerpCone += part->Pt();
2064 //normalize total pT by two times cone are
2065 rhoPerpCone = rhoPerpCone/(2*TMath::Pi()*fBgConeR*fBgConeR);
2068 //__________________________________________________________________
2070 void AliAnalysisTaskJetCorePP::ReadTClonesArray(TString bname, TList *list){
2071 //Convert TClones array of jets to TList
2073 if(!strlen(bname.Data())){
2074 AliError(Form("Jet branch %s not set.", bname.Data()));
2078 TClonesArray *array=0x0;
2079 if(fUseExchContainer){ //take input from exchange containers
2080 if(fJetBranchNameKine.Length()>0 && bname.CompareTo(fJetBranchNameKine.Data())==0){
2081 array = dynamic_cast<TClonesArray*>(GetInputData(1)); //connect exchange container slot 1
2083 if(fJetBranchNameBgKine.Length()>0 && bname.CompareTo(fJetBranchNameBgKine.Data())==0){
2084 array = dynamic_cast<TClonesArray*>(GetInputData(2)); //connect exchange container slot 2
2086 }else{ //take input from AOD
2087 if(fAODOut&&!array){
2088 array = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(bname.Data()));
2090 if(fAODExtension&&!array){
2091 array = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(bname.Data()));
2094 array = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(bname.Data()));
2098 list->Clear(); //clear list beforehand
2102 Printf("########## %s: %d jets",bname.Data(), array->GetEntriesFast());
2105 for(Int_t iJet = 0; iJet < array->GetEntriesFast(); iJet++) {
2106 AliAODJet *jet = dynamic_cast<AliAODJet*>((*array)[iJet]);
2107 if (jet) list->Add(jet);