Add FWHM=A0+A1*sqrt(E) functional form.

Added this functional form primarily for interoperability with other spectroscopy software.

InterSpec/DetectorPeakResponse.h

@@ -123,6 +123,12 @@ class DetectorPeakResponse
      */
     kSqrtEnergyPlusInverse,
 
+    /**
+     FWHM = `A0 + A1*sqrt(E)`
+     
+     For use with other applications - do not recommend actually using.
+     */
+    kConstantPlusSqrtEnergy,
 
     kNumResolutionFnctForm
   };//enum ResolutionFnctForm
@@ -814,6 +820,9 @@ class DetectorPeakResponse
 
   /** On 20230916 updated from version 0 to 1, to account for `m_fixedGeometry` - will still write version 0 if
    `m_geomType == EffGeometryType::FarField`.
+   
+   On 20240410 updated from 1 to 2, to account for `ResolutionFnctForm::kConstantPlusSqrtEnergy` type of FWHM
+   being added.  However, will only write 2 if `m_resolutionForm == ResolutionFnctForm::kConstantPlusSqrtEnergy`.
    */
   static const int sm_xmlSerializationVersion;
 

InterSpec/MakeDrfChart.h

@@ -42,7 +42,7 @@ class MakeDrfChart : public Wt::Chart::WCartesianChart
   /** Fwhm equation type; equivalent to DetectorPeakResponse::ResolutionFnctForm
    but redifined here to not have to include that header.
    */
-  enum class FwhmCoefType{ Gadras, SqrtEnergyPlusInverse, SqrtEqn };
+  enum class FwhmCoefType{ Gadras, SqrtEnergyPlusInverse, SqrtEqn, ConstantPlusSqrtEnergy };
 
   /** Information from the calibration data for a specific peak.  Used to draw
      data points on the chart.

InterSpec/MakeDrfFit.h

@@ -72,6 +72,14 @@ namespace MakeDrfFit
                                        std::vector<float> &coeffs,
                                        std::vector<float> &coeff_uncerts );
 
+  /** Fits `FWHM = A_0 + A_1*sqrt(x)
+   
+   Throws exception if fit fails.
+   */
+  double fit_constant_plus_sqrt_fwhm_lls( const std::deque< std::shared_ptr<const PeakDef> > &peaks,
+                                  std::vector<float> &coeffs,
+                                  std::vector<float> &coeff_uncerts );
+
   double peak_width_chi2( double predicted_sigma, const PeakDef &peak );
 
 

InterSpec/RelActCalcAuto.h

@@ -195,11 +195,16 @@ enum class FwhmForm : int
    */
   Gadras,
 
-  /** The equation FWHM = sqrt( A_0 + A_1*Energy + A_2/Energy ).
+  /** The equation `FWHM = sqrt( A_0 + A_1*Energy + A_2/Energy )`.
    This is what FRAM uses.
    */
   SqrtEnergyPlusInverse,
 
+  /** The equation `FWHM = A_0 + A_1*sqrt(Energy)`
+   This is what some commercial applications use - but should probably only be used for cases of interfacing to the.
+   */
+  ConstantPlusSqrtEnergy,
+
   /** The linear polynomial: e.g. FWHM = sqrt(A_0 + A_1*1*0.001*Energy)
    The "2" refers to how many parameters are being fit
    */

data/common_drfs.tsv

