From 247b63be73a2b4e0610003c67177d910667f7ebb Mon Sep 17 00:00:00 2001
From: David Nathan Brown <dave_brown@lbl.gov>
Date: Mon, 27 Oct 2025 11:01:57 -0700
Subject: [PATCH 01/12] cleanup

---
 MatEnv/DetMaterial.cc | 146 ------------------------------------------
 MatEnv/DetMaterial.hh |  12 ----
 2 files changed, 158 deletions(-)

diff --git a/MatEnv/DetMaterial.cc b/MatEnv/DetMaterial.cc
index d5a03fc2..4048d3a7 100644
--- a/MatEnv/DetMaterial.cc
+++ b/MatEnv/DetMaterial.cc
@@ -74,14 +74,6 @@ 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)));
@@ -101,12 +93,6 @@ namespace MatEnv {
   DetMaterial::~DetMaterial()
   {
     delete _shellCorrectionVector;
-    delete _vecNbOfAtomsPerVolume;
-    delete _vecTau0;
-    delete _vecAlow;
-    delete _vecBlow;
-    delete _vecClow;
-    delete _vecZ;
   }
 
   //
@@ -305,90 +291,7 @@ namespace MatEnv {
       return sh;
     }
 
-  //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;
-      }
-    }
-  //
-  // 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 {
@@ -412,38 +315,6 @@ namespace MatEnv {
       }
     }
 
-
-  //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
   //
