#include "TChain.h"
+#include "TClonesArray.h"
#include "TTree.h"
#include "TList.h"
#include "TMath.h"
#include "AliAODEvent.h"
#include "AliAODHandler.h"
#include "AliAODJet.h"
-
+#include "AliVVertex.h"
#include "AliAnalysisTaskJetCorePP.h"
+#include "AliHeader.h" //KF//
using std::cout;
using std::endl;
fAODIn(0x0),
fAODOut(0x0),
fAODExtension(0x0),
+fMcEvent(0x0),
+fMcHandler(0x0),
fJetBranchName(""),
-fJetBranchNameMC(""),
-fJetBranchNameMCFull(""),
+fJetBranchNameChargMC(""),
+fJetBranchNameKine(""),
+fJetBranchNameFullMC(""),
+fJetBranchNameBg(""),
+fJetBranchNameBgChargMC(""),
+fJetBranchNameBgKine(""),
fListJets(0x0),
fListJetsGen(0x0),
fListJetsGenFull(0x0),
+fListJetsBg(0x0),
+fListJetsBgGen(0x0),
fNonStdFile(""),
fSystem(0), //pp=0 pPb=1
fJetParamR(0.4),
+fBgJetParamR(0.3),
+fBgMaxJetPt(8.0),
+fBgConeR(0.4),
fOfflineTrgMask(AliVEvent::kAny),
fMinContribVtx(1),
fVtxZMin(-10.0),
fTriggerEtaCut(0.9),
fTrackEtaCut(0.9),
fTrackLowPtCut(0.15),
+fUseExchContainer(0),
fOutputList(0x0),
fHistEvtSelection(0x0),
fh2Ntriggers(0x0),
fHJetSpec(0x0),
-fHJetDensity(0x0),
-fHJetDensityA4(0x0),
+fHJetSpecSubUeMedian(0x0),
+fHJetSpecSubUeCone(0x0),
+fHJetPhiCorr(0x0),
+fHJetPhiCorrSubUeMedian(0x0),
+fHJetPhiCorrSubUeCone(0x0),
+fHJetUeMedian(0x0),
+fHJetUeCone(0x0),
+fHRhoUeMedianVsCone(0x0),
+//fHJetDensity(0x0),
+//fHJetDensityA4(0x0),
fhJetPhi(0x0),
fhTriggerPhi(0x0),
fhJetEta(0x0),
fhDphiTriggerJetAccept(0x0),
fhCentrality(0x0),
fhCentralityAccept(0x0),
-fHJetPtRaw(0x0),
-fHLeadingJetPtRaw(0x0),
-fHDphiVsJetPtAll(0x0),
+//fHJetPtRaw(0x0),
+//fHLeadingJetPtRaw(0x0),
+//fHDphiVsJetPtAll(0x0),
fhJetPtGenVsJetPtRec(0x0),
-fhJetPtGen(0x0),
+fhJetPtGenVsJetPtRecSubUeMedian(0x0),
+fhJetPtGenVsJetPtRecSubUeCone(0x0),
+fhJetPtGen(0x0),
+fhJetPtSubUeMedianGen(0x0),
+fhJetPtSubUeConeGen(0x0),
fhJetPtGenChargVsJetPtGenFull(0x0),
fhJetPtGenFull(0x0),
fh2NtriggersGen(0x0),
fHJetSpecGen(0x0),
+fHJetSpecSubUeMedianGen(0x0),
+fHJetSpecSubUeConeGen(0x0),
+fHJetPhiCorrGen(0x0),
+fHJetPhiCorrSubUeMedianGen(0x0),
+fHJetPhiCorrSubUeConeGen(0x0),
+fHJetUeMedianGen(0x0),
+fHJetUeConeGen(0x0),
fhPtTrkTruePrimRec(0x0),
fhPtTrkTruePrimGen(0x0),
fhPtTrkSecOrFakeRec(0x0),
-fIsMC(0),
+fHRhoUeMedianVsConeGen(0x0),
+fhEntriesToMedian(0x0),
+fhEntriesToMedianGen(0x0),
+fhCellAreaToMedian(0x0),
+fhCellAreaToMedianGen(0x0),
+fIsChargedMC(0),
+fIsKine(0),
+fIsFullMC(0),
faGenIndex(0),
faRecIndex(0),
fkAcceptance(2.0*TMath::Pi()*1.8),
fEventNumberRangeHigh(99),
fTriggerPtRangeLow(0.0),
fTriggerPtRangeHigh(10000.0),
-fRandom(0x0)
+fFillRespMx(0),
+fRandom(0x0),
+fnTrials(1000),
+fJetFreeAreaFrac(0.5),
+fnPhi(9),
+fnEta(2),
+fEtaSize(0.7),
+fPhiSize(2*TMath::Pi()/fnPhi),
+fCellArea(fPhiSize*fEtaSize),
+fSafetyMargin(1.1),
+fDoubleBinning(kFALSE)
{
// default Constructor
}
fAODIn(0x0),
fAODOut(0x0),
fAODExtension(0x0),
+fMcEvent(0x0),
+fMcHandler(0x0),
fJetBranchName(""),
-fJetBranchNameMC(""),
-fJetBranchNameMCFull(""),
+fJetBranchNameChargMC(""),
+fJetBranchNameKine(""),
+fJetBranchNameFullMC(""),
+fJetBranchNameBg(""),
+fJetBranchNameBgChargMC(""),
+fJetBranchNameBgKine(""),
fListJets(0x0),
fListJetsGen(0x0),
fListJetsGenFull(0x0),
+fListJetsBg(0x0),
+fListJetsBgGen(0x0),
fNonStdFile(""),
fSystem(0), //pp=0 pPb=1
fJetParamR(0.4),
+fBgJetParamR(0.3),
+fBgMaxJetPt(8.0),
+fBgConeR(0.4),
fOfflineTrgMask(AliVEvent::kAny),
fMinContribVtx(1),
fVtxZMin(-10.0),
fTriggerEtaCut(0.9),
fTrackEtaCut(0.9),
fTrackLowPtCut(0.15),
+fUseExchContainer(0),
fOutputList(0x0),
fHistEvtSelection(0x0),
fh2Ntriggers(0x0),
fHJetSpec(0x0),
-fHJetDensity(0x0),
-fHJetDensityA4(0x0),
+fHJetSpecSubUeMedian(0x0),
+fHJetSpecSubUeCone(0x0),
+fHJetPhiCorr(0x0),
+fHJetPhiCorrSubUeMedian(0x0),
+fHJetPhiCorrSubUeCone(0x0),
+fHJetUeMedian(0x0),
+fHJetUeCone(0x0),
+fHRhoUeMedianVsCone(0x0),
+//fHJetDensity(0x0),
+//fHJetDensityA4(0x0),
fhJetPhi(0x0),
fhTriggerPhi(0x0),
fhJetEta(0x0),
fhDphiTriggerJetAccept(0x0),
fhCentrality(0x0),
fhCentralityAccept(0x0),
-fHJetPtRaw(0x0),
-fHLeadingJetPtRaw(0x0),
-fHDphiVsJetPtAll(0x0),
+//fHJetPtRaw(0x0),
+//fHLeadingJetPtRaw(0x0),
+//fHDphiVsJetPtAll(0x0),
fhJetPtGenVsJetPtRec(0x0),
+fhJetPtGenVsJetPtRecSubUeMedian(0x0),
+fhJetPtGenVsJetPtRecSubUeCone(0x0),
fhJetPtGen(0x0),
+fhJetPtSubUeMedianGen(0x0),
+fhJetPtSubUeConeGen(0x0),
fhJetPtGenChargVsJetPtGenFull(0x0),
fhJetPtGenFull(0x0),
fh2NtriggersGen(0x0),
fHJetSpecGen(0x0),
+fHJetSpecSubUeMedianGen(0x0),
+fHJetSpecSubUeConeGen(0x0),
+fHJetPhiCorrGen(0x0),
+fHJetPhiCorrSubUeMedianGen(0x0),
+fHJetPhiCorrSubUeConeGen(0x0),
+fHJetUeMedianGen(0x0),
+fHJetUeConeGen(0x0),
fhPtTrkTruePrimRec(0x0),
fhPtTrkTruePrimGen(0x0),
fhPtTrkSecOrFakeRec(0x0),
-fIsMC(0),
+fHRhoUeMedianVsConeGen(0x0),
+fhEntriesToMedian(0x0),
+fhEntriesToMedianGen(0x0),
+fhCellAreaToMedian(0x0),
+fhCellAreaToMedianGen(0x0),
+fIsChargedMC(0),
+fIsKine(0),
+fIsFullMC(0),
faGenIndex(0),
faRecIndex(0),
fkAcceptance(2.0*TMath::Pi()*1.8),
fEventNumberRangeHigh(99),
fTriggerPtRangeLow(0.0),
fTriggerPtRangeHigh(10000.0),
-fRandom(0x0)
+fFillRespMx(0),
+fRandom(0x0),
+fnTrials(1000),
+fJetFreeAreaFrac(0.5),
+fnPhi(9),
+fnEta(2),
+fEtaSize(0.7),
+fPhiSize(2*TMath::Pi()/fnPhi),
+fCellArea(fPhiSize*fEtaSize),
+fSafetyMargin(1.1),
+fDoubleBinning(kFALSE)
{
// Constructor
-
DefineOutput(1, TList::Class());
+
+ TString dummy(name);
+ if(dummy.Contains("KINE")){
+ DefineInput(1, TClonesArray::Class());
+ DefineInput(2, TClonesArray::Class());
+ }
}
//--------------------------------------------------------------
fAODIn(a.fAODIn),
fAODOut(a.fAODOut),
fAODExtension(a.fAODExtension),
+fMcEvent(a.fMcEvent),
+fMcHandler(a.fMcHandler),
fJetBranchName(a.fJetBranchName),
-fJetBranchNameMC(a.fJetBranchNameMC),
-fJetBranchNameMCFull(a.fJetBranchNameMCFull),
+fJetBranchNameChargMC(a.fJetBranchNameChargMC),
+fJetBranchNameKine(a.fJetBranchNameKine),
+fJetBranchNameFullMC(a.fJetBranchNameFullMC),
+fJetBranchNameBg(a.fJetBranchNameBg),
+fJetBranchNameBgChargMC(a.fJetBranchNameBgChargMC),
+fJetBranchNameBgKine(a.fJetBranchNameBgKine),
fListJets(a.fListJets),
fListJetsGen(a.fListJetsGen),
fListJetsGenFull(a.fListJetsGenFull),
+fListJetsBg(a.fListJetsBg),
+fListJetsBgGen(a.fListJetsBgGen),
fNonStdFile(a.fNonStdFile),
fSystem(a.fSystem),
fJetParamR(a.fJetParamR),
+fBgJetParamR(a.fBgJetParamR),
+fBgMaxJetPt(a.fBgMaxJetPt),
+fBgConeR(a.fBgConeR),
fOfflineTrgMask(a.fOfflineTrgMask),
fMinContribVtx(a.fMinContribVtx),
fVtxZMin(a.fVtxZMin),
fTriggerEtaCut(a.fTriggerEtaCut),
fTrackEtaCut(a.fTrackEtaCut),
fTrackLowPtCut(a.fTrackLowPtCut),
+fUseExchContainer(a.fUseExchContainer),
fOutputList(a.fOutputList),
fHistEvtSelection(a.fHistEvtSelection),
fh2Ntriggers(a.fh2Ntriggers),
fHJetSpec(a.fHJetSpec),
-fHJetDensity(a.fHJetDensity),
-fHJetDensityA4(a.fHJetDensityA4),
+fHJetSpecSubUeMedian(a.fHJetSpecSubUeMedian),
+fHJetSpecSubUeCone(a.fHJetSpecSubUeCone),
+fHJetPhiCorr(a.fHJetPhiCorr),
+fHJetPhiCorrSubUeMedian(a.fHJetPhiCorrSubUeMedian),
+fHJetPhiCorrSubUeCone(a.fHJetPhiCorrSubUeCone),
+fHJetUeMedian(a.fHJetUeMedian),
+fHJetUeCone(a.fHJetUeCone),
+fHRhoUeMedianVsCone(a.fHRhoUeMedianVsCone),
+//fHJetDensity(a.fHJetDensity),
+//fHJetDensityA4(a.fHJetDensityA4),
fhJetPhi(a.fhJetPhi),
fhTriggerPhi(a.fhTriggerPhi),
fhJetEta(a.fhJetEta),
fhDphiTriggerJetAccept(a.fhDphiTriggerJetAccept),
fhCentrality(a.fhCentrality),
fhCentralityAccept(a.fhCentralityAccept),
-fHJetPtRaw(a.fHJetPtRaw),
-fHLeadingJetPtRaw(a.fHLeadingJetPtRaw),
-fHDphiVsJetPtAll(a.fHDphiVsJetPtAll),
+//fHJetPtRaw(a.fHJetPtRaw),
+//fHLeadingJetPtRaw(a.fHLeadingJetPtRaw),
+//fHDphiVsJetPtAll(a.fHDphiVsJetPtAll),
fhJetPtGenVsJetPtRec(a.fhJetPtGenVsJetPtRec),
+fhJetPtGenVsJetPtRecSubUeMedian(a.fhJetPtGenVsJetPtRecSubUeMedian),
+fhJetPtGenVsJetPtRecSubUeCone(a.fhJetPtGenVsJetPtRecSubUeCone),
fhJetPtGen(a.fhJetPtGen),
+fhJetPtSubUeMedianGen(a.fhJetPtSubUeMedianGen),
+fhJetPtSubUeConeGen(a.fhJetPtSubUeConeGen),
fhJetPtGenChargVsJetPtGenFull(a.fhJetPtGenChargVsJetPtGenFull),
fhJetPtGenFull(a.fhJetPtGenFull),
fh2NtriggersGen(a.fh2NtriggersGen),
fHJetSpecGen(a.fHJetSpecGen),
+fHJetSpecSubUeMedianGen(a.fHJetSpecSubUeMedianGen),
+fHJetSpecSubUeConeGen(a.fHJetSpecSubUeConeGen),
+fHJetPhiCorrGen(a.fHJetPhiCorrGen),
+fHJetPhiCorrSubUeMedianGen(a.fHJetPhiCorrSubUeMedianGen),
+fHJetPhiCorrSubUeConeGen(a.fHJetPhiCorrSubUeConeGen),
+fHJetUeMedianGen(a.fHJetUeMedianGen),
+fHJetUeConeGen(a.fHJetUeConeGen),
fhPtTrkTruePrimRec(a.fhPtTrkTruePrimRec),
fhPtTrkTruePrimGen(a.fhPtTrkTruePrimGen),
fhPtTrkSecOrFakeRec(a.fhPtTrkSecOrFakeRec),
-fIsMC(a.fIsMC),
+fHRhoUeMedianVsConeGen(a.fHRhoUeMedianVsConeGen),
+fhEntriesToMedian(a.fhEntriesToMedian),
+fhEntriesToMedianGen(a.fhEntriesToMedianGen),
+fhCellAreaToMedian(a.fhCellAreaToMedian),
+fhCellAreaToMedianGen(a.fhCellAreaToMedianGen),
+fIsChargedMC(a.fIsChargedMC),
+fIsKine(a.fIsKine),
+fIsFullMC(a.fIsFullMC),
faGenIndex(a.faGenIndex),
faRecIndex(a.faRecIndex),
fkAcceptance(a.fkAcceptance),
fEventNumberRangeHigh(a.fEventNumberRangeHigh),
fTriggerPtRangeLow(a.fTriggerPtRangeLow),
fTriggerPtRangeHigh(a.fTriggerPtRangeHigh),
-fRandom(a.fRandom)
+fFillRespMx(a.fFillRespMx),
+fRandom(a.fRandom),
+fnTrials(a.fnTrials),
+fJetFreeAreaFrac(a.fJetFreeAreaFrac),
+fnPhi(a.fnPhi),
+fnEta(a.fnEta),
+fEtaSize(a.fEtaSize),
+fPhiSize(a.fPhiSize),
+fCellArea(a.fCellArea),
+fSafetyMargin(a.fSafetyMargin),
+fDoubleBinning(a.fDoubleBinning)
{
//Copy Constructor
fRandom->SetSeed(0);
delete fListJets;
delete fListJetsGen;
delete fListJetsGenFull;
+ delete fListJetsBg;
+ delete fListJetsBgGen;
delete fOutputList; // ????
delete fRandom;
}
//Implemented Notify() to read the cross sections
//and number of trials from pyxsec.root
//inspired by AliAnalysisTaskJetSpectrum2::Notify()
- if(!fIsMC) return kFALSE;
+ if(!(fIsChargedMC || fIsKine)) return kFALSE;
+ Float_t xsection = 0;
+ Float_t trials = 1;
+ fAvgTrials = 1;
- fESD = dynamic_cast<AliESDEvent*>(InputEvent());
- if(!fESD){
- if(fDebug>1) AliError("ESD not available");
- fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
- }
+ if(fIsChargedMC){
+ fESD = dynamic_cast<AliESDEvent*>(InputEvent());
+ if(!fESD){
+ if(fDebug>1) AliError("ESD not available");
+ fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
+ }
- fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
+
fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
- if(fNonStdFile.Length()!=0){
- // case that we have an AOD extension we can fetch the jets from the extended output
- AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
- fAODExtension = aodH ? aodH->GetExtension(fNonStdFile.Data()) : 0;
- if(!fAODExtension){
- if(fDebug>1) Printf("AODExtension found for %s",fNonStdFile.Data());
- }
- }
+ if(fNonStdFile.Length()!=0){
+ // case that we have an AOD extension we can fetch the jets from the extended output
+
AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
+ fAODExtension = aodH ? aodH->GetExtension(fNonStdFile.Data()) : 0;
+ if(!fAODExtension){
+ if(fDebug>1) Printf("AODExtension found for %s",fNonStdFile.Data());
+ }
+ }
- TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
- Float_t xsection = 0;
- Float_t ftrials = 1;
+ TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
- fAvgTrials = 1;
- if(tree){
- TFile *curfile = tree->GetCurrentFile();
- if(!curfile) {
- Error("Notify","No current file");
- return kFALSE;
- }
- if(!fh1Xsec || !fh1Trials){
- Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
- return kFALSE;
- }
- AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
- fh1Xsec->Fill("<#sigma>",xsection);
- // construct a poor man average trials
- Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
- if(ftrials>=nEntries && nEntries>0.) fAvgTrials = ftrials/nEntries;
- fh1Trials->Fill("#sum{ntrials}",ftrials);
- }
+ if(tree){
+ TFile *curfile = tree->GetCurrentFile();
+ if(!curfile) {
+ Error("Notify","No current file");
+ return kFALSE;
+ }
+ if(!fh1Xsec || !fh1Trials){
+ Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
+ return kFALSE;
+ }
+ AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,trials);
+ fh1Xsec->Fill("<#sigma>",xsection);
+ // construct a poor man average trials
+ Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
+ if(trials>=nEntries && nEntries>0.) fAvgTrials = trials/nEntries;
+ fh1Trials->Fill("#sum{ntrials}",trials);
+ }
+
+ if(fDebug)Printf("Reading File %s",fInputHandler->GetTree()->GetCurrentFile()->GetName());
+ }
- if(fDebug)Printf("Reading File %s",fInputHandler->GetTree()->GetCurrentFile()->GetName());
return kTRUE;
}
void AliAnalysisTaskJetCorePP::Init()
{
// check for jet branches
- if(!strlen(fJetBranchName.Data())){
- AliError("Jet branch name not set.");
+ if(fJetBranchNameKine.Length()==0){
+ if(!strlen(fJetBranchName.Data())){
+ AliError("Jet branch name not set.");
+ }
}
+ fMcHandler = dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
+
+
}
//--------------------------------------------------------------
void AliAnalysisTaskJetCorePP::UserCreateOutputObjects()
{
- // Create histograms
+ // Create histograms and initilize variables
+ fSafetyMargin = fBgConeR*fBgConeR /(fBgJetParamR*fBgJetParamR);
+
// Called once
- fListJets = new TList(); //reconstructed level
+ fListJets = new TList(); //reconstructed level antikt jets
+ fListJetsBg = new TList(); //reconstructed jets to be removed from bg
- fIsMC = (fJetBranchNameMC.Length()>0) ? kTRUE : kFALSE;
+ fIsChargedMC = (fJetBranchNameChargMC.Length()>0) ? kTRUE : kFALSE;
+ fIsKine = (fJetBranchNameKine.Length()>0) ? kTRUE : kFALSE;
+ fIsFullMC = (fJetBranchNameFullMC.Length()>0) ? kTRUE : kFALSE;
fRandom = new TRandom3(0);
- if(fIsMC){
- fListJetsGen = new TList(); //generator level charged jets
- if(fJetBranchNameMCFull.Length()>0)
+ if(fIsChargedMC || fIsKine){ //full MC or pure kine
+ fListJetsGen = new TList(); //generator level charged antikt jets
+ fListJetsBgGen = new TList(); //generator level jets to be removed from bg
+
+ if(fIsFullMC)
fListJetsGenFull = new TList(); //generator level full jets
}
OpenFile(1);
//trigger pt spectrum (reconstructed)
fh2Ntriggers = new TH2F("fh2Ntriggers","# of triggers",
nBinsCentrality,0.0,100.0,50,0.0,50.0);
- fOutputList->Add(fh2Ntriggers);
+ if(!fIsKine) fOutputList->Add(fh2Ntriggers);
+
+ Int_t bw = (fDoubleBinning==0) ? 1 : 2; //make larger bin width
//Centrality, A, pTjet, pTtrigg, dphi
const Int_t dimSpec = 5;
- const Int_t nBinsSpec[dimSpec] = {nBinsCentrality, 100, 200, 50, TMath::Nint(10*(TMath::Pi()-fkDeltaPhiCut))};
- const Double_t lowBinSpec[dimSpec] = {0.0, 0.0, 0, 0.0, fkDeltaPhiCut};
- const Double_t hiBinSpec[dimSpec] = {100.0, 1.0, 200.0,50.0, TMath::Pi()};
+ const Int_t nBinsSpec[dimSpec] = {nBinsCentrality, 50, bw*110, 50, TMath::Nint(10*(TMath::Pi()-fkDeltaPhiCut))};
+ const Double_t lowBinSpec[dimSpec] = {0.0, 0.0, -20, 0.0, fkDeltaPhiCut};
+ const Double_t hiBinSpec[dimSpec] = {100.0, 1.0, 200.0,50.0, TMath::Pi()};
fHJetSpec = new THnSparseF("fHJetSpec",
"Recoil jet spectrum [cent,A,pTjet,pTtrig,dphi]",
dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
- fOutputList->Add(fHJetSpec);
+ if(!fIsKine) fOutputList->Add(fHJetSpec);
+
+ //background estimated as median of kT jets
+ fHJetSpecSubUeMedian = (THnSparseF*) fHJetSpec->Clone("fHJetSpecSubUeMedian");
+ fHJetSpecSubUeMedian->SetTitle("Recoil jet spectrum [cent,A,pTjet-pTUe,pTtrig,dphi]");
+ if(!fIsKine) fOutputList->Add(fHJetSpecSubUeMedian);
+ //background estimated as weighted median of kT jets ala Cone
+ fHJetSpecSubUeCone = (THnSparseF*) fHJetSpec->Clone("fHJetSpecSubUeCone");
+ fHJetSpecSubUeCone->SetTitle("Recoil jet spectrum [cent,A,pTjet-pTUe,pTtrig,dphi]");
+ if(!fIsKine) fOutputList->Add(fHJetSpecSubUeCone);
+
+
+ //A, pTjet, pTtrigg, dphi
+ const Int_t dimCor = 4;
+ const Int_t nBinsCor[dimCor] = {50, 110, 50, 100};
+ const Double_t lowBinCor[dimCor] = {0.0, -20, 0.0, -0.5*TMath::Pi()};
+ const Double_t hiBinCor[dimCor] = {1.0, 200, 50.0, 1.5*TMath::Pi()};
+ fHJetPhiCorr = new THnSparseF("fHJetPhiCorr",
+ "Recoil jet spectrum [A,pTjet,pTtrig,dphi]",
+ dimCor,nBinsCor,lowBinCor,hiBinCor);
+
+ if(!fIsKine) fOutputList->Add(fHJetPhiCorr); // Dphi distribution jet-triger
+
+ //background estimated as median of kT jets
+ fHJetPhiCorrSubUeMedian = (THnSparseF*) fHJetPhiCorr->Clone("fHJetPhiCorrSubUeMedian");
+ fHJetPhiCorrSubUeMedian->SetTitle("Recoil jet spectrum [A,pTjet-pTUe,pTtrig,dphi]");
+ if(!fIsKine) fOutputList->Add(fHJetPhiCorrSubUeMedian);
+ //background estimated as weighted median of kT jets ala Cone
+ fHJetPhiCorrSubUeCone = (THnSparseF*) fHJetPhiCorr->Clone("fHJetPhiCorrSubUeCone");
+ fHJetPhiCorrSubUeCone->SetTitle("Recoil jet spectrum [A,pTjet-pTUe,pTtrig,dphi]");
+ if(!fIsKine) fOutputList->Add(fHJetPhiCorrSubUeCone);
+
//------------------- HISTOS FOR DIAGNOSTIC ----------------------
+ //A, pTjet, pTjet-pTUe, pTUe, rhoUe bg estimate from kT median
+ const Int_t dimSpecMed = 5;
+ const Int_t nBinsSpecMed[dimSpecMed] = {25, 50, 60, 20, 20};
+ const Double_t lowBinSpecMed[dimSpecMed] = {0.0, 0.0, -20.0, 0.0, 0.0};
+ const Double_t hiBinSpecMed[dimSpecMed] = {1.0, 100.0, 100.0, 10.0, 20.0};
+ fHJetUeMedian = new THnSparseF("fHJetUeMedian",
+ "Recoil jet spectrum [A,pTjet,pTjet-pTUe, pTUe, rhoUe]",
+ dimSpecMed, nBinsSpecMed, lowBinSpecMed, hiBinSpecMed);
+ if(!fIsKine) fOutputList->Add(fHJetUeMedian);
+
+ //A, pTjet, pTjet-pTUe, pTUe, rhoUe bg estimate from kT median Cone
+ fHJetUeCone = (THnSparseF*) fHJetUeMedian->Clone("fHJetUeCone");
+ fHJetUeCone->SetTitle("Recoil jet spectrum [A,pTjet,pTjet-pTUe, pTUe, rhoUe]");
+ if(!fIsKine) fOutputList->Add(fHJetUeCone);
+
+ //rho bacground reconstructed data
+ const Int_t dimRho = 2;
+ const Int_t nBinsRho[dimRho] = {50 , 50};
+ const Double_t lowBinRho[dimRho] = {0.0 , 0.0};
+ const Double_t hiBinRho[dimRho] = {20.0 , 20.0};
+
+ fHRhoUeMedianVsCone = new THnSparseF("hRhoUeMedianVsCone","[Rho Cone, Rho Median]",
+ dimRho, nBinsRho, lowBinRho, hiBinRho);
+ if(!fIsKine) fOutputList->Add(fHRhoUeMedianVsCone);
+
//Jet number density histos [Trk Mult, jet density, pT trigger]
- const Int_t dimJetDens = 3;
+ /*const Int_t dimJetDens = 3;
const Int_t nBinsJetDens[dimJetDens] = {100, 100, 10};
const Double_t lowBinJetDens[dimJetDens] = {0.0, 0.0, 0.0};
const Double_t hiBinJetDens[dimJetDens] = {500.0, 5.0, 50.0};
fOutputList->Add(fHJetDensity);
fOutputList->Add(fHJetDensityA4);
-
+ */
//inclusive azimuthal and pseudorapidity histograms
fhJetPhi = new TH2D("fhJetPhi","Azim dist jets vs pTjet",
fhDphiTriggerJetAccept = new TH1D("fhDphiTriggerJetAccept","Deltaphi trig-jet after cut",50, -TMath::Pi(),TMath::Pi());
fhCentrality = new TH1D("fhCentrality","Centrality",20,0,100);
fhCentralityAccept = new TH1D("fhCentralityAccept","Centrality after cut",20,0,100);
-
- fOutputList->Add(fhJetPhi);
- fOutputList->Add(fhTriggerPhi);
- fOutputList->Add(fhJetEta);
- fOutputList->Add(fhTriggerEta);
+ fhEntriesToMedian = new TH1D("fhEntriesToMedian","fhEntriesToMedian",30,0,30);
+ fhCellAreaToMedian = new TH1D("fhCellAreaToMedian", "fhCellAreaToMedian", 75,0,1.5);
+
+ if(!fIsKine){
+ fOutputList->Add(fhJetPhi);
+ fOutputList->Add(fhTriggerPhi);
+ fOutputList->Add(fhJetEta);
+ fOutputList->Add(fhTriggerEta);
+ }
fOutputList->Add(fhVertexZ);
- fOutputList->Add(fhVertexZAccept);
- fOutputList->Add(fhContribVtx);
- fOutputList->Add(fhContribVtxAccept);
- fOutputList->Add(fhDphiTriggerJet);
- fOutputList->Add(fhDphiTriggerJetAccept);
- fOutputList->Add(fhCentrality);
- fOutputList->Add(fhCentralityAccept);
-
+ fOutputList->Add(fhVertexZAccept);
+ if(!fIsKine){
+ fOutputList->Add(fhContribVtx);
+ fOutputList->Add(fhContribVtxAccept);
+ fOutputList->Add(fhDphiTriggerJet);
+ fOutputList->Add(fhDphiTriggerJetAccept);
+ fOutputList->Add(fhCentrality);
+ fOutputList->Add(fhCentralityAccept);
+ fOutputList->Add(fhEntriesToMedian);
+ fOutputList->Add(fhCellAreaToMedian);
+ }
// raw spectra of INCLUSIVE jets
//Centrality, pTjet, A
- const Int_t dimRaw = 3;
+ /*const Int_t dimRaw = 3;
const Int_t nBinsRaw[dimRaw] = {nBinsCentrality, 50, 100};
const Double_t lowBinRaw[dimRaw] = {0.0, 0.0, 0.0};
const Double_t hiBinRaw[dimRaw] = {100.0, 100, 1.0};
"Dphi vs jet pT [cent,Dphi,pTjet,pTtrigg]",
dimDp,nBinsDp,lowBinDp,hiBinDp);
fOutputList->Add(fHDphiVsJetPtAll);
-
+ */
//analyze MC generator level
- if(fIsMC){
- fhJetPtGenVsJetPtRec = new TH2D("fhJetPtGenVsJetPtRec","JetPtGenVsJetPtRec", 200,0,200, 200,0,200);
- fOutputList->Add(fhJetPtGenVsJetPtRec); //gen MC jet pt spectrum versus reconstructed pt spectrum
-
- fhJetPtGen = new TH1D("fhJetPtGen","Jet Pt (MC charged jets Gen)",200,0,200); //MC generator jet pt spectrum
- fOutputList->Add(fhJetPtGen);
-
- fhJetPtGenChargVsJetPtGenFull = new TH2D("fhJetPtGenChargVsJetPtGenFull","fhJetPtGenChargVsJetPtGenFull", 200,0,200, 200,0,200);
- fOutputList->Add(fhJetPtGenChargVsJetPtGenFull); //gen MC jet pt spectrum versus reconstructed pt spectrum
-
- fhJetPtGenFull = new TH1D("fhJetPtGenFull","Jet Pt (MC Full jets Gen)",200,0,200); //MC generator jet pt spectrum
- fOutputList->Add(fhJetPtGenFull);
-
+ if(fIsChargedMC || fIsKine){
+ if(fFillRespMx){
+ //Fill response matrix only once
+ fhJetPtGenVsJetPtRec = new TH2D("fhJetPtGenVsJetPtRec","JetPtGenVsJetPtRec", bw*100,0,200, bw*100,0,200);
+ fOutputList->Add(fhJetPtGenVsJetPtRec); //gen MC charg jet pt spectrum versus rec charged jet pt spectrum
+ //....
+ fhJetPtGenVsJetPtRecSubUeMedian = new TH2D("fhJetPtGenVsJetPtRecSubUeMedian","fhJetPtGenVsJetPtRecSubUeMedian", bw*110,-20,200, bw*110,-20,200);
+ fOutputList->Add(fhJetPtGenVsJetPtRecSubUeMedian); // with kT median bg subtr
+ //....
+ fhJetPtGenVsJetPtRecSubUeCone=(TH2D*)fhJetPtGenVsJetPtRecSubUeMedian->Clone("fhJetPtGenVsJetPtRecSubUeCone");
+ fhJetPtGenVsJetPtRecSubUeCone->SetTitle("fhJetPtGenVsJetPtRecSubUeCone");
+ fOutputList->Add(fhJetPtGenVsJetPtRecSubUeCone); // with weighted kT median bg subtr
+ //....
+ fhJetPtGen = new TH1D("fhJetPtGen","Jet Pt (MC Gen)",bw*100,0,200); //MC generator charged jet pt spectrum
+ fOutputList->Add(fhJetPtGen);
+ //....
+ fhJetPtSubUeMedianGen = new TH1D("fhJetPtSubUeMedianGen","Jet Pt - UE pT (MC Gen)",bw*110,-20,200);
+ fOutputList->Add(fhJetPtSubUeMedianGen); // with kT median bg subtr
+ //....
+ fhJetPtSubUeConeGen = (TH1D*) fhJetPtSubUeMedianGen->Clone("fhJetPtSubUeConeGen");
+ fOutputList->Add(fhJetPtSubUeConeGen); // with weighted kT median bg subtr
+
+ //....
+ if(fIsFullMC){
+ fhJetPtGenChargVsJetPtGenFull = new TH2D("fhJetPtGenChargVsJetPtGenFull","fhJetPtGenChargVsJetPtGenFull", 100,0,200, 100,0,200);
+ fOutputList->Add(fhJetPtGenChargVsJetPtGenFull); //gen full MC jet pt versus gen charged jet MC pt
+ //....
+ fhJetPtGenFull = new TH1D("fhJetPtGenFull","Jet Pt (MC Full jets Gen)",100,0,200); //MC generator full jet pt spectrum
+ fOutputList->Add(fhJetPtGenFull);
+ }
+ }
+ //....
fh2NtriggersGen = (TH2F*) fh2Ntriggers->Clone("fh2NtriggersGen");
fh2NtriggersGen->SetTitle(Form("%s Gen MC",fh2Ntriggers->GetTitle()));
fOutputList->Add(fh2NtriggersGen);
fHJetSpecGen->SetTitle(Form("%s Gen MC",fHJetSpec->GetTitle()));
fOutputList->Add(fHJetSpecGen);
- //track efficiency/contamination histograms eta versus pT
- fhPtTrkTruePrimRec = new TH2D("fhPtTrkTruePrimRec","PtTrkTruePrimRec",100,0,20,18,-0.9,0.9);
- fOutputList->Add(fhPtTrkTruePrimRec);
+ fHJetSpecSubUeMedianGen = (THnSparseF*) fHJetSpecSubUeMedian->Clone("fHJetSpecSubUeMedianGen");
+ fHJetSpecSubUeMedianGen->SetTitle(Form("%s Gen MC",fHJetSpecSubUeMedian->GetTitle()));
+ fOutputList->Add(fHJetSpecSubUeMedianGen);
+
+ fHJetSpecSubUeConeGen = (THnSparseF*) fHJetSpecSubUeCone->Clone("fHJetSpecSubUeConeGen");
+ fHJetSpecSubUeConeGen->SetTitle(Form("%s Gen MC",fHJetSpecSubUeCone->GetTitle()));
+ fOutputList->Add(fHJetSpecSubUeConeGen);
+ //---
+ fHJetPhiCorrGen = (THnSparseF*) fHJetPhiCorr->Clone("fHJetPhiCorrGen");
+ fHJetPhiCorrGen->SetTitle(Form("%s Gen MC",fHJetPhiCorr->GetTitle()));
+ fOutputList->Add(fHJetPhiCorrGen);
+
+ fHJetPhiCorrSubUeMedianGen = (THnSparseF*) fHJetPhiCorrSubUeMedian->Clone("fHJetPhiCorrSubUeMedianGen");
+ fHJetPhiCorrSubUeMedianGen->SetTitle(Form("%s Gen MC",fHJetPhiCorrSubUeMedian->GetTitle()));
+ fOutputList->Add(fHJetPhiCorrSubUeMedianGen);
+
+ fHJetPhiCorrSubUeConeGen = (THnSparseF*) fHJetPhiCorrSubUeCone->Clone("fHJetPhiCorrSubUeConeGen");
+ fHJetPhiCorrSubUeConeGen->SetTitle(Form("%s Gen MC", fHJetPhiCorrSubUeCone->GetTitle()));
+ fOutputList->Add(fHJetPhiCorrSubUeConeGen);
+
+ //---
+ fHJetUeMedianGen = (THnSparseF*) fHJetUeMedian->Clone("fHJetUeMedianGen");
+ fHJetUeMedianGen->SetTitle(Form("%s Gen MC", fHJetUeMedian->GetTitle()));
+ fOutputList->Add(fHJetUeMedianGen);
+
+ fHJetUeConeGen = (THnSparseF*) fHJetUeCone->Clone("fHJetUeConeGen");
+ fHJetUeConeGen->SetTitle(Form("%s Gen MC", fHJetUeCone->GetTitle()));
+ fOutputList->Add(fHJetUeConeGen);
+
+ if(fFillRespMx){
+ //track efficiency/contamination histograms eta versus pT
+ fhPtTrkTruePrimRec = new TH2D("fhPtTrkTruePrimRec","PtTrkTruePrimRec",100,0,20,18,-0.9,0.9);
+ fOutputList->Add(fhPtTrkTruePrimRec);
- fhPtTrkTruePrimGen = (TH2D*) fhPtTrkTruePrimRec->Clone("fhPtTrkTruePrimGen");
- fhPtTrkTruePrimGen->SetTitle("PtTrkTruePrimGen");
- fOutputList->Add(fhPtTrkTruePrimGen);
+ fhPtTrkTruePrimGen = (TH2D*) fhPtTrkTruePrimRec->Clone("fhPtTrkTruePrimGen");
+ fhPtTrkTruePrimGen->SetTitle("PtTrkTruePrimGen");
+ fOutputList->Add(fhPtTrkTruePrimGen);
- fhPtTrkSecOrFakeRec = (TH2D*) fhPtTrkTruePrimRec->Clone("fhPtTrkSecOrFakeRec");
- fhPtTrkSecOrFakeRec->SetTitle("PtTrkSecOrFakeRec");
- fOutputList->Add(fhPtTrkSecOrFakeRec);
+ fhPtTrkSecOrFakeRec = (TH2D*) fhPtTrkTruePrimRec->Clone("fhPtTrkSecOrFakeRec");
+ fhPtTrkSecOrFakeRec->SetTitle("PtTrkSecOrFakeRec");
+ fOutputList->Add(fhPtTrkSecOrFakeRec);
+ }
+
+ fHRhoUeMedianVsConeGen = (THnSparseF*) fHRhoUeMedianVsCone->Clone("hRhoUeMedianVsConeGen");
+ fHRhoUeMedianVsConeGen->SetTitle(Form("%s Gen MC", fHRhoUeMedianVsCone->GetTitle()));
+ fOutputList->Add(fHRhoUeMedianVsConeGen);
+
+ fhEntriesToMedianGen = (TH1D*) fhEntriesToMedian->Clone("fhEntriesToMedianGen");
+ fhEntriesToMedianGen->SetTitle(Form("%s Gen MC", fhEntriesToMedian->GetTitle()));
+ fOutputList->Add(fhEntriesToMedianGen);
+
+ fhCellAreaToMedianGen = (TH1D*) fhCellAreaToMedian->Clone("fhCellAreaToMedianGen");
+ fhCellAreaToMedianGen->SetTitle(Form("%s Gen MC", fhCellAreaToMedian->GetTitle()));
+ fOutputList->Add(fhCellAreaToMedianGen);
}
//-------------------------------------
// pythia histograms
void AliAnalysisTaskJetCorePP::UserExec(Option_t *)
{
+ //User Exec
+
+
//Event loop
Double_t eventW = 1.0;
Double_t ptHard = 0.0;
Double_t nTrials = 1.0; // Trials for MC trigger
- if(fIsMC) fh1AvgTrials->Fill("#sum{avg ntrials}",fAvgTrials);
+ if(fIsChargedMC || fIsKine) fh1AvgTrials->Fill("#sum{avg ntrials}",fAvgTrials);
if(TMath::Abs((Float_t) fJetParamR)<0.00001){
- AliError("Cone radius is set to zero.");
+ AliError("ANTIKT Cone radius is set to zero.");
return;
}
- if(!strlen(fJetBranchName.Data())){
- AliError("Jet branch name not set.");
+
+ if(TMath::Abs((Float_t) fBgJetParamR)<0.00001){
+ AliError("ANTIKT Cone radius is set to zero.");
return;
}
+ if(!fIsKine){ //real data or full MC
+ if(!strlen(fJetBranchName.Data())){
+ AliError("Name of jet branch not set.");
+ return;
+ }
+ if(!strlen(fJetBranchNameBg.Data())){
+ AliError("Name of jet bg branch not set.");
+ return;
+ }
+ }else{ //Kine
+ if(!strlen(fJetBranchNameBgKine.Data())){
+ AliError("Name of jet bg branch for kine not set.");
+ return;
+ }
+
+ Init();
+ if(fMcHandler){
+ fMcEvent = fMcHandler->MCEvent();
+ }else{
+ if(fDebug > 1) printf("AliAnalysisTaskJetCorePP::Exec() fMcHandler=NULL\n");
+ PostData(1, fOutputList);
+ return;
+ }
+ if(!fMcEvent){
+ if(fDebug > 1) printf("AliAnalysisTaskJetCorePP::Exec() fMcEvent=NULL \n");
+ PostData(1, fOutputList);
+ return;
+ }
+
+ Float_t xsection = 0;
+ Float_t trials = 0;
+
+ AliGenPythiaEventHeader *genPH =
+ dynamic_cast<AliGenPythiaEventHeader*> (fMcEvent->GenEventHeader());
+ if(genPH){
+ xsection = genPH->GetXsection();
+ trials = genPH->Trials();
+ ptHard = genPH->GetPtHard();
+ }
+ fh1Xsec->Fill("<#sigma>",xsection);
+ fh1Trials->Fill("#sum{ntrials}",trials);
+ fh1PtHard->Fill(ptHard,eventW);
+ fh1PtHardNoW->Fill(ptHard,1);
+ fh1PtHardTrials->Fill(ptHard,trials);
+ }
+
+
fESD = dynamic_cast<AliESDEvent*>(InputEvent());
if(!fESD){
if(fDebug>1) AliError("ESD not available");
fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
AliAODEvent* aod = NULL;
// take all other information from the aod we take the tracks from
- if(!aod){
- if(!fESD) aod = fAODIn;
- else aod = fAODOut;
+ if(!aod && !fIsKine){
+ if(!fESD) aod = fAODIn; //not ESD and not kine => input AOD
+ else aod = fAODOut;// ESD or kine
}
+
if(fNonStdFile.Length()!=0){
// case that we have an AOD extension we can fetch the jets from the extended output
AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
// ----------------- EVENT SELECTION --------------------------
fHistEvtSelection->Fill(1); // number of events before event selection
- // physics selection
- AliInputEventHandler* inputHandler = (AliInputEventHandler*)
+ if(!fIsKine){
+ // physics selection
+ AliInputEventHandler* inputHandler = (AliInputEventHandler*)
((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
- if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
- if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
- fHistEvtSelection->Fill(2);
- PostData(1, fOutputList);
- return;
- }
-
+ if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
+ if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
+ fHistEvtSelection->Fill(2);
+ PostData(1, fOutputList);
+ return;
+ }
+
//check AOD pointer
- if(!aod){
- if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
- fHistEvtSelection->Fill(3);
- PostData(1, fOutputList);
- return;
- }
-
- // vertex selection
- AliAODVertex* primVtx = aod->GetPrimaryVertex();
+ if(!aod){
+ if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
+ fHistEvtSelection->Fill(3);
+ PostData(1, fOutputList);
+ return;
+ }
- if(!primVtx){
- if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
- fHistEvtSelection->Fill(3);
- PostData(1, fOutputList);
- return;
- }
+ // vertex selection for reconstructed data
+ AliAODVertex* primVtx = aod->GetPrimaryVertex();
- Int_t nTracksPrim = primVtx->GetNContributors();
- Float_t vtxz = primVtx->GetZ();
- //Input events
- fhContribVtx->Fill(nTracksPrim);
- if( nTracksPrim > 0 ) fhVertexZ->Fill(vtxz);
+ if(!primVtx){
+ if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
+ fHistEvtSelection->Fill(3);
+ PostData(1, fOutputList);
+ return;
+ }
+
+ Int_t nTracksPrim = primVtx->GetNContributors();
+ Float_t vtxz = primVtx->GetZ();
+ //Input events
+ fhContribVtx->Fill(nTracksPrim);
+ if( nTracksPrim > 0 ) fhVertexZ->Fill(vtxz);
- if((nTracksPrim < fMinContribVtx) ||
- (vtxz < fVtxZMin) ||
- (vtxz > fVtxZMax)){
- if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",
+ if((nTracksPrim < fMinContribVtx) ||
+ (vtxz < fVtxZMin) ||
+ (vtxz > fVtxZMax)){
+ if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",
(char*)__FILE__,__LINE__,vtxz);
- fHistEvtSelection->Fill(3);
- PostData(1, fOutputList);
- return;
- }else{
+ fHistEvtSelection->Fill(3);
+ PostData(1, fOutputList);
+ return;
+ }else{
//Accepted events
- fhContribVtxAccept->Fill(nTracksPrim);
- fhVertexZAccept->Fill(vtxz);
+ fhContribVtxAccept->Fill(nTracksPrim);
+ fhVertexZAccept->Fill(vtxz);
+ }
+ }else{ //KINE cut on vertex
+ const AliVVertex *vtxMC = fMcEvent->GetPrimaryVertex();
+ Float_t zVtx = vtxMC->GetZ();
+ fhVertexZ->Fill(zVtx);
+ if(zVtx < fVtxZMin || zVtx > fVtxZMax){ //vertex cut
+ fHistEvtSelection->Fill(3);
+ PostData(1, fOutputList);
+ return;
+ }else{
+ fhVertexZAccept->Fill(zVtx);
+ }
}
//FK// No event class selection imposed
// event class selection (from jet helper task)
}
//-----------------select disjunct event subsamples ----------------
- Int_t eventnum = aod->GetHeader()->GetEventNumberESDFile();
- Int_t lastdigit = eventnum % 10;
- if(!(fEventNumberRangeLow<=lastdigit && lastdigit<=fEventNumberRangeHigh)){
- fHistEvtSelection->Fill(5);
- PostData(1, fOutputList);
- return;
- }
+ if(!fIsKine){ //reconstructed data
+ Int_t eventnum = aod->GetHeader()->GetEventNumberESDFile();
+ Int_t lastdigit = eventnum % 10;
+ if(!(fEventNumberRangeLow<=lastdigit && lastdigit<=fEventNumberRangeHigh)){
+ fHistEvtSelection->Fill(5);
+ PostData(1, fOutputList);
+ return;
+ }
+ }
if(fDebug) std::cout<<" ACCEPTED EVENT "<<endl;
fHistEvtSelection->Fill(0); // accepted events
- // ------------------- end event selection --------------------
-
+ // ==================== end event selection ============================
+
+ Double_t tmpArrayFive[5];
+ Double_t tmpArrayFour[4];
- // fetch RECONSTRUCTED jets
- TClonesArray *aodJets = 0x0;
-
- if(fAODOut && !aodJets){
- aodJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName.Data()));
- }
- if(fAODExtension && !aodJets){
- aodJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName.Data()));
- }
- if(fAODIn && !aodJets){
- aodJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName.Data()));
- }
- //--------- Fill list of RECONSTRUCTED jets --------------
- fListJets->Clear();
- if(aodJets){
- if(fDebug) Printf("########## %s: %d jets",fJetBranchName.Data(), aodJets->GetEntriesFast());
- for(Int_t iJet = 0; iJet < aodJets->GetEntriesFast(); iJet++) {
- AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets)[iJet]);
- if (jet) fListJets->Add(jet);
- }
- }
+ TList particleList; //list of tracks
+ Int_t indexTrigg = -1;
+ Double_t rhoFromCellMedian=0.0, rhoCone=0.0;
+
+ if(!fIsKine){
+ //=============== Reconstructed antikt jets ===============
+ ReadTClonesArray(fJetBranchName.Data() , fListJets);
+ ReadTClonesArray(fJetBranchNameBg.Data(), fListJetsBg);
+
+ //============ Estimate background in reconstructed events ===========
+
+ //Find Hadron trigger and Estimate rho from cone
+ indexTrigg = GetListOfTracks(&particleList); //index of trigger hadron in Particle list
+ EstimateBgCone(fListJets, &particleList, rhoCone);
- //--------- Fill list of MC GENERATOR LEVEL jets --------------
+ //Estimate rho from cell median minus jets
+ EstimateBgRhoMedian(fListJetsBg, &particleList, rhoFromCellMedian,0); //real data
+ }
+ //============== analyze generator level MC ================
TList particleListGen; //list of tracks in MC
- if(fIsMC){ //analyze generator level MC
- // fetch MC generator level jets
- TClonesArray *aodGenJets = NULL;
- TClonesArray *aodGenJetsFull = NULL;
-
- //charged jets
- if(fAODOut&&!aodGenJets){
- aodGenJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchNameMC.Data()));
- }
- if(fAODExtension&&!aodGenJets){
- aodGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchNameMC.Data()));
- }
- if(fAODIn&&!aodGenJets){
- aodGenJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchNameMC.Data()));
- }
- if(!aodGenJets){
- Printf("%s:%d no generated Jet array with name %s in AOD",
- (char*)__FILE__,__LINE__, fJetBranchNameMC.Data());
- PostData(1, fOutputList);
- return;
- }
-
- fListJetsGen->Clear();
+ if(fIsChargedMC || fIsKine){
- //full jets
- if(fJetBranchNameMCFull.Length()>0){
-
- if(fAODOut&&!aodGenJetsFull){
- aodGenJetsFull = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchNameMCFull.Data()));
- }
- if(fAODExtension&&!aodGenJetsFull){
- aodGenJetsFull = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchNameMCFull.Data()));
- }
- if(fAODIn&&!aodGenJetsFull){
- aodGenJetsFull = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchNameMCFull.Data()));
- }
+ if(fIsKine){ //pure kine
+
+ //================= generated charged antikt jets ================
+ ReadTClonesArray(fJetBranchNameKine.Data(), fListJetsGen);
+ ReadTClonesArray(fJetBranchNameBgKine.Data(), fListJetsBgGen);
+ }else{
+ //================= generated charged antikt jets ================
+ ReadTClonesArray(fJetBranchNameChargMC.Data(), fListJetsGen);
+ ReadTClonesArray(fJetBranchNameBgChargMC.Data(), fListJetsBgGen);
- if(!aodGenJetsFull){
- Printf("%s:%d no generated Jet array with name %s in AOD",
- (char*)__FILE__,__LINE__, fJetBranchNameMCFull.Data());
- PostData(1, fOutputList);
- return;
+ if(fIsFullMC){ //generated full jets
+ ReadTClonesArray(fJetBranchNameFullMC.Data(), fListJetsGenFull);
}
-
- fListJetsGenFull->Clear();
}
-
-
+ //========================================================
//serarch for charged trigger at the MC generator level
+
Int_t indexTriggGen = -1;
Double_t ptTriggGen = -1;
- Int_t iCounterGen = 0;
+ Int_t iCounterGen = 0; //number of entries in particleListGen array
Int_t triggersMC[200];//list of trigger candidates
Int_t ntriggersMC = 0; //index in triggers array
- if(fESD){//ESD input
-
- AliMCEvent* mcEvent = MCEvent();
- if(!mcEvent){
- PostData(1, fOutputList);
- return;
- }
-
- AliGenPythiaEventHeader *pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
- if(pythiaGenHeader){
- nTrials = pythiaGenHeader->Trials();
- ptHard = pythiaGenHeader->GetPtHard();
- fh1PtHard->Fill(ptHard,eventW);
- fh1PtHardNoW->Fill(ptHard,1);
- fh1PtHardTrials->Fill(ptHard,nTrials);
- }
+ if(!fIsKine){
+ if(fESD){//ESD input
- for(Int_t it = 0; it < mcEvent->GetNumberOfTracks(); it++){
- if(!mcEvent->IsPhysicalPrimary(it)) continue;
- AliMCParticle* part = (AliMCParticle*) mcEvent->GetTrack(it);
- if(!part) continue;
- if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){
+ AliMCEvent* mcEvent = MCEvent();
+ if(!mcEvent){
+ PostData(1, fOutputList);
+ return;
+ }
- if(fHardest==0 && ntriggersMC<200){//single inclusive trigger
- if(indexTriggGen > -1){//trigger candidater was found
- triggersMC[ntriggersMC] = indexTriggGen;
- ntriggersMC++;
+ AliGenPythiaEventHeader *pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
+ if(pythiaGenHeader){
+ nTrials = pythiaGenHeader->Trials();
+ ptHard = pythiaGenHeader->GetPtHard();
+ fh1PtHard->Fill(ptHard,eventW);
+ fh1PtHardNoW->Fill(ptHard,1);
+ fh1PtHardTrials->Fill(ptHard,nTrials);
+ }
+
+ for(Int_t it = 0; it < mcEvent->GetNumberOfTracks(); it++){
+ if(!mcEvent->IsPhysicalPrimary(it)) continue;
+ AliMCParticle* part = (AliMCParticle*) mcEvent->GetTrack(it);
+ if(!part) continue;
+ if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){
+
+ if(fHardest==0 && ntriggersMC<200){//single inclusive trigger
+ if(indexTriggGen > -1){//trigger candidate was found
+ triggersMC[ntriggersMC] = indexTriggGen;
+ ntriggersMC++;
+ }
}
+
+ iCounterGen++;//index in particleListGen array
+ }
+ }
+ }else{ //AOD input
+ if(!fAODIn){
+ PostData(1, fOutputList);
+ return;
+ }
+ TClonesArray *tca = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(AliAODMCParticle::StdBranchName()));
+ if(!tca){
+ PostData(1, fOutputList);
+ return;
+ }
+ for(Int_t it = 0; it < tca->GetEntriesFast(); it++){
+ AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
+ if(!part) continue;
+ if(!part->IsPhysicalPrimary()) continue;
+ if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){
+
+ if(fHardest==0 && ntriggersMC<200){//single inclusive trigger
+ if(indexTriggGen > -1){ //trigger candidater was found
+ triggersMC[ntriggersMC] = indexTriggGen;
+ ntriggersMC++;
+ }
+ }
+
+ iCounterGen++;//number of entries in particleListGen array
}
-
- iCounterGen++;
}
}
- }else{ //AOD input
- if(!fAODIn){
- PostData(1, fOutputList);
- return;
- }
- TClonesArray *tca = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(AliAODMCParticle::StdBranchName()));
- if(!tca){
- PostData(1, fOutputList);
- return;
- }
- for(Int_t it = 0; it < tca->GetEntriesFast(); it++){
- AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
- if(!part) continue;
- if(!part->IsPhysicalPrimary()) continue;
- if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){
+ }else{ //analyze kine
+ for(Int_t it = 0; it < fMcEvent->GetNumberOfTracks(); it++){
+ if(!fMcEvent->IsPhysicalPrimary(it)) continue;
+ AliMCParticle* part = (AliMCParticle*) fMcEvent->GetTrack(it);
+ if(!part) continue;
+ if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){
+
if(fHardest==0 && ntriggersMC<200){//single inclusive trigger
- if(indexTriggGen > -1){ //trigger candidater was found
+ if(indexTriggGen > -1){//trigger candidate was found
triggersMC[ntriggersMC] = indexTriggGen;
ntriggersMC++;
}
}
-
- iCounterGen++;
+
+ iCounterGen++;//index in particleListGen array
}
}
}
+
+ //============== Estimate bg in generated events ==============
+ Double_t rhoFromCellMedianGen=0.0, rhoConeGen=0.0;
+ //Estimate rho from cone
+ EstimateBgCone(fListJetsGen, &particleListGen, rhoConeGen);
+
+ //Estimate rho from cell median minus jets
+ EstimateBgRhoMedian(fListJetsBgGen, &particleListGen, rhoFromCellMedianGen,1);//mc data
+
+ //============ Generator trigger+jet ==================
if(fHardest==0){ //single inclusive trigger
if(ntriggersMC>0){ //there is at least one trigger
Int_t rnd = fRandom->Integer(ntriggersMC); //0 to ntriggers-1
}
}
- //================== Fill jet list ===================
- if(aodGenJets){ //charged jets
- if(fDebug) Printf("########## %s: %d jets",fJetBranchNameMC.Data(), aodGenJets->GetEntriesFast());
-
- for(Int_t igJet = 0; igJet < aodGenJets->GetEntriesFast(); igJet++) {
- AliAODJet *jetGen = dynamic_cast<AliAODJet*>((*aodGenJets)[igJet]);
- if(jetGen) fListJetsGen->Add(jetGen);
- }
- }
-
- if(aodGenJetsFull){ //full jets
- if(fDebug) Printf("########## %s: %d jets",fJetBranchNameMCFull.Data(), aodGenJetsFull->GetEntriesFast());
-
- for(Int_t igJet = 0; igJet < aodGenJetsFull->GetEntriesFast(); igJet++) {
- AliAODJet *jetGen = dynamic_cast<AliAODJet*>((*aodGenJetsFull)[igJet]);
- if(jetGen) fListJetsGenFull->Add(jetGen);
- }
- }
-
-
- //============ Generator trigger+jet ==================
+ //-----------
Int_t ilowGen = (fHardest==0 || fHardest==1) ? indexTriggGen : 0;
Int_t ihighGen = (fHardest==0 || fHardest==1) ? indexTriggGen+1 : particleListGen.GetEntries();
-
+ Bool_t fillOnceGen = kTRUE;
+ //-----------
+
for(Int_t igen= ilowGen; igen< ihighGen; igen++){ //loop over possible trigger
indexTriggGen = igen; //trigger hadron
fh2NtriggersGen->Fill(centValue, ptTriggGen);
//Count jets and trigger-jet pairs at MC generator level
- if(!aodGenJets) continue;
for(Int_t ij=0; ij<fListJetsGen->GetEntries(); ij++){
AliAODJet* jet = (AliAODJet*)(fListJetsGen->At(ij));
if(!jet) continue;
Double_t etaJetGen = jet->Eta();
+ Double_t areaJetGen = jet->EffectiveAreaCharged();
if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){
+ Double_t ptUeFromCellMedianGen = rhoFromCellMedianGen*areaJetGen;
+ Double_t ptUeConeGen = rhoConeGen*areaJetGen;
+
+ //Rongrong's analysis
+ Double_t dPhiGen = jet->Phi() - triggerGen->Phi();
+ if(dPhiGen>2*TMath::Pi()) dPhiGen -= 2*TMath::Pi();
+ if(dPhiGen<-2*TMath::Pi()) dPhiGen += 2*TMath::Pi();
+ if(dPhiGen<-0.5*TMath::Pi()) dPhiGen += 2*TMath::Pi();
+ if(dPhiGen>1.5*TMath::Pi()) dPhiGen -= 2*TMath::Pi();
+
+ //[A,pTjet,pTtrig,dphi]
+ tmpArrayFour[0] = areaJetGen;
+ tmpArrayFour[1] = jet->Pt();
+ tmpArrayFour[2] = ptTriggGen;
+ tmpArrayFour[3] = dPhiGen;
+
+ fHJetPhiCorrGen->Fill(tmpArrayFour); // Dphi distribution jet-triger
+
+ //[A,pTjet-pTUe,pTtrig,dphi]
+ tmpArrayFour[1] = jet->Pt() - ptUeFromCellMedianGen;
+ fHJetPhiCorrSubUeMedianGen->Fill(tmpArrayFour);
+
+ //[A,pTjet-pTUe,pTtrig,dphi]
+ tmpArrayFour[1] = jet->Pt() - ptUeConeGen;
+ fHJetPhiCorrSubUeConeGen->Fill(tmpArrayFour);
+
+ //Leticia's analysis
Double_t dphi = RelativePhi(triggerGen->Phi(), jet->Phi());
if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue;
//Centrality, A, pT, pTtrigg
- Double_t fillspecgen[] = { centValue,
- jet->EffectiveAreaCharged(),
- jet->Pt(),
- ptTriggGen,
- TMath::Abs((Double_t) dphi)
- };
- fHJetSpecGen->Fill(fillspecgen);
+ tmpArrayFive[0] = centValue;
+ tmpArrayFive[1] = areaJetGen;
+ tmpArrayFive[2] = jet->Pt();
+ tmpArrayFive[3] = ptTriggGen;
+ tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
+ fHJetSpecGen->Fill(tmpArrayFive);
+
+ //Centrality, A, pTjet-pTbgCellMedian, pTtrigg, dphi
+ tmpArrayFive[0] = centValue;
+ tmpArrayFive[1] = areaJetGen;
+ tmpArrayFive[2] = jet->Pt() - ptUeFromCellMedianGen;
+ tmpArrayFive[3] = ptTriggGen;
+ tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
+ fHJetSpecSubUeMedianGen->Fill(tmpArrayFive);
+
+ //Centrality, A, pTjet-pTbgCone, pTtrigg, dphi
+ tmpArrayFive[0] = centValue;
+ tmpArrayFive[1] = areaJetGen;
+ tmpArrayFive[2] = jet->Pt() - ptUeConeGen;
+ tmpArrayFive[3] = ptTriggGen;
+ tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
+ fHJetSpecSubUeConeGen->Fill(tmpArrayFive);
+
+ //Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", kT median
+ tmpArrayFive[0] = areaJetGen;
+ tmpArrayFive[1] = jet->Pt();
+ tmpArrayFive[2] = jet->Pt() - ptUeFromCellMedianGen;
+ tmpArrayFive[3] = ptUeFromCellMedianGen;
+ tmpArrayFive[4] = rhoFromCellMedianGen;
+ fHJetUeMedianGen->Fill(tmpArrayFive);
+
+ //Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", perpendicular Cone
+ tmpArrayFive[0] = areaJetGen;
+ tmpArrayFive[1] = jet->Pt();
+ tmpArrayFive[2] = jet->Pt() - ptUeConeGen;
+ tmpArrayFive[3] = ptUeConeGen;
+ tmpArrayFive[4] = rhoConeGen;
+ fHJetUeConeGen->Fill(tmpArrayFive);
+
+ if(fillOnceGen){
+ Double_t fillRhoGen[] = { rhoConeGen, rhoFromCellMedianGen};
+ fHRhoUeMedianVsConeGen->Fill(fillRhoGen);
+ fillOnceGen = kFALSE;
+ }
}//back to back jet-trigger pair
}//jet loop
}//trigger loop
-
- //================ RESPONSE MATRIX ===============
- if(aodGenJets){
+ if(fFillRespMx && !fIsKine){
+ //================ RESPONSE MATRIX ===============
//Count jets and trigger-jet pairs at MC generator level
for(Int_t ij=0; ij<fListJetsGen->GetEntries(); ij++){
AliAODJet* jet = (AliAODJet*)(fListJetsGen->At(ij));
if(!jet) continue;
Double_t etaJetGen = jet->Eta();
Double_t ptJetGen = jet->Pt();
-
+
if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){
- fhJetPtGen->Fill(ptJetGen); // generator level pt spectum of jets response mx normalization
+ fhJetPtGen->Fill(ptJetGen); // gen level pt spectum of jets response mx normalization
+
+ Double_t areaJetGen = jet->EffectiveAreaCharged();
+ Double_t ptUeFromCellMedianGen = rhoFromCellMedianGen*areaJetGen;
+ Double_t ptUeConeGen = rhoConeGen*areaJetGen;
+
+ fhJetPtSubUeMedianGen->Fill(ptJetGen - ptUeFromCellMedianGen);
+ fhJetPtSubUeConeGen->Fill(ptJetGen - ptUeConeGen);
}
}
if(fListJets->GetEntries()>0 && fListJetsGen->GetEntries()>0){ //at least some reconstructed jets
+
Int_t ng = (Int_t) fListJetsGen->GetEntries();
Int_t nr = (Int_t) fListJets->GetEntries();
// Fill response matrix
for(Int_t ir = 0; ir < nr; ir++){
- AliAODJet *recJet = (AliAODJet*) fListJets->At(ir);
- Double_t etaJetRec = recJet->Eta();
- Double_t ptJetRec = recJet->Pt();
+ AliAODJet *recJet = (AliAODJet*) fListJets->At(ir);
+ Double_t etaJetRec = recJet->Eta();
+ Double_t ptJetRec = recJet->Pt();
+ Double_t areaJetRec = recJet->EffectiveAreaCharged();
//fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
if((fJetEtaMin <= etaJetRec) && (etaJetRec <= fJetEtaMax)){
Double_t ptJetGen = genJet->Pt();
Double_t etaJetGen = genJet->Eta();
- //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
+ //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
if((fJetEtaMin <= etaJetGen) && (etaJetGen <= fJetEtaMax)){
fhJetPtGenVsJetPtRec->Fill(ptJetRec, ptJetGen);
+
+ Double_t areaJetGen = genJet->EffectiveAreaCharged();
+ Double_t ptUeFromCellMedianGen = rhoFromCellMedianGen*areaJetGen;
+ Double_t ptUeConeGen = rhoConeGen*areaJetGen;
+ Double_t ptUeFromCellMedianRec = rhoFromCellMedian*areaJetRec;
+ Double_t ptUeConeRec = rhoCone*areaJetRec;
+ fhJetPtGenVsJetPtRecSubUeMedian->Fill(ptJetRec-ptUeFromCellMedianRec,
+ ptJetGen-ptUeFromCellMedianGen);
+ fhJetPtGenVsJetPtRecSubUeCone->Fill(ptJetRec-ptUeConeRec, ptJetGen-ptUeConeGen);
}
}//ig>=0
}//rec jet in eta acceptance
}//loop over reconstructed jets
}// # of rec jets >0
- }//pointer MC generator jets
-
- //=========================== Full jet charged jet matrix ==========================
- if(aodGenJetsFull){
- //Count full jets and charged-jet pairs at MC generator level
- for(Int_t ij=0; ij<fListJetsGenFull->GetEntries(); ij++){
- AliAODJet* jetFull = (AliAODJet*)(fListJetsGenFull->At(ij));
- if(!jetFull) continue;
- Double_t etaJetFull = jetFull->Eta();
- Double_t ptJetFull = jetFull->Pt();
-
- if((fJetEtaMin<=etaJetFull) && (etaJetFull<=fJetEtaMax)){
- fhJetPtGenFull->Fill(ptJetFull); // generator level pt spectum of full jets
- }
- }
- if(fListJetsGen->GetEntries()>0 && fListJetsGenFull->GetEntries()>0){ //at least some reconstructed jets
- Int_t nful = (Int_t) fListJetsGenFull->GetEntries();
- Int_t nchr = (Int_t) fListJetsGen->GetEntries();
-
- //Find closest MC generator full - charged jet
- if(faGenIndex.GetSize()<nchr) faGenIndex.Set(nchr); //idx of gen FULL jet assoc to gen CHARGED jet
- if(faRecIndex.GetSize()<nful) faRecIndex.Set(nful); //idx of gen CHARGED jet assoc to gen FULL jet
-
- if(fDebug){
- Printf("New Rec List %d gen index Array %d",nchr,faGenIndex.GetSize());
- Printf("New Gen List %d rec index Array %d",nful,faRecIndex.GetSize());
+
+ //=========================== Full jet vs charged jet matrix ==========================
+ if(fIsFullMC){
+ //Count full jets and charged-jet pairs at MC generator level
+ for(Int_t ij=0; ij<fListJetsGenFull->GetEntries(); ij++){
+ AliAODJet* jetFull = (AliAODJet*)(fListJetsGenFull->At(ij));
+ if(!jetFull) continue;
+ Double_t etaJetFull = jetFull->Eta();
+ Double_t ptJetFull = jetFull->Pt();
+
+ if((fJetEtaMin<=etaJetFull) && (etaJetFull<=fJetEtaMax)){
+ fhJetPtGenFull->Fill(ptJetFull); // generator level pt spectum of full jets
+ }
}
- //matching of MC genrator level and reconstructed jets
- AliAnalysisHelperJetTasks::GetClosestJets(fListJetsGenFull,nful,fListJetsGen,nchr,faGenIndex,faRecIndex,fDebug);
-
- // Fill response matrix
- for(Int_t ichr = 0; ichr < nchr; ichr++){ //charged jet loop
- AliAODJet *chJet = (AliAODJet*) fListJetsGen->At(ichr);
- Double_t etaJetCh = chJet->Eta();
- Double_t ptJetCh = chJet->Pt();
- //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
-
- if((fJetEtaMin <= etaJetCh) && (etaJetCh <= fJetEtaMax)){
- Int_t iful = faGenIndex[ichr]; //associated generator level jet
- if(iful >= 0 && iful < nful){
- if(fDebug > 10) Printf("%s:%d iful = %d ichr = %d",(char*)__FILE__,__LINE__,iful,ichr);
- AliAODJet *genJetFull = (AliAODJet*) fListJetsGenFull->At(iful);
- Double_t ptJetFull = genJetFull->Pt();
- Double_t etaJetFull = genJetFull->Eta();
-
+ if(fListJetsGen->GetEntries()>0 && fListJetsGenFull->GetEntries()>0){ //at least some reconstructed jets
+ Int_t nful = (Int_t) fListJetsGenFull->GetEntries();
+ Int_t nchr = (Int_t) fListJetsGen->GetEntries();
+
+ //Find closest MC generator full - charged jet
+ if(faGenIndex.GetSize()<nchr) faGenIndex.Set(nchr); //idx of gen FULL jet assoc to gen CHARGED jet
+ if(faRecIndex.GetSize()<nful) faRecIndex.Set(nful); //idx of gen CHARGED jet assoc to gen FULL jet
+
+ if(fDebug){
+ Printf("New Charg List %d Full index Array %d",nchr,faGenIndex.GetSize());
+ Printf("New Full List %d Charg index Array %d",nful,faRecIndex.GetSize());
+ }
+ //matching of MC genrator level and reconstructed jets
+ AliAnalysisHelperJetTasks::GetClosestJets(fListJetsGenFull,nful,fListJetsGen,nchr,faGenIndex,faRecIndex,fDebug);
+
+ // Fill response matrix
+ for(Int_t ichr = 0; ichr < nchr; ichr++){ //charged jet loop
+ AliAODJet *chJet = (AliAODJet*) fListJetsGen->At(ichr);
+ Double_t etaJetCh = chJet->Eta();
+ Double_t ptJetCh = chJet->Pt();
//fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
- if((fJetEtaMin <= etaJetFull) && (etaJetFull <= fJetEtaMax)){
- fhJetPtGenChargVsJetPtGenFull->Fill(ptJetFull,ptJetCh);
- }
- }//iful>=0
- }//rec jet in eta acceptance
- }//loop over reconstructed jets
- }// # of rec jets >0
- }//pointer MC generator jets
-
-
-
- } //analyze generator level MC
-
+
+ if((fJetEtaMin <= etaJetCh) && (etaJetCh <= fJetEtaMax)){
+ Int_t iful = faGenIndex[ichr]; //associated generator level jet
+ if(iful >= 0 && iful < nful){
+ if(fDebug > 10) Printf("%s:%d iful = %d ichr = %d",(char*)__FILE__,__LINE__,iful,ichr);
+ AliAODJet *genJetFull = (AliAODJet*) fListJetsGenFull->At(iful);
+ Double_t ptJetFull = genJetFull->Pt();
+ Double_t etaJetFull = genJetFull->Eta();
+
+ //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
+ if((fJetEtaMin <= etaJetFull) && (etaJetFull <= fJetEtaMax)){
+ fhJetPtGenChargVsJetPtGenFull->Fill(ptJetFull,ptJetCh);
+ }
+ }//iful>=0
+ }//rec jet in eta acceptance
+ }//loop over reconstructed jets
+ }// # of rec jets >0
+ }//pointer MC generator jets
+ } //fill resp mx only for bin
+ }//analyze generator level MC
+
+
+ if(fIsKine){ //skip reconstructed data analysis in case of kine
+ PostData(1, fOutputList);
+ return;
+ }
//============= RECONSTRUCTED INCLUSIVE JETS ===============
Double_t etaJet = 0.0;
Double_t pTJet = 0.0;
Double_t areaJet = 0.0;
Double_t phiJet = 0.0;
- Int_t indexLeadingJet = -1;
- Double_t pTLeadingJet = -10.0;
- Double_t areaLeadingJet = -10.0;
+ //Int_t indexLeadingJet = -1;
+ //Double_t pTLeadingJet = -10.0;
+ //Double_t areaLeadingJet = -10.0;
for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
if(pTJet==0) continue;
if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
- areaJet = jet->EffectiveAreaCharged();
+ /*areaJet = jet->EffectiveAreaCharged();*/
//Jet Diagnostics---------------------------------
fhJetPhi->Fill(pTJet, RelativePhi(phiJet,0.0)); //phi -pi,pi
fhJetEta->Fill(pTJet, etaJet);
//search for leading jet
- if(pTJet > pTLeadingJet){
+ /*if(pTJet > pTLeadingJet){
indexLeadingJet = ij;
pTLeadingJet = pTJet;
areaLeadingJet = areaJet;
pTJet,
areaJet
};
- fHJetPtRaw->Fill(fillraw);
+ fHJetPtRaw->Fill(fillraw);*/
}
-
+ /*
if(indexLeadingJet > -1){
// raw spectra of LEADING jets
//Centrality, pTjet, A
};
fHLeadingJetPtRaw->Fill(fillleading);
}
-
+ */
// =============== RECONSTRUCTED TRIGGER-JET PAIRS ================
- //Find Hadron trigger
- TList particleList; //list of tracks
- Int_t indexTrigg = GetListOfTracks(&particleList); //index of trigger hadron in Particle list
- if(fIsMC) FillEffHistos(&particleList, &particleListGen); //Fill efficiency histos
-
+ if(fIsChargedMC && fFillRespMx){
+ FillEffHistos(&particleList, &particleListGen); //Fill efficiency histos
+ }
+ Bool_t filledOnce = kTRUE; //fill rho histogram only once per event
//set ranges of the trigger loop
Int_t ilow = (fHardest==0 || fHardest==1) ? indexTrigg : 0;
Int_t ihigh = (fHardest==0 || fHardest==1) ? indexTrigg+1 : particleList.GetEntries();
//---------- make trigger-jet pairs ---------
- Int_t injet4 = 0;
- Int_t injet = 0;
-
+ //Int_t injet4 = 0;
+ //Int_t injet = 0;
+
for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
if(!jet) continue;
if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
areaJet = jet->EffectiveAreaCharged();
- if(areaJet >= 0.07) injet++;
- if(areaJet >= 0.4) injet4++;
+
+ //subtract bg using estinates base on median of kT jets
+ Double_t ptUeFromCellMedian = rhoFromCellMedian*areaJet;
+ Double_t ptUeCone = rhoCone*areaJet;
+
+ //if(areaJet >= 0.07) injet++;
+ //if(areaJet >= 0.4) injet4++;
Double_t dphi = RelativePhi(triggerHadron->Phi(), phiJet);
fhDphiTriggerJet->Fill(dphi); //Input
//Dphi versus jet pT
//Centrality, Dphi=phiTrig-phiJet, pTjet, pTtrigg
- Double_t filldp[] = { centValue,
+ /*Double_t filldp[] = { centValue,
dphi,
pTJet,
triggerHadron->Pt()
};
fHDphiVsJetPtAll->Fill(filldp);
-
+ */
+ //Rongrong's analysis
+
+ Double_t dPhi = phiJet - triggerHadron->Phi();
+ if(dPhi>2*TMath::Pi()) dPhi -= 2*TMath::Pi();
+ if(dPhi<-2*TMath::Pi()) dPhi += 2*TMath::Pi();
+ if(dPhi<-0.5*TMath::Pi()) dPhi += 2*TMath::Pi();
+ if(dPhi>1.5*TMath::Pi()) dPhi -= 2*TMath::Pi();
+
+ //[A,pTjet,pTtrig,dphi]
+ tmpArrayFour[0] = areaJet;
+ tmpArrayFour[1] = pTJet;
+ tmpArrayFour[2] = triggerHadron->Pt();
+ tmpArrayFour[3] = dPhi;
+
+ fHJetPhiCorr->Fill(tmpArrayFour); // Dphi distribution jet-triger
+
+ //[A,pTjet-pTUe,pTtrig,dphi]
+ tmpArrayFour[1] = pTJet - ptUeFromCellMedian;
+ fHJetPhiCorrSubUeMedian->Fill(tmpArrayFour);
+
+ //[A,pTjet-pTUe,pTtrig,dphi]
+ tmpArrayFour[1] = pTJet - ptUeCone;
+ fHJetPhiCorrSubUeCone->Fill(tmpArrayFour);
+
+ //Leticia's analysis
+
// Select back to back trigger - jet pairs
if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue;
fhDphiTriggerJetAccept->Fill(dphi); //Accepted
-
+
+ //NO bg subtraction
+ //Centrality, A, pTjet, pTtrigg, dphi
+ tmpArrayFive[0] = centValue;
+ tmpArrayFive[1] = areaJet;
+ tmpArrayFive[2] = pTJet;
+ tmpArrayFive[3] = triggerHadron->Pt();
+ tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
+ fHJetSpec->Fill(tmpArrayFive);
+
+ //Centrality, A, pTjet-pTbgCellMedian, pTtrigg, dphi
+ tmpArrayFive[0] = centValue;
+ tmpArrayFive[1] = areaJet;
+ tmpArrayFive[2] = pTJet-ptUeFromCellMedian;
+ tmpArrayFive[3] = triggerHadron->Pt();
+ tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
+ fHJetSpecSubUeMedian->Fill(tmpArrayFive);
+
+ //Centrality, A, pTjet-pTbgCone, pTtrigg, dphi
+ tmpArrayFive[0] = centValue;
+ tmpArrayFive[1] = areaJet;
+ tmpArrayFive[2] = pTJet - ptUeCone;
+ tmpArrayFive[3] = triggerHadron->Pt();
+ tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
+ fHJetSpecSubUeCone->Fill(tmpArrayFive);
+
+ //Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", kT median
+ tmpArrayFive[0] = areaJet;
+ tmpArrayFive[1] = pTJet;
+ tmpArrayFive[2] = pTJet - ptUeFromCellMedian;
+ tmpArrayFive[3] = ptUeFromCellMedian;
+ tmpArrayFive[4] = rhoFromCellMedian;
+ fHJetUeMedian->Fill(tmpArrayFive);
- //Centrality, A, pTjet, pTtrigg
- Double_t fillspec[] = { centValue,
- areaJet,
- pTJet,
- triggerHadron->Pt(),
- TMath::Abs((Double_t) dphi)
- };
- fHJetSpec->Fill(fillspec);
+ //Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", Cone median
+ tmpArrayFive[0] = areaJet;
+ tmpArrayFive[1] = pTJet;
+ tmpArrayFive[2] = pTJet - ptUeCone;
+ tmpArrayFive[3] = ptUeCone;
+ tmpArrayFive[4] = rhoCone;
+ fHJetUeCone->Fill(tmpArrayFive);
+
+ if(filledOnce){ //fill for each event only once
+ Double_t fillRho[] = { rhoCone,rhoFromCellMedian};
+ fHRhoUeMedianVsCone->Fill(fillRho);
+ filledOnce = kFALSE;
+ }
}//jet loop
//Fill Jet Density In the Event A>0.07
- if(injet>0){
+ /*if(injet>0){
Double_t filldens[]={ (Double_t) particleList.GetEntries(),
injet/fkAcceptance,
triggerHadron->Pt()
triggerHadron->Pt()
};
fHJetDensityA4->Fill(filldens4);
- }
+ }*/
}
for(Int_t it = 0; it < aodevt->GetNumberOfTracks(); it++){
AliAODTrack *tr = aodevt->GetTrack(it);
- //if((fFilterMask > 0) && !(tr->TestFilterBit(fFilterMask))) continue;
- if((fFilterMask > 0) && !(tr->IsHybridGlobalConstrainedGlobal())) continue;
+ if((fFilterMask > 0) && !(tr->TestFilterBit(fFilterMask))) continue;
+ //if((fFilterMask > 0) && !(tr->IsHybridGlobalConstrainedGlobal())) continue;
if(TMath::Abs((Float_t) tr->Eta()) > fTrackEtaCut) continue;
if(tr->Pt() < fTrackLowPtCut) continue;
list->Add(tr);
if(trk->Charge()==0) return kFALSE;
if(TMath::Abs((Float_t) trk->Eta()) > fTrackEtaCut) return kFALSE;
trkList->Add(trk);
- fhPtTrkTruePrimGen->Fill(trk->Pt(),trk->Eta()); //Efficiency denominator
+ if(fFillRespMx) fhPtTrkTruePrimGen->Fill(trk->Pt(),trk->Eta()); //Efficiency denominator
//Search trigger:
if(fHardest>0){ //leading particle
return;
}
+//________________________________________________________________
+void AliAnalysisTaskJetCorePP::EstimateBgRhoMedian(TList *listJet, TList* listPart, Double_t &rhoMedian, Int_t mode){
+ //Estimate background rho by means of integrating track pT outside identified jet cones
+
+ rhoMedian = 0;
+
+ //identify leading jet in the full track acceptance
+ Int_t idxLeadingJet = -1;
+ Double_t pTleading = 0.0;
+ AliAODJet* jet = NULL;
+
+ for(Int_t ij = 0; ij < listJet->GetEntries(); ij++){ //loop over bg kT jets
+ jet = (AliAODJet*)(listJet->At(ij));
+ if(!jet) continue;
+ if(TMath::Abs(jet->Eta()) > fTrackEtaCut) continue; //track acceptance cut
+ if(jet->Pt() > pTleading){
+ idxLeadingJet = ij;
+ pTleading = jet->Pt();
+ }
+ }
+
+ Int_t njets = 0;
+ static Double_t jphi[100];
+ static Double_t jeta[100];
+ static Double_t jRsquared[100];
+ for(Int_t ij = 0; ij< listJet->GetEntries(); ij++){ //loop over bg kT jets
+ jet = (AliAODJet*)(listJet->At(ij));
+ if(!jet) continue;
+ if(TMath::Abs(jet->Eta()) > fTrackEtaCut) continue;
+
+ if((jet->Pt() > fBgMaxJetPt || ij== idxLeadingJet) && njets<100){
+ jphi[njets] = RelativePhi(jet->Phi(),0.0); //-pi,pi
+ jeta[njets] = jet->Eta();
+ jRsquared[njets] = fSafetyMargin*jet->EffectiveAreaCharged()/TMath::Pi(); //scale jet radius
+ njets++;
+ }
+ }
+
+
+ static Double_t nOutCone[10][4];
+ static Double_t sumPtOutOfCone[10][4];
+
+ for(Int_t ie=0; ie<fnEta; ie++){
+ for(Int_t ip=0; ip<fnPhi; ip++){
+ nOutCone[ip][ie] = 0.0; //initialize counter
+ sumPtOutOfCone[ip][ie] = 0.0;
+ }
+ }
+
+ Double_t rndphi, rndeta;
+ Double_t rndphishift, rndetashift;
+ Double_t dphi, deta;
+ Bool_t bIsInCone;
+
+ for(Int_t it=0; it<fnTrials; it++){
+
+ rndphi = fRandom->Uniform(0, fPhiSize);
+ rndeta = fRandom->Uniform(0, fEtaSize);
+
+ for(Int_t ip=0; ip<fnPhi; ip++){ //move radom position to each cell
+ rndphishift = rndphi + ip*fPhiSize - TMath::Pi();
+ for(Int_t ie=0; ie<fnEta; ie++){
+ rndetashift = rndeta + ie*fEtaSize - fEtaSize;
+
+ bIsInCone = 0; //tag if trial is in the jet cone
+ for(Int_t ij=0; ij<njets; ij++){
+ deta = jeta[ij] - rndetashift;
+ dphi = RelativePhi(rndphishift,jphi[ij]);
+ if((dphi*dphi + deta*deta) < jRsquared[ij]){
+ bIsInCone = 1;
+ break;
+ }
+ }
+ if(!bIsInCone) nOutCone[ip][ie]++;
+ }
+ }
+ }
+
+
+ //Sum particle pT outside jets in cells
+ Int_t npart = listPart->GetEntries();
+ Int_t phicell,etacell;
+ AliVParticle *part = NULL;
+ for(Int_t ip=0; ip < npart; ip++){
+
+ part = (AliVParticle*) listPart->At(ip);
+ if(!part) continue;
+ if( TMath::Abs(part->Eta()) > fEtaSize) continue; //skip part outside +-0.7 in eta
+
+ bIsInCone = 0; //init
+ for(Int_t ij=0; ij<njets; ij++){
+ dphi = RelativePhi(jphi[ij], part->Phi());
+ deta = jeta[ij] - part->Eta();
+ if((dphi*dphi + deta*deta) < jRsquared[ij]){
+ bIsInCone = 1;
+ break;
+ }
+ }
+ if(!bIsInCone){
+ phicell = TMath::Nint(TMath::Floor((RelativePhi(part->Phi(),0.0) + TMath::Pi())/fPhiSize));
+ etacell = TMath::Nint(TMath::Floor((part->Eta()+fEtaSize)/fEtaSize));
+ sumPtOutOfCone[phicell][etacell]+= part->Pt();
+ }
+ }
+
+ static Double_t rhoInCells[20];
+ Double_t relativeArea;
+ Int_t nCells=0;
+ Double_t bufferArea=0.0, bufferPt=0.0; //sum cells where A< fJetFreeAreaFrac
+ for(Int_t ip=0; ip<fnPhi; ip++){
+ for(Int_t ie=0; ie<fnEta; ie++){
+ relativeArea = nOutCone[ip][ie]/fnTrials;
+ //cout<<"RA= "<< relativeArea<<" BA="<<bufferArea<<" BPT="<< bufferPt <<endl;
+
+ bufferArea += relativeArea;
+ bufferPt += sumPtOutOfCone[ip][ie];
+ if(bufferArea > fJetFreeAreaFrac){
+ rhoInCells[nCells] = bufferPt/(bufferArea*fCellArea);
+
+ if(mode==1) fhCellAreaToMedianGen->Fill(bufferArea*fCellArea);
+ else fhCellAreaToMedian->Fill(bufferArea*fCellArea);
+
+ bufferArea = 0.0;
+ bufferPt = 0.0;
+ nCells++;
+ }
+ }
+ }
+
+
+ if(nCells>0){
+ rhoMedian = TMath::Median(nCells, rhoInCells);
+ if(mode==1){ //mc data
+ fhEntriesToMedianGen->Fill(nCells);
+ }else{ //real data
+ fhEntriesToMedian->Fill(nCells);
+ }
+ }
+
+}
+
+//__________________________________________________________________
+void AliAnalysisTaskJetCorePP::EstimateBgCone(TList *listJet, TList* listPart, Double_t &rhoPerpCone){
+
+ rhoPerpCone = 0.0; //init
+
+ if(!listJet) return;
+ Int_t njet = listJet->GetEntries();
+ Int_t npart = listPart->GetEntries();
+ Double_t pTleading = 0.0;
+ Double_t phiLeading = 1000.;
+ Double_t etaLeading = 1000.;
+
+ for(Int_t jsig=0; jsig < njet; jsig++){
+ AliAODJet* signaljet = (AliAODJet*)(listJet->At(jsig));
+ if(!signaljet) continue;
+ if((signaljet->Eta()<fJetEtaMin) && (fJetEtaMax<signaljet->Eta())) continue; //acceptance cut
+ if(signaljet->Pt() >= pTleading){ //replace leading and subleading jet
+ pTleading = signaljet->Pt();
+ phiLeading = signaljet->Phi();
+ etaLeading = signaljet->Eta();
+ }
+ }
+
+ Double_t phileftcone = phiLeading + TMath::Pi()/2;
+ Double_t phirightcone = phiLeading - TMath::Pi()/2;
+ Double_t dp, de;
+
+ for(Int_t ip=0; ip < npart; ip++){
+
+ AliVParticle *part = (AliVParticle*) listPart->At(ip);
+ if(!part){
+ continue;
+ }
+
+ dp = RelativePhi(phileftcone, part->Phi());
+ de = etaLeading - part->Eta();
+ if( sqrt(dp*dp + de*de)< fBgConeR) rhoPerpCone += part->Pt();
+
+ dp = RelativePhi(phirightcone, part->Phi());
+ if( sqrt(dp*dp + de*de)< fBgConeR) rhoPerpCone += part->Pt();
+
+ }
+
+ //normalize total pT by two times cone are
+ rhoPerpCone = rhoPerpCone/(2*TMath::Pi()*fBgConeR*fBgConeR);
+
+}
+//__________________________________________________________________
+
+void AliAnalysisTaskJetCorePP::ReadTClonesArray(TString bname, TList *list){
+ //Convert TClones array of jets to TList
+
+ if(!strlen(bname.Data())){
+ AliError(Form("Jet branch %s not set.", bname.Data()));
+ return;
+ }
+
+ TClonesArray *array=0x0;
+ if(fUseExchContainer){ //take input from exchange containers
+ if(fJetBranchNameKine.Length()>0 && bname.CompareTo(fJetBranchNameKine.Data())==0){
+ array = dynamic_cast<TClonesArray*>(GetInputData(1)); //connect exchange container slot 1
+ }
+ if(fJetBranchNameBgKine.Length()>0 && bname.CompareTo(fJetBranchNameBgKine.Data())==0){
+ array = dynamic_cast<TClonesArray*>(GetInputData(2)); //connect exchange container slot 2
+ }
+ }else{ //take input from AOD
+ if(fAODOut&&!array){
+ array = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(bname.Data()));
+ }
+ if(fAODExtension&&!array){
+ array = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(bname.Data()));
+ }
+ if(fAODIn&&!array){
+ array = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(bname.Data()));
+ }
+ }
+
+ list->Clear(); //clear list beforehand
+
+ if(array){
+ if(fDebug){
+ Printf("########## %s: %d jets",bname.Data(), array->GetEntriesFast());
+ }
+
+ for(Int_t iJet = 0; iJet < array->GetEntriesFast(); iJet++) {
+ AliAODJet *jet = dynamic_cast<AliAODJet*>((*array)[iJet]);
+ if (jet) list->Add(jet);
+ }
+ }
+
+ return;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff