Cut on pt rather x,y components and apply pt cut also on next tracks of a multi track trajectory.

InnerDetector/InDetRecTools/SiCombinatorialTrackFinderTool_xk/SiCombinatorialTrackFinderTool_xk/SiCombinatorialTrackFinder_xk.h

@@ -186,6 +186,8 @@ namespace InDet {
       BooleanProperty m_doFastTracking{this, "doFastTracking", false};
       StringProperty m_fieldmode{this, "MagneticFieldMode", "MapSolenoid", "Mode of magnetic field"};
       DoubleProperty m_qualityCut{this, "TrackQualityCut", 9.3, "Simple track quality cut"};
+      FloatProperty m_minPtCut{this, "MinFinalPtCut", 100, "Cut on the pt of the final track. Must be >0 to avoid NANs when computing eta."};
+      float  m_minPt2Cut=0;
       BooleanProperty m_writeHolesFromPattern{this, "writeHolesFromPattern", false,"Flag to activate writing hole info from the pattern recognition"}; 
       //@}
 
@@ -223,7 +225,7 @@ namespace InDet {
       static void getTrackQualityCuts(SiCombinatorialTrackFinderData_xk& data, const TrackQualityCuts&) ;
 
       Trk::Track* convertToTrack(SiCombinatorialTrackFinderData_xk& data) const;
-      static Trk::Track* convertToNextTrack(SiCombinatorialTrackFinderData_xk& data) ;
+      Trk::Track* convertToNextTrack(SiCombinatorialTrackFinderData_xk& data) const;
  
       void magneticFieldInit();
 

InnerDetector/InDetRecTools/SiCombinatorialTrackFinderTool_xk/src/SiCombinatorialTrackFinder_xk.cxx

@@ -27,6 +27,11 @@
 #include <utility>
 #include <stdexcept>
 
+namespace {
+   template <typename T>
+   inline T sqr(T a) {return a*a; }
+}
+
 ///////////////////////////////////////////////////////////////////
 // Constructor
 ///////////////////////////////////////////////////////////////////
@@ -104,6 +109,8 @@ StatusCode InDet::SiCombinatorialTrackFinder_xk::initialize()
   ATH_CHECK( m_fieldCondObjInputKey.initialize() );
   ATH_CHECK( m_pixelDetElStatus.initialize( !m_pixelDetElStatus.empty() && m_usePIX) );
   ATH_CHECK( m_sctDetElStatus.initialize( !m_sctDetElStatus.empty() && m_useSCT) );
+
+  m_minPt2Cut = sqr(m_minPtCut.value());
   return StatusCode::SUCCESS;
 }
 
@@ -418,7 +425,7 @@ const std::list<Trk::Track*>& InDet::SiCombinatorialTrackFinder_xk::getTracks
   // Trk::Track production
   //
   Trk::Track* t = convertToTrack(data);
-  if (t==nullptr) return data.tracks();
+  if (t==nullptr) return data.tracks(); // @TODO should one check if convertToNextTrack would yield anything ?
   if (m_writeHolesFromPattern) data.addPatternHoleSearchOutcome(t,data.trajectory().getHoleSearchResult());
 
   ++data.findtracks();
@@ -489,7 +496,7 @@ const std::list<Trk::Track*>&  InDet::SiCombinatorialTrackFinder_xk::getTracksWi
     data.trackinfo() = oldinfo;
     data.tools().setMultiTracks(mult,Xi2m);
 
-    if (!t) return data.tracks();
+    if (t==nullptr) return data.tracks(); // @TODO should one check whether the next findTrack call would yield something ?
     if (m_writeHolesFromPattern) data.addPatternHoleSearchOutcome(t,data.trajectory().getHoleSearchResult());
     ++data.findtracks();
     data.tracks().push_back(t);
@@ -801,17 +808,17 @@ InDet::SiCombinatorialTrackFinder_xk::EStat_t InDet::SiCombinatorialTrackFinder_
 
 Trk::Track* InDet::SiCombinatorialTrackFinder_xk::convertToTrack(SiCombinatorialTrackFinderData_xk& data) const
 {
-   Trk::TrackParameters* param = data.trajectory().firstTrackParameters();
-   if (param) {
-      auto momentum = param->momentum();
-      const double  p_limit = momentum.z() * 1e-6;
-      // reject tracks too parallel to z-axis
-      if (std::abs(momentum.x()) < p_limit && std::abs(momentum.y())< p_limit) {
-         ATH_MSG_WARNING( "Reject track with large |eta| : m= " <<  (std::sqrt(momentum.x() * momentum.x() + momentum.y() * momentum.y() + momentum.z()*momentum.z()))
-                          << " z= "  << momentum.z() << " => eta = 1/2 * log( (m+z) / (m-z) )" );
-         return nullptr;
-      }
-   }
+  Trk::TrackParameters* param = data.trajectory().firstTrackParameters();
+  if (param) {
+     auto momentum = param->momentum();
+     const auto  pt2 = momentum.perp2();
+     // reject tracks with small pT
+     // The cut should be large enough otherwise eta computation of such tracks may yield NANs.
+     if (pt2 < m_minPt2Cut) {
+        ATH_MSG_WARNING( "Reject low pT track (pT = " << sqrt(pt2) << " < " << m_minPtCut.value() << ")");
+        return nullptr;
+     }
+  }
   if (!data.simpleTrack()) {
     return new Trk::Track(data.trackinfo(),
                           data.trajectory().convertToTrackStateOnSurface(data.cosmicTrack()),
@@ -836,10 +843,30 @@ Trk::Track* InDet::SiCombinatorialTrackFinder_xk::convertToTrack(SiCombinatorial
 // Next Trk::Track production
 ///////////////////////////////////////////////////////////////////
 
-Trk::Track* InDet::SiCombinatorialTrackFinder_xk::convertToNextTrack(SiCombinatorialTrackFinderData_xk& data) 
+Trk::Track* InDet::SiCombinatorialTrackFinder_xk::convertToNextTrack(SiCombinatorialTrackFinderData_xk& data) const
 {
   auto tsos = data.trajectory().convertToNextTrackStateOnSurface();
   if (tsos.empty()) return nullptr;
+
+  // verify first track parameters
+  const Trk::TrackParameters *param = nullptr;
+  for (const Trk::TrackStateOnSurface *a_tsos : tsos) {
+     const Trk::TrackParameters *param = a_tsos->trackParameters();
+     if (param) {
+        break;
+     }
+  }
+
+  if (param) {
+     auto momentum = param->momentum();
+     const auto  pt2 = momentum.perp2();
+     // reject tracks with small pT
+     // The cut should be large enough otherwise eta computation of such tracks may yield NANs.
+     if (pt2 < m_minPt2Cut) {
+        ATH_MSG_WARNING( "Reject low pT track (pT = " << sqrt(pt2) << " < " << m_minPtCut.value() << ")");
+        return nullptr;
+     }
+  }
   return new Trk::Track(data.trackinfo(),
                         std::move(tsos),
                         data.trajectory().convertToFitQuality());