@@ -100,3 +100,4 @@ BNC SAM-950 3x3 NaI		UrlEncoded	CREATED=1707759637&DIAM=71.426&EFFT=P&EFFX=10*15
 
 Verifinder-SN23N		UrlEncoded	CREATED=1707759245&DIAM=76.166&EFFT=P&EFFX=10*15*20*25*30*35*40*45*50*55*60*70*80*90*100*120*150*180*200*220*240*260*280*300*330*370*400*450*500*550*600*650*700*760*830*910*1E3*1100*1200*1300*1400*1500*1600*1700*1800*1900*2E3*2200*2400*2600*2800*3E3*3500*4E3*4500*5E3*6E3*7E3*8E3*1E4*12000*2E4*5E4*1E5&EFFY=5.17E-4*0.10185*0.37601*0.59765*0.73124*0.63704*0.69512*0.73959*0.77473*0.80303*0.82614*0.86084*0.88483*0.90169*0.91373*0.92888*0.93576*0.93156*0.92333*0.9114*0.89593*0.87718*0.85582*0.83253*0.79549*0.74445*0.70574*0.64453*0.58883*0.53871*0.49372*0.45319*0.42313*0.39435*0.36627*0.34006*0.31639*0.2953*0.27852*0.26482*0.2504*0.23641*0.2241*0.21322*0.20351*0.19472*0.18673*0.17287*0.16122*0.15137*0.14247*0.14091*0.12344*0.10955*0.098185*0.088659*0.073258*0.06071*0.050153*0.037279*0.0371*0.038762*0.047484*0.058621&FWHMC=-3*7*0.65&FWHMT=GAD&HASH=13629241631890916945&HIGHE=2850&LASTUSED=1707759556&LOWE=50&NAME=Verifinder-SN23N&ORIGIN=6&VER=1
 ORTEC RadEagle		UrlEncoded	CREATED=1707767169&DIAM=75.15&EFFT=P&EFFX=10*15*20*25*30*35*40*45*50*55*60*70*80*90*100*120*150*180*200*220*240*260*280*300*330*370*400*450*500*550*600*650*700*760*830*910*1E3*1100*1200*1300*1400*1500*1600*1700*1800*1900*2E3*2200*2400*2600*2800*3E3*3500*4E3*4500*5E3*6E3*7E3*8E3*1E4*12000*2E4*5E4*1E5&EFFY=0*0.001558*0.06278*0.23299*0.41235*0.43239*0.52298*0.59183*0.64475*0.68605*0.71877*0.76631*0.79812*0.82*0.83535*0.85398*0.85804*0.82965*0.79235*0.74592*0.69547*0.64476*0.59627*0.55123*0.49036*0.42333*0.38146*0.32493*0.28106*0.24632*0.21825*0.19525*0.17715*0.15972*0.14336*0.12838*0.11513*0.10343*0.094145*0.086634*0.080092*0.074216*0.069136*0.064625*0.060627*0.057051*0.053841*0.048342*0.043782*0.039963*0.036091*0.034246*0.026753*0.020708*0.015626*0.011197*0.005348*0.00514*0.004973*0.004769*0.004697*0.005048*0.007285*0.009815&FWHMC=7.76*7.41*0.739&FWHMT=GAD&HASH=5209734961709224472&HIGHE=2900&LASTUSED=1707767176&LOWE=25&NAME=RadEagle&ORIGIN=6&VER=1
+Mirion SPIR-Ace		UrlEncoded	CREATED%3d1712789273%26DIAM%3d34.19%26EFFT%3dP%26EFFX%3d10*15*20*25*30*35*40*45*50*55*60*70*80*90*100*120*150*180*200*220*240*260*280*300*330*370*400*450*500*550*600*650*700*760*830*910*1E3*1100*1200*1300*1400*1500*1600*1700*1800*1900*2E3*2200*2400*2600*2800*3E3*3500*4E3*4500*5E3*6E3*7E3*8E3*1E4*12000*2E4*5E4*1E5%26EFFY%3d0*1.01E-4*0.018943*0.11678*0.2565*0.29941*0.38426*0.45071*0.50252*0.54328*0.57572*0.62307*0.65487*0.67666*0.69161*0.70781*0.70747*0.68795*0.66617*0.639*0.60804*0.575*0.54153*0.5088*0.46264*0.40772*0.37247*0.32326*0.28403*0.2525*0.22691*0.20582*0.18823*0.17068*0.1539*0.13838*0.12445*0.11196*0.10175*0.093214*0.086085*0.079968*0.07461*0.06986*0.06562*0.061802*0.058317*0.052185*0.046928*0.042416*0.038*0.037239*0.02868*0.021665*0.015727*0.010421*0.004903*0.00463*0.004425*0.004186*0.004107*0.004489*0.006678*0.00896%26FWHMC%3d-2.02*7.87*0.599%26FWHMT%3dGAD%26HASH%3d3886773262405883703%26LASTUSED%3d1712789295%26NAME%3dSPIR-Ace-NaI%26ORIGIN%3d6%26VER%3d1

src/DetectorPeakResponse.cpp

