File indexing completed on 2024-05-31 04:19:51
0001
0002
0003
0004
0005
0006
0007
0008
0009 #include <memory>
0010 #include <cmath>
0011 #include <tinyxml2.h>
0012
0013
0014 #include "FWCore/Framework/interface/Frameworkfwd.h"
0015 #include "FWCore/Framework/interface/stream/EDProducer.h"
0016
0017 #include "FWCore/Framework/interface/Event.h"
0018 #include "FWCore/Framework/interface/MakerMacros.h"
0019
0020 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0021 #include "FWCore/Utilities/interface/StreamID.h"
0022
0023
0024 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0025 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0026
0027 #include "CommonTools/MVAUtils/interface/TMVAZipReader.h"
0028
0029
0030 #include "DataFormats/TrajectorySeed/interface/TrajectorySeed.h"
0031 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
0032 #include "DataFormats/TrajectorySeed/interface/PropagationDirection.h"
0033 #include "DataFormats/TrajectoryState/interface/PTrajectoryStateOnDet.h"
0034 #include "DataFormats/TrajectoryState/interface/LocalTrajectoryParameters.h"
0035 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0036
0037 #include "RecoMuon/TrackerSeedGenerator/interface/SeedMvaEstimator.h"
0038
0039
0040 bool sortByMvaScore(const std::pair<unsigned, double>& A, const std::pair<unsigned, double>& B) {
0041 return (A.second > B.second);
0042 };
0043
0044 class MuonHLTSeedMVAClassifier : public edm::stream::EDProducer<> {
0045 public:
0046 explicit MuonHLTSeedMVAClassifier(const edm::ParameterSet&);
0047 ~MuonHLTSeedMVAClassifier() override = default;
0048
0049 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0050 bool checkMVAFileConsistency(const std::string& weightsFileFullPath, bool isFromL1) const;
0051
0052 private:
0053 void produce(edm::Event&, const edm::EventSetup&) override;
0054
0055
0056 const edm::EDGetTokenT<TrajectorySeedCollection> seedToken_;
0057 const edm::EDGetTokenT<l1t::MuonBxCollection> l1MuonToken_;
0058 const edm::EDGetTokenT<reco::RecoChargedCandidateCollection> l2MuonToken_;
0059 const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> trackerGeometryToken_;
0060
0061 typedef std::pair<std::unique_ptr<const SeedMvaEstimator>, std::unique_ptr<const SeedMvaEstimator>>
0062 PairSeedMvaEstimator;
0063 PairSeedMvaEstimator mvaEstimator_;
0064
0065 const bool rejectAll_;
0066 const bool isFromL1_;
0067
0068 const edm::FileInPath mvaFileB_;
0069 const edm::FileInPath mvaFileE_;
0070
0071 const std::vector<double> mvaScaleMeanB_;
0072 const std::vector<double> mvaScaleStdB_;
0073 const std::vector<double> mvaScaleMeanE_;
0074 const std::vector<double> mvaScaleStdE_;
0075
0076 const bool doSort_;
0077 const int nSeedsMaxB_;
0078 const int nSeedsMaxE_;
0079
0080 const double etaEdge_;
0081 const double mvaCutB_;
0082 const double mvaCutE_;
0083
0084 const int minL1Qual_;
0085 const double baseScore_;
0086
0087 double getSeedMva(const PairSeedMvaEstimator& pairMvaEstimator,
0088 const TrajectorySeed& seed,
0089 const GlobalVector& global_p,
0090 const l1t::MuonBxCollection& l1Muons,
0091 const reco::RecoChargedCandidateCollection& l2Muons);
0092 };
0093
0094 bool MuonHLTSeedMVAClassifier::checkMVAFileConsistency(const std::string& weightsFileFullPath,
0095 const bool isFromL1) const {
0096 tinyxml2::XMLDocument xmlDoc;
0097 if (reco::details::hasEnding(weightsFileFullPath, ".xml")) {
0098 xmlDoc.LoadFile(weightsFileFullPath.c_str());
0099 } else {
0100 edm::LogError("MuonHLTSeedMVAClassifier") << "unsupported file extension, it should be a .xml file!";
0101 return false;
0102 }
0103 tinyxml2::XMLElement* root = xmlDoc.FirstChildElement("MethodSetup");
0104 if (root == nullptr) {
0105 edm::LogError("MuonHLTSeedMVAClassifier") << "could not retrieve the MethodSetup node from XML file!";
0106 return false;
0107 }
0108
0109 const auto& vars = root->FirstChildElement("Variables");
0110 size_t n = 0;
0111 if (vars != nullptr) {
0112 for (tinyxml2::XMLElement* e = vars->FirstChildElement("Variable"); e != nullptr;
0113 e = e->NextSiblingElement("Variable")) {
0114 ++n;
0115 }
0116 } else {
0117 edm::LogError("MuonHLTSeedMVAClassifier") << "could not retrieve the Variables node from XML file!";
0118 return false;
0119 }
0120
0121 LogTrace("MuonHLTSeedMVAClassifier") << "MVA file:" << weightsFileFullPath.c_str() << " n Var:" << n;
0122 bool condition = (isFromL1 && (n == inputIndexes::kLastL1)) || (!isFromL1 && (n == inputIndexes::kLastL2));
0123 return condition;
0124 }
0125
0126 MuonHLTSeedMVAClassifier::MuonHLTSeedMVAClassifier(const edm::ParameterSet& iConfig)
0127 : seedToken_(consumes<TrajectorySeedCollection>(iConfig.getParameter<edm::InputTag>("src"))),
0128 l1MuonToken_(consumes<l1t::MuonBxCollection>(iConfig.getParameter<edm::InputTag>("L1Muon"))),
0129 l2MuonToken_(consumes<reco::RecoChargedCandidateCollection>(iConfig.getParameter<edm::InputTag>("L2Muon"))),
0130 trackerGeometryToken_(esConsumes<TrackerGeometry, TrackerDigiGeometryRecord>()),
0131 rejectAll_(iConfig.getParameter<bool>("rejectAll")),
0132 isFromL1_(iConfig.getParameter<bool>("isFromL1")),
0133 mvaFileB_(iConfig.getParameter<edm::FileInPath>("mvaFileB")),
0134 mvaFileE_(iConfig.getParameter<edm::FileInPath>("mvaFileE")),
0135 mvaScaleMeanB_(iConfig.getParameter<std::vector<double>>("mvaScaleMeanB")),
0136 mvaScaleStdB_(iConfig.getParameter<std::vector<double>>("mvaScaleStdB")),
0137 mvaScaleMeanE_(iConfig.getParameter<std::vector<double>>("mvaScaleMeanE")),
0138 mvaScaleStdE_(iConfig.getParameter<std::vector<double>>("mvaScaleStdE")),
0139 doSort_(iConfig.getParameter<bool>("doSort")),
0140 nSeedsMaxB_(iConfig.getParameter<int>("nSeedsMaxB")),
0141 nSeedsMaxE_(iConfig.getParameter<int>("nSeedsMaxE")),
0142 etaEdge_(iConfig.getParameter<double>("etaEdge")),
0143 mvaCutB_(iConfig.getParameter<double>("mvaCutB")),
0144 mvaCutE_(iConfig.getParameter<double>("mvaCutE")),
0145 minL1Qual_(iConfig.getParameter<int>("minL1Qual")),
0146 baseScore_(iConfig.getParameter<double>("baseScore")) {
0147 const auto& mvaFileBPath = mvaFileB_.fullPath();
0148 const auto& mvaFileEPath = mvaFileE_.fullPath();
0149
0150 if (!checkMVAFileConsistency(mvaFileBPath, isFromL1_) || !checkMVAFileConsistency(mvaFileEPath, isFromL1_)) {
0151 throw cms::Exception("ConfigurationError") << " MVA files appear to be not consistent with the value of isFromL1 "
0152 "parameter.\n Please check your configuration.";
0153 }
0154
0155 if (!rejectAll_) {
0156 mvaEstimator_ = std::make_pair(
0157 std::make_unique<SeedMvaEstimator>(mvaFileB_, mvaScaleMeanB_, mvaScaleStdB_, isFromL1_, minL1Qual_),
0158 std::make_unique<SeedMvaEstimator>(mvaFileE_, mvaScaleMeanE_, mvaScaleStdE_, isFromL1_, minL1Qual_));
0159 }
0160
0161 produces<TrajectorySeedCollection>();
0162 }
0163
0164
0165 void MuonHLTSeedMVAClassifier::produce(edm::Event& iEvent, edm::EventSetup const& iEventSetup) {
0166 auto result = std::make_unique<TrajectorySeedCollection>();
0167
0168 if (rejectAll_) {
0169 iEvent.put(std::move(result));
0170 return;
0171 }
0172
0173 if (doSort_ && nSeedsMaxB_ <= 0 && nSeedsMaxE_ <= 0) {
0174 iEvent.put(std::move(result));
0175 return;
0176 }
0177
0178 if (!doSort_ && mvaCutB_ > 1. && mvaCutE_ > 1.) {
0179 iEvent.put(std::move(result));
0180 return;
0181 }
0182
0183 const TrajectorySeedCollection& seeds = iEvent.get(seedToken_);
0184 const l1t::MuonBxCollection& l1Muons = iEvent.get(l1MuonToken_);
0185 const reco::RecoChargedCandidateCollection& l2Muons = iEvent.get(l2MuonToken_);
0186 const TrackerGeometry& trkGeom = iEventSetup.getData(trackerGeometryToken_);
0187
0188 std::vector<std::pair<unsigned, double>> pairSeedIdxMvaScoreB = {};
0189 std::vector<std::pair<unsigned, double>> pairSeedIdxMvaScoreE = {};
0190 for (auto& seed : seeds) {
0191 const GlobalVector global_p =
0192 trkGeom.idToDet(seed.startingState().detId())->surface().toGlobal(seed.startingState().parameters().momentum());
0193
0194 bool isB = (std::abs(global_p.eta()) < etaEdge_);
0195
0196 if (doSort_) {
0197 if (isB) {
0198 if (nSeedsMaxB_ <= 0) {
0199 continue;
0200 }
0201 } else {
0202 if (nSeedsMaxE_ <= 0) {
0203 continue;
0204 }
0205 }
0206 } else {
0207 if (isB) {
0208 if (mvaCutB_ > 1.0) {
0209 continue;
0210 } else if (mvaCutB_ <= 0.) {
0211 result->emplace_back(seed);
0212 continue;
0213 }
0214 } else {
0215 if (mvaCutE_ > 1.0) {
0216 continue;
0217 } else if (mvaCutE_ <= 0.) {
0218 result->emplace_back(seed);
0219 continue;
0220 }
0221 }
0222 }
0223
0224 double mva = getSeedMva(mvaEstimator_, seed, global_p, l1Muons, l2Muons);
0225
0226 double score = 1. / (1. + std::exp(-1. * mva));
0227 bool passMva = isB ? score > mvaCutB_ : score > mvaCutE_;
0228 if (!passMva)
0229 continue;
0230
0231 if (doSort_) {
0232 if (isB)
0233 pairSeedIdxMvaScoreB.push_back(std::make_pair(&seed - &seeds.at(0), score));
0234 else
0235 pairSeedIdxMvaScoreE.push_back(std::make_pair(&seed - &seeds.at(0), score));
0236 } else {
0237 result->emplace_back(seed);
0238 }
0239 }
0240
0241 if (doSort_) {
0242 std::sort(pairSeedIdxMvaScoreB.begin(), pairSeedIdxMvaScoreB.end(), sortByMvaScore);
0243 std::sort(pairSeedIdxMvaScoreE.begin(), pairSeedIdxMvaScoreE.end(), sortByMvaScore);
0244
0245 for (auto i = 0U; i < pairSeedIdxMvaScoreB.size(); ++i) {
0246 if ((int)i == nSeedsMaxB_)
0247 break;
0248 const auto& seed(seeds.at(pairSeedIdxMvaScoreB.at(i).first));
0249 result->emplace_back(seed);
0250 }
0251
0252 for (auto i = 0U; i < pairSeedIdxMvaScoreE.size(); ++i) {
0253 if ((int)i == nSeedsMaxE_)
0254 break;
0255 const auto& seed(seeds.at(pairSeedIdxMvaScoreE.at(i).first));
0256 result->emplace_back(seed);
0257 }
0258 }
0259
0260 iEvent.put(std::move(result));
0261 }
0262
0263 double MuonHLTSeedMVAClassifier::getSeedMva(const PairSeedMvaEstimator& pairMvaEstimator,
0264 const TrajectorySeed& seed,
0265 const GlobalVector& global_p,
0266 const l1t::MuonBxCollection& l1Muons,
0267 const reco::RecoChargedCandidateCollection& l2Muons) {
0268 double mva = 0.;
0269 if (std::abs(global_p.eta()) < etaEdge_) {
0270 mva = pairMvaEstimator.first->computeMva(seed, global_p, l1Muons, l2Muons);
0271 } else {
0272 mva = pairMvaEstimator.second->computeMva(seed, global_p, l1Muons, l2Muons);
0273 }
0274
0275 return (mva + baseScore_);
0276 }
0277
0278
0279 void MuonHLTSeedMVAClassifier::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0280 edm::ParameterSetDescription desc;
0281 desc.add<edm::InputTag>("src", edm::InputTag("hltIter2IterL3MuonPixelSeeds", ""));
0282 desc.add<edm::InputTag>("L1Muon", edm::InputTag("hltGtStage2Digis", "Muon"));
0283 desc.add<edm::InputTag>("L2Muon", edm::InputTag("hltL2MuonCandidates", ""));
0284
0285 desc.add<bool>("rejectAll", false);
0286 desc.add<bool>("isFromL1", false);
0287
0288 desc.add<edm::FileInPath>("mvaFileB",
0289 edm::FileInPath("RecoMuon/TrackerSeedGenerator/data/xgb_Run3_Iter2FromL1Seeds_barrel.xml"));
0290 desc.add<edm::FileInPath>("mvaFileE",
0291 edm::FileInPath("RecoMuon/TrackerSeedGenerator/data/xgb_Run3_Iter2FromL1Seeds_endcap.xml"));
0292 desc.add<std::vector<double>>("mvaScaleMeanB", {0., 0., 0., 0., 0., 0., 0., 0.});
0293 desc.add<std::vector<double>>("mvaScaleStdB", {1., 1., 1., 1., 1., 1., 1., 1.});
0294 desc.add<std::vector<double>>("mvaScaleMeanE", {0., 0., 0., 0., 0., 0., 0., 0.});
0295 desc.add<std::vector<double>>("mvaScaleStdE", {1., 1., 1., 1., 1., 1., 1., 1.});
0296
0297 desc.add<bool>("doSort", false);
0298 desc.add<int>("nSeedsMaxB", 1e6);
0299 desc.add<int>("nSeedsMaxE", 1e6);
0300
0301 desc.add<double>("etaEdge", 1.2);
0302 desc.add<double>("mvaCutB", -1.);
0303 desc.add<double>("mvaCutE", -1.);
0304
0305 desc.add<int>("minL1Qual", 7);
0306 desc.add<double>("baseScore", 0.5);
0307
0308 descriptions.add("MuonHLTSeedMVAClassifier", desc);
0309 }
0310
0311
0312 DEFINE_FWK_MODULE(MuonHLTSeedMVAClassifier);
