Fix plotting of double-sided crystal ball distribution not plotting f…

…ull x-range of peak, for large skew.

Finding the coverage-limits of DSCB distribution may fail for really extended skews (i.e., low power).  This commit extends possible range to search for limits, but it still may fail, in which case the full ROI will now be used (as apposed to 25 sigma previously).
This was only a plotting issue - the full range of the ROI was being used to perform the fit - so no previous peak areas/fits should change, they will now just be displayed properly.

Note: the lower coverage limit of the DSCB function is found analytically, but upper limit is found numerically through bisection - finding the upper limit analytically should be implemented at some point.

InterSpec/PeakDists.h

@@ -406,7 +406,9 @@ namespace PeakDists
    
    The area below the retuned limit should be quite accurate, while above may only be good to within 0.02 of requested area (i.e. `0.02*p`).
    
-   Throws error on invalid input.
+   TODO: currently, the upper limit is found via iteration - should fix math so it can analytically be found, like lower limit
+   
+   Throws error on invalid input - or if limit cant be found (which can happen for really larger skews).
    */
   std::pair<double,double> double_sided_crystal_ball_coverage_limits( const double mean, const double sigma,
                                                                      const double left_skew,

src/PeakDef.cpp

@@ -2969,9 +2969,9 @@ std::string PeakDef::gaus_peaks_to_json(const std::vector<std::shared_ptr<const
             }catch( std::exception & )
             {
               // CB dist can have really long tail, causing the coverage limits to fail, because
-              //  of unreasonable values - we'll limit it arbitrarily to 25 sigma
-              vis_limits.first = std::max( p.mean() - 25*p.sigma(), p.lowerX() );
-              vis_limits.second = std::min( p.mean() + 5*p.sigma(), p.upperX() );
+              //  of unreasonable values - in this case we'll use the entire ROI.
+              vis_limits.first = p.lowerX();
+              vis_limits.second = p.upperX();
             }
             break;
           case ExpGaussExp:
@@ -2981,24 +2981,30 @@ std::string PeakDef::gaus_peaks_to_json(const std::vector<std::shared_ptr<const
                                                                   hidden_frac );
             break;
           case DoubleSidedCrystalBall:
-            if( p.coefficient(CoefficientType::SkewPar1) < 5.0
-               || p.coefficient(CoefficientType::SkewPar3) < 5.0)
-              hidden_frac = 1.0E-3; //DSCB doesnt behave well for power-law peaks...
+            // For largely skewed peaks, going out t 1E-6 is unreasonable, so we could limit
+            //  this to 1E-3, to improve chances of success.
+            //if( p.coefficient(CoefficientType::SkewPar1) < 5.0
+            //   || p.coefficient(CoefficientType::SkewPar3) < 5.0)
+            //  hidden_frac = 1.0E-3; //DSCB doesnt behave well for power-law peaks...
 
             try
             {
-              vis_limits = PeakDists::double_sided_crystal_ball_coverage_limits( p.mean(), p.sigma(),
-                                                                              p.coefficient(CoefficientType::SkewPar0),
-                                                                              p.coefficient(CoefficientType::SkewPar1),
-                                                                              p.coefficient(CoefficientType::SkewPar2),
-                                                                              p.coefficient(CoefficientType::SkewPar3),
-                                                                              hidden_frac );
+              const double mean = p.mean();
+              const double sigma = p.sigma();
+              const double left_skew = p.coefficient(CoefficientType::SkewPar0);
+              const double left_n = p.coefficient(CoefficientType::SkewPar1);
+              const double right_skew = p.coefficient(CoefficientType::SkewPar2);
+              const double right_n = p.coefficient(CoefficientType::SkewPar3);
+              const double p = hidden_frac;
+
+              vis_limits = PeakDists::double_sided_crystal_ball_coverage_limits( mean, sigma,
+                                              left_skew, left_n, right_skew, right_n, hidden_frac );
             }catch( std::exception & )
             {
               // DSCB dist can have really long tail, causing the coverage limits to fail, because
-              //  of unreasonable values - we'll limit it arbitrarily to 25 sigma
-              vis_limits.first = std::max( p.mean() - 25*p.sigma(), p.lowerX() );
-              vis_limits.second = std::min( p.mean() + 25*p.sigma(), p.upperX() );
+              //  of unreasonable values - in this case we'll use the entire ROI.
+              vis_limits.first = p.lowerX();
+              vis_limits.second = p.upperX();
             }//try / catch
 
             break;

src/PeakDists.cpp

@@ -685,7 +685,7 @@ namespace PeakDists
       return -0.5*p + (1.0 - indefinite_x( x ));
     };
 