@@ -481,23 +352,6 @@ namespace MatEnv {
         << "  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::nSingleScatter(double mom,double pathlen, double mass) const {
diff --git a/MatEnv/DetMaterial.hh b/MatEnv/DetMaterial.hh
index a334f8a5..477d7ef0 100644
--- a/MatEnv/DetMaterial.hh
+++ b/MatEnv/DetMaterial.hh
@@ -53,10 +53,6 @@ namespace MatEnv {
 
       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
       double energyLoss(double mom,double pathlen,double mass) const;
 
@@ -70,8 +66,6 @@ namespace MatEnv {
         double elrms = energyLossRMS(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; }
@@ -157,12 +151,6 @@ 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

From 00f38960db0f68e4e832b63b0c1471f075a41629 Mon Sep 17 00:00:00 2001
From: David Nathan Brown <dave_brown@lbl.gov>
Date: Mon, 27 Oct 2025 11:20:18 -0700
Subject: [PATCH 02/12] More cleanup

---
 MatEnv/DetMaterial.cc | 51 +------------------------------------------
 MatEnv/DetMaterial.hh | 10 ---------
 2 files changed, 1 insertion(+), 60 deletions(-)

diff --git a/MatEnv/DetMaterial.cc b/MatEnv/DetMaterial.cc
index 4048d3a7..adf9b898 100644
--- a/MatEnv/DetMaterial.cc
+++ b/MatEnv/DetMaterial.cc
@@ -30,10 +30,8 @@ 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 ;
@@ -52,7 +50,6 @@ namespace MatEnv {
     _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),
     _name(detMatName),
     _za(detMtrProp->getZ()/detMtrProp->getA()),
@@ -82,7 +79,6 @@ namespace MatEnv {
     _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) {
@@ -127,52 +123,8 @@ namespace MatEnv {
         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 ;
-
-        return dedx;
-
-      } else
-        return 0.0;
-    }
-
-  //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
+  //Most probable energy loss from https://pdg.lbl.gov/2019/reviews/rpp2018-rev-passage-particles-matter.pdf
   double DetMaterial::energyLoss(double mom, double pathlen, double mass) const {
     if(mom>0.0){
       double beta  = particleBeta(mom,mass) ;
@@ -199,7 +151,6 @@ namespace MatEnv {
       double j = 0.200 ;
       double tau = gamma-1;
 
-
       // most probable energy loss function
       beta2 = beta*beta ;
       gamma2 = gamma*gamma ;
diff --git a/MatEnv/DetMaterial.hh b/MatEnv/DetMaterial.hh
index 477d7ef0..b804637c 100644
--- a/MatEnv/DetMaterial.hh
+++ b/MatEnv/DetMaterial.hh
@@ -46,8 +46,6 @@ 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;
 
@@ -111,11 +109,6 @@ namespace MatEnv {
       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
@@ -125,7 +118,6 @@ namespace MatEnv {
       static double _dgev; // energy characterizing energy loss
       static const double _alpha; // fine structure constant
       double _scatterfrac; // fraction of scattering distribution to include in RMS
-      double _cutOffEnergy; // cut on max energy loss
       dedxtype _elossType;
 
       //
@@ -193,8 +185,6 @@ namespace MatEnv {
       // 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
   };

From bcb0914cc10456015dd13de9b3ac92cafb879a15 Mon Sep 17 00:00:00 2001
From: David Nathan Brown <dave_brown@lbl.gov>
Date: Mon, 27 Oct 2025 11:59:47 -0700
Subject: [PATCH 03/12] more cleanup

---
 MatEnv/DetMaterial.cc | 76 ++++++++++++++-----------------------------
 MatEnv/DetMaterial.hh | 27 +++------------
 2 files changed, 30 insertions(+), 73 deletions(-)

diff --git a/MatEnv/DetMaterial.cc b/MatEnv/DetMaterial.cc
index adf9b898..0bf111a0 100644
--- a/MatEnv/DetMaterial.cc
+++ b/MatEnv/DetMaterial.cc
@@ -31,13 +31,10 @@ using std::ostream;
 namespace MatEnv {
 
   double DetMaterial::_dgev = 0.153536e2;
-  double DetMaterial::_minkappa(1.0e-3);
   // 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);
 
@@ -45,19 +42,16 @@ namespace MatEnv {
   //  Constructor
   //
 
-  double cm(10.0); // temporary hack
+  double cm(10.0); // convert cm to mm
   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),
-    _elossType(loss),
     _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()),
@@ -73,14 +67,10 @@ namespace MatEnv {
       new std::vector< double >(detMtrProp->getShellCorrectionVector());
     // 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) {
-      _elossType = deposit;
-    }
     if (detMtrProp->getState() == "gas" && detMtrProp->getDensity()<0.01) {
       _scatterfrac = 0.999999;
     }
@@ -94,34 +84,26 @@ namespace MatEnv {
   //
   //  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::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
+  }
 
 
   //Most probable energy loss from https://pdg.lbl.gov/2019/reviews/rpp2018-rev-passage-particles-matter.pdf
@@ -296,28 +278,20 @@ namespace MatEnv {
       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
+        << "  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;
     }
 
 
-  double
-    DetMaterial::nSingleScatter(double mom,double pathlen, double mass) const {
-      double beta = particleBeta(mom,mass);
-      return pathlen*_nbar/pow(beta,2);
-    }
-
-
   // 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));
+        return 15*sqrt(radfrac)/(mom*particleBeta(mom,mass));
       } else
         return 1.0;
     }
diff --git a/MatEnv/DetMaterial.hh b/MatEnv/DetMaterial.hh
index b804637c..692d3df4 100644
--- a/MatEnv/DetMaterial.hh
+++ b/MatEnv/DetMaterial.hh
@@ -64,17 +64,13 @@ namespace MatEnv {
         double elrms = energyLossRMS(mom,pathlen,mass);
         return elrms*elrms;
       }
-      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)); }
       //
-      // 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;
+      // 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) {
@@ -113,12 +109,9 @@ namespace MatEnv {
       //
       //  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
       double _scatterfrac; // fraction of scattering distribution to include in RMS
-      dedxtype _elossType;
 
       //
       //  Specific data for this material
@@ -129,7 +122,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
@@ -147,7 +139,6 @@ namespace MatEnv {
 
       // cached values to speed calculations
       double _invx0;
-      double _nbar;
       double _chic2;
       double _chia2_1;
       double _chia2_2;
@@ -159,7 +150,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;}
@@ -176,16 +166,9 @@ 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;}
-      void setDEDXtype(dedxtype elossType) { _elossType = elossType;}
       static constexpr double e_mass_ = 5.10998910E-01; // electron mass in MeVC^2
   };
 }

From 5be71c6f32704a8688bf086e443039d04dc42ee9 Mon Sep 17 00:00:00 2001
From: David Nathan Brown <dave_brown@lbl.gov>
Date: Mon, 27 Oct 2025 12:47:53 -0700
Subject: [PATCH 04/12] More cleanup

---
 MatEnv/DetMaterial.cc |  6 +++---
 MatEnv/DetMaterial.hh | 16 ++++++++--------
 MatEnv/MatDBInfo.cc   |  3 +--
 3 files changed, 12 insertions(+), 13 deletions(-)

diff --git a/MatEnv/DetMaterial.cc b/MatEnv/DetMaterial.cc
index 0bf111a0..82219dbc 100644
--- a/MatEnv/DetMaterial.cc
+++ b/MatEnv/DetMaterial.cc
@@ -43,8 +43,8 @@ namespace MatEnv {
   //
 
   double cm(10.0); // convert cm to mm
-  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
+  DetMaterial::DetMaterial(const char* detMatName, const MtrPropObj* detMtrProp, energylossmode elossmode):
+    _elossmode(elossmode),
     _scatterfrac(0.9999),
     _name(detMatName),
     _za(detMtrProp->getZ()/detMtrProp->getA()),
@@ -106,7 +106,6 @@ namespace MatEnv {
   }
 
 
-  //Most probable energy loss from https://pdg.lbl.gov/2019/reviews/rpp2018-rev-passage-particles-matter.pdf
   double DetMaterial::energyLoss(double mom, double pathlen, double mass) const {
     if(mom>0.0){
       double beta  = particleBeta(mom,mass) ;
@@ -121,6 +120,7 @@ namespace MatEnv {
       return 0.0;
   }
 
+  //Most probable energy loss from https://pdg.lbl.gov/2019/reviews/rpp2018-rev-passage-particles-matter.pdf
   double DetMaterial::energyLossMPV(double mom, double pathlen, double mass) const {
     if(mom>0.0){
       //taking positive lengths
diff --git a/MatEnv/DetMaterial.hh b/MatEnv/DetMaterial.hh
index 692d3df4..84d10c8e 100644
--- a/MatEnv/DetMaterial.hh
+++ b/MatEnv/DetMaterial.hh
@@ -28,11 +28,12 @@
 namespace MatEnv {
   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, as well as discussion about radiative losses
+      enum energylossmode {mpv=0, moyalmean};
       //
       //  Constructor
       // new style
-      DetMaterial(const char* detName, const MtrPropObj* detMtrProp);
+      DetMaterial(const char* detName, const MtrPropObj* detMtrProp, energylossmode);
 
       ~DetMaterial();
       //
@@ -46,10 +47,10 @@ namespace MatEnv {
       bool operator == (const DetMaterial& other) const {
         return _name == other._name; }
 
-      enum energylossmode {mpv=0, moyalmean};
-      energylossmode _elossmode;
 
-      void setEnergyLossMode(energylossmode elossmode) {_elossmode = elossmode;}
+
+
+
 
       //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
       double energyLoss(double mom,double pathlen,double mass) const;
@@ -102,15 +103,14 @@ namespace MatEnv {
       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);}
-      //
+      double kappa(double mom,double pathlen,double mass) const { return eloss_xi(particleBeta(mom,mass),pathlen)/eloss_emax(mom,mass);}
     protected:
       //
       //  Constants used in material calculations
       //
       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
 
       //
diff --git a/MatEnv/MatDBInfo.cc b/MatEnv/MatDBInfo.cc
index 6d5fd2db..bc1e3cea 100644
--- a/MatEnv/MatDBInfo.cc
+++ b/MatEnv/MatDBInfo.cc
@@ -70,8 +70,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, _elossmode ) ;
       that()->_matList[ detMatName ] = theMat;
       return theMat;
     } else {

From fa80eff9ce650836890940ecb67febc94ec84e6d Mon Sep 17 00:00:00 2001
From: David Nathan Brown <dave_brown@lbl.gov>
Date: Mon, 27 Oct 2025 14:02:30 -0700
Subject: [PATCH 05/12] More cleanup. Add explicit radiation energy loss (dummy
 for now)

---
 MatEnv/DetMaterial.cc | 217 ++++++++++++++++++++----------------------
 MatEnv/DetMaterial.hh |  29 +++---
 2 files changed, 115 insertions(+), 131 deletions(-)

diff --git a/MatEnv/DetMaterial.cc b/MatEnv/DetMaterial.cc
index 82219dbc..59630394 100644
--- a/MatEnv/DetMaterial.cc
+++ b/MatEnv/DetMaterial.cc
@@ -105,12 +105,19 @@ namespace MatEnv {
       return 1.0; // 'infinite' scattering
   }
 
+  double DetMaterial::energyLoss(double mom,double pathlen,double mass) const {
+    return ionizationEnergyLoss(mom,pathlen,mass);
+  }
+
+  double DetMaterial::energyLossVar(double mom,double pathlen,double mass) const {
+    return ionizationEnergyLossVar(mom,pathlen,mass);
+  }
 
-  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);
@@ -121,7 +128,7 @@ namespace MatEnv {
   }
 
   //Most probable energy loss from https://pdg.lbl.gov/2019/reviews/rpp2018-rev-passage-particles-matter.pdf
-  double DetMaterial::energyLossMPV(double mom, double pathlen, double mass) const {
+  double DetMaterial::ionizationEnergyLossMPV(double mom, double pathlen, double mass) const {
     if(mom>0.0){
       //taking positive lengths
       pathlen = fabs(pathlen) ;
@@ -156,145 +163,123 @@ namespace MatEnv {
       return 0.0;
   }
 
+  double DetMaterial::radiationEnergyLoss(double mom,double pathlen, double mass) const {
+    return 0;0;
+  }
+  double DetMaterial::radiationEnergyLossVar(double mom,double pathlen, double mass) const {
+    return 0.0;
+  }
+
+
   //////////////////////////////////////////////////////////
 
   //// 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);
+  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.;
       }
-      return delta;
+    } 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;
-        }
-      }
-      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);
+  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;
       }
-      return sh;
     }
-
-
-
-  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;
+  }
 
