File indexing completed on 2024-04-06 12:27:15
0001
0002
0003
0004
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
0050 const edm::EDGetTokenT<reco::TrackCollection> src_;
0051
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
0065 const unsigned int maxSeeds_;
0066
0067 const unsigned int maxHitSeeds_;
0068
0069 const unsigned int maxHitlessSeeds_;
0070
0071 const unsigned int numOfLayersToTry_;
0072
0073 const unsigned int numOfHitsToTry_;
0074
0075 const double fixedErrorRescalingForHitless_;
0076
0077 const std::string estimatorName_;
0078
0079 const double minEtaForTEC_;
0080
0081 const double maxEtaForTOB_;
0082
0083
0084
0085 const unsigned int maxHitlessSeedsIP_;
0086 const unsigned int maxHitlessSeedsMuS_;
0087 const unsigned int maxHitDoubletSeeds_;
0088
0089
0090 const bool getStrategyFromDNN_;
0091
0092 const bool useRegressor_;
0093
0094
0095 std::string dnnModelPath_;
0096 std::unique_ptr<tensorflow::GraphDef> graphDef_;
0097 tensorflow::Session* tf_session_;
0098
0099
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
0111 const std::string dnnMetadataPath_;
0112 pt::ptree metadata_;
0113
0114
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
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
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
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
0158 struct StrategyParameters {
0159 int nHBd, nHLIP, nHLMuS, sf;
0160 };
0161
0162
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
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
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
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 }
0286
0287 void TSGForOIDNN::initializeTensorflow() const {
0288 if (getStrategyFromDNN_ and not useRegressor_) {
0289
0290 StrategyParameters strPars;
0291 bool dnnSuccess = false;
0292 evaluateClassifier(dummyFeatureMap(), tf_session_, metadata_, strPars, dnnSuccess);
0293 }
0294 }
0295
0296
0297
0298
0299 void TSGForOIDNN::produce(edm::StreamID sid, edm::Event& iEvent, edm::EventSetup const& iEventSetup) const {
0300
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
0309 std::unordered_map<std::string, float> feature_map;
0310
0311
0312 StrategyParameters strPars;
0313
0314
0315 Plane::PlanePointer dummyPlane = Plane::build(Plane::PositionType(), Plane::RotationType());
0316
0317
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
0328 auto const& l2TrackCol = iEvent.get(src_);
0329
0330
0331 std::unique_ptr<std::vector<TrajectorySeed> > result(new std::vector<TrajectorySeed>());
0332
0333
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
0342 std::unique_ptr<Propagator> propagatorAlong = SetPropagationDirection(tmpPropagatorAlong, alongMomentum);
0343 std::unique_ptr<Propagator> propagatorOpposite = SetPropagationDirection(tmpPropagatorOpposite, oppositeToMomentum);
0344
0345
0346 edm::ESHandle<Propagator> shpOpposite = iEventSetup.getHandle(t_SHPOpposite_);
0347
0348
0349 LogTrace(theCategory_) << "TSGForOIDNN::produce: Number of L2's: " << l2TrackCol.size();
0350
0351 for (auto const& l2 : l2TrackCol) {
0352
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
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
0377
0378 double l2muonEta = l2.eta();
0379 double absL2muonEta = std::abs(l2muonEta);
0380
0381
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
0390 if (getStrategyFromDNN_) {
0391 bool dnnSuccess = false;
0392
0393
0394 updateFeatureMap(feature_map, l2, tsosAtIP, outerTkStateOutside);
0395
0396 if (useRegressor_) {
0397
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
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
0478 if (absL2muonEta < maxEtaForTOB_) {
0479 layerCount = 0;
0480 createSeeds(tob);
0481 LogTrace("TSGForOIDNN") << "TSGForOIDNN:::produce: NumSeedsMade = " << numSeedsMade
0482 << " , layerCount = " << layerCount;
0483 }
0484
0485
0486 if (absL2muonEta > minEtaForTEC_ && absL2muonEta < maxEtaForTOB_) {
0487 numSeedsMade = 0;
0488 hitlessSeedsMadeIP = 0;
0489 hitlessSeedsMadeMuS = 0;
0490 hitSeedsMade = 0;
0491 hitDoubletSeedsMade = 0;
0492 }
0493
0494
0495 if (l2muonEta > minEtaForTEC_) {
0496 layerCount = 0;
0497 createSeeds(tecPositive);
0498 LogTrace("TSGForOIDNN") << "TSGForOIDNN:::produce: NumSeedsMade = " << numSeedsMade
0499 << " , layerCount = " << layerCount;
0500 }
0501
0502
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 }
0515
0516 edm::LogInfo(theCategory_) << "TSGForOIDNN::produce: number of seeds made: " << result->size();
0517
0518 iEvent.put(std::move(result));
0519 }
0520
0521
0522
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
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
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
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;
0585
0586 std::vector<TrajectoryMeasurement> mymeas =
0587 det.fastMeasurements(detI.second, onLayer, propagatorAlong, estimator);
0588 for (auto const& measurement : mymeas) {
0589 if (measurement.recHit()->isValid())
0590 meas.push_back(measurement);
0591 }
0592 }
0593
0594
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
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
0642
0643
0644
0645
0646
0647 int max_addtnl_layers = 1;
0648 int max_meas = 1;
0649
0650
0651
0652 const TrajectoryStateOnSurface& onLayer(tsos);
0653
0654
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
0663 for (auto const& detI : dets) {
0664 MeasurementDetWithData det = measurementTracker.idToDet(detI.first->geographicalId());
0665
0666 if (det.isNull())
0667 continue;
0668 if (!detI.second.isValid())
0669 continue;
0670
0671
0672 std::vector<TrajectoryMeasurement> mymeas = det.fastMeasurements(detI.second, onLayer, propagatorAlong, estimator);
0673
0674
0675 for (auto const& measurement : mymeas) {
0676 if (measurement.recHit()->isValid())
0677 meas.push_back(measurement);
0678 }
0679 }
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
0686 std::sort(meas.begin(), meas.end(), TrajMeasLessEstim());
0687
0688 unsigned int found = 0;
0689
0690
0691
0692 for (auto const& measurement : meas) {
0693 if (hitDoubletSeedsMade >= maxHitDoubletSeeds)
0694 return;
0695
0696
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;
0702
0703 edm::OwnVector<TrackingRecHit> seedHits;
0704
0705
0706 seedHits.push_back(*measurement.recHit()->hit());
0707 const DetLayer* detLayer = dynamic_cast<const DetLayer*>(&layer);
0708
0709
0710
0711
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
0719 TrajectoryStateOnSurface updatedTSOS_next(updatedTSOS);
0720
0721
0722 for (auto compLayer : compLayers) {
0723 int nmeas = 0;
0724
0725 if (addtnl_layers_scanned >= max_addtnl_layers)
0726 break;
0727 if (found_compatible_on_next_layer > 0)
0728 break;
0729
0730
0731 std::vector<GeometricSearchDet::DetWithState> dets_next;
0732 TrajectoryStateOnSurface onLayer_next(updatedTSOS);
0733
0734 compLayer->compatibleDetsV(onLayer_next, propagatorAlong, estimator, dets_next);
0735
0736
0737 std::vector<TrajectoryMeasurement> meas_next;
0738
0739
0740 for (auto const& detI_next : dets_next) {
0741 MeasurementDetWithData det = measurementTracker.idToDet(detI_next.first->geographicalId());
0742
0743 if (det.isNull())
0744 continue;
0745 if (!detI_next.second.isValid())
0746 continue;
0747
0748
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
0754 if (mea_next.recHit()->isValid())
0755 meas_next.push_back(mea_next);
0756
0757 }
0758 }
0759
0760
0761 std::sort(meas_next.begin(), meas_next.end(), TrajMeasLessEstim());
0762
0763
0764 for (auto const& mea_next : meas_next) {
0765 if (nmeas >= max_meas)
0766 break;
0767
0768
0769 updatedTSOS_next = updator_->update(mea_next.forwardPredictedState(), *mea_next.recHit());
0770
0771 if (not updatedTSOS_next.isValid())
0772 continue;
0773
0774
0775
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 }
0782
0783 addtnl_layers_scanned++;
0784
0785 }
0786
0787 if (found_compatible_on_next_layer == 0)
0788 continue;
0789
0790
0791
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;
0804
0805 }
0806
0807 if (found)
0808 layerCount++;
0809 }
0810
0811
0812
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
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
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
0894 dnnSuccess = false;
0895 } else {
0896 input.matrix<float>()(0, i_feature) = float(feature_map.at(fname));
0897 i_feature++;
0898 }
0899 }
0900
0901
0902 std::vector<tensorflow::Tensor> outputs;
0903
0904
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
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
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
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
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
0955 dnnSuccess = false;
0956 } else {
0957 input.matrix<float>()(0, i_feature) = float(feature_map.at(fname));
0958 i_feature++;
0959 }
0960 }
0961
0962
0963 std::vector<tensorflow::Tensor> outputs_HB;
0964
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
0972 std::vector<tensorflow::Tensor> outputs_HLIP;
0973
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
0981 std::vector<tensorflow::Tensor> outputs_HLMuS;
0982
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
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
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
1001 if (out.nHBd == 0 && out.nHLIP == 0 && out.nHLMuS == 0) {
1002
1003 out.nHBd = 1;
1004 out.nHLIP = 5;
1005 }
1006
1007
1008
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
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);