diff --git a/Fit/Status.hh b/Fit/Status.hh
index 036e907a..f94fa954 100644
--- a/Fit/Status.hh
+++ b/Fit/Status.hh
@@ -16,6 +16,7 @@ namespace KinKal {
     Chisq chisq_; // current chisquared
     std::string comment_; // further information about the status
     bool usable() const { return status_ > unfit && status_ < lowNDOF; }
+    bool unusable() const { return !usable(); }
     bool needsFit() const { return status_ == unfit || status_ == unconverged; }
     Status(unsigned miter,unsigned iter=0,status stat=unfit,const char* comment="") : miter_(miter), iter_(iter), status_(stat), chisq_(NParams()),comment_(comment){}
     static std::string statusName(status stat);
diff --git a/Fit/Track.hh b/Fit/Track.hh
index 2230912a..237f2982 100644
--- a/Fit/Track.hh
+++ b/Fit/Track.hh
@@ -159,7 +159,7 @@ namespace KinKal {
       KKEFFCOL effects_; // effects used in this fit, sorted by time
       HITCOL hits_; // hits used in this fit
       EXINGCOL exings_; // material xings used in this fit
-      DOMAINCOL domains_; // DomainWall domains used in this fit, contiguous and sorted by time
+      DOMAINCOL domains_; // domains used in this fit, contiguous and sorted by time
   };
   // minimal constructor for subclasses. The resulting object has no fit.
   template <class KTRAJ> Track<KTRAJ>::Track(Config const& cfg, BFieldMap const& bfield) :
@@ -248,12 +248,13 @@ namespace KinKal {
   template <class KTRAJ> void Track<KTRAJ>::extend(Config const& cfg, HITCOL& hits, EXINGCOL& exings) {
     // require the existing fit to be usable
     if(!fitStatus().usable())return;
-    // remember the previous config
-    auto const& oldconfig = config_.back();
-    // update the configuration
-    config_.push_back(cfg);
     // configuation check
-    if(config().schedule().size() ==0)throw std::invalid_argument("Invalid configuration: no schedule");
+    if(cfg.schedule().size() ==0)throw std::invalid_argument("Invalid configuration: no schedule");
+    // save the new configuration
+    config_.push_back(cfg);
+    // remember the previous config
+    auto oldciter = config_.rbegin(); oldciter++;
+    auto const& oldconfig = *oldciter;
     // find the range of the added information, and extend as needed
     TimeRange exrange = fittraj_->range();
     if(hits.size() >0 || exings.size() > 0){
@@ -302,7 +303,7 @@ namespace KinKal {
   }
 
   // replace domains when DomainWall correction is added or changed. the traj must also be replaced, so that
-  // the pieces correspond to the new domains
+  // the pieces correspond to the new domains. The new traj is geometrically equivalent, but not parametrically equal.
   template <class KTRAJ> void Track<KTRAJ>::replaceDomains(DOMAINCOL const& domains) {
     // if domains exist, clear them and remove all DomainWall effects
     if(domains_.size() > 0){
@@ -319,7 +320,11 @@ namespace KinKal {
       }
     }
     auto newtraj = std::make_unique<PKTRAJ>();
-    // loop over domains
+    // loop over domains, splitting the overlapping traj pieces at the domain walls, and transforming them to reference the domain's field
+    // This increases the number of traj pieces.
+    // extend the existing traj to the domain range
+    TimeRange drange(domains.begin()->get()->begin(),domains.rbegin()->get()->end());
+    fittraj_->setRange(drange);
     for(auto const& domain : domains) {
       // find the range of existing ptraj pieces that overlaps with this domain's range
       using KTRAJPTR = std::shared_ptr<KTRAJ>;
@@ -327,9 +332,9 @@ namespace KinKal {
       using DKTRAJCITER = typename DKTRAJ::const_iterator;
       DKTRAJCITER first,last;
       fittraj_->pieceRange(domain->range(),first,last);
-      // loop over these pieces
+      // loop over these pieces; first and last can be the same!
       auto olditer = first;
-      while(olditer != last){
+      do {
         auto const& oldpiece = **olditer;
         // copy this piece, translating bnom to this domain's field
         KTRAJ newpiece(oldpiece,domain->bnom(),domain->range().mid());
@@ -340,14 +345,14 @@ namespace KinKal {
           newpiece.range() = TimeRange(tstart,tend);
           newtraj->append(newpiece);
         }
-        ++olditer;
-      }
+        if(olditer != last)++olditer;
+      } while(olditer != last);
     }
     // switch over any existing effects to reference this traj (could be none)
     for (auto& eff : effects_) {
       eff->updateReference(*newtraj);
     }
-    // swap out the fit
+    // swap out the fit trajectory; this will be used as reference for the next iterations
     fittraj_.swap(newtraj);
   }
 
@@ -694,7 +699,10 @@ namespace KinKal {
     }
     KKEFFFWDBND fwdbnds; // bounding sites used for fitting
     KKEFFREVBND revbnds;
-    setBounds(fwdbnds,revbnds);
+    if(!setBounds(fwdbnds,revbnds)){
+      status().status_ = Status::lowNDOF;
+      return;
+    }
     // initialize the fit state where we left off processing
     FitStateArray states;
     TimeRange fitrange(fwdbnds[0]->get()->time(),revbnds[0]->get()->time());
diff --git a/General/TimeRange.cc b/General/TimeRange.cc
index d8bf0aa4..81fe4b98 100644
--- a/General/TimeRange.cc
+++ b/General/TimeRange.cc
@@ -1,5 +1,24 @@
 #include "KinKal/General/TimeRange.hh"
 namespace KinKal {
+
+  bool TimeRange::restrict(TimeRange const& other ) {
+    bool retval = overlaps(other);
+    if(retval){
+      range_[0] = std::max(begin(),other.begin());
+      range_[1] = std::min(end(),other.end());
+    }
+    return retval;
+  }
+
+  void TimeRange::extend(double time, TimeDir tdir) {
+    // require the resulting range to be >0
+    if(tdir == TimeDir::forwards){
+      if(time > range_[0])range_[1] = time;
+    } else {
+      if(time < range_[1])range_[0] = time;
+    }
+  }
+
   std::ostream& operator <<(std::ostream& ost, TimeRange const& trange) {
     ost << " Range [" << trange.begin() << "," << trange.end() << "]";
     return ost;
diff --git a/General/TimeRange.hh b/General/TimeRange.hh
index 00673042..46875a82 100644
--- a/General/TimeRange.hh
+++ b/General/TimeRange.hh
@@ -5,6 +5,8 @@
 #include <ostream>
 #include <array>
 #include <stdexcept>
+#include "KinKal/General/TimeDir.hh"
+
 namespace KinKal {
   class TimeRange {
     public:
@@ -32,14 +34,9 @@ namespace KinKal {
       }
       // restrict the range to the overlap with another range. If there is no overlap
       // leave the object unchanged and return 'false'
-      bool restrict(TimeRange const& other ) {
-        bool retval = overlaps(other);
-        if(retval){
-          range_[0] = std::max(begin(),other.begin());
-          range_[1] = std::min(end(),other.end());
-        }
-        return retval;
-      }
+      bool restrict(TimeRange const& other );
+      // extend in a given direction
+      void extend(double time, TimeDir tdir);
     private:
       std::array<double,2> range_; // range of times
   };
diff --git a/MatEnv/DetMaterial.cc b/MatEnv/DetMaterial.cc
index d5a03fc2..f11d3a0d 100644
--- a/MatEnv/DetMaterial.cc
+++ b/MatEnv/DetMaterial.cc
@@ -30,16 +30,11 @@ using std::ostream;
 //
 namespace MatEnv {
 
-  //double DetMaterial::_msmom = 15.0*MeV;
   double DetMaterial::_dgev = 0.153536e2;
-  double DetMaterial::_minkappa(1.0e-3);
-  //double DetMaterial::_scatterfrac(0.9999); // integrate 99.99% percent of the tail by default, this should be larger
   // if the materials are very thin.
   const double bg2lim = 0.0169;
   const double taulim = 8.4146e-3 ;
   const double twoln10 = 2.0*log(10.);
-  const double betapower = 1.667; // most recent PDG gives beta^-5/3 as dE/dx
-  const int maxnstep = 10; // maximum number of steps through a single material
   // should be from a physics class
   const double DetMaterial::_alpha(1.0/137.036);
 
@@ -47,20 +42,15 @@ namespace MatEnv {
   //  Constructor
   //
 
-  double cm(10.0); // temporary hack
-  DetMaterial::DetMaterial(const char* detMatName, const MtrPropObj* detMtrProp):
-    _elossmode(mpv), //Energy Loss model: choose 'mpv' for the Most Probable Energy Loss, or 'moyalmean' for the mean calculated via the Moyal Distribution approximation, see end of file for more information, as well as discussion about radiative losses
-    _msmom(15.0),
-    _scatterfrac(0.9999),
-    _cutOffEnergy(1000.),
-    _elossType(loss),
+  double cm(10.0); // convert cm to mm
+  DetMaterial::DetMaterial(const char* detMatName, const MtrPropObj* detMtrProp, DetMaterialConfig const& dmconf):
+    _elossmode(dmconf.elossmode_),
     _name(detMatName),
     _za(detMtrProp->getZ()/detMtrProp->getA()),
     _zeff(detMtrProp->getZ()),
     _aeff(detMtrProp->getA()),
     _radthick(detMtrProp->getRadLength()/cm/cm),
     _intLength(detMtrProp->getIntLength()/detMtrProp->getDensity()),
-    _meanion(2.*log(detMtrProp->getMeanExciEnergy()*1.0e6)),
     _eexc(detMtrProp->getMeanExciEnergy()),
     _x0(detMtrProp->getX0density()),
     _x1(detMtrProp->getX1density()),
@@ -74,124 +64,70 @@ namespace MatEnv {
   {
     _shellCorrectionVector =
       new std::vector< double >(detMtrProp->getShellCorrectionVector());
-    _vecNbOfAtomsPerVolume =
-      new std::vector< double >(detMtrProp->getVecNbOfAtomsPerVolume());
-    _vecTau0 =
-      new std::vector< double >(detMtrProp->getVecTau0());
-    _vecAlow = new std::vector< double >(detMtrProp->getVecAlow());
-    _vecBlow = new std::vector< double >(detMtrProp->getVecBlow());
-    _vecClow = new std::vector< double >(detMtrProp->getVecClow());
-    _vecZ = new std::vector< double >(detMtrProp->getVecZ());
     // compute cached values; these are used in detailed scattering models
     _invx0 = _density/_radthick;
-    _nbar = _invx0*1.587e7*pow(_zeff,1.0/3.0)/((_zeff+1)*log(287/sqrt(_zeff)));
     _chic2 = 1.57e1*_zeff*(_zeff+1)/_aeff;
     _chia2_1 = 2.007e-5*pow(_zeff,2.0/3.0);
     _chia2_2 = 3.34*pow(_zeff*_alpha,2);
 
-    if (detMtrProp->getEnergyTcut()>0.0) {
-      _cutOffEnergy = detMtrProp->getEnergyTcut();
-      _elossType = deposit;
-    }
     if (detMtrProp->getState() == "gas" && detMtrProp->getDensity()<0.01) {
-      _scatterfrac = 0.999999;
+      _scatterfrac = dmconf.scatterfrac_gas_;
+    } else {
+      _scatterfrac = dmconf.scatterfrac_solid_;
     }
+    // compute the sum radiated photons*energy for the given cuttoff
+    double ymax = dmconf.ebrehmsfrac_;
+    // energy loss through 1 radiation length of this material
+    _erad = (_density/_radthick)*(4.0*ymax -2.0*pow(ymax,2) + pow(ymax,3)); // energy loss per radiation length per MeV
+    _eradvar = pow(_density/_radthick,2)*(2.0*pow(ymax,2)/3.0 - 4.0*pow(ymax,3)/9.0 + pow(ymax,4)/4.0); // energy variance per (rad length per MeV)^2
   }
 
   DetMaterial::~DetMaterial()
   {
     delete _shellCorrectionVector;
-    delete _vecNbOfAtomsPerVolume;
-    delete _vecTau0;
-    delete _vecAlow;
-    delete _vecBlow;
-    delete _vecClow;
-    delete _vecZ;
   }
 
   //
   //  Multiple scattering function
   //
-  double
-    DetMaterial::scatterAngleRMS(double mom, double pathlen,double mass) const {
-      if(mom>0.0){
-        double beta = particleBeta(mom,mass);
-        // pdg formulat
-        //    double radfrac = fabs(pathlen*_invx0);
-        //    double sigpdg = 0.0136*sqrt(radfrac)*(1.0+0.088*log10(radfrac))/(beta*mom);
-        // old Kalman formula
-        //    double oldsig = 0.011463*sqrt(radfrac)/(mom*particleBeta(mom,mass));
-        // DNB 20/1/2011  Updated to use Dahl-Lynch formula from  NIMB58 (1991)
-        double invmom2 = 1.0/pow(mom,2);
-        double invb2 = 1.0/pow(beta,2);
-        // convert to path in gm/cm^2!!!
-        double path = fabs(pathlen)*_density;
-        double chic2 = _chic2*path*invb2*invmom2;
-        double chia2 = _chia2_1*(1.0 + _chia2_2*invb2)*invmom2;
-        double omega = chic2/chia2;
-        static double vfactor = 0.5/(1-_scatterfrac);
-        double v = vfactor*omega;
-        static double sig2factor = 1.0/(1+_scatterfrac*_scatterfrac);
-        double sig2 = sig2factor*chic2*( (1+v)*log(1+v)/v - 1);
-        // protect against underflow
-        double sigdl = sqrt(std::max(0.0,sig2));
-        // check
-        return sigdl;
-      } else
-        return 1.0; // 'infinite' scattering
-    }
-
-  double
-    DetMaterial::dEdx(double mom,dedxtype type,double mass) const {
-      if(mom>0.0){
-        double Eexc2 = _eexc*_eexc ;
-
-        // New energy loss implementation
-        double Tmax,gamma2,beta2,bg2,rcut,delta,sh,dedx ;
-        double beta  = particleBeta(mom,mass) ;
-        double gamma = particleGamma(mom,mass) ;
-        double tau = gamma-1. ;
-
-        // high energy part , Bethe-Bloch formula
-        beta2 = beta*beta ;
-        gamma2 = gamma*gamma ;
-        bg2 = beta2*gamma2 ;
-
-        double RateMass = e_mass_/ mass;
-
-        Tmax = 2.*e_mass_*bg2
-          /(1.+2.*gamma*RateMass+RateMass*RateMass) ;
-
-        dedx = log(2.*e_mass_*bg2*Tmax/Eexc2);
-        if(type == loss)
-          dedx -= 2.*beta2;
-        else {
-          rcut =  ( _cutOffEnergy< Tmax) ? _cutOffEnergy/Tmax : 1;
-          dedx += log(rcut)-(1.+rcut)*beta2;
-        }
-
-        //// density correction
-        delta = densityCorrection(bg2);
-
-        //// shell correction
-        sh = shellCorrection(bg2, tau);
-
-        dedx -= delta + sh ;
-        dedx *= -_dgev*_density*_za / beta2 ;
+  double DetMaterial::scatterAngleVar(double mom, double pathlen,double mass) const {
+    if(mom>0.0){
+      double beta = particleBeta(mom,mass);
+      // DNB 20/1/2011  Updated to use Dahl-Lynch formula from  NIMB58 (1991)
+      double invmom2 = 1.0/pow(mom,2);
+      double invb2 = 1.0/pow(beta,2);
+      // convert to path in gm/cm^2!!!
+      double path = fabs(pathlen)*_density;
+      double chic2 = _chic2*path*invb2*invmom2;
+      double chia2 = _chia2_1*(1.0 + _chia2_2*invb2)*invmom2;
+      double omega = chic2/chia2;
+      static double vfactor = 0.5/(1-_scatterfrac);
+      double v = vfactor*omega;
+      static double sig2factor = 1.0/(1+_scatterfrac*_scatterfrac);
+      double sig2 = sig2factor*chic2*( (1+v)*log(1+v)/v - 1);
+      // protect against underflow
+      return std::max(0.0,sig2);
+    } else
+      return 1.0; // 'infinite' scattering
+  }
 
-        return dedx;
+  double DetMaterial::energyLoss(double mom,double pathlen,double mass) const {
+    double retval = ionizationEnergyLoss(mom,pathlen,mass);
+    retval += radiationEnergyLoss(mom,pathlen,mass);
+    return retval;
+  }
 
-      } else
-        return 0.0;
-    }
+  double DetMaterial::energyLossVar(double mom,double pathlen,double mass) const {
+    double retval = ionizationEnergyLossVar(mom,pathlen,mass);
+    retval += radiationEnergyLossVar(mom,pathlen,mass);
+    return retval;
+  }
 
-  //Replacement for dEdx-based energy loss function: Most probable energy loss is now used
-  //from https://pdg.lbl.gov/2019/reviews/rpp2018-rev-passage-particles-matter.pdf
-  double DetMaterial::energyLoss(double mom, double pathlen, double mass) const {
+  double DetMaterial::ionizationEnergyLoss(double mom, double pathlen, double mass) const {
     if(mom>0.0){
       double beta  = particleBeta(mom,mass) ;
       double xi = eloss_xi(beta, pathlen);
-      double deltap = energyLossMPV(mom,pathlen,mass);
+      double deltap = ionizationEnergyLossMPV(mom,pathlen,mass);
       //if using mean calculated from the Moyal Dist. Approx: (see end of file for more information)
       if(_elossmode == moyalmean) {
         return moyalMean(deltap, xi);
@@ -201,7 +137,8 @@ namespace MatEnv {
       return 0.0;
   }
 
-  double DetMaterial::energyLossMPV(double mom, double pathlen, double mass) const {
+  //Most probable energy loss from https://pdg.lbl.gov/2019/reviews/rpp2018-rev-passage-particles-matter.pdf
+  double DetMaterial::ionizationEnergyLossMPV(double mom, double pathlen, double mass) const {
     if(mom>0.0){
       //taking positive lengths
       pathlen = fabs(pathlen) ;
@@ -213,7 +150,6 @@ namespace MatEnv {
       double j = 0.200 ;
       double tau = gamma-1;
 
-
       // most probable energy loss function
       beta2 = beta*beta ;
       gamma2 = gamma*gamma ;
@@ -237,285 +173,137 @@ namespace MatEnv {
       return 0.0;
   }
 
+  double DetMaterial::radiationEnergyLoss(double mom,double pathlen, double mass) const {
+    // truncated average bremhsstrahlung radiation energy loss
+    // Based on integrating the 'low-y' approximation of the brehms cross-section from 0 to ymax, using
+    // equation 34.29 of the RPP 'Passage of particles through matter', Phys. Rev. D 110, 030001 (2024)
+    double retval(0.0);
+    static const double small(1e-3);
+    if(fabs(mass-e_mass_)<small){ // crude test to select electrons; should be replaced by a propoer particle enum TODO
+      retval = -mom*pathlen*_erad;
+    }
+    return retval;
+  }
+  double DetMaterial::radiationEnergyLossVar(double mom,double pathlen, double mass) const {
+    double retval(0.0);
+    static const double small(1e-3);
+    if(fabs(mass-e_mass_)<small){ // crude test to select electrons; should be replaced by a propoer particle enum TODO
+      retval = pow(mom*pathlen,2)*_eradvar;
+    }
+    return retval;
+  }
+
+
   //////////////////////////////////////////////////////////
 
   //// Calculate Moyal mean
-  double
-    DetMaterial::moyalMean(double deltap, double xi) const{
-      //getting most probable energy loss, or mpv:
-      double energylossmpv = fabs(deltap);
+  double DetMaterial::moyalMean(double deltap, double xi) const{
+    //getting most probable energy loss, or mpv:
+    double energylossmpv = fabs(deltap);
 
-      //declare moyalsigma for sanity check
-      double moyalsigma = xi;
+    //declare moyalsigma for sanity check
+    double moyalsigma = xi;
 
-      //forming the Moyal Mean
+    //forming the Moyal Mean
 
-      //note: when this is moved to c++20, the eulergamma constant should be replaced by 'egamma_v' in #include <numbers>
-      constexpr static double moyalmeanfactor = 0.57721566490153286 + M_LN2 ; //approximate Euler-Mascheroni (also known as gamma) constant (0.577...), see https://mathworld.wolfram.com/Euler-MascheroniConstant.html, added to log(2). This sum is used for the calculation of the closed-form Moyal mean below
-      double mmean = energylossmpv + moyalsigma * moyalmeanfactor; //formula from https://reference.wolfram.com/language/ref/MoyalDistribution.html, see end of file for more information
+    //note: when this is moved to c++20, the eulergamma constant should be replaced by 'egamma_v' in #include <numbers>
+    constexpr static double moyalmeanfactor = 0.57721566490153286 + M_LN2 ; //approximate Euler-Mascheroni (also known as gamma) constant (0.577...), see https://mathworld.wolfram.com/Euler-MascheroniConstant.html, added to log(2). This sum is used for the calculation of the closed-form Moyal mean below
+    double mmean = energylossmpv + moyalsigma * moyalmeanfactor; //formula from https://reference.wolfram.com/language/ref/MoyalDistribution.html, see end of file for more information
 
-      return -1.0 * mmean;
-    }
+    return -1.0 * mmean;
+  }
 
   ////Calculate density correction for energy loss
-  double
-    DetMaterial::densityCorrection(double bg2) const {
-      // density correction
-      double x = 0;
-      double delta = 0;
-      x = log(bg2)/twoln10 ;
-      if ( x < _x0 ) {
-        if(_delta0 > 0) {
-          delta = _delta0*pow(10.0,2*(x-_x0));
-        }
-        else {
-          delta = 0.;
-        }
-      } else {
-        delta = twoln10*x - _bigc;
-        if ( x < _x1 )
-          delta += _afactor * pow((_x1 - x), _mpower);
-      }
-      return delta;
-    }
-
-  //// Caluclate shell correction for energy loss
-  double
-    DetMaterial::shellCorrection(double bg2, double tau) const {
-      double sh = 0;
-      double x = 1;
-      // shell correction
-      if ( bg2 > bg2lim ) {
-        //sh = 0. ;
-        //x = 1. ;
-        for (int k=0; k<=2; k++) {
-          x *= bg2 ;
-          sh += (*_shellCorrectionVector)[k]/x;
-        }
+  double DetMaterial::densityCorrection(double bg2) const {
+    // density correction
+    double x = 0;
+    double delta = 0;
+    x = log(bg2)/twoln10 ;
+    if ( x < _x0 ) {
+      if(_delta0 > 0) {
+        delta = _delta0*pow(10.0,2*(x-_x0));
       }
       else {
-        //sh = 0. ;
-        //x = 1. ;
-        for (int k=0; k<2; k++) {
-          x *= bg2lim ;
-          sh += (*_shellCorrectionVector)[k]/x;
-        }
-        sh *= log(tau/_taul)/log(taulim/_taul);
+        delta = 0.;
       }
-      return sh;
+    } else {
+      delta = twoln10*x - _bigc;
+      if ( x < _x1 )
+        delta += _afactor * pow((_x1 - x), _mpower);
     }
+    return delta;
+  }
 
-  //below, the old BTrk model 'energyLoss' function based on dE/dx has been renamed (G3 for geant3)
-  //and now 'energyLoss' above refers to the new most probable energy loss method
-
-  double
-    DetMaterial::energyLossG3(double mom, double pathlen,double mass) const {
-      // make sure we take positive lengths!
-      pathlen = fabs(pathlen);
-      double dedx = dEdx(mom,_elossType,mass);
-      // see how far I can step, within tolerance, given this energy loss
-      double maxstep = maxStepdEdx(mom,mass,dedx);
-      // if this is larger than my path, I'm done
-      if(maxstep>pathlen){
-        return dedx*pathlen;
-      } else {
-        // subdivide the material
-        unsigned nstep = std::min(int(pathlen/maxstep) + 1,maxnstep);
-        double step = pathlen/nstep;
-        double energy = particleEnergy(mom,mass);
-        double deltae = step*dedx;
-        double newenergy(energy+deltae);
-        double eloss(deltae);
-        for(unsigned istep=0;istep<nstep-1;istep++){
-          if(newenergy>mass){
-            // compute the new dedx given the new momentum
-            double newmom = particleMomentum(newenergy,mass);
-            deltae = step*dEdx(newmom,_elossType,mass);
-            // compute the loss in this step
-            eloss += deltae;
-            newenergy += deltae;
-          } else {
-            // lost all kinetic energy; stop
-            eloss = mass-energy;
-            break;
-          }
-        }
-        return eloss;
-      }
-    }
-
-
-  double
-    DetMaterial::energyGain(double mom, double pathlen, double mass) const {
-      // make sure we take positive lengths!
-      pathlen = fabs(pathlen);
-      double dedx = dEdx(mom,_elossType,mass);
-      // see how far I can step, within tolerance, given this energy loss
-      double maxstep = maxStepdEdx(mom,mass,dedx);
-      // if this is larger than my path, I'm done
-      if(maxstep>pathlen){
-        return -dedx*pathlen;
-      } else {
-        // subdivide the material
-        unsigned nstep = std::min(int(pathlen/maxstep) + 1,maxnstep);
-        double step = pathlen/nstep;
-        double energy = particleEnergy(mom,mass);
-        double deltae = -step*dedx;
-        // move to the middle of the slice of material
-        double newenergy(energy+deltae);
-        double egain(deltae);
-        for(unsigned istep=0;istep<nstep-1;istep++){
-          // compute the new dedx given the new momentum
-          double newmom = particleMomentum(newenergy,mass);
-          double deltae = -step*dEdx(newmom,_elossType,mass);
-          egain += deltae;
-          newenergy += deltae;
-        }
-        return egain;
+  //// Caluclate shell correction for energy loss
+  double DetMaterial::shellCorrection(double bg2, double tau) const {
+    double sh = 0;
+    double x = 1;
+    // shell correction
+    if ( bg2 > bg2lim ) {
+      //sh = 0. ;
+      //x = 1. ;
+      for (int k=0; k<=2; k++) {
+        x *= bg2 ;
+        sh += (*_shellCorrectionVector)[k]/x;
       }
     }
-  //
-  // calculate the energy deposited in an absorber. That's similiar to
-  // energyLoss, but the delta electron correction in the Bethe Bloch is
-  // switched on, and there is no Bremsstrahlung
-  //
-  double
-    DetMaterial::energyDeposit(double mom, double pathlen, double mass) const {
-      double dedx = dEdx(mom,deposit,mass);
-      return dedx*fabs(pathlen);
-    }
-
-  /********************** end of New Routines **************************/
-
-
-  //this 'energyLossRMS' now refers to the closed-form Moyal distribution RMS, see end of file for more information
-
-  double
-    DetMaterial::energyLossRMS(double mom,double pathlen,double mass) const {
-      if(mom>0.0){
-        //taking positive lengths
-        pathlen = fabs(pathlen) ;
-
-        double beta = particleBeta(mom, mass) ;
-
-        double moyalsigma = eloss_xi(beta, pathlen);
-
-
-        //forming the Moyal RMS
-        constexpr static double pisqrt2 = M_PI/M_SQRT2 ; //constant that is used to calculate the Moyal closed-form RMS: pi/sqrt(2), approx.
-        double mrms = pisqrt2 * moyalsigma ; //from https://reference.wolfram.com/language/ref/MoyalDistribution.html
-
-        return mrms;
-
-      } else {
-        return 0.0 ;
+    else {
+      //sh = 0. ;
+      //x = 1. ;
+      for (int k=0; k<2; k++) {
+        x *= bg2lim ;
+        sh += (*_shellCorrectionVector)[k]/x;
       }
+      sh *= log(tau/_taul)/log(taulim/_taul);
     }
+    return sh;
+  }
 
-
-  //below, the old BTrk model 'energyLossRMS' function, which has been renamed (G3 for geant3)
-  //
-  //  RMS of energy loss.  This is a gaussian approximation, stolen from
-  //  Geant3 (see Phys332)
-  //
-  double
-    DetMaterial::energyLossRMSG3(double mom,double pathlen,double mass) const {
-      //      double beta = particleBeta(mom,mass);
-      //      double emax = eloss_emax(mom,mass);
-      //      double xi = eloss_xi(beta,fabs(pathlen));
-      //      double kappa = xi/emax;
-      //      double gam = sqrt(1.0-0.5*pow(beta,2));
-      //      // formula comes from GFLUCT.F in gphys dnb Jun 4 2004
-      //      //
-      //      // this formula seriously overestimates the rms when kappa<0.001
-      //      // This only really affects electrons
-      //      // as for heavier particles resolution effects already dominate when we get to
-      //      // this range.  I'll truncate
-      //      if(kappa < _minkappa)kappa = _minkappa;
-      //      double elossrms = xi*sqrt(gam/kappa);
-      //      //  cout << "beta = " << beta
-      //      //       << " emax = " << emax
-      //      //       << " xi = " << xi
-      //      //       << " kappa = " << kappa
-      //      //       << " elossrms = " << elossrms << endl;
-      //      // this value is way too big: scale it down for now but this function needs a rewrite FIXME!
-      //      return elossrms;
-      return 0.5*energyLossG3(mom,pathlen,mass);
-    }
-
-
-  //
-  //  Functions needed for energy loss calculation, see reference above
-  //
-  double
-    DetMaterial::eloss_emax(double mom,double mass){
-      double beta = particleBeta(mom,mass);
-      double gamma = particleGamma(mom,mass);
-      double mratio = e_mass_/mass;
-      double emax = 2*e_mass_*pow(beta*gamma,2)/
-        (1+2*gamma*mratio + pow(mratio,2));
-      if(mass <= e_mass_)
-        emax *= 0.5;
-      return emax;
-    }
-
-  double
-    DetMaterial::eloss_xi(double beta,double pathlen) const{
-      return _dgev*_za*_density*fabs(pathlen)/pow(beta,2);
-    }
-
-  void
-    DetMaterial::print(ostream& os) const {
-      os << "Material " << _name << endl;
+  double DetMaterial::ionizationEnergyLossRMS(double mom,double pathlen,double mass) const {
+    if(mom>0.0){
+      //taking positive lengths
+      pathlen = fabs(pathlen) ;
+      double beta = particleBeta(mom, mass) ;
+      double moyalsigma = eloss_xi(beta, pathlen);
+      //forming the Moyal RMS
+      constexpr static double pisqrt2 = M_PI/M_SQRT2 ; //constant that is used to calculate the Moyal closed-form RMS: pi/sqrt(2), approx.
+      double mrms = pisqrt2 * moyalsigma ; //from https://reference.wolfram.com/language/ref/MoyalDistribution.html
+      return mrms;
+    } else {
+      return 0.0 ;
     }
+  }
 
-  void
-    DetMaterial::printAll(ostream& os) const {
-      os << "Material " << _name << " has properties : " << endl
-        << "  Effective Z = " << _zeff << endl
-        << "  Effective A = " << _aeff << endl
-        << "  Density (g/cm^3) = " << _density*cm*cm*cm  << endl
-        << "  Radiation Length (g/cm^2) = " << _radthick*cm*cm << endl
-        << "  Interaction Length (g/cm^2) = " << _intLength << endl
-        //   << "  Mean Ionization energy (MeV) = " << _meanion << endl
-        << "  Mean Ionization energy (MeV) = " << _eexc << endl;
-    }
 
-  double
-    DetMaterial::maxStepdEdx(double mom,double mass,double dEdx,double tol) {
-      // compute betagamma at entrance
-      double betagamma = particleBetaGamma(mom,mass);
-      double energy = particleEnergy(mom,mass);
-      // basic calculation, based on constant dE/dx
-      if(dEdx<0.0){
-        double maxstep = -tol*energy/dEdx;
-        // Modify for steep rise at low momentum
-        if(betagamma<2.0)
-          maxstep *= pow(betagamma,2)/betapower;
-        return maxstep;
-      }
-      else
-        return 1.0e6; // large step
-    }
+  double DetMaterial::eloss_xi(double beta,double pathlen) const{
+    return _dgev*_za*_density*fabs(pathlen)/pow(beta,2);
+  }
 
+  void DetMaterial::print(ostream& os) const {
+    os << "Material " << _name << endl;
+  }
 
-  double
-    DetMaterial::nSingleScatter(double mom,double pathlen, double mass) const {
-      double beta = particleBeta(mom,mass);
-      return pathlen*_nbar/pow(beta,2);
-    }
+  void DetMaterial::printAll(ostream& os) const {
+    os << "Material " << _name << " has properties : " << endl
+      << "  Effective Z = " << _zeff << endl
+      << "  Effective A = " << _aeff << endl
+      << "  Density (g/mm^3) = " << _density*cm*cm*cm  << endl
+      << "  Radiation Length (g/mm^2) = " << _radthick*cm*cm << endl
+      << "  Interaction Length (g/mm^2) = " << _intLength << endl
+      << "  Mean Ionization energy (MeV) = " << _eexc << endl;
+  }
 
 
   // note the scaterring momentum here is hard-coded to the simplified Highland formula,
   // this should only be used as a reference, not as a real estimate of the scattering sigma!!!!!
-  double
-    DetMaterial::highlandSigma(double mom,double pathlen, double mass) const {
-      if(mom>0.0){
-        double radfrac = _invx0*fabs(pathlen);
-        return _msmom*sqrt(radfrac)/(mom*particleBeta(mom,mass));
-      } else
-        return 1.0;
-    }
+  double DetMaterial::highlandSigma(double mom,double pathlen, double mass) const {
+    if(mom>0.0){
+      double radfrac = _invx0*fabs(pathlen);
+      return 15*sqrt(radfrac)/(mom*particleBeta(mom,mass));
+    } else
+      return 1.0;
+  }
+
 
 
   //Information about the Moyal Distribution Approx.:
diff --git a/MatEnv/DetMaterial.hh b/MatEnv/DetMaterial.hh
index a334f8a5..33080aaf 100644
--- a/MatEnv/DetMaterial.hh
+++ b/MatEnv/DetMaterial.hh
@@ -26,13 +26,15 @@
 #include <memory>
 
 namespace MatEnv {
+  struct DetMaterialConfig;
   class DetMaterial{
     public:
-      enum dedxtype {loss=0,deposit};
+ //Energy Loss model: choose 'mpv' for the Most Probable Energy Loss, or 'moyalmean' for the mean calculated via the Moyal Distribution approximation, see end of file for more information
+      enum energylossmode {mpv=0, moyalmean};
       //
       //  Constructor
       // new style
-      DetMaterial(const char* detName, const MtrPropObj* detMtrProp);
+      DetMaterial(const char* detName, const MtrPropObj* detMtrProp, DetMaterialConfig const& dmconf);
 
       ~DetMaterial();
       //
@@ -46,43 +48,32 @@ namespace MatEnv {
       bool operator == (const DetMaterial& other) const {
         return _name == other._name; }
 
-      double dEdx(double mom,dedxtype type,double mass) const;
-
-      enum energylossmode {mpv=0, moyalmean};
-      energylossmode _elossmode;
-
-      void setEnergyLossMode(energylossmode elossmode) {_elossmode = elossmode;}
-
-      double energyLossG3(double mom, double pathlen, double mass) const; // this will be the old BTrk model dE/dx-based energy loss function from Geant3
-
-      double energyLossRMSG3(double mom,double pathlen,double mass) const; //this is the old BTrk model energy loss RMS approx. from Geant3
-
-      //below, 'energyLoss' and 'energyLossRMS' now refer to the MPV-based energy loss (not dE/dx) and closed-form Moyal calculations, see end of DetMaterial.cc for more information on the Moyal distribution and parameters
+      // total energy loss (ionization and radiation) and variance
       double energyLoss(double mom,double pathlen,double mass) const;
+      double energyLossVar(double mom,double pathlen,double mass) const;
+      double energyLossRMS(double mom,double pathlen,double mass) const {
+        return sqrt(energyLossVar(mom,pathlen,mass));
+      }
 
+      // ionization energy loss functions using closed-form Moyal calculations, see end of DetMaterial.cc for more information
+      double ionizationEnergyLoss(double mom,double pathlen,double mass) const;
       // most probable value of energy loss
-      double energyLossMPV(double mom,double pathlen,double mass) const;
-
-
-      double energyLossRMS(double mom,double pathlen,double mass) const;
-
-      double energyLossVar(double mom,double pathlen,double mass) const {
-        double elrms = energyLossRMS(mom,pathlen,mass);
+      double ionizationEnergyLossMPV(double mom,double pathlen,double mass) const;
+      double ionizationEnergyLossRMS(double mom,double pathlen,double mass) const;
+      double ionizationEnergyLossVar(double mom,double pathlen,double mass) const {
+        double elrms = ionizationEnergyLossRMS(mom,pathlen,mass);
         return elrms*elrms;
       }
-      double energyDeposit(double mom, double pathlen,double mass) const;
-      double energyGain(double mom,double pathlen, double mass) const;
-      double nSingleScatter(double mom,double pathlen, double mass) const;
-      // terms used in first-principles single scattering model
-      double aParam(double mom) const { return 2.66e-6*pow(_zeff,0.33333333333333)/mom; }
-      double bParam(double mom) const { return    0.14/(mom*pow(_aeff,0.33333333333333)); }
+      // radiation (brehmsstrahlung) energy loss calculation. This is relevant only for electrons
       //
-      // Single Gaussian approximation, used in Kalman filtering
-      double scatterAngleRMS(double mom,double pathlen,double mass) const;
-      double scatterAngleVar(double mom,double pathlen,double mass) const {
-        double sarms = scatterAngleRMS(mom,pathlen,mass);
-        return sarms*sarms;
+      double radiationEnergyLoss(double mom,double pathlen, double mass) const;
+      double radiationEnergyLossVar(double mom,double pathlen, double mass) const;
+      // Single Gaussian approximation
+      double scatterAngleVar(double mom,double pathlen,double mass) const;
+      double scatterAngleRMS(double mom,double pathlen,double mass) const {
+        return sqrt(scatterAngleVar(mom,pathlen,mass));
       }
+
       double highlandSigma(double mom,double pathlen, double mass) const;
 
       static double particleEnergy(double mom,double mass) {
@@ -109,31 +100,20 @@ namespace MatEnv {
       //
       //  functions used to compute energy loss
       //
-      static double eloss_emax(double mom,double mass) ;
       double eloss_xi(double beta,double pathlen) const;
       double densityCorrection(double bg2) const;
       double shellCorrection(double bg2, double tau) const;
       double moyalMean(double deltap, double xi) const;
-      double kappa(double mom,double pathlen,double mass) const {
-        return eloss_xi(particleBeta(mom,mass),pathlen)/eloss_emax(mom,mass);}
-      //
-      // return the maximum step one can make through this material
-      // for a given momentum and particle type without dE/dx changing
-      // by more than the given tolerance (fraction).  This is an _approximate_
-      // function, based on a crude model of dE/dx.
-      static double maxStepdEdx(double mom,double mass, double dEdx,double tol=0.05);
     protected:
       //
       //  Constants used in material calculations
       //
-      double _msmom; // constant in Highland scattering formula
-      static double _minkappa; // ionization randomization parameter
       static double _dgev; // energy characterizing energy loss
       static const double _alpha; // fine structure constant
+      energylossmode _elossmode;
       double _scatterfrac; // fraction of scattering distribution to include in RMS
-      double _cutOffEnergy; // cut on max energy loss
-      dedxtype _elossType;
-
+      double _erad; // bremsstrahlong fractional energy loss scale
+      double _eradvar; // bremsstrahlong fractional energy loss variance scale
       //
       //  Specific data for this material
       //
@@ -143,7 +123,6 @@ namespace MatEnv {
       double _aeff; // effective Z of our material
       double _radthick; // radiation thickness in g/cm**2
       double _intLength; // ineraction length from MatMtrObj in g/cm**2
-      double _meanion; // mean ionization energy loss
       double _eexc; // mean ionization energy loss for new e_loss routine
       double _x0; /*  The following specify parameters for energy loss. see
                       Sternheimer etal,'Atomic Data and
@@ -157,17 +136,10 @@ namespace MatEnv {
 
       int _noem;
       std::vector< double >* _shellCorrectionVector;
-      std::vector< double >* _vecNbOfAtomsPerVolume;
-      std::vector< double >* _vecTau0;
-      std::vector< double >* _vecAlow;
-      std::vector< double >* _vecBlow;
-      std::vector< double >* _vecClow;
-      std::vector< double >* _vecZ;
       double _taul;
 
       // cached values to speed calculations
       double _invx0;
-      double _nbar;
       double _chic2;
       double _chia2_1;
       double _chia2_2;
@@ -179,7 +151,6 @@ namespace MatEnv {
       double aeff() const { return _aeff;}
       double radiationLength()const {return _radthick;}
       double intLength()const {return _intLength;}
-      double meanIon()const {return _meanion;}
       double eexc() const { return _eexc; }
       double X0()const {return _x0;}
       double X1()const {return _x1;}
@@ -196,20 +167,19 @@ namespace MatEnv {
       void print(std::ostream& os) const;
       void printAll(std::ostream& os ) const;
 
-      // parameters used in ionization energy loss
-      static double energyLossScale() { return _dgev; }
-      static void setEnergyLossScale(double dgev) { _dgev = dgev; }
       // parameters used in ionization energy loss randomization
-      static double minKappa() { return _minkappa; }
-      static void setMinimumKappa(double minkappa) { _minkappa = minkappa; }
       // scattering parameter
       double scatterFraction() const { return _scatterfrac;}
-      void setScatterFraction(double scatterfrac) {_scatterfrac = scatterfrac;}
-      double cutOffEnergy() const { return _cutOffEnergy;}
-      void setCutOffEnergy(double cutOffEnergy) {_cutOffEnergy = cutOffEnergy; _elossType = deposit; }
-      void setDEDXtype(dedxtype elossType) { _elossType = elossType;}
       static constexpr double e_mass_ = 5.10998910E-01; // electron mass in MeVC^2
   };
+
+  struct DetMaterialConfig {
+    DetMaterial::energylossmode elossmode_ = DetMaterial::moyalmean; // ionization energy loss mode
+    double scatterfrac_solid_ = 0.995; // scattering angle fraction cutoff for computing RMS in solids
+    double scatterfrac_gas_ = 0.999; // scattering angle fraction cutoff for computing RMS in gases
+    double ebrehmsfrac_ = 0.04; // electron Brehmsstrahlung energy fraction cutoff
+  };
+
 }
 #endif
 
diff --git a/MatEnv/MatDBInfo.cc b/MatEnv/MatDBInfo.cc
index 6d5fd2db..421cf16c 100644
--- a/MatEnv/MatDBInfo.cc
+++ b/MatEnv/MatDBInfo.cc
@@ -20,10 +20,13 @@
 #include <map>
 namespace MatEnv {
 
-  MatDBInfo::MatDBInfo(FileFinderInterface const& interface, DetMaterial::energylossmode elossmode ) :
-    _genMatFactory(RecoMatFactory::getInstance(interface)), _elossmode(elossmode)
+  MatDBInfo::MatDBInfo(FileFinderInterface const& interface, DetMaterialConfig const& dmconf ) :
+    _genMatFactory(RecoMatFactory::getInstance(interface)), _dmconf(dmconf)
   {;}
 
+  MatDBInfo::MatDBInfo(FileFinderInterface const& interface) : MatDBInfo(interface, DetMaterialConfig()){}
+
+
   MatDBInfo::~MatDBInfo() {;}
 
   void
@@ -70,8 +73,7 @@ namespace MatEnv {
 
     genMtrProp = _genMatFactory->GetMtrProperties(db_name);
     if(genMtrProp != 0){
-      theMat = std::make_shared<DetMaterial> ( detMatName.c_str(), genMtrProp ) ;
-      theMat->setEnergyLossMode(_elossmode);
+      theMat = std::make_shared<DetMaterial> ( detMatName.c_str(), genMtrProp, _dmconf ) ;
       that()->_matList[ detMatName ] = theMat;
       return theMat;
     } else {
diff --git a/MatEnv/MatDBInfo.hh b/MatEnv/MatDBInfo.hh
index 78bb5b7c..5aa7720f 100644
--- a/MatEnv/MatDBInfo.hh
+++ b/MatEnv/MatDBInfo.hh
@@ -37,7 +37,8 @@ namespace MatEnv {
 
   class MatDBInfo : public MaterialInfo {
     public:
-      MatDBInfo(FileFinderInterface const& interface, DetMaterial::energylossmode elossmode);
+      MatDBInfo(FileFinderInterface const& interface);
+      MatDBInfo(FileFinderInterface const& interface, DetMaterialConfig const& dmconf);
       virtual ~MatDBInfo();
       //  Find the material, given the name
       const std::shared_ptr<DetMaterial> findDetMaterial( const std::string& matName ) const override;
@@ -57,7 +58,7 @@ namespace MatEnv {
       MatDBInfo* that() const {
         return const_cast<MatDBInfo*>(this);
       }
-      DetMaterial::energylossmode _elossmode;
+      DetMaterialConfig _dmconf;
   };
 
 }
diff --git a/Tests/MatEnv_unit.cc b/Tests/MatEnv_unit.cc
index 3df06520..30e66d75 100644
--- a/Tests/MatEnv_unit.cc
+++ b/Tests/MatEnv_unit.cc
@@ -85,7 +85,13 @@ int main(int argc, char **argv) {
   cout << "Test for particle " << pname  << " mass " << pmass << endl;
   cout << "Searching for material " << matname << endl;
   MatEnv::SimpleFileFinder sfinder;
-  MatDBInfo matdbinfo(sfinder,MatEnv::DetMaterial::moyalmean);
+  DetMaterialConfig dmconf;
+  dmconf.elossmode_ = MatEnv::DetMaterial::moyalmean;
+  dmconf.scatterfrac_solid_ = 0.995;
+  dmconf.scatterfrac_gas_ = 0.999;
+  dmconf.ebrehmsfrac_ = 0.04;
+
+  MatDBInfo matdbinfo(sfinder,dmconf);
   const std::shared_ptr<DetMaterial> dmat = matdbinfo.findDetMaterial(matname);
   if(dmat != 0){
     cout << "Found DetMaterial " << dmat->name() << endl;
diff --git a/Tests/ToyMC.hh b/Tests/ToyMC.hh
index 90b43694..d0891748 100644
--- a/Tests/ToyMC.hh
+++ b/Tests/ToyMC.hh
@@ -41,7 +41,7 @@ namespace KKTest {
       using PCA = PiecewiseClosestApproach<KTRAJ,SensorLine>;
       // create from aseed
       ToyMC(BFieldMap const& bfield, double mom, int icharge, double zrange, int iseed, unsigned nhits, bool simmat, bool scinthit, double ambigdoca ,double simmass) :
-        bfield_(bfield), matdb_(sfinder_,MatEnv::DetMaterial::moyalmean), // use the moyal based eloss model
+        bfield_(bfield), matdb_(sfinder_),
         mom_(mom), icharge_(icharge),
         tr_(iseed), nhits_(nhits), simmat_(simmat), scinthit_(scinthit), ambigdoca_(ambigdoca), simmass_(simmass),
         sprop_(0.8*CLHEP::c_light), sdrift_(0.065),
diff --git a/Trajectory/PiecewiseTrajectory.hh b/Trajectory/PiecewiseTrajectory.hh
index 4d9cb365..08236fec 100644
--- a/Trajectory/PiecewiseTrajectory.hh
+++ b/Trajectory/PiecewiseTrajectory.hh
@@ -83,9 +83,9 @@ namespace KinKal {
       pieces_.erase(jpiece.base(),pieces_.end());
     } else if(trange.begin() > pieces_.front()->range().end() || trange.end() < pieces_.back()->range().begin())
       throw std::invalid_argument("PiecewiseTrajectory::setRange; Invalid Range");
-    // update piece range
-    pieces_.front()->range().restrict(trange);
-    pieces_.back()->range().restrict(trange);
+    // update end piece range
+    pieces_.front()->range().extend(trange.begin(),TimeDir::backwards);
+    pieces_.back()->range().extend(trange.end(),TimeDir::forwards);
   }
 
   template <class KTRAJ> PiecewiseTrajectory<KTRAJ>::PiecewiseTrajectory(KTRAJ const& piece) : pieces_(1,std::make_shared<KTRAJ>(piece))