TrackerSeedGenerator/plugins/MuonHLTSeedMVAClassifierPhase2.cc

0001 // Package:    RecoMuon_TrackerSeedGenerator
0002 // Class:      MuonHLTSeedMVAClassifierPhase2
0003
0004 // Original Author:  OH Minseok, Sungwon Kim, Won Jun
0005 //         Created:  Mon, 08 Jun 2020 06:20:44 GMT
0006
0007 // system include files
0008 #include <memory>
0009 #include <cmath>
0010
0011 // user include files
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 // Geometry
0023 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0024 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0025
0026 // TrajectorySeed
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 // -- for L1TkMu propagation
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 // class declaration
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   // ----------member data ---------------------------
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 // constructors and destructor
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 // member functions
0141 // ------------ method called on each new Event  ------------
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   // -- sort seeds by MVA score and chooes top nSeedsMax_B_ / nSeedsMax_E_
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   // -- simple fitering based on Mva threshold
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 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
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 //define this as a plug-in
0367 DEFINE_FWK_MODULE(MuonHLTSeedMVAClassifierPhase2);