File indexing completed on 2024-12-20 03:14:08
0001
0002
0003
0004
0005
0006
0007
0008 #include <memory>
0009 #include <cmath>
0010
0011
0012 #include "FWCore/Framework/interface/Frameworkfwd.h"
0013 #include "FWCore/Framework/interface/stream/EDProducer.h"
0014 #include "FWCore/Framework/interface/Event.h"
0015 #include "FWCore/Framework/interface/MakerMacros.h"
0016 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0018 #include "FWCore/Utilities/interface/StreamID.h"
0019 #include "FWCore/Utilities/interface/ESGetToken.h"
0020 #include "FWCore/Utilities/interface/ESInputTag.h"
0021
0022
0023 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0024 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0025
0026
0027 #include "DataFormats/TrajectorySeed/interface/TrajectorySeed.h"
0028 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
0029 #include "DataFormats/TrajectorySeed/interface/PropagationDirection.h"
0030 #include "DataFormats/TrajectoryState/interface/PTrajectoryStateOnDet.h"
0031 #include "DataFormats/TrajectoryState/interface/LocalTrajectoryParameters.h"
0032 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0033
0034
0035 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0036 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0037 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0038 #include "RecoTracker/TkDetLayers/interface/GeometricSearchTracker.h"
0039 #include "RecoTracker/TkDetLayers/interface/GeometricSearchTrackerBuilder.h"
0040 #include "Geometry/TrackerNumberingBuilder/interface/GeometricDet.h"
0041 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0042 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0043 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0044
0045 #include "RecoMuon/TrackerSeedGenerator/interface/SeedMvaEstimatorPhase2.h"
0046
0047
0048 bool sortByMvaScorePhase2(const std::pair<unsigned, double>& A, const std::pair<unsigned, double>& B) {
0049 return (A.second > B.second);
0050 };
0051
0052 class MuonHLTSeedMVAClassifierPhase2 : public edm::stream::EDProducer<> {
0053 public:
0054 explicit MuonHLTSeedMVAClassifierPhase2(const edm::ParameterSet&);
0055 ~MuonHLTSeedMVAClassifierPhase2() override = default;
0056
0057 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0058
0059 private:
0060 void produce(edm::Event&, const edm::EventSetup&) override;
0061
0062
0063 const edm::EDGetTokenT<TrajectorySeedCollection> t_Seed_;
0064 const edm::EDGetTokenT<l1t::TrackerMuonCollection> t_L1TkMu_;
0065
0066 typedef std::pair<std::unique_ptr<const SeedMvaEstimatorPhase2>, std::unique_ptr<const SeedMvaEstimatorPhase2>>
0067 pairSeedMvaEstimator;
0068 pairSeedMvaEstimator mvaEstimator_;
0069
0070 const edm::FileInPath mvaFile_B_0_;
0071 const edm::FileInPath mvaFile_E_0_;
0072
0073 const std::vector<double> mvaScaleMean_B_;
0074 const std::vector<double> mvaScaleStd_B_;
0075 const std::vector<double> mvaScaleMean_E_;
0076 const std::vector<double> mvaScaleStd_E_;
0077
0078 const double etaEdge_;
0079 const double mvaCut_B_;
0080 const double mvaCut_E_;
0081
0082 const bool doSort_;
0083 const int nSeedsMax_B_;
0084 const int nSeedsMax_E_;
0085
0086 const double baseScore_;
0087
0088 const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> trackerTopologyESToken_;
0089 const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> trackerGeometryESToken_;
0090 const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magFieldESToken_;
0091 const edm::ESGetToken<GeometricDet, IdealGeometryRecord> geomDetESToken_;
0092 const edm::ESGetToken<Propagator, TrackingComponentsRecord> propagatorESToken_;
0093
0094 double getSeedMva(const pairSeedMvaEstimator& pairMvaEstimator,
0095 const TrajectorySeed& seed,
0096 const GlobalVector& global_p,
0097 const GlobalPoint& global_x,
0098 const edm::Handle<l1t::TrackerMuonCollection>& h_L1TkMu,
0099 const edm::ESHandle<MagneticField>& magfieldH,
0100 const Propagator& propagatorAlong,
0101 const GeometricSearchTracker& geomTracker);
0102 };
0103
0104
0105 MuonHLTSeedMVAClassifierPhase2::MuonHLTSeedMVAClassifierPhase2(const edm::ParameterSet& iConfig)
0106 : t_Seed_(consumes<TrajectorySeedCollection>(iConfig.getParameter<edm::InputTag>("src"))),
0107 t_L1TkMu_(consumes<l1t::TrackerMuonCollection>(iConfig.getParameter<edm::InputTag>("L1TkMu"))),
0108
0109 mvaFile_B_0_(iConfig.getParameter<edm::FileInPath>("mvaFile_B_0")),
0110 mvaFile_E_0_(iConfig.getParameter<edm::FileInPath>("mvaFile_E_0")),
0111
0112 mvaScaleMean_B_(iConfig.getParameter<std::vector<double>>("mvaScaleMean_B")),
0113 mvaScaleStd_B_(iConfig.getParameter<std::vector<double>>("mvaScaleStd_B")),
0114 mvaScaleMean_E_(iConfig.getParameter<std::vector<double>>("mvaScaleMean_E")),
0115 mvaScaleStd_E_(iConfig.getParameter<std::vector<double>>("mvaScaleStd_E")),
0116
0117 etaEdge_(iConfig.getParameter<double>("etaEdge")),
0118 mvaCut_B_(iConfig.getParameter<double>("mvaCut_B")),
0119 mvaCut_E_(iConfig.getParameter<double>("mvaCut_E")),
0120
0121 doSort_(iConfig.getParameter<bool>("doSort")),
0122 nSeedsMax_B_(iConfig.getParameter<int>("nSeedsMax_B")),
0123 nSeedsMax_E_(iConfig.getParameter<int>("nSeedsMax_E")),
0124
0125 baseScore_(iConfig.getParameter<double>("baseScore")),
0126
0127 trackerTopologyESToken_(esConsumes<TrackerTopology, TrackerTopologyRcd>()),
0128 trackerGeometryESToken_(esConsumes<TrackerGeometry, TrackerDigiGeometryRecord>()),
0129 magFieldESToken_(esConsumes<MagneticField, IdealMagneticFieldRecord>()),
0130 geomDetESToken_(esConsumes<GeometricDet, IdealGeometryRecord>()),
0131 propagatorESToken_(
0132 esConsumes<Propagator, TrackingComponentsRecord>(edm::ESInputTag("", "PropagatorWithMaterialParabolicMf"))) {
0133 produces<TrajectorySeedCollection>();
0134
0135 mvaEstimator_ =
0136 std::make_pair(std::make_unique<SeedMvaEstimatorPhase2>(mvaFile_B_0_, mvaScaleMean_B_, mvaScaleStd_B_),
0137 std::make_unique<SeedMvaEstimatorPhase2>(mvaFile_E_0_, mvaScaleMean_E_, mvaScaleStd_E_));
0138 }
0139
0140
0141
0142 void MuonHLTSeedMVAClassifierPhase2::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0143 auto result = std::make_unique<TrajectorySeedCollection>();
0144
0145 edm::Handle<l1t::TrackerMuonCollection> h_L1TkMu;
0146 bool hasL1TkMu = iEvent.getByToken(t_L1TkMu_, h_L1TkMu);
0147
0148 edm::Handle<TrajectorySeedCollection> h_Seed;
0149 bool hasSeed = iEvent.getByToken(t_Seed_, h_Seed);
0150
0151 if (!(hasL1TkMu && hasSeed)) {
0152 edm::LogError("SeedClassifierError") << "Error! Cannot find L1TkMuon or TrajectorySeed\n"
0153 << "hasL1TkMu : " << hasL1TkMu << "\n"
0154 << "hasSeed : " << hasSeed << "\n";
0155 return;
0156 }
0157
0158 if (h_L1TkMu->empty() or h_Seed->empty()) {
0159 if (!h_Seed->empty()) {
0160 edm::LogInfo("SeedClassifierError") << "Empty L1TkMu collection" << '\n';
0161 }
0162 if (!h_L1TkMu->empty()) {
0163 edm::LogInfo("SeedClassifierError") << "Empty Muon Pixel seeds collection" << '\n';
0164 } else {
0165 edm::LogInfo("SeedClassifierError") << "Empty L1TkMu and Muon Pixel seeds collections" << '\n';
0166 }
0167 iEvent.put(std::move(result));
0168 return;
0169 }
0170
0171 edm::ESHandle<TrackerGeometry> trkGeom = iSetup.getHandle(trackerGeometryESToken_);
0172 edm::ESHandle<GeometricDet> geomDet = iSetup.getHandle(geomDetESToken_);
0173 edm::ESHandle<TrackerTopology> trkTopo = iSetup.getHandle(trackerTopologyESToken_);
0174
0175 GeometricSearchTrackerBuilder builder;
0176 GeometricSearchTracker geomTracker = *(builder.build(&(*geomDet), &(*trkGeom), &(*trkTopo)));
0177
0178 edm::ESHandle<MagneticField> magfieldH = iSetup.getHandle(magFieldESToken_);
0179 edm::ESHandle<Propagator> propagatorAlongH = iSetup.getHandle(propagatorESToken_);
0180 std::unique_ptr<Propagator> propagatorAlong = SetPropagationDirection(*propagatorAlongH, alongMomentum);
0181
0182
0183 if (doSort_) {
0184 std::vector<std::pair<unsigned, double>> pairSeedIdxMvaScore_B = {};
0185 std::vector<std::pair<unsigned, double>> pairSeedIdxMvaScore_E = {};
0186
0187 for (auto i = 0U; i < h_Seed->size(); ++i) {
0188 const auto& seed(h_Seed->at(i));
0189
0190 GlobalVector global_p = trkGeom->idToDet(seed.startingState().detId())
0191 ->surface()
0192 .toGlobal(seed.startingState().parameters().momentum());
0193 GlobalPoint global_x = trkGeom->idToDet(seed.startingState().detId())
0194 ->surface()
0195 .toGlobal(seed.startingState().parameters().position());
0196
0197 bool isB = (std::abs(global_p.eta()) < etaEdge_);
0198
0199 if (isB && nSeedsMax_B_ < 0) {
0200 result->emplace_back(seed);
0201 continue;
0202 }
0203
0204 if (!isB && nSeedsMax_E_ < 0) {
0205 result->emplace_back(seed);
0206 continue;
0207 }
0208
0209 double mva = getSeedMva(
0210 mvaEstimator_, seed, global_p, global_x, h_L1TkMu, magfieldH, *(propagatorAlong.get()), geomTracker);
0211
0212 double logistic = 1 / (1 + std::exp(-mva));
0213
0214 if (isB)
0215 pairSeedIdxMvaScore_B.push_back(make_pair(i, logistic));
0216 else
0217 pairSeedIdxMvaScore_E.push_back(make_pair(i, logistic));
0218 }
0219
0220 std::partial_sort(pairSeedIdxMvaScore_B.begin(),
0221 std::min(pairSeedIdxMvaScore_B.begin() + nSeedsMax_B_, pairSeedIdxMvaScore_B.end()),
0222 pairSeedIdxMvaScore_B.end(),
0223 sortByMvaScorePhase2);
0224 std::partial_sort(pairSeedIdxMvaScore_E.begin(),
0225 std::min(pairSeedIdxMvaScore_E.begin() + nSeedsMax_E_, pairSeedIdxMvaScore_E.end()),
0226 pairSeedIdxMvaScore_E.end(),
0227 sortByMvaScorePhase2);
0228
0229 for (size_t i = 0; i < std::min(pairSeedIdxMvaScore_B.size(), static_cast<size_t>(nSeedsMax_B_)); ++i) {
0230 const auto& seed(h_Seed->at(pairSeedIdxMvaScore_B.at(i).first));
0231 result->emplace_back(seed);
0232 }
0233
0234 for (size_t i = 0; i < std::min(pairSeedIdxMvaScore_E.size(), static_cast<size_t>(nSeedsMax_E_)); ++i) {
0235 const auto& seed(h_Seed->at(pairSeedIdxMvaScore_E.at(i).first));
0236 result->emplace_back(seed);
0237 }
0238 }
0239
0240
0241 else {
0242 for (auto i = 0U; i < h_Seed->size(); ++i) {
0243 const auto& seed(h_Seed->at(i));
0244
0245 GlobalVector global_p = trkGeom->idToDet(seed.startingState().detId())
0246 ->surface()
0247 .toGlobal(seed.startingState().parameters().momentum());
0248 GlobalPoint global_x = trkGeom->idToDet(seed.startingState().detId())
0249 ->surface()
0250 .toGlobal(seed.startingState().parameters().position());
0251
0252 bool isB = (std::abs(global_p.eta()) < etaEdge_);
0253
0254 if (isB && mvaCut_B_ <= 0.) {
0255 result->emplace_back(seed);
0256 continue;
0257 }
0258
0259 if (!isB && mvaCut_E_ <= 0.) {
0260 result->emplace_back(seed);
0261 continue;
0262 }
0263
0264 double mva = getSeedMva(
0265 mvaEstimator_, seed, global_p, global_x, h_L1TkMu, magfieldH, *(propagatorAlong.get()), geomTracker);
0266
0267 double logistic = 1 / (1 + std::exp(-mva));
0268
0269 bool passMva = ((isB && (logistic > mvaCut_B_)) || (!isB && (logistic > mvaCut_E_)));
0270
0271 if (passMva)
0272 result->emplace_back(seed);
0273 }
0274 }
0275
0276 iEvent.put(std::move(result));
0277 }
0278
0279 double MuonHLTSeedMVAClassifierPhase2::getSeedMva(const pairSeedMvaEstimator& pairMvaEstimator,
0280 const TrajectorySeed& seed,
0281 const GlobalVector& global_p,
0282 const GlobalPoint& global_x,
0283 const edm::Handle<l1t::TrackerMuonCollection>& h_L1TkMu,
0284 const edm::ESHandle<MagneticField>& magfieldH,
0285 const Propagator& propagatorAlong,
0286 const GeometricSearchTracker& geomTracker) {
0287 double mva = 0.;
0288
0289 if (std::abs(global_p.eta()) < etaEdge_) {
0290 mva =
0291 pairMvaEstimator.first->computeMva(seed, global_p, global_x, h_L1TkMu, magfieldH, propagatorAlong, geomTracker);
0292 } else {
0293 mva = pairMvaEstimator.second->computeMva(
0294 seed, global_p, global_x, h_L1TkMu, magfieldH, propagatorAlong, geomTracker);
0295 }
0296
0297 return (baseScore_ + mva);
0298 }
0299
0300
0301 void MuonHLTSeedMVAClassifierPhase2::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0302 edm::ParameterSetDescription desc;
0303 desc.add<edm::InputTag>("src", edm::InputTag("hltIter2Phase2L3FromL1TkMuonPixelSeeds", ""));
0304 desc.add<edm::InputTag>("L1TkMu", edm::InputTag("L1TkMuons", "Muon"));
0305 desc.add<edm::FileInPath>("mvaFile_B_0",
0306 edm::FileInPath("RecoMuon/TrackerSeedGenerator/data/xgb_Phase2_Iter2FromL1_barrel_v0.xml"));
0307 desc.add<edm::FileInPath>("mvaFile_E_0",
0308 edm::FileInPath("RecoMuon/TrackerSeedGenerator/data/xgb_Phase2_Iter2FromL1_endcap_v0.xml"));
0309 desc.add<std::vector<double>>("mvaScaleMean_B",
0310 {0.00033113700731766336,
0311 1.6825601468762878e-06,
0312 1.790932122524803e-06,
0313 0.010534608406382916,
0314 0.005969459957330139,
0315 0.0009605022254971113,
0316 0.04384189672781466,
0317 7.846741237608237e-05,
0318 0.40725050850004824,
0319 0.41125151617410227,
0320 0.39815551065544846});
0321 desc.add<std::vector<double>>("mvaScaleStd_B",
0322 {0.0006042948363798624,
0323 2.445644111872427e-06,
0324 3.454992543447134e-06,
0325 0.09401581628887255,
0326 0.7978806947573766,
0327 0.4932933044535928,
0328 0.04180518265631776,
0329 0.058296511682094855,
0330 0.4071857009373577,
0331 0.41337782307392973,
0332 0.4101160349549534});
0333 desc.add<std::vector<double>>("mvaScaleMean_E",
0334 {0.00022658482374555603,
0335 5.358921973784045e-07,
0336 1.010003713549798e-06,
0337 0.0007886873612224615,
0338 0.001197730548842408,
0339 -0.0030252353426003594,
0340 0.07151944804171254,
0341 -0.0006940626775109026,
0342 0.20535152195939896,
0343 0.2966816533783824,
0344 0.28798220230180455});
0345 desc.add<std::vector<double>>("mvaScaleStd_E",
0346 {0.0003857726789049956,
0347 1.4853721474087994e-06,
0348 6.982997036736564e-06,
0349 0.04071340757666084,
0350 0.5897606560095399,
0351 0.33052121398064654,
0352 0.05589386786541949,
0353 0.08806273533388546,
0354 0.3254586902665612,
0355 0.3293354496231377,
0356 0.3179899794578072});
0357 desc.add<bool>("doSort", true);
0358 desc.add<int>("nSeedsMax_B", 20);
0359 desc.add<int>("nSeedsMax_E", 20);
0360 desc.add<double>("mvaCut_B", 0.);
0361 desc.add<double>("mvaCut_E", 0.);
0362 desc.add<double>("etaEdge", 1.2);
0363 desc.add<double>("baseScore", 0.5);
0364 descriptions.add("MuonHLTSeedMVAClassifierPhase2", desc);
0365 }
0366
0367 DEFINE_FWK_MODULE(MuonHLTSeedMVAClassifierPhase2);