Merge branch 'master' of github.com:sandialabs/InterSpec

InterSpec/BatchInfoLog.h

@@ -178,6 +178,8 @@ namespace BatchInfoLog
                        const std::string &filename,
                        const BatchPeak::BatchPeakFitResult * const peak_fit );
 
+  void add_energy_cal_json( nlohmann::basic_json<> &data,
+                           const std::shared_ptr<const SpecUtils::EnergyCalibration> &cal );
 
   void add_activity_fit_options_to_json( nlohmann::basic_json<> &data,
                                         const BatchActivity::BatchActivityFitOptions &options );

InterSpec/BatchPeak.h

@@ -138,6 +138,12 @@ namespace BatchPeak
 
     bool success;
     std::vector<std::string> warnings;
+
+    /** The original energy calibration of the spectrum, before re-fitting it (if done). */
+    std::shared_ptr<const SpecUtils::EnergyCalibration> original_energy_cal;
+
+    /** The energy calibration after fitting for it - will only be non-null if energy calibration was performed */
+    std::shared_ptr<const SpecUtils::EnergyCalibration> refit_energy_cal;
   };//struct BatchPeakFitResult
 
 

InterSpec_resources/static_text/ShieldSourceFitLog/std_peak_fit_log.tmplt.txt

@@ -1,5 +1,7 @@
 {% if exists("foreground") %}
+Foreground StartTime: {{ foreground.StartTime }}
 Foreground LiveTime: {% if foreground.LiveTime_s %}{{ foreground.LiveTime_s }} seconds{% endif %} {% if foreground.LiveTime %}({{ foreground.LiveTime }}){% endif %}
+Foreground RealTime: {% if foreground.RealTime_s %}{{ foreground.RealTime_s }} seconds{% endif %} {% if foreground.RealTime %}({{ foreground.RealTime }}){% endif %}
 {% endif %}
 
 {% if HasWarnings %}

InterSpec_resources/static_text/ShieldSourceFitLog/std_peak_fit_summary.tmplt.csv

@@ -9,11 +9,11 @@ Setup
     {% if exists("ExemplarSampleNumbers") %}Exemplar samples: {{ ExemplarSampleNumbers }}{% endif %}
 
 
-Filename, Mean, Amplitude, FWHM, LowerEnergy, UpperEnergy, LowerChannel, UpperChannel, ContinuumIndex, SourceName, SourceType, SourceEnergy, PeakColor, UserLabel, Mean Uncert, Amplitude Uncert, FWHM Uncert, Chi2Dof
+Filename, Mean, Amplitude, FWHM, LowerEnergy, UpperEnergy, LowerChannel, UpperChannel, ContinuumIndex, ContinuumType, SourceName, SourceType, SourceEnergy, PeakColor, UserLabel, Mean Uncert, Amplitude Uncert, FWHM Uncert, Chi2Dof, StartTime, LiveTime (s), RealTime(s),
 ## for file in Files
 ## if existsIn(file,"FitPeaks")
 ## for peak in file.FitPeaks.Peaks
