Merge pull request #45283 from AuroraPerego/removeTrajectories

TICLv5: do not use trajectories to compute the track path length

Author: cmsbuild

RecoHGCal/TICL/plugins/TICLCandidateProducer.cc

@@ -15,7 +15,6 @@
 #include "FWCore/Framework/interface/ConsumesCollector.h"
 #include "DataFormats/Common/interface/OrphanHandle.h"
 
-#include "DataFormats/BeamSpot/interface/BeamSpot.h"
 #include "DataFormats/CaloRecHit/interface/CaloCluster.h"
 #include "DataFormats/HGCalReco/interface/Common.h"
 #include "DataFormats/HGCalReco/interface/MtdHostCollection.h"
@@ -27,6 +26,8 @@
 #include "DataFormats/HGCalReco/interface/TICLCandidate.h"
 #include "DataFormats/TrackReco/interface/TrackFwd.h"
 #include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"
+#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
+#include "DataFormats/GeometryVector/interface/GlobalVector.h"
 
 #include "RecoHGCal/TICL/interface/TICLInterpretationAlgoBase.h"
 #include "RecoHGCal/TICL/plugins/TICLInterpretationPluginFactory.h"
@@ -41,9 +42,7 @@
 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateClosestToBeamLine.h"
-#include "TrackingTools/PatternTools/interface/Trajectory.h"
-#include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
-#include "TrackingTools/PatternTools/interface/TSCBLBuilderWithPropagator.h"
+#include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
 
 #include "MagneticField/Engine/interface/MagneticField.h"
 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