-  //
-  //  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::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 ;
     }
+  }
 
-  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::eloss_xi(double beta,double pathlen) const{
+    return _dgev*_za*_density*fabs(pathlen)/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;
-    }
+  void DetMaterial::print(ostream& os) const {
+    os << "Material " << _name << endl;
+  }
+
+  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 15*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 84d10c8e..b7fb0f7d 100644
--- a/MatEnv/DetMaterial.hh
+++ b/MatEnv/DetMaterial.hh
@@ -47,25 +47,26 @@ namespace MatEnv {
       bool operator == (const DetMaterial& other) const {
         return _name == other._name; }
 
-
-
-
-
-
-      //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;
       }
+      // radiation (brehmsstrahlung) energy loss calculation. This is relevant only for electrons
       //
+      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 {
@@ -98,12 +99,10 @@ 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);}
     protected:
       //
       //  Constants used in material calculations

From 6acbf840472531e35bdd0b4c8b0cd15426e9b008 Mon Sep 17 00:00:00 2001
From: David Nathan Brown <dave_brown@lbl.gov>
Date: Thu, 30 Oct 2025 11:04:25 -0700
Subject: [PATCH 06/12] Partial implementaiton of radiation energy loss; not
 yet normalized

---
 MatEnv/DetMaterial.cc | 46 +++++++++++++++++++++++++++++++++++++------
 MatEnv/DetMaterial.hh |  5 ++++-
 2 files changed, 44 insertions(+), 7 deletions(-)