-{{ file.Filename }}, {{ printFixed(peak.PeakMean,2) }}, {{ printCompact(peak.PeakAmplitude,6) }}, {{ printFixed(peak.PeakFwhm,2) }}, {{ printFixed(peak.LowerEnergy,2) }}, {{ printFixed(peak.UpperEnergy,2) }}, {% if peak.HasChannelRange %}{{ peak.LowerChannelInt }}{% endif %}, {% if peak.HasChannelRange %}{{ peak.UpperChannelInt }}{% endif %}, {{ peak.ContinuumIndex }}, {{ at(file.FitPeaks.Continua,peak.ContinuumIndex).ContinuumType }}, {{ peak.SourceName }}, {{ peak.SourceType }}, {{ peak.SourceEnergy }}, {{ peak.PeakColor }}, {{ peak.PeakUserLabel }}, {{ printFixed(peak.PeakMeanUncert,3) }}, {{ printFixed(peak.PeakAmplitudeUncert,1) }}, {{ printFixed(peak.PeakFwhmUncert,3) }}, {% if peak.HasChi2Dof %}{{ printFixed(peak.Chi2Dof, 4) }}{% endif %}
+{{ file.Filename }}, {{ printFixed(peak.PeakMean,2) }}, {{ printCompact(peak.PeakAmplitude,6) }}, {{ printFixed(peak.PeakFwhm,2) }}, {{ printFixed(peak.LowerEnergy,2) }}, {{ printFixed(peak.UpperEnergy,2) }}, {% if peak.HasChannelRange %}{{ peak.LowerChannelInt }}{% endif %}, {% if peak.HasChannelRange %}{{ peak.UpperChannelInt }}{% endif %}, {{ peak.ContinuumIndex }}, {{ at(file.FitPeaks.Continua,peak.ContinuumIndex).ContinuumType }}, {{ peak.SourceName }}, {{ peak.SourceType }}, {{ peak.SourceEnergy }}, {{ peak.PeakColor }}, {{ peak.PeakUserLabel }}, {{ printFixed(peak.PeakMeanUncert,3) }}, {{ printFixed(peak.PeakAmplitudeUncert,1) }}, {{ printFixed(peak.PeakFwhmUncert,3) }}, {% if peak.HasChi2Dof %}{{ printFixed(peak.Chi2Dof, 4) }}{% endif %}, {{file.foreground.StartTime}}, {{file.foreground.LiveTime_s}}, {{file.foreground.RealTime_s}}
 ## endfor
 ## endif
 ## endfor

cmake/FetchInterSpecDeps.txt

@@ -170,7 +170,6 @@ if( USE_REL_ACT_TOOL )
     GIT_REPOSITORY https://gitlab.com/libeigen/eigen.git
     #GIT_TAG        3147391d946bb4b6c68edd901f2add6ac1f31f8c # release-3.4.0
     GIT_TAG        9df21dc8b4b576a7aa5c0094daa8d7e8b8be60f0 #Updated 3.4 release, to pickup some CMake fixes
-    GIT_SHALLOW TRUE
   )
 
   FetchContent_MakeAvailable( eigen )

src/BatchInfoLog.cpp

