Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-24 22:51:38 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!