diff --git a/MatEnv/DetMaterial.cc b/MatEnv/DetMaterial.cc
index 59630394..ba97dfb9 100644
--- a/MatEnv/DetMaterial.cc
+++ b/MatEnv/DetMaterial.cc
@@ -45,7 +45,8 @@ namespace MatEnv {
   double cm(10.0); // convert cm to mm
   DetMaterial::DetMaterial(const char* detMatName, const MtrPropObj* detMtrProp, energylossmode elossmode):
     _elossmode(elossmode),
-    _scatterfrac(0.9999),
+    _scatterfrac(0.995), // should come from configuration TODO
+    _ymax(0.1), // should come from configuration TODO
     _name(detMatName),
     _za(detMtrProp->getZ()/detMtrProp->getA()),
     _zeff(detMtrProp->getZ()),
@@ -72,8 +73,12 @@ namespace MatEnv {
     _chia2_2 = 3.34*pow(_zeff*_alpha,2);
 
     if (detMtrProp->getState() == "gas" && detMtrProp->getDensity()<0.01) {
-      _scatterfrac = 0.999999;
+      _scatterfrac = 0.999; // should come from configuration TODO
     }
+    // calculate the brehms energy loss scale based on the cutoff fraction ymax = max(E_gamma/E_electron)
+    // The cutoff describes where the energy loss is so extreme as to prevent subsequent hits from being
+    // included in the track fit.
+    _e_lpm = 7.7e6*_radthick/_density;
   }
 
   DetMaterial::~DetMaterial()
@@ -106,11 +111,13 @@ namespace MatEnv {
   }
 
   double DetMaterial::energyLoss(double mom,double pathlen,double mass) const {
-    return ionizationEnergyLoss(mom,pathlen,mass);
+    double retval = ionizationEnergyLoss(mom,pathlen,mass); // + radiationEnergyLoss(mom,pathlen,mass);
+    return retval;
   }
 
   double DetMaterial::energyLossVar(double mom,double pathlen,double mass) const {
-    return ionizationEnergyLossVar(mom,pathlen,mass);
+    double retval = ionizationEnergyLossVar(mom,pathlen,mass); //+ radiationEnergyLossVar(mom,pathlen,mass);
+    return retval;
   }
 
   double DetMaterial::ionizationEnergyLoss(double mom, double pathlen, double mass) const {
@@ -163,11 +170,37 @@ namespace MatEnv {
       return 0.0;
   }
 
+  // radiation energy loss normalization
+  double DetMaterial::radiationNormalization(double mom) const {
+    double ymin = _e_lpm/mom; // electron mass is negligible
+    return 4*(log(_ymax/ymin) - _ymax)/3.0 + 0.5*(pow(_ymax,2));
+  }
+
   double DetMaterial::radiationEnergyLoss(double mom,double pathlen, double mass) const {
-    return 0;0;
+    // trancated average bremhsstrahlung radiation energy loss
+    // Based on integrating the 'low-y' approximation of the brehms cross-section from 0 to ymax, using
+    // formula 34.29 of the PDG 'Passage of particles through matter', Phys. Rev. D 110, 030001 (2024)
+    // lower cutoff to brehmsstrahlung due to coherence effect. from equation 34.33 in RPG
+    // calculation based on equations 33.29 and 34.33 of the RPG 'pass
+    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
+      // normalization for brehms calculations
+      double norm = radiationNormalization(mom);
+      retval = -(4.0*_ymax-2.0*pow(_ymax,2) + pow(_ymax,3))/(3*norm);
+    }
+    return retval;
   }
   double DetMaterial::radiationEnergyLossVar(double mom,double pathlen, double mass) const {
-    return 0.0;
+    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
+      double norm = radiationNormalization(mom);
+      double meany = radiationEnergyLoss(mom,pathlen,mass);
+      double meany2 = (2.0*pow(_ymax,2)/3.0 - 4.0*pow(_ymax,3)/9.0 + pow(_ymax,4)/4.0)/norm;
+      retval = meany2 - pow(meany,2);
+    }
+    return retval;
   }
 
 