-    if( (p <= 1.0E-11) || (p > 0.999)  ) // for gaussian, 7 sigma would be 1.279812544E-12
+    if( (p <= 1.0E-11) || (p > (1.0 - 1.0E-11))  ) // for gaussian, 7 sigma would be 1.279812544E-12
       throw runtime_error( "double_sided_crystal_ball_coverage_limits: invalid p" );
 
     auto x_from_eqn = [mean, sigma, norm, left_skew, left_n, right_skew, right_n,
@@ -726,27 +726,94 @@ namespace PeakDists
 
       // TODO: if value is in right tail, we will fail, in which case we will resort to an iterative solution
       double lower_x = -999, upper_x = -999;
-      try{ lower_x = x_from_eqn( 0.5*p ); }catch( std::exception & ){ cout << "Failed to find lower answer by eqn from p=" << p << endl; }
-      try{ upper_x = x_from_eqn( 1.0 - 0.5*p ); }catch( std::exception & ){ cout << "Failed to find upper answer by eqn from p=" << p << endl; }
+      try
+      {
+        lower_x = x_from_eqn( 0.5*p );
+      }catch( std::exception & )
+      {
+        cout << "Failed to find lower answer by eqn from p=" << p << endl;
+      }
+
+      try
+      {
+        upper_x = x_from_eqn( 1.0 - 0.5*p );
+      }catch( std::exception & )
+      {
+        cout << "Failed to find upper answer by eqn from p=" << p << endl;
+      }
 
 
-      if( (lower_x < 989) || IsNan(lower_x) || IsInf(lower_x) )
+      if( (lower_x < -998) || IsNan(lower_x) || IsInf(lower_x) )
       {
-        boost::uintmax_t max_iter = 100;
-        const double low_low_limit = mean - 200*sigma;  //200 sigma is arbitrary
-        const double low_up_limit = mean;
+        // With huge skew and/or small probabilities, the energy-value we are looking for can be
+        //  really far away from mean, so we will search by doubling the search range.
+        double low_low_limit;
+        bool found_lower = false;
+        const size_t max_range_doublings = 40;
+        double current_dx = 20.0 * std::max(sigma, 0.1); //20 sigma starting is arbitrary
+        for( size_t i = 0; !found_lower && (i < max_range_doublings); ++i, current_dx *= 2 )
+        {
+          low_low_limit = mean - current_dx;
+          const double y = lower_fcn( low_low_limit );
+          found_lower = (y < 0.0);
+        }//for( look for x where pdf has gone below limit )
+
+        // If we didnt find the lower limit, `boost::math::tools::bisect(...)` will throw exception
+        //assert( found_lower );
+
+#if( PERFORM_DEVELOPER_CHECKS )
+        if( !found_lower )
+        {
+          const string msg = "Failed to find lower limit by " + std::to_string(mean - current_dx)
+          + " for double_sided_crystal_ball_coverage_limits( "
+          + std::to_string(mean) + ", " + std::to_string(sigma) + ", " + std::to_string(left_skew)
+          + ", " + std::to_string(left_n) + ", " + std::to_string(right_skew) + ", "
+          + std::to_string(right_n) + ", " + std::to_string(p) + " )";
+
+          log_developer_error( __func__, msg.c_str() );
+        }//if( !found_upper )
+#endif
+
+        boost::uintmax_t max_iter = 1000;
         const pair<double,double> lower_val = boost::math::tools::bisect( lower_fcn, low_low_limit,
-                                                                         low_up_limit, term_condition, max_iter );
+                                                                         mean, term_condition, max_iter );
         lower_x = 0.5*(lower_val.first + lower_val.second);
       }//if( lower_x < 989 || IsNan(lower_x) || IsInf(lower_x) )
 
 
-      if( (upper_x < 989) || IsNan(upper_x) || IsInf(upper_x) )
+      if( (upper_x < -998) || IsNan(upper_x) || IsInf(upper_x) )
       {
-        boost::uintmax_t max_iter = 100;
-        const double up_low_limit = mean;
-        const double up_up_limit = mean + 200*sigma;
-        const pair<double,double> upper_val = boost::math::tools::bisect( upper_fcn, up_low_limit,
+        // With huge skew and/or small probabilities, the energy-value we are looking for can be
+        //  really huge, so we will search by doubling the search range.
+        double up_up_limit;
+        bool found_upper = false;
+        const size_t max_range_doublings = 40;
+        double current_dx = 20.0 * std::max(sigma, 0.1);
+        for( size_t i = 0; !found_upper && (i < max_range_doublings); ++i, current_dx *= 2 )
+        {
+          up_up_limit = mean + current_dx;
+          const double y = upper_fcn( up_up_limit );
+          found_upper = (y < 0.0);
+        }//for( look for x where pdf has gone below limit )
+
+        // If we didnt find the upper limit, `boost::math::tools::bisect(...)` will throw exception
+        //assert( found_upper );
+
+#if( PERFORM_DEVELOPER_CHECKS )
+        if( !found_upper )
+        {
+          const string msg = "Failed to find upper limit by " + std::to_string(mean + current_dx)
+          + " for double_sided_crystal_ball_coverage_limits( "
+          + std::to_string(mean) + ", " + std::to_string(sigma) + ", " + std::to_string(left_skew)
+          + ", " + std::to_string(left_n) + ", " + std::to_string(right_skew) + ", "
+          + std::to_string(right_n) + ", " + std::to_string(p) + " )";
+
+          log_developer_error( __func__, msg.c_str() );
+        }//if( !found_upper )
+#endif
+
+        boost::uintmax_t max_iter = 1000;
+        const pair<double,double> upper_val = boost::math::tools::bisect( upper_fcn, mean,
                                                                          up_up_limit, term_condition, max_iter );
         upper_x = 0.5*(upper_val.first + upper_val.second);
       }//if( upper_x < 989 || IsNan(upper_x) || IsInf(upper_x) )
@@ -759,7 +826,8 @@ namespace PeakDists
       return pair<double,double>( lower_x, upper_x );
     }catch( std::exception &e )
     {
-      throw runtime_error( "double_sided_crystal_ball_coverage_limits: failed to find limit: " + string(e.what()) );
+      const string excmsg = e.what();
+      throw runtime_error( "double_sided_crystal_ball_coverage_limits: failed to find limit: " + excmsg );
     }//try / catch
 
     assert( 0 );

