Implement bethemean energy loss option. Enable optional per-material …

MatEnv/DetMaterial.cc

@@ -44,7 +44,9 @@ namespace MatEnv {
 
   double cm(10.0); // convert cm to mm
   DetMaterial::DetMaterial(const char* detMatName, const MtrPropObj* detMtrProp, DetMaterialConfig const& dmconf):
-    _elossmode(dmconf.elossmode_),
+    // per-material eloss mode if the material specifies one (>=0), else the global default
+    _elossmode(detMtrProp->getElossMode() >= 0 ?
+        (energylossmode)detMtrProp->getElossMode() : dmconf.elossmode_),
     _name(detMatName),
     _za(detMtrProp->getZ()/detMtrProp->getA()),
     _zeff(detMtrProp->getZ()),
@@ -131,12 +133,45 @@ namespace MatEnv {
       //if using mean calculated from the Moyal Dist. Approx: (see end of file for more information)
       if(_elossmode == moyalmean) {
         return moyalMean(deltap, xi);
+      } else if(_elossmode == bethemean) {
+        //unrestricted Bethe-Bloch mean (full loss including energy carried off by delta-rays), see end of file
+        return ionizationEnergyLossBetheMean(mom, pathlen, mass);
       } else
         return deltap;
     } else
       return 0.0;
   }
 
+  //Unrestricted Bethe-Bloch MEAN ionization energy loss (PDG RPP eq. 34.5), NOT the most-probable value.
+  //With xi == (K/2)(Z/A)(rho x / beta^2) as used throughout this class, the standard mean stopping power
+  //  -dE/dx = K (Z/A)(1/beta^2) [ 0.5 ln(2 me c^2 beta^2 gamma^2 Tmax / I^2) - beta^2 - delta/2 - C/Z ]
+  //multiplied by the traversed grammage rho*x becomes
+  //  <dE> = xi [ ln(2 me c^2 beta^2 gamma^2 / I) + ln(Tmax / I) - 2 beta^2 - delta - 2 sh ],
+  //where Tmax is the maximum single-collision energy transfer, delta the density-effect correction and sh the
+  //shell correction (the same terms ionizationEnergyLossMPV uses). Unlike the MPV this is additive in path
+  //length, so integrating it over sub-steps of a decelerating track yields an unbiased mean.
+  double DetMaterial::ionizationEnergyLossBetheMean(double mom, double pathlen, double mass) const {
+    if(mom>0.0){
+      //taking positive lengths
+      pathlen = fabs(pathlen) ;
+      double beta  = particleBeta(mom,mass) ;
+      double gamma = particleGamma(mom,mass) ;
+      double beta2 = beta*beta ;
+      double bg2   = beta2*gamma*gamma ;                 // (beta*gamma)^2
+      double tau   = gamma - 1 ;
+      double xi    = eloss_xi(beta, pathlen);            // (K/2)(Z/A)(rho x / beta^2)
+      // maximum kinetic energy transferable to a free electron in a single collision (PDG RPP eq. 34.4)
+      double Tmax  = 2.*e_mass_*bg2 / (1. + 2.*gamma*e_mass_/mass + (e_mass_/mass)*(e_mass_/mass)) ;
+      // density-effect and shell corrections (identical to ionizationEnergyLossMPV)
+      double delta = densityCorrection(bg2);
+      double sh    = shellCorrection(bg2, tau);
+      double meanloss = log(2.*e_mass_*bg2/_eexc) + log(Tmax/_eexc) - 2.*beta2 - delta - 2.*sh ;
+      meanloss *= -xi ; // sign convention: energy loss is returned as a NEGATIVE energy change
+      return meanloss;
+    } else
+      return 0.0;
+  }
+
   //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){
