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::ESHandle<TrackerGeometry> trkGeom = iSetup.getHandle(trackerGeometryESToken_);
0146   edm::ESHandle<GeometricDet> geomDet = iSetup.getHandle(geomDetESToken_);
0147   edm::ESHandle<TrackerTopology> trkTopo = iSetup.getHandle(trackerTopologyESToken_);
0148
0149   GeometricSearchTrackerBuilder builder;
0150   GeometricSearchTracker geomTracker = *(builder.build(&(*geomDet), &(*trkGeom), &(*trkTopo)));
0151
0152   edm::Handle<l1t::TrackerMuonCollection> h_L1TkMu;
0153   bool hasL1TkMu = iEvent.getByToken(t_L1TkMu_, h_L1TkMu);
0154
0155   edm::Handle<TrajectorySeedCollection> h_Seed;
0156   bool hasSeed = iEvent.getByToken(t_Seed_, h_Seed);
0157
0158   edm::ESHandle<MagneticField> magfieldH = iSetup.getHandle(magFieldESToken_);
0159   edm::ESHandle<Propagator> propagatorAlongH = iSetup.getHandle(propagatorESToken_);
0160   std::unique_ptr<Propagator> propagatorAlong = SetPropagationDirection(*propagatorAlongH, alongMomentum);
0161
0162   if (!(hasL1TkMu && hasSeed)) {
0163     edm::LogError("SeedClassifierError") << "Error! Cannot find L1TkMuon or TrajectorySeed\n"
0164                                          << "hasL1TkMu : " << hasL1TkMu << "\n"
0165                                          << "hasSeed : " << hasSeed << "\n";
0166     return;
0167   }
0168
0169   // -- sort seeds by MVA score and chooes top nSeedsMax_B_ / nSeedsMax_E_
0170   if (doSort_) {
0171     std::vector<std::pair<unsigned, double>> pairSeedIdxMvaScore_B = {};
0172     std::vector<std::pair<unsigned, double>> pairSeedIdxMvaScore_E = {};
0173
0174     for (auto i = 0U; i < h_Seed->size(); ++i) {
0175       const auto& seed(h_Seed->at(i));
0176
0177       GlobalVector global_p = trkGeom->idToDet(seed.startingState().detId())
0178                                   ->surface()
0179                                   .toGlobal(seed.startingState().parameters().momentum());
0180       GlobalPoint global_x = trkGeom->idToDet(seed.startingState().detId())
0181                                  ->surface()
0182                                  .toGlobal(seed.startingState().parameters().position());
0183
0184       bool isB = (std::abs(global_p.eta()) < etaEdge_);
0185
0186       if (isB && nSeedsMax_B_ < 0) {
0187         result->emplace_back(seed);
0188         continue;
0189       }
0190
0191       if (!isB && nSeedsMax_E_ < 0) {
0192         result->emplace_back(seed);
0193         continue;
0194       }
0195
0196       double mva = getSeedMva(
0197           mvaEstimator_, seed, global_p, global_x, h_L1TkMu, magfieldH, *(propagatorAlong.get()), geomTracker);
0198
0199       double logistic = 1 / (1 + std::exp(-mva));
0200
0201       if (isB)
0202         pairSeedIdxMvaScore_B.push_back(make_pair(i, logistic));
0203       else
0204         pairSeedIdxMvaScore_E.push_back(make_pair(i, logistic));
0205     }
0206
0207     std::sort(pairSeedIdxMvaScore_B.begin(), pairSeedIdxMvaScore_B.end(), sortByMvaScorePhase2);
0208     std::sort(pairSeedIdxMvaScore_E.begin(), pairSeedIdxMvaScore_E.end(), sortByMvaScorePhase2);
0209
0210     for (auto i = 0U; i < pairSeedIdxMvaScore_B.size(); ++i) {
0211       if ((int)i == nSeedsMax_B_)
0212         break;
0213       const auto& seed(h_Seed->at(pairSeedIdxMvaScore_B.at(i).first));
0214       result->emplace_back(seed);
0215     }
0216
0217     for (auto i = 0U; i < pairSeedIdxMvaScore_E.size(); ++i) {
0218       if ((int)i == nSeedsMax_E_)
0219         break;
0220       const auto& seed(h_Seed->at(pairSeedIdxMvaScore_E.at(i).first));
0221       result->emplace_back(seed);
0222     }
0223   }
0224
0225   // -- simple fitering based on Mva threshold
0226   else {
0227     for (auto i = 0U; i < h_Seed->size(); ++i) {
0228       const auto& seed(h_Seed->at(i));
0229
0230       GlobalVector global_p = trkGeom->idToDet(seed.startingState().detId())
0231                                   ->surface()
0232                                   .toGlobal(seed.startingState().parameters().momentum());
0233       GlobalPoint global_x = trkGeom->idToDet(seed.startingState().detId())
0234                                  ->surface()
0235                                  .toGlobal(seed.startingState().parameters().position());
0236
0237       bool isB = (std::abs(global_p.eta()) < etaEdge_);
0238
0239       if (isB && mvaCut_B_ <= 0.) {
0240         result->emplace_back(seed);
0241         continue;
0242       }
0243
0244       if (!isB && mvaCut_E_ <= 0.) {
0245         result->emplace_back(seed);
0246         continue;
0247       }
0248
0249       double mva = getSeedMva(
0250           mvaEstimator_, seed, global_p, global_x, h_L1TkMu, magfieldH, *(propagatorAlong.get()), geomTracker);
0251
0252       double logistic = 1 / (1 + std::exp(-mva));
0253
0254       bool passMva = ((isB && (logistic > mvaCut_B_)) || (!isB && (logistic > mvaCut_E_)));
0255
0256       if (passMva)
0257         result->emplace_back(seed);
0258     }
0259   }
0260
0261   iEvent.put(std::move(result));
0262 }
0263
0264 double MuonHLTSeedMVAClassifierPhase2::getSeedMva(const pairSeedMvaEstimator& pairMvaEstimator,
0265                                                   const TrajectorySeed& seed,
0266                                                   const GlobalVector& global_p,
0267                                                   const GlobalPoint& global_x,
0268                                                   const edm::Handle<l1t::TrackerMuonCollection>& h_L1TkMu,
0269                                                   const edm::ESHandle<MagneticField>& magfieldH,
0270                                                   const Propagator& propagatorAlong,
0271                                                   const GeometricSearchTracker& geomTracker) {
0272   double mva = 0.;
0273
0274   if (std::abs(global_p.eta()) < etaEdge_) {
0275     mva =
0276         pairMvaEstimator.first->computeMva(seed, global_p, global_x, h_L1TkMu, magfieldH, propagatorAlong, geomTracker);
0277   } else {
0278     mva = pairMvaEstimator.second->computeMva(
0279         seed, global_p, global_x, h_L1TkMu, magfieldH, propagatorAlong, geomTracker);
0280   }
0281
0282   return (baseScore_ + mva);
0283 }
0284
0285 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0286 void MuonHLTSeedMVAClassifierPhase2::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0287   edm::ParameterSetDescription desc;
0288   desc.add<edm::InputTag>("src", edm::InputTag("hltIter2Phase2L3FromL1TkMuonPixelSeeds", ""));
0289   desc.add<edm::InputTag>("L1TkMu", edm::InputTag("L1TkMuons", "Muon"));
0290   desc.add<edm::FileInPath>("mvaFile_B_0",
0291                             edm::FileInPath("RecoMuon/TrackerSeedGenerator/data/xgb_Phase2_Iter2FromL1_barrel_v0.xml"));
0292   desc.add<edm::FileInPath>("mvaFile_E_0",
0293                             edm::FileInPath("RecoMuon/TrackerSeedGenerator/data/xgb_Phase2_Iter2FromL1_endcap_v0.xml"));
0294   desc.add<std::vector<double>>("mvaScaleMean_B",
0295                                 {0.00033113700731766336,
0296                                  1.6825601468762878e-06,
0297                                  1.790932122524803e-06,
0298                                  0.010534608406382916,
0299                                  0.005969459957330139,
0300                                  0.0009605022254971113,
0301                                  0.04384189672781466,
0302                                  7.846741237608237e-05,
0303                                  0.40725050850004824,
0304                                  0.41125151617410227,
0305                                  0.39815551065544846});
0306   desc.add<std::vector<double>>("mvaScaleStd_B",
0307                                 {0.0006042948363798624,
0308                                  2.445644111872427e-06,
0309                                  3.454992543447134e-06,
0310                                  0.09401581628887255,
0311                                  0.7978806947573766,
0312                                  0.4932933044535928,
0313                                  0.04180518265631776,
0314                                  0.058296511682094855,
0315                                  0.4071857009373577,
0316                                  0.41337782307392973,
0317                                  0.4101160349549534});
0318   desc.add<std::vector<double>>("mvaScaleMean_E",
0319                                 {0.00022658482374555603,
0320                                  5.358921973784045e-07,
0321                                  1.010003713549798e-06,
0322                                  0.0007886873612224615,
0323                                  0.001197730548842408,
0324                                  -0.0030252353426003594,
0325                                  0.07151944804171254,
0326                                  -0.0006940626775109026,
0327                                  0.20535152195939896,
0328                                  0.2966816533783824,
0329                                  0.28798220230180455});
0330   desc.add<std::vector<double>>("mvaScaleStd_E",
0331                                 {0.0003857726789049956,
0332                                  1.4853721474087994e-06,
0333                                  6.982997036736564e-06,
0334                                  0.04071340757666084,
0335                                  0.5897606560095399,
0336                                  0.33052121398064654,
0337                                  0.05589386786541949,
0338                                  0.08806273533388546,
0339                                  0.3254586902665612,
0340                                  0.3293354496231377,
0341                                  0.3179899794578072});
0342   desc.add<bool>("doSort", true);
0343   desc.add<int>("nSeedsMax_B", 20);
0344   desc.add<int>("nSeedsMax_E", 20);
0345   desc.add<double>("mvaCut_B", 0.);
0346   desc.add<double>("mvaCut_E", 0.);
0347   desc.add<double>("etaEdge", 1.2);
0348   desc.add<double>("baseScore", 0.5);
0349   descriptions.add("MuonHLTSeedMVAClassifierPhase2", desc);
0350 }
0351 //define this as a plug-in
0352 DEFINE_FWK_MODULE(MuonHLTSeedMVAClassifierPhase2);