@@ -66,7 +66,7 @@
 using namespace std;
 using SpecUtils::Measurement;
 
-const int DetectorPeakResponse::sm_xmlSerializationVersion = 1;
+const int DetectorPeakResponse::sm_xmlSerializationVersion = 2;
 
 namespace
 {
@@ -623,6 +623,8 @@ bool DetectorPeakResponse::hasResolutionInfo() const
       return (m_resolutionCoeffs.size() == 3);
     case kSqrtEnergyPlusInverse:
       return (m_resolutionCoeffs.size() == 3);
+    case kConstantPlusSqrtEnergy:
+      return (m_resolutionCoeffs.size() == 2);
     case kSqrtPolynomial:
       return (m_resolutionCoeffs.size() > 1);
     case kNumResolutionFnctForm:
@@ -1359,6 +1361,9 @@ std::string DetectorPeakResponse::toAppUrl() const
     case kSqrtEnergyPlusInverse:
       parts["FWHMT"] = "FRAM";
       break;
+    case kConstantPlusSqrtEnergy:
+      parts["FWHMT"] = "GENIE";
+      break;
     case ResolutionFnctForm::kSqrtPolynomial:
       parts["FWHMT"] = "SQRTPOLY";
       break;
@@ -1669,6 +1674,9 @@ void DetectorPeakResponse::fromAppUrl( std::string url_query )
     }else if( parts["FWHMT"] == "FRAM" )
     {
       resolutionForm = ResolutionFnctForm::kSqrtEnergyPlusInverse;
+    }else if( parts["FWHMT"] == "GENIE" )
+    {
+      resolutionForm = ResolutionFnctForm::kConstantPlusSqrtEnergy;
     }else if( parts["FWHMT"] == "SQRTPOLY" )
     {
       resolutionForm = ResolutionFnctForm::kSqrtPolynomial;
@@ -1685,22 +1693,29 @@ void DetectorPeakResponse::fromAppUrl( std::string url_query )
     if( resolutionUncerts.size() && (resolutionUncerts.size() != resolutionCoeffs.size()) )
        throw runtime_error( "fromAppUrl: FWHMC and FWHMU different lengths" );
 
+    bool invalidNumCoef = false;
     switch( resolutionForm )
     {
       case ResolutionFnctForm::kGadrasResolutionFcn:
       case ResolutionFnctForm::kSqrtEnergyPlusInverse:
-       if( resolutionCoeffs.size() != 3 )
-         throw runtime_error( "fromAppUrl: invalid number resolution coefs" );
+        invalidNumCoef = (resolutionCoeffs.size() != 3);
+
         break;
 
+      case ResolutionFnctForm::kConstantPlusSqrtEnergy:
+        invalidNumCoef = (resolutionCoeffs.size() != 2);
+         break;
+
       case ResolutionFnctForm::kSqrtPolynomial:
-        if( resolutionCoeffs.size() < 2 )
-          throw runtime_error( "fromAppUrl: not enough resolution coefs" );
+        invalidNumCoef = (resolutionCoeffs.size() < 2);
         break;
 
       case ResolutionFnctForm::kNumResolutionFnctForm:
         break;
     }//switch( resolutionForm )
+
+    if( invalidNumCoef )
+      throw runtime_error( "fromAppUrl: invalid number resolution coefs" );
   }//if( parts.count("FWHMT") )
 
   if( parts.count("ORIGIN") )
