From 86b6c2652ceb7ea2896c5c4774d7150e8a8e4c40 Mon Sep 17 00:00:00 2001 From: Martin Date: Wed, 13 Oct 2021 15:54:43 -0500 Subject: [PATCH 1/2] DataFormat for MuonRechitCluster --- .../MuonReco/interface/MuonRecHitCluster.h | 66 +++++ DataFormats/MuonReco/src/MuonRecHitCluster.cc | 26 ++ DataFormats/MuonReco/src/classes.h | 1 + DataFormats/MuonReco/src/classes_def.xml | 6 + .../plugins/BuildFile.xml | 9 + .../plugins/RechitClusterProducer.cc | 268 ++++++++++++++++++ .../python/rechitCluster_cff.py | 18 ++ 7 files changed, 394 insertions(+) create mode 100644 DataFormats/MuonReco/interface/MuonRecHitCluster.h create mode 100644 DataFormats/MuonReco/src/MuonRecHitCluster.cc create mode 100644 RecoMuon/MuonRechitClusterProducer/plugins/BuildFile.xml create mode 100644 RecoMuon/MuonRechitClusterProducer/plugins/RechitClusterProducer.cc create mode 100644 RecoMuon/MuonRechitClusterProducer/python/rechitCluster_cff.py diff --git a/DataFormats/MuonReco/interface/MuonRecHitCluster.h b/DataFormats/MuonReco/interface/MuonRecHitCluster.h new file mode 100644 index 0000000000000..6383cd1b4e9f3 --- /dev/null +++ b/DataFormats/MuonReco/interface/MuonRecHitCluster.h @@ -0,0 +1,66 @@ +#ifndef DataFormats_MuonReco_MuonRecHitCluster_h +#define DataFormats_MuonReco_MuonRecHitCluster_h + +#include +#include "DataFormats/Common/interface/SortedCollection.h" +#include "DataFormats/Math/interface/Vector3D.h" + +namespace reco { + + class MuonRecHitCluster { + public: + //default constructor + MuonRecHitCluster() = default; + + MuonRecHitCluster(const math::RhoEtaPhiVectorF position, + const int size, + const int nStation, + const float avgStation, + const float time, + const float timeSpread, + const int nME11, + const int nME12, + const int nME41, + const int nME42, + const int nMB1, + const int nMB2); + + // + ~MuonRecHitCluster() = default; + + float eta() const { return position_.Eta(); } + float phi() const { return position_.Phi(); } + float x() const { return position_.X(); } + float y() const { return position_.Y(); } + float z() const { return position_.Z(); } + float r() const { return position_.Rho(); } + int size() const { return size_; } + int nStation() const { return nStation_; } + float avgStation() const { return avgStation_; } + int nMB1() const { return nMB1_; } + int nMB2() const { return nMB2_; } + int nME11() const { return nME11_; } + int nME12() const { return nME12_; } + int nME41() const { return nME41_; } + int nME42() const { return nME42_; } + float time() const { return time_; } + float timeSpread() const { return timeSpread_; } + + private: + math::RhoEtaPhiVectorF position_; + int size_; + int nStation_; + float avgStation_; + float time_; + float timeSpread_; + int nME11_; + int nME12_; + int nME41_; + int nME42_; + int nMB1_; + int nMB2_; + }; + + typedef std::vector MuonRecHitClusterCollection; +} // namespace reco +#endif diff --git a/DataFormats/MuonReco/src/MuonRecHitCluster.cc b/DataFormats/MuonReco/src/MuonRecHitCluster.cc new file mode 100644 index 0000000000000..22e96f004e789 --- /dev/null +++ b/DataFormats/MuonReco/src/MuonRecHitCluster.cc @@ -0,0 +1,26 @@ +#include "DataFormats/MuonReco/interface/MuonRecHitCluster.h" + +reco::MuonRecHitCluster::MuonRecHitCluster(const math::RhoEtaPhiVectorF position, + const int size, + const int nStation, + const float avgStation, + const float time, + const float timeSpread, + const int nME11, + const int nME12, + const int nME41, + const int nME42, + const int nMB1, + const int nMB2) + : position_(position), + size_(size), + nStation_(nStation), + avgStation_(avgStation), + time_(time), + timeSpread_(timeSpread), + nME11_(nME11), + nME12_(nME12), + nME41_(nME41), + nME42_(nME42), + nMB1_(nMB1), + nMB2_(nMB2) {} diff --git a/DataFormats/MuonReco/src/classes.h b/DataFormats/MuonReco/src/classes.h index a9b8eb5c5bb17..bfbdebc1a5655 100644 --- a/DataFormats/MuonReco/src/classes.h +++ b/DataFormats/MuonReco/src/classes.h @@ -20,6 +20,7 @@ #include "DataFormats/MuonReco/interface/MuonQuality.h" #include "DataFormats/MuonReco/interface/MuonCosmicCompatibility.h" #include "DataFormats/MuonReco/interface/MuonShower.h" +#include "DataFormats/MuonReco/interface/MuonRecHitCluster.h" #include "DataFormats/MuonReco/interface/MuonToMuonMap.h" #include "DataFormats/TrackReco/interface/Track.h" #include "DataFormats/Common/interface/AssociationMap.h" diff --git a/DataFormats/MuonReco/src/classes_def.xml b/DataFormats/MuonReco/src/classes_def.xml index 1184cdbbd1e46..77e7419873b8a 100644 --- a/DataFormats/MuonReco/src/classes_def.xml +++ b/DataFormats/MuonReco/src/classes_def.xml @@ -164,6 +164,12 @@ initial version number of a class which has never been stored before. + + + + + + diff --git a/RecoMuon/MuonRechitClusterProducer/plugins/BuildFile.xml b/RecoMuon/MuonRechitClusterProducer/plugins/BuildFile.xml new file mode 100644 index 0000000000000..8e4d6eaea8968 --- /dev/null +++ b/RecoMuon/MuonRechitClusterProducer/plugins/BuildFile.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/RecoMuon/MuonRechitClusterProducer/plugins/RechitClusterProducer.cc b/RecoMuon/MuonRechitClusterProducer/plugins/RechitClusterProducer.cc new file mode 100644 index 0000000000000..5751491aa67b5 --- /dev/null +++ b/RecoMuon/MuonRechitClusterProducer/plugins/RechitClusterProducer.cc @@ -0,0 +1,268 @@ +#include +#include + +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/global/EDProducer.h" + +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Utilities/interface/StreamID.h" + +#include "DataFormats/Common/interface/Ptr.h" +#include "DataFormats/Math/interface/deltaPhi.h" +#include "DataFormats/MuonReco/interface/MuonRecHitCluster.h" +#include "Geometry/CSCGeometry/interface/CSCGeometry.h" +#include "Geometry/Records/interface/MuonGeometryRecord.h" +#include "DataFormats/MuonDetId/interface/CSCDetId.h" + +#include "fastjet/JetDefinition.hh" +#include "fastjet/ClusterSequence.hh" +#include "fastjet/PseudoJet.hh" + +#include "Geometry/CSCGeometry/interface/CSCGeometry.h" +#include "Geometry/DTGeometry/interface/DTGeometry.h" +#include "DataFormats/MuonDetId/interface/CSCDetId.h" +#include "DataFormats/Math/interface/Vector3D.h" +#include "DataFormats/Math/interface/PtEtaPhiMass.h" +#include "DataFormats/DTRecHit/interface/DTRecHitCollection.h" +#include "DataFormats/CSCRecHit/interface/CSCRecHit2DCollection.h" +#include "DataFormats/MuonReco/interface/MuonRecHitCluster.h" + +typedef std::vector RecHitClusterCollection; + +template +class RechitClusterProducerT : public edm::global::EDProducer<> { + typedef typename Trait::RecHitRef RechitRef; + typedef typename Trait::RecHitRefVector RecHitRefVector; + +public: + explicit RechitClusterProducerT(const edm::ParameterSet&); + ~RechitClusterProducerT() override = default; + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + +private: + void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override; + const edm::ESGetToken GeometryToken_; + edm::EDGetTokenT InputToken_; + typedef std::vector RecHitClusterCollection; + + const double rParam_; // distance paramter + const int nRechitMin_; // min number of rechits + const int nStationThres_; // min number of rechits to count towards nStation +}; + +template +RechitClusterProducerT::RechitClusterProducerT(const edm::ParameterSet& iConfig) + : GeometryToken_(esConsumes()), + InputToken_(consumes(iConfig.getParameter("recHitLabel"))), + rParam_(iConfig.getParameter("rParam")), + nRechitMin_(iConfig.getParameter("nRechitMin")), + nStationThres_(iConfig.getParameter("nStationThres")) { + produces(); +} + +// +// member functions +// + +// ------------ method called to produce the data ------------ +template +void RechitClusterProducerT::produce(edm::StreamID, edm::Event& ev, const edm::EventSetup& iSetup) const { + auto const& geo = iSetup.getData(GeometryToken_); + auto const& rechits = ev.get(InputToken_); + + std::vector fjInput; + RecHitRefVector inputs_; + + fastjet::JetDefinition jet_def(fastjet::cambridge_algorithm, rParam_); + + inputs_.clear(); + fjInput.clear(); + + Trait RecHitTrait; + + int recIt = 0; + for (auto const& rechit : rechits) { + LocalPoint RecHitLocalPosition = rechit.localPosition(); + auto detid = RecHitTrait.detid(rechit); + auto thischamber = geo.chamber(detid); + if (thischamber) { + GlobalPoint globalPosition = thischamber->toGlobal(RecHitLocalPosition); + float x = globalPosition.x(); + float y = globalPosition.y(); + float z = globalPosition.z(); + RechitRef ref = RechitRef(&rechits, recIt); + inputs_.push_back(ref); + fjInput.push_back(fastjet::PseudoJet(x, y, z, globalPosition.mag())); + fjInput.back().set_user_index(recIt); + } + recIt++; + } + fastjet::ClusterSequence clus_seq(fjInput, jet_def); + + //keep all the clusters + double ptmin = 0.0; + std::vector fjJets = fastjet::sorted_by_pt(clus_seq.inclusive_jets(ptmin)); + + auto clusters = std::make_unique(); + if (!fjJets.empty()) { + for (unsigned int ijet = 0; ijet < fjJets.size(); ++ijet) { + // get the fastjet jet + const fastjet::PseudoJet& fjJet = fjJets[ijet]; + + // skip if the cluster has too few rechits + if (int(fjJet.constituents().size()) < nRechitMin_) + continue; + // get the constituents from fastjet + RecHitRefVector rechits; + for (unsigned int i = 0; i < fjJet.constituents().size(); i++) { + auto index = fjJet.constituents()[i].user_index(); + if (index >= 0 && static_cast(index) < inputs_.size()) { + rechits.push_back(inputs_[index]); + } + } + + //Derive cluster properties + int nStation = 0; + int totStation = 0; + float avgStation = 0.0; + std::map station_count_map; + for (auto& rechit : rechits) { + station_count_map[RecHitTrait.station(*rechit)]++; + } + //station statistics + std::map::iterator it; + for (auto const& [station, count] : station_count_map) { + if (count >= nStationThres_) { + nStation++; + avgStation += (station) * (count); + totStation += (count); + } + } + if (totStation != 0) { + avgStation = avgStation / totStation; + } + float invN = 1.f / rechits.size(); + // cluster position is the average position of the constituent rechits + float jetX = fjJet.px() * invN; + float jetY = fjJet.py() * invN; + float jetZ = fjJet.pz() * invN; + + math::RhoEtaPhiVectorF position( + std::sqrt(jetX * jetX + jetY * jetY), etaFromXYZ(jetX, jetY, jetZ), std::atan2(jetY, jetX)); + Trait::emplace_back(clusters.get(), position, nStation, avgStation, rechits); + } + } + ev.put(std::move(clusters)); +} + +// ------------ method fills 'descriptions' with the allowed parameters for the module ------------ +template +void RechitClusterProducerT::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("nRechitMin", 50); + desc.add("rParam", 0.4); + desc.add("nStationThres", 10); + desc.add("recHitLabel", edm::InputTag(Trait::recHitLabel())); + descriptions.add(Trait::producerName(), desc); +} +struct DTRecHitTrait { + using GeometryType = DTGeometry; + using InputType = DTRecHitCollection; + using RecHitRef = edm::Ref; + using RecHitRefVector = edm::RefVector; + static std::string recHitLabel() { return "dt1DRecHits"; } + static std::string producerName() { return "DTrechitClusterProducer"; } + + static int station(const DTRecHit1DPair& dtRechit) { return detid(dtRechit).station(); } + static DTChamberId detid(const DTRecHit1DPair& dtRechit) { + DetId geoid = dtRechit.geographicalId(); + DTChamberId dtdetid = DTChamberId(geoid); + return dtdetid; + } + static void emplace_back(RecHitClusterCollection* clusters, + math::RhoEtaPhiVectorF const& position, + int nStation, + float avgStation, + RecHitRefVector const& rechits) { + // compute nMB1, nMB2 + int nMB1 = 0; + int nMB2 = 0; + for (auto& rechit : rechits) { + DetId geoid = rechit->geographicalId(); + DTChamberId dtdetid = DTChamberId(geoid); + if (dtdetid.station() == 1) + nMB1++; + if (dtdetid.station() == 2) + nMB2++; + } + //set time, timespread, nME11,nME12 to 0 + reco::MuonRecHitCluster cls(position, rechits.size(), nStation, avgStation, 0.0, 0.0, 0, 0, 0, 0, nMB1, nMB2); + clusters->emplace_back(cls); + } +}; + +struct CSCRecHitTrait { + using GeometryType = CSCGeometry; + using InputType = CSCRecHit2DCollection; + using RecHitRef = edm::Ref; + using RecHitRefVector = edm::RefVector; + static std::string recHitLabel() { return "csc2DRecHits"; } + static std::string producerName() { return "CSCrechitClusterProducer"; } + + static int station(const CSCRecHit2D& cscRechit) { return CSCDetId::station(detid(cscRechit)); } + static CSCDetId detid(const CSCRecHit2D& cscRechit) { return cscRechit.cscDetId(); } + static void emplace_back(RecHitClusterCollection* clusters, + math::RhoEtaPhiVectorF const& position, + int nStation, + float avgStation, + RecHitRefVector const& rechits) { + int nME11 = 0; + int nME12 = 0; + int nME41 = 0; + int nME42 = 0; + float timeSpread = 0.0; + float time = 0.0; + float time_strip = 0.0; // for timeSpread calculation + for (auto& rechit : rechits) { + CSCDetId cscdetid = rechit->cscDetId(); + int stationRing = (CSCDetId::station(cscdetid) * 10 + CSCDetId::ring(cscdetid)); + if (CSCDetId::ring(cscdetid) == 4) + stationRing = (CSCDetId::station(cscdetid) * 10 + 1); // ME1/a has ring==4 + if (stationRing == 11) + nME11++; + if (stationRing == 12) + nME12++; + if (stationRing == 41) + nME41++; + if (stationRing == 42) + nME42++; + time += (rechit->tpeak() + rechit->wireTime()); + time_strip += rechit->tpeak(); + } + float invN = 1.f / rechits.size(); + time = (time / 2.f) * invN; + time_strip = time_strip * invN; + + //derive cluster statistics + for (auto& rechit : rechits) { + timeSpread += (rechit->tpeak() - time_strip) * (rechit->tpeak() - time_strip); + } + timeSpread = std::sqrt(timeSpread * invN); + + //set nMB1,nMB2 to 0 + reco::MuonRecHitCluster cls( + position, rechits.size(), nStation, avgStation, time, timeSpread, nME11, nME12, nME41, nME42, 0, 0); + clusters->emplace_back(cls); + } +}; + +using DTrechitClusterProducer = RechitClusterProducerT; +using CSCrechitClusterProducer = RechitClusterProducerT; +DEFINE_FWK_MODULE(DTrechitClusterProducer); +DEFINE_FWK_MODULE(CSCrechitClusterProducer); diff --git a/RecoMuon/MuonRechitClusterProducer/python/rechitCluster_cff.py b/RecoMuon/MuonRechitClusterProducer/python/rechitCluster_cff.py new file mode 100644 index 0000000000000..f4dc69e7f77c8 --- /dev/null +++ b/RecoMuon/MuonRechitClusterProducer/python/rechitCluster_cff.py @@ -0,0 +1,18 @@ +import FWCore.ParameterSet.Config as cms + +import RecoMuon.MuonRechitClusterProducer.CSCrechitClusterProducer_cfi as CSCcluster +import RecoMuon.MuonRechitClusterProducer.DTrechitClusterProducer_cfi as DTcluster + +ca4CSCrechitClusters= CSCcluster.CSCrechitClusterProducer.clone( + recHitLabel = "csc2DRecHits", + nRechitMin = 50, + rParam = 0.4, + nStationThres = 10, +) +ca4DTrechitClusters = DTcluster.DTrechitClusterProducer.clone( + recHitLabel = "dt1DRecHits", + nRechitMin = 50, + rParam = 0.4, + nStationThres = 10, +) + From 5b73c3439ce20091997914fc4e2ef0a80816272f Mon Sep 17 00:00:00 2001 From: Martin Date: Tue, 23 Nov 2021 14:55:24 -0600 Subject: [PATCH 2/2] add HLT filter --- .../MuonReco/interface/MuonRecHitCluster.h | 1 - HLTrigger/special/plugins/BuildFile.xml | 1 + .../plugins/HLTMuonRecHitClusterFilter.cc | 106 +++++++++++ .../plugins/RechitClusterProducer.cc | 173 +++++++++--------- .../python/rechitCluster_cff.py | 8 +- 5 files changed, 197 insertions(+), 92 deletions(-) create mode 100644 HLTrigger/special/plugins/HLTMuonRecHitClusterFilter.cc diff --git a/DataFormats/MuonReco/interface/MuonRecHitCluster.h b/DataFormats/MuonReco/interface/MuonRecHitCluster.h index 6383cd1b4e9f3..a792e11cd79a9 100644 --- a/DataFormats/MuonReco/interface/MuonRecHitCluster.h +++ b/DataFormats/MuonReco/interface/MuonRecHitCluster.h @@ -2,7 +2,6 @@ #define DataFormats_MuonReco_MuonRecHitCluster_h #include -#include "DataFormats/Common/interface/SortedCollection.h" #include "DataFormats/Math/interface/Vector3D.h" namespace reco { diff --git a/HLTrigger/special/plugins/BuildFile.xml b/HLTrigger/special/plugins/BuildFile.xml index 335f4c61ac7bf..8ddc513c461f2 100644 --- a/HLTrigger/special/plugins/BuildFile.xml +++ b/HLTrigger/special/plugins/BuildFile.xml @@ -25,6 +25,7 @@ + diff --git a/HLTrigger/special/plugins/HLTMuonRecHitClusterFilter.cc b/HLTrigger/special/plugins/HLTMuonRecHitClusterFilter.cc new file mode 100644 index 0000000000000..86f0c0a8db3f5 --- /dev/null +++ b/HLTrigger/special/plugins/HLTMuonRecHitClusterFilter.cc @@ -0,0 +1,106 @@ +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/Framework/interface/global/EDFilter.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" +#include "FWCore/Utilities/interface/EDGetToken.h" +#include "FWCore/Utilities/interface/InputTag.h" +#include "FWCore/Utilities/interface/StreamID.h" +#include "DataFormats/MuonReco/interface/MuonRecHitCluster.h" + +class HLTMuonRecHitClusterFilter : public edm::global::EDFilter<> { +public: + explicit HLTMuonRecHitClusterFilter(const edm::ParameterSet&); + ~HLTMuonRecHitClusterFilter() override = default; + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + +private: + bool filter(edm::StreamID, edm::Event&, edm::EventSetup const&) const override; + const edm::EDGetTokenT cluster_token_; + const int min_N_; + const int min_Size_; + const int min_SizeMinusMB1_; + const int max_nMB1_; + const int max_nMB2_; + const int max_nME11_; + const int max_nME12_; + const int max_nME41_; + const int max_nME42_; + const int min_Nstation_; + const double min_AvgStation_; + const double min_Time_; + const double max_Time_; + const double min_Eta_; + const double max_Eta_; + const double max_TimeSpread_; +}; +// +// constructors and destructor +// +HLTMuonRecHitClusterFilter::HLTMuonRecHitClusterFilter(const edm::ParameterSet& iConfig) + : cluster_token_(consumes(iConfig.getParameter("ClusterTag"))), + min_N_(iConfig.getParameter("MinN")), + min_Size_(iConfig.getParameter("MinSize")), + min_SizeMinusMB1_(iConfig.getParameter("MinSizeMinusMB1")), + max_nMB1_(iConfig.getParameter("Max_nMB1")), + max_nMB2_(iConfig.getParameter("Max_nMB2")), + max_nME11_(iConfig.getParameter("Max_nME11")), + max_nME12_(iConfig.getParameter("Max_nME12")), + max_nME41_(iConfig.getParameter("Max_nME41")), + max_nME42_(iConfig.getParameter("Max_nME42")), + min_Nstation_(iConfig.getParameter("MinNstation")), + min_AvgStation_(iConfig.getParameter("MinAvgStation")), + min_Time_(iConfig.getParameter("MinTime")), + max_Time_(iConfig.getParameter("MaxTime")), + min_Eta_(iConfig.getParameter("MinEta")), + max_Eta_(iConfig.getParameter("MaxEta")), + max_TimeSpread_(iConfig.getParameter("MaxTimeSpread")) {} + +void HLTMuonRecHitClusterFilter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("ClusterTag", edm::InputTag("hltCSCrechitClusters")); + desc.add("MinN", 1); + desc.add("MinSize", 50); + desc.add("MinSizeMinusMB1", 0); + desc.add("Max_nMB1", 0); + desc.add("Max_nMB2", 0); + desc.add("Max_nME11", 0); + desc.add("Max_nME12", 0); + desc.add("Max_nME41", 0); + desc.add("Max_nME42", 0); + desc.add("MinNstation", 0); + desc.add("MinAvgStation", 0.0); + desc.add("MinTime", -999); + desc.add("MaxTime", 999); + desc.add("MinEta", -1.0); + desc.add("MaxEta", -1.0); + desc.add("MaxTimeSpread", 999); + descriptions.addWithDefaultLabel(desc); +} + +// +// member functions +// + +// ------------ method called on each new Event ------------ +bool HLTMuonRecHitClusterFilter::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const { + int nClusterPassed = 0; + + auto const& rechitClusters = iEvent.get(cluster_token_); + + for (auto const& cluster : rechitClusters) { + if ((cluster.size() >= min_Size_) && ((cluster.size() - cluster.nMB1()) >= min_SizeMinusMB1_) && + (cluster.nMB1() <= max_nMB1_) && (cluster.nMB2() <= max_nMB2_) && (cluster.nME11() <= max_nME11_) && + (cluster.nME12() <= max_nME12_) && (cluster.nME41() <= max_nME41_) && (cluster.nME42() <= max_nME42_) && + (cluster.nStation() >= min_Nstation_) && (cluster.avgStation() >= min_AvgStation_) && + ((min_Eta_ < 0.0) || (std::abs(cluster.eta()) >= min_Eta_)) && + ((max_Eta_ < 0.0) || (std::abs(cluster.eta()) <= max_Eta_)) && (cluster.time() > min_Time_) && + (cluster.time() <= max_Time_) && (cluster.timeSpread() <= max_TimeSpread_)) { + nClusterPassed++; + } + } + + return (nClusterPassed >= min_N_); +} + +DEFINE_FWK_MODULE(HLTMuonRecHitClusterFilter); diff --git a/RecoMuon/MuonRechitClusterProducer/plugins/RechitClusterProducer.cc b/RecoMuon/MuonRechitClusterProducer/plugins/RechitClusterProducer.cc index 5751491aa67b5..2b06cbda6d25b 100644 --- a/RecoMuon/MuonRechitClusterProducer/plugins/RechitClusterProducer.cc +++ b/RecoMuon/MuonRechitClusterProducer/plugins/RechitClusterProducer.cc @@ -1,37 +1,44 @@ +#include #include -#include - -#include "FWCore/Framework/interface/Frameworkfwd.h" -#include "FWCore/Framework/interface/global/EDProducer.h" - +#include +#include +#include +#include "DataFormats/CSCRecHit/interface/CSCRecHit2D.h" +#include "DataFormats/CSCRecHit/interface/CSCRecHit2DCollection.h" +#include "DataFormats/Common/interface/RangeMap.h" +#include "DataFormats/Common/interface/Ref.h" +#include "DataFormats/Common/interface/RefVector.h" +#include "DataFormats/Common/interface/RefVectorIterator.h" +#include "DataFormats/DTRecHit/interface/DTRecHit1DPair.h" +#include "DataFormats/DTRecHit/interface/DTRecHitCollection.h" +#include "DataFormats/DetId/interface/DetId.h" +#include "DataFormats/GeometryVector/interface/GlobalPoint.h" +#include "DataFormats/GeometryVector/interface/LocalPoint.h" +#include "DataFormats/Math/interface/PtEtaPhiMass.h" +#include "DataFormats/Math/interface/Vector3D.h" +#include "DataFormats/MuonDetId/interface/CSCDetId.h" +#include "DataFormats/MuonDetId/interface/DTChamberId.h" +#include "DataFormats/MuonReco/interface/MuonRecHitCluster.h" +#include "FWCore/Framework/interface/ESHandle.h" #include "FWCore/Framework/interface/Event.h" -#include "FWCore/Framework/interface/EventSetup.h" #include "FWCore/Framework/interface/MakerMacros.h" -#include "FWCore/MessageLogger/interface/MessageLogger.h" - +#include "FWCore/Framework/interface/global/EDProducer.h" +#include "FWCore/Framework/interface/maker/WorkerMaker.h" #include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescriptionFiller.h" +#include "FWCore/Utilities/interface/EDGetToken.h" +#include "FWCore/Utilities/interface/ESGetToken.h" +#include "FWCore/Utilities/interface/InputTag.h" #include "FWCore/Utilities/interface/StreamID.h" - -#include "DataFormats/Common/interface/Ptr.h" -#include "DataFormats/Math/interface/deltaPhi.h" -#include "DataFormats/MuonReco/interface/MuonRecHitCluster.h" +#include "FWCore/Utilities/interface/typelookup.h" #include "Geometry/CSCGeometry/interface/CSCGeometry.h" +#include "Geometry/DTGeometry/interface/DTGeometry.h" #include "Geometry/Records/interface/MuonGeometryRecord.h" -#include "DataFormats/MuonDetId/interface/CSCDetId.h" - -#include "fastjet/JetDefinition.hh" #include "fastjet/ClusterSequence.hh" +#include "fastjet/JetDefinition.hh" #include "fastjet/PseudoJet.hh" -#include "Geometry/CSCGeometry/interface/CSCGeometry.h" -#include "Geometry/DTGeometry/interface/DTGeometry.h" -#include "DataFormats/MuonDetId/interface/CSCDetId.h" -#include "DataFormats/Math/interface/Vector3D.h" -#include "DataFormats/Math/interface/PtEtaPhiMass.h" -#include "DataFormats/DTRecHit/interface/DTRecHitCollection.h" -#include "DataFormats/CSCRecHit/interface/CSCRecHit2DCollection.h" -#include "DataFormats/MuonReco/interface/MuonRecHitCluster.h" - typedef std::vector RecHitClusterCollection; template @@ -47,8 +54,8 @@ class RechitClusterProducerT : public edm::global::EDProducer<> { private: void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override; - const edm::ESGetToken GeometryToken_; - edm::EDGetTokenT InputToken_; + const edm::ESGetToken geometryToken_; + edm::EDGetTokenT inputToken_; typedef std::vector RecHitClusterCollection; const double rParam_; // distance paramter @@ -58,8 +65,8 @@ class RechitClusterProducerT : public edm::global::EDProducer<> { template RechitClusterProducerT::RechitClusterProducerT(const edm::ParameterSet& iConfig) - : GeometryToken_(esConsumes()), - InputToken_(consumes(iConfig.getParameter("recHitLabel"))), + : geometryToken_(esConsumes()), + inputToken_(consumes(iConfig.getParameter("recHitLabel"))), rParam_(iConfig.getParameter("rParam")), nRechitMin_(iConfig.getParameter("nRechitMin")), nStationThres_(iConfig.getParameter("nStationThres")) { @@ -73,31 +80,28 @@ RechitClusterProducerT::RechitClusterProducerT(const edm::ParameterSet& i // ------------ method called to produce the data ------------ template void RechitClusterProducerT::produce(edm::StreamID, edm::Event& ev, const edm::EventSetup& iSetup) const { - auto const& geo = iSetup.getData(GeometryToken_); - auto const& rechits = ev.get(InputToken_); + auto const& geo = iSetup.getData(geometryToken_); + auto const& rechits = ev.get(inputToken_); std::vector fjInput; - RecHitRefVector inputs_; + RecHitRefVector inputs; fastjet::JetDefinition jet_def(fastjet::cambridge_algorithm, rParam_); - inputs_.clear(); - fjInput.clear(); - - Trait RecHitTrait; + Trait recHitTrait; int recIt = 0; for (auto const& rechit : rechits) { - LocalPoint RecHitLocalPosition = rechit.localPosition(); - auto detid = RecHitTrait.detid(rechit); + LocalPoint recHitLocalPosition = rechit.localPosition(); + auto detid = recHitTrait.detid(rechit); auto thischamber = geo.chamber(detid); if (thischamber) { - GlobalPoint globalPosition = thischamber->toGlobal(RecHitLocalPosition); + GlobalPoint globalPosition = thischamber->toGlobal(recHitLocalPosition); float x = globalPosition.x(); float y = globalPosition.y(); float z = globalPosition.z(); RechitRef ref = RechitRef(&rechits, recIt); - inputs_.push_back(ref); + inputs.push_back(ref); fjInput.push_back(fastjet::PseudoJet(x, y, z, globalPosition.mag())); fjInput.back().set_user_index(recIt); } @@ -107,56 +111,51 @@ void RechitClusterProducerT::produce(edm::StreamID, edm::Event& ev, const //keep all the clusters double ptmin = 0.0; - std::vector fjJets = fastjet::sorted_by_pt(clus_seq.inclusive_jets(ptmin)); + std::vector fjJets = clus_seq.inclusive_jets(ptmin); auto clusters = std::make_unique(); - if (!fjJets.empty()) { - for (unsigned int ijet = 0; ijet < fjJets.size(); ++ijet) { - // get the fastjet jet - const fastjet::PseudoJet& fjJet = fjJets[ijet]; - - // skip if the cluster has too few rechits - if (int(fjJet.constituents().size()) < nRechitMin_) - continue; - // get the constituents from fastjet - RecHitRefVector rechits; - for (unsigned int i = 0; i < fjJet.constituents().size(); i++) { - auto index = fjJet.constituents()[i].user_index(); - if (index >= 0 && static_cast(index) < inputs_.size()) { - rechits.push_back(inputs_[index]); - } + for (auto const& fjJet : fjJets) { + // skip if the cluster has too few rechits + if (int(fjJet.constituents().size()) < nRechitMin_) + continue; + // get the constituents from fastjet + RecHitRefVector rechits; + for (auto const& constituent : fjJet.constituents()) { + auto index = constituent.user_index(); + if (index >= 0 && static_cast(index) < inputs.size()) { + rechits.push_back(inputs[index]); } + } - //Derive cluster properties - int nStation = 0; - int totStation = 0; - float avgStation = 0.0; - std::map station_count_map; - for (auto& rechit : rechits) { - station_count_map[RecHitTrait.station(*rechit)]++; - } - //station statistics - std::map::iterator it; - for (auto const& [station, count] : station_count_map) { - if (count >= nStationThres_) { - nStation++; - avgStation += (station) * (count); - totStation += (count); - } - } - if (totStation != 0) { - avgStation = avgStation / totStation; + //Derive cluster properties + int nStation = 0; + int totStation = 0; + float avgStation = 0.0; + std::map station_count_map; + for (auto const& rechit : rechits) { + station_count_map[recHitTrait.station(*rechit)]++; + } + //station statistics + std::map::iterator it; + for (auto const& [station, count] : station_count_map) { + if (count >= nStationThres_) { + nStation++; + avgStation += station * count; + totStation += count; } - float invN = 1.f / rechits.size(); - // cluster position is the average position of the constituent rechits - float jetX = fjJet.px() * invN; - float jetY = fjJet.py() * invN; - float jetZ = fjJet.pz() * invN; - - math::RhoEtaPhiVectorF position( - std::sqrt(jetX * jetX + jetY * jetY), etaFromXYZ(jetX, jetY, jetZ), std::atan2(jetY, jetX)); - Trait::emplace_back(clusters.get(), position, nStation, avgStation, rechits); } + if (totStation != 0) { + avgStation = avgStation / totStation; + } + float invN = 1.f / rechits.size(); + // cluster position is the average position of the constituent rechits + float jetX = fjJet.px() * invN; + float jetY = fjJet.py() * invN; + float jetZ = fjJet.pz() * invN; + + math::RhoEtaPhiVectorF position( + std::sqrt(jetX * jetX + jetY * jetY), etaFromXYZ(jetX, jetY, jetZ), std::atan2(jetY, jetX)); + Trait::emplace_back(clusters.get(), position, nStation, avgStation, rechits); } ev.put(std::move(clusters)); } @@ -177,7 +176,7 @@ struct DTRecHitTrait { using RecHitRef = edm::Ref; using RecHitRefVector = edm::RefVector; static std::string recHitLabel() { return "dt1DRecHits"; } - static std::string producerName() { return "DTrechitClusterProducer"; } + static std::string producerName() { return "dtRechitClusterProducer"; } static int station(const DTRecHit1DPair& dtRechit) { return detid(dtRechit).station(); } static DTChamberId detid(const DTRecHit1DPair& dtRechit) { @@ -193,7 +192,7 @@ struct DTRecHitTrait { // compute nMB1, nMB2 int nMB1 = 0; int nMB2 = 0; - for (auto& rechit : rechits) { + for (auto const& rechit : rechits) { DetId geoid = rechit->geographicalId(); DTChamberId dtdetid = DTChamberId(geoid); if (dtdetid.station() == 1) @@ -213,7 +212,7 @@ struct CSCRecHitTrait { using RecHitRef = edm::Ref; using RecHitRefVector = edm::RefVector; static std::string recHitLabel() { return "csc2DRecHits"; } - static std::string producerName() { return "CSCrechitClusterProducer"; } + static std::string producerName() { return "cscRechitClusterProducer"; } static int station(const CSCRecHit2D& cscRechit) { return CSCDetId::station(detid(cscRechit)); } static CSCDetId detid(const CSCRecHit2D& cscRechit) { return cscRechit.cscDetId(); } @@ -229,7 +228,7 @@ struct CSCRecHitTrait { float timeSpread = 0.0; float time = 0.0; float time_strip = 0.0; // for timeSpread calculation - for (auto& rechit : rechits) { + for (auto const& rechit : rechits) { CSCDetId cscdetid = rechit->cscDetId(); int stationRing = (CSCDetId::station(cscdetid) * 10 + CSCDetId::ring(cscdetid)); if (CSCDetId::ring(cscdetid) == 4) diff --git a/RecoMuon/MuonRechitClusterProducer/python/rechitCluster_cff.py b/RecoMuon/MuonRechitClusterProducer/python/rechitCluster_cff.py index f4dc69e7f77c8..bdc5ebd08aca1 100644 --- a/RecoMuon/MuonRechitClusterProducer/python/rechitCluster_cff.py +++ b/RecoMuon/MuonRechitClusterProducer/python/rechitCluster_cff.py @@ -1,15 +1,15 @@ import FWCore.ParameterSet.Config as cms -import RecoMuon.MuonRechitClusterProducer.CSCrechitClusterProducer_cfi as CSCcluster -import RecoMuon.MuonRechitClusterProducer.DTrechitClusterProducer_cfi as DTcluster +import RecoMuon.MuonRechitClusterProducer.cscRechitClusterProducer_cfi as CSCcluster +import RecoMuon.MuonRechitClusterProducer.dtRechitClusterProducer_cfi as DTcluster -ca4CSCrechitClusters= CSCcluster.CSCrechitClusterProducer.clone( +ca4CSCrechitClusters= CSCcluster.cscRechitClusterProducer.clone( recHitLabel = "csc2DRecHits", nRechitMin = 50, rParam = 0.4, nStationThres = 10, ) -ca4DTrechitClusters = DTcluster.DTrechitClusterProducer.clone( +ca4DTrechitClusters = DTcluster.dtRechitClusterProducer.clone( recHitLabel = "dt1DRecHits", nRechitMin = 50, rParam = 0.4,