KinKal stability knobs to allow bfcorr=true for CentralHelix extrap to CRV

BFieldGeom/inc/BFieldConfig.hh

@@ -64,9 +64,12 @@ namespace mu2e {
 
         bool flipBFieldMaps() const { return flipBFieldMaps_; }
 
+        // smallest |B| (T) KinKal may use for field-corrected transport; 0 disables its low-field protection
+        double minUsableField() const { return minUsableField_; }
+
        private:
         BFieldConfig()
-            : scaleFactor_(1.), writeBinaries_(false), verbosityLevel_(1), flipBFieldMaps_(false) {}
+            : scaleFactor_(1.), writeBinaries_(false), verbosityLevel_(1), flipBFieldMaps_(false), minUsableField_(0.) {}
 
         // G4BL, PARAM or possible future types.
         BFMapType mapType_;
@@ -91,6 +94,7 @@ namespace mu2e {
         bool writeBinaries_;
         int verbosityLevel_;
         bool flipBFieldMaps_;
+        double minUsableField_;
     };
 
 }  // namespace mu2e

GeometryService/src/BFieldConfigMaker.cc

@@ -47,6 +47,7 @@ namespace mu2e {
         bfconf_->flipBFieldMaps_ = config.getBool("bfield.flipMaps", false);
 
         bfconf_->scaleFactor_ = config.getDouble("bfield.scaleFactor", 1.0);
+        bfconf_->minUsableField_ = config.getDouble("bfield.minUsableField", 0.0);
 
         bfconf_->dsFieldForm_ = BFieldConfig::DSFieldModel(
             config.getInt("detSolFieldForm", BFieldConfig::dsModelUniform));

Mu2eG4/geom/bfgeom_v01_minfield.txt

@@ -0,0 +1,17 @@
+//
+// bfgeom_v01 with KinKal low-field protection enabled at the null-field threshold.
+//
+// bfield.minUsableField is the smallest |B| (T) KinKal uses for field-corrected transport. Below it, and
+// outside the maps, bfcorr extrapolation hands off to geometric free-particle continuation instead of
+// building domains on a field it can't parameterize. 1e-6 T is KinKal's null-field value
+// (BFieldMap::zeroField), so this only removes the degenerate samples. Reconstruction only: Geant4 ignores it.
+
+#include "Offline/Mu2eG4/geom/bfgeom_v01.txt"
+
+double bfield.minUsableField = 1.0e-6;
+
+//
+// This tells emacs to view this file in c++ mode.
+// Local Variables:
+// mode:c++
+// End:

Mu2eG4/geom/bfgeom_v01_minfield0p3.txt

@@ -0,0 +1,16 @@
+//
+// bfgeom_v01 with KinKal low-field protection at 0.3 T.
+//
+// bfield.minUsableField is the smallest |B| (T) KinKal uses for field-corrected transport; see
+// bfgeom_v01_minfield.txt. 0.3 T hands off to geometric continuation while real bending remains, so it
+// trades extrapolation accuracy for stability. Reconstruction only: Geant4 ignores it.
+
+#include "Offline/Mu2eG4/geom/bfgeom_v01.txt"
+
+double bfield.minUsableField = 0.3;
+
+//
+// This tells emacs to view this file in c++ mode.
+// Local Variables:
+// mode:c++
+// End:

Mu2eKinKal/fcl/prolog.fcl

@@ -53,6 +53,7 @@ Mu2eKinKal : {
      DivergenceGap : 20 # mm
      BFieldCorrection : false
      BCorrTolerance : 1e-2 # momemntum fraction
+     MinDomainTimeStep : 0.0 # (ns); 0 leaves the BField domain step unbounded
      ProcessEnds : false
      MetaIterationSettings : [
        # annealing temp, strawhit updater algorithm
@@ -89,6 +90,7 @@ Mu2eKinKal : {
     BFieldCorrection : true
     ProcessEnds : false
     BCorrTolerance : 1e-2 # relative momemntum precision
+    MinDomainTimeStep : 0.0 # (ns); 0 leaves the BField domain step unbounded
     MetaIterationSettings : [
       # annealing temp, strawhit updater algorithm
       [ 2.0, "CADSHU" ],
@@ -128,6 +130,7 @@ Mu2eKinKal : {
     BFieldCorrection : false
     ProcessEnds : false
     BCorrTolerance : 1e-4 # momemntum fraction
+    MinDomainTimeStep : 0.0 # (ns); 0 leaves the BField domain step unbounded
     MetaIterationSettings : [
       # annealing temp, strawhit updater algorithms
       [ 2.0, "CADSHU" ],
@@ -194,6 +197,7 @@ Mu2eKinKal : {
     DivergenceGap : 10 # mm
     BFieldCorrection : true
     BCorrTolerance : 1e-2 # momemntum fraction
+    MinDomainTimeStep : 0.0 # (ns); 0 leaves the BField domain step unbounded
     ProcessEnds : false
     MetaIterationSettings : [
       # annealing temp, strawhit updater algorithm
@@ -227,6 +231,7 @@ Mu2eKinKal : {
     BFieldCorrection : false
     ProcessEnds : false
     BCorrTolerance : 1e-4 # momemntum fraction
+    MinDomainTimeStep : 0.0 # (ns); 0 leaves the BField domain step unbounded
     MetaIterationSettings : []
     CADSHUSettings : []
     StrawXingUpdaterSettings : []
@@ -244,6 +249,7 @@ Mu2eKinKal : {
     BFieldCorrection : true
     ProcessEnds : false
     BCorrTolerance : 1e-2 # relative momemntum precision
+    MinDomainTimeStep : 0.0 # (ns); 0 leaves the BField domain step unbounded
     MetaIterationSettings : [
       # annealing temp, strawhit updater algorithm
       [ 2.0, "CADSHU" ],
@@ -283,6 +289,7 @@ Mu2eKinKal : {
     BFieldCorrection : true
     ProcessEnds : true
     BCorrTolerance : 1e-4 # momemntum fraction
+    MinDomainTimeStep : 0.0 # (ns); 0 leaves the BField domain step unbounded
     MetaIterationSettings : [
       # annealing temp, strawhit updater algorithms
       [ 2.0, "CADSHU" ],
@@ -358,6 +365,7 @@ Mu2eKinKal : {
     DivergenceGap : 100 # mm
     ProcessEnds : true
     BCorrTolerance : 1e-4
+    MinDomainTimeStep : 0.0 # (ns); 0 leaves the BField domain step unbounded
     MetaIterationSettings : [ [ 0.0, "" ] ]
     CADSHUSettings : [ ]
     BkgANNSHUSettings : [ ]

Mu2eKinKal/inc/KKBField.hh

@@ -17,7 +17,8 @@ namespace mu2e
       using Grad = ROOT::Math::SMatrix<double,3>; // field gradient: ie dBi/d(x,y,z)
       // construct from BField object and system translator.
       // This should be a single BField map valid in the detector system, to avoid making continuous translations FIXME!
-      KKBField(BFieldManager const& bfmgr, DetectorSystem const& det) : bfmgr_(bfmgr), det_(det) {}
+      // The minimum usable field comes from the BField geometry (bfield.minUsableField), so every module in a job shares it.
+      KKBField(BFieldManager const& bfmgr, DetectorSystem const& det);
       virtual ~KKBField() {}
       // KinKal BField interface
       // return value of the field at a poin

Mu2eKinKal/inc/KKExtrap.hh

@@ -19,6 +19,7 @@
 #include "Offline/Mu2eKinKal/inc/KKShellXing.hh"
 #include "Offline/KinKalGeom/inc/KKMaterial.hh"
 #include "KinKal/Geometry/ParticleTrajectoryIntersect.hh"
+#include "cetlib_except/exception.h"
 #include "Offline/GeometryService/inc/GeomHandle.hh"
 #include <unordered_set>
 #include <map>
@@ -41,6 +42,8 @@ namespace mu2e {
   class KKExtrap {
     public:
       explicit KKExtrap(KKExtrapConfig const& exconfig);
+      // throw at configuration if the fit config's minimum domain step can't fit this extrapolation's maximum step
+      void checkConfig(KinKal::Config const& config) const;
       // extrapolation functions; these are templated on the type of trajectory
       template <class KTRAJ> void extrapolate(KKTrack<KTRAJ>& ktrk) const;
       template <class KTRAJ> bool extrapolateIPA(KKTrack<KTRAJ>& ktrk,TimeDir trkdir) const;
@@ -80,6 +83,12 @@ namespace mu2e {
     toCRV_(extrapconfig.ToCRV())
   {}
 
+  inline void KKExtrap::checkConfig(KinKal::Config const& config) const {
+    if(config.bfcorr_ && config.mindtstep_ > maxdtstep_)
+      throw cet::exception("RECO") << "mu2e::KKExtrap: MinDomainTimeStep " << config.mindtstep_
+        << " ns exceeds ExtrapolationSettings.MaxDtStep " << maxdtstep_ << " ns" << std::endl;
+  }
+
   template <class KTRAJ> void KKExtrap::extrapolate(KKTrack<KTRAJ>& ktrk) const {
     GeomHandle<mu2e::KinKalGeom> kkg_h;
     // define the time direction according to the fit direction inside the tracker

Mu2eKinKal/inc/KKFitSettings.hh

@@ -4,6 +4,7 @@
 // Structs for configuring the Mu2e KinKal fit
 //
 #include "fhiclcpp/types/Atom.h"
+#include "fhiclcpp/types/OptionalAtom.h"
 #include "fhiclcpp/types/Sequence.h"
 #include "fhiclcpp/types/OptionalSequence.h"
 #include "fhiclcpp/types/Tuple.h"
@@ -35,6 +36,8 @@ namespace mu2e {
       fhicl::Atom<bool> bfieldCorr { Name("BFieldCorrection"), Comment("Apply correction for BField inhomogeneity") };
       fhicl::Atom<bool> ends { Name("ProcessEnds"), Comment("Process purely passive sites at the time range ends") };
       fhicl::Atom<float> btol { Name("BCorrTolerance"), Comment("Tolerance on BField correction momentum fractional accuracy (dimensionless)") };
+      fhicl::Atom<double> minDomainTimeStep { Name("MinDomainTimeStep"), Comment("Minimum BField domain time step (ns); 0 leaves it unbounded"), 0.0 };
+      fhicl::OptionalAtom<double> domainMargin { Name("DomainMargin"), Comment("Maximum time a BField domain may extend beyond the active range (ns, >= 0); absent leaves domains unconfined") };
       // Updater settings
       using MetaIterationSettings = fhicl::Sequence<fhicl::Tuple<float,std::string>>;
       MetaIterationSettings miConfig { Name("MetaIterationSettings"), Comment("Temperature (dimensionless), StrawHitUpdater algorithm") };

Mu2eKinKal/src/CentralHelixFit_module.cc

@@ -233,6 +233,8 @@ namespace mu2e {
       }
       // setup extrapolation
       if(settings().extrapSettings())extrap_ = make_unique<KKExtrap>(*settings().extrapSettings());
+      // extrapolation runs under the extension config when there is one
+      if(extrap_) extrap_->checkConfig(exconfig_.schedule().size() > 0 ? exconfig_ : config_);
 
       // surfaces to sample; this interface is deprecatecd and should be replaced with extrapolation TODO
       for(auto const& sidname : settings().modSettings().sampleSurfaces()) {
@@ -361,9 +363,13 @@ namespace mu2e {
             goodfit = goodFit(*ktrk);
           }
           // extrapolate as required
-          if(goodfit && extrap_)extrap_->extrapolate(*ktrk);
+          if(goodfit && extrap_){
+            extrap_->extrapolate(*ktrk);
+            goodfit = goodFit(*ktrk); // extrapolation can mark the fit failing
+          }
           if(print_>1)ktrk->printFit(std::cout,print_);
-          if(goodfit || saveall_){
+          // a fit that failed before its trajectory was built has nothing to save
+          if((goodfit || saveall_) && ktrk->hasTraj()){
             TrkFitFlag fitflag;
             fitflag.merge(fitflag_);
             if(goodfit)

Mu2eKinKal/src/KKBField.cc

@@ -1,9 +1,14 @@
 #include "Offline/Mu2eKinKal/inc/KKBField.hh"
+#include "Offline/BFieldGeom/inc/BFieldConfig.hh"
+#include "Offline/GeometryService/inc/GeomHandle.hh"
 #include "cetlib_except/exception.h"
 namespace mu2e {
   using Grad = ROOT::Math::SMatrix<double,3>;
   using SVEC3 = KinKal::SVEC3;
 
+  KKBField::KKBField(BFieldManager const& bfmgr, DetectorSystem const& det) :
+    KinKal::BFieldMap(GeomHandle<BFieldConfig>()->minUsableField()), bfmgr_(bfmgr), det_(det) {}
+
   VEC3 KKBField::fieldVect(VEC3 const& position) const {
     // change coordinates to mu2e; the map should be native in detector coordinates FIXME!
     CLHEP::Hep3Vector vpoint(position.x(),position.y(),position.z());

Mu2eKinKal/src/KKFitSettings.cc

@@ -2,6 +2,7 @@
 #include "Offline/Mu2eKinKal/inc/WireHitState.hh"
 #include "Offline/Mu2eKinKal/inc/StrawHitUpdaters.hh"
 #include "Offline/GeneralUtilities/inc/splitLine.hh"
+#include "cetlib_except/exception.h"
 #include <iostream>
 
 namespace mu2e {
@@ -22,6 +23,12 @@ namespace mu2e {
       config.bfcorr_ = fitconfig.bfieldCorr();
       config.ends_ = fitconfig.ends();
       config.tol_ = fitconfig.btol();
+      config.mindtstep_ = fitconfig.minDomainTimeStep();
+      double domainmargin(0.0);
+      if(fitconfig.domainMargin(domainmargin)){ // absent keeps KinKal's unconfined default
+        if(domainmargin < 0.0) throw cet::exception("RECO") << "mu2e::Mu2eKinKal::makeConfig: DomainMargin must be >= 0, got " << domainmargin << std::endl;
+        config.domainmargin_ = domainmargin;
+      }
       // create the updaters requested
       std::vector<CADSHU::Config> cadshusettings;
       std::vector<DriftANNSHU::Config> driftannshusettings;

Mu2eKinKal/src/KKTruthSeed_module.cc

@@ -145,6 +145,8 @@ namespace mu2e {
     produces<KalSeedCollection>();
     produces<KalSeedMCAssns>();
     if(settings().extrapSettings()) extrap_ = std::make_unique<KKExtrap>(*settings().extrapSettings());
+    // the truth track extrapolates with the extension config's bfcorr_ and mindtstep_
+    if(extrap_) extrap_->checkConfig(exconfig_);
     if(truthSeedParamConstraints_.size() != KinKal::NParams())
       throw cet::exception("RECO") << "mu2e::KKTruthSeed: TruthSeedParameterConstraints must have "
         << KinKal::NParams() << " entries" << std::endl;
@@ -219,6 +221,14 @@ namespace mu2e {
       auto const htime = fpiece.range().mid();
       VEC3 tpos, tmom; double ttime = 0.0;
       if(!trueStateAt(mctrajs, det, fpiece.position3(htime), mass, tpos, tmom, ttime)) continue;
+      // a helix needs a usable local field: outside the map, or where the helix can't represent the true state there
+      // (null field, circle centered on the axis), the constructor would refuse it
+      if constexpr(!isLine){
+        VEC4 tpos4(tpos.X(),tpos.Y(),tpos.Z(),ttime);
+        MOM4 tmom4(tmom.X(),tmom.Y(),tmom.Z(),mass);
+        if(!kkbf_->usable(tpos) ||
+            !KTRAJ::constructible(KinKal::ParticleState(tpos4,tmom4,tcharge),VEC3(0.0,0.0,kkbf_->fieldVect(tpos).Z()))) continue;
+      }
       // the KTRAJ nominal BField: the local field for a helix; a fixed non-zero z for a field-free line
       // (KinematicLine requires a non-zero bnom for interface consistency).
       KTRAJ tpiece = [&]{
@@ -289,6 +299,12 @@ namespace mu2e {
     // without re-running the full extension schedule (which diverges a hit-less truth-pinned fit).
     truthConfig.bfcorr_ = exconfig_.bfcorr_;
     truthConfig.tol_    = exconfig_.tol_;
+    truthConfig.mindtstep_ = exconfig_.mindtstep_; // floors the extrapolation step too
+    // a field-corrected truth fit needs the reco domains, as in RegrowLoopHelix
+    if(truthConfig.bfcorr_ && truthDomains.empty()){
+      if(debug_ > 0) std::cout << "KKTruthSeed: no BField domains from " << recoseed.domainBounds().size() << " saved bounds; skipping" << std::endl;
+      return std::nullopt;
+    }
     auto ktrk_truth = std::make_unique<KKTRK>(truthConfig,*kkbf_,recoseed.particle(),truthTraj,
         nostrawhits,nostrawxings,nocalohits,truthParamHits,truthDomains);
     if(!ktrk_truth->fitStatus().usable()){

Mu2eKinKal/src/KinematicLineFit_module.cc

@@ -219,6 +219,8 @@ namespace mu2e {
       }
       // setup extrapolation
       if(settings().extrapSettings())extrap_ = make_unique<KKExtrap>(*settings().extrapSettings());
+      // extrapolation runs under the extension config when there is one
+      if(extrap_) extrap_->checkConfig(exconfig_.schedule().size() > 0 ? exconfig_ : config_);
       if(print_ > 0) std::cout << config_;
     }
 
@@ -316,7 +318,8 @@ namespace mu2e {
           // extrapolate as required
           if(goodfit && extrap_) extrap_->extrapolate(*kktrk);
           bool save = goodFit(*kktrk);
-          if(save || saveall_){
+          // a fit that failed before its trajectory was built has nothing to save
+          if((save || saveall_) && kktrk->hasTraj()){
             TrkFitFlag fitflag(hptr->status());
             fitflag.merge(TrkFitFlag::KKLine);
             sampleFit(*kktrk);

Mu2eKinKal/src/LoopHelixFit_module.cc

@@ -225,6 +225,8 @@ namespace mu2e {
       }
       // setup extrapolation
       if(settings().extrapSettings())extrap_ = make_unique<KKExtrap>(*settings().extrapSettings());
+      // extrapolation runs under the extension config when there is one
+      if(extrap_) extrap_->checkConfig(exconfig_.schedule().size() > 0 ? exconfig_ : config_);
 
       // setup optional fit finalization; this just updates the internals, not the fit result itself
       if(settings().finalSettings()){
@@ -366,7 +368,8 @@ namespace mu2e {
     }
     if(print_>0) printf("[LoopHelixFit::%s] After extending the fit : goodFit = %o, fitcon = %.4f, nHits = %2lu, %lu calo-hits\n",
         __func__, goodfit, ktrk->fitStatus().chisq_.probability(), ktrk->strawHits().size(), ktrk->caloHits().size());
-    if((!goodfit) && (! saveall_)) ktrk.reset();
+    // a fit that failed before its trajectory was built has nothing to save
+    if(!ktrk->hasTraj() || ((!goodfit) && (! saveall_))) ktrk.reset();
     return ktrk;
   }
 

RecoDataProducts/inc/TrkStrawHitSeed.hh

@@ -172,7 +172,7 @@ namespace mu2e {
     float         _stoca =0;    // reference sensor time of closest approach (TOCA)
     float         _rdoca, _rdocavar =0;   // reference (biased) DOCA from the track to the wire, signed by the angular momentum WRT the wire and the measurement end (and variance)
     float         _rdt, _rtocavar =0;   // reference (biased) time difference (and variance) at POCA
-    float         _udoca, _udocavar =0; // unbiaed DOCA (and variance)
+    float         _udoca, _udocavar =0; // unbiased DOCA (and variance) WRT the final fit trajectory; the hit-state updaters (ANNs) used the unbiased DOCA from the previous fit cycle.
     float         _udt, _utocavar =0;   //unbiased dt and variance
     float         _rupos =0; // reference POCA position along the straw WRT the straw middle
     float         _uupos =0; // unbiased POCA position along the straw WRT the straw middle

scripts/build/SConscript

@@ -41,4 +41,3 @@ exe_env.AlwaysBuild(env.Alias('INCPATH', [], incpath))
 # mode:python
 # End:
 # vi:syntax=python
-