@@ -308,7 +343,8 @@ 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.
+  //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 energy-loss parameter returned by ionizationEnergyLoss is toggled with the DetMaterial class member '_elossmode': 'mpv' (most probable value), 'moyalmean' (closed-form Moyal mean of the restricted loss), or 'bethemean'.
+  //'bethemean' returns the unrestricted Bethe-Bloch MEAN (ionizationEnergyLossBetheMean): the full mean energy loss including the energy carried away by energetic delta-rays. The Moyal/MPV values track only the locally-deposited (restricted) loss, whose dE/dx saturates at the Fermi plateau; the unrestricted mean keeps rising (the relativistic rise), so for a thick absorber it is the correct estimate of the momentum the particle actually loses. It is additive in path length (integrable over sub-steps), whereas the MPV is not. See PDG RPP 'Passage of particles through matter' eqs. 34.4-34.5.
   //reference for Moyal dist.: Theory of Ionization Fluctuation by J. E. Moyal, Phil. Mag. 46 (1955) 263
   //more useful references: https://reference.wolfram.com/language/ref/MoyalDistribution.html, http://www.stat.rice.edu/~dobelman/textfiles/DistributionsHandbook.pdf, and https://arxiv.org/pdf/1702.06655.pdf
 

MatEnv/DetMaterial.hh