target/testing/test_PeakDists.cpp

@@ -848,7 +848,7 @@ BOOST_AUTO_TEST_CASE( DoubleSidedCrystalBall )
 
     double prob = 0.000000573303;
     pair<double,double> limits = double_sided_crystal_ball_coverage_limits( mean, sigma, alpha_low,
-                                                                n_low, alpha_high, n_high, prob );
+                                                                           n_low, alpha_high, n_high, prob );
     double fraction = 1.0 - double_sided_crystal_ball_integral( mean, sigma, alpha_low, n_low,
                                                                alpha_high, n_high,
                                                                limits.first, limits.second );
@@ -857,35 +857,62 @@ BOOST_AUTO_TEST_CASE( DoubleSidedCrystalBall )
 
     prob = 1.0E-3;
     limits = double_sided_crystal_ball_coverage_limits( mean, sigma, alpha_low,
-                                                                n_low, alpha_high, n_high, prob );
+                                                       n_low, alpha_high, n_high, prob );
     fraction = 1.0 - double_sided_crystal_ball_integral( mean, sigma, alpha_low, n_low,
                                                         alpha_high, n_high,
                                                         limits.first, limits.second );
     BOOST_CHECK_CLOSE( fraction, prob, 0.01 );
 
     prob = 1.0E-2;
     limits = double_sided_crystal_ball_coverage_limits( mean, sigma, alpha_low,
-                                                                n_low, alpha_high, n_high, prob );
+                                                       n_low, alpha_high, n_high, prob );
     fraction = 1.0 - double_sided_crystal_ball_integral( mean, sigma, alpha_low, n_low,
                                                         alpha_high, n_high,
                                                         limits.first, limits.second );
     BOOST_CHECK_CLOSE( fraction, prob, 0.01 );
 
     prob = 1.0E-1;
     limits = double_sided_crystal_ball_coverage_limits( mean, sigma, alpha_low,
-                                                                n_low, alpha_high, n_high, prob );
+                                                       n_low, alpha_high, n_high, prob );
     fraction = 1.0 - double_sided_crystal_ball_integral( mean, sigma, alpha_low, n_low,
                                                         alpha_high, n_high,
                                                         limits.first, limits.second );
     BOOST_CHECK_CLOSE( fraction, prob, 0.01 );
 
     prob = 0.317310507863;
     limits = double_sided_crystal_ball_coverage_limits( mean, sigma, alpha_low,
-                                                                n_low, alpha_high, n_high, prob );
+                                                       n_low, alpha_high, n_high, prob );
     fraction = 1.0 - double_sided_crystal_ball_integral( mean, sigma, alpha_low, n_low,
                                                         alpha_high, n_high,
                                                         limits.first, limits.second );
     BOOST_CHECK_CLOSE( fraction, prob, 0.1 );
+
+
+    mean = 200;
+    sigma = 1.25;
+    alpha_low = 1.0;
+    n_low = 2.3;
+    alpha_high = 2.0;
+    n_high = 1.5;
+    prob = 1.0E-3;
+
+    limits = double_sided_crystal_ball_coverage_limits( mean, sigma, alpha_low,
+                                                       n_low, alpha_high, n_high, prob );
+    //cout << "Limits are {" << limits.first << ", " << limits.second << "}" << endl; //{-215.032, 14154.5}
+
+
+    /*
+     //This next case will fail.
+    mean = 5577.6396484375;
+    sigma = 1.2235283851623535;
+    alpha_low = 0.79716122150421142;
+    n_low = 5.5732345581054688;
+    alpha_high = 1.9554067850112915;
+    n_high = 1.0500000715255737;
+    prob = 1.0E-6;
+    limits = double_sided_crystal_ball_coverage_limits( mean, sigma, alpha_low,
+                                                       n_low, alpha_high, n_high, prob );
+     */
   }