Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:27:15

0001 /**
0002   \class    TSGForOIDNN
0003   \brief    Create L3MuonTrajectorySeeds from L2 Muons in an outside-in manner
0004   \author   Dmitry Kondratyev, Arnab Purohit, Jan-Frederik Schulte (Purdue University, West Lafayette, USA)
0005  */
0006 
0007 #include "DataFormats/TrackReco/interface/Track.h"
0008 #include "FWCore/Framework/interface/global/EDProducer.h"
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/Framework/interface/EventSetup.h"
0011 #include "FWCore/Framework/interface/ESHandle.h"
0012 #include "FWCore/Framework/interface/MakerMacros.h"
0013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
0016 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0017 #include "MagneticField/Engine/interface/MagneticField.h"
0018 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0019 #include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h"
0020 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0021 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
0022 #include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
0023 #include "TrackingTools/MeasurementDet/interface/MeasurementDet.h"
0024 #include "TrackingTools/PatternTools/interface/TrajMeasLessEstim.h"
0025 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0026 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0027 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0028 #include "TrackingTools/GeomPropagators/interface/StateOnTrackerBound.h"
0029 #include "PhysicsTools/TensorFlow/interface/TensorFlow.h"
0030 #include "TrackingTools/DetLayers/interface/NavigationSchool.h"
0031 #include "RecoTracker/Record/interface/NavigationSchoolRecord.h"
0032 #include "DataFormats/Math/interface/deltaR.h"
0033 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0034 #include <boost/property_tree/ptree.hpp>
0035 #include <boost/property_tree/json_parser.hpp>
0036 #include <boost/range/adaptor/reversed.hpp>
0037 #include <memory>
0038 namespace pt = boost::property_tree;
0039 
0040 class TSGForOIDNN : public edm::global::EDProducer<> {
0041 public:
0042   explicit TSGForOIDNN(const edm::ParameterSet& iConfig);
0043   ~TSGForOIDNN() override;
0044   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0045   void produce(edm::StreamID sid, edm::Event& iEvent, const edm::EventSetup& iSetup) const override;
0046   void beginJob() override;
0047 
0048 private:
0049   /// Labels for input collections
0050   const edm::EDGetTokenT<reco::TrackCollection> src_;
0051   /// Tokens for ESHandle
0052   const edm::ESGetToken<Chi2MeasurementEstimatorBase, TrackingComponentsRecord> t_estimatorH_;
0053   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> t_magfieldH_;
0054   const edm::ESGetToken<Propagator, TrackingComponentsRecord> t_propagatorAlongH_;
0055   const edm::ESGetToken<Propagator, TrackingComponentsRecord> t_propagatorOppositeH_;
0056   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> t_tmpTkGeometryH_;
0057   const edm::ESGetToken<GlobalTrackingGeometry, GlobalTrackingGeometryRecord> t_geometryH_;
0058   const edm::ESGetToken<NavigationSchool, NavigationSchoolRecord> t_navSchool_;
0059   const edm::ESGetToken<Propagator, TrackingComponentsRecord> t_SHPOpposite_;
0060   const std::unique_ptr<TrajectoryStateUpdator> updator_;
0061   const edm::EDGetTokenT<MeasurementTrackerEvent> measurementTrackerTag_;
0062   const std::string theCategory_;
0063 
0064   /// Maximum number of seeds for each L2
0065   const unsigned int maxSeeds_;
0066   /// Maximum number of hitbased seeds for each L2
0067   const unsigned int maxHitSeeds_;
0068   /// Maximum number of hitless seeds for each L2
0069   const unsigned int maxHitlessSeeds_;
0070   /// How many layers to try
0071   const unsigned int numOfLayersToTry_;
0072   /// How many hits to try per layer
0073   const unsigned int numOfHitsToTry_;
0074   /// Rescale L2 parameter uncertainties (fixed error vs pT, eta)
0075   const double fixedErrorRescalingForHitless_;
0076   /// Estimator used to find dets and TrajectoryMeasurements
0077   const std::string estimatorName_;
0078   /// Minimum eta value to activate searching in the TEC
0079   const double minEtaForTEC_;
0080   /// Maximum eta value to activate searching in the TOB
0081   const double maxEtaForTOB_;
0082 
0083   /// IP refers to TSOS at interaction point,
0084   /// MuS refers to TSOS at muon system
0085   const unsigned int maxHitlessSeedsIP_;
0086   const unsigned int maxHitlessSeedsMuS_;
0087   const unsigned int maxHitDoubletSeeds_;
0088 
0089   /// Get number of seeds to use from DNN output instead of "max..Seeds" parameters
0090   const bool getStrategyFromDNN_;
0091   /// Whether to use DNN regressor (if false, will use classifier)
0092   const bool useRegressor_;
0093 
0094   /// Settings for classifier
0095   std::string dnnModelPath_;
0096   std::unique_ptr<tensorflow::GraphDef> graphDef_;
0097   tensorflow::Session* tf_session_;
0098 
0099   /// Settings for regressor
0100   std::string dnnModelPath_HB_;
0101   std::string dnnModelPath_HLIP_;
0102   std::string dnnModelPath_HLMuS_;
0103   std::unique_ptr<tensorflow::GraphDef> graphDef_HB_;
0104   tensorflow::Session* tf_session_HB_;
0105   std::unique_ptr<tensorflow::GraphDef> graphDef_HLIP_;
0106   tensorflow::Session* tf_session_HLIP_;
0107   std::unique_ptr<tensorflow::GraphDef> graphDef_HLMuS_;
0108   tensorflow::Session* tf_session_HLMuS_;
0109 
0110   /// DNN metadata
0111   const std::string dnnMetadataPath_;
0112   pt::ptree metadata_;
0113 
0114   /// Create seeds without hits on a given layer (TOB or TEC)
0115   void makeSeedsWithoutHits(const GeometricSearchDet& layer,
0116                             const TrajectoryStateOnSurface& tsos,
0117                             const Propagator& propagatorAlong,
0118                             const Chi2MeasurementEstimatorBase& estimator,
0119                             double errorSF,
0120                             unsigned int& hitlessSeedsMade,
0121                             unsigned int& numSeedsMade,
0122                             std::vector<TrajectorySeed>& out) const;
0123 
0124   /// Find hits on a given layer (TOB or TEC) and create seeds from updated TSOS with hit
0125   void makeSeedsFromHits(const GeometricSearchDet& layer,
0126                          const TrajectoryStateOnSurface& tsos,
0127                          const Propagator& propagatorAlong,
0128                          const Chi2MeasurementEstimatorBase& estimator,
0129                          const MeasurementTrackerEvent& measurementTracker,
0130                          unsigned int& hitSeedsMade,
0131                          unsigned int& numSeedsMade,
0132                          const unsigned int& maxHitSeeds,
0133                          unsigned int& layerCount,
0134                          std::vector<TrajectorySeed>& out) const;
0135 
0136   /// Similar to makeSeedsFromHits, but seed is created only if there are compatible hits on two adjacent layers
0137   void makeSeedsFromHitDoublets(const GeometricSearchDet& layer,
0138                                 const TrajectoryStateOnSurface& tsos,
0139                                 const Propagator& propagatorAlong,
0140                                 const Chi2MeasurementEstimatorBase& estimator,
0141                                 const MeasurementTrackerEvent& measurementTracker,
0142                                 const NavigationSchool& navSchool,
0143                                 unsigned int& hitDoubletSeedsMade,
0144                                 unsigned int& numSeedsMade,
0145                                 const unsigned int& maxHitDoubletSeeds,
0146                                 unsigned int& layerCount,
0147                                 std::vector<TrajectorySeed>& out) const;
0148 
0149   /// Update dictionary of inputs for DNN
0150   void updateFeatureMap(std::unordered_map<std::string, float>& the_map,
0151                         const reco::Track& l2,
0152                         const TrajectoryStateOnSurface& tsos_IP,
0153                         const TrajectoryStateOnSurface& tsos_MuS) const;
0154 
0155   void initializeTensorflow() const;
0156 
0157   /// Container for DNN outupts
0158   struct StrategyParameters {
0159     int nHBd, nHLIP, nHLMuS, sf;
0160   };
0161 
0162   /// Evaluate DNN classifier
0163   void evaluateClassifier(const std::unordered_map<std::string, float>& feature_map,
0164                           tensorflow::Session* session,
0165                           const pt::ptree& metadata,
0166                           StrategyParameters& out,
0167                           bool& dnnSuccess) const;
0168 
0169   /// Evaluate DNN regressor
0170   void evaluateRegressor(const std::unordered_map<std::string, float>& feature_map,
0171                          tensorflow::Session* session_HB,
0172                          const pt::ptree& metadata_HB,
0173                          tensorflow::Session* session_HLIP,
0174                          const pt::ptree& metadata_HLIP,
0175                          tensorflow::Session* session_HLMuS,
0176                          const pt::ptree& metadata_HLMuS,
0177                          StrategyParameters& out,
0178                          bool& dnnSuccess) const;
0179 };
0180 
0181 TSGForOIDNN::TSGForOIDNN(const edm::ParameterSet& iConfig)
0182     : src_(consumes(iConfig.getParameter<edm::InputTag>("src"))),
0183       t_estimatorH_(esConsumes(edm::ESInputTag("", iConfig.getParameter<std::string>("estimator")))),
0184       t_magfieldH_(esConsumes()),
0185       t_propagatorAlongH_(esConsumes(edm::ESInputTag("", iConfig.getParameter<std::string>("propagatorName")))),
0186       t_propagatorOppositeH_(esConsumes(edm::ESInputTag("", iConfig.getParameter<std::string>("propagatorName")))),
0187       t_tmpTkGeometryH_(esConsumes()),
0188       t_geometryH_(esConsumes()),
0189       t_navSchool_(esConsumes(edm::ESInputTag("", "SimpleNavigationSchool"))),
0190       t_SHPOpposite_(esConsumes(edm::ESInputTag("", "hltESPSteppingHelixPropagatorOpposite"))),
0191       updator_(new KFUpdator()),
0192       measurementTrackerTag_(consumes(iConfig.getParameter<edm::InputTag>("MeasurementTrackerEvent"))),
0193       theCategory_(std::string("Muon|RecoMuon|TSGForOIDNN")),
0194       maxSeeds_(iConfig.getParameter<uint32_t>("maxSeeds")),
0195       maxHitSeeds_(iConfig.getParameter<uint32_t>("maxHitSeeds")),
0196       maxHitlessSeeds_(iConfig.getParameter<uint32_t>("maxHitlessSeeds")),
0197       numOfLayersToTry_(iConfig.getParameter<int32_t>("layersToTry")),
0198       numOfHitsToTry_(iConfig.getParameter<int32_t>("hitsToTry")),
0199       fixedErrorRescalingForHitless_(iConfig.getParameter<double>("fixedErrorRescaleFactorForHitless")),
0200       minEtaForTEC_(iConfig.getParameter<double>("minEtaForTEC")),
0201       maxEtaForTOB_(iConfig.getParameter<double>("maxEtaForTOB")),
0202       maxHitlessSeedsIP_(iConfig.getParameter<uint32_t>("maxHitlessSeedsIP")),
0203       maxHitlessSeedsMuS_(iConfig.getParameter<uint32_t>("maxHitlessSeedsMuS")),
0204       maxHitDoubletSeeds_(iConfig.getParameter<uint32_t>("maxHitDoubletSeeds")),
0205       getStrategyFromDNN_(iConfig.getParameter<bool>("getStrategyFromDNN")),
0206       useRegressor_(iConfig.getParameter<bool>("useRegressor")),
0207       dnnMetadataPath_(iConfig.getParameter<std::string>("dnnMetadataPath")) {
0208   if (getStrategyFromDNN_) {
0209     edm::FileInPath dnnMetadataPath(dnnMetadataPath_);
0210     pt::read_json(dnnMetadataPath.fullPath(), metadata_);
0211     tensorflow::setLogging("3");
0212 
0213     if (useRegressor_) {
0214       // use regressor
0215       dnnModelPath_HB_ = metadata_.get<std::string>("HB.dnnmodel_path");
0216       edm::FileInPath dnnPath_HB(dnnModelPath_HB_);
0217       graphDef_HB_ = std::unique_ptr<tensorflow::GraphDef>(tensorflow::loadGraphDef(dnnPath_HB.fullPath()));
0218       tf_session_HB_ = tensorflow::createSession(graphDef_HB_.get());
0219 
0220       dnnModelPath_HLIP_ = metadata_.get<std::string>("HLIP.dnnmodel_path");
0221       edm::FileInPath dnnPath_HLIP(dnnModelPath_HLIP_);
0222       graphDef_HLIP_ = std::unique_ptr<tensorflow::GraphDef>(tensorflow::loadGraphDef(dnnPath_HLIP.fullPath()));
0223       tf_session_HLIP_ = tensorflow::createSession(graphDef_HLIP_.get());
0224 
0225       dnnModelPath_HLMuS_ = metadata_.get<std::string>("HLMuS.dnnmodel_path");
0226       edm::FileInPath dnnPath_HLMuS(dnnModelPath_HLMuS_);
0227       graphDef_HLMuS_ = std::unique_ptr<tensorflow::GraphDef>(tensorflow::loadGraphDef(dnnPath_HLMuS.fullPath()));
0228       tf_session_HLMuS_ = tensorflow::createSession(graphDef_HLMuS_.get());
0229     } else {
0230       // use classifier (default)
0231       dnnModelPath_ = metadata_.get<std::string>("dnnmodel_path");
0232       edm::FileInPath dnnPath(dnnModelPath_);
0233       graphDef_ = std::unique_ptr<tensorflow::GraphDef>(tensorflow::loadGraphDef(dnnPath.fullPath()));
0234       tf_session_ = tensorflow::createSession(graphDef_.get());
0235     }
0236   }
0237   produces<std::vector<TrajectorySeed> >();
0238 }
0239 
0240 TSGForOIDNN::~TSGForOIDNN() {
0241   if (getStrategyFromDNN_) {
0242     if (useRegressor_) {
0243       tensorflow::closeSession(tf_session_HB_);
0244       tensorflow::closeSession(tf_session_HLIP_);
0245       tensorflow::closeSession(tf_session_HLMuS_);
0246     } else {
0247       tensorflow::closeSession(tf_session_);
0248     }
0249   }
0250 }
0251 
0252 void TSGForOIDNN::beginJob() { initializeTensorflow(); }
0253 
0254 namespace {
0255   std::unordered_map<std::string, float> dummyFeatureMap() {
0256     std::unordered_map<std::string, float> the_map;
0257     the_map["pt"] = 2.0;
0258     the_map["eta"] = 0.5;
0259     the_map["phi"] = 1.0;
0260     the_map["validHits"] = 0.;
0261     the_map["tsos_IP_eta"] = -999;
0262     the_map["tsos_IP_phi"] = -999;
0263     the_map["tsos_IP_pt"] = -999;
0264     the_map["tsos_IP_pt_eta"] = -999;
0265     the_map["tsos_IP_pt_phi"] = -999;
0266     the_map["err0_IP"] = -999;
0267     the_map["err1_IP"] = -999;
0268     the_map["err2_IP"] = -999;
0269     the_map["err3_IP"] = -999;
0270     the_map["err4_IP"] = -999;
0271     the_map["tsos_IP_valid"] = 0.0;
0272     the_map["tsos_MuS_eta"] = -999;
0273     the_map["tsos_MuS_phi"] = -999;
0274     the_map["tsos_MuS_pt"] = -999;
0275     the_map["tsos_MuS_pt_eta"] = -999;
0276     the_map["tsos_MuS_pt_phi"] = -999;
0277     the_map["err0_MuS"] = -999;
0278     the_map["err1_MuS"] = -999;
0279     the_map["err2_MuS"] = -999;
0280     the_map["err3_MuS"] = -999;
0281     the_map["err4_MuS"] = -999;
0282     the_map["tsos_MuS_valid"] = 0.0;
0283     return the_map;
0284   }
0285 }  // namespace
0286 
0287 void TSGForOIDNN::initializeTensorflow() const {
0288   if (getStrategyFromDNN_ and not useRegressor_) {
0289     // Container for DNN outputs
0290     StrategyParameters strPars;
0291     bool dnnSuccess = false;
0292     evaluateClassifier(dummyFeatureMap(), tf_session_, metadata_, strPars, dnnSuccess);
0293   }
0294 }
0295 
0296 //
0297 // Produce seeds
0298 //
0299 void TSGForOIDNN::produce(edm::StreamID sid, edm::Event& iEvent, edm::EventSetup const& iEventSetup) const {
0300   // Initialize variables
0301   unsigned int numSeedsMade = 0;
0302   unsigned int layerCount = 0;
0303   unsigned int hitlessSeedsMadeIP = 0;
0304   unsigned int hitlessSeedsMadeMuS = 0;
0305   unsigned int hitSeedsMade = 0;
0306   unsigned int hitDoubletSeedsMade = 0;
0307 
0308   // Container for DNN inputs
0309   std::unordered_map<std::string, float> feature_map;
0310 
0311   // Container for DNN outputs
0312   StrategyParameters strPars;
0313 
0314   // Surface used to make a TSOS at the PCA to the beamline
0315   Plane::PlanePointer dummyPlane = Plane::build(Plane::PositionType(), Plane::RotationType());
0316 
0317   // Get setup objects
0318   const MagneticField& magfield = iEventSetup.getData(t_magfieldH_);
0319   const Chi2MeasurementEstimatorBase& estimator = iEventSetup.getData(t_estimatorH_);
0320   const Propagator& tmpPropagatorAlong = iEventSetup.getData(t_propagatorAlongH_);
0321   const Propagator& tmpPropagatorOpposite = iEventSetup.getData(t_propagatorOppositeH_);
0322   const TrackerGeometry& tmpTkGeometry = iEventSetup.getData(t_tmpTkGeometryH_);
0323   const GlobalTrackingGeometry& geometry = iEventSetup.getData(t_geometryH_);
0324   const NavigationSchool& navSchool = iEventSetup.getData(t_navSchool_);
0325   auto const& measurementTracker = iEvent.get(measurementTrackerTag_);
0326 
0327   // Read L2 track collection
0328   auto const& l2TrackCol = iEvent.get(src_);
0329 
0330   // The product
0331   std::unique_ptr<std::vector<TrajectorySeed> > result(new std::vector<TrajectorySeed>());
0332 
0333   // Get vector of Detector layers
0334   auto const* gsTracker = measurementTracker.geometricSearchTracker();
0335   std::vector<BarrelDetLayer const*> const& tob = gsTracker->tobLayers();
0336   std::vector<ForwardDetLayer const*> const& tecPositive =
0337       tmpTkGeometry.isThere(GeomDetEnumerators::P2OTEC) ? gsTracker->posTidLayers() : gsTracker->posTecLayers();
0338   std::vector<ForwardDetLayer const*> const& tecNegative =
0339       tmpTkGeometry.isThere(GeomDetEnumerators::P2OTEC) ? gsTracker->negTidLayers() : gsTracker->negTecLayers();
0340 
0341   // Get suitable propagators
0342   std::unique_ptr<Propagator> propagatorAlong = SetPropagationDirection(tmpPropagatorAlong, alongMomentum);
0343   std::unique_ptr<Propagator> propagatorOpposite = SetPropagationDirection(tmpPropagatorOpposite, oppositeToMomentum);
0344 
0345   // Stepping Helix Propagator for propogation from muon system to tracker
0346   edm::ESHandle<Propagator> shpOpposite = iEventSetup.getHandle(t_SHPOpposite_);
0347 
0348   // Loop over the L2's and make seeds for all of them
0349   LogTrace(theCategory_) << "TSGForOIDNN::produce: Number of L2's: " << l2TrackCol.size();
0350 
0351   for (auto const& l2 : l2TrackCol) {
0352     // Container of Seeds
0353     std::vector<TrajectorySeed> out;
0354     LogTrace("TSGForOIDNN") << "TSGForOIDNN::produce: L2 muon pT, eta, phi --> " << l2.pt() << " , " << l2.eta()
0355                             << " , " << l2.phi();
0356 
0357     FreeTrajectoryState fts = trajectoryStateTransform::initialFreeState(l2, &magfield);
0358 
0359     dummyPlane->move(fts.position() - dummyPlane->position());
0360     TrajectoryStateOnSurface tsosAtIP = TrajectoryStateOnSurface(fts, *dummyPlane);
0361     LogTrace("TSGForOIDNN") << "TSGForOIDNN::produce: Created TSOSatIP: " << tsosAtIP;
0362 
0363     // Get the TSOS on the innermost layer of the L2
0364     TrajectoryStateOnSurface tsosAtMuonSystem = trajectoryStateTransform::innerStateOnSurface(l2, geometry, &magfield);
0365     LogTrace("TSGForOIDNN") << "TSGForOIDNN::produce: Created TSOSatMuonSystem: " << tsosAtMuonSystem;
0366 
0367     LogTrace("TSGForOIDNN")
0368         << "TSGForOIDNN::produce: Check the error of the L2 parameter and use hit seeds if big errors";
0369 
0370     StateOnTrackerBound fromInside(propagatorAlong.get());
0371     TrajectoryStateOnSurface outerTkStateInside = fromInside(fts);
0372 
0373     StateOnTrackerBound fromOutside(&*shpOpposite);
0374     TrajectoryStateOnSurface outerTkStateOutside = fromOutside(tsosAtMuonSystem);
0375 
0376     // Check if the two positions (using updated and not-updated TSOS) agree withing certain extent.
0377     // If both TSOSs agree, use only the one at vertex, as it uses more information. If they do not agree, search for seeds based on both.
0378     double l2muonEta = l2.eta();
0379     double absL2muonEta = std::abs(l2muonEta);
0380 
0381     // make non-const copies of parameters, so they can be overriden for individual L2 muons
0382     unsigned int maxHitSeeds = maxHitSeeds_;
0383     unsigned int maxHitDoubletSeeds = maxHitDoubletSeeds_;
0384     unsigned int maxHitlessSeedsIP = maxHitlessSeedsIP_;
0385     unsigned int maxHitlessSeedsMuS = maxHitlessSeedsMuS_;
0386 
0387     float errorSFHitless = fixedErrorRescalingForHitless_;
0388 
0389     // update strategy parameters by evaluating DNN
0390     if (getStrategyFromDNN_) {
0391       bool dnnSuccess = false;
0392 
0393       // Update feature map with parameters of the current muon
0394       updateFeatureMap(feature_map, l2, tsosAtIP, outerTkStateOutside);
0395 
0396       if (useRegressor_) {
0397         // Use regressor
0398         evaluateRegressor(feature_map,
0399                           tf_session_HB_,
0400                           metadata_.get_child("HB"),
0401                           tf_session_HLIP_,
0402                           metadata_.get_child("HLIP"),
0403                           tf_session_HLMuS_,
0404                           metadata_.get_child("HLMuS"),
0405                           strPars,
0406                           dnnSuccess);
0407       } else {
0408         // Use classifier
0409         evaluateClassifier(feature_map, tf_session_, metadata_, strPars, dnnSuccess);
0410       }
0411       if (!dnnSuccess)
0412         break;
0413 
0414       maxHitSeeds = 0;
0415       maxHitDoubletSeeds = strPars.nHBd;
0416       maxHitlessSeedsIP = strPars.nHLIP;
0417       maxHitlessSeedsMuS = strPars.nHLMuS;
0418       errorSFHitless = strPars.sf;
0419     }
0420 
0421     numSeedsMade = 0;
0422     hitlessSeedsMadeIP = 0;
0423     hitlessSeedsMadeMuS = 0;
0424     hitSeedsMade = 0;
0425     hitDoubletSeedsMade = 0;
0426 
0427     auto createSeeds = [&](auto const& layers) {
0428       for (auto const& layer : boost::adaptors::reverse(layers)) {
0429         if (hitlessSeedsMadeIP < maxHitlessSeedsIP && numSeedsMade < maxSeeds_)
0430           makeSeedsWithoutHits(*layer,
0431                                tsosAtIP,
0432                                *(propagatorAlong.get()),
0433                                estimator,
0434                                errorSFHitless,
0435                                hitlessSeedsMadeIP,
0436                                numSeedsMade,
0437                                out);
0438 
0439         if (outerTkStateInside.isValid() && outerTkStateOutside.isValid() && hitlessSeedsMadeMuS < maxHitlessSeedsMuS &&
0440             numSeedsMade < maxSeeds_)
0441           makeSeedsWithoutHits(*layer,
0442                                outerTkStateOutside,
0443                                *(propagatorOpposite.get()),
0444                                estimator,
0445                                errorSFHitless,
0446                                hitlessSeedsMadeMuS,
0447                                numSeedsMade,
0448                                out);
0449 
0450         if (hitSeedsMade < maxHitSeeds && numSeedsMade < maxSeeds_)
0451           makeSeedsFromHits(*layer,
0452                             tsosAtIP,
0453                             *(propagatorAlong.get()),
0454                             estimator,
0455                             measurementTracker,
0456                             hitSeedsMade,
0457                             numSeedsMade,
0458                             maxHitSeeds,
0459                             layerCount,
0460                             out);
0461 
0462         if (hitDoubletSeedsMade < maxHitDoubletSeeds && numSeedsMade < maxSeeds_)
0463           makeSeedsFromHitDoublets(*layer,
0464                                    tsosAtIP,
0465                                    *(propagatorAlong.get()),
0466                                    estimator,
0467                                    measurementTracker,
0468                                    navSchool,
0469                                    hitDoubletSeedsMade,
0470                                    numSeedsMade,
0471                                    maxHitDoubletSeeds,
0472                                    layerCount,
0473                                    out);
0474       };
0475     };
0476 
0477     // BARREL
0478     if (absL2muonEta < maxEtaForTOB_) {
0479       layerCount = 0;
0480       createSeeds(tob);
0481       LogTrace("TSGForOIDNN") << "TSGForOIDNN:::produce: NumSeedsMade = " << numSeedsMade
0482                               << " , layerCount = " << layerCount;
0483     }
0484 
0485     // Reset number of seeds if in overlap region
0486     if (absL2muonEta > minEtaForTEC_ && absL2muonEta < maxEtaForTOB_) {
0487       numSeedsMade = 0;
0488       hitlessSeedsMadeIP = 0;
0489       hitlessSeedsMadeMuS = 0;
0490       hitSeedsMade = 0;
0491       hitDoubletSeedsMade = 0;
0492     }
0493 
0494     // ENDCAP+
0495     if (l2muonEta > minEtaForTEC_) {
0496       layerCount = 0;
0497       createSeeds(tecPositive);
0498       LogTrace("TSGForOIDNN") << "TSGForOIDNN:::produce: NumSeedsMade = " << numSeedsMade
0499                               << " , layerCount = " << layerCount;
0500     }
0501 
0502     // ENDCAP-
0503     if (l2muonEta < -minEtaForTEC_) {
0504       layerCount = 0;
0505       createSeeds(tecNegative);
0506       LogTrace("TSGForOIDNN") << "TSGForOIDNN:::produce: NumSeedsMade = " << numSeedsMade
0507                               << " , layerCount = " << layerCount;
0508     }
0509 
0510     for (std::vector<TrajectorySeed>::iterator it = out.begin(); it != out.end(); ++it) {
0511       result->push_back(*it);
0512     }
0513 
0514   }  // L2Collection
0515 
0516   edm::LogInfo(theCategory_) << "TSGForOIDNN::produce: number of seeds made: " << result->size();
0517 
0518   iEvent.put(std::move(result));
0519 }
0520 
0521 //
0522 // Create seeds without hits on a given layer (TOB or TEC)
0523 //
0524 void TSGForOIDNN::makeSeedsWithoutHits(const GeometricSearchDet& layer,
0525                                        const TrajectoryStateOnSurface& tsos,
0526                                        const Propagator& propagatorAlong,
0527                                        const Chi2MeasurementEstimatorBase& estimator,
0528                                        double errorSF,
0529                                        unsigned int& hitlessSeedsMade,
0530                                        unsigned int& numSeedsMade,
0531                                        std::vector<TrajectorySeed>& out) const {
0532   // create hitless seeds
0533   LogTrace("TSGForOIDNN") << "TSGForOIDNN::makeSeedsWithoutHits: Start hitless";
0534   std::vector<GeometricSearchDet::DetWithState> dets;
0535   layer.compatibleDetsV(tsos, propagatorAlong, estimator, dets);
0536   if (!dets.empty()) {
0537     auto const& detOnLayer = dets.front().first;
0538     auto& tsosOnLayer = dets.front().second;
0539     LogTrace("TSGForOIDNN") << "TSGForOIDNN::makeSeedsWithoutHits: tsosOnLayer " << tsosOnLayer;
0540     if (!tsosOnLayer.isValid()) {
0541       edm::LogInfo(theCategory_) << "ERROR!: Hitless TSOS is not valid!";
0542     } else {
0543       tsosOnLayer.rescaleError(errorSF);
0544       PTrajectoryStateOnDet const& ptsod =
0545           trajectoryStateTransform::persistentState(tsosOnLayer, detOnLayer->geographicalId().rawId());
0546       TrajectorySeed::RecHitContainer rHC;
0547       out.push_back(TrajectorySeed(ptsod, rHC, oppositeToMomentum));
0548       LogTrace("TSGForOIDNN") << "TSGForOIDNN::makeSeedsWithoutHits: TSOS (Hitless) done ";
0549       hitlessSeedsMade++;
0550       numSeedsMade++;
0551     }
0552   }
0553 }
0554 
0555 //
0556 // Find hits on a given layer (TOB or TEC) and create seeds from updated TSOS with hit
0557 //
0558 void TSGForOIDNN::makeSeedsFromHits(const GeometricSearchDet& layer,
0559                                     const TrajectoryStateOnSurface& tsos,
0560                                     const Propagator& propagatorAlong,
0561                                     const Chi2MeasurementEstimatorBase& estimator,
0562                                     const MeasurementTrackerEvent& measurementTracker,
0563                                     unsigned int& hitSeedsMade,
0564                                     unsigned int& numSeedsMade,
0565                                     const unsigned int& maxHitSeeds,
0566                                     unsigned int& layerCount,
0567                                     std::vector<TrajectorySeed>& out) const {
0568   if (layerCount > numOfLayersToTry_)
0569     return;
0570 
0571   const TrajectoryStateOnSurface& onLayer(tsos);
0572 
0573   std::vector<GeometricSearchDet::DetWithState> dets;
0574   layer.compatibleDetsV(onLayer, propagatorAlong, estimator, dets);
0575 
0576   // Find Measurements on each DetWithState
0577   LogTrace("TSGForOIDNN") << "TSGForOIDNN::makeSeedsFromHits: Find measurements on each detWithState  " << dets.size();
0578   std::vector<TrajectoryMeasurement> meas;
0579   for (auto const& detI : dets) {
0580     MeasurementDetWithData det = measurementTracker.idToDet(detI.first->geographicalId());
0581     if (det.isNull())
0582       continue;
0583     if (!detI.second.isValid())
0584       continue;  // Skip if TSOS is not valid
0585 
0586     std::vector<TrajectoryMeasurement> mymeas =
0587         det.fastMeasurements(detI.second, onLayer, propagatorAlong, estimator);  // Second TSOS is not used
0588     for (auto const& measurement : mymeas) {
0589       if (measurement.recHit()->isValid())
0590         meas.push_back(measurement);  // Only save those which are valid
0591     }
0592   }
0593 
0594   // Update TSOS using TMs after sorting, then create Trajectory Seed and put into vector
0595   LogTrace("TSGForOIDNN") << "TSGForOIDNN::makeSeedsFromHits: Update TSOS using TMs after sorting, then create "
0596                              "Trajectory Seed, number of TM = "
0597                           << meas.size();
0598   std::sort(meas.begin(), meas.end(), TrajMeasLessEstim());
0599 
0600   unsigned int found = 0;
0601   for (auto const& measurement : meas) {
0602     if (hitSeedsMade >= maxHitSeeds)
0603       return;
0604     TrajectoryStateOnSurface updatedTSOS = updator_->update(measurement.forwardPredictedState(), *measurement.recHit());
0605     LogTrace("TSGForOIDNN") << "TSGForOIDNN::makeSeedsFromHits: TSOS for TM " << found;
0606     if (not updatedTSOS.isValid())
0607       continue;
0608 
0609     edm::OwnVector<TrackingRecHit> seedHits;
0610     seedHits.push_back(*measurement.recHit()->hit());
0611     PTrajectoryStateOnDet const& pstate =
0612         trajectoryStateTransform::persistentState(updatedTSOS, measurement.recHit()->geographicalId().rawId());
0613     LogTrace("TSGForOIDNN") << "TSGForOIDNN::makeSeedsFromHits: Number of seedHits: " << seedHits.size();
0614     TrajectorySeed seed(pstate, std::move(seedHits), oppositeToMomentum);
0615     out.push_back(seed);
0616     found++;
0617     numSeedsMade++;
0618     hitSeedsMade++;
0619     if (found == numOfHitsToTry_)
0620       break;
0621   }
0622 
0623   if (found)
0624     layerCount++;
0625 }
0626 
0627 //
0628 // Find hits compatible with L2 trajectory on two adjacent layers; if found, create a seed using both hits
0629 //
0630 void TSGForOIDNN::makeSeedsFromHitDoublets(const GeometricSearchDet& layer,
0631                                            const TrajectoryStateOnSurface& tsos,
0632                                            const Propagator& propagatorAlong,
0633                                            const Chi2MeasurementEstimatorBase& estimator,
0634                                            const MeasurementTrackerEvent& measurementTracker,
0635                                            const NavigationSchool& navSchool,
0636                                            unsigned int& hitDoubletSeedsMade,
0637                                            unsigned int& numSeedsMade,
0638                                            const unsigned int& maxHitDoubletSeeds,
0639                                            unsigned int& layerCount,
0640                                            std::vector<TrajectorySeed>& out) const {
0641   // This method is similar to makeSeedsFromHits, but the seed is created
0642   // only when in addition to a hit on a given layer, there are more compatible hits
0643   // on next layers (going from outside inwards), compatible with updated TSOS.
0644   // If that's the case, multiple compatible hits are used to create a single seed.
0645 
0646   // Configured to only check the immideately adjacent layer and add one more hit
0647   int max_addtnl_layers = 1;  // max number of additional layers to scan
0648   int max_meas = 1;           // number of measurements to consider on each additional layer
0649 
0650   // // // First, regular procedure to find a compatible hit - like in makeSeedsFromHits // // //
0651 
0652   const TrajectoryStateOnSurface& onLayer(tsos);
0653 
0654   // Find dets compatible with original TSOS
0655   std::vector<GeometricSearchDet::DetWithState> dets;
0656   layer.compatibleDetsV(onLayer, propagatorAlong, estimator, dets);
0657 
0658   LogTrace("TSGForOIDNN") << "TSGForOIDNN::makeSeedsFromHitDoublets: Find measurements on each detWithState  "
0659                           << dets.size();
0660   std::vector<TrajectoryMeasurement> meas;
0661 
0662   // Loop over dets
0663   for (auto const& detI : dets) {
0664     MeasurementDetWithData det = measurementTracker.idToDet(detI.first->geographicalId());
0665 
0666     if (det.isNull())
0667       continue;  // skip if det does not exist
0668     if (!detI.second.isValid())
0669       continue;  // skip if TSOS is invalid
0670 
0671     // Find measurements on this det
0672     std::vector<TrajectoryMeasurement> mymeas = det.fastMeasurements(detI.second, onLayer, propagatorAlong, estimator);
0673 
0674     // Save valid measurements
0675     for (auto const& measurement : mymeas) {
0676       if (measurement.recHit()->isValid())
0677         meas.push_back(measurement);
0678     }  // end loop over meas
0679   }    // end loop over dets
0680 
0681   LogTrace("TSGForOIDNN") << "TSGForOIDNN::makeSeedsFromHitDoublets: Update TSOS using TMs after sorting, then create "
0682                              "Trajectory Seed, number of TM = "
0683                           << meas.size();
0684 
0685   // sort valid measurements found on the first layer
0686   std::sort(meas.begin(), meas.end(), TrajMeasLessEstim());
0687 
0688   unsigned int found = 0;
0689 
0690   // Loop over all valid measurements compatible with original TSOS
0691   //for (std::vector<TrajectoryMeasurement>::const_iterator mea = meas.begin(); mea != meas.end(); ++mea) {
0692   for (auto const& measurement : meas) {
0693     if (hitDoubletSeedsMade >= maxHitDoubletSeeds)
0694       return;  // abort if enough seeds created
0695 
0696     // Update TSOS with measurement on first considered layer
0697     TrajectoryStateOnSurface updatedTSOS = updator_->update(measurement.forwardPredictedState(), *measurement.recHit());
0698 
0699     LogTrace("TSGForOIDNN") << "TSGForOIDNN::makeSeedsFromHitDoublets: TSOS for TM " << found;
0700     if (not updatedTSOS.isValid())
0701       continue;  // Skip if updated TSOS is invalid
0702 
0703     edm::OwnVector<TrackingRecHit> seedHits;
0704 
0705     // Save hit on first layer
0706     seedHits.push_back(*measurement.recHit()->hit());
0707     const DetLayer* detLayer = dynamic_cast<const DetLayer*>(&layer);
0708 
0709     // // // Now for this measurement we will loop over additional layers and try to update the TSOS again // // //
0710 
0711     // find layers compatible with updated TSOS
0712     auto const& compLayers = navSchool.nextLayers(*detLayer, *updatedTSOS.freeState(), alongMomentum);
0713 
0714     int addtnl_layers_scanned = 0;
0715     int found_compatible_on_next_layer = 0;
0716     int det_id = 0;
0717 
0718     // Copy updated TSOS - we will update it again with a measurement from the next layer, if we find it
0719     TrajectoryStateOnSurface updatedTSOS_next(updatedTSOS);
0720 
0721     // loop over layers compatible with updated TSOS
0722     for (auto compLayer : compLayers) {
0723       int nmeas = 0;
0724 
0725       if (addtnl_layers_scanned >= max_addtnl_layers)
0726         break;  // break if we already looped over enough layers
0727       if (found_compatible_on_next_layer > 0)
0728         break;  // break if we already found additional hit
0729 
0730       // find dets compatible with updated TSOS
0731       std::vector<GeometricSearchDet::DetWithState> dets_next;
0732       TrajectoryStateOnSurface onLayer_next(updatedTSOS);
0733 
0734       compLayer->compatibleDetsV(onLayer_next, propagatorAlong, estimator, dets_next);
0735 
0736       //if (!detWithState.size()) continue;
0737       std::vector<TrajectoryMeasurement> meas_next;
0738 
0739       // find measurements on dets_next and save the valid ones
0740       for (auto const& detI_next : dets_next) {
0741         MeasurementDetWithData det = measurementTracker.idToDet(detI_next.first->geographicalId());
0742 
0743         if (det.isNull())
0744           continue;  // skip if det does not exist
0745         if (!detI_next.second.isValid())
0746           continue;  // skip if TSOS is invalid
0747 
0748         // Find measurements on this det
0749         std::vector<TrajectoryMeasurement> mymeas_next =
0750             det.fastMeasurements(detI_next.second, onLayer_next, propagatorAlong, estimator);
0751 
0752         for (auto const& mea_next : mymeas_next) {
0753           // save valid measurements
0754           if (mea_next.recHit()->isValid())
0755             meas_next.push_back(mea_next);
0756 
0757         }  // end loop over mymeas_next
0758       }    // end loop over dets_next
0759 
0760       // sort valid measurements found on this layer
0761       std::sort(meas_next.begin(), meas_next.end(), TrajMeasLessEstim());
0762 
0763       // loop over valid measurements compatible with updated TSOS (TSOS updated with a hit on the first layer)
0764       for (auto const& mea_next : meas_next) {
0765         if (nmeas >= max_meas)
0766           break;  // skip if we already found enough hits
0767 
0768         // try to update TSOS again, with an additional hit
0769         updatedTSOS_next = updator_->update(mea_next.forwardPredictedState(), *mea_next.recHit());
0770 
0771         if (not updatedTSOS_next.isValid())
0772           continue;  // skip if TSOS updated with additional hit is not valid
0773 
0774         // If there was a compatible hit on this layer, we end up here.
0775         // An additional compatible hit is saved.
0776         seedHits.push_back(*mea_next.recHit()->hit());
0777         det_id = mea_next.recHit()->geographicalId().rawId();
0778         nmeas++;
0779         found_compatible_on_next_layer++;
0780 
0781       }  // end loop over meas_next
0782 
0783       addtnl_layers_scanned++;
0784 
0785     }  // end loop over compLayers (additional layers scanned after the original layer)
0786 
0787     if (found_compatible_on_next_layer == 0)
0788       continue;
0789     // only consider the hit if there was a compatible hit on one of the additional scanned layers
0790 
0791     // Create a seed from two saved hits
0792     PTrajectoryStateOnDet const& pstate = trajectoryStateTransform::persistentState(updatedTSOS_next, det_id);
0793     TrajectorySeed seed(pstate, std::move(seedHits), oppositeToMomentum);
0794 
0795     LogTrace("TSGForOIDNN") << "TSGForOIDNN::makeSeedsFromHitDoublets: Number of seedHits: " << seedHits.size();
0796     out.push_back(seed);
0797 
0798     found++;
0799     numSeedsMade++;
0800     hitDoubletSeedsMade++;
0801 
0802     if (found == numOfHitsToTry_)
0803       break;  // break if enough measurements scanned
0804 
0805   }  // end loop over measurements compatible with original TSOS
0806 
0807   if (found)
0808     layerCount++;
0809 }
0810 
0811 //
0812 // Update the dictionary of variables to use as input features for DNN
0813 //
0814 void TSGForOIDNN::updateFeatureMap(std::unordered_map<std::string, float>& the_map,
0815                                    const reco::Track& l2,
0816                                    const TrajectoryStateOnSurface& tsos_IP,
0817                                    const TrajectoryStateOnSurface& tsos_MuS) const {
0818   the_map["pt"] = l2.pt();
0819   the_map["eta"] = l2.eta();
0820   the_map["phi"] = l2.phi();
0821   the_map["validHits"] = l2.found();
0822   if (tsos_IP.isValid()) {
0823     the_map["tsos_IP_eta"] = tsos_IP.globalPosition().eta();
0824     the_map["tsos_IP_phi"] = tsos_IP.globalPosition().phi();
0825     the_map["tsos_IP_pt"] = tsos_IP.globalMomentum().perp();
0826     the_map["tsos_IP_pt_eta"] = tsos_IP.globalMomentum().eta();
0827     the_map["tsos_IP_pt_phi"] = tsos_IP.globalMomentum().phi();
0828     const AlgebraicSymMatrix55& matrix_IP = tsos_IP.curvilinearError().matrix();
0829     the_map["err0_IP"] = sqrt(matrix_IP[0][0]);
0830     the_map["err1_IP"] = sqrt(matrix_IP[1][1]);
0831     the_map["err2_IP"] = sqrt(matrix_IP[2][2]);
0832     the_map["err3_IP"] = sqrt(matrix_IP[3][3]);
0833     the_map["err4_IP"] = sqrt(matrix_IP[4][4]);
0834     the_map["tsos_IP_valid"] = 1.0;
0835   } else {
0836     the_map["tsos_IP_eta"] = -999;
0837     the_map["tsos_IP_phi"] = -999;
0838     the_map["tsos_IP_pt"] = -999;
0839     the_map["tsos_IP_pt_eta"] = -999;
0840     the_map["tsos_IP_pt_phi"] = -999;
0841     the_map["err0_IP"] = -999;
0842     the_map["err1_IP"] = -999;
0843     the_map["err2_IP"] = -999;
0844     the_map["err3_IP"] = -999;
0845     the_map["err4_IP"] = -999;
0846     the_map["tsos_IP_valid"] = 0.0;
0847   }
0848   if (tsos_MuS.isValid()) {
0849     the_map["tsos_MuS_eta"] = tsos_MuS.globalPosition().eta();
0850     the_map["tsos_MuS_phi"] = tsos_MuS.globalPosition().phi();
0851     the_map["tsos_MuS_pt"] = tsos_MuS.globalMomentum().perp();
0852     the_map["tsos_MuS_pt_eta"] = tsos_MuS.globalMomentum().eta();
0853     the_map["tsos_MuS_pt_phi"] = tsos_MuS.globalMomentum().phi();
0854     const AlgebraicSymMatrix55& matrix_MuS = tsos_MuS.curvilinearError().matrix();
0855     the_map["err0_MuS"] = sqrt(matrix_MuS[0][0]);
0856     the_map["err1_MuS"] = sqrt(matrix_MuS[1][1]);
0857     the_map["err2_MuS"] = sqrt(matrix_MuS[2][2]);
0858     the_map["err3_MuS"] = sqrt(matrix_MuS[3][3]);
0859     the_map["err4_MuS"] = sqrt(matrix_MuS[4][4]);
0860     the_map["tsos_MuS_valid"] = 1.0;
0861   } else {
0862     the_map["tsos_MuS_eta"] = -999;
0863     the_map["tsos_MuS_phi"] = -999;
0864     the_map["tsos_MuS_pt"] = -999;
0865     the_map["tsos_MuS_pt_eta"] = -999;
0866     the_map["tsos_MuS_pt_phi"] = -999;
0867     the_map["err0_MuS"] = -999;
0868     the_map["err1_MuS"] = -999;
0869     the_map["err2_MuS"] = -999;
0870     the_map["err3_MuS"] = -999;
0871     the_map["err4_MuS"] = -999;
0872     the_map["tsos_MuS_valid"] = 0.0;
0873   }
0874 }
0875 
0876 //
0877 // Obtain seeding strategy parameters by evaluating DNN classifier for a given muon
0878 //
0879 void TSGForOIDNN::evaluateClassifier(const std::unordered_map<std::string, float>& feature_map,
0880                                      tensorflow::Session* session,
0881                                      const pt::ptree& metadata,
0882                                      StrategyParameters& out,
0883                                      bool& dnnSuccess) const {
0884   int n_features = metadata.get<int>("n_features", 0);
0885 
0886   // Prepare tensor for DNN inputs
0887   tensorflow::Tensor input(tensorflow::DT_FLOAT, {1, n_features});
0888   std::string fname;
0889   int i_feature = 0;
0890   for (const pt::ptree::value_type& feature : metadata.get_child("feature_names")) {
0891     fname = feature.second.data();
0892     if (feature_map.find(fname) == feature_map.end()) {
0893       // don't evaluate DNN if any input feature is missing
0894       dnnSuccess = false;
0895     } else {
0896       input.matrix<float>()(0, i_feature) = float(feature_map.at(fname));
0897       i_feature++;
0898     }
0899   }
0900 
0901   // Prepare tensor for DNN outputs
0902   std::vector<tensorflow::Tensor> outputs;
0903 
0904   // Evaluate DNN and put results in output tensor
0905   std::string inputLayer = metadata.get<std::string>("input_layer");
0906   std::string outputLayer = metadata.get<std::string>("output_layer");
0907 
0908   tensorflow::run(session, {{inputLayer, input}}, {outputLayer}, &outputs);
0909   tensorflow::Tensor out_tensor = outputs[0];
0910   tensorflow::TTypes<float, 1>::Matrix dnn_outputs = out_tensor.matrix<float>();
0911 
0912   // Find output with largest prediction
0913   int imax = 0;
0914   float out_max = 0;
0915   for (long long int i = 0; i < out_tensor.dim_size(1); i++) {
0916     float ith_output = dnn_outputs(0, i);
0917     if (ith_output > out_max) {
0918       imax = i;
0919       out_max = ith_output;
0920     }
0921   }
0922 
0923   // Decode output
0924   const std::string label = "output_labels.label_" + std::to_string(imax);
0925   out.nHBd = metadata.get<int>(label + ".nHBd");
0926   out.nHLIP = metadata.get<int>(label + ".nHLIP");
0927   out.nHLMuS = metadata.get<int>(label + ".nHLMuS");
0928   out.sf = metadata.get<int>(label + ".SF");
0929 
0930   dnnSuccess = true;
0931 }
0932 
0933 //
0934 // Obtain seeding strategy parameters by evaluating DNN regressor for a given muon
0935 //
0936 void TSGForOIDNN::evaluateRegressor(const std::unordered_map<std::string, float>& feature_map,
0937                                     tensorflow::Session* session_HB,
0938                                     const pt::ptree& metadata_HB,
0939                                     tensorflow::Session* session_HLIP,
0940                                     const pt::ptree& metadata_HLIP,
0941                                     tensorflow::Session* session_HLMuS,
0942                                     const pt::ptree& metadata_HLMuS,
0943                                     StrategyParameters& out,
0944                                     bool& dnnSuccess) const {
0945   int n_features = metadata_HB.get<int>("n_features", 0);
0946 
0947   // Prepare tensor for DNN inputs
0948   tensorflow::Tensor input(tensorflow::DT_FLOAT, {1, n_features});
0949   std::string fname;
0950   int i_feature = 0;
0951   for (const pt::ptree::value_type& feature : metadata_HB.get_child("feature_names")) {
0952     fname = feature.second.data();
0953     if (feature_map.find(fname) == feature_map.end()) {
0954       // don't evaluate DNN if any input feature is missing
0955       dnnSuccess = false;
0956     } else {
0957       input.matrix<float>()(0, i_feature) = float(feature_map.at(fname));
0958       i_feature++;
0959     }
0960   }
0961 
0962   // Prepare tensor for DNN outputs
0963   std::vector<tensorflow::Tensor> outputs_HB;
0964   // Evaluate DNN and put results in output tensor
0965   std::string inputLayer_HB = metadata_HB.get<std::string>("input_layer");
0966   std::string outputLayer_HB = metadata_HB.get<std::string>("output_layer");
0967   tensorflow::run(session_HB, {{inputLayer_HB, input}}, {outputLayer_HB}, &outputs_HB);
0968   tensorflow::Tensor out_tensor_HB = outputs_HB[0];
0969   tensorflow::TTypes<float, 1>::Matrix dnn_outputs_HB = out_tensor_HB.matrix<float>();
0970 
0971   // Prepare tensor for DNN outputs
0972   std::vector<tensorflow::Tensor> outputs_HLIP;
0973   // Evaluate DNN and put results in output tensor
0974   std::string inputLayer_HLIP = metadata_HLIP.get<std::string>("input_layer");
0975   std::string outputLayer_HLIP = metadata_HLIP.get<std::string>("output_layer");
0976   tensorflow::run(session_HLIP, {{inputLayer_HLIP, input}}, {outputLayer_HLIP}, &outputs_HLIP);
0977   tensorflow::Tensor out_tensor_HLIP = outputs_HLIP[0];
0978   tensorflow::TTypes<float, 1>::Matrix dnn_outputs_HLIP = out_tensor_HLIP.matrix<float>();
0979 
0980   // Prepare tensor for DNN outputs
0981   std::vector<tensorflow::Tensor> outputs_HLMuS;
0982   // Evaluate DNN and put results in output tensor
0983   std::string inputLayer_HLMuS = metadata_HLMuS.get<std::string>("input_layer");
0984   std::string outputLayer_HLMuS = metadata_HLMuS.get<std::string>("output_layer");
0985   tensorflow::run(session_HLMuS, {{inputLayer_HLMuS, input}}, {outputLayer_HLMuS}, &outputs_HLMuS);
0986   tensorflow::Tensor out_tensor_HLMuS = outputs_HLMuS[0];
0987   tensorflow::TTypes<float, 1>::Matrix dnn_outputs_HLMuS = out_tensor_HLMuS.matrix<float>();
0988 
0989   // Decode output
0990   out.nHBd = round(dnn_outputs_HB(0, 0));
0991   out.nHLIP = round(dnn_outputs_HLIP(0, 0));
0992   out.sf = round(dnn_outputs_HLIP(0, 1));
0993   out.nHLMuS = round(dnn_outputs_HLMuS(0, 0));
0994 
0995   // Prevent prediction of negative number of seeds or too many seeds
0996   out.nHBd = std::clamp(out.nHBd, 0, 10);
0997   out.nHLIP = std::clamp(out.nHLIP, 0, 10);
0998   out.nHLMuS = std::clamp(out.nHLMuS, 0, 10);
0999 
1000   // Prevent prediction of 0 seeds in total
1001   if (out.nHBd == 0 && out.nHLIP == 0 && out.nHLMuS == 0) {
1002     // default strategy, similar to Run 2
1003     out.nHBd = 1;
1004     out.nHLIP = 5;
1005   }
1006 
1007   // Prevent extreme predictions for scale factors
1008   // (on average SF=2 was found to be optimal)
1009   if (out.sf <= 0)
1010     out.sf = 2;
1011   if (out.sf > 10)
1012     out.sf = 10;
1013 
1014   dnnSuccess = true;
1015 }
1016 
1017 //
1018 // Default values of configuration parameters
1019 //
1020 void TSGForOIDNN::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
1021   edm::ParameterSetDescription desc;
1022   desc.add<edm::InputTag>("src", edm::InputTag("hltL2Muons", "UpdatedAtVtx"));
1023   desc.add<int>("layersToTry", 2);
1024   desc.add<double>("fixedErrorRescaleFactorForHitless", 2.0);
1025   desc.add<int>("hitsToTry", 1);
1026   desc.add<edm::InputTag>("MeasurementTrackerEvent", edm::InputTag("hltSiStripClusters"));
1027   desc.add<std::string>("estimator", "hltESPChi2MeasurementEstimator100");
1028   desc.add<double>("maxEtaForTOB", 1.8);
1029   desc.add<double>("minEtaForTEC", 0.7);
1030   desc.addUntracked<bool>("debug", false);
1031   desc.add<unsigned int>("maxSeeds", 20);
1032   desc.add<unsigned int>("maxHitlessSeeds", 5);
1033   desc.add<unsigned int>("maxHitSeeds", 1);
1034   desc.add<std::string>("propagatorName", "PropagatorWithMaterialParabolicMf");
1035   desc.add<unsigned int>("maxHitlessSeedsIP", 5);
1036   desc.add<unsigned int>("maxHitlessSeedsMuS", 0);
1037   desc.add<unsigned int>("maxHitDoubletSeeds", 0);
1038   desc.add<bool>("getStrategyFromDNN", false);
1039   desc.add<bool>("useRegressor", false);
1040   desc.add<std::string>("dnnMetadataPath", "");
1041   descriptions.add("tsgForOIDNN", desc);
1042 }
1043 
1044 DEFINE_FWK_MODULE(TSGForOIDNN);