Skip to content

Commit

Permalink
Merge pull request #44518 from stahlleiton/PFHFCentrality_CMSSW_14_1_X
Browse files Browse the repository at this point in the history 
Add PF HF sum energy to reco::Centrality
  • Loading branch information
@cmsbuild
cmsbuild committed Apr 3, 2024
2 parents f34cb55 + b328705 commit 6cda5fd
Show file tree
Hide file tree
Showing 14 changed files with 107 additions and 31 deletions.
7 changes: 7 additions & 0 deletions DQM/Physics/src/CentralityDQM.cc
View file
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
View file
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
View file
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
View file
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
View file
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
View file
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
View file
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
View file
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
View file
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
View file
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_)) {
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
View file
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
View file
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
View file
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
View file
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

0 comments on commit 6cda5fd

Please sign in to comment.