@@ -1331,55 +1331,70 @@ void add_basic_src_details( const GammaInteractionCalc::SourceDetails &src,
     spec_obj["DerivedIsBackground"] = static_cast<bool>(spec.derived_data_properties()
                                                         & static_cast<uint32_t>(SpecUtils::Measurement::DerivedDataProperties::IsBackground));
 
-    auto cal = spec.energy_calibration();
-    auto &cal_obj = spec_obj["EnergyCal"];
+    add_energy_cal_json( spec_obj["EnergyCal"], spec.energy_calibration() );
 
+    if( spec_file )
+    {
+      spec_obj["InstrumentModel"] = spec_file->instrument_model();
+      spec_obj["SerialNumber"] = spec_file->instrument_id();
+      spec_obj["Manufacturer"] = spec_file->manufacturer();
+      spec_obj["InstrumentType"] = spec_file->instrument_type();
+      spec_obj["DetectorType"] = SpecUtils::detectorTypeToString( spec_file->detector_type() );
+      spec_obj["NumberRecordsInFile"] = static_cast<int>( spec_file->num_measurements() );
+      spec_obj["RemarksInFile"] = spec_file->remarks();
+      spec_obj["ParseWarningsForFile"] = spec_file->parse_warnings();
+
+      spec_obj["HasInstrumentRid"] = !!spec_file->detectors_analysis();
+      if( spec_file->detectors_analysis() )
+        spec_obj["InstrumentRidSummary"] = riidAnaSummary( spec_file );
+    }//if( spec_file )
+  }//add_hist_to_json(...)
+
+
+  void add_energy_cal_json( nlohmann::basic_json<> &cal_obj,
+                           const std::shared_ptr<const SpecUtils::EnergyCalibration> &cal )
+  {
     cal_obj["NumChannels"] = cal ? static_cast<int>(cal->num_channels()) : 0;
     const SpecUtils::EnergyCalType cal_type = cal ? cal->type() : SpecUtils::EnergyCalType::InvalidEquationType;
     switch( cal_type )
     {
       case SpecUtils::EnergyCalType::Polynomial:
       case SpecUtils::EnergyCalType::UnspecifiedUsingDefaultPolynomial:
         cal_obj["Type"] = "Polynomial";
-        cal_obj["Coefficients"] = vector<double>{ begin(cal->coefficients()), end(cal->coefficients()) };
         break;
 
       case SpecUtils::EnergyCalType::FullRangeFraction:
         cal_obj["Type"] = "FullRangeFraction";
-        cal_obj["Coefficients"] = vector<double>{ begin(cal->coefficients()), end(cal->coefficients()) };
         break;
+
       case SpecUtils::EnergyCalType::LowerChannelEdge:
         cal_obj["Type"] = "LowerChannelEdge";
         break;
+
       case SpecUtils::EnergyCalType::InvalidEquationType:
         cal_obj["Type"] = "Invalid";
         break;
     }//switch( cal->type() )
 
-    if( cal && !cal->deviation_pairs().empty() )
-    {
-      vector<vector<double>> pairs;
-      for( const pair<float,float> &p : cal->deviation_pairs() )
-        pairs.push_back( { static_cast<double>(p.first), static_cast<double>(p.second) } );
-      cal_obj["DeviationPairs"] = pairs;
-    }//if( dev pairs )
+    if( cal_type == SpecUtils::EnergyCalType::InvalidEquationType )
+      return;
 
-    if( spec_file )
+    cal_obj["LowerEnergy"] = cal->lower_energy();
+    cal_obj["UpperEnergy"] = cal->upper_energy();
+
+    if( cal->type() != SpecUtils::EnergyCalType::LowerChannelEdge )
     {
-      spec_obj["InstrumentModel"] = spec_file->instrument_model();
-      spec_obj["SerialNumber"] = spec_file->instrument_id();
-      spec_obj["Manufacturer"] = spec_file->manufacturer();
-      spec_obj["InstrumentType"] = spec_file->instrument_type();
-      spec_obj["DetectorType"] = SpecUtils::detectorTypeToString( spec_file->detector_type() );
-      spec_obj["NumberRecordsInFile"] = static_cast<int>( spec_file->num_measurements() );
-      spec_obj["RemarksInFile"] = spec_file->remarks();
-      spec_obj["ParseWarningsForFile"] = spec_file->parse_warnings();
+      cal_obj["Coefficients"] = vector<double>{ begin(cal->coefficients()), end(cal->coefficients()) };
 
-      spec_obj["HasInstrumentRid"] = !!spec_file->detectors_analysis();
-      if( spec_file->detectors_analysis() )
-        spec_obj["InstrumentRidSummary"] = riidAnaSummary( spec_file );
-    }//if( spec_file )
-  }//add_hist_to_json(...)
+      if( !cal->deviation_pairs().empty() )
+      {
+        vector<vector<double>> pairs;
+        for( const pair<float,float> &p : cal->deviation_pairs() )
+          pairs.push_back( { static_cast<double>(p.first), static_cast<double>(p.second) } );
+        cal_obj["DeviationPairs"] = pairs;
+      }//if( dev pairs )
+    }//if( polynomial or FRF )
+  }//void add_energy_cal_json(...)
 
 
   void add_activity_fit_options_to_json( nlohmann::json &data,
@@ -1480,6 +1495,20 @@ void add_basic_src_details( const GammaInteractionCalc::SourceDetails &src,
     }//if( fit_results.spectrum )
 
 
+    const bool successful_ene_refit = (fit_results.original_energy_cal && fit_results.refit_energy_cal);
+    data["EnergyCalIsRefit"] = successful_ene_refit;
+
+    if( successful_ene_refit )
+    {
+      auto &ene_cal_json = data["EnergyCalRefitResult"];
+
+      if( fit_results.original_energy_cal )
+        add_energy_cal_json( ene_cal_json["Original"], fit_results.original_energy_cal );
+
+      if( fit_results.refit_energy_cal )
+        add_energy_cal_json( ene_cal_json["Refit"], fit_results.refit_energy_cal );
+    }//if( successful_ene_refit )
+
 
     auto add_peaks = []( nlohmann::json &json, deque<shared_ptr<const PeakDef>> peaks, const shared_ptr<const SpecUtils::Measurement> &spectrum ) {
 

src/BatchPeak.cpp

@@ -162,7 +162,7 @@ void propagate_energy_cal( const shared_ptr<const SpecUtils::EnergyCalibration>
   }//if( updated_cal )
 }//propagate_energy_cal(...)
 
-void fit_energy_cal_from_fit_peaks( shared_ptr<SpecUtils::Measurement> &raw, vector<PeakDef> peaks, const size_t num_coefs )
+void fit_energy_cal_from_fit_peaks( shared_ptr<SpecUtils::Measurement> &raw, vector<PeakDef> peaks )
 {
   if( !raw || raw->num_gamma_channels() < 16 )
     throw runtime_error( "update_gain_from_peak: invalid input spectrum" );
@@ -172,15 +172,39 @@ void fit_energy_cal_from_fit_peaks( shared_ptr<SpecUtils::Measurement> &raw, vec
   shared_ptr<const SpecUtils::EnergyCalibration> orig_cal = raw->energy_calibration();
   assert( orig_cal && orig_cal->valid() && (orig_cal->coefficients().size() > 1) );
 
+  const SpecUtils::EnergyCalType energy_cal_type = (orig_cal && orig_cal->valid()) 
+                                                    ? orig_cal->type()
+                                                    : SpecUtils::EnergyCalType::InvalidEquationType;
+  switch( energy_cal_type )
+  {
+    case SpecUtils::EnergyCalType::Polynomial:
+    case SpecUtils::EnergyCalType::UnspecifiedUsingDefaultPolynomial:
+    case SpecUtils::EnergyCalType::FullRangeFraction:
+      if( orig_cal->coefficients().size() < 2 )
+        throw std::logic_error( "Somehow the energy calibration has less than two coefficients." );
+      break;
+
+    case SpecUtils::EnergyCalType::LowerChannelEdge:
+    case SpecUtils::EnergyCalType::InvalidEquationType:
+      assert( 0 );
+      throw std::logic_error( "The original calibration must be either polynomial or full-range-fraction." );
+      break;
+  }//switch( exemplar_cal->type() )
+
+
   for( const PeakDef &peak : peaks )
   {
     if( !peak.hasSourceGammaAssigned() )
     {
       cerr << "Warning: peak at " << peak.mean() << " keV doesnt have a source assigned, so will"
       << " not be used for energy calibration" << endl;
-      continue;;
+      continue;
     }
 
+
+    if( !peak.useForEnergyCalibration() )  //This defaults to true, so if this is false, then user selected it
+      continue;
+
     EnergyCal::RecalPeakInfo recalpeak;
     recalpeak.peakMean = peak.mean();
     recalpeak.peakMeanUncert = peak.meanUncert();
@@ -190,22 +214,68 @@ void fit_energy_cal_from_fit_peaks( shared_ptr<SpecUtils::Measurement> &raw, vec
     peakinfos.push_back( recalpeak );
   }//for( const double peak_energy : peak_energies )
 
+  std::sort( begin(peakinfos), end(peakinfos), 
+    []( const EnergyCal::RecalPeakInfo &lhs, const EnergyCal::RecalPeakInfo &rhs ) -> bool {
+    return lhs.peakMean < rhs.peakMean;
+  } );
+
+  if( peakinfos.empty() )
+    throw runtime_error( "No peaks selected for use in energy calibration." );
+
   const size_t nchannels = raw->num_gamma_channels();
 
+  const vector<float> &orig_coefs = orig_cal->coefficients();
   const vector<pair<float,float>> &dev_pairs = orig_cal->deviation_pairs();
+  const size_t num_coefs = orig_coefs.size();
+  assert( num_coefs >= 2 );
 
   vector<float> fit_coefs_uncert;
-  vector<float> fit_coefs = orig_cal->coefficients();
+  vector<float> fit_coefs = orig_coefs;
+  vector<bool> fitfor( num_coefs, false );
+
+  // If we only have a couple peaks, then we cant fit for like 4 coefficients; we'll
+  //  just hardcode a rough heuristic for what coefficients to fit for, because I think
+  //  bothering the user with this level of detail is probably a bit too much.
+  //
+  //  But also note that we have already fit peaks, so we must be reasonably close
+  //  right now anyway, so we can be a bit more liberal in terms of fitting more
+  //  coefficients than a user might normally do during an interactive session.
+  //
+  //  We dont want to update both gain and offset from like the two Co60 peaks, so
+  //  we'll count the effective number of peaks, requiring them to be separated a bit.
+  //  We will require the separation to be max( 100, min(0.1*total-energy-range, 200)) keV,
+  //  with the 0.1, 100 and 200, all being entirely arbitrary, but hopefully reasonable.
+  size_t num_effective_peaks = 1;
+  const double energy_range = (orig_cal->upper_energy() - orig_cal->lower_energy());
+  const double sep_dist = std::max( 100.0, std::min(0.1*energy_range, 200.0) );
+  size_t last_peak = 0;
+  for( size_t peak_index = 1; peak_index < peakinfos.size(); ++peak_index )
+  {
+    const double delta_energy = peakinfos[peak_index].peakMean - peakinfos[last_peak].peakMean;
+    assert( delta_energy >= 0.0 );
+    if( delta_energy >= sep_dist )
+    {
+      num_effective_peaks += 1;
+      last_peak = peak_index;
+    }
+  }//for( loop over peaks to count how many are separated by at least `sep_dist` )
 
-  vector<bool> fitfor( orig_cal->coefficients().size(), false );
-  if( fitfor.size() < num_coefs )
-    fitfor.resize( num_coefs );
-  for( size_t i = 0; i < num_coefs; ++i )
-    fitfor[i] = true;
+  if( num_effective_peaks == 1 )
+  {
+    fitfor[1] = true; // Only fit gain
+  }else
+  {
+    for( size_t index = 0; (index < num_effective_peaks) && (index < fitfor.size()); ++index )
+      fitfor[index] = true;
+    // TODO: we could consider not fitting for offset if `peakinfos[0].peakMean` is less than 
+    //       ~200 keV or something, but for the moment we wont, because we know we are close
+    //       in energy calibration, and we know the peaks the user is interested in, so overfitting
+    //       isnt as large of a concern as not lining up the ROI edges as much
+  }
 
   shared_ptr<SpecUtils::EnergyCalibration> updated_cal = make_shared<SpecUtils::EnergyCalibration>();
 
-  switch( orig_cal->type() )
+  switch( energy_cal_type )
   {
     case SpecUtils::EnergyCalType::Polynomial:
     case SpecUtils::EnergyCalType::UnspecifiedUsingDefaultPolynomial:
@@ -226,6 +296,7 @@ void fit_energy_cal_from_fit_peaks( shared_ptr<SpecUtils::Measurement> &raw, vec
 
   raw->set_energy_calibration( updated_cal );
 
+#if( PERFORM_DEVELOPER_CHECKS )
   cout << "Updated energy calibration using ROIs from exemplar.\n\tCoefficients:\n";
   assert( fit_coefs.size() == orig_cal->coefficients().size() );
   for( size_t i = 0; i < fit_coefs.size(); ++i )
@@ -244,6 +315,7 @@ void fit_energy_cal_from_fit_peaks( shared_ptr<SpecUtils::Measurement> &raw, vec
   }
 
   cout << endl << endl;
+#endif
 }//void fit_energy_cal_from_fit_peaks(...)
 
 
@@ -1052,6 +1124,8 @@ BatchPeak::BatchPeakFitResult fit_peaks_in_file( const std::string &exemplar_fil
       candidate_peaks.push_back( peak );
     }//for( const auto &p : exemplar_peaks )
 
+    results.original_energy_cal = spec ? spec->energy_calibration() : nullptr;
+
     if( options.refit_energy_cal )
     {
       // We will refit the energy calibration - maybe a few times - to really hone in on things
@@ -1068,7 +1142,16 @@ BatchPeak::BatchPeakFitResult fit_peaks_in_file( const std::string &exemplar_fil
         vector<PeakDef> peaks = fitPeaksInRange( lower_energy, uppper_energy, ncausalitysigma,
                                                 stat_threshold, hypothesis_threshold,
                                                 energy_cal_peaks, spec, {}, isRefit );
-        fit_energy_cal_from_fit_peaks( spec, peaks, 4 );
+        try
+        {
+          fit_energy_cal_from_fit_peaks( spec, peaks );
+        }catch( std::exception &e )
+        {
+          const string msg = "Energy calibration not performed: " + string(e.what());
+          auto pos = std::find( begin(results.warnings), end(results.warnings), msg );
+          if( pos == end(results.warnings) )
+            results.warnings.push_back( msg );
+        }
       }//for( size_t i = 0; i < 1; ++i )
 
       // Propagate the updated energy cal to the result file
@@ -1087,6 +1170,8 @@ BatchPeak::BatchPeakFitResult fit_peaks_in_file( const std::string &exemplar_fil
       {
         results.warnings.push_back( "Failed to fit an appropriate energy calibration in '" + filename + "'." );
       }
+
+      results.refit_energy_cal = spec ? spec->energy_calibration() : nullptr;
     }//if( options.refit_energy_cal )
 
     vector<PeakDef> fit_peaks = fitPeaksInRange( lower_energy, uppper_energy, ncausalitysigma,

src/ExportSpecFile.cpp

@@ -2987,12 +2987,14 @@ std::shared_ptr<const SpecMeas> ExportSpecFileTool::generateFileToSave()
       samples.erase( sample );
 
     // Now lets map sample numbers to 1 through N
+    samples.clear();
     int new_sample_number = 0;
     vector<pair<int,int>> old_to_new_samplenum;
     for( const int old_sample_number : answer->sample_numbers() )
     {
       ++new_sample_number;
       old_to_new_samplenum.emplace_back( old_sample_number, new_sample_number );
+      samples.insert( new_sample_number );
     }//for( loop over sample numbers )
 
     answer->change_sample_numbers( old_to_new_samplenum );

target/electron/build_linux_app_from_docker.sh

@@ -21,7 +21,7 @@ trap 'handle_error $LINENO' ERR
 #  git clone --recursive [email protected]:sandialabs/InterSpec.git ./InterSpec_code
 #  mkdir build_electron
 #  mkdir build_working_dir
-#  docker run --rm -it -v `pwd`/InterSpec_code:/interspec -v `pwd`/build_electron/:/build_app -v `pwd`/build_working_dir:/build_working_dir quay.io/pypa/manylinux_2_28_x86_64:latest /interspec/target/electron/Docker/build_linux_app_from_docker.sh /interspec /build_app /build_working_dir
+#  docker run --rm -it -v `pwd`/InterSpec_code:/interspec -v `pwd`/build_electron/:/build_app -v `pwd`/build_working_dir:/build_working_dir quay.io/pypa/manylinux_2_28_x86_64:latest /interspec/target/electron/build_linux_app_from_docker.sh /interspec /build_app /build_working_dir
 #  #Then results will then be in build_working_dir/InterSpec-linux-x64
 
 if [ $# -ne 3 ]; then
@@ -44,7 +44,7 @@ export GIT_HASH=$(git -C ${InterSpecCodePath} rev-parse HEAD)
 echo "GIT_HASH = ${GIT_HASH}"
 
 echo "Will install npm and global packages"
-yum update
+yum update -y
 yum install -y npm zip
 # Use the 'n' package to install a fairly modern version of npm, by default we are on like version 6 or 8, which is too old to run some of the packages we need
 npm install -g n