Stability controls for the bfcorr BField domain walk

Fit/Config.cc

@@ -8,6 +8,8 @@ namespace KinKal {
       << " diverge dpar chisq " << kkconfig.pdchisq_
       << " diverge traj gap (mm) " << kkconfig.divgap_
       << " fractional momentum tolerance " << kkconfig.tol_
+      << " min domain step (ns) " << kkconfig.mindtstep_
+      << " domain margin (ns) " << kkconfig.domainmargin_
       << " min NDOF " << kkconfig.minndof_
       << " BField correction " << kkconfig.bfcorr_
       << " with " << kkconfig.schedule().size()

Fit/Config.hh

@@ -11,6 +11,7 @@
 #include <vector>
 #include <memory>
 #include <algorithm>
+#include <limits>
 #include <ostream>
 #include <istream>
 
@@ -29,6 +30,8 @@ namespace KinKal {
     double pdchisq_ = 1.0e6; // maximum allowed parameter change (units of chisqred) WRT previous reference
     double divgap_ = 1.0e2; // maximum average gap of trajectory before calling it diverged (mm)
     double tol_ = 1.0e-4; // tolerance on fractional momentum accuracy due to BField domain steps
+    double mindtstep_ = 0.0; // ns: minimum domain range. >0 bounds the domain count at (walk window)/mindtstep_; 0 leaves it unbounded
+    double domainmargin_ = std::numeric_limits<double>::max(); // ns: max time a domain may extend beyond the active range; finite confines the walk and its field sampling to range +/- margin
     unsigned minndof_ = 5; // minimum number of DOFs to continue fit
     bool bfcorr_ = true; // whether to make BFieldMap corrections in the fit
     bool ends_ = true; // process the passive effects at each end of the track after schedule completion

Fit/Status.cc

@@ -21,6 +21,8 @@ namespace KinKal {
         return "OutsideBFieldMap ";
       case Status::failed:
         return "Failed ";
+      case Status::incompatiblepiece:
+        return "IncompatiblePiece ";
     }
   }
 

Fit/Status.hh

@@ -9,7 +9,7 @@
 namespace KinKal {
 // struct to define fit status
   struct Status {
-    enum status {unfit=-1,converged,unconverged,lowNDOF,gapdiverged,paramsdiverged,chisqdiverged,outsidemap,failed}; // fit status
+    enum status {unfit=-1,converged,unconverged,lowNDOF,gapdiverged,paramsdiverged,chisqdiverged,outsidemap,failed,incompatiblepiece}; // fit status
     unsigned miter_; // meta-iteration number;
     unsigned iter_; // iteration number;
     status status_; // current status

General/BFieldMap.hh

@@ -11,6 +11,7 @@
 #include <algorithm>
 #include <cstdarg>
 #include <cmath>
+#include <optional>
 #include <ostream>
 
 namespace KinKal {
@@ -27,17 +28,41 @@ namespace KinKal {
       virtual bool inRange(VEC3 const& position) const = 0;
       virtual ~BFieldMap(){}
       virtual void print(std::ostream& os ) const = 0;
-      BFieldMap(){}
+      // smallest |B| (T) usable for field-corrected transport; 0 (default) disables low-field
+      // protection, leaving an unusable sample a hard failure for the caller
+      BFieldMap(double minfield=0.0) : minfield_(minfield) {}
+      double minField() const { return minfield_; }
+      bool protecting() const { return minfield_ > 0.0; } // is low-field protection enabled?
+      // The field here, if this position can support field-corrected transport; nothing if it cannot.
+      std::optional<VEC3> usableField(VEC3 const& position) const {
+        if(!inRange(position)) return std::nullopt;   // outside: fieldDeriv is undefined here
+        VEC3 bf = fieldVect(position);
+        if(protecting() && bf.R() < std::max(minfield_,zeroField())) return std::nullopt;
+        return bf;
+      }
+      // predicate form, for the callers that only decide and never use the field
+      bool usable(VEC3 const& position) const { return usableField(position).has_value(); }
       // disallow copy and equivalence
       BFieldMap(BFieldMap const& ) = delete;
       BFieldMap& operator =(BFieldMap const& ) = delete;
       // speed of light in units to convert Tesla to mm (bending radius)
       static double constexpr cbar() { return CLHEP::c_light/1000.0; }
+      // |B| below this (T) is treated as physically zero when deciding to hand off extrapolation;
+      // fit paths must not sample here
+      static double constexpr zeroField() { return 1.0e-6; }
+      static bool isZeroField(VEC3 const& bvec) { return bvec.R() < zeroField(); }
       // templated interface for interacting with kinematic trajectory classes
       // how far can you go along the given kinematic trajectory till BField inhomogeneity makes the momentum accuracy out of (fractional) tolerance
       template<class KTRAJ> double rangeInTolerance(KTRAJ const& ktraj, double tstart, double tol) const;
+      // the domain step at tstart: rangeInTolerance with a floor applied. Callers must have confirmed
+      // the position is in range first -- rangeInTolerance samples fieldDeriv, undefined outside it.
+      template<class KTRAJ> double domainStep(KTRAJ const& ktraj, double tstart, double tol, double mindtstep) const {
+        return std::max(rangeInTolerance(ktraj,tstart,tol),mindtstep);
+      }
       // integrate the residual magentic force over the given kinematic trajectory and range due to the difference between the true field and the nominal field in the
       template<class KTRAJ> VEC3 integrate(KTRAJ const& ktraj, TimeRange const& trange) const;
+    private:
+      double minfield_; // smallest usable |B| (T); 0 disables low-field protection
   };
 
   template<class KTRAJ> VEC3 BFieldMap::integrate(KTRAJ const& ktraj, TimeRange const& trange) const {

Trajectory/CentralHelix.cc

@@ -32,6 +32,12 @@ namespace KinKal {
   CentralHelix::CentralHelix(VEC4 const &pos0, MOM4 const &mom0, int charge, VEC3 const &bnom,
       TimeRange const &trange) : trange_(trange), mass_(mom0.M()), bnom_(bnom)
   {
+    // A null nominal field has no valid parameterization here: radToMom below is 0, so omega comes out
+    // as signed zero and momentum() as 0/0 = NaN, while charge() -- which reads the sign of omega --
+    // becomes meaningless. Refuse rather than return a silently degenerate object. Callers that can meet
+    // a field-free region test BFieldMap::usableField() first; this is the backstop behind them.
+    if(bnom_.R() < BFieldMap::zeroField())
+      throw std::invalid_argument("CentralHelix::CentralHelix; null BNom");
     // Transform into the system where Z is along the Bfield.  This is a pure rotation about the origin
     VEC4 pos(pos0);
     MOM4 mom(mom0);
@@ -85,6 +91,9 @@ namespace KinKal {
   }
 
   void CentralHelix::resetBNom(VEC3 const& bnom) {
+    // same degeneracy as construction: refuse to move an existing piece onto a null field
+    if(bnom.R() < BFieldMap::zeroField())
+      throw std::invalid_argument("CentralHelix::resetBNom; null BNom");
     bnom_ = bnom;
     setTransforms();
   }

Trajectory/ParticleTrajectory.hh

@@ -18,15 +18,19 @@ namespace KinKal {
       // construct from an initial piece, which also provides kinematic information
       ParticleTrajectory(KTRAJ const& piece) : PTTRAJ(piece) {}
       ParticleTrajectory() : PTTRAJ() {}
+      // does this piece describe the same particle (mass and charge) as the existing ones? Callers that
+      // expect incompatibility as a routine outcome test this instead of catching append's throw
+      bool compatible(KTRAJ const& newpiece) const {
+        return PTTRAJ::pieces().size() == 0 ||
+          (fabs(newpiece.mass()-mass()) <= 1e-6 && newpiece.charge() == charge());
+      }
       //  append and prepend to check mass and charge consistency
       void append(KTRAJ const& newpiece, bool allowremove=false)  {
-        if(PTTRAJ::pieces().size() > 0){
-          if(fabs(newpiece.mass()-mass())>1e-6 || newpiece.charge() != charge()) throw std::invalid_argument("Invalid particle parameters");
-        }
+        if(!compatible(newpiece)) throw std::invalid_argument("Invalid particle parameters");
         PTTRAJ::append(newpiece,allowremove);
       }
       void prepend(KTRAJ const& newpiece, bool allowremove=false)  {
-        if(fabs(newpiece.mass()-mass())>1e-6 || newpiece.charge() != charge()) throw std::invalid_argument("Invalid particle parameters");
+        if(!compatible(newpiece)) throw std::invalid_argument("Invalid particle parameters");
         PTTRAJ::prepend(newpiece,allowremove);
       }
       // kinematic interface