File indexing completed on 2022-04-07 05:53:37
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
0047 private:
0048
0049 const edm::EDGetTokenT<reco::TrackCollection> src_;
0050
0051 const edm::ESGetToken<Chi2MeasurementEstimatorBase, TrackingComponentsRecord> t_estimatorH_;
0052 const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> t_magfieldH_;
0053 const edm::ESGetToken<Propagator, TrackingComponentsRecord> t_propagatorAlongH_;
0054 const edm::ESGetToken<Propagator, TrackingComponentsRecord> t_propagatorOppositeH_;
0055 const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> t_tmpTkGeometryH_;
0056 const edm::ESGetToken<GlobalTrackingGeometry, GlobalTrackingGeometryRecord> t_geometryH_;
0057 const edm::ESGetToken<NavigationSchool, NavigationSchoolRecord> t_navSchool_;
0058 const edm::ESGetToken<Propagator, TrackingComponentsRecord> t_SHPOpposite_;
0059 const std::unique_ptr<TrajectoryStateUpdator> updator_;
0060 const edm::EDGetTokenT<MeasurementTrackerEvent> measurementTrackerTag_;
0061 const std::string theCategory_;
0062
0063
0064 const unsigned int maxSeeds_;
0065
0066 const unsigned int maxHitSeeds_;
0067
0068 const unsigned int maxHitlessSeeds_;
0069
0070 const unsigned int numOfLayersToTry_;
0071
0072 const unsigned int numOfHitsToTry_;
0073
0074 const double fixedErrorRescalingForHitless_;
0075
0076 const std::string estimatorName_;
0077
0078 const double minEtaForTEC_;
0079
0080 const double maxEtaForTOB_;
0081
0082
0083
0084 const unsigned int maxHitlessSeedsIP_;
0085 const unsigned int maxHitlessSeedsMuS_;
0086 const unsigned int maxHitDoubletSeeds_;
0087
0088
0089 const bool getStrategyFromDNN_;
0090
0091 const bool useRegressor_;
0092
0093
0094 std::string dnnModelPath_;
0095 std::unique_ptr<tensorflow::GraphDef> graphDef_;
0096 tensorflow::Session* tf_session_;
0097
0098
0099 std::string dnnModelPath_HB_;
0100 std::string dnnModelPath_HLIP_;
0101 std::string dnnModelPath_HLMuS_;
0102 std::unique_ptr<tensorflow::GraphDef> graphDef_HB_;
0103 tensorflow::Session* tf_session_HB_;
0104 std::unique_ptr<tensorflow::GraphDef> graphDef_HLIP_;
0105 tensorflow::Session* tf_session_HLIP_;
0106 std::unique_ptr<tensorflow::GraphDef> graphDef_HLMuS_;
0107 tensorflow::Session* tf_session_HLMuS_;
0108
0109
0110 const std::string dnnMetadataPath_;
0111 pt::ptree metadata_;
0112
0113
0114 void makeSeedsWithoutHits(const GeometricSearchDet& layer,
0115 const TrajectoryStateOnSurface& tsos,
0116 const Propagator& propagatorAlong,
0117 const Chi2MeasurementEstimatorBase& estimator,
0118 double errorSF,
0119 unsigned int& hitlessSeedsMade,
0120 unsigned int& numSeedsMade,
0121 std::vector<TrajectorySeed>& out) const;
0122
0123
0124 void makeSeedsFromHits(const GeometricSearchDet& layer,
0125 const TrajectoryStateOnSurface& tsos,
0126 const Propagator& propagatorAlong,
0127 const Chi2MeasurementEstimatorBase& estimator,
0128 const MeasurementTrackerEvent& measurementTracker,
0129 unsigned int& hitSeedsMade,
0130 unsigned int& numSeedsMade,
0131 const unsigned int& maxHitSeeds,
0132 unsigned int& layerCount,
0133 std::vector<TrajectorySeed>& out) const;
0134
0135
0136 void makeSeedsFromHitDoublets(const GeometricSearchDet& layer,
0137 const TrajectoryStateOnSurface& tsos,
0138 const Propagator& propagatorAlong,
0139 const Chi2MeasurementEstimatorBase& estimator,
0140 const MeasurementTrackerEvent& measurementTracker,
0141 const NavigationSchool& navSchool,
0142 unsigned int& hitDoubletSeedsMade,
0143 unsigned int& numSeedsMade,
0144 const unsigned int& maxHitDoubletSeeds,
0145 unsigned int& layerCount,
0146 std::vector<TrajectorySeed>& out) const;
0147
0148
0149 void updateFeatureMap(std::unordered_map<std::string, float>& the_map,
0150 const reco::Track& l2,
0151 const TrajectoryStateOnSurface& tsos_IP,
0152 const TrajectoryStateOnSurface& tsos_MuS) const;
0153
0154
0155 struct StrategyParameters {
0156 int nHBd, nHLIP, nHLMuS, sf;
0157 };
0158
0159
0160 void evaluateClassifier(const std::unordered_map<std::string, float>& feature_map,
0161 tensorflow::Session* session,
0162 const pt::ptree& metadata,
0163 StrategyParameters& out,
0164 bool& dnnSuccess) const;
0165
0166
0167 void evaluateRegressor(const std::unordered_map<std::string, float>& feature_map,
0168 tensorflow::Session* session_HB,
0169 const pt::ptree& metadata_HB,
0170 tensorflow::Session* session_HLIP,
0171 const pt::ptree& metadata_HLIP,
0172 tensorflow::Session* session_HLMuS,
0173 const pt::ptree& metadata_HLMuS,
0174 StrategyParameters& out,
0175 bool& dnnSuccess) const;
0176 };
0177
0178 TSGForOIDNN::TSGForOIDNN(const edm::ParameterSet& iConfig)
0179 : src_(consumes(iConfig.getParameter<edm::InputTag>("src"))),
0180 t_estimatorH_(esConsumes(edm::ESInputTag("", iConfig.getParameter<std::string>("estimator")))),
0181 t_magfieldH_(esConsumes()),
0182 t_propagatorAlongH_(esConsumes(edm::ESInputTag("", iConfig.getParameter<std::string>("propagatorName")))),
0183 t_propagatorOppositeH_(esConsumes(edm::ESInputTag("", iConfig.getParameter<std::string>("propagatorName")))),
0184 t_tmpTkGeometryH_(esConsumes()),
0185 t_geometryH_(esConsumes()),
0186 t_navSchool_(esConsumes(edm::ESInputTag("", "SimpleNavigationSchool"))),
0187 t_SHPOpposite_(esConsumes(edm::ESInputTag("", "hltESPSteppingHelixPropagatorOpposite"))),
0188 updator_(new KFUpdator()),
0189 measurementTrackerTag_(consumes(iConfig.getParameter<edm::InputTag>("MeasurementTrackerEvent"))),
0190 theCategory_(std::string("Muon|RecoMuon|TSGForOIDNN")),
0191 maxSeeds_(iConfig.getParameter<uint32_t>("maxSeeds")),
0192 maxHitSeeds_(iConfig.getParameter<uint32_t>("maxHitSeeds")),
0193 maxHitlessSeeds_(iConfig.getParameter<uint32_t>("maxHitlessSeeds")),
0194 numOfLayersToTry_(iConfig.getParameter<int32_t>("layersToTry")),
0195 numOfHitsToTry_(iConfig.getParameter<int32_t>("hitsToTry")),
0196 fixedErrorRescalingForHitless_(iConfig.getParameter<double>("fixedErrorRescaleFactorForHitless")),
0197 minEtaForTEC_(iConfig.getParameter<double>("minEtaForTEC")),
0198 maxEtaForTOB_(iConfig.getParameter<double>("maxEtaForTOB")),
0199 maxHitlessSeedsIP_(iConfig.getParameter<uint32_t>("maxHitlessSeedsIP")),
0200 maxHitlessSeedsMuS_(iConfig.getParameter<uint32_t>("maxHitlessSeedsMuS")),
0201 maxHitDoubletSeeds_(iConfig.getParameter<uint32_t>("maxHitDoubletSeeds")),
0202 getStrategyFromDNN_(iConfig.getParameter<bool>("getStrategyFromDNN")),
0203 useRegressor_(iConfig.getParameter<bool>("useRegressor")),
0204 dnnMetadataPath_(iConfig.getParameter<std::string>("dnnMetadataPath")) {
0205 if (getStrategyFromDNN_) {
0206 edm::FileInPath dnnMetadataPath(dnnMetadataPath_);
0207 pt::read_json(dnnMetadataPath.fullPath(), metadata_);
0208 tensorflow::setLogging("3");
0209
0210 if (useRegressor_) {
0211
0212 dnnModelPath_HB_ = metadata_.get<std::string>("HB.dnnmodel_path");
0213 edm::FileInPath dnnPath_HB(dnnModelPath_HB_);
0214 graphDef_HB_ = std::unique_ptr<tensorflow::GraphDef>(tensorflow::loadGraphDef(dnnPath_HB.fullPath()));
0215 tf_session_HB_ = tensorflow::createSession(graphDef_HB_.get());
0216
0217 dnnModelPath_HLIP_ = metadata_.get<std::string>("HLIP.dnnmodel_path");
0218 edm::FileInPath dnnPath_HLIP(dnnModelPath_HLIP_);
0219 graphDef_HLIP_ = std::unique_ptr<tensorflow::GraphDef>(tensorflow::loadGraphDef(dnnPath_HLIP.fullPath()));
0220 tf_session_HLIP_ = tensorflow::createSession(graphDef_HLIP_.get());
0221
0222 dnnModelPath_HLMuS_ = metadata_.get<std::string>("HLMuS.dnnmodel_path");
0223 edm::FileInPath dnnPath_HLMuS(dnnModelPath_HLMuS_);
0224 graphDef_HLMuS_ = std::unique_ptr<tensorflow::GraphDef>(tensorflow::loadGraphDef(dnnPath_HLMuS.fullPath()));
0225 tf_session_HLMuS_ = tensorflow::createSession(graphDef_HLMuS_.get());
0226 } else {
0227
0228 dnnModelPath_ = metadata_.get<std::string>("dnnmodel_path");
0229 edm::FileInPath dnnPath(dnnModelPath_);
0230 graphDef_ = std::unique_ptr<tensorflow::GraphDef>(tensorflow::loadGraphDef(dnnPath.fullPath()));
0231 tf_session_ = tensorflow::createSession(graphDef_.get());
0232 }
0233 }
0234 produces<std::vector<TrajectorySeed> >();
0235 }
0236
0237 TSGForOIDNN::~TSGForOIDNN() {
0238 if (getStrategyFromDNN_) {
0239 if (useRegressor_) {
0240 tensorflow::closeSession(tf_session_HB_);
0241 tensorflow::closeSession(tf_session_HLIP_);
0242 tensorflow::closeSession(tf_session_HLMuS_);
0243 } else {
0244 tensorflow::closeSession(tf_session_);
0245 }
0246 }
0247 }
0248
0249
0250
0251
0252 void TSGForOIDNN::produce(edm::StreamID sid, edm::Event& iEvent, edm::EventSetup const& iEventSetup) const {
0253
0254 unsigned int numSeedsMade = 0;
0255 unsigned int layerCount = 0;
0256 unsigned int hitlessSeedsMadeIP = 0;
0257 unsigned int hitlessSeedsMadeMuS = 0;
0258 unsigned int hitSeedsMade = 0;
0259 unsigned int hitDoubletSeedsMade = 0;
0260
0261
0262 std::unordered_map<std::string, float> feature_map;
0263
0264
0265 StrategyParameters strPars;
0266
0267
0268 Plane::PlanePointer dummyPlane = Plane::build(Plane::PositionType(), Plane::RotationType());
0269
0270
0271 const MagneticField& magfield = iEventSetup.getData(t_magfieldH_);
0272 const Chi2MeasurementEstimatorBase& estimator = iEventSetup.getData(t_estimatorH_);
0273 const Propagator& tmpPropagatorAlong = iEventSetup.getData(t_propagatorAlongH_);
0274 const Propagator& tmpPropagatorOpposite = iEventSetup.getData(t_propagatorOppositeH_);
0275 const TrackerGeometry& tmpTkGeometry = iEventSetup.getData(t_tmpTkGeometryH_);
0276 const GlobalTrackingGeometry& geometry = iEventSetup.getData(t_geometryH_);
0277 const NavigationSchool& navSchool = iEventSetup.getData(t_navSchool_);
0278 auto const& measurementTracker = iEvent.get(measurementTrackerTag_);
0279
0280
0281 auto const& l2TrackCol = iEvent.get(src_);
0282
0283
0284 std::unique_ptr<std::vector<TrajectorySeed> > result(new std::vector<TrajectorySeed>());
0285
0286
0287 auto const* gsTracker = measurementTracker.geometricSearchTracker();
0288 std::vector<BarrelDetLayer const*> const& tob = gsTracker->tobLayers();
0289 std::vector<ForwardDetLayer const*> const& tecPositive =
0290 tmpTkGeometry.isThere(GeomDetEnumerators::P2OTEC) ? gsTracker->posTidLayers() : gsTracker->posTecLayers();
0291 std::vector<ForwardDetLayer const*> const& tecNegative =
0292 tmpTkGeometry.isThere(GeomDetEnumerators::P2OTEC) ? gsTracker->negTidLayers() : gsTracker->negTecLayers();
0293
0294
0295 std::unique_ptr<Propagator> propagatorAlong = SetPropagationDirection(tmpPropagatorAlong, alongMomentum);
0296 std::unique_ptr<Propagator> propagatorOpposite = SetPropagationDirection(tmpPropagatorOpposite, oppositeToMomentum);
0297
0298
0299 edm::ESHandle<Propagator> shpOpposite = iEventSetup.getHandle(t_SHPOpposite_);
0300
0301
0302 LogTrace(theCategory_) << "TSGForOIDNN::produce: Number of L2's: " << l2TrackCol.size();
0303
0304 for (auto const& l2 : l2TrackCol) {
0305
0306 std::vector<TrajectorySeed> out;
0307 LogTrace("TSGForOIDNN") << "TSGForOIDNN::produce: L2 muon pT, eta, phi --> " << l2.pt() << " , " << l2.eta()
0308 << " , " << l2.phi();
0309
0310 FreeTrajectoryState fts = trajectoryStateTransform::initialFreeState(l2, &magfield);
0311
0312 dummyPlane->move(fts.position() - dummyPlane->position());
0313 TrajectoryStateOnSurface tsosAtIP = TrajectoryStateOnSurface(fts, *dummyPlane);
0314 LogTrace("TSGForOIDNN") << "TSGForOIDNN::produce: Created TSOSatIP: " << tsosAtIP;
0315
0316
0317 TrajectoryStateOnSurface tsosAtMuonSystem = trajectoryStateTransform::innerStateOnSurface(l2, geometry, &magfield);
0318 LogTrace("TSGForOIDNN") << "TSGForOIDNN::produce: Created TSOSatMuonSystem: " << tsosAtMuonSystem;
0319
0320 LogTrace("TSGForOIDNN")
0321 << "TSGForOIDNN::produce: Check the error of the L2 parameter and use hit seeds if big errors";
0322
0323 StateOnTrackerBound fromInside(propagatorAlong.get());
0324 TrajectoryStateOnSurface outerTkStateInside = fromInside(fts);
0325
0326 StateOnTrackerBound fromOutside(&*shpOpposite);
0327 TrajectoryStateOnSurface outerTkStateOutside = fromOutside(tsosAtMuonSystem);
0328
0329
0330
0331 double l2muonEta = l2.eta();
0332 double absL2muonEta = std::abs(l2muonEta);
0333
0334
0335 unsigned int maxHitSeeds = maxHitSeeds_;
0336 unsigned int maxHitDoubletSeeds = maxHitDoubletSeeds_;
0337 unsigned int maxHitlessSeedsIP = maxHitlessSeedsIP_;
0338 unsigned int maxHitlessSeedsMuS = maxHitlessSeedsMuS_;
0339
0340 float errorSFHitless = fixedErrorRescalingForHitless_;
0341
0342
0343 if (getStrategyFromDNN_) {
0344 bool dnnSuccess = false;
0345
0346
0347 updateFeatureMap(feature_map, l2, tsosAtIP, outerTkStateOutside);
0348
0349 if (useRegressor_) {
0350
0351 evaluateRegressor(feature_map,
0352 tf_session_HB_,
0353 metadata_.get_child("HB"),
0354 tf_session_HLIP_,
0355 metadata_.get_child("HLIP"),
0356 tf_session_HLMuS_,
0357 metadata_.get_child("HLMuS"),
0358 strPars,
0359 dnnSuccess);
0360 } else {
0361
0362 evaluateClassifier(feature_map, tf_session_, metadata_, strPars, dnnSuccess);
0363 }
0364 if (!dnnSuccess)
0365 break;
0366
0367 maxHitSeeds = 0;
0368 maxHitDoubletSeeds = strPars.nHBd;
0369 maxHitlessSeedsIP = strPars.nHLIP;
0370 maxHitlessSeedsMuS = strPars.nHLMuS;
0371 errorSFHitless = strPars.sf;
0372 }
0373
0374 numSeedsMade = 0;
0375 hitlessSeedsMadeIP = 0;
0376 hitlessSeedsMadeMuS = 0;
0377 hitSeedsMade = 0;
0378 hitDoubletSeedsMade = 0;
0379
0380 auto createSeeds = [&](auto const& layers) {
0381 for (auto const& layer : boost::adaptors::reverse(layers)) {
0382 if (hitlessSeedsMadeIP < maxHitlessSeedsIP && numSeedsMade < maxSeeds_)
0383 makeSeedsWithoutHits(*layer,
0384 tsosAtIP,
0385 *(propagatorAlong.get()),
0386 estimator,
0387 errorSFHitless,
0388 hitlessSeedsMadeIP,
0389 numSeedsMade,
0390 out);
0391
0392 if (outerTkStateInside.isValid() && outerTkStateOutside.isValid() && hitlessSeedsMadeMuS < maxHitlessSeedsMuS &&
0393 numSeedsMade < maxSeeds_)
0394 makeSeedsWithoutHits(*layer,
0395 outerTkStateOutside,
0396 *(propagatorOpposite.get()),
0397 estimator,
0398 errorSFHitless,
0399 hitlessSeedsMadeMuS,
0400 numSeedsMade,
0401 out);
0402
0403 if (hitSeedsMade < maxHitSeeds && numSeedsMade < maxSeeds_)
0404 makeSeedsFromHits(*layer,
0405 tsosAtIP,
0406 *(propagatorAlong.get()),
0407 estimator,
0408 measurementTracker,
0409 hitSeedsMade,
0410 numSeedsMade,
0411 maxHitSeeds,
0412 layerCount,
0413 out);
0414
0415 if (hitDoubletSeedsMade < maxHitDoubletSeeds && numSeedsMade < maxSeeds_)
0416 makeSeedsFromHitDoublets(*layer,
0417 tsosAtIP,
0418 *(propagatorAlong.get()),
0419 estimator,
0420 measurementTracker,
0421 navSchool,
0422 hitDoubletSeedsMade,
0423 numSeedsMade,
0424 maxHitDoubletSeeds,
0425 layerCount,
0426 out);
0427 };
0428 };
0429
0430
0431 if (absL2muonEta < maxEtaForTOB_) {
0432 layerCount = 0;
0433 createSeeds(tob);
0434 LogTrace("TSGForOIDNN") << "TSGForOIDNN:::produce: NumSeedsMade = " << numSeedsMade
0435 << " , layerCount = " << layerCount;
0436 }
0437
0438
0439 if (absL2muonEta > minEtaForTEC_ && absL2muonEta < maxEtaForTOB_) {
0440 numSeedsMade = 0;
0441 hitlessSeedsMadeIP = 0;
0442 hitlessSeedsMadeMuS = 0;
0443 hitSeedsMade = 0;
0444 hitDoubletSeedsMade = 0;
0445 }
0446
0447
0448 if (l2muonEta > minEtaForTEC_) {
0449 layerCount = 0;
0450 createSeeds(tecPositive);
0451 LogTrace("TSGForOIDNN") << "TSGForOIDNN:::produce: NumSeedsMade = " << numSeedsMade
0452 << " , layerCount = " << layerCount;
0453 }
0454
0455
0456 if (l2muonEta < -minEtaForTEC_) {
0457 layerCount = 0;
0458 createSeeds(tecNegative);
0459 LogTrace("TSGForOIDNN") << "TSGForOIDNN:::produce: NumSeedsMade = " << numSeedsMade
0460 << " , layerCount = " << layerCount;
0461 }
0462
0463 for (std::vector<TrajectorySeed>::iterator it = out.begin(); it != out.end(); ++it) {
0464 result->push_back(*it);
0465 }
0466
0467 }
0468
0469 edm::LogInfo(theCategory_) << "TSGForOIDNN::produce: number of seeds made: " << result->size();
0470
0471 iEvent.put(std::move(result));
0472 }
0473
0474
0475
0476
0477 void TSGForOIDNN::makeSeedsWithoutHits(const GeometricSearchDet& layer,
0478 const TrajectoryStateOnSurface& tsos,
0479 const Propagator& propagatorAlong,
0480 const Chi2MeasurementEstimatorBase& estimator,
0481 double errorSF,
0482 unsigned int& hitlessSeedsMade,
0483 unsigned int& numSeedsMade,
0484 std::vector<TrajectorySeed>& out) const {
0485
0486 LogTrace("TSGForOIDNN") << "TSGForOIDNN::makeSeedsWithoutHits: Start hitless";
0487 std::vector<GeometricSearchDet::DetWithState> dets;
0488 layer.compatibleDetsV(tsos, propagatorAlong, estimator, dets);
0489 if (!dets.empty()) {
0490 auto const& detOnLayer = dets.front().first;
0491 auto& tsosOnLayer = dets.front().second;
0492 LogTrace("TSGForOIDNN") << "TSGForOIDNN::makeSeedsWithoutHits: tsosOnLayer " << tsosOnLayer;
0493 if (!tsosOnLayer.isValid()) {
0494 edm::LogInfo(theCategory_) << "ERROR!: Hitless TSOS is not valid!";
0495 } else {
0496 tsosOnLayer.rescaleError(errorSF);
0497 PTrajectoryStateOnDet const& ptsod =
0498 trajectoryStateTransform::persistentState(tsosOnLayer, detOnLayer->geographicalId().rawId());
0499 TrajectorySeed::RecHitContainer rHC;
0500 out.push_back(TrajectorySeed(ptsod, rHC, oppositeToMomentum));
0501 LogTrace("TSGForOIDNN") << "TSGForOIDNN::makeSeedsWithoutHits: TSOS (Hitless) done ";
0502 hitlessSeedsMade++;
0503 numSeedsMade++;
0504 }
0505 }
0506 }
0507
0508
0509
0510
0511 void TSGForOIDNN::makeSeedsFromHits(const GeometricSearchDet& layer,
0512 const TrajectoryStateOnSurface& tsos,
0513 const Propagator& propagatorAlong,
0514 const Chi2MeasurementEstimatorBase& estimator,
0515 const MeasurementTrackerEvent& measurementTracker,
0516 unsigned int& hitSeedsMade,
0517 unsigned int& numSeedsMade,
0518 const unsigned int& maxHitSeeds,
0519 unsigned int& layerCount,
0520 std::vector<TrajectorySeed>& out) const {
0521 if (layerCount > numOfLayersToTry_)
0522 return;
0523
0524 const TrajectoryStateOnSurface& onLayer(tsos);
0525
0526 std::vector<GeometricSearchDet::DetWithState> dets;
0527 layer.compatibleDetsV(onLayer, propagatorAlong, estimator, dets);
0528
0529
0530 LogTrace("TSGForOIDNN") << "TSGForOIDNN::makeSeedsFromHits: Find measurements on each detWithState " << dets.size();
0531 std::vector<TrajectoryMeasurement> meas;
0532 for (auto const& detI : dets) {
0533 MeasurementDetWithData det = measurementTracker.idToDet(detI.first->geographicalId());
0534 if (det.isNull())
0535 continue;
0536 if (!detI.second.isValid())
0537 continue;
0538
0539 std::vector<TrajectoryMeasurement> mymeas =
0540 det.fastMeasurements(detI.second, onLayer, propagatorAlong, estimator);
0541 for (auto const& measurement : mymeas) {
0542 if (measurement.recHit()->isValid())
0543 meas.push_back(measurement);
0544 }
0545 }
0546
0547
0548 LogTrace("TSGForOIDNN") << "TSGForOIDNN::makeSeedsFromHits: Update TSOS using TMs after sorting, then create "
0549 "Trajectory Seed, number of TM = "
0550 << meas.size();
0551 std::sort(meas.begin(), meas.end(), TrajMeasLessEstim());
0552
0553 unsigned int found = 0;
0554 for (auto const& measurement : meas) {
0555 if (hitSeedsMade >= maxHitSeeds)
0556 return;
0557 TrajectoryStateOnSurface updatedTSOS = updator_->update(measurement.forwardPredictedState(), *measurement.recHit());
0558 LogTrace("TSGForOIDNN") << "TSGForOIDNN::makeSeedsFromHits: TSOS for TM " << found;
0559 if (not updatedTSOS.isValid())
0560 continue;
0561
0562 edm::OwnVector<TrackingRecHit> seedHits;
0563 seedHits.push_back(*measurement.recHit()->hit());
0564 PTrajectoryStateOnDet const& pstate =
0565 trajectoryStateTransform::persistentState(updatedTSOS, measurement.recHit()->geographicalId().rawId());
0566 LogTrace("TSGForOIDNN") << "TSGForOIDNN::makeSeedsFromHits: Number of seedHits: " << seedHits.size();
0567 TrajectorySeed seed(pstate, std::move(seedHits), oppositeToMomentum);
0568 out.push_back(seed);
0569 found++;
0570 numSeedsMade++;
0571 hitSeedsMade++;
0572 if (found == numOfHitsToTry_)
0573 break;
0574 }
0575
0576 if (found)
0577 layerCount++;
0578 }
0579
0580
0581
0582
0583 void TSGForOIDNN::makeSeedsFromHitDoublets(const GeometricSearchDet& layer,
0584 const TrajectoryStateOnSurface& tsos,
0585 const Propagator& propagatorAlong,
0586 const Chi2MeasurementEstimatorBase& estimator,
0587 const MeasurementTrackerEvent& measurementTracker,
0588 const NavigationSchool& navSchool,
0589 unsigned int& hitDoubletSeedsMade,
0590 unsigned int& numSeedsMade,
0591 const unsigned int& maxHitDoubletSeeds,
0592 unsigned int& layerCount,
0593 std::vector<TrajectorySeed>& out) const {
0594
0595
0596
0597
0598
0599
0600 int max_addtnl_layers = 1;
0601 int max_meas = 1;
0602
0603
0604
0605 const TrajectoryStateOnSurface& onLayer(tsos);
0606
0607
0608 std::vector<GeometricSearchDet::DetWithState> dets;
0609 layer.compatibleDetsV(onLayer, propagatorAlong, estimator, dets);
0610
0611 LogTrace("TSGForOIDNN") << "TSGForOIDNN::makeSeedsFromHitDoublets: Find measurements on each detWithState "
0612 << dets.size();
0613 std::vector<TrajectoryMeasurement> meas;
0614
0615
0616 for (auto const& detI : dets) {
0617 MeasurementDetWithData det = measurementTracker.idToDet(detI.first->geographicalId());
0618
0619 if (det.isNull())
0620 continue;
0621 if (!detI.second.isValid())
0622 continue;
0623
0624
0625 std::vector<TrajectoryMeasurement> mymeas = det.fastMeasurements(detI.second, onLayer, propagatorAlong, estimator);
0626
0627
0628 for (auto const& measurement : mymeas) {
0629 if (measurement.recHit()->isValid())
0630 meas.push_back(measurement);
0631 }
0632 }
0633
0634 LogTrace("TSGForOIDNN") << "TSGForOIDNN::makeSeedsFromHitDoublets: Update TSOS using TMs after sorting, then create "
0635 "Trajectory Seed, number of TM = "
0636 << meas.size();
0637
0638
0639 std::sort(meas.begin(), meas.end(), TrajMeasLessEstim());
0640
0641 unsigned int found = 0;
0642 int hit_num = 0;
0643
0644
0645
0646 for (auto const& measurement : meas) {
0647 if (hitDoubletSeedsMade >= maxHitDoubletSeeds)
0648 return;
0649
0650 hit_num++;
0651
0652
0653 TrajectoryStateOnSurface updatedTSOS = updator_->update(measurement.forwardPredictedState(), *measurement.recHit());
0654
0655 LogTrace("TSGForOIDNN") << "TSGForOIDNN::makeSeedsFromHitDoublets: TSOS for TM " << found;
0656 if (not updatedTSOS.isValid())
0657 continue;
0658
0659 edm::OwnVector<TrackingRecHit> seedHits;
0660
0661
0662 seedHits.push_back(*measurement.recHit()->hit());
0663 const DetLayer* detLayer = dynamic_cast<const DetLayer*>(&layer);
0664
0665
0666
0667
0668 auto const& compLayers = navSchool.nextLayers(*detLayer, *updatedTSOS.freeState(), alongMomentum);
0669
0670 int addtnl_layers_scanned = 0;
0671 int found_compatible_on_next_layer = 0;
0672 int det_id = 0;
0673
0674
0675 TrajectoryStateOnSurface updatedTSOS_next(updatedTSOS);
0676
0677
0678 for (auto compLayer : compLayers) {
0679 int nmeas = 0;
0680
0681 if (addtnl_layers_scanned >= max_addtnl_layers)
0682 break;
0683 if (found_compatible_on_next_layer > 0)
0684 break;
0685
0686
0687 std::vector<GeometricSearchDet::DetWithState> dets_next;
0688 TrajectoryStateOnSurface onLayer_next(updatedTSOS);
0689
0690 compLayer->compatibleDetsV(onLayer_next, propagatorAlong, estimator, dets_next);
0691
0692
0693 std::vector<TrajectoryMeasurement> meas_next;
0694
0695
0696 for (auto const& detI_next : dets_next) {
0697 MeasurementDetWithData det = measurementTracker.idToDet(detI_next.first->geographicalId());
0698
0699 if (det.isNull())
0700 continue;
0701 if (!detI_next.second.isValid())
0702 continue;
0703
0704
0705 std::vector<TrajectoryMeasurement> mymeas_next =
0706 det.fastMeasurements(detI_next.second, onLayer_next, propagatorAlong, estimator);
0707
0708 for (auto const& mea_next : mymeas_next) {
0709
0710 if (mea_next.recHit()->isValid())
0711 meas_next.push_back(mea_next);
0712
0713 }
0714 }
0715
0716
0717 std::sort(meas_next.begin(), meas_next.end(), TrajMeasLessEstim());
0718
0719
0720 for (auto const& mea_next : meas_next) {
0721 if (nmeas >= max_meas)
0722 break;
0723
0724
0725 updatedTSOS_next = updator_->update(mea_next.forwardPredictedState(), *mea_next.recHit());
0726
0727 if (not updatedTSOS_next.isValid())
0728 continue;
0729
0730
0731
0732 seedHits.push_back(*mea_next.recHit()->hit());
0733 det_id = mea_next.recHit()->geographicalId().rawId();
0734 nmeas++;
0735 found_compatible_on_next_layer++;
0736
0737 }
0738
0739 addtnl_layers_scanned++;
0740
0741 }
0742
0743 if (found_compatible_on_next_layer == 0)
0744 continue;
0745
0746
0747
0748 PTrajectoryStateOnDet const& pstate = trajectoryStateTransform::persistentState(updatedTSOS_next, det_id);
0749 TrajectorySeed seed(pstate, std::move(seedHits), oppositeToMomentum);
0750
0751 LogTrace("TSGForOIDNN") << "TSGForOIDNN::makeSeedsFromHitDoublets: Number of seedHits: " << seedHits.size();
0752 out.push_back(seed);
0753
0754 found++;
0755 numSeedsMade++;
0756 hitDoubletSeedsMade++;
0757
0758 if (found == numOfHitsToTry_)
0759 break;
0760
0761 }
0762
0763 if (found)
0764 layerCount++;
0765 }
0766
0767
0768
0769
0770 void TSGForOIDNN::updateFeatureMap(std::unordered_map<std::string, float>& the_map,
0771 const reco::Track& l2,
0772 const TrajectoryStateOnSurface& tsos_IP,
0773 const TrajectoryStateOnSurface& tsos_MuS) const {
0774 the_map["pt"] = l2.pt();
0775 the_map["eta"] = l2.eta();
0776 the_map["phi"] = l2.phi();
0777 the_map["validHits"] = l2.found();
0778 if (tsos_IP.isValid()) {
0779 the_map["tsos_IP_eta"] = tsos_IP.globalPosition().eta();
0780 the_map["tsos_IP_phi"] = tsos_IP.globalPosition().phi();
0781 the_map["tsos_IP_pt"] = tsos_IP.globalMomentum().perp();
0782 the_map["tsos_IP_pt_eta"] = tsos_IP.globalMomentum().eta();
0783 the_map["tsos_IP_pt_phi"] = tsos_IP.globalMomentum().phi();
0784 const AlgebraicSymMatrix55& matrix_IP = tsos_IP.curvilinearError().matrix();
0785 the_map["err0_IP"] = sqrt(matrix_IP[0][0]);
0786 the_map["err1_IP"] = sqrt(matrix_IP[1][1]);
0787 the_map["err2_IP"] = sqrt(matrix_IP[2][2]);
0788 the_map["err3_IP"] = sqrt(matrix_IP[3][3]);
0789 the_map["err4_IP"] = sqrt(matrix_IP[4][4]);
0790 the_map["tsos_IP_valid"] = 1.0;
0791 } else {
0792 the_map["tsos_IP_eta"] = -999;
0793 the_map["tsos_IP_phi"] = -999;
0794 the_map["tsos_IP_pt"] = -999;
0795 the_map["tsos_IP_pt_eta"] = -999;
0796 the_map["tsos_IP_pt_phi"] = -999;
0797 the_map["err0_IP"] = -999;
0798 the_map["err1_IP"] = -999;
0799 the_map["err2_IP"] = -999;
0800 the_map["err3_IP"] = -999;
0801 the_map["err4_IP"] = -999;
0802 the_map["tsos_IP_valid"] = 0.0;
0803 }
0804 if (tsos_MuS.isValid()) {
0805 the_map["tsos_MuS_eta"] = tsos_MuS.globalPosition().eta();
0806 the_map["tsos_MuS_phi"] = tsos_MuS.globalPosition().phi();
0807 the_map["tsos_MuS_pt"] = tsos_MuS.globalMomentum().perp();
0808 the_map["tsos_MuS_pt_eta"] = tsos_MuS.globalMomentum().eta();
0809 the_map["tsos_MuS_pt_phi"] = tsos_MuS.globalMomentum().phi();
0810 const AlgebraicSymMatrix55& matrix_MuS = tsos_MuS.curvilinearError().matrix();
0811 the_map["err0_MuS"] = sqrt(matrix_MuS[0][0]);
0812 the_map["err1_MuS"] = sqrt(matrix_MuS[1][1]);
0813 the_map["err2_MuS"] = sqrt(matrix_MuS[2][2]);
0814 the_map["err3_MuS"] = sqrt(matrix_MuS[3][3]);
0815 the_map["err4_MuS"] = sqrt(matrix_MuS[4][4]);
0816 the_map["tsos_MuS_valid"] = 1.0;
0817 } else {
0818 the_map["tsos_MuS_eta"] = -999;
0819 the_map["tsos_MuS_phi"] = -999;
0820 the_map["tsos_MuS_pt"] = -999;
0821 the_map["tsos_MuS_pt_eta"] = -999;
0822 the_map["tsos_MuS_pt_phi"] = -999;
0823 the_map["err0_MuS"] = -999;
0824 the_map["err1_MuS"] = -999;
0825 the_map["err2_MuS"] = -999;
0826 the_map["err3_MuS"] = -999;
0827 the_map["err4_MuS"] = -999;
0828 the_map["tsos_MuS_valid"] = 0.0;
0829 }
0830 }
0831
0832
0833
0834
0835 void TSGForOIDNN::evaluateClassifier(const std::unordered_map<std::string, float>& feature_map,
0836 tensorflow::Session* session,
0837 const pt::ptree& metadata,
0838 StrategyParameters& out,
0839 bool& dnnSuccess) const {
0840 int n_features = metadata.get<int>("n_features", 0);
0841
0842
0843 tensorflow::Tensor input(tensorflow::DT_FLOAT, {1, n_features});
0844 std::string fname;
0845 int i_feature = 0;
0846 for (const pt::ptree::value_type& feature : metadata.get_child("feature_names")) {
0847 fname = feature.second.data();
0848 if (feature_map.find(fname) == feature_map.end()) {
0849
0850 dnnSuccess = false;
0851 } else {
0852 input.matrix<float>()(0, i_feature) = float(feature_map.at(fname));
0853 i_feature++;
0854 }
0855 }
0856
0857
0858 std::vector<tensorflow::Tensor> outputs;
0859
0860
0861 std::string inputLayer = metadata.get<std::string>("input_layer");
0862 std::string outputLayer = metadata.get<std::string>("output_layer");
0863
0864 tensorflow::run(session, {{inputLayer, input}}, {outputLayer}, &outputs);
0865 tensorflow::Tensor out_tensor = outputs[0];
0866 tensorflow::TTypes<float, 1>::Matrix dnn_outputs = out_tensor.matrix<float>();
0867
0868
0869 int imax = -1;
0870 float out_max = 0;
0871 for (long long int i = 0; i < out_tensor.dim_size(1); i++) {
0872 float ith_output = dnn_outputs(0, i);
0873 if (ith_output > out_max) {
0874 imax = i;
0875 out_max = ith_output;
0876 }
0877 }
0878
0879
0880 const std::string label = "output_labels.label_" + std::to_string(imax);
0881 out.nHBd = metadata.get<int>(label + ".nHBd");
0882 out.nHLIP = metadata.get<int>(label + ".nHLIP");
0883 out.nHLMuS = metadata.get<int>(label + ".nHLMuS");
0884 out.sf = metadata.get<int>(label + ".SF");
0885
0886 dnnSuccess = true;
0887 }
0888
0889
0890
0891
0892 void TSGForOIDNN::evaluateRegressor(const std::unordered_map<std::string, float>& feature_map,
0893 tensorflow::Session* session_HB,
0894 const pt::ptree& metadata_HB,
0895 tensorflow::Session* session_HLIP,
0896 const pt::ptree& metadata_HLIP,
0897 tensorflow::Session* session_HLMuS,
0898 const pt::ptree& metadata_HLMuS,
0899 StrategyParameters& out,
0900 bool& dnnSuccess) const {
0901 int n_features = metadata_HB.get<int>("n_features", 0);
0902
0903
0904 tensorflow::Tensor input(tensorflow::DT_FLOAT, {1, n_features});
0905 std::string fname;
0906 int i_feature = 0;
0907 for (const pt::ptree::value_type& feature : metadata_HB.get_child("feature_names")) {
0908 fname = feature.second.data();
0909 if (feature_map.find(fname) == feature_map.end()) {
0910
0911 dnnSuccess = false;
0912 } else {
0913 input.matrix<float>()(0, i_feature) = float(feature_map.at(fname));
0914 i_feature++;
0915 }
0916 }
0917
0918
0919 std::vector<tensorflow::Tensor> outputs_HB;
0920
0921 std::string inputLayer_HB = metadata_HB.get<std::string>("input_layer");
0922 std::string outputLayer_HB = metadata_HB.get<std::string>("output_layer");
0923 tensorflow::run(session_HB, {{inputLayer_HB, input}}, {outputLayer_HB}, &outputs_HB);
0924 tensorflow::Tensor out_tensor_HB = outputs_HB[0];
0925 tensorflow::TTypes<float, 1>::Matrix dnn_outputs_HB = out_tensor_HB.matrix<float>();
0926
0927
0928 std::vector<tensorflow::Tensor> outputs_HLIP;
0929
0930 std::string inputLayer_HLIP = metadata_HLIP.get<std::string>("input_layer");
0931 std::string outputLayer_HLIP = metadata_HLIP.get<std::string>("output_layer");
0932 tensorflow::run(session_HLIP, {{inputLayer_HLIP, input}}, {outputLayer_HLIP}, &outputs_HLIP);
0933 tensorflow::Tensor out_tensor_HLIP = outputs_HLIP[0];
0934 tensorflow::TTypes<float, 1>::Matrix dnn_outputs_HLIP = out_tensor_HLIP.matrix<float>();
0935
0936
0937 std::vector<tensorflow::Tensor> outputs_HLMuS;
0938
0939 std::string inputLayer_HLMuS = metadata_HLMuS.get<std::string>("input_layer");
0940 std::string outputLayer_HLMuS = metadata_HLMuS.get<std::string>("output_layer");
0941 tensorflow::run(session_HLMuS, {{inputLayer_HLMuS, input}}, {outputLayer_HLMuS}, &outputs_HLMuS);
0942 tensorflow::Tensor out_tensor_HLMuS = outputs_HLMuS[0];
0943 tensorflow::TTypes<float, 1>::Matrix dnn_outputs_HLMuS = out_tensor_HLMuS.matrix<float>();
0944
0945
0946 out.nHBd = round(dnn_outputs_HB(0, 0));
0947 out.nHLIP = round(dnn_outputs_HLIP(0, 0));
0948 out.sf = round(dnn_outputs_HLIP(0, 1));
0949 out.nHLMuS = round(dnn_outputs_HLMuS(0, 0));
0950
0951
0952 out.nHBd = std::clamp(out.nHBd, 0, 10);
0953 out.nHLIP = std::clamp(out.nHLIP, 0, 10);
0954 out.nHLMuS = std::clamp(out.nHLMuS, 0, 10);
0955
0956
0957 if (out.nHBd == 0 && out.nHLIP == 0 && out.nHLMuS == 0) {
0958
0959 out.nHBd = 1;
0960 out.nHLIP = 5;
0961 }
0962
0963
0964
0965 if (out.sf <= 0)
0966 out.sf = 2;
0967 if (out.sf > 10)
0968 out.sf = 10;
0969
0970 dnnSuccess = true;
0971 }
0972
0973
0974
0975
0976 void TSGForOIDNN::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0977 edm::ParameterSetDescription desc;
0978 desc.add<edm::InputTag>("src", edm::InputTag("hltL2Muons", "UpdatedAtVtx"));
0979 desc.add<int>("layersToTry", 2);
0980 desc.add<double>("fixedErrorRescaleFactorForHitless", 2.0);
0981 desc.add<int>("hitsToTry", 1);
0982 desc.add<edm::InputTag>("MeasurementTrackerEvent", edm::InputTag("hltSiStripClusters"));
0983 desc.add<std::string>("estimator", "hltESPChi2MeasurementEstimator100");
0984 desc.add<double>("maxEtaForTOB", 1.8);
0985 desc.add<double>("minEtaForTEC", 0.7);
0986 desc.addUntracked<bool>("debug", false);
0987 desc.add<unsigned int>("maxSeeds", 20);
0988 desc.add<unsigned int>("maxHitlessSeeds", 5);
0989 desc.add<unsigned int>("maxHitSeeds", 1);
0990 desc.add<std::string>("propagatorName", "PropagatorWithMaterialParabolicMf");
0991 desc.add<unsigned int>("maxHitlessSeedsIP", 5);
0992 desc.add<unsigned int>("maxHitlessSeedsMuS", 0);
0993 desc.add<unsigned int>("maxHitDoubletSeeds", 0);
0994 desc.add<bool>("getStrategyFromDNN", false);
0995 desc.add<bool>("useRegressor", false);
0996 desc.add<std::string>("dnnMetadataPath", "");
0997 descriptions.add("tsgForOIDNN", desc);
0998 }
0999
1000 DEFINE_FWK_MODULE(TSGForOIDNN);