@@ -282,6 +315,7 @@ namespace MatEnv {
   }
 
 
+
   //Information about the Moyal Distribution Approx.:
 
   //The Moyal distribution is an approximation for the ionization energy loss distribution. Unlike the Landau distribution is provides a closed-form energy loss mean and RMS. Code above uses the closed-form Moyal RMS for RMS, and allows the option of choosing the closed-form Moyal mean for the total energy loss parameter, which utilizes the most probable energy loss function. The options for either most probable energy loss and moyal distribution mean is toggled with the DetMaterial class member '_elossmode' with the options 'mpv' or 'moyalmean' respectively.
diff --git a/MatEnv/DetMaterial.hh b/MatEnv/DetMaterial.hh
index b7fb0f7d..0af1da5a 100644
--- a/MatEnv/DetMaterial.hh
+++ b/MatEnv/DetMaterial.hh
@@ -111,7 +111,8 @@ namespace MatEnv {
       static const double _alpha; // fine structure constant
       energylossmode _elossmode;
       double _scatterfrac; // fraction of scattering distribution to include in RMS
-
+      double _e_lpm; // brehmsstrahlung energy lower limit
+      double _ymax; // bremsstrahlong fractional energy upper limit
       //
       //  Specific data for this material
       //
@@ -141,6 +142,8 @@ namespace MatEnv {
       double _chic2;
       double _chia2_1;
       double _chia2_2;
+      // helper function
+      double radiationNormalization(double mom) const;
 
     public:
       // baseic accessors

From 3e08f734a8281d23943f62d92440e6b6ff983511 Mon Sep 17 00:00:00 2001
From: David Nathan Brown <dave_brown@lbl.gov>
Date: Thu, 30 Oct 2025 14:19:24 -0700
Subject: [PATCH 07/12] Add bounds check

---
 Fit/Track.hh | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/Fit/Track.hh b/Fit/Track.hh
index 2230912a..799b9964 100644
--- a/Fit/Track.hh
+++ b/Fit/Track.hh
@@ -694,7 +694,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());

From 6092f239a28b8614fab72cb8ba6271b4f8db33b2 Mon Sep 17 00:00:00 2001
From: David Nathan Brown <dave_brown@lbl.gov>
Date: Thu, 30 Oct 2025 16:56:27 -0700
Subject: [PATCH 08/12] Working version of radiation energy loss

---
 MatEnv/DetMaterial.cc | 37 +++++++++++++------------------------
 MatEnv/DetMaterial.hh |  6 ++----
 2 files changed, 15 insertions(+), 28 deletions(-)

diff --git a/MatEnv/DetMaterial.cc b/MatEnv/DetMaterial.cc
index ba97dfb9..4fa3985c 100644
--- a/MatEnv/DetMaterial.cc
+++ b/MatEnv/DetMaterial.cc
@@ -46,7 +46,6 @@ namespace MatEnv {
   DetMaterial::DetMaterial(const char* detMatName, const MtrPropObj* detMtrProp, energylossmode elossmode):
     _elossmode(elossmode),
     _scatterfrac(0.995), // should come from configuration TODO
-    _ymax(0.1), // should come from configuration TODO
     _name(detMatName),
     _za(detMtrProp->getZ()/detMtrProp->getA()),
     _zeff(detMtrProp->getZ()),
@@ -75,10 +74,11 @@ namespace MatEnv {
     if (detMtrProp->getState() == "gas" && detMtrProp->getDensity()<0.01) {
       _scatterfrac = 0.999; // should come from configuration TODO
     }
-    // calculate the brehms energy loss scale based on the cutoff fraction ymax = max(E_gamma/E_electron)
-    // The cutoff describes where the energy loss is so extreme as to prevent subsequent hits from being
-    // included in the track fit.
-    _e_lpm = 7.7e6*_radthick/_density;
+    // compute the sum radiated photons*energy for the given cuttoff
+    double ymax(0.04);  // should come from configuration TODO
+    // 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()
@@ -111,12 +111,14 @@ namespace MatEnv {
   }
 
   double DetMaterial::energyLoss(double mom,double pathlen,double mass) const {
-    double retval = ionizationEnergyLoss(mom,pathlen,mass); // + radiationEnergyLoss(mom,pathlen,mass);
+    double retval = ionizationEnergyLoss(mom,pathlen,mass);
+    retval += radiationEnergyLoss(mom,pathlen,mass);
     return retval;
   }
 
   double DetMaterial::energyLossVar(double mom,double pathlen,double mass) const {
-    double retval = ionizationEnergyLossVar(mom,pathlen,mass); //+ radiationEnergyLossVar(mom,pathlen,mass);
+    double retval = ionizationEnergyLossVar(mom,pathlen,mass);
+    retval += radiationEnergyLossVar(mom,pathlen,mass);
     return retval;
   }
 
@@ -170,24 +172,14 @@ namespace MatEnv {
       return 0.0;
   }
 
-  // radiation energy loss normalization
-  double DetMaterial::radiationNormalization(double mom) const {
-    double ymin = _e_lpm/mom; // electron mass is negligible
-    return 4*(log(_ymax/ymin) - _ymax)/3.0 + 0.5*(pow(_ymax,2));
-  }
-
   double DetMaterial::radiationEnergyLoss(double mom,double pathlen, double mass) const {
-    // trancated average bremhsstrahlung radiation energy loss
+    // truncated average bremhsstrahlung radiation energy loss
     // Based on integrating the 'low-y' approximation of the brehms cross-section from 0 to ymax, using
-    // formula 34.29 of the PDG 'Passage of particles through matter', Phys. Rev. D 110, 030001 (2024)
-    // lower cutoff to brehmsstrahlung due to coherence effect. from equation 34.33 in RPG
-    // calculation based on equations 33.29 and 34.33 of the RPG 'pass
+    // 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
-      // normalization for brehms calculations
-      double norm = radiationNormalization(mom);
-      retval = -(4.0*_ymax-2.0*pow(_ymax,2) + pow(_ymax,3))/(3*norm);
+      retval = -mom*pathlen*_erad;
     }
     return retval;
   }
@@ -195,10 +187,7 @@ namespace MatEnv {
     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
-      double norm = radiationNormalization(mom);
-      double meany = radiationEnergyLoss(mom,pathlen,mass);
-      double meany2 = (2.0*pow(_ymax,2)/3.0 - 4.0*pow(_ymax,3)/9.0 + pow(_ymax,4)/4.0)/norm;
-      retval = meany2 - pow(meany,2);
+      retval = pow(mom*pathlen,2)*_eradvar;
     }
     return retval;
   }
diff --git a/MatEnv/DetMaterial.hh b/MatEnv/DetMaterial.hh
index 0af1da5a..0a4127c9 100644
--- a/MatEnv/DetMaterial.hh
+++ b/MatEnv/DetMaterial.hh
@@ -111,8 +111,8 @@ namespace MatEnv {
       static const double _alpha; // fine structure constant
       energylossmode _elossmode;
       double _scatterfrac; // fraction of scattering distribution to include in RMS
-      double _e_lpm; // brehmsstrahlung energy lower limit
-      double _ymax; // bremsstrahlong fractional energy upper limit
+      double _erad; // bremsstrahlong fractional energy loss scale
+      double _eradvar; // bremsstrahlong fractional energy loss variance scale
       //
       //  Specific data for this material
       //
@@ -142,8 +142,6 @@ namespace MatEnv {
       double _chic2;
       double _chia2_1;
       double _chia2_2;
-      // helper function
-      double radiationNormalization(double mom) const;
 
     public:
       // baseic accessors

From 3c7405d495b57bee4d4afe732ff607d1206c4d69 Mon Sep 17 00:00:00 2001
From: David Nathan Brown <dave_brown@lbl.gov>
Date: Thu, 30 Oct 2025 22:43:46 -0700
Subject: [PATCH 09/12] Fix longstanding memory issue

---
 Fit/Status.hh |  1 +
 Fit/Track.hh  | 11 ++++++-----
 2 files changed, 7 insertions(+), 5 deletions(-)

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 799b9964..f819cd08 100644
--- a/Fit/Track.hh
+++ b/Fit/Track.hh
@@ -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){

From a9504a3e3b69ff9cbfd4430da9a813be751b94bd Mon Sep 17 00:00:00 2001
From: David Nathan Brown <dave_brown@lbl.gov>
Date: Fri, 31 Oct 2025 16:01:16 -0700
Subject: [PATCH 10/12] Fix rare range bug

---
 Fit/Track.hh                      | 20 ++++++++++++--------
 General/TimeRange.cc              | 19 +++++++++++++++++++
 General/TimeRange.hh              | 13 +++++--------
 Trajectory/PiecewiseTrajectory.hh |  6 +++---
 4 files changed, 39 insertions(+), 19 deletions(-)

diff --git a/Fit/Track.hh b/Fit/Track.hh
index f819cd08..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) :
@@ -303,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){
@@ -320,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>;
@@ -328,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());
@@ -341,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);
   }
 
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/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))