@@ -29,8 +29,10 @@ namespace MatEnv {
   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
-      enum energylossmode {mpv=0, moyalmean};
+ //Energy Loss model: choose 'mpv' for the Most Probable Energy Loss, 'moyalmean' for the (restricted) mean
+ //calculated via the Moyal Distribution approximation, or 'bethemean' for the unrestricted Bethe-Bloch mean
+ //(the full mean loss including energy carried off by delta-rays); see end of file for more information
+      enum energylossmode {mpv=0, moyalmean, bethemean};
       //
       //  Constructor
       // new style
@@ -59,6 +61,10 @@ namespace MatEnv {
       double ionizationEnergyLoss(double mom,double pathlen,double mass) const;
       // most probable value of energy loss
       double ionizationEnergyLossMPV(double mom,double pathlen,double mass) const;
+      // unrestricted Bethe-Bloch MEAN ionization energy loss (the full mean, including the energy carried off by
+      // energetic delta-rays; PDG RPP eq. 34.5). Unlike ionizationEnergyLossMPV it is additive in path length, so
+      // integrating it over sub-steps of a decelerating track gives an unbiased mean. Selected by elossmode==bethemean.
+      double ionizationEnergyLossBetheMean(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);
@@ -167,6 +173,9 @@ namespace MatEnv {
       void print(std::ostream& os) const;
       void printAll(std::ostream& os ) const;
 
+      // ionization energy loss mode in effect for this material (per-material if the
+      // material specified one, otherwise the global DetMaterialConfig default)
+      energylossmode elossMode() const { return _elossmode; }
       // parameters used in ionization energy loss randomization
       // scattering parameter
       double scatterFraction() const { return _scatterfrac;}

MatEnv/MatMaterialList.cc

@@ -24,6 +24,8 @@
 #include <stdlib.h>
 #include <assert.h>
 #include <algorithm>
+#include <cctype>
+#include <stdexcept>
 #include "KinKal/MatEnv/BbrCollectionUtils.hh"
 
 //----------------------
@@ -96,7 +98,6 @@ namespace MatEnv {
     double refindex = 0.;
     double temperature = 0.;
     double pressure = 0.;
-    double Tcut = 0.;
     materials >> name;
     while( !materials.eof())
     {
@@ -146,8 +147,31 @@ namespace MatEnv {
 	matObj->setTemperature(temperature);
 	matObj->setPressure(pressure);
 	matObj->setState(state);
-	if (iss>>Tcut) {
-	  matObj->setTcut(Tcut);
+	// Optional trailing columns (order-independent): a numeric token is the Tcut
+	// (restricted energy-loss cut, as before), a keyword token selects the ionization
+	// energy loss mode for this material ("mpv"/"moyalmean"/"bethemean", matching
+	// DetMaterial::energylossmode). Anything else is ignored.
+	std::string token;
+	while (iss >> token) {
+	  char* endptr = nullptr;
+	  double tokval = strtod(token.c_str(), &endptr);
+	  if (endptr != token.c_str() && *endptr == '\0') {
+	    matObj->setTcut(tokval);
+	  } else {
+	    std::string key = token;
+	    std::transform(key.begin(), key.end(), key.begin(),
+		[](unsigned char c){ return std::tolower(c); });
+	    if (key == "mpv") {
+	      matObj->setElossMode(0);
+	    } else if (key == "moyalmean") {
+	      matObj->setElossMode(1);
+	    } else if (key == "bethemean") {
+	      matObj->setElossMode(2);
+	    } else {
+	      throw std::invalid_argument("MatMaterialList: unrecognized trailing token '" + token
+		+ "' for material " + name);
+	    }
+	  }
 	}
 
 	_vector.push_back(matObj);
@@ -226,6 +250,10 @@ namespace MatEnv {
       out <<" "<< radlen <<" "<< intlen <<" "<< refindex <<" "<< temperature 
 	<<" "<< pressure <<" "<< state;
       if (energyCut>0.0) { out <<" "<< energyCut; }
+      int elossMode = matObj->getElossMode();
+      if (elossMode==0) { out <<" mpv"; }
+      else if (elossMode==1) { out <<" moyalmean"; }
+      else if (elossMode==2) { out <<" bethemean"; }
 
       out << std::endl;
     }

MatEnv/MatMaterialObj.cc

@@ -48,7 +48,8 @@ namespace MatEnv {
     _matTemperature(0),
     _matPressure(0),
     _matState(" "),
-    _matTcut(0)
+    _matTcut(0),
+    _elossMode(-1)
   {
   }
 
@@ -107,6 +108,7 @@ namespace MatEnv {
     _matPressure = matcp.getPressure();
     _matState = matcp.getState();
     _matTcut = matcp.getTcut();
+    _elossMode = matcp.getElossMode();
   }
 
   MatMaterialObj& MatMaterialObj::operator= (const MatMaterialObj& matrl)
@@ -129,8 +131,9 @@ namespace MatEnv {
     _refIndex = matrl.getRefIndex();
     _matTemperature = matrl.getTemperature();
     _matPressure = matrl.getPressure();
-    _matState = matrl.getState();            
+    _matState = matrl.getState();
     _matTcut = matrl.getTcut();
+    _elossMode = matrl.getElossMode();
 
     return *this;
   }
@@ -172,6 +175,7 @@ namespace MatEnv {
       << " RefIndex: " << getRefIndex()  << "  Temperature: " << getTemperature() 
       << "  Pressure: " << getPressure() << "  State: " << getState();
     if (_matTcut>0.0) { cout<< " Tcut: "<<getTcut(); }
+    if (_elossMode>=0) { cout<< " ElossMode: "<<getElossMode(); }
     cout << endl;
   }
 

MatEnv/MatMaterialObj.hh

@@ -78,8 +78,11 @@ namespace MatEnv {
       double getRefIndex() const { return _refIndex; };  
       double getTemperature() const { return _matTemperature; };           
       double getPressure() const { return _matPressure; };  
-      std::string getState() const { return std::string(_matState); };              
+      std::string getState() const { return std::string(_matState); };
       double getTcut() const { return _matTcut; };
+      // ionization energy loss mode for this material (matches DetMaterial::energylossmode:
+      // 0=mpv, 1=moyalmean, 2=bethemean); -1 means unspecified -> fall back to the global default
+      int getElossMode() const { return _elossMode; };
 
       void setName(const std::string Name) {_matName=Name; };
       void setDensity(double Density) {_matDensity=Density; };
@@ -98,6 +101,7 @@ namespace MatEnv {
       void setPressure(double Pressure) {_matPressure=Pressure; };
       void setState(const std::string State) {_matState=State; };
       void setTcut(double Tcut) {_matTcut=Tcut; };
+      void setElossMode(int ElossMode) {_elossMode=ElossMode; };
 
 
     private:
@@ -123,6 +127,7 @@ namespace MatEnv {
       double _matPressure;          // Material Pressure
       std::string _matState;        // Material State (Gas, liquid, ...)
       double _matTcut;              // Maximum energy transfer allowed per interaction step (this means that for the dE/dx is used the restricted energy loss rate parameterization instead of the Bethe-Bloch)
+      int _elossMode = -1;          // per-material ionization energy loss mode (DetMaterial::energylossmode); -1 = unspecified (use global default)
 
       friend bool testCdb(const MatMaterialObj*, const MatMaterialObj*);
   };

MatEnv/MaterialsList.data

@@ -19,6 +19,9 @@ HDPE      0.96        0.   0.   -2      4  Hydrogen 0    2 Carbon 0    -10 -20 -
 HeCF4  0.0005331 0. 0. 2 0.334 Helium 0 0.666 CF4 1 -10 -20 -30 20.0 1.0 gas
 Kapton       1.4300    0.   0.   -4    37    Carbon     0 6 Oxygen 0 2 Nitrogen 0 6 Hydrogen 0    -10     -20     -30     20.0    1.0    solid
 Mylar       1.4    0.   0.   -3    10    Carbon     0 4 Oxygen 0 8 Hydrogen 0    -10     -20     -30     20.0    1.0    solid
+# Concrete (MARS composition, matches Offline/TrackerConditions/data/MaterialsList.data). Tagged
+# bethemean: a thick passive absorber where the unrestricted Bethe-Bloch mean is the correct estimator.
+CONCRETE_MARS 2.35 0.0 0.0 +9 0.006 Hydrogen 0 0.030 Carbon 0 0.500 Oxygen 0 0.010 Sodium 0 0.030 Aluminum 0 0.200 Silicon 0 0.010 Potasium 0 0.200 Calcium 0 0.014 Iron 0 -10 -20 -30 20.0 1.0 solid bethemean
 # specific mu2e materials
 straw-wire    19.300000    0.   0.   +1    100.0e-2    Tungsten   0    -10 -20 -30     20.0    1.0    solid
 # 80:20 (by volume) Argon:CO2

MatEnv/MtrPropObj.cc

@@ -57,6 +57,7 @@ namespace MatEnv {
     _dEdxFactor(0),
     _meanExciEnergy(0),
     _energyTcut(0),
+    _elossMode(-1),
     _zeff( 0 ),
     _aeff( 0 ),
     _temp( 0 ),
@@ -94,6 +95,7 @@ namespace MatEnv {
     _radLength( theMaterial->getRadLength() ),
     _intLength( theMaterial->getIntLength() ),
     _energyTcut( theMaterial->getTcut() ),
+    _elossMode( theMaterial->getElossMode() ),
     _zeff( theMaterial->getZeff() ),
     _aeff( theMaterial->getAeff() ),
     _temp( theMaterial->getTemperature()+STP_Temperature ),
@@ -190,6 +192,7 @@ namespace MatEnv {
     _radLength = matcp.getRadLength();
     _intLength = matcp.getIntLength();
     _energyTcut = matcp.getEnergyTcut();
+    _elossMode = matcp.getElossMode();
 
   }
 
@@ -212,6 +215,7 @@ namespace MatEnv {
       _radLength = matrl.getRadLength();
       _intLength = matrl.getIntLength();
       _energyTcut = matrl.getEnergyTcut();
+      _elossMode = matrl.getElossMode();
       *_state = matrl.getState();
 
       return *this;

MatEnv/MtrPropObj.hh

@@ -92,6 +92,8 @@ namespace MatEnv {
       double getMeanExciEnergy() const { return _meanExciEnergy; }
       const std::vector< double >& getShellCorrectionVector() const;
       double getEnergyTcut() const { return _energyTcut; }
+      // per-material ionization energy loss mode (DetMaterial::energylossmode); -1 = unspecified
+      int getElossMode() const { return _elossMode; }
 
       double getRadLength() const { return _radLength; }
       double getIntLength() const { return _intLength; }
@@ -140,6 +142,7 @@ namespace MatEnv {
       double _meanExciEnergy;
       std::vector< double >* _shellCorrectionVector;
       double _energyTcut;
+      int _elossMode = -1;          // per-material ionization energy loss mode; -1 = unspecified (use global default)
 
       double _zeff;
       double _aeff;

Tests/MatEnv_unit.cc

@@ -95,6 +95,48 @@ int main(int argc, char **argv) {
   const std::shared_ptr<DetMaterial> dmat = matdbinfo.findDetMaterial(matname);
   if(dmat != 0){
     cout << "Found DetMaterial " << dmat->name() << endl;
+
+    // ---- issue #221: per-material unrestricted (Bethe) energy-loss mean testbench ----
+    // Demonstrates (a) per-material eloss selection (CONCRETE_MARS is tagged 'bethemean' in the
+    // data file, overriding the global moyalmean default; an untagged material falls back to the
+    // global default), and (b) that for a thick slab the unrestricted Bethe mean is the largest
+    // loss estimator (relativistic rise), matching PDG RPP eqs. 34.4-34.5. Validate the table by
+    // hand: <dE> = xi*[ln(2 me bg2/I) + ln(Tmax/I) - 2 beta^2 - delta - 2 shell], xi=(K/2)(Z/A)(rho x/beta^2).
+    auto modeName = [](DetMaterial::energylossmode m)->const char*{
+      switch(m){
+        case DetMaterial::mpv: return "mpv";
+        case DetMaterial::moyalmean: return "moyalmean";
+        case DetMaterial::bethemean: return "bethemean";
+        default: return "unknown";
+      } };
+    cout << "\n==== per-material ionization energy loss (issue #221) ====" << endl;
+    cout << "global default eloss mode (dmconf): " << modeName(dmconf.elossmode_) << endl;
+    cout << "requested material '" << dmat->name() << "' resolved eloss mode: "
+         << modeName(dmat->elossMode()) << endl;
+    const std::shared_ptr<DetMaterial> dmc = matdbinfo.findDetMaterial("CONCRETE_MARS");
+    const std::shared_ptr<DetMaterial> dmu = matdbinfo.findDetMaterial("Target"); // untagged Al
+    if(dmc) cout << "  CONCRETE_MARS (tagged bethemean) -> elossMode() = " << modeName(dmc->elossMode()) << endl;
+    if(dmu) cout << "  Target (untagged Al)             -> elossMode() = " << modeName(dmu->elossMode())
+                 << " (== global default, proving the fallback)" << endl;
+    if(dmc){
+      double slab = 447.0; // mm, ExtShield crossbar concrete thickness
+      cout << "\nthick concrete slab " << slab << " mm, particle " << pname
+           << " (mass " << pmass << " MeV) -- energy loss in MeV (more negative = larger loss):" << endl;
+      printf("  %10s %12s %12s %12s %14s\n","p[MeV/c]","MPV","moyalmean","bethemean","dispatched");
+      double testmom[6] = {100.,300.,1000.,3000.,10000.,30000.};
+      for(double tp : testmom){
+        double beta = DetMaterial::particleBeta(tp,pmass);
+        double xi   = dmc->eloss_xi(beta, slab);
+        double mpvv = dmc->ionizationEnergyLossMPV(tp, slab, pmass);
+        double moyv = dmc->moyalMean(mpvv, xi);
+        double betv = dmc->ionizationEnergyLossBetheMean(tp, slab, pmass);
+        double disp = dmc->ionizationEnergyLoss(tp, slab, pmass); // routes to bethemean via the tag
+        printf("  %10.0f %12.4f %12.4f %12.4f %14.4f\n", tp, mpvv, moyv, betv, disp);
+      }
+    }
+    cout << "==========================================================\n" << endl;
+    // ---- end issue #221 testbench ----
+
     unsigned nstep(100);
     double momstep = (momend-momstart)/(nstep-1);
     TGraph* geloss = new TGraph(nstep);