Fix critical bugs in KinKal integration code

Mu2eKinKal/inc/Chi2SHU_updateCluster.hh

@@ -33,7 +33,7 @@ namespace mu2e {
     // in marginal fits removing even a few hits can leave the fit underconstrained, resulting in a zero determinant.  For now, do nothing
     // with these clusters
     double determinant;
-    if(!uparams.covariance().Det(determinant) || determinant < std::numeric_limits<float>::min()){
+    if(!uparams.covariance().Det(determinant) || determinant < std::numeric_limits<double>::min()){
       if(diag_ > 2)std::cout << "Negative unbiased covar determinant = " << determinant << std::endl;
       return;
     }

Mu2eKinKal/inc/DriftANNSHU.hh

@@ -34,7 +34,7 @@ namespace mu2e {
       double dtmvacut_ =0; // cut value for using dt constraint
       KKSHFlag flag_ = KKSHFlag(KKSHFlag::tot); // constrain time with TOT by default
       WHSMask freeze_; // states to freeze
-      int diag_; // diag print level
+      int diag_ = 0; // diag print level
   };
 }
 #endif

Mu2eKinKal/inc/ExtrapolateIPA.hh

@@ -17,8 +17,8 @@ namespace mu2e {
     public:
       using CylPtr = std::shared_ptr<KinKal::Cylinder>;
       ExtrapolateIPA() : maxDt_(-1.0), dptol_(1e10), intertol_(1e10),
-      zmin_(std::numeric_limits<float>::max()),
-      zmax_(-std::numeric_limits<float>::max()),debug_(0) {}
+      zmin_(std::numeric_limits<double>::max()),
+      zmax_(std::numeric_limits<double>::lowest()),debug_(0) {}
 
       ExtrapolateIPA(double maxdt, double dptol,double intertol, CylPtr const& ipa, int debug=0) :
         maxDt_(maxdt), dptol_(dptol), intertol_(intertol), ipa_(ipa),

Mu2eKinKal/inc/ExtrapolateST.hh

@@ -24,10 +24,10 @@ namespace mu2e {
   class ExtrapolateST {
     public:
       ExtrapolateST() : maxDt_(-1.0), dptol_(1e10), intertol_(1e10),
-      zmin_(std::numeric_limits<float>::max()),
-      zmax_(-std::numeric_limits<float>::max()),
-      rmin_(std::numeric_limits<float>::max()),
-      rmax_(-std::numeric_limits<float>::max()),
+      zmin_(std::numeric_limits<double>::max()),
+      zmax_(std::numeric_limits<double>::lowest()),
+      rmin_(std::numeric_limits<double>::max()),
+      rmax_(std::numeric_limits<double>::lowest()),
       debug_(0){}
 
       ExtrapolateST(double maxdt, double dptol, double intertol, StoppingTarget const& stoptarg, int debug=0) :

Mu2eKinKal/inc/ExtrapolateTCRV.hh

@@ -22,8 +22,8 @@ namespace mu2e {
     public:
       ExtrapolateTCRV() : maxDt_(-1.0), dptol_(1e10), intertol_(1e10), minv_(1e-5),
       step_(0),
-      ymin_(std::numeric_limits<float>::max()),
-      ymax_(-std::numeric_limits<float>::max()),
+      ymin_(std::numeric_limits<double>::max()),
+      ymax_(std::numeric_limits<double>::lowest()),
       debug_(0){}
 
       ExtrapolateTCRV(double maxdt, double dptol, double intertol, double minv, TestCRV const& tcrv, int debug=0) :

Mu2eKinKal/inc/KKFit.hh

@@ -146,10 +146,10 @@ namespace mu2e {
       float maxStrawHitDoca_, maxStrawHitDt_, maxStrawDoca_, maxStrawDocaCon_, maxStrawUposBuff_;
       int maxDStraw_; // maximum distance from the track a strawhit can be to consider it for adding.
       // cached info computed from the tracker, used in hit adding; these must be lazy-evaluated as the tracker doesn't exist on construction
-      mutable double strawradius_;
-      mutable double ymin_, ymax_, umax_; // panel-level info
-      mutable double rmin_, rmax_; // plane-level info
-      mutable double spitch_;
+      mutable double strawradius_ = 0;
+      mutable double ymin_ = 0, ymax_ = 0, umax_ = 0; // panel-level info
+      mutable double rmin_ = 0, rmax_ = 0; // plane-level info
+      mutable double spitch_ = 0;
       mutable bool needstrackerinfo_ = true;
       // extrapolation and sampling options
       SurfaceMap::SurfacePairCollection sample_; // surfaces to sample the fit
@@ -476,7 +476,7 @@ namespace mu2e {
                 // require consistency with this track passing through this straw
                 double du = fabs((pca.sensorPoca().Vect()-VEC3(straw.wirePosition(0.0))).Dot(VEC3(straw.wireDirection(0.0))));
                 double doca = fabs(pca.doca());
-                double dsig = std::max(0.0,doca-strawradius_)/sqrt(pca.docaVar());
+                double dsig = std::max(0.0,doca-strawradius_)/sqrt(std::max(pca.docaVar(),std::numeric_limits<double>::min()));
                 if(doca < maxStrawDoca_ && dsig < maxStrawDocaCon_ && du < straw.halfLength() + maxStrawUposBuff_){
                   addexings.push_back(std::make_shared<KKSTRAWXING>(shptr,static_cast<CA>(pca),smat,straw));
                   oldstraws.insert(straw.id());
@@ -533,7 +533,7 @@ namespace mu2e {
                       // require consistency with this track passing through this straw
                       double du = fabs((pca.sensorPoca().Vect()-VEC3(straw.wirePosition(0.0))).Dot(VEC3(straw.wireDirection(0.0))));
                       double doca = fabs(pca.doca());
-                      double dsig = std::max(0.0,doca-strawradius_)/sqrt(pca.docaVar());
+                      double dsig = std::max(0.0,doca-strawradius_)/sqrt(std::max(pca.docaVar(),std::numeric_limits<double>::min()));
                       if(doca < maxStrawDoca_ && dsig < maxStrawDocaCon_ && du < straw.halfLength() + maxStrawUposBuff_){
                         addexings.push_back(std::make_shared<KKSTRAWXING>(shptr,static_cast<CA>(pca),smat,straw));
                         oldstraws.insert(straw.id());
@@ -625,8 +625,8 @@ namespace mu2e {
 
 
   template <class KTRAJ> TimeRange KKFit<KTRAJ>::range(KKSTRAWHITCOL const& strawhits, KKCALOHITCOL const& calohits, KKSTRAWXINGCOL const& strawxings) const{
-    double tmin = std::numeric_limits<float>::max();
-    double tmax = -tmin;
+    double tmin = std::numeric_limits<double>::max();
+    double tmax = std::numeric_limits<double>::lowest();
     for( auto const& strawhit : strawhits) {
       tmin = std::min(tmin,strawhit->time());
       tmax = std::max(tmax,strawhit->time());
@@ -763,8 +763,8 @@ namespace mu2e {
         }
       } else if (savetraj_ == detector ) {
         // only save segments inside the tracker volume. Find the limits for that. Start with the times of active hits
-        double tmin = std::numeric_limits<float>::max();
-        double tmax = -tmin;
+        double tmin = std::numeric_limits<double>::max();
+        double tmax = std::numeric_limits<double>::lowest();
         for(auto const& kkshp : kktrk.strawHits()){
           if(kkshp->active()){
             tmin = std::min(tmin,kkshp->time());

Mu2eKinKal/inc/KKShellXing.hh

@@ -101,7 +101,8 @@ namespace mu2e {
   }
 
   template <class KTRAJ,class SURF> double KKShellXing<KTRAJ,SURF>::transitTime() const {
-    double pathlen = thick_/(inter_.norm_.Dot(inter_.pdir_));
+    double dotprod = std::max(1e-6,fabs(inter_.norm_.Dot(inter_.pdir_)));
+    double pathlen = thick_/dotprod;
     return pathlen/reftrajptr_->speed();
   }
 

Mu2eKinKal/inc/KKStrawHit.hh

@@ -54,8 +54,8 @@ namespace mu2e {
       KKStrawHit(KKStrawHit<KTRAJ> const& rhs):
           bfield_(rhs.bfield()),
           whstate_(rhs.hitState()),
-          dVar_(driftVariance()),
-          dDdT_(driftVelocity()),
+          dVar_(rhs.dVar_),
+          dDdT_(rhs.dDdT_),
           ca_(rhs.closestApproach()),
           resids_(rhs.refResiduals()),
           chit_(rhs.hit()),

Mu2eKinKal/inc/KKStrawHitCluster.hh

@@ -143,7 +143,7 @@ namespace mu2e {
     // return time just before the first hit's time.  This insures hit clusters are updated before individual hits
     // This insures the weights subtracted correspond to the reference fit, and that any changes made to the
     // hits in the cluster get propagated to the residuals and weights before the next fit
-    double mintime(std::numeric_limits<float>::max());
+    double mintime(std::numeric_limits<double>::max());
     for(auto const& hit : hits_){
       mintime = std::min(hit->time(),mintime);
     }

Mu2eKinKal/src/Chi2SHU.cc

@@ -79,8 +79,8 @@ namespace mu2e {
 
   std::ostream& operator <<(std::ostream& os, ClusterState const& cstate ) {
     os << "ClusterState " << cstate.chi2_ << " states: ";
-    for(auto whstate : cstate.hitstates_) std::cout << "  " << whstate.state_;
-    std::cout << std::endl;
+    for(auto whstate : cstate.hitstates_) os << "  " << whstate.state_;
+    os << std::endl;
     return os;
   }
 

Mu2eKinKal/src/DriftANNSHU.cc

@@ -23,7 +23,7 @@ namespace mu2e {
     diag_ = std::get<7>(config);
     if(diag_ > 0)
       std::cout << "DriftANNSHU LR sign weights " << std::get<0>(config) << " cut " << signmvacut_
-        << " cluster weights " << std::get<0>(config) << " cut " << clustermvacut_ << " dt cut " << dtmvacut_
+        << " cluster weights " << std::get<2>(config) << " cut " << clustermvacut_ << " dt cut " << dtmvacut_
         << " freezing " << freeze_ << " flags " << flag << " diag " << diag_ << std::endl;
   }
 

Mu2eKinKal/src/KKFitUtilities.cc

@@ -42,8 +42,8 @@ namespace mu2e {
 
     KinKal::TimeRange timeBounds(ComboHitCollection const& chits) {
       if(chits.size() == 0) return KinKal::TimeRange();
-      double tmin = std::numeric_limits<float>::max();
-      double tmax = std::numeric_limits<float>::min();
+      double tmin = std::numeric_limits<double>::max();
+      double tmax = std::numeric_limits<double>::lowest();
       for( auto const& hit : chits) {
         // filter out hits already flagged as bad TODO
         tmin = std::min(tmin,(double)hit.correctedTime());
@@ -54,8 +54,8 @@ namespace mu2e {
 
     double zMid(ComboHitCollection const& chits) {
       if(chits.size() == 0) return 0.0;
-      double zmin = std::numeric_limits<float>::max();
-      double zmax = std::numeric_limits<float>::min();
+      double zmin = std::numeric_limits<double>::max();
+      double zmax = std::numeric_limits<double>::lowest();
       for( auto const& hit : chits) {
         // filter out hits already flagged as bad TODO
         double zpos = hit.pos().z();

Mu2eKinKal/src/KKStrawMaterial.cc

@@ -31,7 +31,7 @@ namespace mu2e {
   KKStrawMaterial::PathCalc KKStrawMaterial::averagePathLengths(double& wallpath, double& gaspath, double& wirepath) const {
     gaspath = avggaspath_;
     wallpath = avgwallpath_;
-    wallpath = 0.0;
+    wirepath = 0.0;
     return KKStrawMaterial::average;
   }
 

Mu2eKinKal/src/KinematicLineFit_module.cc

@@ -337,7 +337,6 @@ namespace mu2e {
             kkseedcol->push_back(kkseed);
             kkseedcol->back()._status.merge(TrkFitFlag::KKLine);
             // fill assns with the cosmic seed
-            auto hptr = art::Ptr<CosmicTrackSeed>(hseedcol_h,iseed);
             auto kseedptr = art::Ptr<KalSeed>(KalSeedCollectionPID,kkseedcol->size()-1,KalSeedCollectionGetter);
             kkseedassns->addSingle(kseedptr,hptr);
             // save (unpersistable) KKTrk in the event

Mu2eKinKal/src/LoopHelixFit_module.cc

@@ -478,7 +478,7 @@ namespace mu2e {
       // require physical consistency: fit can succeed but the result can have changed charge or helicity. Test at the t0 segment
       auto t0 = Mu2eKinKal::zTime(ktrk.fitTraj(),0.0,ktrk.fitTraj().range().mid());
       auto const& t0seg = ktrk.fitTraj().nearestPiece(t0);
-      bool retval = ktrk.fitStatus().usable() && t0seg.parameterSign()*seed.parameterSign() > 0 && t0seg.helicity()*seed.helicity() > 0;
+      retval = t0seg.parameterSign()*seed.parameterSign() > 0 && t0seg.helicity()*seed.helicity() > 0;
       // also check that the fit is inside the physical detector volume.  Test where the StrawHits are
       if(retval){
         for(auto const& shptr : ktrk.strawHits()) {