diff --git a/InterSpec/DetectorPeakResponse.h b/InterSpec/DetectorPeakResponse.h index 07cb624f..25c326b5 100644 --- a/InterSpec/DetectorPeakResponse.h +++ b/InterSpec/DetectorPeakResponse.h @@ -324,7 +324,7 @@ class DetectorPeakResponse const float upperEnergy, const EffGeometryType geometry_type ); - /** Makes a callable funtion from the provided mathematical formula + /** Makes a callable function from the provided mathematical formula @param fcn Mathematical formula to evaluate. Example: "exp( 1.2 + 3.2*ln(x) + -2.1*ln(x)^2 )" @param isMeV Whether units are in MeV, or keV @@ -417,6 +417,16 @@ class DetectorPeakResponse std::vector &credits, std::vector> &drfs ); + /** Parses detector efficiency function originating from GammaQuant. + These detector efficiency functions are defined in a column format. + A major limitation is that no CSV cell may have a new-line in it. + + Throws std::exception on format error, or now det eff functions found + */ + static void parseGammaQuantRelEffDrfCsv( std::istream &input, + std::vector> &drfs, + std::vector &warnings ); + /** Creates a DetectorPeakResponse from "App URL" data. Takes in just the "query" portion of the URL (i.e., the data after the '?' character), that diff --git a/InterSpec/SpecMeasManager.h b/InterSpec/SpecMeasManager.h index 96b343e6..6a8f78b6 100644 --- a/InterSpec/SpecMeasManager.h +++ b/InterSpec/SpecMeasManager.h @@ -196,7 +196,16 @@ class SpecMeasManager : public Wt::WObject bool handleMultipleDrfCsv( std::istream &input, const std::string &displayName, const std::string &fileLocation ); - + + /** Reads in GammaQuant CSV of detector efficiencies. + + The first cell must be "Detector ID", then rows through "Coefficient h". + Each column is a different detector, with the first column being the labels ("Detector ID", "Calibration Geometry", + "Comments", ..., "Coefficient h"). + */ + bool handleGammaQuantDrfCsv( std::istream &input, + const std::string &displayName, + const std::string &fileLocation ); /** Reads a CALp file from input stream and then either applies the CALp to current spectra, or prompts the user how to apply it. Function may return asynchronously to the CALp being applied, as the application may prompt user for options/input. diff --git a/src/DetectorPeakResponse.cpp b/src/DetectorPeakResponse.cpp index 4811daab..f646d601 100644 --- a/src/DetectorPeakResponse.cpp +++ b/src/DetectorPeakResponse.cpp @@ -1272,6 +1272,293 @@ void DetectorPeakResponse::parseMultipleRelEffDrfCsv( std::istream &input, }//parseMultipleRelEffDrfCsv(...) + +void DetectorPeakResponse::parseGammaQuantRelEffDrfCsv( std::istream &input, + std::vector> &drfs, + std::vector &warnings ) +{ + drfs.clear(); + warnings.clear(); + + const std::streampos start_pos = input.tellg(); + + try + { + string line; + + // Allow for some empty lines at the top + while( SpecUtils::safe_get_line( input, line, 128*1024 ) ) + { + SpecUtils::trim( line ); + if( line.empty() || (line[0] == '#') ) + continue; + break; + } + + if( line.empty() ) + throw runtime_error( "No content" ); + + size_t row_num = 0; + vector detector_id; + split_escaped_csv( detector_id, line ); + const size_t num_col = detector_id.size(); + if( (num_col < 2) || !SpecUtils::icontains(detector_id[0], "Detector ID") ) + throw runtime_error( "First column must be row headers, and start with \"Detector ID\"" ); + + // Define a utility lambda to help get each line + auto get_row = [num_col, &line, &row_num, &input, &warnings]( const bool is_empty, const string &expected_title ) + -> vector { + row_num += 1; + if( !SpecUtils::safe_get_line( input, line, 128*1024 ) ) + throw runtime_error( "Not enough rows - only " + std::to_string(row_num) + " rows." ); + + // TODO: check if last cell has newline in it, and if so, get the next line, and append to current. + + SpecUtils::trim( line ); + vector cols; + split_escaped_csv( cols, line ); + + if( !is_empty && (cols.size() != num_col) ) + throw runtime_error( "Inconsistent number of columns on line " + std::to_string(row_num) + + " - expected " + std::to_string(num_col) + " and received " + + std::to_string(cols.size()) ); + + SpecUtils::ireplace_all(cols[0], "\n", " "); + SpecUtils::ireplace_all(cols[0], "\r", " "); + SpecUtils::ireplace_all(cols[0], "\t", " "); + while( SpecUtils::icontains(cols[0], " ") ) + SpecUtils::ireplace_all(cols[0], " ", " "); + + if( !SpecUtils::icontains(cols[0], expected_title) ) + warnings.push_back( "Warning: row " + std::to_string(row_num) + + " doesnt have expected header '" + expected_title + "'" ); + + return cols; + };//get_row lambda + + + const vector calibration_geometry = get_row(false, "Calibration Geometry"); + const vector comments = get_row(false, "Comments"); + const vector is_far_field = get_row(false, "Far-Field Point Source Eff Cal?"); + get_row(true, ""); + const vector det_radius = get_row(false, "Detector Radius - a (cm)"); + const vector distance = get_row(false, "Source to Detector Face - d (cm)"); + get_row(false, "Source Radius - r (cm)"); //We dont use this info, yet + const vector detector_setbacks = get_row(false, "Detector Setback (cm)"); //We dont use this info, yet + get_row(false, "Detector Area - (cm2)"); + get_row(false, "Source to Crystal Face - d (cm)"); + get_row(false, "Relative Efficiency (%)"); + get_row(true, ""); + get_row(true, "EFFICIENCY CALIBRATION AND CURVE FIT INFORMATION"); + const vector lower_energy = get_row(false, "Lower Energy Cutoff (keV)"); + const vector upper_energy = get_row(false, "Upper Energy Cutoff (keV)"); + get_row(true, ""); + get_row(true, "y = exp (a + b(lnx) + c(lnx)2 + d(lnx)3 + e(lnx)4 + f(lnx)5 + g(lnx)6 + h(lnx)7)"); + const vector keV_or_MeV = get_row(false, "keV or MeV"); + const vector coef_a = get_row(false, "Coefficient a"); + const vector coef_b = get_row(false, "Coefficient b"); + const vector coef_c = get_row(false, "Coefficient c"); + const vector coef_d = get_row(false, "Coefficient d"); + const vector coef_e = get_row(false, "Coefficient e"); + const vector coef_f = get_row(false, "Coefficient f"); + vector coef_g, coef_h; + try + { + coef_g = get_row(false, "Coefficient g"); + coef_h = get_row(false, "Coefficient h"); + }catch( std::exception & ) + { + //G and H arent used that much anyway + } + + for( size_t col = 1; col < num_col; ++col ) + { + const string &name = detector_id[col]; + if( name.empty() ) + continue; + + const string &comment = comments[col]; + + const string &far_field = is_far_field[col]; + if( far_field.empty() ) + continue; + + const bool is_fixed_geom = SpecUtils::iequals_ascii( far_field, "No" ); + const bool is_far_field_val = SpecUtils::iequals_ascii( far_field, "Yes" ); + + if( !is_fixed_geom && !is_far_field_val ) + { + warnings.push_back( "Far-Field row of column " + std::to_string(col) + + " should be either Yes or No, was '" + far_field + "'" ); + continue; + } + + double det_rad = 0, src_to_det_face = 0, det_setback = 0, low_energy = 0, up_energy = 0; + + const string &det_rad_str = det_radius[col]; + if( !SpecUtils::parse_double(det_rad_str.c_str(), det_rad_str.size(), det_rad) ) + { + warnings.push_back( "Failed to parse detector radius, '" + det_rad_str + + "', to a number, skipping detector '" + name + "'" ); + continue; + } + det_rad *= PhysicalUnits::cm; + if( det_rad < 0.0 ) + { + warnings.push_back( "Detector radius negative ('" + det_rad_str + + "') - skipping detector '" + name + "'" ); + continue; + } + + const string &dist = distance[col]; + if( !SpecUtils::parse_double(dist.c_str(), dist.size(), src_to_det_face) ) + { + warnings.push_back( "Failed to parse distance to detector, '" + dist + + "', to a number, skipping detector '" + name + "'" ); + continue; + } + src_to_det_face *= PhysicalUnits::cm; + + if( src_to_det_face < 0.0 ) + { + warnings.push_back( "Src to face distance, '" + dist + + "', is negative - skipping detector '" + name + "'" ); + continue; + } + + const string &setback = detector_setbacks[col]; + if( !SpecUtils::parse_double(setback.c_str(), setback.size(), det_setback) ) + { + warnings.push_back( "Failed to parse detector setback, '" + setback + + "', to a number, skipping detector '" + name + "'" ); + continue; + } + + + const string &low_energy_str = lower_energy[col]; + if( !SpecUtils::parse_double(low_energy_str.c_str(), low_energy_str.size(), low_energy) ) + { + low_energy = 0; + warnings.push_back( "Failed to parse lower energy, '" + low_energy_str + + "', to a number, for detector '" + name + "'" ); + //continue; + } + + const string &up_energy_str = upper_energy[col]; + if( !SpecUtils::parse_double(up_energy_str.c_str(), up_energy_str.size(), up_energy) ) + { + up_energy = 0; + low_energy = 0; + warnings.push_back( "Failed to parse upper energy, '" + up_energy_str + + "', to a number, for detector '" + name + "'" ); + //continue; + } + + const string &keV_or_MeV_str = keV_or_MeV[col]; + const bool is_keV = SpecUtils::iequals_ascii(keV_or_MeV_str, "keV"); + const bool is_MeV = SpecUtils::iequals_ascii(keV_or_MeV_str, "MeV"); + if( !is_keV && !is_MeV ) + { + warnings.push_back( "The keV/MeV field for detector '" + name + "', is invalid, '" + + keV_or_MeV_str+ "' (must be either keV or MeV_ - skipping detector" ); + continue; + } + const double energy_units = is_keV ? PhysicalUnits::keV : PhysicalUnits::MeV; + + vector coef_strs; + coef_strs.push_back( coef_a[col] ); + coef_strs.push_back( coef_b[col] ); + coef_strs.push_back( coef_c[col] ); + coef_strs.push_back( coef_d[col] ); + coef_strs.push_back( coef_e[col] ); + coef_strs.push_back( coef_f[col] ); + if( (coef_g.size() > col) && (!coef_g[col].empty()) ) + { + coef_strs.push_back( coef_g[col] ); + if( (coef_h.size() > col) && (!coef_h[col].empty()) ) + coef_strs.push_back( coef_h[col] ); + }//if( we have coef g or h ) + + bool success_parse_all_coef = true; + size_t last_non_zero_coef = 0; + vector coef_values; + for( size_t coef_index = 0; coef_index < coef_strs.size(); ++coef_index ) + { + const string &str = coef_strs[coef_index]; + float coef; + if( !SpecUtils::parse_float(str.c_str(), str.size(), coef) ) + { + success_parse_all_coef = false; + break; + } + + coef_values.push_back( coef ); + + if( coef != 0.0f ) + last_non_zero_coef = coef_index; + }//for( loop over coefficients ) + + if( !success_parse_all_coef ) + { + warnings.push_back( "Failed to parse all eqn coefficients for detector '" + name + "' - skipping detector" ); + continue; + } + + if( last_non_zero_coef < 1 ) + { + warnings.push_back( "No non-zero eqn coefficients for detector '" + name + "' - skipping detector" ); + continue; + } + + coef_values.resize( last_non_zero_coef + 1 ); + + EffGeometryType geometry_type = EffGeometryType::FarField; + if( is_fixed_geom ) + { + if( SpecUtils::icontains(name, "Bq-cm2") ) + { + geometry_type = EffGeometryType::FixedGeomActPerCm2; + }else if( SpecUtils::icontains(name, "OnContact") || SpecUtils::icontains(comment, "OnContact") ) + { + geometry_type = EffGeometryType::FixedGeomTotalAct; + }else + { + warnings.push_back( "Couldnt determine fixed-geometry type for detector '" + name + "' - skipping detector" ); + continue; + } + }//if( is_fixed_geom ) + + try + { + auto det = make_shared(); + det->fromExpOfLogPowerSeriesAbsEff( coef_values, {}, src_to_det_face, 2*det_rad, + energy_units, low_energy, up_energy, geometry_type ); + if( !det->isValid() ) + throw runtime_error( "Not valid" ); + + det->setName( name ); + det->setDescription( comment ); + + drfs.push_back( det ); + }catch( std::exception &e ) + { + warnings.push_back( "Detector '" + name + "' was invalid ('" + string(e.what()) + + "') - skipping detector" ); + } + }//for( size_t col = 1; col < num_col; ++col ) + + if( drfs.empty() ) + throw runtime_error( "No valid detector efficiencies found." ); + }catch( std::exception &e ) + { + cerr << "Failed to parse GammaQuant Det Effs: " << e.what() << endl; + input.seekg( start_pos, ios::beg ); + input.clear( ios::failbit ); + throw; + } +}//void DetectorPeakResponse::parseGammaQuantRelEffDrfCsv(...) + + std::shared_ptr DetectorPeakResponse::parseFromAppUrl( const std::string &url_query ) { auto drf = make_shared(); diff --git a/src/SpecMeasManager.cpp b/src/SpecMeasManager.cpp index 3f4571e3..e1d82099 100644 --- a/src/SpecMeasManager.cpp +++ b/src/SpecMeasManager.cpp @@ -2326,6 +2326,13 @@ bool SpecMeasManager::handleNonSpectrumFile( const std::string &displayName, return true; } + // Check if this is TSV/CSV file containing multiple DRFs + if( header_contains( "Detector ID" ) + && handleGammaQuantDrfCsv(infile, displayName, fileLocation) ) + { + delete dialog; + return true; + } // Check if this is a PeakEasy CALp file if( currdata @@ -2510,6 +2517,34 @@ bool SpecMeasManager::handleMultipleDrfCsv( std::istream &input, }//bool handleMultipleDrfCsv( std::istream &input, SimpleDialog *dialog ) +bool SpecMeasManager::handleGammaQuantDrfCsv( std::istream &input, + const std::string &displayName, + const std::string &fileLocation ) +{ + vector warnings; + vector> drfs; + + try + { + DetectorPeakResponse::parseGammaQuantRelEffDrfCsv( input, drfs, warnings ); + }catch( std::exception &e ) + { + return false; + } + + if( drfs.empty() ) + return false; + + std::function saveDrfFile; + WString dialogmsg = WString::tr( (drfs.size()==1) ? "smm-file-is-drf" : "smm-file-is-multi-drf"); + + string creditsHtml; + DrfSelect::createChooseDrfDialog( drfs, dialogmsg, creditsHtml, saveDrfFile ); + + return true; +}//handleGammaQuantDrfCsv(...) + + bool SpecMeasManager::handleCALpFile( std::istream &infile, SimpleDialog *dialog, bool autoApply ) { WGridLayout *stretcher = nullptr;