diff --git a/HLTrigger/Configuration/python/HLT_75e33/modules/hltIter2Phase2L3FromL1TkMuonCkfTrackCandidates_cfi.py b/HLTrigger/Configuration/python/HLT_75e33/modules/hltIter2Phase2L3FromL1TkMuonCkfTrackCandidates_cfi.py index 1f760e7e18906..9fca8db69ad90 100644 --- a/HLTrigger/Configuration/python/HLT_75e33/modules/hltIter2Phase2L3FromL1TkMuonCkfTrackCandidates_cfi.py +++ b/HLTrigger/Configuration/python/HLT_75e33/modules/hltIter2Phase2L3FromL1TkMuonCkfTrackCandidates_cfi.py @@ -17,6 +17,6 @@ doSeedingRegionRebuilding = cms.bool(False), maxNSeeds = cms.uint32(100000), maxSeedsBeforeCleaning = cms.uint32(1000), - src = cms.InputTag("hltIter2Phase2L3FromL1TkMuonPixelSeeds"), + src = cms.InputTag("hltIter2Phase2L3FromL1TkMuonPixelSeedsFiltered"), useHitsSplitting = cms.bool(False) ) diff --git a/HLTrigger/Configuration/python/HLT_75e33/modules/hltIter2Phase2L3FromL1TkMuonPixelSeedsFiltered_cfi.py b/HLTrigger/Configuration/python/HLT_75e33/modules/hltIter2Phase2L3FromL1TkMuonPixelSeedsFiltered_cfi.py new file mode 100644 index 0000000000000..8d3074e5fdf1d --- /dev/null +++ b/HLTrigger/Configuration/python/HLT_75e33/modules/hltIter2Phase2L3FromL1TkMuonPixelSeedsFiltered_cfi.py @@ -0,0 +1,19 @@ +import FWCore.ParameterSet.Config as cms + +hltIter2Phase2L3FromL1TkMuonPixelSeedsFiltered = cms.EDProducer("MuonHLTSeedMVAClassifierPhase2", + src = cms.InputTag("hltIter2Phase2L3FromL1TkMuonPixelSeeds"), + L1TkMu = cms.InputTag("l1tTkMuonsGmt"), + mvaFile_B_0 = cms.FileInPath("RecoMuon/TrackerSeedGenerator/data/xgb_Phase2_Iter2FromL1_barrel_v0.xml"), + mvaFile_E_0 = cms.FileInPath("RecoMuon/TrackerSeedGenerator/data/xgb_Phase2_Iter2FromL1_endcap_v0.xml"), + mvaScaleMean_B = cms.vdouble(0.00033113700731766336, 1.6825601468762878e-06, 1.790932122524803e-06, 0.010534608406382916, 0.005969459957330139, 0.0009605022254971113, 0.04384189672781466, 7.846741237608237e-05, 0.40725050850004824, 0.41125151617410227, 0.39815551065544846), + mvaScaleStd_B = cms.vdouble(0.0006042948363798624, 2.445644111872427e-06, 3.454992543447134e-06, 0.09401581628887255, 0.7978806947573766, 0.4932933044535928, 0.04180518265631776, 0.058296511682094855, 0.4071857009373577, 0.41337782307392973, 0.4101160349549534), + mvaScaleMean_E = cms.vdouble(0.00022658482374555603, 5.358921973784045e-07, 1.010003713549798e-06, 0.0007886873612224615, 0.001197730548842408, -0.0030252353426003594, 0.07151944804171254, -0.0006940626775109026, 0.20535152195939896, 0.2966816533783824, 0.28798220230180455), + mvaScaleStd_E = cms.vdouble(0.0003857726789049956, 1.4853721474087994e-06, 6.982997036736564e-06, 0.04071340757666084, 0.5897606560095399, 0.33052121398064654, 0.05589386786541949, 0.08806273533388546, 0.3254586902665612, 0.3293354496231377, 0.3179899794578072), + doSort = cms.bool(True), + nSeedsMax_B = cms.int32(20), + nSeedsMax_E = cms.int32(20), + etaEdge = cms.double(1.2), + mvaCut_B = cms.double(0.), + mvaCut_E = cms.double(0.), + baseScore = cms.double(0.5) +) \ No newline at end of file diff --git a/HLTrigger/Configuration/python/HLT_75e33/paths/HLT_IsoMu24_FromL1TkMuon_cfi.py b/HLTrigger/Configuration/python/HLT_75e33/paths/HLT_IsoMu24_FromL1TkMuon_cfi.py index 6e2c3c6252d05..8b40790591a90 100644 --- a/HLTrigger/Configuration/python/HLT_75e33/paths/HLT_IsoMu24_FromL1TkMuon_cfi.py +++ b/HLTrigger/Configuration/python/HLT_75e33/paths/HLT_IsoMu24_FromL1TkMuon_cfi.py @@ -42,6 +42,7 @@ from ..modules.hltIter2Phase2L3FromL1TkMuonPixelHitTriplets_cfi import * from ..modules.hltIter2Phase2L3FromL1TkMuonPixelLayerTriplets_cfi import * from ..modules.hltIter2Phase2L3FromL1TkMuonPixelSeeds_cfi import * +from ..modules.hltIter2Phase2L3FromL1TkMuonPixelSeedsFiltered_cfi import * from ..modules.hltIter2Phase2L3FromL1TkMuonTrackCutClassifier_cfi import * from ..modules.hltIter2Phase2L3FromL1TkMuonTrackSelectionHighPurity_cfi import * from ..modules.hltL2MuonFromL1TkMuonCandidates_cfi import * @@ -165,6 +166,7 @@ hltIter2Phase2L3FromL1TkMuonPixelHitTriplets, hltIter2Phase2L3FromL1TkMuonPixelLayerTriplets, hltIter2Phase2L3FromL1TkMuonPixelSeeds, + hltIter2Phase2L3FromL1TkMuonPixelSeedsFiltered, hltIter2Phase2L3FromL1TkMuonTrackCutClassifier, hltIter2Phase2L3FromL1TkMuonTrackSelectionHighPurity, hltL2MuonFromL1TkMuonCandidates, diff --git a/HLTrigger/Configuration/python/HLT_75e33/paths/HLT_Mu17_TrkIsoVVL_Mu8_TrkIsoVVL_DZ_FromL1TkMuon_cfi.py b/HLTrigger/Configuration/python/HLT_75e33/paths/HLT_Mu17_TrkIsoVVL_Mu8_TrkIsoVVL_DZ_FromL1TkMuon_cfi.py index 5df2b67f573d5..3e155a837a301 100644 --- a/HLTrigger/Configuration/python/HLT_75e33/paths/HLT_Mu17_TrkIsoVVL_Mu8_TrkIsoVVL_DZ_FromL1TkMuon_cfi.py +++ b/HLTrigger/Configuration/python/HLT_75e33/paths/HLT_Mu17_TrkIsoVVL_Mu8_TrkIsoVVL_DZ_FromL1TkMuon_cfi.py @@ -28,6 +28,7 @@ from ..modules.hltIter2Phase2L3FromL1TkMuonPixelHitTriplets_cfi import * from ..modules.hltIter2Phase2L3FromL1TkMuonPixelLayerTriplets_cfi import * from ..modules.hltIter2Phase2L3FromL1TkMuonPixelSeeds_cfi import * +from ..modules.hltIter2Phase2L3FromL1TkMuonPixelSeedsFiltered_cfi import * from ..modules.hltIter2Phase2L3FromL1TkMuonTrackCutClassifier_cfi import * from ..modules.hltIter2Phase2L3FromL1TkMuonTrackSelectionHighPurity_cfi import * from ..modules.hltL2MuonFromL1TkMuonCandidates_cfi import * @@ -127,6 +128,7 @@ hltIter2Phase2L3FromL1TkMuonPixelHitTriplets, hltIter2Phase2L3FromL1TkMuonPixelLayerTriplets, hltIter2Phase2L3FromL1TkMuonPixelSeeds, + hltIter2Phase2L3FromL1TkMuonPixelSeedsFiltered, hltIter2Phase2L3FromL1TkMuonTrackCutClassifier, hltIter2Phase2L3FromL1TkMuonTrackSelectionHighPurity, hltL2MuonFromL1TkMuonCandidates, diff --git a/HLTrigger/Configuration/python/HLT_75e33/paths/HLT_Mu37_Mu27_FromL1TkMuon_cfi.py b/HLTrigger/Configuration/python/HLT_75e33/paths/HLT_Mu37_Mu27_FromL1TkMuon_cfi.py index 173c303fba3e0..b1d2826a2fb59 100644 --- a/HLTrigger/Configuration/python/HLT_75e33/paths/HLT_Mu37_Mu27_FromL1TkMuon_cfi.py +++ b/HLTrigger/Configuration/python/HLT_75e33/paths/HLT_Mu37_Mu27_FromL1TkMuon_cfi.py @@ -26,6 +26,7 @@ from ..modules.hltIter2Phase2L3FromL1TkMuonPixelHitTriplets_cfi import * from ..modules.hltIter2Phase2L3FromL1TkMuonPixelLayerTriplets_cfi import * from ..modules.hltIter2Phase2L3FromL1TkMuonPixelSeeds_cfi import * +from ..modules.hltIter2Phase2L3FromL1TkMuonPixelSeedsFiltered_cfi import * from ..modules.hltIter2Phase2L3FromL1TkMuonTrackCutClassifier_cfi import * from ..modules.hltIter2Phase2L3FromL1TkMuonTrackSelectionHighPurity_cfi import * from ..modules.hltL2MuonFromL1TkMuonCandidates_cfi import * @@ -99,6 +100,7 @@ hltIter2Phase2L3FromL1TkMuonPixelHitTriplets, hltIter2Phase2L3FromL1TkMuonPixelLayerTriplets, hltIter2Phase2L3FromL1TkMuonPixelSeeds, + hltIter2Phase2L3FromL1TkMuonPixelSeedsFiltered, hltIter2Phase2L3FromL1TkMuonTrackCutClassifier, hltIter2Phase2L3FromL1TkMuonTrackSelectionHighPurity, hltL2MuonFromL1TkMuonCandidates, diff --git a/HLTrigger/Configuration/python/HLT_75e33/paths/HLT_Mu50_FromL1TkMuon_cfi.py b/HLTrigger/Configuration/python/HLT_75e33/paths/HLT_Mu50_FromL1TkMuon_cfi.py index b1e0e2e182f30..8a2c891741c53 100644 --- a/HLTrigger/Configuration/python/HLT_75e33/paths/HLT_Mu50_FromL1TkMuon_cfi.py +++ b/HLTrigger/Configuration/python/HLT_75e33/paths/HLT_Mu50_FromL1TkMuon_cfi.py @@ -23,6 +23,7 @@ from ..modules.hltIter2Phase2L3FromL1TkMuonPixelHitTriplets_cfi import * from ..modules.hltIter2Phase2L3FromL1TkMuonPixelLayerTriplets_cfi import * from ..modules.hltIter2Phase2L3FromL1TkMuonPixelSeeds_cfi import * +from ..modules.hltIter2Phase2L3FromL1TkMuonPixelSeedsFiltered_cfi import * from ..modules.hltIter2Phase2L3FromL1TkMuonTrackCutClassifier_cfi import * from ..modules.hltIter2Phase2L3FromL1TkMuonTrackSelectionHighPurity_cfi import * from ..modules.hltL2MuonFromL1TkMuonCandidates_cfi import * @@ -93,6 +94,7 @@ hltIter2Phase2L3FromL1TkMuonPixelHitTriplets, hltIter2Phase2L3FromL1TkMuonPixelLayerTriplets, hltIter2Phase2L3FromL1TkMuonPixelSeeds, + hltIter2Phase2L3FromL1TkMuonPixelSeedsFiltered, hltIter2Phase2L3FromL1TkMuonTrackCutClassifier, hltIter2Phase2L3FromL1TkMuonTrackSelectionHighPurity, hltL2MuonFromL1TkMuonCandidates, diff --git a/HLTrigger/Configuration/python/HLT_75e33/paths/HLT_TriMu_10_5_5_DZ_FromL1TkMuon_cfi.py b/HLTrigger/Configuration/python/HLT_75e33/paths/HLT_TriMu_10_5_5_DZ_FromL1TkMuon_cfi.py index 9cb6cf27bb6d2..6fc1501f2f987 100644 --- a/HLTrigger/Configuration/python/HLT_75e33/paths/HLT_TriMu_10_5_5_DZ_FromL1TkMuon_cfi.py +++ b/HLTrigger/Configuration/python/HLT_75e33/paths/HLT_TriMu_10_5_5_DZ_FromL1TkMuon_cfi.py @@ -28,6 +28,7 @@ from ..modules.hltIter2Phase2L3FromL1TkMuonPixelHitTriplets_cfi import * from ..modules.hltIter2Phase2L3FromL1TkMuonPixelLayerTriplets_cfi import * from ..modules.hltIter2Phase2L3FromL1TkMuonPixelSeeds_cfi import * +from ..modules.hltIter2Phase2L3FromL1TkMuonPixelSeedsFiltered_cfi import * from ..modules.hltIter2Phase2L3FromL1TkMuonTrackCutClassifier_cfi import * from ..modules.hltIter2Phase2L3FromL1TkMuonTrackSelectionHighPurity_cfi import * from ..modules.hltL2MuonFromL1TkMuonCandidates_cfi import * @@ -103,6 +104,7 @@ hltIter2Phase2L3FromL1TkMuonPixelHitTriplets, hltIter2Phase2L3FromL1TkMuonPixelLayerTriplets, hltIter2Phase2L3FromL1TkMuonPixelSeeds, + hltIter2Phase2L3FromL1TkMuonPixelSeedsFiltered, hltIter2Phase2L3FromL1TkMuonTrackCutClassifier, hltIter2Phase2L3FromL1TkMuonTrackSelectionHighPurity, hltL2MuonFromL1TkMuonCandidates, diff --git a/RecoMuon/TrackerSeedGenerator/BuildFile.xml b/RecoMuon/TrackerSeedGenerator/BuildFile.xml index 485ac85a0dd2e..b151d85bd4ce5 100644 --- a/RecoMuon/TrackerSeedGenerator/BuildFile.xml +++ b/RecoMuon/TrackerSeedGenerator/BuildFile.xml @@ -1,32 +1,60 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/RecoMuon/TrackerSeedGenerator/interface/SeedMvaEstimatorPhase2.h b/RecoMuon/TrackerSeedGenerator/interface/SeedMvaEstimatorPhase2.h new file mode 100644 index 0000000000000..6dcce3a799ad2 --- /dev/null +++ b/RecoMuon/TrackerSeedGenerator/interface/SeedMvaEstimatorPhase2.h @@ -0,0 +1,88 @@ +#ifndef RecoMuon_TrackerSeedGenerator_SeedMvaEstimatorPhase2_h +#define RecoMuon_TrackerSeedGenerator_SeedMvaEstimatorPhase2_h + +#include "FWCore/Framework/interface/ConsumesCollector.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Utilities/interface/ESGetToken.h" +#include "FWCore/Utilities/interface/ESInputTag.h" +#include "DataFormats/Common/interface/Handle.h" +#include "DataFormats/TrackerCommon/interface/TrackerTopology.h" +#include "DataFormats/TrajectorySeed/interface/TrajectorySeed.h" +#include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h" +#include "DataFormats/TrajectorySeed/interface/PropagationDirection.h" +#include "DataFormats/TrajectoryState/interface/PTrajectoryStateOnDet.h" +#include "DataFormats/TrajectoryState/interface/LocalTrajectoryParameters.h" +#include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h" +#include "DataFormats/RecoCandidate/interface/RecoChargedCandidateFwd.h" +#include "DataFormats/L1TCorrelator/interface/TkMuon.h" +#include "DataFormats/L1TCorrelator/interface/TkMuonFwd.h" +#include "DataFormats/L1TMuonPhase2/interface/TrackerMuon.h" +#include "TrackingTools/GeomPropagators/interface/Propagator.h" +#include "TrackingTools/Records/interface/TrackingComponentsRecord.h" +#include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h" +#include "RecoTracker/TkDetLayers/interface/GeometricSearchTracker.h" +#include "RecoTracker/TkDetLayers/interface/GeometricSearchTrackerBuilder.h" +#include "Geometry/TrackerNumberingBuilder/interface/GeometricDet.h" +#include "Geometry/Records/interface/IdealGeometryRecord.h" +#include "Geometry/Records/interface/TrackerTopologyRcd.h" + +#include +#include + +typedef pair LayerTSOS; +typedef pair LayerHit; + +class GBRForest; + +namespace edm { + class FileInPath; +} + +class SeedMvaEstimatorPhase2 { +public: + SeedMvaEstimatorPhase2(const edm::FileInPath& weightsfile, + const std::vector& scale_mean, + const std::vector& scale_std); + ~SeedMvaEstimatorPhase2(); + + double computeMva(const TrajectorySeed&, + const GlobalVector&, + const GlobalPoint&, + const edm::Handle&, + const edm::ESHandle&, + const Propagator&, + const GeometricSearchTracker&) const; + +private: + std::unique_ptr gbrForest_; + const std::vector scale_mean_; + const std::vector scale_std_; + + vector getTsosOnPixels(const TTTrack&, + const edm::ESHandle&, + const Propagator&, + const GeometricSearchTracker&) const; + + vector > getHitTsosPairs(const TrajectorySeed&, + const edm::Handle&, + const edm::ESHandle&, + const Propagator&, + const GeometricSearchTracker&) const; + + void getL1TTVariables(const TrajectorySeed&, + const GlobalVector&, + const GlobalPoint&, + const edm::Handle&, + float&, + float&) const; + void getHitL1TkVariables(const TrajectorySeed&, + const edm::Handle&, + const edm::ESHandle&, + const Propagator&, + const GeometricSearchTracker&, + float&, + float&, + float&) const; +}; +#endif diff --git a/RecoMuon/TrackerSeedGenerator/plugins/BuildFile.xml b/RecoMuon/TrackerSeedGenerator/plugins/BuildFile.xml index 5528edd4e9959..7a47fd226d191 100644 --- a/RecoMuon/TrackerSeedGenerator/plugins/BuildFile.xml +++ b/RecoMuon/TrackerSeedGenerator/plugins/BuildFile.xml @@ -1,7 +1,17 @@ + + + + + + + + + + @@ -10,30 +20,42 @@ - + + + + + - + - + + + + + + + + + diff --git a/RecoMuon/TrackerSeedGenerator/plugins/MuonHLTSeedMVAClassifierPhase2.cc b/RecoMuon/TrackerSeedGenerator/plugins/MuonHLTSeedMVAClassifierPhase2.cc new file mode 100644 index 0000000000000..bb0ecd31ea163 --- /dev/null +++ b/RecoMuon/TrackerSeedGenerator/plugins/MuonHLTSeedMVAClassifierPhase2.cc @@ -0,0 +1,352 @@ +// Package: RecoMuon_TrackerSeedGenerator +// Class: MuonHLTSeedMVAClassifierPhase2 + +// Original Author: OH Minseok, Sungwon Kim, Won Jun +// Created: Mon, 08 Jun 2020 06:20:44 GMT + +// system include files +#include +#include + +// user include files +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "FWCore/Framework/interface/Event.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 "FWCore/Utilities/interface/ESGetToken.h" +#include "FWCore/Utilities/interface/ESInputTag.h" + +// Geometry +#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h" +#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" + +// TrajectorySeed +#include "DataFormats/TrajectorySeed/interface/TrajectorySeed.h" +#include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h" +#include "DataFormats/TrajectorySeed/interface/PropagationDirection.h" +#include "DataFormats/TrajectoryState/interface/PTrajectoryStateOnDet.h" +#include "DataFormats/TrajectoryState/interface/LocalTrajectoryParameters.h" +#include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h" + +// -- for L1TkMu propagation +#include "TrackingTools/GeomPropagators/interface/Propagator.h" +#include "TrackingTools/Records/interface/TrackingComponentsRecord.h" +#include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h" +#include "RecoTracker/TkDetLayers/interface/GeometricSearchTracker.h" +#include "RecoTracker/TkDetLayers/interface/GeometricSearchTrackerBuilder.h" +#include "Geometry/TrackerNumberingBuilder/interface/GeometricDet.h" +#include "Geometry/Records/interface/IdealGeometryRecord.h" +#include "DataFormats/TrackerCommon/interface/TrackerTopology.h" +#include "Geometry/Records/interface/TrackerTopologyRcd.h" + +#include "RecoMuon/TrackerSeedGenerator/interface/SeedMvaEstimatorPhase2.h" + +// class declaration +bool sortByMvaScorePhase2(const std::pair& A, const std::pair& B) { + return (A.second > B.second); +}; + +class MuonHLTSeedMVAClassifierPhase2 : public edm::stream::EDProducer<> { +public: + explicit MuonHLTSeedMVAClassifierPhase2(const edm::ParameterSet&); + ~MuonHLTSeedMVAClassifierPhase2() override = default; + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + +private: + void produce(edm::Event&, const edm::EventSetup&) override; + + // ----------member data --------------------------- + const edm::EDGetTokenT t_Seed_; + const edm::EDGetTokenT t_L1TkMu_; + + typedef std::pair, std::unique_ptr> + pairSeedMvaEstimator; + pairSeedMvaEstimator mvaEstimator_; + + const edm::FileInPath mvaFile_B_0_; + const edm::FileInPath mvaFile_E_0_; + + const std::vector mvaScaleMean_B_; + const std::vector mvaScaleStd_B_; + const std::vector mvaScaleMean_E_; + const std::vector mvaScaleStd_E_; + + const double etaEdge_; + const double mvaCut_B_; + const double mvaCut_E_; + + const bool doSort_; + const int nSeedsMax_B_; + const int nSeedsMax_E_; + + const double baseScore_; + + const edm::ESGetToken trackerTopologyESToken_; + const edm::ESGetToken trackerGeometryESToken_; + const edm::ESGetToken magFieldESToken_; + const edm::ESGetToken geomDetESToken_; + const edm::ESGetToken propagatorESToken_; + + double getSeedMva(const pairSeedMvaEstimator& pairMvaEstimator, + const TrajectorySeed& seed, + const GlobalVector& global_p, + const GlobalPoint& global_x, + const edm::Handle& h_L1TkMu, + const edm::ESHandle& magfieldH, + const Propagator& propagatorAlong, + const GeometricSearchTracker& geomTracker); +}; + +// constructors and destructor +MuonHLTSeedMVAClassifierPhase2::MuonHLTSeedMVAClassifierPhase2(const edm::ParameterSet& iConfig) + : t_Seed_(consumes(iConfig.getParameter("src"))), + t_L1TkMu_(consumes(iConfig.getParameter("L1TkMu"))), + + mvaFile_B_0_(iConfig.getParameter("mvaFile_B_0")), + mvaFile_E_0_(iConfig.getParameter("mvaFile_E_0")), + + mvaScaleMean_B_(iConfig.getParameter>("mvaScaleMean_B")), + mvaScaleStd_B_(iConfig.getParameter>("mvaScaleStd_B")), + mvaScaleMean_E_(iConfig.getParameter>("mvaScaleMean_E")), + mvaScaleStd_E_(iConfig.getParameter>("mvaScaleStd_E")), + + etaEdge_(iConfig.getParameter("etaEdge")), + mvaCut_B_(iConfig.getParameter("mvaCut_B")), + mvaCut_E_(iConfig.getParameter("mvaCut_E")), + + doSort_(iConfig.getParameter("doSort")), + nSeedsMax_B_(iConfig.getParameter("nSeedsMax_B")), + nSeedsMax_E_(iConfig.getParameter("nSeedsMax_E")), + + baseScore_(iConfig.getParameter("baseScore")), + + trackerTopologyESToken_(esConsumes()), + trackerGeometryESToken_(esConsumes()), + magFieldESToken_(esConsumes()), + geomDetESToken_(esConsumes()), + propagatorESToken_( + esConsumes(edm::ESInputTag("", "PropagatorWithMaterialParabolicMf"))) { + produces(); + + mvaEstimator_ = + std::make_pair(std::make_unique(mvaFile_B_0_, mvaScaleMean_B_, mvaScaleStd_B_), + std::make_unique(mvaFile_E_0_, mvaScaleMean_E_, mvaScaleStd_E_)); +} + +// member functions +// ------------ method called on each new Event ------------ +void MuonHLTSeedMVAClassifierPhase2::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { + auto result = std::make_unique(); + + edm::ESHandle trkGeom = iSetup.getHandle(trackerGeometryESToken_); + edm::ESHandle geomDet = iSetup.getHandle(geomDetESToken_); + edm::ESHandle trkTopo = iSetup.getHandle(trackerTopologyESToken_); + + GeometricSearchTrackerBuilder builder; + GeometricSearchTracker geomTracker = *(builder.build(&(*geomDet), &(*trkGeom), &(*trkTopo))); + + edm::Handle h_L1TkMu; + bool hasL1TkMu = iEvent.getByToken(t_L1TkMu_, h_L1TkMu); + + edm::Handle h_Seed; + bool hasSeed = iEvent.getByToken(t_Seed_, h_Seed); + + edm::ESHandle magfieldH = iSetup.getHandle(magFieldESToken_); + edm::ESHandle propagatorAlongH = iSetup.getHandle(propagatorESToken_); + std::unique_ptr propagatorAlong = SetPropagationDirection(*propagatorAlongH, alongMomentum); + + if (!(hasL1TkMu && hasSeed)) { + edm::LogError("SeedClassifierError") << "Error! Cannot find L1TkMuon or TrajectorySeed\n" + << "hasL1TkMu : " << hasL1TkMu << "\n" + << "hasSeed : " << hasSeed << "\n"; + return; + } + + // -- sort seeds by MVA score and chooes top nSeedsMax_B_ / nSeedsMax_E_ + if (doSort_) { + std::vector> pairSeedIdxMvaScore_B = {}; + std::vector> pairSeedIdxMvaScore_E = {}; + + for (auto i = 0U; i < h_Seed->size(); ++i) { + const auto& seed(h_Seed->at(i)); + + GlobalVector global_p = trkGeom->idToDet(seed.startingState().detId()) + ->surface() + .toGlobal(seed.startingState().parameters().momentum()); + GlobalPoint global_x = trkGeom->idToDet(seed.startingState().detId()) + ->surface() + .toGlobal(seed.startingState().parameters().position()); + + bool isB = (std::abs(global_p.eta()) < etaEdge_); + + if (isB && nSeedsMax_B_ < 0) { + result->emplace_back(seed); + continue; + } + + if (!isB && nSeedsMax_E_ < 0) { + result->emplace_back(seed); + continue; + } + + double mva = getSeedMva( + mvaEstimator_, seed, global_p, global_x, h_L1TkMu, magfieldH, *(propagatorAlong.get()), geomTracker); + + double logistic = 1 / (1 + std::exp(-mva)); + + if (isB) + pairSeedIdxMvaScore_B.push_back(make_pair(i, logistic)); + else + pairSeedIdxMvaScore_E.push_back(make_pair(i, logistic)); + } + + std::sort(pairSeedIdxMvaScore_B.begin(), pairSeedIdxMvaScore_B.end(), sortByMvaScorePhase2); + std::sort(pairSeedIdxMvaScore_E.begin(), pairSeedIdxMvaScore_E.end(), sortByMvaScorePhase2); + + for (auto i = 0U; i < pairSeedIdxMvaScore_B.size(); ++i) { + if ((int)i == nSeedsMax_B_) + break; + const auto& seed(h_Seed->at(pairSeedIdxMvaScore_B.at(i).first)); + result->emplace_back(seed); + } + + for (auto i = 0U; i < pairSeedIdxMvaScore_E.size(); ++i) { + if ((int)i == nSeedsMax_E_) + break; + const auto& seed(h_Seed->at(pairSeedIdxMvaScore_E.at(i).first)); + result->emplace_back(seed); + } + } + + // -- simple fitering based on Mva threshold + else { + for (auto i = 0U; i < h_Seed->size(); ++i) { + const auto& seed(h_Seed->at(i)); + + GlobalVector global_p = trkGeom->idToDet(seed.startingState().detId()) + ->surface() + .toGlobal(seed.startingState().parameters().momentum()); + GlobalPoint global_x = trkGeom->idToDet(seed.startingState().detId()) + ->surface() + .toGlobal(seed.startingState().parameters().position()); + + bool isB = (std::abs(global_p.eta()) < etaEdge_); + + if (isB && mvaCut_B_ <= 0.) { + result->emplace_back(seed); + continue; + } + + if (!isB && mvaCut_E_ <= 0.) { + result->emplace_back(seed); + continue; + } + + double mva = getSeedMva( + mvaEstimator_, seed, global_p, global_x, h_L1TkMu, magfieldH, *(propagatorAlong.get()), geomTracker); + + double logistic = 1 / (1 + std::exp(-mva)); + + bool passMva = ((isB && (logistic > mvaCut_B_)) || (!isB && (logistic > mvaCut_E_))); + + if (passMva) + result->emplace_back(seed); + } + } + + iEvent.put(std::move(result)); +} + +double MuonHLTSeedMVAClassifierPhase2::getSeedMva(const pairSeedMvaEstimator& pairMvaEstimator, + const TrajectorySeed& seed, + const GlobalVector& global_p, + const GlobalPoint& global_x, + const edm::Handle& h_L1TkMu, + const edm::ESHandle& magfieldH, + const Propagator& propagatorAlong, + const GeometricSearchTracker& geomTracker) { + double mva = 0.; + + if (std::abs(global_p.eta()) < etaEdge_) { + mva = + pairMvaEstimator.first->computeMva(seed, global_p, global_x, h_L1TkMu, magfieldH, propagatorAlong, geomTracker); + } else { + mva = pairMvaEstimator.second->computeMva( + seed, global_p, global_x, h_L1TkMu, magfieldH, propagatorAlong, geomTracker); + } + + return (baseScore_ + mva); +} + +// ------------ method fills 'descriptions' with the allowed parameters for the module ------------ +void MuonHLTSeedMVAClassifierPhase2::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("src", edm::InputTag("hltIter2Phase2L3FromL1TkMuonPixelSeeds", "")); + desc.add("L1TkMu", edm::InputTag("L1TkMuons", "Muon")); + desc.add("mvaFile_B_0", + edm::FileInPath("RecoMuon/TrackerSeedGenerator/data/xgb_Phase2_Iter2FromL1_barrel_v0.xml")); + desc.add("mvaFile_E_0", + edm::FileInPath("RecoMuon/TrackerSeedGenerator/data/xgb_Phase2_Iter2FromL1_endcap_v0.xml")); + desc.add>("mvaScaleMean_B", + {0.00033113700731766336, + 1.6825601468762878e-06, + 1.790932122524803e-06, + 0.010534608406382916, + 0.005969459957330139, + 0.0009605022254971113, + 0.04384189672781466, + 7.846741237608237e-05, + 0.40725050850004824, + 0.41125151617410227, + 0.39815551065544846}); + desc.add>("mvaScaleStd_B", + {0.0006042948363798624, + 2.445644111872427e-06, + 3.454992543447134e-06, + 0.09401581628887255, + 0.7978806947573766, + 0.4932933044535928, + 0.04180518265631776, + 0.058296511682094855, + 0.4071857009373577, + 0.41337782307392973, + 0.4101160349549534}); + desc.add>("mvaScaleMean_E", + {0.00022658482374555603, + 5.358921973784045e-07, + 1.010003713549798e-06, + 0.0007886873612224615, + 0.001197730548842408, + -0.0030252353426003594, + 0.07151944804171254, + -0.0006940626775109026, + 0.20535152195939896, + 0.2966816533783824, + 0.28798220230180455}); + desc.add>("mvaScaleStd_E", + {0.0003857726789049956, + 1.4853721474087994e-06, + 6.982997036736564e-06, + 0.04071340757666084, + 0.5897606560095399, + 0.33052121398064654, + 0.05589386786541949, + 0.08806273533388546, + 0.3254586902665612, + 0.3293354496231377, + 0.3179899794578072}); + desc.add("doSort", true); + desc.add("nSeedsMax_B", 20); + desc.add("nSeedsMax_E", 20); + desc.add("mvaCut_B", 0.); + desc.add("mvaCut_E", 0.); + desc.add("etaEdge", 1.2); + desc.add("baseScore", 0.5); + descriptions.add("MuonHLTSeedMVAClassifierPhase2", desc); +} +//define this as a plug-in +DEFINE_FWK_MODULE(MuonHLTSeedMVAClassifierPhase2); \ No newline at end of file diff --git a/RecoMuon/TrackerSeedGenerator/src/SeedMvaEstimatorPhase2.cc b/RecoMuon/TrackerSeedGenerator/src/SeedMvaEstimatorPhase2.cc new file mode 100644 index 0000000000000..1571c90c12f65 --- /dev/null +++ b/RecoMuon/TrackerSeedGenerator/src/SeedMvaEstimatorPhase2.cc @@ -0,0 +1,309 @@ +#include "RecoMuon/TrackerSeedGenerator/interface/SeedMvaEstimatorPhase2.h" + +#include "DataFormats/Math/interface/deltaPhi.h" +#include "DataFormats/Math/interface/deltaR.h" + +#include "FWCore/ParameterSet/interface/FileInPath.h" +#include "CommonTools/MVAUtils/interface/GBRForestTools.h" + +#include + +using namespace std; + +namespace Phase2 { + enum inputIndexesPhase2 { + kTsosErr0, // 0 + kTsosErr2, // 1 + kTsosErr5, // 2 + kTsosDxdz, // 3 + kTsosDydz, // 4 + kTsosQbp, // 5 + kDRdRL1TkMuSeedP, // 6 + kDRdPhiL1TkMuSeedP, // 7 + kExpd2HitL1Tk1, // 8 + kExpd2HitL1Tk2, // 9 + kExpd2HitL1Tk3, // 10 + kLast // 11 + }; +} + +SeedMvaEstimatorPhase2::SeedMvaEstimatorPhase2(const edm::FileInPath& weightsfile, + const std::vector& scale_mean, + const std::vector& scale_std) + : scale_mean_(scale_mean), scale_std_(scale_std) { + gbrForest_ = createGBRForest(weightsfile); +} + +SeedMvaEstimatorPhase2::~SeedMvaEstimatorPhase2() {} + +void SeedMvaEstimatorPhase2::getL1TTVariables(const TrajectorySeed& seed, + const GlobalVector& global_p, + const GlobalPoint& global_x, + const edm::Handle& h_L1TkMu, + float& DRL1TkMu, + float& DPhiL1TkMu) const { + for (auto L1TkMu = h_L1TkMu->begin(); L1TkMu != h_L1TkMu->end(); ++L1TkMu) { + auto TkRef = L1TkMu->trkPtr(); + float DRL1TkMu_tmp = reco::deltaR(*TkRef, global_p); + float DPhiL1TkMu_tmp = reco::deltaPhi(TkRef->phi(), global_p.phi()); + if (DRL1TkMu_tmp < DRL1TkMu) { + DRL1TkMu = DRL1TkMu_tmp; + DPhiL1TkMu = DPhiL1TkMu_tmp; + } + } +} + +vector SeedMvaEstimatorPhase2::getTsosOnPixels(const TTTrack& l1tk, + const edm::ESHandle& magfieldH, + const Propagator& propagatorAlong, + const GeometricSearchTracker& geomTracker) const { + vector v_tsos = {}; + + std::vector const& bpix = geomTracker.pixelBarrelLayers(); + std::vector const& nfpix = geomTracker.negPixelForwardLayers(); + std::vector const& pfpix = geomTracker.posPixelForwardLayers(); + + int chargeTk = l1tk.rInv() > 0. ? 1 : -1; + GlobalPoint gpos = l1tk.POCA(); + GlobalVector gmom = l1tk.momentum(); + + FreeTrajectoryState fts = FreeTrajectoryState(gpos, gmom, chargeTk, magfieldH.product()); + + for (auto it = bpix.begin(); it != bpix.end(); ++it) { + TrajectoryStateOnSurface tsos = propagatorAlong.propagate(fts, (**it).specificSurface()); + if (!tsos.isValid()) + continue; + + auto z0 = std::abs(tsos.globalPosition().z() - (**it).specificSurface().position().z()); + auto deltaZ = 0.5f * (**it).surface().bounds().length(); + deltaZ += 0.5f * (**it).surface().bounds().thickness() * std::abs(tsos.globalDirection().z()) / + tsos.globalDirection().perp(); + bool is_compatible = (z0 < deltaZ); + + if (is_compatible) { + v_tsos.push_back(make_pair((const DetLayer*)(*it), tsos)); + } + } + for (auto it = nfpix.begin(); it != nfpix.end(); ++it) { + TrajectoryStateOnSurface tsos = propagatorAlong.propagate(fts, (**it).specificSurface()); + if (!tsos.isValid()) + continue; + + auto r2 = tsos.localPosition().perp2(); + float deltaR = 0.5f * (**it).surface().bounds().thickness() * tsos.localDirection().perp() / + std::abs(tsos.localDirection().z()); + auto ri2 = std::max((**it).specificSurface().innerRadius() - deltaR, 0.f); + ri2 *= ri2; + auto ro2 = (**it).specificSurface().outerRadius() + deltaR; + ro2 *= ro2; + bool is_compatible = ((r2 > ri2) && (r2 < ro2)); + + if (is_compatible) { + v_tsos.push_back(make_pair((const DetLayer*)(*it), tsos)); + } + } + for (auto it = pfpix.begin(); it != pfpix.end(); ++it) { + TrajectoryStateOnSurface tsos = propagatorAlong.propagate(fts, (**it).specificSurface()); + if (!tsos.isValid()) + continue; + + auto r2 = tsos.localPosition().perp2(); + float deltaR = 0.5f * (**it).surface().bounds().thickness() * tsos.localDirection().perp() / + std::abs(tsos.localDirection().z()); + auto ri2 = std::max((**it).specificSurface().innerRadius() - deltaR, 0.f); + ri2 *= ri2; + auto ro2 = (**it).specificSurface().outerRadius() + deltaR; + ro2 *= ro2; + bool is_compatible = ((r2 > ri2) && (r2 < ro2)); + if (is_compatible) { + v_tsos.push_back(make_pair((const DetLayer*)(*it), tsos)); + } + } + return v_tsos; +} + +// FIX HERE +vector > SeedMvaEstimatorPhase2::getHitTsosPairs( + const TrajectorySeed& seed, + const edm::Handle& L1TkMuonHandle, + const edm::ESHandle& magfieldH, + const Propagator& propagatorAlong, + const GeometricSearchTracker& geomTracker) const { + vector > out = {}; + float av_dr_min = 999.; + + // -- loop on L1TkMu + for (auto L1TkMu = L1TkMuonHandle->begin(); L1TkMu != L1TkMuonHandle->end(); ++L1TkMu) { + const vector v_tsos = getTsosOnPixels(*L1TkMu->trkPtr(), magfieldH, propagatorAlong, geomTracker); + + // -- loop on recHits + vector v_tsos_skip(v_tsos.size(), 0); + vector > hitTsosPair = {}; + int ihit = 0; + for (const auto& hit : seed.recHits()) { + // -- look for closest tsos by absolute distance + int the_tsos = -99999; + float dr_min = 20.; + for (auto i = 0U; i < v_tsos.size(); ++i) { + if (v_tsos_skip.at(i)) + continue; + float dr = (v_tsos.at(i).second.globalPosition() - hit.globalPosition()).mag(); + if (dr < dr_min) { + dr_min = dr; + the_tsos = i; + v_tsos_skip.at(i) = 1; + } + } + if (the_tsos > -1) { + const DetLayer* thelayer = geomTracker.idToLayer(hit.geographicalId()); + hitTsosPair.push_back(make_pair(make_pair(thelayer, &hit), v_tsos.at(the_tsos))); + } + ihit++; + } // loop on recHits + + // -- find tsos for all recHits? + if ((int)hitTsosPair.size() == ihit) { + float av_dr = 0.; + for (auto it = hitTsosPair.begin(); it != hitTsosPair.end(); ++it) { + auto hit = it->first.second; + auto tsos = it->second.second; + av_dr += (hit->globalPosition() - tsos.globalPosition()).mag(); + } + av_dr = av_dr > 0 ? av_dr / (float)hitTsosPair.size() : 1.e6; + if (av_dr < av_dr_min) { + av_dr_min = av_dr; + out = hitTsosPair; + } + } + } // loop on L1TkMu + return out; +} + +void SeedMvaEstimatorPhase2::getHitL1TkVariables(const TrajectorySeed& seed, + const edm::Handle& L1TkMuonHandle, + const edm::ESHandle& magfieldH, + const Propagator& propagatorAlong, + const GeometricSearchTracker& geomTracker, + float& expd2HitL1Tk1, + float& expd2HitL1Tk2, + float& expd2HitL1Tk3) const { + vector > hitTsosPair = + getHitTsosPairs(seed, L1TkMuonHandle, magfieldH, propagatorAlong, geomTracker); + + bool found = (!hitTsosPair.empty()); + + if (found) { + float l1x1 = 99999.; + float l1y1 = 99999.; + float l1z1 = 99999.; + float hitx1 = 99999.; + float hity1 = 99999.; + float hitz1 = 99999.; + + float l1x2 = 99999.; + float l1y2 = 99999.; + float l1z2 = 99999.; + float hitx2 = 99999.; + float hity2 = 99999.; + float hitz2 = 99999.; + + float l1x3 = 99999.; + float l1y3 = 99999.; + float l1z3 = 99999.; + float hitx3 = 99999.; + float hity3 = 99999.; + float hitz3 = 99999.; + + if (hitTsosPair.size() > 1) { + auto hit1 = hitTsosPair.at(0).first.second; + auto tsos1 = hitTsosPair.at(0).second.second; + + l1x1 = tsos1.globalPosition().x(); + l1y1 = tsos1.globalPosition().y(); + l1z1 = tsos1.globalPosition().z(); + hitx1 = hit1->globalPosition().x(); + hity1 = hit1->globalPosition().y(); + hitz1 = hit1->globalPosition().z(); + + auto hit2 = hitTsosPair.at(1).first.second; + auto tsos2 = hitTsosPair.at(1).second.second; + + l1x2 = tsos2.globalPosition().x(); + l1y2 = tsos2.globalPosition().y(); + l1z2 = tsos2.globalPosition().z(); + hitx2 = hit2->globalPosition().x(); + hity2 = hit2->globalPosition().y(); + hitz2 = hit2->globalPosition().z(); + } + if (hitTsosPair.size() > 2) { + auto hit3 = hitTsosPair.at(2).first.second; + auto tsos3 = hitTsosPair.at(2).second.second; + + l1x3 = tsos3.globalPosition().x(); + l1y3 = tsos3.globalPosition().y(); + l1z3 = tsos3.globalPosition().z(); + hitx3 = hit3->globalPosition().x(); + hity3 = hit3->globalPosition().y(); + hitz3 = hit3->globalPosition().z(); + } + + // If hit == 99999. >> cannot find hit info >> distance btw hit~tsos is large >> expd2HitL1Tk ~ 0 + if (hitx1 != 99999.) { + expd2HitL1Tk1 = exp(-(pow((l1x1 - hitx1), 2) + pow((l1y1 - hity1), 2) + pow((l1z1 - hitz1), 2))); + } else { + expd2HitL1Tk1 = 0.; + } + + if (hitx2 != 99999.) { + expd2HitL1Tk2 = exp(-(pow((l1x2 - hitx2), 2) + pow((l1y2 - hity2), 2) + pow((l1z2 - hitz2), 2))); + } else { + expd2HitL1Tk2 = 0.; + } + + if (hitx3 != 99999.) { + expd2HitL1Tk3 = exp(-(pow((l1x3 - hitx3), 2) + pow((l1y3 - hity3), 2) + pow((l1z3 - hitz3), 2))); + } else { + expd2HitL1Tk3 = 0.; + } + } +} + +double SeedMvaEstimatorPhase2::computeMva(const TrajectorySeed& seed, + const GlobalVector& global_p, + const GlobalPoint& global_x, + const edm::Handle& h_L1TkMu, + const edm::ESHandle& magfieldH, + const Propagator& propagatorAlong, + const GeometricSearchTracker& geomTracker) const { + float Phase2var[Phase2::kLast]{}; + + Phase2var[Phase2::kTsosErr0] = seed.startingState().error(0); + Phase2var[Phase2::kTsosErr2] = seed.startingState().error(2); + Phase2var[Phase2::kTsosErr5] = seed.startingState().error(5); + Phase2var[Phase2::kTsosDxdz] = seed.startingState().parameters().dxdz(); + Phase2var[Phase2::kTsosDydz] = seed.startingState().parameters().dydz(); + Phase2var[Phase2::kTsosQbp] = seed.startingState().parameters().qbp(); + + float initDRdPhi = 99999.; + + float DRL1TkMuSeedP = initDRdPhi; + float DPhiL1TkMuSeedP = initDRdPhi; + getL1TTVariables(seed, global_p, global_x, h_L1TkMu, DRL1TkMuSeedP, DPhiL1TkMuSeedP); + Phase2var[Phase2::kDRdRL1TkMuSeedP] = DRL1TkMuSeedP; + Phase2var[Phase2::kDRdPhiL1TkMuSeedP] = DPhiL1TkMuSeedP; + + float expd2HitL1Tk1 = initDRdPhi; + float expd2HitL1Tk2 = initDRdPhi; + float expd2HitL1Tk3 = initDRdPhi; + getHitL1TkVariables( + seed, h_L1TkMu, magfieldH, propagatorAlong, geomTracker, expd2HitL1Tk1, expd2HitL1Tk2, expd2HitL1Tk3); + Phase2var[Phase2::kExpd2HitL1Tk1] = expd2HitL1Tk1; + Phase2var[Phase2::kExpd2HitL1Tk2] = expd2HitL1Tk2; + Phase2var[Phase2::kExpd2HitL1Tk3] = expd2HitL1Tk3; + + for (int iv = 0; iv < Phase2::kLast; ++iv) { + Phase2var[iv] = (Phase2var[iv] - scale_mean_.at(iv)) / scale_std_.at(iv); + } + + return gbrForest_->GetResponse(Phase2var); +}