#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(""),
fJetBranchNameChargMC(""),
+fJetBranchNameKine(""),
fJetBranchNameFullMC(""),
fJetBranchNameBg(""),
fJetBranchNameBgChargMC(""),
+fJetBranchNameBgKine(""),
fListJets(0x0),
fListJetsGen(0x0),
fListJetsGenFull(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),
fHJetSpecSubUeMedian(0x0),
-fHJetSpecSubUeCMS(0x0),
+fHJetSpecSubUeCone(0x0),
+fHJetPhiCorr(0x0),
+fHJetPhiCorrSubUeMedian(0x0),
+fHJetPhiCorrSubUeCone(0x0),
fHJetUeMedian(0x0),
-fHJetUeCMS(0x0),
-fHRhoUeMedianVsCMS(0x0),
+fHJetUeCone(0x0),
+fHRhoUeMedianVsCone(0x0),
//fHJetDensity(0x0),
//fHJetDensityA4(0x0),
fhJetPhi(0x0),
//fHDphiVsJetPtAll(0x0),
fhJetPtGenVsJetPtRec(0x0),
fhJetPtGenVsJetPtRecSubUeMedian(0x0),
-fhJetPtGenVsJetPtRecSubUeCMS(0x0),
+fhJetPtGenVsJetPtRecSubUeCone(0x0),
fhJetPtGen(0x0),
fhJetPtSubUeMedianGen(0x0),
-fhJetPtSubUeCMSGen(0x0),
+fhJetPtSubUeConeGen(0x0),
fhJetPtGenChargVsJetPtGenFull(0x0),
fhJetPtGenFull(0x0),
fh2NtriggersGen(0x0),
fHJetSpecGen(0x0),
fHJetSpecSubUeMedianGen(0x0),
-fHJetSpecSubUeCMSGen(0x0),
+fHJetSpecSubUeConeGen(0x0),
+fHJetPhiCorrGen(0x0),
+fHJetPhiCorrSubUeMedianGen(0x0),
+fHJetPhiCorrSubUeConeGen(0x0),
fHJetUeMedianGen(0x0),
-fHJetUeCMSGen(0x0),
+fHJetUeConeGen(0x0),
fhPtTrkTruePrimRec(0x0),
fhPtTrkTruePrimGen(0x0),
fhPtTrkSecOrFakeRec(0x0),
-fHRhoUeMedianVsCMSGen(0x0),
+fHRhoUeMedianVsConeGen(0x0),
+fhEntriesToMedian(0x0),
+fhEntriesToMedianGen(0x0),
+fhCellAreaToMedian(0x0),
+fhCellAreaToMedianGen(0x0),
fIsChargedMC(0),
+fIsKine(0),
fIsFullMC(0),
faGenIndex(0),
faRecIndex(0),
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(""),
fJetBranchNameChargMC(""),
+fJetBranchNameKine(""),
fJetBranchNameFullMC(""),
fJetBranchNameBg(""),
fJetBranchNameBgChargMC(""),
+fJetBranchNameBgKine(""),
fListJets(0x0),
fListJetsGen(0x0),
fListJetsGenFull(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),
fHJetSpecSubUeMedian(0x0),
-fHJetSpecSubUeCMS(0x0),
+fHJetSpecSubUeCone(0x0),
+fHJetPhiCorr(0x0),
+fHJetPhiCorrSubUeMedian(0x0),
+fHJetPhiCorrSubUeCone(0x0),
fHJetUeMedian(0x0),
-fHJetUeCMS(0x0),
-fHRhoUeMedianVsCMS(0x0),
+fHJetUeCone(0x0),
+fHRhoUeMedianVsCone(0x0),
//fHJetDensity(0x0),
//fHJetDensityA4(0x0),
fhJetPhi(0x0),
//fHDphiVsJetPtAll(0x0),
fhJetPtGenVsJetPtRec(0x0),
fhJetPtGenVsJetPtRecSubUeMedian(0x0),
-fhJetPtGenVsJetPtRecSubUeCMS(0x0),
+fhJetPtGenVsJetPtRecSubUeCone(0x0),
fhJetPtGen(0x0),
fhJetPtSubUeMedianGen(0x0),
-fhJetPtSubUeCMSGen(0x0),
+fhJetPtSubUeConeGen(0x0),
fhJetPtGenChargVsJetPtGenFull(0x0),
fhJetPtGenFull(0x0),
fh2NtriggersGen(0x0),
fHJetSpecGen(0x0),
fHJetSpecSubUeMedianGen(0x0),
-fHJetSpecSubUeCMSGen(0x0),
+fHJetSpecSubUeConeGen(0x0),
+fHJetPhiCorrGen(0x0),
+fHJetPhiCorrSubUeMedianGen(0x0),
+fHJetPhiCorrSubUeConeGen(0x0),
fHJetUeMedianGen(0x0),
-fHJetUeCMSGen(0x0),
+fHJetUeConeGen(0x0),
fhPtTrkTruePrimRec(0x0),
fhPtTrkTruePrimGen(0x0),
fhPtTrkSecOrFakeRec(0x0),
-fHRhoUeMedianVsCMSGen(0x0),
+fHRhoUeMedianVsConeGen(0x0),
+fhEntriesToMedian(0x0),
+fhEntriesToMedianGen(0x0),
+fhCellAreaToMedian(0x0),
+fhCellAreaToMedianGen(0x0),
fIsChargedMC(0),
+fIsKine(0),
fIsFullMC(0),
faGenIndex(0),
faRecIndex(0),
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),
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),
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),
fHJetSpecSubUeMedian(a.fHJetSpecSubUeMedian),
-fHJetSpecSubUeCMS(a.fHJetSpecSubUeCMS),
+fHJetSpecSubUeCone(a.fHJetSpecSubUeCone),
+fHJetPhiCorr(a.fHJetPhiCorr),
+fHJetPhiCorrSubUeMedian(a.fHJetPhiCorrSubUeMedian),
+fHJetPhiCorrSubUeCone(a.fHJetPhiCorrSubUeCone),
fHJetUeMedian(a.fHJetUeMedian),
-fHJetUeCMS(a.fHJetUeCMS),
-fHRhoUeMedianVsCMS(a.fHRhoUeMedianVsCMS),
+fHJetUeCone(a.fHJetUeCone),
+fHRhoUeMedianVsCone(a.fHRhoUeMedianVsCone),
//fHJetDensity(a.fHJetDensity),
//fHJetDensityA4(a.fHJetDensityA4),
fhJetPhi(a.fhJetPhi),
//fHDphiVsJetPtAll(a.fHDphiVsJetPtAll),
fhJetPtGenVsJetPtRec(a.fhJetPtGenVsJetPtRec),
fhJetPtGenVsJetPtRecSubUeMedian(a.fhJetPtGenVsJetPtRecSubUeMedian),
-fhJetPtGenVsJetPtRecSubUeCMS(a.fhJetPtGenVsJetPtRecSubUeCMS),
+fhJetPtGenVsJetPtRecSubUeCone(a.fhJetPtGenVsJetPtRecSubUeCone),
fhJetPtGen(a.fhJetPtGen),
fhJetPtSubUeMedianGen(a.fhJetPtSubUeMedianGen),
-fhJetPtSubUeCMSGen(a.fhJetPtSubUeCMSGen),
+fhJetPtSubUeConeGen(a.fhJetPtSubUeConeGen),
fhJetPtGenChargVsJetPtGenFull(a.fhJetPtGenChargVsJetPtGenFull),
fhJetPtGenFull(a.fhJetPtGenFull),
fh2NtriggersGen(a.fh2NtriggersGen),
fHJetSpecGen(a.fHJetSpecGen),
fHJetSpecSubUeMedianGen(a.fHJetSpecSubUeMedianGen),
-fHJetSpecSubUeCMSGen(a.fHJetSpecSubUeCMSGen),
+fHJetSpecSubUeConeGen(a.fHJetSpecSubUeConeGen),
+fHJetPhiCorrGen(a.fHJetPhiCorrGen),
+fHJetPhiCorrSubUeMedianGen(a.fHJetPhiCorrSubUeMedianGen),
+fHJetPhiCorrSubUeConeGen(a.fHJetPhiCorrSubUeConeGen),
fHJetUeMedianGen(a.fHJetUeMedianGen),
-fHJetUeCMSGen(a.fHJetUeCMSGen),
+fHJetUeConeGen(a.fHJetUeConeGen),
fhPtTrkTruePrimRec(a.fhPtTrkTruePrimRec),
fhPtTrkTruePrimGen(a.fhPtTrkTruePrimGen),
fhPtTrkSecOrFakeRec(a.fhPtTrkSecOrFakeRec),
-fHRhoUeMedianVsCMSGen(a.fHRhoUeMedianVsCMSGen),
+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),
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);
//Implemented Notify() to read the cross sections
//and number of trials from pyxsec.root
//inspired by AliAnalysisTaskJetSpectrum2::Notify()
- if(!fIsChargedMC) 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 antikt jets
- fListJetsBg = new TList(); //reconstructed level bg kT jets
+ fListJetsBg = new TList(); //reconstructed jets to be removed from bg
fIsChargedMC = (fJetBranchNameChargMC.Length()>0) ? kTRUE : kFALSE;
+ fIsKine = (fJetBranchNameKine.Length()>0) ? kTRUE : kFALSE;
fIsFullMC = (fJetBranchNameFullMC.Length()>0) ? kTRUE : kFALSE;
fRandom = new TRandom3(0);
- if(fIsChargedMC){
+ if(fIsChargedMC || fIsKine){ //full MC or pure kine
fListJetsGen = new TList(); //generator level charged antikt jets
- fListJetsBgGen = new TList(); //generator level charged bg kT jets
+ fListJetsBgGen = new TList(); //generator level jets to be removed from bg
if(fIsFullMC)
fListJetsGenFull = new TList(); //generator level full jets
//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, 220, 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()};
+ 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]");
- fOutputList->Add(fHJetSpecSubUeMedian);
- //background estimated as weighted median of kT jets ala CMS
- fHJetSpecSubUeCMS = (THnSparseF*) fHJetSpec->Clone("fHJetSpecSubUeCMS");
- fHJetSpecSubUeCMS->SetTitle("Recoil jet spectrum [cent,A,pTjet-pTUe,pTtrig,dphi]");
- fOutputList->Add(fHJetSpecSubUeCMS);
+ 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] = {50, 50, 50, 50, 50};
- const Double_t lowBinSpecMed[dimSpecMed] = {0.0, 0.0, -50.0, 0.0, 0.0};
- const Double_t hiBinSpecMed[dimSpecMed] = {1.0, 100.0, 100.0, 20.0, 20.0};
+ 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);
- fOutputList->Add(fHJetUeMedian);
+ if(!fIsKine) fOutputList->Add(fHJetUeMedian);
- //A, pTjet, pTjet-pTUe, pTUe, rhoUe bg estimate from kT median CMS
- fHJetUeCMS = (THnSparseF*) fHJetUeMedian->Clone("fHJetUeCMS");
- fHJetUeCMS->SetTitle("Recoil jet spectrum [A,pTjet,pTjet-pTUe, pTUe, rhoUe]");
- fOutputList->Add(fHJetUeCMS);
+ //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] = {200 , 200};
+ 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};
- fHRhoUeMedianVsCMS = new THnSparseF("hRhoUeMedianVsCMS","[Rho CMS, Rho Median]",
+ fHRhoUeMedianVsCone = new THnSparseF("hRhoUeMedianVsCone","[Rho Cone, Rho Median]",
dimRho, nBinsRho, lowBinRho, hiBinRho);
- fOutputList->Add(fHRhoUeMedianVsCMS);
+ if(!fIsKine) fOutputList->Add(fHRhoUeMedianVsCone);
//Jet number density histos [Trk Mult, jet density, pT trigger]
/*const Int_t dimJetDens = 3;
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;
*/
//analyze MC generator level
- if(fIsChargedMC){
- fhJetPtGenVsJetPtRec = new TH2D("fhJetPtGenVsJetPtRec","JetPtGenVsJetPtRec", 200,0,200, 200,0,200);
- fOutputList->Add(fhJetPtGenVsJetPtRec); //gen MC charg jet pt spectrum versus rec charged jet pt spectrum
- //....
- fhJetPtGenVsJetPtRecSubUeMedian = new TH2D("fhJetPtGenVsJetPtRecSubUeMedian","fhJetPtGenVsJetPtRecSubUeMedian", 220,-20,200, 220,-20,200);
- fOutputList->Add(fhJetPtGenVsJetPtRecSubUeMedian); // with kT median bg subtr
- //....
- fhJetPtGenVsJetPtRecSubUeCMS=(TH2D*)fhJetPtGenVsJetPtRecSubUeMedian->Clone("fhJetPtGenVsJetPtRecSubUeCMS");
- fhJetPtGenVsJetPtRecSubUeCMS->SetTitle("fhJetPtGenVsJetPtRecSubUeCMS");
- fOutputList->Add(fhJetPtGenVsJetPtRecSubUeCMS); // with weighted kT median bg subtr
- //....
- fhJetPtGen = new TH1D("fhJetPtGen","Jet Pt (MC Gen)",200,0,200); //MC generator charged jet pt spectrum
- fOutputList->Add(fhJetPtGen);
- //....
- fhJetPtSubUeMedianGen = new TH1D("fhJetPtSubUeMedianGen","Jet Pt - UE pT (MC Gen)",220,-20,200);
- fOutputList->Add(fhJetPtSubUeMedianGen); // with kT median bg subtr
- //....
- fhJetPtSubUeCMSGen = (TH1D*) fhJetPtSubUeMedianGen->Clone("fhJetPtSubUeCMSGen");
- fOutputList->Add(fhJetPtSubUeCMSGen); // with weighted kT median bg subtr
- //....
- if(fIsFullMC){
- fhJetPtGenChargVsJetPtGenFull = new TH2D("fhJetPtGenChargVsJetPtGenFull","fhJetPtGenChargVsJetPtGenFull", 200,0,200, 200,0,200);
- fOutputList->Add(fhJetPtGenChargVsJetPtGenFull); //gen full MC jet pt versus gen charged jet MC pt
+ 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
//....
- fhJetPtGenFull = new TH1D("fhJetPtGenFull","Jet Pt (MC Full jets Gen)",200,0,200); //MC generator full jet pt spectrum
- fOutputList->Add(fhJetPtGenFull);
+ 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");
fHJetSpecSubUeMedianGen->SetTitle(Form("%s Gen MC",fHJetSpecSubUeMedian->GetTitle()));
fOutputList->Add(fHJetSpecSubUeMedianGen);
- fHJetSpecSubUeCMSGen = (THnSparseF*) fHJetSpecSubUeCMS->Clone("fHJetSpecSubUeCMSGen");
- fHJetSpecSubUeCMSGen->SetTitle(Form("%s Gen MC",fHJetSpecSubUeCMS->GetTitle()));
- fOutputList->Add(fHJetSpecSubUeCMSGen);
+ 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);
- fHJetUeCMSGen = (THnSparseF*) fHJetUeCMS->Clone("fHJetUeCMSGen");
- fHJetUeCMSGen->SetTitle(Form("%s Gen MC", fHJetUeCMS->GetTitle()));
- fOutputList->Add(fHJetUeCMSGen);
+ fHJetUeConeGen = (THnSparseF*) fHJetUeCone->Clone("fHJetUeConeGen");
+ fHJetUeConeGen->SetTitle(Form("%s Gen MC", fHJetUeCone->GetTitle()));
+ fOutputList->Add(fHJetUeConeGen);
- //track efficiency/contamination histograms eta versus pT
- fhPtTrkTruePrimRec = new TH2D("fhPtTrkTruePrimRec","PtTrkTruePrimRec",100,0,20,18,-0.9,0.9);
- fOutputList->Add(fhPtTrkTruePrimRec);
+ 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);
- fHRhoUeMedianVsCMSGen = (THnSparseF*) fHRhoUeMedianVsCMS->Clone("hRhoUeMedianVsCMSGen");
- fHRhoUeMedianVsCMSGen->SetTitle(Form("%s Gen MC", fHRhoUeMedianVsCMS->GetTitle()));
- fOutputList->Add(fHRhoUeMedianVsCMSGen);
+ 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(fIsChargedMC) 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(!strlen(fJetBranchNameBg.Data())){
- AliError("Jet bg branch name not set.");
- 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);
}
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 ============================
Double_t tmpArrayFive[5];
-
- //=============== Recnstricted antikt and kt jets ===============
- ReadTClonesArray(fJetBranchName.Data() , fListJets);
- ReadTClonesArray(fJetBranchNameBg.Data(), fListJetsBg);
+ Double_t tmpArrayFour[4];
+
- //============ Estimate background in reconstructed events ===========
- Double_t rhoFromKtMedian=0.0, rhoAlaCMS=0.0;
- EstimateBgRhoAlaCMS(fListJetsBg, fListJets, rhoFromKtMedian, rhoAlaCMS);
- fListJetsBg->Clear(); //this list is further not needed
+ 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);
+ //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(fIsChargedMC){
- //================= generated charged antikt and kt jets ================
- ReadTClonesArray(fJetBranchNameChargMC.Data(), fListJetsGen);
- ReadTClonesArray(fJetBranchNameBgChargMC.Data(), fListJetsBgGen);
- if(fIsFullMC){ //generated full jets
- ReadTClonesArray(fJetBranchNameFullMC.Data(), fListJetsGenFull);
- }
+ if(fIsChargedMC || fIsKine){
- //============== Estimate bg in generated events ==============
- Double_t rhoFromKtMedianGen=0.0, rhoAlaCMSGen=0.0;
- EstimateBgRhoAlaCMS(fListJetsBgGen, fListJetsGen, rhoFromKtMedianGen, rhoAlaCMSGen);
- fListJetsBgGen->Clear();
+ 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(fIsFullMC){ //generated full jets
+ ReadTClonesArray(fJetBranchNameFullMC.Data(), fListJetsGenFull);
+ }
+ }
//========================================================
//serarch for charged trigger at the MC generator level
Int_t triggersMC[200];//list of trigger candidates
Int_t ntriggersMC = 0; //index in triggers array
- if(fESD){//ESD input
+ if(!fIsKine){
+ if(fESD){//ESD input
- AliMCEvent* mcEvent = MCEvent();
- if(!mcEvent){
- PostData(1, fOutputList);
- return;
- }
+ 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);
+ }
- 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
+ }
+ }
}
+ }else{ //analyze kine
- for(Int_t it = 0; it < mcEvent->GetNumberOfTracks(); it++){
- if(!mcEvent->IsPhysicalPrimary(it)) continue;
- AliMCParticle* part = (AliMCParticle*) mcEvent->GetTrack(it);
+ 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 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)){
+ }
+
+ //============== Estimate bg in generated events ==============
+ Double_t rhoFromCellMedianGen=0.0, rhoConeGen=0.0;
- if(fHardest==0 && ntriggersMC<200){//single inclusive trigger
- if(indexTriggGen > -1){ //trigger candidater was found
- triggersMC[ntriggersMC] = indexTriggGen;
- ntriggersMC++;
- }
- }
+ //Estimate rho from cone
+ EstimateBgCone(fListJetsGen, &particleListGen, rhoConeGen);
- iCounterGen++;//number of entries in particleListGen array
- }
- }
- }
+ //Estimate rho from cell median minus jets
+ EstimateBgRhoMedian(fListJetsBgGen, &particleListGen, rhoFromCellMedianGen,1);//mc data
//============ Generator trigger+jet ==================
if(fHardest==0){ //single inclusive trigger
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;
tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
fHJetSpecGen->Fill(tmpArrayFive);
-
- Double_t ptUeFromKtMedianGen = rhoFromKtMedianGen*areaJetGen;
- Double_t ptUeAlaCMSGen = rhoAlaCMSGen*areaJetGen;
-
- //Centrality, A, pTjet-pTbgKtMedian, pTtrigg, dphi
+ //Centrality, A, pTjet-pTbgCellMedian, pTtrigg, dphi
tmpArrayFive[0] = centValue;
tmpArrayFive[1] = areaJetGen;
- tmpArrayFive[2] = jet->Pt() - ptUeFromKtMedianGen;
+ tmpArrayFive[2] = jet->Pt() - ptUeFromCellMedianGen;
tmpArrayFive[3] = ptTriggGen;
tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
fHJetSpecSubUeMedianGen->Fill(tmpArrayFive);
- //Centrality, A, pTjet-pTbgCMS, pTtrigg, dphi
+ //Centrality, A, pTjet-pTbgCone, pTtrigg, dphi
tmpArrayFive[0] = centValue;
tmpArrayFive[1] = areaJetGen;
- tmpArrayFive[2] = jet->Pt() - ptUeAlaCMSGen;
+ tmpArrayFive[2] = jet->Pt() - ptUeConeGen;
tmpArrayFive[3] = ptTriggGen;
tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
- fHJetSpecSubUeCMSGen->Fill(tmpArrayFive);
+ fHJetSpecSubUeConeGen->Fill(tmpArrayFive);
//Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", kT median
tmpArrayFive[0] = areaJetGen;
tmpArrayFive[1] = jet->Pt();
- tmpArrayFive[2] = jet->Pt() - ptUeFromKtMedianGen;
- tmpArrayFive[3] = ptUeFromKtMedianGen;
- tmpArrayFive[4] = rhoFromKtMedianGen;
+ tmpArrayFive[2] = jet->Pt() - ptUeFromCellMedianGen;
+ tmpArrayFive[3] = ptUeFromCellMedianGen;
+ tmpArrayFive[4] = rhoFromCellMedianGen;
fHJetUeMedianGen->Fill(tmpArrayFive);
- //Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", CMS median
+ //Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", perpendicular Cone
tmpArrayFive[0] = areaJetGen;
tmpArrayFive[1] = jet->Pt();
- tmpArrayFive[2] = jet->Pt() - ptUeAlaCMSGen;
- tmpArrayFive[3] = ptUeAlaCMSGen;
- tmpArrayFive[4] = rhoAlaCMSGen;
- fHJetUeCMSGen->Fill(tmpArrayFive);
+ tmpArrayFive[2] = jet->Pt() - ptUeConeGen;
+ tmpArrayFive[3] = ptUeConeGen;
+ tmpArrayFive[4] = rhoConeGen;
+ fHJetUeConeGen->Fill(tmpArrayFive);
if(fillOnceGen){
- Double_t fillRhoGen[] = { rhoAlaCMSGen,rhoFromKtMedianGen};
- fHRhoUeMedianVsCMSGen->Fill(fillRhoGen);
+ Double_t fillRhoGen[] = { rhoConeGen, rhoFromCellMedianGen};
+ fHRhoUeMedianVsConeGen->Fill(fillRhoGen);
fillOnceGen = kFALSE;
}
}//back to back jet-trigger pair
}//jet loop
}//trigger loop
- //================ 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(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); // gen level pt spectum of jets response mx normalization
-
- Double_t areaJetGen = jet->EffectiveAreaCharged();
- Double_t ptUeFromKtMedianGen = rhoFromKtMedianGen*areaJetGen;
- Double_t ptUeAlaCMSGen = rhoAlaCMSGen*areaJetGen;
-
- fhJetPtSubUeMedianGen->Fill(ptJetGen - ptUeFromKtMedianGen);
- fhJetPtSubUeCMSGen->Fill(ptJetGen - ptUeAlaCMSGen);
- }
- }
- 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();
-
- //Find closest MC generator - reconstructed jets
- if(faGenIndex.GetSize()<nr) faGenIndex.Set(nr); //idx of gen jet assoc to rec jet
- if(faRecIndex.GetSize()<ng) faRecIndex.Set(ng); //idx of rec jet assoc to gen jet
+ if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){
+ fhJetPtGen->Fill(ptJetGen); // gen level pt spectum of jets response mx normalization
- if(fDebug){
- Printf("New Rec List %d gen index Array %d",nr,faGenIndex.GetSize());
- Printf("New Gen List %d rec index Array %d",ng,faRecIndex.GetSize());
- }
- //matching of MC genrator level and reconstructed jets
- AliAnalysisHelperJetTasks::GetClosestJets(fListJetsGen,ng,fListJets,nr,faGenIndex,faRecIndex,fDebug);
-
- // 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();
- Double_t areaJetRec = recJet->EffectiveAreaCharged();
- //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
-
- if((fJetEtaMin <= etaJetRec) && (etaJetRec <= fJetEtaMax)){
- Int_t ig = faGenIndex[ir]; //associated generator level jet
- if(ig >= 0 && ig < ng){
- if(fDebug > 10) Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
- AliAODJet *genJet = (AliAODJet*) fListJetsGen->At(ig);
- Double_t ptJetGen = genJet->Pt();
- Double_t etaJetGen = genJet->Eta();
+ Double_t areaJetGen = jet->EffectiveAreaCharged();
+ Double_t ptUeFromCellMedianGen = rhoFromCellMedianGen*areaJetGen;
+ Double_t ptUeConeGen = rhoConeGen*areaJetGen;
- //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 ptUeFromKtMedianGen = rhoFromKtMedianGen*areaJetGen;
- Double_t ptUeAlaCMSGen = rhoAlaCMSGen*areaJetGen;
- Double_t ptUeFromKtMedianRec = rhoFromKtMedian*areaJetRec;
- Double_t ptUeAlaCMSRec = rhoAlaCMS*areaJetRec;
- fhJetPtGenVsJetPtRecSubUeMedian->Fill(ptJetRec-ptUeFromKtMedianRec,
- ptJetGen-ptUeFromKtMedianGen);
- fhJetPtGenVsJetPtRecSubUeCMS->Fill(ptJetRec-ptUeAlaCMSRec, ptJetGen-ptUeAlaCMSGen);
- }
- }//ig>=0
- }//rec jet in eta acceptance
- }//loop over reconstructed jets
- }// # of rec jets >0
-
- //=========================== 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
+ fhJetPtSubUeMedianGen->Fill(ptJetGen - ptUeFromCellMedianGen);
+ fhJetPtSubUeConeGen->Fill(ptJetGen - ptUeConeGen);
}
}
- 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();
+ 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();
- //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
+ //Find closest MC generator - reconstructed jets
+ if(faGenIndex.GetSize()<nr) faGenIndex.Set(nr); //idx of gen jet assoc to rec jet
+ if(faRecIndex.GetSize()<ng) faRecIndex.Set(ng); //idx of rec jet assoc to gen 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());
+ Printf("New Rec List %d gen index Array %d",nr,faGenIndex.GetSize());
+ Printf("New Gen List %d rec index Array %d",ng,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();
+ AliAnalysisHelperJetTasks::GetClosestJets(fListJetsGen,ng,fListJets,nr,faGenIndex,faRecIndex,fDebug);
+
+ // 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();
+ Double_t areaJetRec = recJet->EffectiveAreaCharged();
//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((fJetEtaMin <= etaJetRec) && (etaJetRec <= fJetEtaMax)){
+ Int_t ig = faGenIndex[ir]; //associated generator level jet
+ if(ig >= 0 && ig < ng){
+ if(fDebug > 10) Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
+ AliAODJet *genJet = (AliAODJet*) fListJetsGen->At(ig);
+ Double_t ptJetGen = genJet->Pt();
+ Double_t etaJetGen = genJet->Eta();
//fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
- if((fJetEtaMin <= etaJetFull) && (etaJetFull <= fJetEtaMax)){
- fhJetPtGenChargVsJetPtGenFull->Fill(ptJetFull,ptJetCh);
+ 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);
}
- }//iful>=0
+ }//ig>=0
}//rec jet in eta acceptance
}//loop over reconstructed jets
}// # of rec jets >0
- }//pointer MC generator jets
-
- } //analyze generator level MC
-
-
-
-
+
+ //=========================== 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
+ }
+ }
+ 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 <= 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));
*/
// =============== 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(fIsChargedMC) 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;
//---------- make trigger-jet pairs ---------
//Int_t injet4 = 0;
//Int_t injet = 0;
- //Double_t rhoFromKtMedian=0.0, rhoAlaCMS=0.0;
- //EstimateBgRhoAlaCMS(fListJetsBg, fListJets, rhoFromKtMedian, rhoAlaCMS);
for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
areaJet = jet->EffectiveAreaCharged();
+
+ //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++;
};
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
tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
fHJetSpec->Fill(tmpArrayFive);
- //subtract bg using estinates base on median of kT jets
- Double_t ptUeFromKtMedian = rhoFromKtMedian*areaJet;
- Double_t ptUeAlaCMS = rhoAlaCMS*areaJet;
-
- //Centrality, A, pTjet-pTbgKtMedian, pTtrigg, dphi
+ //Centrality, A, pTjet-pTbgCellMedian, pTtrigg, dphi
tmpArrayFive[0] = centValue;
tmpArrayFive[1] = areaJet;
- tmpArrayFive[2] = pTJet-ptUeFromKtMedian;
+ tmpArrayFive[2] = pTJet-ptUeFromCellMedian;
tmpArrayFive[3] = triggerHadron->Pt();
tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
fHJetSpecSubUeMedian->Fill(tmpArrayFive);
- //Centrality, A, pTjet-pTbgCMS, pTtrigg, dphi
+ //Centrality, A, pTjet-pTbgCone, pTtrigg, dphi
tmpArrayFive[0] = centValue;
tmpArrayFive[1] = areaJet;
- tmpArrayFive[2] = pTJet - ptUeAlaCMS;
+ tmpArrayFive[2] = pTJet - ptUeCone;
tmpArrayFive[3] = triggerHadron->Pt();
tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
- fHJetSpecSubUeCMS->Fill(tmpArrayFive);
+ fHJetSpecSubUeCone->Fill(tmpArrayFive);
//Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", kT median
tmpArrayFive[0] = areaJet;
tmpArrayFive[1] = pTJet;
- tmpArrayFive[2] = pTJet - ptUeFromKtMedian;
- tmpArrayFive[3] = ptUeFromKtMedian;
- tmpArrayFive[4] = rhoFromKtMedian;
+ tmpArrayFive[2] = pTJet - ptUeFromCellMedian;
+ tmpArrayFive[3] = ptUeFromCellMedian;
+ tmpArrayFive[4] = rhoFromCellMedian;
fHJetUeMedian->Fill(tmpArrayFive);
- //Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", CMS median
+ //Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", Cone median
tmpArrayFive[0] = areaJet;
tmpArrayFive[1] = pTJet;
- tmpArrayFive[2] = pTJet - ptUeAlaCMS;
- tmpArrayFive[3] = ptUeAlaCMS;
- tmpArrayFive[4] = rhoAlaCMS;
- fHJetUeCMS->Fill(tmpArrayFive);
+ tmpArrayFive[2] = pTJet - ptUeCone;
+ tmpArrayFive[3] = ptUeCone;
+ tmpArrayFive[4] = rhoCone;
+ fHJetUeCone->Fill(tmpArrayFive);
if(filledOnce){ //fill for each event only once
- Double_t fillRho[] = { rhoAlaCMS,rhoFromKtMedian};
- fHRhoUeMedianVsCMS->Fill(fillRho);
+ Double_t fillRho[] = { rhoCone,rhoFromCellMedian};
+ fHRhoUeMedianVsCone->Fill(fillRho);
filledOnce = kFALSE;
}
}//jet loop
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::EstimateBgRhoAlaCMS(TList *listJetBg, TList *listJet, Double_t &rhoMedian, Double_t& rhoImprovedCMS){
- //CMS method to estimate bg
- //adopted from Ruedinger
- // http://portal.nersc.gov/project/alice/htmldoc/src/AliAnalysisTaskChargedJetsPA.cxx.html#nP_gSD
- //AliAnalysisTaskChargedJetsPA::GetKTBackgroundDensity
- if(!listJetBg) return;
-
- static Double_t tmpRhoImprovedCMS[1024];
- Double_t tmpCoveredArea = 0.0;
-
- // Initialize
- rhoMedian = 0.0;
- rhoImprovedCMS = 0.0;
- Int_t rhoImprovedCMSJetCount = 0;
-
- //Find two leading signal jets in the acceptance
- const Int_t kSigJets = 2;
- Int_t idxLeadingJets[kSigJets]={-1,-1};
- Double_t pTleading=-1., pTsubleading=-1.;
- if(listJet){
- for(Int_t jsig=0; jsig < listJet->GetEntries(); 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
- pTsubleading = pTleading;
- idxLeadingJets[1]=idxLeadingJets[0];
- pTleading = signaljet->Pt();
- idxLeadingJets[0]=jsig;
- }else if( signaljet->Pt() > pTsubleading){ //replace subleading jet
- pTsubleading = signaljet->Pt();
- idxLeadingJets[1]=jsig;
- }
+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++;
}
}
- //____
- for(Int_t ibg = 0; ibg< listJetBg->GetEntries(); ibg++){ //loop over bg kT jets
-
- AliAODJet* bgjet = (AliAODJet*)(listJetBg->At(ibg));
- if(!bgjet) continue;
- if((bgjet->Eta()<fJetEtaMin) && (fJetEtaMax<bgjet->Eta())) continue; //acceptance cut
-
- Int_t nbgtracks = (bgjet->GetRefTracks())->GetLast()+1;
-
- // Search for overlap with signal jets
- Bool_t isOverlapping = kFALSE;
- if(listJet){
- for(Int_t j=0;j<kSigJets;j++){ //skip jets wich share particles with the two leading jets
- if(idxLeadingJets[j]<0) continue;
- AliAODJet* signaljet = (AliAODJet*)(listJet->At(idxLeadingJets[j]));
- if(!signaljet) continue;
- Int_t nsignaltracks = (signaljet->GetRefTracks())->GetLast()+1;
+
- for(Int_t itrkbg=0; itrkbg < nbgtracks; itrkbg++ ){
- AliVParticle *trkbg = dynamic_cast<AliVParticle*> ( bgjet->GetTrack(itrkbg));
- if(!trkbg) continue;
- for(Int_t itrksig=0; itrksig < nsignaltracks; itrksig++ ){
- AliVParticle *trksig = dynamic_cast<AliVParticle*> ( signaljet->GetTrack(itrksig));
- if(!trksig) continue;
- if(trkbg->GetLabel() == trksig->GetLabel()){
- isOverlapping = kTRUE; // check phi and eta of the tracks id it is the same
- break;
- }
- }//loop over antikt jet tracks
- if(isOverlapping) break;
- }//loop over jet bg tracks
- if(isOverlapping) break;
- }//loop over 2 leading antikt jets
+ 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;
}
-
- if(!bgjet->EffectiveAreaCharged()) continue;
+ }
- if(bgjet->Pt() > 0.150) tmpCoveredArea += bgjet->EffectiveAreaCharged();
+ 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]++;
+ }
+ }
+ }
+
- if(bgjet->Pt() > 0.150 && !isOverlapping){
- Double_t tmpRho = bgjet->Pt() / bgjet->EffectiveAreaCharged();
- tmpRhoImprovedCMS[rhoImprovedCMSJetCount] = tmpRho;
- rhoImprovedCMSJetCount++;
+ //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;
+ }
}
- }//loop over bg jets
+ 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();
+ }
+ }
- if(rhoImprovedCMSJetCount > 0){
- rhoMedian = TMath::Median(rhoImprovedCMSJetCount, tmpRhoImprovedCMS);
- rhoImprovedCMS = rhoMedian * tmpCoveredArea/fkAcceptance; //rescale to full akceptance
+ 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){
}
TClonesArray *array=0x0;
-
- 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()));
+ 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(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);
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff