Merge pull request #235 from aglowacki/master

Polar updates

Author: aglowacki

CMakeLists.txt

@@ -373,8 +373,6 @@ ELSE()
              )
 ENDIF()
 
-set(HDF5_LIB_LEN 0)
-list(LENGTH HDF5_LIBRARIES HDF5_LIB_LEN)
 
 #--------------- start pyxrf-maps lib -----------------
 IF (BUILD_WITH_PYBIND11)
@@ -385,24 +383,16 @@ IF (BUILD_WITH_PYBIND11)
     pybind11_add_module(pyxrfmaps src/pybindings/main.cpp)
     IF (BUILD_WITH_ZMQ)
       IF (WIN32)
-        IF(${HDF5_LIB_LEN} LESS 1)
-          target_link_libraries(pyxrfmaps PRIVATE libxrf_fit libxrf_io netCDF::netcdf libzmq-static ws2_32.lib rpcrt4.lib iphlpapi.lib)
-        ELSE()
           target_link_libraries(pyxrfmaps PRIVATE libxrf_fit libxrf_io netCDF::netcdf hdf5::hdf5-shared libzmq-static ws2_32.lib rpcrt4.lib iphlpapi.lib)
-        ENDIF()
       ELSEIF (UNIX)
-        IF(${HDF5_LIB_LEN} LESS 1)
-        target_link_libraries(pyxrfmaps PRIVATE libxrf_fit libxrf_io netCDF::netcdf libzmq-static)
-        ELSE()
-          target_link_libraries(pyxrfmaps PRIVATE libxrf_fit libxrf_io netCDF::netcdf hdf5::hdf5-shared libzmq-static)
-        ENDIF()
+          target_link_libraries(pyxrfmaps PRIVATE libxrf_fit libxrf_io netCDF::netcdf libzmq-static)
       ENDIF()
     ELSE()
-      IF(${HDF5_LIB_LEN} LESS 1)
-        target_link_libraries(pyxrfmaps PRIVATE libxrf_fit libxrf_io netCDF::netcdf)
-      ELSE()
-        target_link_libraries(pyxrfmaps PRIVATE libxrf_fit libxrf_io netCDF::netcdf hdf5::hdf5-shared)
-      ENDIF()  
+        IF (WIN32)
+          target_link_libraries(pyxrfmaps PRIVATE libxrf_fit libxrf_io netCDF::netcdf hdf5::hdf5-shared)
+      ELSEIF (UNIX)
+          target_link_libraries(pyxrfmaps PRIVATE libxrf_fit libxrf_io netCDF::netcdf libzmq-static)
+      ENDIF()
     ENDIF()
 ENDIF()
 
@@ -527,15 +517,13 @@ ELSEIF (UNIX)
   ENDIF()
 ENDIF()
 
-# Building on theta doesn't need hdf5 and it is empty os need a check 
-IF(${HDF5_LIB_LEN} LESS 1)
-  target_link_libraries (libxrf_fit PRIVATE NLopt::nlopt ${CMAKE_THREAD_LIBS_INIT} )
-  target_link_libraries(libxrf_io PRIVATE libxrf_fit netCDF::netcdf yaml-cpp::yaml-cpp JsonCpp::JsonCpp NLopt::nlopt ${CMAKE_THREAD_LIBS_INIT} )
-  target_link_libraries (xrf_maps PRIVATE libxrf_io libxrf_fit netCDF::netcdf yaml-cpp::yaml-cpp JsonCpp::JsonCpp NLopt::nlopt ${CMAKE_THREAD_LIBS_INIT} )
-ELSE()
-  target_link_libraries(libxrf_fit PRIVATE NLopt::nlopt ${CMAKE_THREAD_LIBS_INIT} )
+target_link_libraries(libxrf_fit PRIVATE NLopt::nlopt ${CMAKE_THREAD_LIBS_INIT} )
+IF (WIN32)
   target_link_libraries(libxrf_io PRIVATE libxrf_fit netCDF::netcdf hdf5::hdf5-shared yaml-cpp::yaml-cpp JsonCpp::JsonCpp NLopt::nlopt ${CMAKE_THREAD_LIBS_INIT} )
   target_link_libraries (xrf_maps PRIVATE libxrf_io libxrf_fit netCDF::netcdf hdf5::hdf5-shared yaml-cpp::yaml-cpp JsonCpp::JsonCpp NLopt::nlopt ${CMAKE_THREAD_LIBS_INIT} )
+ELSEIF (UNIX)
+  target_link_libraries(libxrf_io PRIVATE libxrf_fit netCDF::netcdf yaml-cpp::yaml-cpp JsonCpp::JsonCpp NLopt::nlopt ${CMAKE_THREAD_LIBS_INIT} )
+  target_link_libraries (xrf_maps PRIVATE libxrf_io libxrf_fit netCDF::netcdf yaml-cpp::yaml-cpp JsonCpp::JsonCpp NLopt::nlopt ${CMAKE_THREAD_LIBS_INIT} )
 ENDIF()
 
 IF (BUILD_WITH_QT)

src/core/defines.h

@@ -337,6 +337,8 @@ const std::string STR_LIMITED_LO = "limit_lo";
 const std::string STR_LIMITED_HI = "limit_hi";
 const std::string STR_FIT = "fit";
 
+const std::string STR_CUSTOM_EL_REGION = "CUSTOM_ELEMENT_REGION";
+
 // ROI
 const std::string STR_MAPS_ROIS = "MAPS_ROIS";
 const std::string STR_MAP_ROI_NAME = "Name";

src/data_struct/fit_element_map.cpp

@@ -81,9 +81,9 @@ Fit_Element_Map<T_real>::Fit_Element_Map(std::string name, Element_Info<T_real>*
     else
     {
         std::string str_shell = _full_name.substr(idx);
-        if (str_shell == "L")
+        if (str_shell == "M")
         {
-            _shell_type = Electron_Shell::L_SHELL;
+            _shell_type = Electron_Shell::M_SHELL;
         }
         else if (str_shell == "L")
         {
@@ -118,6 +118,33 @@ Fit_Element_Map<T_real>::Fit_Element_Map(std::string name, Element_Info<T_real>*
 
 //-----------------------------------------------------------------------------
 
+template<typename T_real>
+Fit_Element_Map<T_real>::Fit_Element_Map(std::string name, T_real center, T_real width_multi)
+{
+    _full_name = name;
+
+    _element_info = nullptr;
+
+    _pileup_element_info = nullptr;
+
+    _width_multi = width_multi;
+
+    _center = center;
+
+    _width = int( std::sqrt( std::pow(ENERGY_RES_OFFSET, 2) + std::pow( (_center * ENERGY_RES_SQRT), 2)  ) );
+    _width /= 2000.0;
+    _width *= _width_multi;
+
+    size_t num_ratios = 1;
+    _energy_ratios.push_back(Element_Energy_Ratio<T_real>(_center, 1.0, 0.0, Element_Param_Type::Ka1_Line));
+    for(size_t i=0; i<num_ratios; i++)
+    {
+        _energy_ratio_custom_multipliers.push_back(1.0);
+    }
+}
+
+//-----------------------------------------------------------------------------
+
 template<typename T_real>
 Fit_Element_Map<T_real>::~Fit_Element_Map()
 {
@@ -135,6 +162,10 @@ void Fit_Element_Map<T_real>::init_energy_ratio_for_detector_element(const Eleme
 
     if(_element_info == nullptr)
     {
+        _width = int( std::sqrt( std::pow(ENERGY_RES_OFFSET, 2) + std::pow( (_center * ENERGY_RES_SQRT), 2)  ) );
+        _width /= 2000.0;
+        _width *= _width_multi;
+
         //Don't log because we can have non elements such as Compton that doesn't have element properties
         //logE<<"Element info was not properly loaded. Variable is null! Can't initialize energy ratios!\n";
         return;

src/data_struct/fit_element_map.h

@@ -159,6 +159,8 @@ class DLL_EXPORT Fit_Element_Map
 public:
     Fit_Element_Map(std::string name, Element_Info<T_real>* element_info);
 
+    Fit_Element_Map(std::string name, T_real center, T_real width_multi);
+
     ~Fit_Element_Map();
 
     const T_real& center() const { return _center; }
@@ -194,6 +196,8 @@ class DLL_EXPORT Fit_Element_Map
     std::unordered_map<std::string, T_real> generate_roi_centers_per_shell();
 
 	bool check_binding_energy(T_real incident_energy, size_t energy_ratio_idx) const;
+
+    bool has_element_info() { return !(_element_info==nullptr);}
 protected:
 
     void generate_energy_ratio(T_real energy, T_real ratio, Element_Param_Type et, const Element_Info<T_real> * const detector_element);

src/io/file/aps/aps_fit_params_import.h

@@ -468,6 +468,31 @@ DLL_EXPORT bool load_parameters_override(std::string path, Params_Override<T_rea
                             fit_map->set_width_multi(factor);
                         }
                     }
+                    else if (tag == STR_CUSTOM_EL_REGION)
+                    {
+                        std::string element_symb;
+                        std::string energy_value;
+                        std::string width_value;
+
+                        std::getline(strstream, element_symb, ',');
+                        element_symb.erase(std::remove_if(element_symb.begin(), element_symb.end(), ::isspace), element_symb.end());
+
+                        T_real center = 1.0;
+                        std::getline(strstream, energy_value, ',');
+                        center = parse_input_real<T_real>(energy_value);
+    
+                        T_real width_multi = 1.0;
+                        std::getline(strstream, width_value, ',');
+                        width_multi = parse_input_real<T_real>(width_value);
+
+                        Fit_Element_Map<T_real>* fit_map;
+                        if (params_override->elements_to_fit.count(element_symb) < 1)
+                        {
+                            
+                            fit_map = new Fit_Element_Map<T_real>(element_symb, center, width_multi);
+                            params_override->elements_to_fit[element_symb] = fit_map;
+                        }
+                    }
                     else if (tag == "BRANCHING_FAMILY_ADJUSTMENT_L" || tag == "BRANCHING_RATIO_ADJUSTMENT_L" || tag == "BRANCHING_RATIO_ADJUSTMENT_K" || tag == "BRANCHING_RATIO_ADJUSTMENT_M")
                     {
                         unsigned int cnt = 0;
@@ -954,12 +979,24 @@ DLL_EXPORT bool save_parameters_override(std::string path, Params_Override<T_rea
         for (const auto& itr : params_override->elements_to_fit)
         {
             // if not pileup
-            if (itr.second->pileup_element() == nullptr && (itr.first != "COHERENT_SCT_AMPLITUDE" || itr.first != "COMPTON_AMPLITUDE"))
+            if (itr.second->pileup_element() == nullptr && itr.second->has_element_info())
             {
                 out_stream << itr.first << " , ";
             }
         }
         out_stream << "\n";
+        for (const auto& itr : params_override->elements_to_fit)
+        {
+            // if not pileup
+            if (itr.second->pileup_element() == nullptr && false == itr.second->has_element_info())
+            {
+                if(itr.first == STR_COHERENT_SCT_AMPLITUDE || itr.first == STR_COMPTON_AMPLITUDE)
+                {
+                    continue;
+                }
+                out_stream << STR_CUSTOM_EL_REGION<<": "<< itr.first << "," << itr.second->center() << ","  << itr.second->width_multi() << "\n";
+            }
+        }
         out_stream << "   list the element combinations you want to fit for pileup, e.g., Si_Si, Si_Si_Si, Si_Cl, etc\n";
         out_stream << "ELEMENTS_WITH_PILEUP: ";
         for (const auto& itr : params_override->elements_to_fit)

src/io/file/hdf5_io.h

@@ -1658,6 +1658,13 @@ class DLL_EXPORT HDF5_IO
                     ocr_array.setOnes();
                 }
             }
+            else
+            {
+                elt_array.setOnes();
+                ert_array.setOnes();
+                icr_array.setOnes();
+                ocr_array.setOnes();
+            }
             err = _read_h5d<T_real2>(energy_dset_id, H5S_ALL, energy_space_id, H5P_DEFAULT, energy_array.data());
             if(err < 0)
             {
@@ -1694,6 +1701,9 @@ class DLL_EXPORT HDF5_IO
             int left_pol_idx = 0;
             int right_pol_idx = 0;
 
+            int scaler_idx = 0;
+            int scaler_cntr = 0;
+            int scaler_inc = params_override->polarity_pattern.size();
             for(size_t i = 0; i < dims3[0]; i++)
             {
                 offset3[0] = i;
@@ -1774,11 +1784,21 @@ class DLL_EXPORT HDF5_IO
                 (*spec_vol)[polarity][idx].input_counts(icr_array[i]);
                 (*spec_vol)[polarity][idx].output_counts(ocr_array[i]);
 
-                energy_map.values(polarity, idx) = energy_array[i];
-                i0_map.values(polarity, idx) = i0_array[i];
+                
+                energy_map.values(0, idx) = energy_array[scaler_idx];
+                energy_map.values(1, idx) = energy_array[scaler_idx];
+                i0_map.values(0, idx) = i0_array[scaler_idx];
+                i0_map.values(1, idx) = i0_array[scaler_idx];
 
-                dtf_map.values(polarity, idx) = dtf_array[i];
-                dtp_map.values(polarity, idx) = dtp_array[i];
+                scaler_cntr ++;
+                if(scaler_cntr > scaler_inc)
+                {
+                    scaler_cntr = 0;
+                    scaler_idx ++;
+                }   
+                
+                dtf_map.values(0, idx) = dtf_array[i];
+                dtp_map.values(0, idx) = dtp_array[i];
             }
 
             scan_info.scaler_maps.push_back(dtf_map);

src/io/file/hl_file_io.h

@@ -688,22 +688,24 @@ DLL_EXPORT bool load_and_integrate_spectra_volume(std::string dataset_directory,
                         size_t spec_size = io::file::NetCDF_IO<T_real>::inst()->load_scalers_line(full_filename, "tetra2_", i, mda_io.get_scan_info(), params_override);
                     }
                 }
-                T_real val = mda_io.get_scan_info()->scaler_avg_value(STR_US_IC);
-                if(val > (T_real)0.0)
+                if(params_override != nullptr)
                 {
-                    params_override->US_IC = val;
-                }
-                val = mda_io.get_scan_info()->scaler_avg_value(STR_DS_IC);
-                if(val > (T_real)0.0)
-                {
-                    params_override->DS_IC = val;
-                }
-                val = mda_io.get_scan_info()->scaler_avg_value(STR_US_FM);
-                if(val > (T_real)0.0)
-                {
-                    params_override->US_FM = val;
+                    T_real val = mda_io.get_scan_info()->scaler_avg_value(STR_US_IC);
+                    if(val > (T_real)0.0)
+                    {
+                        params_override->US_IC = val;
+                    }
+                    val = mda_io.get_scan_info()->scaler_avg_value(STR_DS_IC);
+                    if(val > (T_real)0.0)
+                    {
+                        params_override->DS_IC = val;
+                    }
+                    val = mda_io.get_scan_info()->scaler_avg_value(STR_US_FM);
+                    if(val > (T_real)0.0)
+                    {
+                        params_override->US_FM = val;
+                    }
                 }
-                
             }
             if (hasNetcdf)
             {
@@ -1232,13 +1234,14 @@ DLL_EXPORT bool load_spectra_volume(std::string dataset_directory,
     }
 
     // try to load spectra from mda file
-    if (false == mda_io.load_spectra_volume(dataset_directory + "mda" + DIR_END_CHAR + dataset_file, detector_num, spectra_volume, (hasNetcdf || hasBnpNetcdf || hasHdf || hasXspress), hasNetcdf))
+    auto mda_ret_val = mda_io.load_spectra_volume(dataset_directory + "mda" + DIR_END_CHAR + dataset_file, detector_num, spectra_volume, (hasNetcdf || hasBnpNetcdf || hasHdf || hasXspress), hasNetcdf);
+    if (Load_Status::Failed == mda_ret_val)
     {
         scan_type = STR_SCAN_TYPE_2D_MAP;
         logE << "Load spectra " << dataset_directory + "mda" + DIR_END_CHAR + dataset_file << "\n";
         return false;
     }
-    else
+    else if (Load_Status::Half_need_spectra == mda_ret_val)
     {
         // tetramm is scalers only
         if(hasTetraMM)