@@ -71,8 +70,7 @@ class TICLCandidateProducer : public edm::stream::EDProducer<> {
   template <typename F>
   void assignTimeToCandidates(std::vector<TICLCandidate> &resultCandidates,
                               edm::Handle<std::vector<reco::Track>> track_h,
-                              edm::Handle<MtdHostCollection> inputTiming_h,
-                              TrajTrackAssociationCollection trjtrks,
+                              MtdHostCollection::ConstView &inputTimingView,
                               F func) const;
 
   std::unique_ptr<TICLInterpretationAlgoBase<reco::Track>> generalInterpretationAlgo_;
@@ -88,7 +86,6 @@ class TICLCandidateProducer : public edm::stream::EDProducer<> {
   std::vector<edm::EDGetTokenT<std::vector<float>>> original_masks_tokens_;
 
   const edm::EDGetTokenT<std::vector<reco::Track>> tracks_token_;
-  const edm::EDGetTokenT<TrajTrackAssociationCollection> trajTrackAssToken_;
   edm::EDGetTokenT<MtdHostCollection> inputTimingToken_;
 
   const edm::EDGetTokenT<std::vector<reco::Muon>> muons_token_;
@@ -101,7 +98,7 @@ class TICLCandidateProducer : public edm::stream::EDProducer<> {
   const std::string detector_;
   const std::string propName_;
   const edm::ESGetToken<Propagator, TrackingComponentsRecord> propagator_token_;
-  const edm::EDGetTokenT<reco::BeamSpot> bsToken_;
+  const edm::ESGetToken<GlobalTrackingGeometry, GlobalTrackingGeometryRecord> trackingGeometry_token_;
 
   std::once_flag initializeGeometry_;
   const HGCalDDDConstants *hgcons_;
@@ -111,15 +108,16 @@ class TICLCandidateProducer : public edm::stream::EDProducer<> {
   edm::ESGetToken<HGCalDDDConstants, IdealGeometryRecord> hdc_token_;
   edm::ESHandle<MagneticField> bfield_;
   edm::ESHandle<Propagator> propagator_;
+  edm::ESHandle<GlobalTrackingGeometry> trackingGeometry_;
   static constexpr float c_light_ = CLHEP::c_light * CLHEP::ns / CLHEP::cm;
+  static constexpr float timeRes = 0.02f;
 };
 
 TICLCandidateProducer::TICLCandidateProducer(const edm::ParameterSet &ps)
     : clusters_token_(consumes<std::vector<reco::CaloCluster>>(ps.getParameter<edm::InputTag>("layer_clusters"))),
       clustersTime_token_(
           consumes<edm::ValueMap<std::pair<float, float>>>(ps.getParameter<edm::InputTag>("layer_clustersTime"))),
       tracks_token_(consumes<std::vector<reco::Track>>(ps.getParameter<edm::InputTag>("tracks"))),
-      trajTrackAssToken_(consumes<TrajTrackAssociationCollection>(ps.getParameter<edm::InputTag>("trjtrkAss"))),
       inputTimingToken_(consumes<MtdHostCollection>(ps.getParameter<edm::InputTag>("timingSoA"))),
       muons_token_(consumes<std::vector<reco::Muon>>(ps.getParameter<edm::InputTag>("muons"))),
       useMTDTiming_(ps.getParameter<bool>("useMTDTiming")),
@@ -131,7 +129,8 @@ TICLCandidateProducer::TICLCandidateProducer(const edm::ParameterSet &ps)
       propName_(ps.getParameter<std::string>("propagator")),
       propagator_token_(
           esConsumes<Propagator, TrackingComponentsRecord, edm::Transition::BeginRun>(edm::ESInputTag("", propName_))),
-      bsToken_(consumes<reco::BeamSpot>(ps.getParameter<edm::InputTag>("beamspot"))),
+      trackingGeometry_token_(
+          esConsumes<GlobalTrackingGeometry, GlobalTrackingGeometryRecord, edm::Transition::BeginRun>()),
       cutTk_(ps.getParameter<std::string>("cutTk")) {
   // These are the CLUE3DEM Tracksters put in the event by the TracksterLinksProducer with the superclustering algorithm
   for (auto const &tag : ps.getParameter<std::vector<edm::InputTag>>("egamma_tracksters_collections")) {
@@ -192,6 +191,8 @@ void TICLCandidateProducer::beginRun(edm::Run const &iEvent, edm::EventSetup con
   bfield_ = es.getHandle(bfield_token_);
   propagator_ = es.getHandle(propagator_token_);
   generalInterpretationAlgo_->initialize(hgcons_, rhtools_, bfield_, propagator_);
+
+  trackingGeometry_ = es.getHandle(trackingGeometry_token_);
 };
 
 void filterTracks(edm::Handle<std::vector<reco::Track>> tkH,
@@ -238,15 +239,13 @@ void TICLCandidateProducer::produce(edm::Event &evt, const edm::EventSetup &es)
   evt.getByToken(tracks_token_, tracks_h);
   const auto &tracks = *tracks_h;
 
-  const auto &trjtrks = evt.get(trajTrackAssToken_);
-
   edm::Handle<MtdHostCollection> inputTiming_h;
+  MtdHostCollection::ConstView inputTimingView;
   if (useMTDTiming_) {
     evt.getByToken(inputTimingToken_, inputTiming_h);
+    inputTimingView = (*inputTiming_h).const_view();
   }
 
-  const auto &bs = evt.get(bsToken_);
-
   auto const bFieldProd = bfield_.product();
   const Propagator *propagator = propagator_.product();
 
@@ -348,68 +347,55 @@ void TICLCandidateProducer::produce(edm::Event &evt, const edm::EventSetup &es)
     }
   }
 
-  auto getPathLength =
-      [&](const reco::Track track, float zVal, const Trajectory &traj, TrajectoryStateClosestToBeamLine &tscbl) {
-        TrajectoryStateOnSurface stateForProjectionToBeamLineOnSurface =
-            traj.closestMeasurement(GlobalPoint(bs.x0(), bs.y0(), bs.z0())).updatedState();
+  auto getPathLength = [&](const reco::Track &track, float zVal) {
+    const auto &fts_inn = trajectoryStateTransform::innerFreeState(track, bFieldProd);
+    const auto &fts_out = trajectoryStateTransform::outerFreeState(track, bFieldProd);
+    const auto &surf_inn = trajectoryStateTransform::innerStateOnSurface(track, *trackingGeometry_, bFieldProd);
+    const auto &surf_out = trajectoryStateTransform::outerStateOnSurface(track, *trackingGeometry_, bFieldProd);
 
-        if (!stateForProjectionToBeamLineOnSurface.isValid()) {
-          edm::LogError("CannotPropagateToBeamLine")
-              << "the state on the closest measurement is not valid. skipping track.";
-          return 0.f;
-        }
-        const FreeTrajectoryState &stateForProjectionToBeamLine = *stateForProjectionToBeamLineOnSurface.freeState();
-
-        TSCBLBuilderWithPropagator tscblBuilder(*propagator);
-        tscbl = tscblBuilder(stateForProjectionToBeamLine, bs);
-
-        if (tscbl.isValid()) {
-          auto const &tscblPCA = tscbl.trackStateAtPCA();
-          auto const &innSurface = traj.direction() == alongMomentum ? traj.firstMeasurement().updatedState().surface()
-                                                                     : traj.lastMeasurement().updatedState().surface();
-          auto const &extSurface = traj.direction() == alongMomentum ? traj.lastMeasurement().updatedState().surface()
-                                                                     : traj.firstMeasurement().updatedState().surface();
-          float pathlength = propagator->propagateWithPath(tscblPCA, innSurface).second;
-
-          if (pathlength) {
-            const auto &fts_inn = trajectoryStateTransform::innerFreeState((track), bFieldProd);
-            const auto &t_inn_out = propagator->propagateWithPath(fts_inn, extSurface);
-
-            if (t_inn_out.first.isValid()) {
-              pathlength += t_inn_out.second;
-
-              std::pair<float, float> rMinMax = hgcons_->rangeR(zVal, true);
-
-              int iSide = int(track.eta() > 0);
-              float zSide = (iSide == 0) ? (-1. * zVal) : zVal;
-              const auto &disk = std::make_unique<GeomDet>(
-                  Disk::build(Disk::PositionType(0, 0, zSide),
-                              Disk::RotationType(),
-                              SimpleDiskBounds(rMinMax.first, rMinMax.second, zSide - 0.5, zSide + 0.5))
-                      .get());
-              const auto &fts_out = trajectoryStateTransform::outerFreeState((track), bFieldProd);
-              const auto &tsos = propagator->propagateWithPath(fts_out, disk->surface());
-
-              if (tsos.first.isValid()) {
-                pathlength += tsos.second;
-                return pathlength;
-              }
-            }
-          }
+    Basic3DVector<float> pos(track.referencePoint());
+    Basic3DVector<float> mom(track.momentum());
+    FreeTrajectoryState stateAtBeamspot{GlobalPoint(pos), GlobalVector(mom), track.charge(), bFieldProd};
+
+    float pathlength = propagator->propagateWithPath(stateAtBeamspot, surf_inn.surface()).second;
+
+    if (pathlength) {
+      const auto &t_inn_out = propagator->propagateWithPath(fts_inn, surf_out.surface());
+
+      if (t_inn_out.first.isValid()) {
+        pathlength += t_inn_out.second;
+
+        std::pair<float, float> rMinMax = hgcons_->rangeR(zVal, true);
+
+        int iSide = int(track.eta() > 0);
+        float zSide = (iSide == 0) ? (-1. * zVal) : zVal;
+        const auto &disk = std::make_unique<GeomDet>(
+            Disk::build(Disk::PositionType(0, 0, zSide),
+                        Disk::RotationType(),
+                        SimpleDiskBounds(rMinMax.first, rMinMax.second, zSide - 0.5, zSide + 0.5))
+                .get());
+        const auto &tsos = propagator->propagateWithPath(fts_out, disk->surface());
+
+        if (tsos.first.isValid()) {
+          pathlength += tsos.second;
+          return pathlength;
         }
-        return 0.f;
-      };
+      }
+    }
+    edm::LogWarning("TICLCandidateProducer")
+        << "Not able to use the track to compute the path length. A straight line will be used instead.";
+    return 0.f;
+  };
 
-  assignTimeToCandidates(*resultCandidates, tracks_h, inputTiming_h, trjtrks, getPathLength);
+  assignTimeToCandidates(*resultCandidates, tracks_h, inputTimingView, getPathLength);
 
   evt.put(std::move(resultCandidates));
 }
 
 template <typename F>
 void TICLCandidateProducer::assignTimeToCandidates(std::vector<TICLCandidate> &resultCandidates,
                                                    edm::Handle<std::vector<reco::Track>> track_h,
-                                                   edm::Handle<MtdHostCollection> inputTiming_h,
-                                                   TrajTrackAssociationCollection trjtrks,
+                                                   MtdHostCollection::ConstView &inputTimingView,
                                                    F func) const {
   for (auto &cand : resultCandidates) {
     float beta = 1;
@@ -426,30 +412,18 @@ void TICLCandidateProducer::assignTimeToCandidates(std::vector<TICLCandidate> &r
         auto path = std::sqrt(x * x + y * y + z * z);
         if (cand.trackPtr().get() != nullptr) {
           const auto &trackIndex = cand.trackPtr().get() - (edm::Ptr<reco::Track>(track_h, 0)).get();
-          if (useMTDTiming_) {
-            auto const &inputTimingView = (*inputTiming_h).const_view();
-            if (inputTimingView.timeErr()[trackIndex] > 0) {
-              const auto xMtd = inputTimingView.posInMTD_x()[trackIndex];
-              const auto yMtd = inputTimingView.posInMTD_y()[trackIndex];
-              const auto zMtd = inputTimingView.posInMTD_z()[trackIndex];
-
-              beta = inputTimingView.beta()[trackIndex];
-              path = std::sqrt((x - xMtd) * (x - xMtd) + (y - yMtd) * (y - yMtd) + (z - zMtd) * (z - zMtd)) +
-                     inputTimingView.pathLength()[trackIndex];
-            }
+          if (useMTDTiming_ and inputTimingView.timeErr()[trackIndex] > 0) {
+            const auto xMtd = inputTimingView.posInMTD_x()[trackIndex];
+            const auto yMtd = inputTimingView.posInMTD_y()[trackIndex];
+            const auto zMtd = inputTimingView.posInMTD_z()[trackIndex];
+
+            beta = inputTimingView.beta()[trackIndex];
+            path = std::sqrt((x - xMtd) * (x - xMtd) + (y - yMtd) * (y - yMtd) + (z - zMtd) * (z - zMtd)) +
+                   inputTimingView.pathLength()[trackIndex];
           } else {
-            const auto &trackIndex = cand.trackPtr().get() - (edm::Ptr<reco::Track>(track_h, 0)).get();
-            for (const auto &trj : trjtrks) {
-              if (trj.val != edm::Ref<std::vector<reco::Track>>(track_h, trackIndex))
-                continue;
-              const Trajectory &traj = *trj.key;
-              TrajectoryStateClosestToBeamLine tscbl;
-
-              float pathLength = func(*(cand.trackPtr().get()), z, traj, tscbl);
-              if (pathLength) {
-                path = pathLength;
-                break;
-              }
+            float pathLength = func(*(cand.trackPtr().get()), z);
+            if (pathLength) {
+              path = pathLength;
             }
           }
         }
@@ -460,13 +434,14 @@ void TICLCandidateProducer::assignTimeToCandidates(std::vector<TICLCandidate> &r
     if (invTimeErr > 0) {
       time = time / invTimeErr;
       // FIXME_ set a liminf of 0.02 ns on the ts error (based on residuals)
-      timeErr = sqrt(1.f / invTimeErr) > 0.02 ? sqrt(1.f / invTimeErr) : 0.02;
+      timeErr = sqrt(1.f / invTimeErr);
+      if (timeErr < timeRes)
+        timeErr = timeRes;
       cand.setTime(time, timeErr);
     }
 
     if (useMTDTiming_ and cand.charge()) {
       // Check MTD timing availability
-      auto const &inputTimingView = (*inputTiming_h).const_view();
       const auto &trackIndex = cand.trackPtr().get() - (edm::Ptr<reco::Track>(track_h, 0)).get();
       const bool assocQuality = inputTimingView.MVAquality()[trackIndex] > timingQualityThreshold_;
       if (assocQuality) {
@@ -507,12 +482,10 @@ void TICLCandidateProducer::fillDescriptions(edm::ConfigurationDescriptions &des
   desc.add<edm::InputTag>("layer_clusters", edm::InputTag("hgcalMergeLayerClusters"));
   desc.add<edm::InputTag>("layer_clustersTime", edm::InputTag("hgcalMergeLayerClusters", "timeLayerCluster"));
   desc.add<edm::InputTag>("tracks", edm::InputTag("generalTracks"));
-  desc.add<edm::InputTag>("trjtrkAss", edm::InputTag("generalTracks"));
   desc.add<edm::InputTag>("timingSoA", edm::InputTag("mtdSoA"));
   desc.add<edm::InputTag>("muons", edm::InputTag("muons1stStep"));
   desc.add<std::string>("detector", "HGCAL");
   desc.add<std::string>("propagator", "PropagatorWithMaterial");
-  desc.add<edm::InputTag>("beamspot", edm::InputTag("offlineBeamSpot"));
   desc.add<bool>("useMTDTiming", true);
   desc.add<bool>("useTimingAverage", true);
   desc.add<double>("timingQualityThreshold", 0.5);