From 599387474839fac491848ef84992b78c6cb2a52f Mon Sep 17 00:00:00 2001
From: David Nathan Brown <dave_brown@lbl.gov>
Date: Fri, 31 Oct 2025 16:01:43 -0700
Subject: [PATCH 11/12] Move DetMaterial configuration parameters into a struct

---
 MatEnv/DetMaterial.cc | 11 ++++++-----
 MatEnv/DetMaterial.hh | 13 +++++++++++--
 MatEnv/MatDBInfo.cc   |  9 ++++++---
 MatEnv/MatDBInfo.hh   |  5 +++--
 Tests/MatEnv_unit.cc  |  8 +++++++-
 Tests/ToyMC.hh        |  2 +-
 6 files changed, 34 insertions(+), 14 deletions(-)

diff --git a/MatEnv/DetMaterial.cc b/MatEnv/DetMaterial.cc
index 4fa3985c..f11d3a0d 100644
--- a/MatEnv/DetMaterial.cc
+++ b/MatEnv/DetMaterial.cc
@@ -43,9 +43,8 @@ namespace MatEnv {
   //
 
   double cm(10.0); // convert cm to mm
-  DetMaterial::DetMaterial(const char* detMatName, const MtrPropObj* detMtrProp, energylossmode elossmode):
-    _elossmode(elossmode),
-    _scatterfrac(0.995), // should come from configuration TODO
+  DetMaterial::DetMaterial(const char* detMatName, const MtrPropObj* detMtrProp, DetMaterialConfig const& dmconf):
+    _elossmode(dmconf.elossmode_),
     _name(detMatName),
     _za(detMtrProp->getZ()/detMtrProp->getA()),
     _zeff(detMtrProp->getZ()),
@@ -72,10 +71,12 @@ namespace MatEnv {
     _chia2_2 = 3.34*pow(_zeff*_alpha,2);
 
     if (detMtrProp->getState() == "gas" && detMtrProp->getDensity()<0.01) {
-      _scatterfrac = 0.999; // should come from configuration TODO
+      _scatterfrac = dmconf.scatterfrac_gas_;
+    } else {
+      _scatterfrac = dmconf.scatterfrac_solid_;
     }
     // compute the sum radiated photons*energy for the given cuttoff
-    double ymax(0.04);  // should come from configuration TODO
+    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
diff --git a/MatEnv/DetMaterial.hh b/MatEnv/DetMaterial.hh
index 0a4127c9..153f5f8f 100644
--- a/MatEnv/DetMaterial.hh
+++ b/MatEnv/DetMaterial.hh
@@ -26,14 +26,15 @@
 #include <memory>
 
 namespace MatEnv {
+  class DetMaterialConfig;
   class DetMaterial{
     public:
- //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
+ //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, energylossmode);
+      DetMaterial(const char* detName, const MtrPropObj* detMtrProp, DetMaterialConfig const& dmconf);
 
       ~DetMaterial();
       //
@@ -171,6 +172,14 @@ namespace MatEnv {
       double scatterFraction() const { return _scatterfrac;}
       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 bc1e3cea..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,7 +73,7 @@ namespace MatEnv {
 
     genMtrProp = _genMatFactory->GetMtrProperties(db_name);
     if(genMtrProp != 0){
-      theMat = std::make_shared<DetMaterial> ( detMatName.c_str(), genMtrProp, _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),

From 9d56111d770f77f9fef00a9f74cc590a64bf71d0 Mon Sep 17 00:00:00 2001
From: David Nathan Brown <dave_brown@lbl.gov>
Date: Sun, 2 Nov 2025 21:29:19 -0800
Subject: [PATCH 12/12] Fix name

---
 MatEnv/DetMaterial.hh | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/MatEnv/DetMaterial.hh b/MatEnv/DetMaterial.hh
index 153f5f8f..33080aaf 100644
--- a/MatEnv/DetMaterial.hh
+++ b/MatEnv/DetMaterial.hh
@@ -26,7 +26,7 @@
 #include <memory>
 
 namespace MatEnv {
-  class DetMaterialConfig;
+  struct DetMaterialConfig;
   class DetMaterial{
     public:
  //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