Skip to content

Commit

Permalink
Fix simple detection limit bugs.
Browse files Browse the repository at this point in the history
The right-clicked energy wasnt always being selected in Simple MDA tool, due to using the wrong/incorrect index to select gamma.
Also, fixed a few bugs that could cause crash if no result was found.
  • Loading branch information
wcjohns committed Jul 26, 2024
1 parent a070967 commit 6bd0c8a
Show file tree
Hide file tree
Showing 4 changed files with 45 additions and 30 deletions.
1 change: 1 addition & 0 deletions InterSpec_resources/app_text/DetectionLimitSimple.xml
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,7 @@
<message id="dls-qr-tool-state-txt">Current state of the Simple MDA tool.</message>

<message id="dls-det-act-with-range">Detected: {1}, range [{2}, {3}], @{4} CL</message>
<message id="dls-det-no-result">Sorry, no result</message>
<message id="dls-det-counts-with-range">Excess counts: {1}, range [{2}, {3}], @{4} CL</message>
<message id="dls-det-act-less-zero">Activity &lt; 0 {1} (observed significantly fewer counts than expected).</message>
<message id="dls-det-counts-less-zero">Excess counts &lt; 0 {2} (observed significantly fewer counts than expected).</message>
Expand Down
4 changes: 2 additions & 2 deletions src/DetectionLimitCalc.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -753,8 +753,8 @@ pair<size_t,size_t> round_roi_to_channels( shared_ptr<const SpecUtils::Measureme
if( !cal || !cal->valid() )
throw runtime_error( "mda_counts_calc: invalid energy calibration" );

const float peak_region_lower_ch = cal->channel_for_energy( roi_lower_energy );
const float peak_region_upper_ch = cal->channel_for_energy( roi_upper_energy );
const double peak_region_lower_ch = std::max(0.0, cal->channel_for_energy( roi_lower_energy ) );
const double peak_region_upper_ch = std::max(0.0, cal->channel_for_energy( roi_upper_energy ) );

//if( (peak_region_lower_ch - num_lower_side_channels) < 0.0 )
// throw runtime_error( "mda_counts_calc: lower energy goes off spectrum" );
Expand Down
65 changes: 39 additions & 26 deletions src/DetectionLimitSimple.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -824,6 +824,7 @@ void DetectionLimitSimple::setNuclide( const SandiaDecay::Nuclide *nuc,
m_currentEnergy = energy;

m_nucEnterController->setNuclideText( nuc ? nuc->symbol : string() );
assert( m_currentNuclide == nuc );
m_currentNuclide = nuc;

if( (age >= 0.0) && !m_nucEnterController->nuclideAgeStr().empty() )
Expand Down Expand Up @@ -998,6 +999,8 @@ void DetectionLimitSimple::handleNuclideChanged()
? det->peakResolutionSigma( static_cast<float>(energyToSelect) )
: -1.0;

const float fwhm = std::max( 0.1f, PeakSearchGuiUtils::estimate_FWHM_of_foreground(energyToSelect) );

size_t transition_index = 0;
const SandiaDecay::Transition *transition = nullptr;
PeakDef::SourceGammaType sourceGammaType = PeakDef::SourceGammaType::NormalGamma;
Expand All @@ -1015,36 +1018,38 @@ void DetectionLimitSimple::handleNuclideChanged()
}


double minDE = DBL_MAX;
double min_scale_delta_e = DBL_MAX;
int currentIndex = -1;

for( size_t i = 0; i < photons.size(); ++i )
{
double energy = photons[i].energy;
const double intensity = photons[i].numPerSecond / dummy_activity;

const double dE = fabs( energyToSelect - energy );
energy = floor(10000.0*energy + 0.5)/10000.0;

if( dE < minDE )
{
minDE = dE;
currentIndex = static_cast<int>( i );
}

char text[128] = { '\0' };
if( i < xrays.size() )
{
snprintf( text, sizeof(text), "%.4f keV xray I=%.1e", energy, intensity );
}else
{
snprintf( text, sizeof(text), "%.4f keV I=%.1e", energy, intensity );
}//if( i < xrays.size() )

if( intensity > std::numeric_limits<double>::epsilon() )
{
const double delta_e = fabs( energyToSelect - energy );
const double scale_delta_e = (0.1*fwhm + delta_e) / intensity;

energy = floor(10000.0*energy + 0.5)/10000.0;

char text[128] = { '\0' };
if( i < xrays.size() )
{
snprintf( text, sizeof(text), "%.4f keV xray I=%.1e", energy, intensity );
}else
{
snprintf( text, sizeof(text), "%.4f keV I=%.1e", energy, intensity );
}//if( i < xrays.size() )

m_photoPeakEnergiesAndBr.push_back( make_pair(energy, intensity) );
m_photoPeakEnergy->addItem( text );

if( scale_delta_e < min_scale_delta_e )
{
min_scale_delta_e = scale_delta_e;
currentIndex = static_cast<int>( m_photoPeakEnergiesAndBr.size() - 1 );
}
}//if( intensity > 0.0 )
}//for each( const double energy, energies )

Expand All @@ -1068,7 +1073,7 @@ void DetectionLimitSimple::handleGammaChanged()
const int gamma_index = m_photoPeakEnergy->currentIndex();
assert( gamma_index < static_cast<int>(m_photoPeakEnergiesAndBr.size()) );

if( (gamma_index >= 0) && (gamma_index < static_cast<int>(m_photoPeakEnergiesAndBr.size())) )
if( (gamma_index >= 0) && (gamma_index < static_cast<int>(m_photoPeakEnergiesAndBr.size())) )
{
m_currentEnergy = m_photoPeakEnergiesAndBr[gamma_index].first;
}else
Expand Down Expand Up @@ -1326,14 +1331,19 @@ void DetectionLimitSimple::updateSpectrumDecorationsAndResultText()
const bool use_curie = use_curie_units();
const DetectorPeakResponse::EffGeometryType det_geom = drf ? drf->geometryType() : DetectorPeakResponse::EffGeometryType::FarField;

assert( (result->input.num_lower_side_channels != 0)
assert( !result
|| (result->input.num_lower_side_channels != 0)
|| (result->input.num_upper_side_channels != 0)
|| (result->input.num_lower_side_channels == result->input.num_upper_side_channels) );

const bool assertedIsBackground = ((result->input.num_lower_side_channels == 0)
const bool assertedIsBackground = (result
&& (result->input.num_lower_side_channels == 0)
&& (result->input.num_upper_side_channels == 0));

if( result->source_counts > result->decision_threshold )
if( !result )
{
result_txt = WString::tr("dls-det-no-result");
}else if( result->source_counts > result->decision_threshold )
{
assert( !assertedIsBackground );

Expand Down Expand Up @@ -1396,7 +1406,10 @@ void DetectionLimitSimple::updateSpectrumDecorationsAndResultText()

WString full_result_txt;

if( assertedIsBackground )
if( !result )
{
full_result_txt = result_txt;
}else if( assertedIsBackground )
{
full_result_txt = WString( "{1}" );
}else
Expand All @@ -1405,14 +1418,14 @@ void DetectionLimitSimple::updateSpectrumDecorationsAndResultText()
full_result_txt.arg(result_txt);
}//if( assertedIsBackground ) / else

if( gammas_per_bq > 0.0 )
if( result && (gammas_per_bq > 0.0) )
{
const double detection_act = result->detection_limit / gammas_per_bq;
const string act = PhysicalUnits::printToBestActivityUnits( detection_act, 2, use_curie )
+ DetectorPeakResponse::det_eff_geom_type_postfix( det_geom );

full_result_txt.arg( WString::tr("dls-min-detectable-act").arg(act) );
}else
}else if( result )
{
const string counts = SpecUtils::printCompact(result->detection_limit, 4);
full_result_txt.arg( WString::tr("dls-min-detectable-counts").arg( counts ) );
Expand Down
5 changes: 3 additions & 2 deletions src/DetectionLimitTool.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -1488,7 +1488,8 @@ void DetectionLimitTool::update_spectrum_for_currie_result( D3SpectrumDisplayDiv

// If number of of side channels is zero, then we expect the following
// variables to all be zero
assert( (result->input.num_lower_side_channels != 0)
assert( !result
|| (result->input.num_lower_side_channels != 0)
|| ((result->first_lower_continuum_channel == 0)
&& (result->last_lower_continuum_channel == 0)
&& (result->lower_continuum_counts_sum == 0)
Expand All @@ -1497,7 +1498,7 @@ void DetectionLimitTool::update_spectrum_for_currie_result( D3SpectrumDisplayDiv
&& (result->upper_continuum_counts_sum == 0))
);

const bool assertedNoSignal = (result->input.num_lower_side_channels == 0);
const bool assertedNoSignal = (input.num_lower_side_channels == 0);

// We will prefer to set the highlight regions from the results, but if we dont have results
// we'll do it from the input.
Expand Down

0 comments on commit 6bd0c8a

Please sign in to comment.