Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add PF HF sum energy to reco::Centrality #44518

Merged
merged 3 commits into from
Apr 3, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 7 additions & 0 deletions DQM/Physics/src/CentralityDQM.cc
Original file line number Diff line number Diff line change
Expand Up @@ -84,6 +84,9 @@ void CentralityDQM::bookHistograms(DQMStore::IBooker& bei, edm::Run const&, edm:
h_hiZDC = bei.book1D("h_hiZDC", "h_hiZDC", 600, 0, 6000);
h_hiZDCplus = bei.book1D("h_hiZDCplus", "h_hiZDCplus", 600, 0, 6000);
h_hiZDCminus = bei.book1D("h_hiZDCminus", "h_hiZDCminus", 600, 0, 6000);
h_hiPF = bei.book1D("h_hiPF", "h_hiPF", 900, 0, 9000);
h_hiPFplus = bei.book1D("h_hiPFplus", "h_hiPFplus", 900, 0, 9000);
h_hiPFminus = bei.book1D("h_hiPFminus", "h_hiPFminus", 900, 0, 9000);

h_vertex_x = bei.book1D("h_vertex_x", "h_vertex_x", 400, -4, 4);
h_vertex_y = bei.book1D("h_vertex_y", "h_vertex_y", 400, -4, 4);
Expand Down Expand Up @@ -155,6 +158,10 @@ void CentralityDQM::analyze(const edm::Event& iEvent, const edm::EventSetup& iSe
h_hiEB->Fill(cent->EtEBSum());
h_hiET->Fill(cent->EtMidRapiditySum());

h_hiPF->Fill(cent->EtPFhfSum());
h_hiPFplus->Fill(cent->EtPFhfSumPlus());
h_hiPFminus->Fill(cent->EtPFhfSumMinus());

edm::Handle<std::vector<reco::Vertex> > vertex;
iEvent.getByToken(vertexToken, vertex);
h_vertex_x->Fill(vertex->begin()->x());
Expand Down
3 changes: 3 additions & 0 deletions DQM/Physics/src/CentralityDQM.h
Original file line number Diff line number Diff line change
Expand Up @@ -76,6 +76,9 @@ class CentralityDQM : public DQMEDAnalyzer {
MonitorElement* h_hiZDC;
MonitorElement* h_hiZDCplus;
MonitorElement* h_hiZDCminus;
MonitorElement* h_hiPF;
MonitorElement* h_hiPFplus;
MonitorElement* h_hiPFminus;

MonitorElement* h_vertex_x;
MonitorElement* h_vertex_y;
Expand Down
7 changes: 7 additions & 0 deletions DQM/Physics/src/CentralitypADQM.cc
Original file line number Diff line number Diff line change
Expand Up @@ -77,6 +77,9 @@ void CentralitypADQM::bookHistograms(DQMStore::IBooker& bei, edm::Run const&, ed
h_hiZDC = bei.book1D("h_hiZDC", "h_hiZDC", 600, 0, 6000);
h_hiZDCplus = bei.book1D("h_hiZDCplus", "h_hiZDCplus", 600, 0, 6000);
h_hiZDCminus = bei.book1D("h_hiZDCminus", "h_hiZDCminus", 600, 0, 6000);
h_hiPF = bei.book1D("h_hiPF", "h_hiPF", 900, 0, 9000);
h_hiPFplus = bei.book1D("h_hiPFplus", "h_hiPFplus", 900, 0, 9000);
h_hiPFminus = bei.book1D("h_hiPFminus", "h_hiPFminus", 900, 0, 9000);

h_vertex_x = bei.book1D("h_vertex_x", "h_vertex_x", 400, -4, 4);
h_vertex_y = bei.book1D("h_vertex_y", "h_vertex_y", 400, -4, 4);
Expand Down Expand Up @@ -118,6 +121,10 @@ void CentralitypADQM::analyze(const edm::Event& iEvent, const edm::EventSetup& i
h_hiEB->Fill(cent->EtEBSum());
h_hiET->Fill(cent->EtMidRapiditySum());

h_hiPF->Fill(cent->EtPFhfSum());
h_hiPFplus->Fill(cent->EtPFhfSumPlus());
h_hiPFminus->Fill(cent->EtPFhfSumMinus());

edm::Handle<std::vector<reco::Vertex> > vertex;
iEvent.getByToken(vertexToken, vertex);
h_vertex_x->Fill(vertex->begin()->x());
Expand Down
3 changes: 3 additions & 0 deletions DQM/Physics/src/CentralitypADQM.h
Original file line number Diff line number Diff line change
Expand Up @@ -73,6 +73,9 @@ class CentralitypADQM : public DQMEDAnalyzer {
MonitorElement* h_hiZDC;
MonitorElement* h_hiZDCplus;
MonitorElement* h_hiZDCminus;
MonitorElement* h_hiPF;
MonitorElement* h_hiPFplus;
MonitorElement* h_hiPFminus;

MonitorElement* h_vertex_x;
MonitorElement* h_vertex_y;
Expand Down
6 changes: 3 additions & 3 deletions DQMOffline/JetMET/src/JetAnalyzer_HeavyIons.cc
Original file line number Diff line number Diff line change
Expand Up @@ -642,16 +642,16 @@ void JetAnalyzer_HeavyIons::analyze(const edm::Event &mEvent, const edm::EventSe
edm::Handle<reco::Centrality> cent;
mEvent.getByToken(centralityToken, cent); //_centralitytag comes from the cfg

if (!cent.isValid())
return;

mHF->Fill(cent->EtHFtowerSum());
Float_t HF_energy = cent->EtHFtowerSum();

//for later when centrality gets added to RelVal
//edm::Handle<int> cbin;
//mEvent.getByToken(centralityBinToken, cbin);

if (!cent.isValid())
return;

/*int hibin = -999;
if(cbin.isValid()){
hibin = *cbin;
Expand Down
5 changes: 5 additions & 0 deletions DataFormats/HeavyIonEvent/interface/Centrality.h
Original file line number Diff line number Diff line change
Expand Up @@ -55,6 +55,9 @@ namespace reco {
double zdcSumPlus() const { return zdcSumPlus_; }
double zdcSumMinus() const { return zdcSumMinus_; }
double EtMidRapiditySum() const { return etMidRapiditySum_; }
double EtPFhfSum() const { return etPFhfSumPlus_ + etPFhfSumMinus_; }
double EtPFhfSumPlus() const { return etPFhfSumPlus_; }
double EtPFhfSumMinus() const { return etPFhfSumMinus_; }

protected:
double value_;
Expand Down Expand Up @@ -84,6 +87,8 @@ namespace reco {
double zdcSumPlus_;
double zdcSumMinus_;
double etMidRapiditySum_;
double etPFhfSumPlus_;
double etPFhfSumMinus_;
double ntracksPtCut_;
double ntracksEtaCut_;
double ntracksEtaPtCut_;
Expand Down
2 changes: 2 additions & 0 deletions DataFormats/HeavyIonEvent/src/Centrality.cc
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,8 @@ Centrality::Centrality(double d, std::string label)
zdcSumPlus_(0),
zdcSumMinus_(0),
etMidRapiditySum_(0),
etPFhfSumPlus_(0),
etPFhfSumMinus_(0),
ntracksPtCut_(0),
ntracksEtaCut_(0),
ntracksEtaPtCut_(0),
Expand Down
3 changes: 2 additions & 1 deletion DataFormats/HeavyIonEvent/src/classes_def.xml
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,8 @@
<version ClassVersion="11" checksum="393681411"/>
</class>
<class name="edm::Wrapper<reco::EvtPlane>"/>
<class name="reco::Centrality" ClassVersion="11">
<class name="reco::Centrality" ClassVersion="12">
<version ClassVersion="12" checksum="926006939"/>
<version ClassVersion="11" checksum="688182903"/>
<version ClassVersion="10" checksum="3770933457"/>
</class>
Expand Down
34 changes: 28 additions & 6 deletions RecoHI/HiCentralityAlgos/plugins/CentralityBinProducer.cc
Original file line number Diff line number Diff line change
Expand Up @@ -59,7 +59,10 @@ class CentralityBinProducer : public edm::stream::EDProducer<> {
EE = 11,
ZDChitsPlus = 12,
ZDChitsMinus = 13,
Missing = 14
PFhf = 14,
PFhfPlus = 15,
PFhfMinus = 16,
Missing = 17
};

public:
Expand Down Expand Up @@ -131,12 +134,18 @@ CentralityBinProducer::CentralityBinProducer(const edm::ParameterSet& iConfig) :
varType_ = ZDChitsPlus;
if (centralityVariable_ == "ZDChitsMinus")
varType_ = ZDChitsMinus;
if (centralityVariable_ == "PFhf")
varType_ = PFhf;
if (centralityVariable_ == "PFhfPlus")
varType_ = PFhfPlus;
if (centralityVariable_ == "PFhfMinus")
varType_ = PFhfMinus;
if (varType_ == Missing) {
std::string errorMessage = "Requested Centrality variable does not exist : " + centralityVariable_ + "\n" +
"Supported variables are: \n" +
"HFtowers HFtowersPlus HFtowersMinus HFtowersTrunc HFtowersPlusTrunc HFtowersMinusTrunc "
"HFhits PixelHits PixelTracks Tracks EB EE" +
"\n";
std::string errorMessage =
"Requested Centrality variable does not exist : " + centralityVariable_ + "\n" + "Supported variables are: \n" +
"HFtowers HFtowersPlus HFtowersMinus HFtowersTrunc HFtowersPlusTrunc HFtowersMinusTrunc "
"HFhits PixelHits PixelTracks Tracks EB EE ZDChitsPlus ZDChitsMinus PFhf PFhfPlus PFhfMinus" +
"\n";
throw cms::Exception("Configuration", errorMessage);
}

Expand Down Expand Up @@ -207,6 +216,15 @@ void CentralityBinProducer::produce(edm::Event& iEvent, const edm::EventSetup& i
case ZDChitsMinus:
value = chandle_->zdcSumMinus();
break;
case PFhf:
value = chandle_->EtPFhfSum();
break;
case PFhfPlus:
value = chandle_->EtPFhfSumPlus();
break;
case PFhfMinus:
value = chandle_->EtPFhfSumMinus();
break;
default:
throw cms::Exception("CentralityBinProducer", "Centrality variable not recognized.");
}
Expand Down Expand Up @@ -236,6 +254,10 @@ void CentralityBinProducer::beginRun(edm::Run const& iRun, const edm::EventSetup
varType_ = ZDChitsMinus;
if (centralityVariable_ == "ZDChitsMinus")
varType_ = ZDChitsPlus;
if (centralityVariable_ == "PFhfPlus")
varType_ = PFhfMinus;
if (centralityVariable_ == "PFhfMinus")
varType_ = PFhfPlus;
}
prevRun_ = iRun.run();

Expand Down
23 changes: 23 additions & 0 deletions RecoHI/HiCentralityAlgos/plugins/CentralityProducer.cc
Original file line number Diff line number Diff line change
Expand Up @@ -73,6 +73,7 @@ namespace reco {
const bool doPixelCut_;
const bool produceTracks_;
const bool producePixelTracks_;
const bool producePF_;

const double midRapidityRange_;
const double trackPtCut_;
Expand All @@ -91,6 +92,7 @@ namespace reco {
const edm::EDGetTokenT<SiPixelRecHitCollection> srcPixelhits_;
const edm::EDGetTokenT<TrackCollection> srcTracks_;
const edm::EDGetTokenT<TrackCollection> srcPixelTracks_;
const edm::EDGetTokenT<reco::CandidateView> srcPF_;
const edm::EDGetTokenT<VertexCollection> srcVertex_;
const edm::EDGetTokenT<Centrality> reuseTag_;

Expand Down Expand Up @@ -121,6 +123,7 @@ namespace reco {
doPixelCut_(iConfig.getParameter<bool>("doPixelCut")),
produceTracks_(iConfig.getParameter<bool>("produceTracks")),
producePixelTracks_(iConfig.getParameter<bool>("producePixelTracks")),
producePF_(iConfig.getParameter<bool>("producePF")),
midRapidityRange_(iConfig.getParameter<double>("midRapidityRange")),
trackPtCut_(iConfig.getParameter<double>("trackPtCut")),
trackEtaCut_(iConfig.getParameter<double>("trackEtaCut")),
Expand All @@ -147,6 +150,8 @@ namespace reco {
srcPixelTracks_(producePixelTracks_
? consumes<TrackCollection>(iConfig.getParameter<edm::InputTag>("srcPixelTracks"))
: edm::EDGetTokenT<TrackCollection>()),
srcPF_(producePF_ ? consumes<reco::CandidateView>(iConfig.getParameter<edm::InputTag>("srcPF"))
: edm::EDGetTokenT<reco::CandidateView>()),
srcVertex_((produceTracks_ || producePixelTracks_)
? consumes<VertexCollection>(iConfig.getParameter<edm::InputTag>("srcVertex"))
: edm::EDGetTokenT<VertexCollection>()),
Expand Down Expand Up @@ -245,6 +250,22 @@ namespace reco {
}
}

if (producePF_) {
creco->etPFhfSumPlus_ = 0;
creco->etPFhfSumMinus_ = 0;
for (const auto& pf : iEvent.get(srcPF_)) {
stahlleiton marked this conversation as resolved.
Show resolved Hide resolved
if (pf.pdgId() != 1 && pf.pdgId() != 2)
continue;
if (pf.eta() > 0)
creco->etPFhfSumPlus_ += pf.pt();
else
creco->etPFhfSumMinus_ += pf.pt();
}
} else if (reuseAny_) {
creco->etPFhfSumMinus_ = inputCentrality->EtPFhfSumMinus();
creco->etPFhfSumPlus_ = inputCentrality->EtPFhfSumPlus();
}

if (produceEcalhits_) {
creco->etEESumPlus_ = 0;
creco->etEESumMinus_ = 0;
Expand Down Expand Up @@ -497,6 +518,7 @@ namespace reco {
desc.add<bool>("producePixelhits", true);
desc.add<bool>("produceTracks", true);
desc.add<bool>("producePixelTracks", true);
desc.add<bool>("producePF", true);
desc.add<bool>("reUseCentrality", false);
desc.add<edm::InputTag>("srcHFhits", edm::InputTag("hfreco"));
desc.add<edm::InputTag>("srcTowers", edm::InputTag("towerMaker"));
Expand All @@ -508,6 +530,7 @@ namespace reco {
desc.add<edm::InputTag>("srcVertex", edm::InputTag("hiSelectedVertex"));
desc.add<edm::InputTag>("srcReUse", edm::InputTag("hiCentrality"));
desc.add<edm::InputTag>("srcPixelTracks", edm::InputTag("hiPixel3PrimTracks"));
desc.add<edm::InputTag>("srcPF", edm::InputTag("particleFlow"));
desc.add<bool>("doPixelCut", true);
desc.add<bool>("useQuality", true);
desc.add<string>("trackQuality", "highPurity");
Expand Down
2 changes: 2 additions & 0 deletions RecoHI/HiCentralityAlgos/python/HiCentrality_cfi.py
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@
producePixelhits = cms.bool(True),
produceTracks = cms.bool(True),
producePixelTracks = cms.bool(True),
producePF = cms.bool(True),
reUseCentrality = cms.bool(False),

srcHFhits = cms.InputTag("hfreco"),
Expand All @@ -22,6 +23,7 @@
srcVertex= cms.InputTag("hiSelectedVertex"),
srcReUse = cms.InputTag("hiCentrality"),
srcPixelTracks = cms.InputTag("hiPixel3PrimTracks"),
srcPF = cms.InputTag("particleFlow"),

doPixelCut = cms.bool(True),
useQuality = cms.bool(True),
Expand Down
2 changes: 2 additions & 0 deletions RecoHI/HiCentralityAlgos/python/pACentrality_cfi.py
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@
producePixelhits = cms.bool(True),
produceTracks = cms.bool(True),
producePixelTracks = cms.bool(True),
producePF = cms.bool(True),
reUseCentrality = cms.bool(False),

srcHFhits = cms.InputTag("hfreco"),
Expand All @@ -22,6 +23,7 @@
srcVertex= cms.InputTag("offlinePrimaryVertices"),
srcReUse = cms.InputTag("pACentrality"),
srcPixelTracks = cms.InputTag("pixelTracks"),
srcPF = cms.InputTag("particleFlow"),

doPixelCut = cms.bool(True),
useQuality = cms.bool(True),
Expand Down
31 changes: 15 additions & 16 deletions RecoHI/HiJetAlgos/plugins/HiFJGridEmptyAreaCalculator.cc
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,7 @@ HiFJGridEmptyAreaCalculator::HiFJGridEmptyAreaCalculator(const edm::ParameterSet
ntotalJet_ = 0;
nyJet_ = 0;

pfCandsToken_ = consumes<reco::PFCandidateCollection>(iConfig.getParameter<edm::InputTag>("pfCandSource"));
pfCandsToken_ = consumes<reco::CandidateView>(iConfig.getParameter<edm::InputTag>("pfCandSource"));
mapEtaToken_ = consumes<std::vector<double>>(iConfig.getParameter<edm::InputTag>("mapEtaEdges"));
mapRhoToken_ = consumes<std::vector<double>>(iConfig.getParameter<edm::InputTag>("mapToRho"));
mapRhoMToken_ = consumes<std::vector<double>>(iConfig.getParameter<edm::InputTag>("mapToRhoM"));
Expand Down Expand Up @@ -171,11 +171,10 @@ void HiFJGridEmptyAreaCalculator::produce(edm::Event& iEvent, const edm::EventSe
void HiFJGridEmptyAreaCalculator::calculateGridRho(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
vector<vector<double>> scalarPt(ny_, vector<double>(nphi_, 0.0));

edm::Handle<reco::PFCandidateCollection> pfCands;
edm::Handle<reco::CandidateView> pfCands;
iEvent.getByToken(pfCandsToken_, pfCands);
const reco::PFCandidateCollection* pfCandidateColl = pfCands.product();
for (unsigned icand = 0; icand < pfCandidateColl->size(); icand++) {
const reco::PFCandidate pfCandidate = pfCandidateColl->at(icand);
const auto& pfCandidateColl = pfCands.product();
for (const auto& pfCandidate : *pfCandidateColl) {
//use ony the particles within the eta range
if (pfCandidate.eta() < ymin_ || pfCandidate.eta() > ymax_)
continue;
Expand Down Expand Up @@ -231,17 +230,17 @@ void HiFJGridEmptyAreaCalculator::calculateAreaFractionOfJets(const edm::Event&
//calculate jet kt area fraction inside boundary by grid
totalInboundArea_ = 0;

for (auto jet = jets->begin(); jet != jets->end(); ++jet) {
if (jet->eta() < etaminJet_ || jet->eta() > etamaxJet_)
for (const auto& jet : *jets) {
if (jet.eta() < etaminJet_ || jet.eta() > etamaxJet_)
continue;

double areaKt = jet->jetArea();
setupGridJet(&*jet);
double areaKt = jet.jetArea();
setupGridJet(&jet);
std::vector<std::pair<int, int>> pfIndicesJet;
std::vector<std::pair<int, int>> pfIndicesJetInbound;
int nConstitJet = 0;
int nConstitJetInbound = 0;
for (auto daughter : jet->getJetConstituentsQuick()) {
for (const auto& daughter : jet.getJetConstituentsQuick()) {
auto pfCandidate = static_cast<const reco::PFCandidate*>(daughter);

int jeta = tileIndexEtaJet(&*pfCandidate);
Expand Down Expand Up @@ -323,7 +322,7 @@ void HiFJGridEmptyAreaCalculator::setupGrid(double etamin, double etamax) {

//----------------------------------------------------------------------
// retrieve the grid tile index for a given PseudoJet
int HiFJGridEmptyAreaCalculator::tileIndexPhi(const reco::PFCandidate* pfCand) {
int HiFJGridEmptyAreaCalculator::tileIndexPhi(const reco::Candidate* pfCand) {
// directly taking int does not work for values between -1 and 0
// so use floor instead
// double iy_double = (p.rap() - _ymin) / _dy;
Expand All @@ -345,7 +344,7 @@ int HiFJGridEmptyAreaCalculator::tileIndexPhi(const reco::PFCandidate* pfCand) {

//----------------------------------------------------------------------
// retrieve the grid tile index for a given PseudoJet
int HiFJGridEmptyAreaCalculator::tileIndexEta(const reco::PFCandidate* pfCand) {
int HiFJGridEmptyAreaCalculator::tileIndexEta(const reco::Candidate* pfCand) {
// directly taking int does not work for values between -1 and 0
// so use floor instead
// double iy_double = (p.rap() - _ymin) / _dy;
Expand Down Expand Up @@ -393,7 +392,7 @@ void HiFJGridEmptyAreaCalculator::setupGridJet(const reco::Jet* jet) {

//----------------------------------------------------------------------
// retrieve the grid tile index for a given PseudoJet
int HiFJGridEmptyAreaCalculator::tileIndexEtaJet(const reco::PFCandidate* pfCand) {
int HiFJGridEmptyAreaCalculator::tileIndexEtaJet(const reco::Candidate* pfCand) {
// directly taking int does not work for values between -1 and 0
// so use floor instead
// double iy_double = (p.rap() - _ymin) / _dy;
Expand Down Expand Up @@ -422,9 +421,9 @@ int HiFJGridEmptyAreaCalculator::numJetGridCells(std::vector<std::pair<int, int>
int highestJEta = lowestJEta;

//for each fixed phi value calculate the number of grids in eta
for (unsigned int iconst = 1; iconst < indices.size(); iconst++) {
int jphi = indices[iconst].first;
int jeta = indices[iconst].second;
for (const auto& index : indices) {
int jphi = index.first;
int jeta = index.second;
if (jphi == lowestJPhi) {
if (jeta < lowestJEta)
lowestJEta = jeta;
Expand Down
10 changes: 5 additions & 5 deletions RecoHI/HiJetAlgos/plugins/HiFJGridEmptyAreaCalculator.h
Original file line number Diff line number Diff line change
Expand Up @@ -46,10 +46,10 @@ class HiFJGridEmptyAreaCalculator : public edm::stream::EDProducer<> {
void setupGridJet(const reco::Jet* jet);

/// retrieve the grid cell index for a given PseudoJet
int tileIndexJet(const reco::PFCandidate* pfCand);
int tileIndexEta(const reco::PFCandidate* pfCand);
int tileIndexEtaJet(const reco::PFCandidate* pfCand);
int tileIndexPhi(const reco::PFCandidate* pfCand);
int tileIndexJet(const reco::Candidate* pfCand);
int tileIndexEta(const reco::Candidate* pfCand);
int tileIndexEtaJet(const reco::Candidate* pfCand);
int tileIndexPhi(const reco::Candidate* pfCand);

///number of grid cells that overlap with jet constituents filling in the in between area
int numJetGridCells(std::vector<std::pair<int, int>>& indices);
Expand Down Expand Up @@ -106,7 +106,7 @@ class HiFJGridEmptyAreaCalculator : public edm::stream::EDProducer<> {

/// input tokens
edm::EDGetTokenT<edm::View<reco::Jet>> jetsToken_;
edm::EDGetTokenT<reco::PFCandidateCollection> pfCandsToken_;
edm::EDGetTokenT<reco::CandidateView> pfCandsToken_;
edm::EDGetTokenT<std::vector<double>> mapEtaToken_;
edm::EDGetTokenT<std::vector<double>> mapRhoToken_;
edm::EDGetTokenT<std::vector<double>> mapRhoMToken_;
Expand Down