@@ -2234,6 +2249,11 @@ void DetectorPeakResponse::setFwhmCoefficients( const std::vector<float> &coefs,
       if( coefs.size() != 3 )
         throw runtime_error( "setFwhmCoefficients: sqrt(A0+A1*E+A2/E) equation must have three coefficients." );
       break;
+
+    case ResolutionFnctForm::kConstantPlusSqrtEnergy:
+      if( coefs.size() != 2 )
+        throw runtime_error( "setFwhmCoefficients: A0 + A1*sqrt(E) equation must have two coefficients." );
+      break;
 
     case ResolutionFnctForm::kNumResolutionFnctForm:
       if( !coefs.empty() )
@@ -2259,8 +2279,15 @@ void DetectorPeakResponse::toXml( ::rapidxml::xml_node<char> *parent,
   parent->append_node( base_node );
 
   // We will write XML version 0, if m_geomType is not far-field (this was only change between 0 and 1)
-  static_assert( sm_xmlSerializationVersion == 1, "Update DetectorPeakResponse sm_xmlSerializationVersion");
-  snprintf( buffer, sizeof(buffer), "%i", ((m_geomType != EffGeometryType::FarField) ? sm_xmlSerializationVersion : 0) );
+  static_assert( sm_xmlSerializationVersion == 2, "Update DetectorPeakResponse sm_xmlSerializationVersion");
+
+  if( m_resolutionForm == ResolutionFnctForm::kConstantPlusSqrtEnergy )
+  {
+    snprintf( buffer, sizeof(buffer), "%i", sm_xmlSerializationVersion );
+  }else
+  {
+    snprintf( buffer, sizeof(buffer), "%i", ((m_geomType != EffGeometryType::FarField) ? 1 : 0) );
+  }
 
   const char *value = doc->allocate_string( buffer );
   xml_attribute<char> *attr = doc->allocate_attribute( "version", value );
@@ -2308,10 +2335,11 @@ void DetectorPeakResponse::toXml( ::rapidxml::xml_node<char> *parent,
 
   switch( m_resolutionForm )
   {
-    case kGadrasResolutionFcn:   val = "GadrasResolutionFcn";   break;
-    case kSqrtEnergyPlusInverse: val = "SqrtEnergyPlusInverse"; break;
-    case kSqrtPolynomial:        val = "SqrtPolynomial";        break;
-    case kNumResolutionFnctForm: val = "Undefined";             break;
+    case kGadrasResolutionFcn:    val = "GadrasResolutionFcn";   break;
+    case kSqrtEnergyPlusInverse:  val = "SqrtEnergyPlusInverse"; break;
+    case kConstantPlusSqrtEnergy: val = "ConstantPlusSqrtEnergy"; break;
+    case kSqrtPolynomial:         val = "SqrtPolynomial";        break;
+    case kNumResolutionFnctForm:  val = "Undefined";             break;
   }//switch( m_resolutionForm )
 
   node = doc->allocate_node( node_element, "ResolutionForm", val );
@@ -2586,6 +2614,8 @@ void DetectorPeakResponse::fromXml( const ::rapidxml::xml_node<char> *parent )
     m_resolutionForm = kGadrasResolutionFcn;
   else if( compare(node->value(),node->value_size(),"SqrtEnergyPlusInverse",21,false) )
     m_resolutionForm = kSqrtEnergyPlusInverse;
+  else if( compare(node->value(),node->value_size(),"ConstantPlusSqrtEnergy",22,false) )
+    m_resolutionForm = kConstantPlusSqrtEnergy;
   else if( compare(node->value(),node->value_size(),"SqrtPolynomial",14,false) )
     m_resolutionForm = kSqrtPolynomial;
   else if( compare(node->value(),node->value_size(),"Undefined",9,false) )
@@ -3271,6 +3301,16 @@ float DetectorPeakResponse::peakResolutionFWHM( float energy,
       return sqrt(pars[0] + pars[1]*energy + pars[2]/energy);
     }//case kSqrtEnergyPlusInverse:
 
+    case kConstantPlusSqrtEnergy:
+    {
+      if( pars.size() != 2 )
+        throw std::runtime_error( "DetectorPeakResponse::peakResolutionSigma():"
+                                 " pars not defined" );
+      energy /= PhysicalUnits::keV;
+
+      return pars[0] + pars[1]*sqrt(energy);
+    }//case kConstantPlusSqrtEnergy:
+
     case kSqrtPolynomial:
     {
       if( pars.size() < 1 )

src/DrfSelect.cpp

@@ -4009,8 +4009,14 @@ std::shared_ptr<DetectorPeakResponse> DrfSelect::parseRelEffCsvFile( const std::
       {
         if( SpecUtils::icontains( line, "Full width half maximum (FWHM) follows equation" ) )
         {
-          const bool isSqrt = SpecUtils::icontains( line, "sqrt(" );  //Else contains "GadrasEqn"
-          const auto form = isSqrt ? DetectorPeakResponse::ResolutionFnctForm::kSqrtPolynomial : DetectorPeakResponse::ResolutionFnctForm::kGadrasResolutionFcn;
+          const bool isConstPlusSqrt = SpecUtils::icontains( line, "A0 + A1*sqrt" );
+          const bool isSqrt = !isConstPlusSqrt && SpecUtils::icontains( line, "sqrt(" );  //Else contains "GadrasEqn"
+          auto form = DetectorPeakResponse::ResolutionFnctForm::kGadrasResolutionFcn;
+          if( isConstPlusSqrt )
+            form = DetectorPeakResponse::ResolutionFnctForm::kConstantPlusSqrtEnergy;
+          else if( isSqrt )
+            form = DetectorPeakResponse::ResolutionFnctForm::kSqrtPolynomial;
+
           int nlinecheck = 0;
           while( SpecUtils::safe_get_line(csvfile, line, 2048)
                 && !SpecUtils::icontains(line, "Values")
@@ -4029,7 +4035,7 @@ std::shared_ptr<DetectorPeakResponse> DrfSelect::parseRelEffCsvFile( const std::
               for( size_t i = 1; i < fwhm_fields.size(); ++i )
                 coefs.push_back( stof(fwhm_fields[i]) );
               if( !coefs.empty() )
-                det->setFwhmCoefficients( coefs,form );
+                det->setFwhmCoefficients( coefs, form );
             }catch(...)
             {
             }

src/MakeDrf.cpp

@@ -1498,6 +1498,10 @@ MakeDrf::MakeDrf( InterSpec *viewer, MaterialDB *materialDB,
         m_fwhmEqnType->addItem( "sqrt(A + B*E + C/E)" );
         break;
 
+      case DetectorPeakResponse::kConstantPlusSqrtEnergy:
+        m_fwhmEqnType->addItem( "A + B*sqrt(E)" );
+        break;
+
       case DetectorPeakResponse::kSqrtPolynomial:
         m_fwhmEqnType->addItem( "Sqrt Power Series" );
         break;
@@ -2505,6 +2509,7 @@ void MakeDrf::handleFwhmTypeChanged()
   {
     case DetectorPeakResponse::kGadrasResolutionFcn:
     case DetectorPeakResponse::kSqrtEnergyPlusInverse:
+    case DetectorPeakResponse::kConstantPlusSqrtEnergy:
     case DetectorPeakResponse::kNumResolutionFnctForm:
       m_sqrtEqnOrder->hide();
       valid_fwhm_form = true;
@@ -2629,6 +2634,7 @@ void MakeDrf::fitFwhmEqn( std::vector< std::shared_ptr<const PeakDef> > peaks,
     case DetectorPeakResponse::kGadrasResolutionFcn:
     case DetectorPeakResponse::kNumResolutionFnctForm:
     case DetectorPeakResponse::kSqrtEnergyPlusInverse:
+    case DetectorPeakResponse::kConstantPlusSqrtEnergy:
       break;
 
     case DetectorPeakResponse::kSqrtPolynomial:
@@ -2709,6 +2715,11 @@ void MakeDrf::updateFwhmEqn( std::vector<float> coefs,
       valid_fcntform = true;
       break;
 
+    case DetectorPeakResponse::ResolutionFnctForm::kConstantPlusSqrtEnergy:
+      eqnType = MakeDrfChart::FwhmCoefType::ConstantPlusSqrtEnergy;
+      valid_fcntform = true;
+      break;
+
     case DetectorPeakResponse::ResolutionFnctForm::kSqrtPolynomial:
       eqnType = MakeDrfChart::FwhmCoefType::SqrtEqn;
       valid_fcntform = true;
@@ -3131,6 +3142,10 @@ shared_ptr<DetectorPeakResponse> MakeDrf::assembleDrf( const string &name, const
         drf->setFwhmCoefficients( m_fwhmCoefs, DetectorPeakResponse::ResolutionFnctForm::kSqrtEnergyPlusInverse );
         break;
 
+      case DetectorPeakResponse::ResolutionFnctForm::kConstantPlusSqrtEnergy:
+        drf->setFwhmCoefficients( m_fwhmCoefs, DetectorPeakResponse::ResolutionFnctForm::kConstantPlusSqrtEnergy );
+        break;
+
       case DetectorPeakResponse::ResolutionFnctForm::kSqrtPolynomial:
         drf->setFwhmCoefficients( m_fwhmCoefs, DetectorPeakResponse::ResolutionFnctForm::kSqrtPolynomial );
         break;
@@ -3361,6 +3376,18 @@ void MakeDrf::writeCsvSummary( std::ostream &out,
         break;
       }//case DetectorPeakResponse::kSqrtEnergyPlusInverse:
 
+      case DetectorPeakResponse::kConstantPlusSqrtEnergy:
+      {
+        out << "FWHM = A0 + A1*sqrt( energy )" << endline
+        << "# Energy in keV" << endline
+        << "# ";
+        for( size_t i = 0; i < fwhmCoefs.size(); ++i )
+          out << ",A" << i;
+        out << endline;
+
+        break;
+      }//case DetectorPeakResponse::kConstantPlusSqrtEnergy:
+
       case DetectorPeakResponse::kNumResolutionFnctForm:
         assert( 0 );
         break;
@@ -3614,6 +3641,18 @@ void MakeDrf::writeRefSheet( std::ostream &output, std::string drfname, std::str
         break;
       }//case DetectorPeakResponse::kSqrtEnergyPlusInverse:
 
+      case DetectorPeakResponse::kConstantPlusSqrtEnergy:
+      {
+        assert( m_fwhmCoefs.size() == 2 );
+        const float A0 = (m_fwhmCoefs.size()) > 0 ? m_fwhmCoefs[0] : 0.0f;
+        const float A1 = (m_fwhmCoefs.size()) > 1 ? m_fwhmCoefs[1] : 0.0f;
+
+        char buffer[128] = { '\0' };
+        snprintf( buffer, sizeof(buffer), "FWHM(keV) = %.4g + %.4g*sqrt(x)", A0, A1 );
+        fwhm_eqn = buffer;
+        break;
+      }//case DetectorPeakResponse::kConstantPlusSqrtEnergy:
+
       case DetectorPeakResponse::kSqrtPolynomial:
       {
         fwhm_eqn = "FWHM(keV) = sqrt(";
@@ -3749,6 +3788,10 @@ void MakeDrf::writeRefSheet( std::ostream &output, std::string drfname, std::str
       case DetectorPeakResponse::kSqrtEnergyPlusInverse:
         fwhmEqnType = MakeDrfChart::FwhmCoefType::SqrtEnergyPlusInverse;
         break;
+
+      case DetectorPeakResponse::kConstantPlusSqrtEnergy:
+        fwhmEqnType = MakeDrfChart::FwhmCoefType::ConstantPlusSqrtEnergy;
+        break;
 
       case DetectorPeakResponse::kSqrtPolynomial:
         fwhmEqnType = MakeDrfChart::FwhmCoefType::SqrtEqn;

src/MakeDrfChart.cpp

@@ -507,6 +507,7 @@ void MakeDrfChart::updateFwhmEquationToModel()
 
   if( ((m_fwhmEqnType==FwhmCoefType::Gadras) && m_fwhmCoefs.size() != 3)
      || ((m_fwhmEqnType==FwhmCoefType::SqrtEnergyPlusInverse) && m_fwhmCoefs.size() != 3)
+     || ((m_fwhmEqnType==FwhmCoefType::ConstantPlusSqrtEnergy) && m_fwhmCoefs.size() != 2)
      || (m_fwhmEqnType==FwhmCoefType::SqrtEqn && m_fwhmCoefs.size() < 1) )
   {
     if( m->data(0, sm_equation_fwhm_col).empty() )  //
@@ -528,6 +529,10 @@ void MakeDrfChart::updateFwhmEquationToModel()
     case FwhmCoefType::SqrtEnergyPlusInverse:
       eqnType = DetectorPeakResponse::ResolutionFnctForm::kSqrtEnergyPlusInverse;
       break;
+
+    case FwhmCoefType::ConstantPlusSqrtEnergy:
+      eqnType = DetectorPeakResponse::ResolutionFnctForm::kConstantPlusSqrtEnergy;
+      break;
 
     case FwhmCoefType::SqrtEqn:
       eqnType = DetectorPeakResponse::ResolutionFnctForm::kSqrtPolynomial;