From f21d7f914b05def2f3372750e86c8cc923563a43 Mon Sep 17 00:00:00 2001
From: cclaessens <claessens.christine@hotmail.de>
Date: Fri, 12 Jun 2020 16:45:01 +0200
Subject: [PATCH 01/49] can use helium instead of krypton as second gas

---
 mermithid/misc/ComplexLineShapeUtilities.py   |  2 ++
 .../processors/misc/KrComplexLineShape.py     | 28 +++++++++----------
 2 files changed, 16 insertions(+), 14 deletions(-)

diff --git a/mermithid/misc/ComplexLineShapeUtilities.py b/mermithid/misc/ComplexLineShapeUtilities.py
index 074013d0..9e7d48f0 100644
--- a/mermithid/misc/ComplexLineShapeUtilities.py
+++ b/mermithid/misc/ComplexLineShapeUtilities.py
@@ -58,6 +58,8 @@ def aseev_func_tail(energy_loss_array, gas_type):
         A2, omeg2, eps2 = 0.195, 14.13, 10.60
     elif gas_type=="Kr":
         A2, omeg2, eps2 = 0.4019, 22.31, 16.725
+    elif gas_type=="He":
+        A2, omeg2, eps2 = 0.1187, 33.40, 10.43
     return A2*omeg2**2./(omeg2**2.+4*(energy_loss_array-eps2)**2.)
 
 #convert oscillator strength into energy loss spectrum
diff --git a/mermithid/processors/misc/KrComplexLineShape.py b/mermithid/processors/misc/KrComplexLineShape.py
index bba80328..aab299f1 100644
--- a/mermithid/processors/misc/KrComplexLineShape.py
+++ b/mermithid/processors/misc/KrComplexLineShape.py
@@ -84,6 +84,7 @@ def InternalRun(self):
         kr17kev_in_hz = guess*(bins[1]-bins[0])+bins[0]
         #self.B_field = B(17.8, kr17kev_in_hz + 0)
         if self.fix_scatter_proportion == True:
+            logger.info('Scatter proportion fixed')
             self.results = self.fit_data_1(freq_bins, data_hist_freq)
         else:
             self.results = self.fit_data(freq_bins, data_hist_freq)
@@ -127,7 +128,6 @@ def normalize(self, f):
     def single_scatter_f(self, gas_type):
         energy_loss_array = self.std_eV_array()
         f = 0 * energy_loss_array
-
         input_filename = self.path_to_osc_strengths_files + gas_type + "OscillatorStrength.txt"
         energy_fOsc = ComplexLineShapeUtilities.read_oscillator_str_file(input_filename)
         fData = interpolate.interp1d(energy_fOsc[0], energy_fOsc[1], kind='linear')
@@ -151,7 +151,7 @@ def another_scatter(self, input_spectrum, gas_type):
         return f_normed
 
     # Convolves the scatter functions and saves
-    # the results to a .npy file.    
+    # the results to a .npy file.
     def generate_scatter_convolution_file(self):
         t = time.time()
         scatter_spectra_single_gas = {}
@@ -179,7 +179,7 @@ def generate_scatter_convolution_file(self):
                 total_scatter = self.normalize(signal.convolve(H2_scatter, Kr_scatter, mode='same'))
                 scatter_spectra['{}_{}'.format(self.gases[0], self.gases[1])]['{}_{}'.format(str(j).zfill(2), str(i-j).zfill(2))] = total_scatter
         np.save(
-        self.path_to_osc_strengths_files+'scatter_spectra_file/scatter_spectra.npy', 
+        self.path_to_osc_strengths_files+'scatter_spectra_file/scatter_spectra.npy',
         scatter_spectra
         )
         elapsed = time.time() - t
@@ -191,13 +191,13 @@ def generate_scatter_convolution_file(self):
     # If not, this function calls generate_scatter_convolution_file.
     # This function also checks to make sure that the scatter file have the correct
     # number of entries and correct number of points in the SELA, and if not, it generates a fresh file.
-    # When the variable regenerate is set as True, it generates a fresh file   
+    # When the variable regenerate is set as True, it generates a fresh file
     def check_existence_of_scatter_file(self, regenerate = False):
         gases = self.gases
         if regenerate == True:
             logger.info('generate fresh scatter file')
             self.generate_scatter_convolution_file()
-        else: 
+        else:
             stuff_in_dir = os.listdir(self.path_to_osc_strengths_files)
             if 'scatter_spectra_file' not in stuff_in_dir:
                 logger.info('Scatter spectra folder not found, generating')
@@ -209,7 +209,7 @@ def check_existence_of_scatter_file(self, regenerate = False):
                 strippeddirs = [s.strip('\n') for s in directory]
                 if 'scatter_spectra.npy' not in strippeddirs:
                     self.generate_scatter_convolution_file()
-                test_file = self.path_to_osc_strengths_files+'scatter_spectra_file/scatter_spectra.npy' 
+                test_file = self.path_to_osc_strengths_files+'scatter_spectra_file/scatter_spectra.npy'
                 test_dict = np.load(test_file, allow_pickle = True)
                 if list(test_dict.item().keys())[0] != '{}_{}'.format(gases[0], gases[1]):
                     logger.info('first entry not matching, generating fresh files')
@@ -294,12 +294,12 @@ def spectrum_func(self, x_keV, *p0):
         amplitude = p0[2]
         prob_parameter = p0[3]
         scatter_proportion = p0[4]
-    
+
         line_pos_eV = line_pos_keV*1000.
         x_eV_minus_line = x_eV - line_pos_eV
         zero_idx = np.r_[np.where(x_eV_minus_line< en_loss_array_min)[0],np.where(x_eV_minus_line>en_loss_array_max)[0]]
         nonzero_idx = [i for i in range(len(x_keV)) if i not in zero_idx]
-    
+
         full_spectrum = self.make_spectrum(FWHM_G_eV, prob_parameter, scatter_proportion)
         full_spectrum_rev = ComplexLineShapeUtilities.flip_array(full_spectrum)
         f_intermediate[nonzero_idx] = np.interp(x_eV_minus_line[nonzero_idx],en_array_rev,full_spectrum_rev)
@@ -340,8 +340,8 @@ def fit_data(self, freq_bins, data_hist_freq, print_params=True):
         amplitude_guess = np.sum(data_hist)/2
         prob_parameter_guess = 0.5
         scatter_proportion_guess = 0.9
-        p0_guess = [FWHM_guess, line_pos_guess, amplitude_guess, prob_parameter_guess, scatter_proportion_guess] 
-        p0_bounds = ([FWHM_eV_min, line_pos_keV_min, amplitude_min, prob_parameter_min, scatter_proportion_min],  
+        p0_guess = [FWHM_guess, line_pos_guess, amplitude_guess, prob_parameter_guess, scatter_proportion_guess]
+        p0_bounds = ([FWHM_eV_min, line_pos_keV_min, amplitude_min, prob_parameter_min, scatter_proportion_min],
                     [FWHM_eV_max, line_pos_keV_max, amplitude_max, prob_parameter_max, scatter_proportion_max])
         # Actually do the fitting
         params , cov = curve_fit(self.spectrum_func, bins_keV_nonzero, data_hist_nonzero, sigma=data_hist_err, p0=p0_guess, bounds=p0_bounds)
@@ -362,7 +362,7 @@ def fit_data(self, freq_bins, data_hist_freq, print_params=True):
         prob_parameter_fit_err = perr[3]
         scatter_proportion_fit_err = perr[4]
         total_counts_fit_err = amplitude_fit_err
-    
+
         fit = self.spectrum_func(bins_keV,*params)
 
         line_pos_Hz_fit , line_pos_Hz_fit_err = ComplexLineShapeUtilities.energy_guess_to_frequency(line_pos_keV_fit, line_pos_keV_fit_err, self.B_field)
@@ -506,8 +506,8 @@ def fit_data_1(self, freq_bins, data_hist_freq):
         line_pos_guess = bins_keV[np.argmax(data_hist)]
         amplitude_guess = np.sum(data_hist)/2
         prob_parameter_guess = 0.5
-        p0_guess = [FWHM_guess, line_pos_guess, amplitude_guess, prob_parameter_guess] 
-        p0_bounds = ([FWHM_eV_min, line_pos_keV_min, amplitude_min, prob_parameter_min],  
+        p0_guess = [FWHM_guess, line_pos_guess, amplitude_guess, prob_parameter_guess]
+        p0_bounds = ([FWHM_eV_min, line_pos_keV_min, amplitude_min, prob_parameter_min],
                     [FWHM_eV_max, line_pos_keV_max, amplitude_max, prob_parameter_max])
         # Actually do the fitting
         params , cov = curve_fit(self.spectrum_func_1, bins_keV_nonzero, data_hist_nonzero, sigma=data_hist_err, p0=p0_guess, bounds=p0_bounds)
@@ -526,7 +526,7 @@ def fit_data_1(self, freq_bins, data_hist_freq):
         amplitude_fit_err = perr[2]
         prob_parameter_fit_err = perr[3]
         total_counts_fit_err = amplitude_fit_err
-    
+
         fit = self.spectrum_func_1(bins_keV,*params)
 
         line_pos_Hz_fit , line_pos_Hz_fit_err = ComplexLineShapeUtilities.energy_guess_to_frequency(line_pos_keV_fit, line_pos_keV_fit_err, self.B_field)

From 6ec42ecd18ecf2fc2ee9a8fdfcea00bc3078685f Mon Sep 17 00:00:00 2001
From: cclaessens <claessens.christine@hotmail.de>
Date: Fri, 12 Jun 2020 18:19:11 +0200
Subject: [PATCH 02/49] FakeDataGenerator using KrComplexLineShape

Made lineshape base shape configurable in KrComplexLineShape.py.
FakeDataGenerator.py now imports and configures KrComplexLineShape.py
instead of using functions in KrLineshapeFunctions.
---
 mermithid/misc/FakeTritiumDataFunctions.py    | 10 +++--
 .../TritiumSpectrum/FakeDataGenerator.py      | 42 +++++++++++++++----
 .../processors/misc/KrComplexLineShape.py     | 10 +++--
 3 files changed, 47 insertions(+), 15 deletions(-)

diff --git a/mermithid/misc/FakeTritiumDataFunctions.py b/mermithid/misc/FakeTritiumDataFunctions.py
index bbf8edf2..24a9bd52 100644
--- a/mermithid/misc/FakeTritiumDataFunctions.py
+++ b/mermithid/misc/FakeTritiumDataFunctions.py
@@ -20,6 +20,7 @@
 from mermithid.misc.ConversionFunctions import *
 from mermithid.misc.KrLineshapeFunctions import *
 
+
 """
 Constants and functions used by processors/TritiumSpectrum/FakeDataGenerator.py
 """
@@ -245,7 +246,7 @@ def convolved_bkgd_rate(K, Kmin, Kmax, lineshape, ls_params, min_energy, max_ene
 
 #Convolution of signal and lineshape using scipy.signal.convolve
 def convolved_spectral_rate_arrays(K, Q, mnu, Kmin,
-                                   lineshape, ls_params, min_energy, max_energy):
+                                   lineshape, ls_params, min_energy, max_energy, complexLineShape):
     """K is an array-like object
     """
     logger.info('Using scipy convolve')
@@ -263,7 +264,8 @@ def convolved_spectral_rate_arrays(K, Q, mnu, Kmin,
     elif lineshape=='simplified_scattering' or lineshape=='simplified':
         lineshape_rates = simplified_ls(K_lineshape, 0, ls_params[0], ls_params[1], ls_params[2], ls_params[3], ls_params[4], ls_params[5])
     elif lineshape=='detailed_scattering' or lineshape=='detailed':
-        lineshape_rates = spectrum_func(K_lineshape/1000., ls_params[0], 0, ls_params[1], 1)
+
+        lineshape_rates = complexLineShape.spectrum_func_1(K_lineshape/1000., ls_params[0], 0, ls_params[1], 1)
 
     beta_rates = np.zeros(len(K))
     for i,ke in enumerate(K):
@@ -278,7 +280,7 @@ def convolved_spectral_rate_arrays(K, Q, mnu, Kmin,
 
 
 #Convolution of background and lineshape using scipy.signal.convolve
-def convolved_bkgd_rate_arrays(K, Kmin, Kmax, lineshape, ls_params, min_energy, max_energy):
+def convolved_bkgd_rate_arrays(K, Kmin, Kmax, lineshape, ls_params, min_energy, max_energy, complexLineShape):
     """K is an array-like object
     """
     energy_half_range = max(max_energy, abs(min_energy))
@@ -294,7 +296,7 @@ def convolved_bkgd_rate_arrays(K, Kmin, Kmax, lineshape, ls_params, min_energy,
     elif lineshape=='simplified_scattering' or lineshape=='simplified':
         lineshape_rates = simplified_ls(K_lineshape, 0, ls_params[0], ls_params[1], ls_params[2], ls_params[3], ls_params[4], ls_params[5])
     elif lineshape=='detailed_scattering' or lineshape=='detailed':
-        lineshape_rates = spectrum_func(K_lineshape/1000., ls_params[0], 0, ls_params[1], 1)
+        lineshape_rates = complexLineShape.spectrum_func_1(K_lineshape/1000., ls_params[0], 0, ls_params[1], 1)
 
     bkgd_rates = np.full(len(K), bkgd_rate())
     if len(K) < len(K_lineshape):
diff --git a/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py b/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
index 397532df..41a9f568 100644
--- a/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
+++ b/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
@@ -16,7 +16,7 @@
 from morpho.utilities import morphologging, reader
 from morpho.processors import BaseProcessor
 from mermithid.misc.FakeTritiumDataFunctions import *
-
+from mermithid.processors.misc.KrComplexLineShape import KrComplexLineShape
 logger = morphologging.getLogger(__name__)
 
 
@@ -52,7 +52,7 @@ def InternalConfigure(self, params):
         - scattering_sigma [eV]: lineshape parameter - 0-th peak gaussian broadening standard deviation
         - NScatters: lineshape parameter - number of scatters included in lineshape
         - simplified_scattering_path: path to simplified lineshape parameters
-        #- detailed_scattering_path: path to H2 scattering files for detailed lineshape
+        - path_to_detailed_scatter_spectra_dir: path to oscillator and or scatter_spectra_file
         - efficiency_path: path to efficiency vs. frequency (and uncertainties)
         - use_lineshape (boolean): determines whether tritium spectrum is smeared by lineshape.
           If False, it will only be smeared with a Gaussian
@@ -99,7 +99,7 @@ def InternalConfigure(self, params):
 
         #paths
         self.simplified_scattering_path = reader.read_param(params, 'simplified_scattering_path', '/host/input_data/simplified_scattering_params.txt')
-        #self.detailed_scattering_path = reader.read_param(params, 'detailed_scattering_path', None)
+        self.detailed_scatter_spectra_path = reader.read_param(params, 'path_to_detailed_scatter_spectra_dir', '.')
         self.efficiency_path = reader.read_param(params, 'efficiency_path', '')
 
         #options
@@ -108,6 +108,8 @@ def InternalConfigure(self, params):
         self.apply_efficiency = reader.read_param(params, 'apply_efficiency', False)
         self.return_frequency = reader.read_param(params, 'return_frequency', True)
 
+        # nwill be replaced with complex lineshape object if detailed lineshape is used
+        self.complexLineShape = None
 
         # get file content if needed
         # get efficiency dictionary
@@ -125,9 +127,35 @@ def InternalConfigure(self, params):
                                                         self.survival_prob,
                                                         self.NScatters)
             elif self.lineshape=='detailed':
-                if not os.path.exists('./scatter_spectra_files'):
-                    raise IOError('./scatter_spectra_files does not exist')
+                if 'scatter_spectra_file' in self.detailed_scatter_spectra_path:
+                    full_path = self.detailed_scatter_spectra_path
+                    self.detailed_scatter_spectra_path, _ = os.path.split(full_path)
+                logger.info('Path to scatter_spectra_file: {}'.format(self.detailed_scatter_spectra_path))
+
+                if not os.path.exists(os.path.join(self.detailed_scatter_spectra_path, scatter_spectra_file)):
+                    raise IOError('{} does not exist'.format(full_path))
                 self.SimpParams = [self.scattering_sigma*2*math.sqrt(2*math.log(2)), self.survival_prob]
+
+                # complex lineshape object
+                # Setup and configure lineshape processor
+                complexLineShape_config = {
+                    'gases': ["H2","He"],
+                    'max_scatters': 20,
+                    'fix_scatter_proportion': True,
+                    # When fix_scatter_proportion is True, set the scatter proportion for gas1 below
+                    'gas1_scatter_proportion': 0.8,
+                    # This is an important parameter which determines how finely resolved
+                    # the scatter calculations are. 10000 seems to produce a stable fit, with minimal slowdown
+                    'num_points_in_std_array': 10000,
+                    'B_field': self.B_field,
+                    'base_shape': 'lorentzian', # needs to be replaced by dirac
+                    'path_to_osc_strengths_files': self.detailed_scatter_spectra_path
+                }
+
+                # create complex lineshape object and configure
+                self.complexLineShape = KrComplexLineShape("complexLineShape")
+                logger.info('Configuring complex lineshape')
+                self.complexLineShape.Configure(complexLineShape_config)
             else:
                 raise ValueError("'detailed_or_simplified' is neither 'detailed' nor 'simplified'")
 
@@ -257,7 +285,7 @@ def generate_unbinned_data(self, Q_mean, mass, ROIbound, S, B_1kev, nsteps=10**4
         time0 = time.time()
 
         if array_method == True:
-            ratesS = convolved_spectral_rate_arrays(self.Koptions, Q_mean, mass, Kmin, lineshape, params, min_energy, max_energy)
+            ratesS = convolved_spectral_rate_arrays(self.Koptions, Q_mean, mass, Kmin, lineshape, params, min_energy, max_energy, self.complexLineShape)
         else:
             ratesS = [convolved_spectral_rate(K, Q_mean, mass, Kmin, lineshape, params, min_energy, max_energy) for K in self.Koptions]
 
@@ -269,7 +297,7 @@ def generate_unbinned_data(self, Q_mean, mass, ROIbound, S, B_1kev, nsteps=10**4
 
         # background
         if array_method == True:
-            ratesB = convolved_bkgd_rate_arrays(self.Koptions, Kmin, Kmax, lineshape, params, min_energy, max_energy)
+            ratesB = convolved_bkgd_rate_arrays(self.Koptions, Kmin, Kmax, lineshape, params, min_energy, max_energy, self.complexLineShape)
         else:
             ratesB = [convolved_bkgd_rate(K, Kmin, Kmax, lineshape, params, min_energy, max_energy) for K in self.Koptions]
 
diff --git a/mermithid/processors/misc/KrComplexLineShape.py b/mermithid/processors/misc/KrComplexLineShape.py
index aab299f1..d4095af4 100644
--- a/mermithid/processors/misc/KrComplexLineShape.py
+++ b/mermithid/processors/misc/KrComplexLineShape.py
@@ -17,6 +17,7 @@
 B_field: can be put in hand or found by position of the peak of the frequency histogram.
 shake_spectrum_parameters_json_path: path to json file storing shake spectrum parameters.
 path_to_osc_strength_files: path to oscillator strength files.
+base_shape: "shake" or "lorentzian"
 '''
 
 from __future__ import absolute_import
@@ -59,9 +60,10 @@ def InternalConfigure(self, params):
         self.RF_ROI_MIN = reader.read_param(params, 'RF_ROI_MIN', 25850000000.0)
         self.B_field = reader.read_param(params, 'B_field', 0.957810722501)
         self.shake_spectrum_parameters_json_path = reader.read_param(params, 'shake_spectrum_parameters_json_path', 'shake_spectrum_parameters.json')
+        self.base_shape = reader.read_param(params, 'base_shape', 'shake')
         self.path_to_osc_strengths_files = reader.read_param(params, 'path_to_osc_strengths_files', '/host/')
 
-        if not os.path.exists(self.shake_spectrum_parameters_json_path):
+        if self.base_shape=='shake' and not os.path.exists(self.shake_spectrum_parameters_json_path):
             raise IOError('Shake spectrum path does not exist')
         if not os.path.exists(self.path_to_osc_strengths_files):
             raise IOError('Path to osc strengths files does not exist')
@@ -104,7 +106,7 @@ def std_eV_array(self):
     # A lorentzian line centered at 0 eV, with 2.83 eV width on the SELA
     def std_lorenztian_17keV(self):
         x_array = self.std_eV_array()
-        ans = lorentzian(x_array,0,kr_line_width)
+        ans = ComplexLineShapeUtilities.lorentzian(x_array,0,ComplexLineShapeUtilities.kr_17keV_line_width)
         return ans
 
     # A gaussian centered at 0 eV with variable width, on the SELA
@@ -421,7 +423,7 @@ def make_spectrum_1(self, gauss_FWHM_eV, prob_parameter, emitted_peak='shake'):
         p = self.scatter_proportion
         q = 1 - p
         scatter_spectra = np.load(
-        current_path + 'scatter_spectra_file/scatter_spectra.npy', allow_pickle = True
+        os.path.join(current_path,'scatter_spectra_file/scatter_spectra.npy'), allow_pickle = True
         )
         en_array = self.std_eV_array()
         current_full_spectrum = np.zeros(len(en_array))
@@ -476,7 +478,7 @@ def spectrum_func_1(self, x_keV, *p0):
         zero_idx = np.r_[np.where(x_eV_minus_line< en_loss_array_min)[0],np.where(x_eV_minus_line>en_loss_array_max)[0]]
         nonzero_idx = [i for i in range(len(x_keV)) if i not in zero_idx]
 
-        full_spectrum = self.make_spectrum_1(FWHM_G_eV, prob_parameter,)
+        full_spectrum = self.make_spectrum_1(FWHM_G_eV, prob_parameter,emitted_peak=self.base_shape)
         full_spectrum_rev = ComplexLineShapeUtilities.flip_array(full_spectrum)
         f_intermediate[nonzero_idx] = np.interp(x_eV_minus_line[nonzero_idx],en_array_rev,full_spectrum_rev)
         f[nonzero_idx] += amplitude*f_intermediate[nonzero_idx]/np.sum(f_intermediate[nonzero_idx])

From a6ec33bb53b5669d7cdaeffcc1c33c6017694b70 Mon Sep 17 00:00:00 2001
From: cclaessens <claessens.christine@hotmail.de>
Date: Fri, 12 Jun 2020 18:26:50 +0200
Subject: [PATCH 03/49] added scatter_proportion configurable

---
 .../processors/TritiumSpectrum/FakeDataGenerator.py      | 9 +++++----
 tests/Fake_data_generator_test.py                        | 2 ++
 2 files changed, 7 insertions(+), 4 deletions(-)

diff --git a/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py b/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
index 41a9f568..81a70535 100644
--- a/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
+++ b/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
@@ -51,6 +51,7 @@ def InternalConfigure(self, params):
         - survival_prob: lineshape parameter - ratio of n+t/nth peak
         - scattering_sigma [eV]: lineshape parameter - 0-th peak gaussian broadening standard deviation
         - NScatters: lineshape parameter - number of scatters included in lineshape
+        - scatter_proportion: fraction of hydrogen in complex lineshape
         - simplified_scattering_path: path to simplified lineshape parameters
         - path_to_detailed_scatter_spectra_dir: path to oscillator and or scatter_spectra_file
         - efficiency_path: path to efficiency vs. frequency (and uncertainties)
@@ -91,11 +92,11 @@ def InternalConfigure(self, params):
         self.err_from_B = reader.read_param(params, 'err_from_B', 0.) #In eV, kinetic energy error from f_c --> K conversion
 
 
-        #Simplified scattering model parameters
+        #Scattering model parameters
         self.survival_prob = reader.read_param(params, 'survival_prob', 0.77)
         self.scattering_sigma = reader.read_param(params, 'scattering_sigma', 18.6)
         self.NScatters = reader.read_param(params, 'NScatters', 20)
-
+        self.scatter_proportion = reader.read_param(params, 'scatter_proportion', 0.8)
 
         #paths
         self.simplified_scattering_path = reader.read_param(params, 'simplified_scattering_path', '/host/input_data/simplified_scattering_params.txt')
@@ -140,10 +141,10 @@ def InternalConfigure(self, params):
                 # Setup and configure lineshape processor
                 complexLineShape_config = {
                     'gases': ["H2","He"],
-                    'max_scatters': 20,
+                    'max_scatters': self.NScatters,
                     'fix_scatter_proportion': True,
                     # When fix_scatter_proportion is True, set the scatter proportion for gas1 below
-                    'gas1_scatter_proportion': 0.8,
+                    'gas1_scatter_proportion': self.scatter_proportion,
                     # This is an important parameter which determines how finely resolved
                     # the scatter calculations are. 10000 seems to produce a stable fit, with minimal slowdown
                     'num_points_in_std_array': 10000,
diff --git a/tests/Fake_data_generator_test.py b/tests/Fake_data_generator_test.py
index 22a41025..2f0d0044 100644
--- a/tests/Fake_data_generator_test.py
+++ b/tests/Fake_data_generator_test.py
@@ -20,11 +20,13 @@ def test_data_generation(self):
             "apply_efficiency": False,
             "efficiency_path": "./combined_energy_corrected_eff_at_quad_trap_frequencies.json",
             "simplified_lineshape_path": None,
+            "path_to_detailed_scatter_spectra_dir": '/host',
             "detailed_or_simplified_lineshape": "detailed", #"simplified" or "detailed"
             "use_lineshape": False, # if False only gaussian smearing is applied
             "return_frequency": True,
             "scattering_sigma": 18.6, # only used if use_lineshape = True
             "scattering_prob": 0.77, # only used if use_lineshape = True
+            "scatter_proportion": 0.8, # only used if use_lineshape = True and lineshape = detailed
             "B_field": 0.9578186017836624,
             "S": 4500, # number of tritium events
             "n_steps": 1000, # stepsize for pseudo continuous data is: (Kmax_eff-Kmin_eff)/nsteps

From 7130957eeb1d0984bc84b2947284f2dab8bd8384 Mon Sep 17 00:00:00 2001
From: cclaessens <claessens.christine@hotmail.de>
Date: Fri, 12 Jun 2020 18:32:06 +0200
Subject: [PATCH 04/49] small fixes

---
 mermithid/processors/TritiumSpectrum/FakeDataGenerator.py | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py b/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
index 81a70535..aaeb3123 100644
--- a/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
+++ b/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
@@ -128,13 +128,19 @@ def InternalConfigure(self, params):
                                                         self.survival_prob,
                                                         self.NScatters)
             elif self.lineshape=='detailed':
+                # check path exists
                 if 'scatter_spectra_file' in self.detailed_scatter_spectra_path:
                     full_path = self.detailed_scatter_spectra_path
                     self.detailed_scatter_spectra_path, _ = os.path.split(full_path)
+                else:
+                    full_path = os.path.join(self.detailed_scatter_spectra_path, 'scatter_spectra_file')
+
                 logger.info('Path to scatter_spectra_file: {}'.format(self.detailed_scatter_spectra_path))
 
-                if not os.path.exists(os.path.join(self.detailed_scatter_spectra_path, scatter_spectra_file)):
+                if not os.path.exists(full_path):
                     raise IOError('{} does not exist'.format(full_path))
+
+                # lineshape params
                 self.SimpParams = [self.scattering_sigma*2*math.sqrt(2*math.log(2)), self.survival_prob]
 
                 # complex lineshape object

From 8897a184b3601ef65a608e58629aa3763d29bf39 Mon Sep 17 00:00:00 2001
From: cclaessens <claessens.christine@hotmail.de>
Date: Fri, 12 Jun 2020 19:28:53 +0200
Subject: [PATCH 05/49] scatter_spectra.npy file does not need to pre-exist.
 KrLineshapeFunctions no longer used.

---
 mermithid/misc/FakeTritiumDataFunctions.py               | 6 +++---
 .../processors/TritiumSpectrum/FakeDataGenerator.py      | 9 ++++-----
 mermithid/processors/misc/KrComplexLineShape.py          | 6 +++---
 3 files changed, 10 insertions(+), 11 deletions(-)

diff --git a/mermithid/misc/FakeTritiumDataFunctions.py b/mermithid/misc/FakeTritiumDataFunctions.py
index 24a9bd52..6e565ed9 100644
--- a/mermithid/misc/FakeTritiumDataFunctions.py
+++ b/mermithid/misc/FakeTritiumDataFunctions.py
@@ -18,7 +18,7 @@
 
 from mermithid.misc.Constants import *
 from mermithid.misc.ConversionFunctions import *
-from mermithid.misc.KrLineshapeFunctions import *
+#from mermithid.misc.KrLineshapeFunctions import *
 
 
 """
@@ -171,8 +171,8 @@ def bkgd_rate():
 
 
 ##Lineshape option: In this case, simply a Gaussian
-#def gaussian(x,a):
-#    return 1/((2.*np.pi)**0.5*a[0])*(np.exp(-0.5*((x-a[1])/a[0])**2))
+def gaussian(x,a):
+    return 1/((2.*np.pi)**0.5*a[0])*(np.exp(-0.5*((x-a[1])/a[0])**2))
 
 
 #Normalized simplified linesahape with scattering
diff --git a/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py b/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
index aaeb3123..c08b5954 100644
--- a/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
+++ b/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
@@ -137,13 +137,11 @@ def InternalConfigure(self, params):
 
                 logger.info('Path to scatter_spectra_file: {}'.format(self.detailed_scatter_spectra_path))
 
-                if not os.path.exists(full_path):
-                    raise IOError('{} does not exist'.format(full_path))
 
                 # lineshape params
                 self.SimpParams = [self.scattering_sigma*2*math.sqrt(2*math.log(2)), self.survival_prob]
 
-                # complex lineshape object
+
                 # Setup and configure lineshape processor
                 complexLineShape_config = {
                     'gases': ["H2","He"],
@@ -158,11 +156,12 @@ def InternalConfigure(self, params):
                     'base_shape': 'lorentzian', # needs to be replaced by dirac
                     'path_to_osc_strengths_files': self.detailed_scatter_spectra_path
                 }
-
-                # create complex lineshape object and configure
+                logger.info('Setting up complex lineshape object')
                 self.complexLineShape = KrComplexLineShape("complexLineShape")
                 logger.info('Configuring complex lineshape')
                 self.complexLineShape.Configure(complexLineShape_config)
+                logger.info('Checking existence of scatter spectra files')
+                self.complexLineShape.check_existence_of_scatter_file()
             else:
                 raise ValueError("'detailed_or_simplified' is neither 'detailed' nor 'simplified'")
 
diff --git a/mermithid/processors/misc/KrComplexLineShape.py b/mermithid/processors/misc/KrComplexLineShape.py
index d4095af4..0c45ce64 100644
--- a/mermithid/processors/misc/KrComplexLineShape.py
+++ b/mermithid/processors/misc/KrComplexLineShape.py
@@ -130,7 +130,7 @@ def normalize(self, f):
     def single_scatter_f(self, gas_type):
         energy_loss_array = self.std_eV_array()
         f = 0 * energy_loss_array
-        input_filename = self.path_to_osc_strengths_files + gas_type + "OscillatorStrength.txt"
+        input_filename = os.path.join(self.path_to_osc_strengths_files,  "{}OscillatorStrength.txt".format(gas_type))
         energy_fOsc = ComplexLineShapeUtilities.read_oscillator_str_file(input_filename)
         fData = interpolate.interp1d(energy_fOsc[0], energy_fOsc[1], kind='linear')
         for i in range(len(energy_loss_array)):
@@ -181,7 +181,7 @@ def generate_scatter_convolution_file(self):
                 total_scatter = self.normalize(signal.convolve(H2_scatter, Kr_scatter, mode='same'))
                 scatter_spectra['{}_{}'.format(self.gases[0], self.gases[1])]['{}_{}'.format(str(j).zfill(2), str(i-j).zfill(2))] = total_scatter
         np.save(
-        self.path_to_osc_strengths_files+'scatter_spectra_file/scatter_spectra.npy',
+        os.path.join(self.path_to_osc_strengths_files, 'scatter_spectra_file/scatter_spectra.npy'),
         scatter_spectra
         )
         elapsed = time.time() - t
@@ -211,7 +211,7 @@ def check_existence_of_scatter_file(self, regenerate = False):
                 strippeddirs = [s.strip('\n') for s in directory]
                 if 'scatter_spectra.npy' not in strippeddirs:
                     self.generate_scatter_convolution_file()
-                test_file = self.path_to_osc_strengths_files+'scatter_spectra_file/scatter_spectra.npy'
+                test_file = os.path.join(self.path_to_osc_strengths_files, 'scatter_spectra_file/scatter_spectra.npy')
                 test_dict = np.load(test_file, allow_pickle = True)
                 if list(test_dict.item().keys())[0] != '{}_{}'.format(gases[0], gases[1]):
                     logger.info('first entry not matching, generating fresh files')

From 754024b2e9c4c1d1c6a44df22a07e0069096566d Mon Sep 17 00:00:00 2001
From: Talia Weiss <>
Date: Fri, 12 Jun 2020 12:55:50 -0500
Subject: [PATCH 06/49] Added delta base shape option; removed extra input to
 spectrum_func_1

---
 mermithid/misc/FakeTritiumDataFunctions.py    |  4 ++--
 .../TritiumSpectrum/FakeDataGenerator.py      |  4 ++--
 .../processors/misc/KrComplexLineShape.py     | 23 +++++++++++++++++--
 3 files changed, 25 insertions(+), 6 deletions(-)

diff --git a/mermithid/misc/FakeTritiumDataFunctions.py b/mermithid/misc/FakeTritiumDataFunctions.py
index 24a9bd52..7dfbf064 100644
--- a/mermithid/misc/FakeTritiumDataFunctions.py
+++ b/mermithid/misc/FakeTritiumDataFunctions.py
@@ -265,7 +265,7 @@ def convolved_spectral_rate_arrays(K, Q, mnu, Kmin,
         lineshape_rates = simplified_ls(K_lineshape, 0, ls_params[0], ls_params[1], ls_params[2], ls_params[3], ls_params[4], ls_params[5])
     elif lineshape=='detailed_scattering' or lineshape=='detailed':
 
-        lineshape_rates = complexLineShape.spectrum_func_1(K_lineshape/1000., ls_params[0], 0, ls_params[1], 1)
+        lineshape_rates = complexLineShape.spectrum_func_1(K_lineshape/1000., ls_params[0], 0, ls_params[1])
 
     beta_rates = np.zeros(len(K))
     for i,ke in enumerate(K):
@@ -296,7 +296,7 @@ def convolved_bkgd_rate_arrays(K, Kmin, Kmax, lineshape, ls_params, min_energy,
     elif lineshape=='simplified_scattering' or lineshape=='simplified':
         lineshape_rates = simplified_ls(K_lineshape, 0, ls_params[0], ls_params[1], ls_params[2], ls_params[3], ls_params[4], ls_params[5])
     elif lineshape=='detailed_scattering' or lineshape=='detailed':
-        lineshape_rates = complexLineShape.spectrum_func_1(K_lineshape/1000., ls_params[0], 0, ls_params[1], 1)
+        lineshape_rates = complexLineShape.spectrum_func_1(K_lineshape/1000., ls_params[0], 0, ls_params[1])
 
     bkgd_rates = np.full(len(K), bkgd_rate())
     if len(K) < len(K_lineshape):
diff --git a/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py b/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
index aaeb3123..661dfb98 100644
--- a/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
+++ b/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
@@ -96,7 +96,7 @@ def InternalConfigure(self, params):
         self.survival_prob = reader.read_param(params, 'survival_prob', 0.77)
         self.scattering_sigma = reader.read_param(params, 'scattering_sigma', 18.6)
         self.NScatters = reader.read_param(params, 'NScatters', 20)
-        self.scatter_proportion = reader.read_param(params, 'scatter_proportion', 0.8)
+        self.scatter_proportion = reader.read_param(params, 'scatter_proportion', 1.0)
 
         #paths
         self.simplified_scattering_path = reader.read_param(params, 'simplified_scattering_path', '/host/input_data/simplified_scattering_params.txt')
@@ -155,7 +155,7 @@ def InternalConfigure(self, params):
                     # the scatter calculations are. 10000 seems to produce a stable fit, with minimal slowdown
                     'num_points_in_std_array': 10000,
                     'B_field': self.B_field,
-                    'base_shape': 'lorentzian', # needs to be replaced by dirac
+                    'base_shape': 'dirac',
                     'path_to_osc_strengths_files': self.detailed_scatter_spectra_path
                 }
 
diff --git a/mermithid/processors/misc/KrComplexLineShape.py b/mermithid/processors/misc/KrComplexLineShape.py
index d4095af4..1619c4d0 100644
--- a/mermithid/processors/misc/KrComplexLineShape.py
+++ b/mermithid/processors/misc/KrComplexLineShape.py
@@ -17,7 +17,7 @@
 B_field: can be put in hand or found by position of the peak of the frequency histogram.
 shake_spectrum_parameters_json_path: path to json file storing shake spectrum parameters.
 path_to_osc_strength_files: path to oscillator strength files.
-base_shape: "shake" or "lorentzian"
+base_shape: "shake", "lorentzian" or "dirac"
 '''
 
 from __future__ import absolute_import
@@ -109,6 +109,21 @@ def std_lorenztian_17keV(self):
         ans = ComplexLineShapeUtilities.lorentzian(x_array,0,ComplexLineShapeUtilities.kr_17keV_line_width)
         return ans
 
+    #A Dirac delta functin
+    def std_dirac():
+        x_array = std_eV_array()
+        ans = np.zeros(len(x_array))
+        min_x = np.min(np.abs(x_array))
+        ans[np.abs(x_array)==min_x] = 1.
+        logger.warning('Spectrum will be shifted by lineshape by {} eV'.format(min_x))
+        if min_x > 0.1:
+            logger.warning('Lineshape will shift spectrum by > 0.1 eV')
+        if min_x > 1.:
+            logger.warning('Lineshape will shift spectrum by > 1 eV')
+            raise ValueError('problem with std_eV_array()')
+        return ans
+
+
     # A gaussian centered at 0 eV with variable width, on the SELA
     def std_gaussian(self, sigma):
         x_array = self.std_eV_array()
@@ -251,6 +266,8 @@ def make_spectrum(self, gauss_FWHM_eV, prob_parameter, scatter_proportion, emitt
             current_working_spectrum = self.std_lorenztian_17keV()
         elif emitted_peak == 'shake':
             current_working_spectrum = self.shakeSpectrumClassInstance.shake_spectrum()
+        elif emitted_peak == 'dirac':
+            current_working_spectrum = std_dirac()
         current_working_spectrum = self.convolve_gaussian(current_working_spectrum, gauss_FWHM_eV)
         zeroth_order_peak = current_working_spectrum
         current_full_spectrum += current_working_spectrum
@@ -431,6 +448,8 @@ def make_spectrum_1(self, gauss_FWHM_eV, prob_parameter, emitted_peak='shake'):
             current_working_spectrum = self.std_lorenztian_17keV()
         elif emitted_peak == 'shake':
             current_working_spectrum = self.shakeSpectrumClassInstance.shake_spectrum()
+        elif emitted_peak == 'dirac':
+            current_working_spectrum = std_dirac()
         current_working_spectrum = self.convolve_gaussian(current_working_spectrum, gauss_FWHM_eV)
         zeroth_order_peak = current_working_spectrum
         current_full_spectrum += current_working_spectrum
@@ -570,4 +589,4 @@ def fit_data_1(self, freq_bins, data_hist_freq):
         'amplitude_fit_err': amplitude_fit_err,
         'data_hist_freq': data_hist_freq
         }
-        return dictionary_of_fit_results
\ No newline at end of file
+        return dictionary_of_fit_results

From 0c8e97f74c9b7189d0eb1ea19c7f3431d8066e9a Mon Sep 17 00:00:00 2001
From: cclaessens <claessens.christine@hotmail.de>
Date: Fri, 12 Jun 2020 20:14:58 +0200
Subject: [PATCH 07/49] some fixes

---
 mermithid/misc/FakeTritiumDataFunctions.py      | 4 ++--
 mermithid/processors/misc/KrComplexLineShape.py | 8 ++++----
 2 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/mermithid/misc/FakeTritiumDataFunctions.py b/mermithid/misc/FakeTritiumDataFunctions.py
index 38b0dc50..343ee319 100644
--- a/mermithid/misc/FakeTritiumDataFunctions.py
+++ b/mermithid/misc/FakeTritiumDataFunctions.py
@@ -265,7 +265,7 @@ def convolved_spectral_rate_arrays(K, Q, mnu, Kmin,
         lineshape_rates = simplified_ls(K_lineshape, 0, ls_params[0], ls_params[1], ls_params[2], ls_params[3], ls_params[4], ls_params[5])
     elif lineshape=='detailed_scattering' or lineshape=='detailed':
 
-        lineshape_rates = complexLineShape.spectrum_func_1(K_lineshape/1000., ls_params[0], 0, ls_params[1])
+        lineshape_rates = complexLineShape.spectrum_func_1(K_lineshape/1000., ls_params[0], 0, 1, ls_params[1])
 
     beta_rates = np.zeros(len(K))
     for i,ke in enumerate(K):
@@ -296,7 +296,7 @@ def convolved_bkgd_rate_arrays(K, Kmin, Kmax, lineshape, ls_params, min_energy,
     elif lineshape=='simplified_scattering' or lineshape=='simplified':
         lineshape_rates = simplified_ls(K_lineshape, 0, ls_params[0], ls_params[1], ls_params[2], ls_params[3], ls_params[4], ls_params[5])
     elif lineshape=='detailed_scattering' or lineshape=='detailed':
-        lineshape_rates = complexLineShape.spectrum_func_1(K_lineshape/1000., ls_params[0], 0, ls_params[1])
+        lineshape_rates = complexLineShape.spectrum_func_1(K_lineshape/1000., ls_params[0], 0, 1, ls_params[1])
 
     bkgd_rates = np.full(len(K), bkgd_rate())
     if len(K) < len(K_lineshape):
diff --git a/mermithid/processors/misc/KrComplexLineShape.py b/mermithid/processors/misc/KrComplexLineShape.py
index 0de43030..bc3db21b 100644
--- a/mermithid/processors/misc/KrComplexLineShape.py
+++ b/mermithid/processors/misc/KrComplexLineShape.py
@@ -110,8 +110,8 @@ def std_lorenztian_17keV(self):
         return ans
 
     #A Dirac delta functin
-    def std_dirac():
-        x_array = std_eV_array()
+    def std_dirac(self):
+        x_array = self.std_eV_array()
         ans = np.zeros(len(x_array))
         min_x = np.min(np.abs(x_array))
         ans[np.abs(x_array)==min_x] = 1.
@@ -218,7 +218,7 @@ def check_existence_of_scatter_file(self, regenerate = False):
             stuff_in_dir = os.listdir(self.path_to_osc_strengths_files)
             if 'scatter_spectra_file' not in stuff_in_dir:
                 logger.info('Scatter spectra folder not found, generating')
-                os.mkdir(self.path_to_osc_strengths_files+'scatter_spectra_file')
+                os.mkdir(os.path.join(self.path_to_osc_strengths_files,'scatter_spectra_file'))
                 time.sleep(2)
                 self.generate_scatter_convolution_file()
             else:
@@ -449,7 +449,7 @@ def make_spectrum_1(self, gauss_FWHM_eV, prob_parameter, emitted_peak='shake'):
         elif emitted_peak == 'shake':
             current_working_spectrum = self.shakeSpectrumClassInstance.shake_spectrum()
         elif emitted_peak == 'dirac':
-            current_working_spectrum = std_dirac()
+            current_working_spectrum = self.std_dirac()
         current_working_spectrum = self.convolve_gaussian(current_working_spectrum, gauss_FWHM_eV)
         zeroth_order_peak = current_working_spectrum
         current_full_spectrum += current_working_spectrum

From c10eccdfd87ef9c64f1fc93170a6c457b959aaca Mon Sep 17 00:00:00 2001
From: Talia Weiss <>
Date: Fri, 12 Jun 2020 14:16:58 -0500
Subject: [PATCH 08/49] Added option to normalize complex lineshape

---
 mermithid/misc/FakeTritiumDataFunctions.py    |  4 +--
 .../TritiumSpectrum/FakeDataGenerator.py      |  1 +
 .../processors/misc/KrComplexLineShape.py     | 33 ++++++++++++++++---
 3 files changed, 32 insertions(+), 6 deletions(-)

diff --git a/mermithid/misc/FakeTritiumDataFunctions.py b/mermithid/misc/FakeTritiumDataFunctions.py
index 38b0dc50..343ee319 100644
--- a/mermithid/misc/FakeTritiumDataFunctions.py
+++ b/mermithid/misc/FakeTritiumDataFunctions.py
@@ -265,7 +265,7 @@ def convolved_spectral_rate_arrays(K, Q, mnu, Kmin,
         lineshape_rates = simplified_ls(K_lineshape, 0, ls_params[0], ls_params[1], ls_params[2], ls_params[3], ls_params[4], ls_params[5])
     elif lineshape=='detailed_scattering' or lineshape=='detailed':
 
-        lineshape_rates = complexLineShape.spectrum_func_1(K_lineshape/1000., ls_params[0], 0, ls_params[1])
+        lineshape_rates = complexLineShape.spectrum_func_1(K_lineshape/1000., ls_params[0], 0, 1, ls_params[1])
 
     beta_rates = np.zeros(len(K))
     for i,ke in enumerate(K):
@@ -296,7 +296,7 @@ def convolved_bkgd_rate_arrays(K, Kmin, Kmax, lineshape, ls_params, min_energy,
     elif lineshape=='simplified_scattering' or lineshape=='simplified':
         lineshape_rates = simplified_ls(K_lineshape, 0, ls_params[0], ls_params[1], ls_params[2], ls_params[3], ls_params[4], ls_params[5])
     elif lineshape=='detailed_scattering' or lineshape=='detailed':
-        lineshape_rates = complexLineShape.spectrum_func_1(K_lineshape/1000., ls_params[0], 0, ls_params[1])
+        lineshape_rates = complexLineShape.spectrum_func_1(K_lineshape/1000., ls_params[0], 0, 1, ls_params[1])
 
     bkgd_rates = np.full(len(K), bkgd_rate())
     if len(K) < len(K_lineshape):
diff --git a/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py b/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
index f498b379..60110453 100644
--- a/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
+++ b/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
@@ -154,6 +154,7 @@ def InternalConfigure(self, params):
                     'num_points_in_std_array': 10000,
                     'B_field': self.B_field,
                     'base_shape': 'dirac',
+                    'normalize_lineshape': True,
                     'path_to_osc_strengths_files': self.detailed_scatter_spectra_path
                 }
                 logger.info('Setting up complex lineshape object')
diff --git a/mermithid/processors/misc/KrComplexLineShape.py b/mermithid/processors/misc/KrComplexLineShape.py
index 0de43030..d1f1e851 100644
--- a/mermithid/processors/misc/KrComplexLineShape.py
+++ b/mermithid/processors/misc/KrComplexLineShape.py
@@ -54,6 +54,7 @@ def InternalConfigure(self, params):
         self.fix_scatter_proportion = reader.read_param(params, 'fix_scatter_proportion', True)
         if self.fix_scatter_proportion == True:
             self.scatter_proportion = reader.read_param(params, 'gas1_scatter_proportion', 0.8)
+        self.normalize_lineshape = reader.read_params('normalize_lineshape', False)
         # This is an important parameter which determines how finely resolved
         # the scatter calculations are. 10000 seems to produce a stable fit, with minimal slowdown
         self.num_points_in_std_array = reader.read_param(params, 'num_points_in_std_array', 10000)
@@ -267,10 +268,12 @@ def make_spectrum(self, gauss_FWHM_eV, prob_parameter, scatter_proportion, emitt
         elif emitted_peak == 'shake':
             current_working_spectrum = self.shakeSpectrumClassInstance.shake_spectrum()
         elif emitted_peak == 'dirac':
-            current_working_spectrum = std_dirac()
+            current_working_spectrum = self.std_dirac()
         current_working_spectrum = self.convolve_gaussian(current_working_spectrum, gauss_FWHM_eV)
         zeroth_order_peak = current_working_spectrum
         current_full_spectrum += current_working_spectrum
+        if self.normalize_lineshape==True:
+            norm = 1
         for n in range(1, max_comprehensive_scatters + 1):
             for r in range(0, n + 1):
                 current_working_spectrum = \
@@ -280,6 +283,9 @@ def make_spectrum(self, gauss_FWHM_eV, prob_parameter, scatter_proportion, emitt
                 self.normalize(signal.convolve(zeroth_order_peak, current_working_spectrum, mode='same'))
                 current_full_spectrum += current_working_spectrum*comb(n, r)\
                 *(prob_parameter*p)**(r)*(prob_parameter*q)**(n-r)
+                if self.normalize_lineshape==True:
+                    norm += comb(n, r)\
+                    *(prob_parameter*p)**(r)*(prob_parameter*q)**(n-r)
 
         for n in range(max_comprehensive_scatters + 1, max_scatters + 1):
             current_working_spectrum = \
@@ -288,6 +294,8 @@ def make_spectrum(self, gauss_FWHM_eV, prob_parameter, scatter_proportion, emitt
             current_working_spectrum = \
             self.normalize(signal.convolve(zeroth_order_peak, current_working_spectrum, mode='same'))
             current_full_spectrum += current_working_spectrum*(prob_parameter*p)**(n)
+            if self.normalize_lineshape==True:
+                norm += (prob_parameter*p)**(n)
 
             current_working_spectrum = \
             scatter_spectra.item()['{}_{}'.format(gases[0], gases[1])]\
@@ -295,7 +303,12 @@ def make_spectrum(self, gauss_FWHM_eV, prob_parameter, scatter_proportion, emitt
             current_working_spectrum = \
             self.normalize(signal.convolve(zeroth_order_peak, current_working_spectrum, mode='same'))
             current_full_spectrum += current_working_spectrum*(prob_parameter*q)**(n)
-        return current_full_spectrum
+            if self.normalize_lineshape==True:
+                norm += (prob_parameter*q)**(n)
+        if self.normalize_lineshape==True:
+            return current_full_spectrum/norm
+        else:
+            return current_full_spectrum
 
     # Produces a spectrum in real energy that can now be evaluated off of the SELA.
     #def spectrum_func(x_keV,FWHM_G_eV,line_pos_keV,scatter_prob,amplitude):
@@ -449,10 +462,12 @@ def make_spectrum_1(self, gauss_FWHM_eV, prob_parameter, emitted_peak='shake'):
         elif emitted_peak == 'shake':
             current_working_spectrum = self.shakeSpectrumClassInstance.shake_spectrum()
         elif emitted_peak == 'dirac':
-            current_working_spectrum = std_dirac()
+            current_working_spectrum = self.std_dirac()
         current_working_spectrum = self.convolve_gaussian(current_working_spectrum, gauss_FWHM_eV)
         zeroth_order_peak = current_working_spectrum
         current_full_spectrum += current_working_spectrum
+        if self.normalize_lineshape==True:
+            norm = 1
         for n in range(1, max_comprehensive_scatters + 1):
             for r in range(0, n + 1):
                 current_working_spectrum = \
@@ -462,6 +477,9 @@ def make_spectrum_1(self, gauss_FWHM_eV, prob_parameter, emitted_peak='shake'):
                 self.normalize(signal.convolve(zeroth_order_peak, current_working_spectrum, mode='same'))
                 current_full_spectrum += current_working_spectrum*comb(n, r)\
                 *(prob_parameter*p)**(r)*(prob_parameter*q)**(n-r)
+                if self.normalize_lineshape==True:
+                    norm += comb(n, r)\
+                    *(prob_parameter*p)**(r)*(prob_parameter*q)**(n-r)
         for n in range(max_comprehensive_scatters + 1, max_scatters + 1):
             current_working_spectrum = \
             scatter_spectra.item()['{}_{}'.format(gases[0], gases[1])]\
@@ -469,6 +487,8 @@ def make_spectrum_1(self, gauss_FWHM_eV, prob_parameter, emitted_peak='shake'):
             current_working_spectrum = \
             self.normalize(signal.convolve(zeroth_order_peak, current_working_spectrum, mode='same'))
             current_full_spectrum += current_working_spectrum*(prob_parameter*p)**(n)
+            if self.normalize_lineshape==True:
+                norm += (prob_parameter*p)**(n)
 
             current_working_spectrum = \
             scatter_spectra.item()['{}_{}'.format(gases[0], gases[1])]\
@@ -476,7 +496,12 @@ def make_spectrum_1(self, gauss_FWHM_eV, prob_parameter, emitted_peak='shake'):
             current_working_spectrum = \
             self.normalize(signal.convolve(zeroth_order_peak, current_working_spectrum, mode='same'))
             current_full_spectrum += current_working_spectrum*(prob_parameter*q)**(n)
-        return current_full_spectrum
+            if self.normalize_lineshape==True:
+                norm += (prob_parameter*q)**(n)
+        if self.normalize_lineshape==True:
+            return current_full_spectrum/norm
+        else:
+            return current_full_spectrum
 
     def spectrum_func_1(self, x_keV, *p0):
         x_eV = x_keV*1000.

From f9846394afd4d71adf84999db920ce39e7315773 Mon Sep 17 00:00:00 2001
From: Talia Weiss <>
Date: Fri, 12 Jun 2020 14:42:01 -0500
Subject: [PATCH 09/49] Fixed errors

---
 mermithid/processors/misc/KrComplexLineShape.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/mermithid/processors/misc/KrComplexLineShape.py b/mermithid/processors/misc/KrComplexLineShape.py
index 2e238581..99872568 100644
--- a/mermithid/processors/misc/KrComplexLineShape.py
+++ b/mermithid/processors/misc/KrComplexLineShape.py
@@ -54,7 +54,7 @@ def InternalConfigure(self, params):
         self.fix_scatter_proportion = reader.read_param(params, 'fix_scatter_proportion', True)
         if self.fix_scatter_proportion == True:
             self.scatter_proportion = reader.read_param(params, 'gas1_scatter_proportion', 0.8)
-        self.normalize_lineshape = reader.read_params('normalize_lineshape', False)
+        self.normalize_lineshape = reader.read_param(params, 'normalize_lineshape', False)
         # This is an important parameter which determines how finely resolved
         # the scatter calculations are. 10000 seems to produce a stable fit, with minimal slowdown
         self.num_points_in_std_array = reader.read_param(params, 'num_points_in_std_array', 10000)

From e733607476f8117f2c8852ee73df663ccc9d40fb Mon Sep 17 00:00:00 2001
From: cclaessens <claessens.christine@hotmail.de>
Date: Sat, 13 Jun 2020 20:22:35 +0200
Subject: [PATCH 10/49] iminuit binned data fitter

---
 mermithid/processors/misc/BinnedDataFitter.py | 152 ++++++++++++++++++
 mermithid/processors/misc/__init__.py         |   1 +
 tests/Misc_test.py                            | 114 ++++++++++---
 3 files changed, 246 insertions(+), 21 deletions(-)
 create mode 100644 mermithid/processors/misc/BinnedDataFitter.py

diff --git a/mermithid/processors/misc/BinnedDataFitter.py b/mermithid/processors/misc/BinnedDataFitter.py
new file mode 100644
index 00000000..7d5f01ce
--- /dev/null
+++ b/mermithid/processors/misc/BinnedDataFitter.py
@@ -0,0 +1,152 @@
+'''
+Author: C. Claessens
+Date:6/12/2020
+Description: Configure with arbitrary histogram and pdf.
+Minimizes poisson loglikelihood using iMinuit.
+
+'''
+
+from __future__ import absolute_import
+
+import numpy as np
+import scipy
+import sys
+import json
+from iminuit import Minuit
+from morpho.utilities import morphologging, reader
+from morpho.processors import BaseProcessor
+
+logger = morphologging.getLogger(__name__)
+
+
+
+__all__ = []
+__all__.append(__name__)
+
+class BinnedDataFitter(BaseProcessor):
+    '''
+    Processor that
+    Args:
+        variables
+        parameter_names
+        initial_values
+        limits
+        fixed
+        bins
+        print_level
+        constrained_parameter_indices
+        constrained_parameter_means
+        constrained_parameter_widths
+
+    Inputs:
+        data:
+    Output:
+        result: dictionary containing
+    '''
+    def InternalConfigure(self, params):
+        '''
+        Configure
+        '''
+
+        # Read other parameters
+        self.namedata = reader.read_param(params, 'variables', "required")
+        self.parameter_names = reader.read_param(params, 'parameter_names', ['A', 'mu', 'sigma'])
+        self.initial_values = reader.read_param(params, 'initial_values', [1, 0, 1])
+        self.limits = reader.read_param(params, 'limits', [[None, None], [None, None], [None, None]])
+        self.fixes = reader.read_param(params, 'fixed', [False, False, False])
+        self.bins = reader.read_param(params, 'bins', np.linspace(-2, 2, 100))
+        self.print_level = reader.read_param(params, 'print_level', 1)
+
+        self.constrained_parameters = reader.read_param(params, 'constrained_parameter_indices', [])
+        self.constrained_means = reader.read_param(params, 'constrained_parameter_means', [])
+        self.constrained_widths = reader.read_param(params, 'constrained_parameter_widths', [])
+
+        # derived configurations
+        self.bin_centers = self.bins[0:-1]+0.5*(self.bins[1]-self.bins[0])
+        self.parameter_errors = [max([0.1, 0.1*p]) for p in self.initial_values]
+
+        return True
+
+    def InternalRun(self):
+        logger.info('namedata: {}'.format(self.namedata))
+
+        self.hist, _ = np.histogram(self.data[self.namedata], self.bins)
+
+
+        # Now minimize neg log likelihood using iMinuit
+        logger.info('This is the plan:')
+        logger.info('Fitting data consisting of {} elements'.format(np.sum(self.hist)))
+        logger.info('Fit parameters: {}'.format(self.parameter_names))
+        logger.info('Initial values: {}'.format(self.initial_values))
+        logger.info('Initial error: {}'.format(self.parameter_errors))
+        logger.info('Fixed in fit: {}'.format(self.fixes))
+        logger.info('Print level: {}'.format(self.print_level))
+        logger.info('Constrained parameters: {}'.format([self.parameter_names[i] for i in self.constrained_parameters]))
+
+        m_binned = Minuit.from_array_func(self.negPoissonLogLikelihood,
+                                      self.initial_values,
+                                      error=self.parameter_errors,
+                                      errordef = 0.5, limit = self.limits,
+                                      name=self.parameter_names,
+                                      fix=self.fixes, print_level=self.print_level)
+
+
+        # minimze
+        m_binned.migrad(resume=False)
+        result_array = m_binned.np_values()
+        error_array = m_binned.np_errors()
+
+        # save results
+        self.results = {}
+        self.results['param_values'] = result_array
+        self.results['param_errors'] = error_array
+        for i, k in enumerate(self.parameter_names):
+            self.results[k] = {'value': result_array[i], 'error': error_array[i]}
+
+
+        return True
+
+    def PDF(self, x, A, mu, sigma):
+        """
+        Overwrite by whatever PDF
+        """
+        f = A*(1/(sigma*np.sqrt(2*np.pi)))*np.exp(-(((x-mu)/sigma)**2.)/2.)
+        return f
+
+
+    def PoissonLogLikelihood(self, params):
+
+        # binned data
+        hist = self.hist
+
+        # expectation
+        pdf_return = self.PDF(self.bin_centers, *params)
+        if np.shape(pdf_return)[0] == 2:
+            expectation, expectation_error = pdf_return
+        else:
+            expectation = pdf_return
+
+        # exclude bins where expectation is <= zero or nan
+        index = np.where(expectation>0)
+
+        # poisson log likelihoood
+        ll = (hist[index]*np.log(expectation[index]) - expectation[index]).sum()
+
+        # extended ll: poisson total number of events
+        N = np.nansum(expectation)
+        extended_ll = -N+np.sum(hist)*np.log(N)+ll
+        return extended_ll
+
+
+    def negPoissonLogLikelihood(self, params):
+
+        # neg log likelihood
+        neg_ll = - self.PoissonLogLikelihood(params)
+
+        # constrained parameters
+        if len(self.constrained_parameters) > 0:
+            for i in range(len(self.constrained_parameters)):
+                i_param = self.constrained_parameters[i]
+                neg_ll += 0.5 * ((params[i_param] - self.constrained_means[i])/ self.constrained_widths[i])**2
+
+        return neg_ll
diff --git a/mermithid/processors/misc/__init__.py b/mermithid/processors/misc/__init__.py
index 3c9983d2..1a407864 100644
--- a/mermithid/processors/misc/__init__.py
+++ b/mermithid/processors/misc/__init__.py
@@ -5,3 +5,4 @@
 
 from .FrequencyEnergyConversionProcessor import FrequencyEnergyConversionProcessor
 from .FrequencyShifter import FrequencyShifter
+from .BinnedDataFitter import BinnedDataFitter
diff --git a/tests/Misc_test.py b/tests/Misc_test.py
index 54e0edcb..8f26c9bf 100644
--- a/tests/Misc_test.py
+++ b/tests/Misc_test.py
@@ -8,27 +8,99 @@
 
 import unittest
 
-from morpho.utilities import morphologging, read_param
+from morpho.utilities import morphologging, parser
 logger = morphologging.getLogger(__name__)
 
+import matplotlib.pyplot as plt
+import numpy as np
 
-class TritiumTests(unittest.TestCase):
-    from mermithid.processors.misc import FrequencyEnergyConversionProcessor
-    
-    freq_data = [27.9925*10**9, 27.0094*10**9,
-                 26.4195*10**9, 26.4169*10**9,
-                 26.3460*10**9, 26.3457*10**9]
-    logger.info("Will convert the following frequencies: %s"%freq_data)
-    logger.debug("At 1 T, these correspond to kinetic energies (in keV) of " +
-          "[0, 18.6, 30.424, 30.477, 31.934, 31.942]")
-                 
-    freq_energy_dict = {
-        "B": 1
-    }
-
-    freq_proc = FrequencyEnergyConversionProcessor("freq_energy_processor")
-    freq_proc.Configure(freq_energy_dict)
-    freq_proc.frequencies = freq_data
-    freq_proc.Run()
-
-    logger.info("Resulting energies: %s"%freq_proc.energies)
+class MiscTest(unittest.TestCase):
+
+    def test_FreqConversionTest(self):
+        from mermithid.processors.misc import FrequencyEnergyConversionProcessor
+
+        freq_data = [27.9925*10**9, 27.0094*10**9,
+                     26.4195*10**9, 26.4169*10**9,
+                     26.3460*10**9, 26.3457*10**9]
+        logger.info("Will convert the following frequencies: %s"%freq_data)
+        logger.debug("At 1 T, these correspond to kinetic energies (in keV) of " +
+              "[0, 18.6, 30.424, 30.477, 31.934, 31.942]")
+
+        freq_energy_dict = {
+            "B": 1
+        }
+
+        freq_proc = FrequencyEnergyConversionProcessor("freq_energy_processor")
+        freq_proc.Configure(freq_energy_dict)
+        freq_proc.frequencies = freq_data
+        freq_proc.Run()
+
+        logger.info("Resulting energies: %s"%freq_proc.energies)
+
+
+    def test_BinnedDataFitter(self):
+        from mermithid.processors.misc import BinnedDataFitter
+
+        logger.info('iMinuit fit test')
+        config_dict = {
+            'variables': 'K',
+            'bins': np.linspace(-3, 3, 60),
+            'parameter_names': ['A', 'mu', 'sigma'],
+            'initial_values': [100, 0, 1],
+            'limits': [[0, None], [None, None], [0, None]],
+            'constrained_parameter_indices': [1],
+            'constrained_parameter_means': [0],
+            'constrained_parameter_widths': [0.1]
+
+            }
+
+        random_data = {'K': np.random.randn(1000)*0.5+1}
+        negll_fitter = BinnedDataFitter('iminuit_processor')
+        negll_fitter.Configure(config_dict)
+        negll_fitter.data = random_data
+        negll_fitter.Run()
+
+        results = negll_fitter.results
+
+        for k in results.keys():
+            logger.info('{}: {}'.format(k, results[k]))
+
+        x = negll_fitter.bin_centers
+        hist = negll_fitter.hist
+        result_list = results['param_values']
+        error_list = results['param_errors']
+
+        plt.rcParams.update({'font.size': 20})
+        plt.figure(figsize=(10,10))
+        plt.subplot(211)
+        plt.errorbar(x, hist, np.sqrt(hist), drawstyle='steps-mid', label='Binned data')
+        plt.plot(x, negll_fitter.PDF(x, *result_list), label='fit')
+        plt.xlabel(negll_fitter.namedata)
+        plt.ylabel('Counts')
+
+        plt.subplot(212)
+        plt.scatter(x, (hist-negll_fitter.PDF(x, *result_list))/np.sqrt(hist), label='Pearson residuals')
+        plt.xlabel(negll_fitter.namedata)
+        plt.ylabel('Residuals')
+
+
+        plt.savefig('iminit_fit.png')
+
+
+
+
+if __name__ == '__main__':
+
+    args = parser.parse_args(False)
+
+
+    logger = morphologging.getLogger('morpho',
+                                     level=args.verbosity,
+                                     stderr_lb=args.stderr_verbosity,
+                                     propagate=False)
+    logger = morphologging.getLogger(__name__,
+                                     level=args.verbosity,
+                                     stderr_lb=args.stderr_verbosity,
+                                     propagate=False)
+
+    unittest.main()
\ No newline at end of file

From b5b2db127066f60370dc62a0e4b5cf794a783266 Mon Sep 17 00:00:00 2001
From: Talia Weiss <>
Date: Sun, 14 Jun 2020 16:40:53 -0500
Subject: [PATCH 11/49] Fixed wrong var name for survival prob

---
 tests/Fake_data_generator_test.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/Fake_data_generator_test.py b/tests/Fake_data_generator_test.py
index 2f0d0044..fc809250 100644
--- a/tests/Fake_data_generator_test.py
+++ b/tests/Fake_data_generator_test.py
@@ -25,7 +25,7 @@ def test_data_generation(self):
             "use_lineshape": False, # if False only gaussian smearing is applied
             "return_frequency": True,
             "scattering_sigma": 18.6, # only used if use_lineshape = True
-            "scattering_prob": 0.77, # only used if use_lineshape = True
+            "survival_prob": 0.77, # only used if use_lineshape = True
             "scatter_proportion": 0.8, # only used if use_lineshape = True and lineshape = detailed
             "B_field": 0.9578186017836624,
             "S": 4500, # number of tritium events

From ea6af4ef31ee78459aa880aff524675b053be01d Mon Sep 17 00:00:00 2001
From: Xueying-Huyan <huyanxy@gmail.com>
Date: Sun, 14 Jun 2020 16:32:11 -0700
Subject: [PATCH 12/49] added Ar scatter

---
 mermithid/misc/ComplexLineShapeUtilities.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/mermithid/misc/ComplexLineShapeUtilities.py b/mermithid/misc/ComplexLineShapeUtilities.py
index 9e7d48f0..4fe05dc7 100644
--- a/mermithid/misc/ComplexLineShapeUtilities.py
+++ b/mermithid/misc/ComplexLineShapeUtilities.py
@@ -60,6 +60,8 @@ def aseev_func_tail(energy_loss_array, gas_type):
         A2, omeg2, eps2 = 0.4019, 22.31, 16.725
     elif gas_type=="He":
         A2, omeg2, eps2 = 0.1187, 33.40, 10.43
+    elif gas_type=="Ar":
+        A2, omeg2, eps2 = 0.3344, 21.91, 21.14
     return A2*omeg2**2./(omeg2**2.+4*(energy_loss_array-eps2)**2.)
 
 #convert oscillator strength into energy loss spectrum

From 9002d63f3032fcb52f6768bb5a2009eb32d3e690 Mon Sep 17 00:00:00 2001
From: cclaessens <claessens.christine@hotmail.de>
Date: Wed, 17 Jun 2020 15:59:08 +0200
Subject: [PATCH 13/49] info only printed when print level is 1

---
 mermithid/processors/misc/BinnedDataFitter.py | 61 +++++++++++--------
 tests/Misc_test.py                            | 24 +++++---
 2 files changed, 52 insertions(+), 33 deletions(-)

diff --git a/mermithid/processors/misc/BinnedDataFitter.py b/mermithid/processors/misc/BinnedDataFitter.py
index 7d5f01ce..2639f3fd 100644
--- a/mermithid/processors/misc/BinnedDataFitter.py
+++ b/mermithid/processors/misc/BinnedDataFitter.py
@@ -72,29 +72,7 @@ def InternalRun(self):
 
         self.hist, _ = np.histogram(self.data[self.namedata], self.bins)
 
-
-        # Now minimize neg log likelihood using iMinuit
-        logger.info('This is the plan:')
-        logger.info('Fitting data consisting of {} elements'.format(np.sum(self.hist)))
-        logger.info('Fit parameters: {}'.format(self.parameter_names))
-        logger.info('Initial values: {}'.format(self.initial_values))
-        logger.info('Initial error: {}'.format(self.parameter_errors))
-        logger.info('Fixed in fit: {}'.format(self.fixes))
-        logger.info('Print level: {}'.format(self.print_level))
-        logger.info('Constrained parameters: {}'.format([self.parameter_names[i] for i in self.constrained_parameters]))
-
-        m_binned = Minuit.from_array_func(self.negPoissonLogLikelihood,
-                                      self.initial_values,
-                                      error=self.parameter_errors,
-                                      errordef = 0.5, limit = self.limits,
-                                      name=self.parameter_names,
-                                      fix=self.fixes, print_level=self.print_level)
-
-
-        # minimze
-        m_binned.migrad(resume=False)
-        result_array = m_binned.np_values()
-        error_array = m_binned.np_errors()
+        result_array, error_array = self.fit()
 
         # save results
         self.results = {}
@@ -103,9 +81,10 @@ def InternalRun(self):
         for i, k in enumerate(self.parameter_names):
             self.results[k] = {'value': result_array[i], 'error': error_array[i]}
 
-
         return True
 
+
+
     def PDF(self, x, A, mu, sigma):
         """
         Overwrite by whatever PDF
@@ -114,6 +93,40 @@ def PDF(self, x, A, mu, sigma):
         return f
 
 
+
+    def fit(self):
+        # Now minimize neg log likelihood using iMinuit
+        if self.print_level == 1:
+            logger.info('This is the plan:')
+            logger.info('Fitting data consisting of {} elements'.format(np.sum(self.hist)))
+            logger.info('Fit parameters: {}'.format(self.parameter_names))
+            logger.info('Initial values: {}'.format(self.initial_values))
+            logger.info('Initial error: {}'.format(self.parameter_errors))
+            logger.info('Fixed in fit: {}'.format(self.fixes))
+            logger.info('Constrained parameters: {}'.format([self.parameter_names[i] for i in self.constrained_parameters]))
+
+        m_binned = Minuit.from_array_func(self.negPoissonLogLikelihood,
+                                      self.initial_values,
+                                      error=self.parameter_errors,
+                                      errordef = 0.5, limit = self.limits,
+                                      name=self.parameter_names,
+                                      fix=self.fixes,
+                                      print_level=self.print_level,
+                                      throw_nan=True
+                                      )
+
+
+        # minimze
+        m_binned.migrad(resume=False)
+        self.param_states = m_binned.get_param_states()
+        self.m_binned = m_binned
+
+        # results
+        result_array = m_binned.np_values()
+        error_array = m_binned.np_errors()
+        return result_array, error_array
+
+
     def PoissonLogLikelihood(self, params):
 
         # binned data
diff --git a/tests/Misc_test.py b/tests/Misc_test.py
index 8f26c9bf..a311a0b6 100644
--- a/tests/Misc_test.py
+++ b/tests/Misc_test.py
@@ -44,17 +44,17 @@ def test_BinnedDataFitter(self):
         logger.info('iMinuit fit test')
         config_dict = {
             'variables': 'K',
-            'bins': np.linspace(-3, 3, 60),
+            'bins': np.linspace(-3, 3, 100),
             'parameter_names': ['A', 'mu', 'sigma'],
             'initial_values': [100, 0, 1],
             'limits': [[0, None], [None, None], [0, None]],
-            'constrained_parameter_indices': [1],
-            'constrained_parameter_means': [0],
-            'constrained_parameter_widths': [0.1]
+            'constrained_parameter_indices': [],
+            'constrained_parameter_means': [0.5],
+            'constrained_parameter_widths': [1]
 
             }
 
-        random_data = {'K': np.random.randn(1000)*0.5+1}
+        random_data = {'K': np.random.randn(10000)*0.5+1}
         negll_fitter = BinnedDataFitter('iminuit_processor')
         negll_fitter.Configure(config_dict)
         negll_fitter.data = random_data
@@ -65,23 +65,29 @@ def test_BinnedDataFitter(self):
         for k in results.keys():
             logger.info('{}: {}'.format(k, results[k]))
 
-        x = negll_fitter.bin_centers
-        hist = negll_fitter.hist
+
         result_list = results['param_values']
         error_list = results['param_errors']
+        x = negll_fitter.bin_centers
+        hist = negll_fitter.hist
+        hist_fit = negll_fitter.PDF(x, *result_list)
+        residuals = (hist-hist_fit)/np.sqrt(hist_fit)
 
         plt.rcParams.update({'font.size': 20})
         plt.figure(figsize=(10,10))
         plt.subplot(211)
         plt.errorbar(x, hist, np.sqrt(hist), drawstyle='steps-mid', label='Binned data')
-        plt.plot(x, negll_fitter.PDF(x, *result_list), label='fit')
+        plt.plot(x, negll_fitter.PDF(x, *result_list), label='Fit')
         plt.xlabel(negll_fitter.namedata)
         plt.ylabel('Counts')
+        plt.legend()
 
         plt.subplot(212)
-        plt.scatter(x, (hist-negll_fitter.PDF(x, *result_list))/np.sqrt(hist), label='Pearson residuals')
+        plt.scatter(x, residuals, label='Pearson residuals')
+        plt.axhline(np.nanmean(residuals))
         plt.xlabel(negll_fitter.namedata)
         plt.ylabel('Residuals')
+        plt.legend()
 
 
         plt.savefig('iminit_fit.png')

From 8ee8a2f7409a8ce16425f7115612819033fbd5cb Mon Sep 17 00:00:00 2001
From: Talia Weiss <>
Date: Fri, 19 Jun 2020 15:30:06 -0500
Subject: [PATCH 14/49] Removed lineshape normalization option (was not
 working)

---
 .../TritiumSpectrum/FakeDataGenerator.py      |  2 --
 .../processors/misc/KrComplexLineShape.py     | 34 ++++---------------
 2 files changed, 6 insertions(+), 30 deletions(-)

diff --git a/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py b/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
index 60110453..39bba576 100644
--- a/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
+++ b/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
@@ -141,7 +141,6 @@ def InternalConfigure(self, params):
                 # lineshape params
                 self.SimpParams = [self.scattering_sigma*2*math.sqrt(2*math.log(2)), self.survival_prob]
 
-
                 # Setup and configure lineshape processor
                 complexLineShape_config = {
                     'gases': ["H2","He"],
@@ -154,7 +153,6 @@ def InternalConfigure(self, params):
                     'num_points_in_std_array': 10000,
                     'B_field': self.B_field,
                     'base_shape': 'dirac',
-                    'normalize_lineshape': True,
                     'path_to_osc_strengths_files': self.detailed_scatter_spectra_path
                 }
                 logger.info('Setting up complex lineshape object')
diff --git a/mermithid/processors/misc/KrComplexLineShape.py b/mermithid/processors/misc/KrComplexLineShape.py
index 99872568..b57c91e1 100644
--- a/mermithid/processors/misc/KrComplexLineShape.py
+++ b/mermithid/processors/misc/KrComplexLineShape.py
@@ -54,7 +54,7 @@ def InternalConfigure(self, params):
         self.fix_scatter_proportion = reader.read_param(params, 'fix_scatter_proportion', True)
         if self.fix_scatter_proportion == True:
             self.scatter_proportion = reader.read_param(params, 'gas1_scatter_proportion', 0.8)
-        self.normalize_lineshape = reader.read_param(params, 'normalize_lineshape', False)
+            logger.info('Using an H2 scatter proportion of {} with gases {}'.format(self.gases, self.scatter_proportion))
         # This is an important parameter which determines how finely resolved
         # the scatter calculations are. 10000 seems to produce a stable fit, with minimal slowdown
         self.num_points_in_std_array = reader.read_param(params, 'num_points_in_std_array', 10000)
@@ -272,8 +272,6 @@ def make_spectrum(self, gauss_FWHM_eV, prob_parameter, scatter_proportion, emitt
         current_working_spectrum = self.convolve_gaussian(current_working_spectrum, gauss_FWHM_eV)
         zeroth_order_peak = current_working_spectrum
         current_full_spectrum += current_working_spectrum
-        if self.normalize_lineshape==True:
-            norm = 1
         for n in range(1, max_comprehensive_scatters + 1):
             for r in range(0, n + 1):
                 current_working_spectrum = \
@@ -283,10 +281,7 @@ def make_spectrum(self, gauss_FWHM_eV, prob_parameter, scatter_proportion, emitt
                 self.normalize(signal.convolve(zeroth_order_peak, current_working_spectrum, mode='same'))
                 current_full_spectrum += current_working_spectrum*comb(n, r)\
                 *(prob_parameter*p)**(r)*(prob_parameter*q)**(n-r)
-                if self.normalize_lineshape==True:
-                    norm += comb(n, r)\
-                    *(prob_parameter*p)**(r)*(prob_parameter*q)**(n-r)
-
+                
         for n in range(max_comprehensive_scatters + 1, max_scatters + 1):
             current_working_spectrum = \
             scatter_spectra.item()['{}_{}'.format(gases[0], gases[1])]\
@@ -294,8 +289,6 @@ def make_spectrum(self, gauss_FWHM_eV, prob_parameter, scatter_proportion, emitt
             current_working_spectrum = \
             self.normalize(signal.convolve(zeroth_order_peak, current_working_spectrum, mode='same'))
             current_full_spectrum += current_working_spectrum*(prob_parameter*p)**(n)
-            if self.normalize_lineshape==True:
-                norm += (prob_parameter*p)**(n)
 
             current_working_spectrum = \
             scatter_spectra.item()['{}_{}'.format(gases[0], gases[1])]\
@@ -303,12 +296,8 @@ def make_spectrum(self, gauss_FWHM_eV, prob_parameter, scatter_proportion, emitt
             current_working_spectrum = \
             self.normalize(signal.convolve(zeroth_order_peak, current_working_spectrum, mode='same'))
             current_full_spectrum += current_working_spectrum*(prob_parameter*q)**(n)
-            if self.normalize_lineshape==True:
-                norm += (prob_parameter*q)**(n)
-        if self.normalize_lineshape==True:
-            return current_full_spectrum/norm
-        else:
-            return current_full_spectrum
+
+        return current_full_spectrum
 
     # Produces a spectrum in real energy that can now be evaluated off of the SELA.
     #def spectrum_func(x_keV,FWHM_G_eV,line_pos_keV,scatter_prob,amplitude):
@@ -466,8 +455,6 @@ def make_spectrum_1(self, gauss_FWHM_eV, prob_parameter, emitted_peak='shake'):
         current_working_spectrum = self.convolve_gaussian(current_working_spectrum, gauss_FWHM_eV)
         zeroth_order_peak = current_working_spectrum
         current_full_spectrum += current_working_spectrum
-        if self.normalize_lineshape==True:
-            norm = 1
         for n in range(1, max_comprehensive_scatters + 1):
             for r in range(0, n + 1):
                 current_working_spectrum = \
@@ -477,9 +464,6 @@ def make_spectrum_1(self, gauss_FWHM_eV, prob_parameter, emitted_peak='shake'):
                 self.normalize(signal.convolve(zeroth_order_peak, current_working_spectrum, mode='same'))
                 current_full_spectrum += current_working_spectrum*comb(n, r)\
                 *(prob_parameter*p)**(r)*(prob_parameter*q)**(n-r)
-                if self.normalize_lineshape==True:
-                    norm += comb(n, r)\
-                    *(prob_parameter*p)**(r)*(prob_parameter*q)**(n-r)
         for n in range(max_comprehensive_scatters + 1, max_scatters + 1):
             current_working_spectrum = \
             scatter_spectra.item()['{}_{}'.format(gases[0], gases[1])]\
@@ -487,8 +471,6 @@ def make_spectrum_1(self, gauss_FWHM_eV, prob_parameter, emitted_peak='shake'):
             current_working_spectrum = \
             self.normalize(signal.convolve(zeroth_order_peak, current_working_spectrum, mode='same'))
             current_full_spectrum += current_working_spectrum*(prob_parameter*p)**(n)
-            if self.normalize_lineshape==True:
-                norm += (prob_parameter*p)**(n)
 
             current_working_spectrum = \
             scatter_spectra.item()['{}_{}'.format(gases[0], gases[1])]\
@@ -496,12 +478,8 @@ def make_spectrum_1(self, gauss_FWHM_eV, prob_parameter, emitted_peak='shake'):
             current_working_spectrum = \
             self.normalize(signal.convolve(zeroth_order_peak, current_working_spectrum, mode='same'))
             current_full_spectrum += current_working_spectrum*(prob_parameter*q)**(n)
-            if self.normalize_lineshape==True:
-                norm += (prob_parameter*q)**(n)
-        if self.normalize_lineshape==True:
-            return current_full_spectrum/norm
-        else:
-            return current_full_spectrum
+
+        return current_full_spectrum
 
     def spectrum_func_1(self, x_keV, *p0):
         x_eV = x_keV*1000.

From e89af4a76a051824abcd4cff2dcc51251bb7522e Mon Sep 17 00:00:00 2001
From: Noah Oblath <nsoblath@gmail.com>
Date: Wed, 24 Jun 2020 06:06:16 -0700
Subject: [PATCH 15/49] Install iminuit in the Dockerfile (temporary; will
 eventually be done in p8compute_dependencies)

---
 Dockerfile | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/Dockerfile b/Dockerfile
index 1d5b63fb..dbde3af5 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -18,6 +18,10 @@ RUN mkdir -p $MERMITHID_BUILD_PREFIX &&\
     echo 'export PYTHONPATH=$MERMITHID_BUILD_PREFIX/$(python -m site --user-site | sed "s%$(python -m site --user-base)%%"):$PYTHONPATH' >> setup.sh &&\
     /bin/true
 
+RUN source $COMMON_BUILD_PREFIX/setup.sh &&\
+    pip install iminuit &&\
+    /bin/true
+
 ########################
 FROM mermithid_common as mermithid_done
 

From d427805aa31fdf2a72f3e0d81ad5ba36ad34707d Mon Sep 17 00:00:00 2001
From: cclaessens <claessens.christine@hotmail.de>
Date: Wed, 8 Jul 2020 17:36:15 +0200
Subject: [PATCH 16/49] returning lists

---
 mermithid/processors/IO/MultiChannelCicadaReader.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/mermithid/processors/IO/MultiChannelCicadaReader.py b/mermithid/processors/IO/MultiChannelCicadaReader.py
index e3b3bc3c..69241b17 100644
--- a/mermithid/processors/IO/MultiChannelCicadaReader.py
+++ b/mermithid/processors/IO/MultiChannelCicadaReader.py
@@ -99,10 +99,10 @@ def Reader(self):
                 index = np.where((true_frequencies>=self.transition_freqs[i][0]) &(true_frequencies<self.transition_freqs[i][1]))
                 all_frequencies.extend(list(true_frequencies[index]))
 
-                self.data[self.channel_ids[i]]['TrueStartFrequenciesCut'] = true_frequencies[index]
+                self.data[self.channel_ids[i]]['TrueStartFrequenciesCut'] = list(true_frequencies[index])
                 self.data[self.channel_ids[i]]['FrequencyRange'] = [self.transition_freqs[i][0], self.transition_freqs[i][1]]
 
-            self.data[self.frequency_variable_name] = all_frequencies
+            self.data[self.frequency_variable_name] = list(all_frequencies)
         return True
 
     def Writer(self):

From 24c54be85aafe6067e04ae04be3237e23097ad6f Mon Sep 17 00:00:00 2001
From: cclaessens <claessens.christine@hotmail.de>
Date: Sat, 11 Jul 2020 14:20:46 +0200
Subject: [PATCH 17/49] added results logging for print_level 1

---
 mermithid/processors/misc/BinnedDataFitter.py | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/mermithid/processors/misc/BinnedDataFitter.py b/mermithid/processors/misc/BinnedDataFitter.py
index 2639f3fd..edd0b0e2 100644
--- a/mermithid/processors/misc/BinnedDataFitter.py
+++ b/mermithid/processors/misc/BinnedDataFitter.py
@@ -124,6 +124,10 @@ def fit(self):
         # results
         result_array = m_binned.np_values()
         error_array = m_binned.np_errors()
+
+        if self.print_level == 1:
+            logger.info('Fit results: {}'.format(result_array))
+            logger.info('Errors: {}'.format(error_array))
         return result_array, error_array
 
 

From f13d44bed587910f24977187d584b10656516ae2 Mon Sep 17 00:00:00 2001
From: cclaessens <claessens.christine@hotmail.de>
Date: Sat, 11 Jul 2020 14:24:41 +0200
Subject: [PATCH 18/49] option to tilt efficiency curve

---
 mermithid/processors/TritiumSpectrum/FakeDataGenerator.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py b/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
index 60110453..638b05b3 100644
--- a/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
+++ b/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
@@ -90,6 +90,7 @@ def InternalConfigure(self, params):
         self.B =self.A_b*self.runtime*(self.Kmax-self.Kmin) #Background poisson rate
         self.poisson_stats = reader.read_param(params, 'poisson_stats', True)
         self.err_from_B = reader.read_param(params, 'err_from_B', 0.) #In eV, kinetic energy error from f_c --> K conversion
+        self.efficiency_tilt = reader.read_param(params, 'efficiency_tilt', 0.) # per keV
 
 
         #Scattering model parameters
@@ -280,6 +281,11 @@ def generate_unbinned_data(self, Q_mean, mass, ROIbound, S, B_1kev, nsteps=10**4
         if efficiency_dict is not None:
             logger.info('Evaluating efficiencies')
             efficiency_mean, efficiency_error = efficiency_from_interpolation(self.Koptions, efficiency_dict, B_field)
+            slope = self.efficiency_tilt/(1e3)
+            logger.info('Tilted efficiencies: {}'.format(self.efficiency_tilt))
+            efficiency_mean *= (1+slope*(self.Koptions-17.83e3))
+            efficiency_error *= (1+slope*(self.Koptions-17.83e3))
+
             logger.info("Sampling efficiencies given means and uncertainties")
             efficiency = np.random.normal(efficiency_mean, efficiency_error)
             eff_negative = (efficiency<0.)

From 7b8cebe21c7946edcabe6e5b16d9d32a68947285 Mon Sep 17 00:00:00 2001
From: cclaessens <claessens.christine@hotmail.de>
Date: Tue, 14 Jul 2020 21:10:19 +0200
Subject: [PATCH 19/49] added return True to remove configure error

---
 mermithid/processors/misc/KrComplexLineShape.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/mermithid/processors/misc/KrComplexLineShape.py b/mermithid/processors/misc/KrComplexLineShape.py
index b57c91e1..b273ad90 100644
--- a/mermithid/processors/misc/KrComplexLineShape.py
+++ b/mermithid/processors/misc/KrComplexLineShape.py
@@ -69,6 +69,8 @@ def InternalConfigure(self, params):
         if not os.path.exists(self.path_to_osc_strengths_files):
             raise IOError('Path to osc strengths files does not exist')
 
+        return True
+
     def InternalRun(self):
 
         # Read shake parameters from JSON file
@@ -281,7 +283,7 @@ def make_spectrum(self, gauss_FWHM_eV, prob_parameter, scatter_proportion, emitt
                 self.normalize(signal.convolve(zeroth_order_peak, current_working_spectrum, mode='same'))
                 current_full_spectrum += current_working_spectrum*comb(n, r)\
                 *(prob_parameter*p)**(r)*(prob_parameter*q)**(n-r)
-                
+
         for n in range(max_comprehensive_scatters + 1, max_scatters + 1):
             current_working_spectrum = \
             scatter_spectra.item()['{}_{}'.format(gases[0], gases[1])]\

From 40e6060b916383b27190cd590e28ec958268e3b5 Mon Sep 17 00:00:00 2001
From: Talia Weiss <>
Date: Fri, 17 Jul 2020 14:02:07 -0400
Subject: [PATCH 20/49] Added option to set maxf for testing

---
 .../TritiumSpectrum/FakeDataGenerator.py        | 17 ++++++++++++-----
 1 file changed, 12 insertions(+), 5 deletions(-)

diff --git a/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py b/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
index 39bba576..ea5979a0 100644
--- a/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
+++ b/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
@@ -67,7 +67,8 @@ def InternalConfigure(self, params):
         self.m = reader.read_param(params, 'neutrino_mass', 0.2) #Neutrino mass (eV)
         self.Kmin = reader.read_param(params, 'Kmin', self.Q-self.m-2300)  #Energy corresponding to lower bound of frequency ROI (eV)
         self.Kmax = reader.read_param(params, 'Kmax', self.Q-self.m+1000)   #Same, for upper bound (eV)
-        self.minf = reader.read_param(params, 'minf', 25.8e+9) #Minimum frequency
+        self.minf = reader.read_param(params, 'minf', 25813125000.0) #Minimum frequency
+        self.maxf = reader.read_param(params, 'maxf', None)
         if self.Kmax <= self.Kmin:
             logger.error("Kmax <= Kmin!")
             return False
@@ -175,7 +176,10 @@ def InternalConfigure(self, params):
     def InternalRun(self):
 
         if self.return_frequency:
-            ROIbound = [self.minf]
+            if self.maxf == None:
+                ROIbound = [self.minf]
+            else:
+                ROIbound = [self.minf, self.maxf]
         else:
             ROIbound = [self.Kmin, self.Kmax]
 
@@ -246,10 +250,13 @@ def generate_unbinned_data(self, Q_mean, mass, ROIbound, S, B_1kev, nsteps=10**4
 
         if self.return_frequency:
             minf = ROIbound[0]
-            if efficiency_dict is not None:
-                maxf = max(efficiency_dict['frequencies'])
+            if len(ROIbound)==2:
+                maxf = ROIbound[1]
             else:
-                maxf = max(self.load_efficiency_curve()['frequencies'])
+                if efficiency_dict is not None:
+                    maxf = max(efficiency_dict['frequencies'])
+                else:
+                    maxf = max(self.load_efficiency_curve()['frequencies'])
             Kmax, Kmin = Energy(minf, B_field), Energy(maxf, B_field)
         else:
             Kmin, Kmax = ROIbound[0], ROIbound[1]

From a2a91aa6efa77b5fb81982614e59b432fcd26405 Mon Sep 17 00:00:00 2001
From: Talia Weiss <>
Date: Fri, 17 Jul 2020 14:43:17 -0400
Subject: [PATCH 21/49] Cosmetic changes; added to description of
 num_points_in_std_array

---
 mermithid/misc/FakeTritiumDataFunctions.py                | 6 ++----
 mermithid/processors/TritiumSpectrum/FakeDataGenerator.py | 5 +++--
 2 files changed, 5 insertions(+), 6 deletions(-)

diff --git a/mermithid/misc/FakeTritiumDataFunctions.py b/mermithid/misc/FakeTritiumDataFunctions.py
index 343ee319..68308718 100644
--- a/mermithid/misc/FakeTritiumDataFunctions.py
+++ b/mermithid/misc/FakeTritiumDataFunctions.py
@@ -18,8 +18,6 @@
 
 from mermithid.misc.Constants import *
 from mermithid.misc.ConversionFunctions import *
-#from mermithid.misc.KrLineshapeFunctions import *
-
 
 """
 Constants and functions used by processors/TritiumSpectrum/FakeDataGenerator.py
@@ -170,12 +168,12 @@ def bkgd_rate():
     return 1.
 
 
-##Lineshape option: In this case, simply a Gaussian
+#Lineshape option: In this case, simply a Gaussian
 def gaussian(x,a):
     return 1/((2.*np.pi)**0.5*a[0])*(np.exp(-0.5*((x-a[1])/a[0])**2))
 
 
-#Normalized simplified linesahape with scattering
+#Normalized simplified lineshape with scattering
 def simplified_ls(K, Kcenter, FWHM, prob, p0, p1, p2, p3):
     sig0 = FWHM/float(2*np.sqrt(2*np.log(2)))
     shape = gaussian(K, [sig0, Kcenter])
diff --git a/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py b/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
index ea5979a0..c6c122f3 100644
--- a/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
+++ b/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
@@ -110,7 +110,7 @@ def InternalConfigure(self, params):
         self.apply_efficiency = reader.read_param(params, 'apply_efficiency', False)
         self.return_frequency = reader.read_param(params, 'return_frequency', True)
 
-        # nwill be replaced with complex lineshape object if detailed lineshape is used
+        # will be replaced with complex lineshape object if detailed lineshape is used
         self.complexLineShape = None
 
         # get file content if needed
@@ -150,7 +150,8 @@ def InternalConfigure(self, params):
                     # When fix_scatter_proportion is True, set the scatter proportion for gas1 below
                     'gas1_scatter_proportion': self.scatter_proportion,
                     # This is an important parameter which determines how finely resolved
-                    # the scatter calculations are. 10000 seems to produce a stable fit, with minimal slowdown
+                    # the scatter calculations are. 10000 seems to produce a stable fit with minimal slowdown, for ~4000 fake events. The parameter may need to
+                    # be increased for larger datasets.
                     'num_points_in_std_array': 10000,
                     'B_field': self.B_field,
                     'base_shape': 'dirac',

From afa7fa91a7c6f8e79862d591df475204b69f741f Mon Sep 17 00:00:00 2001
From: cclaessens <claessens.christine@hotmail.de>
Date: Sun, 16 Aug 2020 10:18:36 +0200
Subject: [PATCH 22/49] added error message for expectation < 0

also some small changes adressing other review comments
---
 mermithid/processors/misc/BinnedDataFitter.py | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/mermithid/processors/misc/BinnedDataFitter.py b/mermithid/processors/misc/BinnedDataFitter.py
index edd0b0e2..e5c1a433 100644
--- a/mermithid/processors/misc/BinnedDataFitter.py
+++ b/mermithid/processors/misc/BinnedDataFitter.py
@@ -118,7 +118,7 @@ def fit(self):
 
         # minimze
         m_binned.migrad(resume=False)
-        self.param_states = m_binned.get_param_states()
+        self.param_states = m_binned.params
         self.m_binned = m_binned
 
         # results
@@ -133,8 +133,6 @@ def fit(self):
 
     def PoissonLogLikelihood(self, params):
 
-        # binned data
-        hist = self.hist
 
         # expectation
         pdf_return = self.PDF(self.bin_centers, *params)
@@ -143,15 +141,18 @@ def PoissonLogLikelihood(self, params):
         else:
             expectation = pdf_return
 
+        if np.min(expectation) < 0:
+            logger.error('Expectation contains negative numbers. They will be excluded but something could be horribly wrong.')
+
         # exclude bins where expectation is <= zero or nan
         index = np.where(expectation>0)
 
         # poisson log likelihoood
-        ll = (hist[index]*np.log(expectation[index]) - expectation[index]).sum()
+        ll = (self.hist[index]*np.log(expectation[index]) - expectation[index]).sum()
 
         # extended ll: poisson total number of events
         N = np.nansum(expectation)
-        extended_ll = -N+np.sum(hist)*np.log(N)+ll
+        extended_ll = -N+np.sum(self.hist)*np.log(N)+ll
         return extended_ll
 
 
@@ -162,8 +163,7 @@ def negPoissonLogLikelihood(self, params):
 
         # constrained parameters
         if len(self.constrained_parameters) > 0:
-            for i in range(len(self.constrained_parameters)):
-                i_param = self.constrained_parameters[i]
-                neg_ll += 0.5 * ((params[i_param] - self.constrained_means[i])/ self.constrained_widths[i])**2
+            for i, param in enumerate(self.constrained_parameters):
+                neg_ll += 0.5 * ((params[param] - self.constrained_means[i])/ self.constrained_widths[i])**2
 
         return neg_ll

From d1b654cf6877629e408106f4c0c6d52ae26973db Mon Sep 17 00:00:00 2001
From: cclaessens <claessens.christine@hotmail.de>
Date: Sun, 16 Aug 2020 10:39:48 +0200
Subject: [PATCH 23/49] moved BinnedDataFitter to Fitters

---
 .../{misc => Fitters}/BinnedDataFitter.py     |  0
 mermithid/processors/Fitters/__init__.py      |  7 ++
 mermithid/processors/__init__.py              |  1 +
 mermithid/processors/misc/__init__.py         |  1 -
 tests/Fitters_test.py                         | 90 +++++++++++++++++++
 tests/Misc_test.py                            | 55 ------------
 6 files changed, 98 insertions(+), 56 deletions(-)
 rename mermithid/processors/{misc => Fitters}/BinnedDataFitter.py (100%)
 create mode 100644 mermithid/processors/Fitters/__init__.py
 create mode 100644 tests/Fitters_test.py

diff --git a/mermithid/processors/misc/BinnedDataFitter.py b/mermithid/processors/Fitters/BinnedDataFitter.py
similarity index 100%
rename from mermithid/processors/misc/BinnedDataFitter.py
rename to mermithid/processors/Fitters/BinnedDataFitter.py
diff --git a/mermithid/processors/Fitters/__init__.py b/mermithid/processors/Fitters/__init__.py
new file mode 100644
index 00000000..634545ab
--- /dev/null
+++ b/mermithid/processors/Fitters/__init__.py
@@ -0,0 +1,7 @@
+'''
+'''
+
+from __future__ import absolute_import
+
+from .BinnedDataFitter import BinnedDataFitter
+
diff --git a/mermithid/processors/__init__.py b/mermithid/processors/__init__.py
index 430c7bf3..bc325afd 100644
--- a/mermithid/processors/__init__.py
+++ b/mermithid/processors/__init__.py
@@ -7,3 +7,4 @@
 from . import TritiumSpectrum
 from . import plots
 from . import misc
+from . import Fitters
diff --git a/mermithid/processors/misc/__init__.py b/mermithid/processors/misc/__init__.py
index 1a407864..3c9983d2 100644
--- a/mermithid/processors/misc/__init__.py
+++ b/mermithid/processors/misc/__init__.py
@@ -5,4 +5,3 @@
 
 from .FrequencyEnergyConversionProcessor import FrequencyEnergyConversionProcessor
 from .FrequencyShifter import FrequencyShifter
-from .BinnedDataFitter import BinnedDataFitter
diff --git a/tests/Fitters_test.py b/tests/Fitters_test.py
new file mode 100644
index 00000000..6fa2849f
--- /dev/null
+++ b/tests/Fitters_test.py
@@ -0,0 +1,90 @@
+"""
+Script to test fit processors
+Author: C. Claessens
+Date: August 16, 2020
+
+"""
+
+import unittest
+
+from morpho.utilities import morphologging, parser
+logger = morphologging.getLogger(__name__)
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+class FittersTest(unittest.TestCase):
+
+
+    def test_BinnedDataFitter(self):
+        from mermithid.processors.Fitters import BinnedDataFitter
+
+        logger.info('iMinuit fit test')
+        config_dict = {
+            'variables': 'K',
+            'bins': np.linspace(-3, 3, 100),
+            'parameter_names': ['A', 'mu', 'sigma'],
+            'initial_values': [100, 0, 1],
+            'limits': [[0, None], [None, None], [0, None]],
+            'constrained_parameter_indices': [],
+            'constrained_parameter_means': [0.5],
+            'constrained_parameter_widths': [1]
+
+            }
+
+        random_data = {'K': np.random.randn(10000)*0.5+1}
+        negll_fitter = BinnedDataFitter('iminuit_processor')
+        negll_fitter.Configure(config_dict)
+        negll_fitter.data = random_data
+        negll_fitter.Run()
+
+        results = negll_fitter.results
+
+        for k in results.keys():
+            logger.info('{}: {}'.format(k, results[k]))
+
+
+        result_list = results['param_values']
+        error_list = results['param_errors']
+        x = negll_fitter.bin_centers
+        hist = negll_fitter.hist
+        hist_fit = negll_fitter.PDF(x, *result_list)
+        residuals = (hist-hist_fit)/np.sqrt(hist_fit)
+
+        plt.rcParams.update({'font.size': 20})
+        plt.figure(figsize=(10,10))
+        plt.subplot(211)
+        plt.errorbar(x, hist, np.sqrt(hist), drawstyle='steps-mid', label='Binned data')
+        plt.plot(x, negll_fitter.PDF(x, *result_list), label='Fit')
+        plt.xlabel(negll_fitter.namedata)
+        plt.ylabel('Counts')
+        plt.legend()
+
+        plt.subplot(212)
+        plt.scatter(x, residuals, label='Pearson residuals')
+        plt.axhline(np.nanmean(residuals))
+        plt.xlabel(negll_fitter.namedata)
+        plt.ylabel('Residuals')
+        plt.legend()
+
+
+        plt.savefig('iminit_fit.png')
+
+
+
+
+if __name__ == '__main__':
+
+    args = parser.parse_args(False)
+
+
+    logger = morphologging.getLogger('morpho',
+                                     level=args.verbosity,
+                                     stderr_lb=args.stderr_verbosity,
+                                     propagate=False)
+    logger = morphologging.getLogger(__name__,
+                                     level=args.verbosity,
+                                     stderr_lb=args.stderr_verbosity,
+                                     propagate=False)
+
+    unittest.main()
\ No newline at end of file
diff --git a/tests/Misc_test.py b/tests/Misc_test.py
index a311a0b6..d1f05019 100644
--- a/tests/Misc_test.py
+++ b/tests/Misc_test.py
@@ -38,61 +38,6 @@ def test_FreqConversionTest(self):
         logger.info("Resulting energies: %s"%freq_proc.energies)
 
 
-    def test_BinnedDataFitter(self):
-        from mermithid.processors.misc import BinnedDataFitter
-
-        logger.info('iMinuit fit test')
-        config_dict = {
-            'variables': 'K',
-            'bins': np.linspace(-3, 3, 100),
-            'parameter_names': ['A', 'mu', 'sigma'],
-            'initial_values': [100, 0, 1],
-            'limits': [[0, None], [None, None], [0, None]],
-            'constrained_parameter_indices': [],
-            'constrained_parameter_means': [0.5],
-            'constrained_parameter_widths': [1]
-
-            }
-
-        random_data = {'K': np.random.randn(10000)*0.5+1}
-        negll_fitter = BinnedDataFitter('iminuit_processor')
-        negll_fitter.Configure(config_dict)
-        negll_fitter.data = random_data
-        negll_fitter.Run()
-
-        results = negll_fitter.results
-
-        for k in results.keys():
-            logger.info('{}: {}'.format(k, results[k]))
-
-
-        result_list = results['param_values']
-        error_list = results['param_errors']
-        x = negll_fitter.bin_centers
-        hist = negll_fitter.hist
-        hist_fit = negll_fitter.PDF(x, *result_list)
-        residuals = (hist-hist_fit)/np.sqrt(hist_fit)
-
-        plt.rcParams.update({'font.size': 20})
-        plt.figure(figsize=(10,10))
-        plt.subplot(211)
-        plt.errorbar(x, hist, np.sqrt(hist), drawstyle='steps-mid', label='Binned data')
-        plt.plot(x, negll_fitter.PDF(x, *result_list), label='Fit')
-        plt.xlabel(negll_fitter.namedata)
-        plt.ylabel('Counts')
-        plt.legend()
-
-        plt.subplot(212)
-        plt.scatter(x, residuals, label='Pearson residuals')
-        plt.axhline(np.nanmean(residuals))
-        plt.xlabel(negll_fitter.namedata)
-        plt.ylabel('Residuals')
-        plt.legend()
-
-
-        plt.savefig('iminit_fit.png')
-
-
 
 
 if __name__ == '__main__':

From 227216f54d9e24bca0420bdcaac22f09e8e7264d Mon Sep 17 00:00:00 2001
From: cclaessens <claessens.christine@hotmail.de>
Date: Sun, 16 Aug 2020 10:44:57 +0200
Subject: [PATCH 24/49] reverted changes in fake data generator

---
 mermithid/processors/TritiumSpectrum/FakeDataGenerator.py | 6 ------
 1 file changed, 6 deletions(-)

diff --git a/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py b/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
index 8a20a325..c6c122f3 100644
--- a/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
+++ b/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
@@ -91,7 +91,6 @@ def InternalConfigure(self, params):
         self.B =self.A_b*self.runtime*(self.Kmax-self.Kmin) #Background poisson rate
         self.poisson_stats = reader.read_param(params, 'poisson_stats', True)
         self.err_from_B = reader.read_param(params, 'err_from_B', 0.) #In eV, kinetic energy error from f_c --> K conversion
-        self.efficiency_tilt = reader.read_param(params, 'efficiency_tilt', 0.) # per keV
 
 
         #Scattering model parameters
@@ -287,11 +286,6 @@ def generate_unbinned_data(self, Q_mean, mass, ROIbound, S, B_1kev, nsteps=10**4
         if efficiency_dict is not None:
             logger.info('Evaluating efficiencies')
             efficiency_mean, efficiency_error = efficiency_from_interpolation(self.Koptions, efficiency_dict, B_field)
-            slope = self.efficiency_tilt/(1e3)
-            logger.info('Tilted efficiencies: {}'.format(self.efficiency_tilt))
-            efficiency_mean *= (1+slope*(self.Koptions-17.83e3))
-            efficiency_error *= (1+slope*(self.Koptions-17.83e3))
-
             logger.info("Sampling efficiencies given means and uncertainties")
             efficiency = np.random.normal(efficiency_mean, efficiency_error)
             eff_negative = (efficiency<0.)

From 6c9bb9c4badb84154f65017c6e085ddff1f5de3a Mon Sep 17 00:00:00 2001
From: cclaessens <claessens.christine@hotmail.de>
Date: Sun, 16 Aug 2020 10:52:12 +0200
Subject: [PATCH 25/49] added new test script to test.sh

---
 tests/test.sh | 1 +
 1 file changed, 1 insertion(+)

diff --git a/tests/test.sh b/tests/test.sh
index 85c649e3..8483f480 100644
--- a/tests/test.sh
+++ b/tests/test.sh
@@ -5,3 +5,4 @@ python Misc_test.py
 python Tritium_test.py
 python TritiumAndEfficiencyBinner_test.py
 python Fake_data_generator_test.py
+python Fitters_test.py

From fbe046cf96e3291a205c5663f573beb8eb9d7e92 Mon Sep 17 00:00:00 2001
From: cclaessens <claessens.christine@hotmail.de>
Date: Sun, 16 Aug 2020 13:56:43 +0200
Subject: [PATCH 26/49] added molecular final states in fake data generator

---
 mermithid/misc/FakeTritiumDataFunctions.py    | 75 ++++++++++++++++---
 .../TritiumSpectrum/FakeDataGenerator.py      | 23 +++++-
 tests/Fake_data_generator_test.py             |  4 +-
 3 files changed, 89 insertions(+), 13 deletions(-)

diff --git a/mermithid/misc/FakeTritiumDataFunctions.py b/mermithid/misc/FakeTritiumDataFunctions.py
index 68308718..705748c7 100644
--- a/mermithid/misc/FakeTritiumDataFunctions.py
+++ b/mermithid/misc/FakeTritiumDataFunctions.py
@@ -7,12 +7,15 @@
 from __future__ import absolute_import
 
 import numpy as np
+import json
 
 from scipy.special import gamma
 from scipy import integrate
 from scipy.interpolate import interp1d
 from scipy.signal import convolve
 
+import matplotlib.pyplot as plt
+
 from morpho.utilities import morphologging
 logger = morphologging.getLogger(__name__)
 
@@ -146,13 +149,66 @@ def spectral_rate_in_window(K, Q, mnu, Kmin):
     else:
         return 0.
 
+
+def beta_rates(K, Q, mnu, index):
+    beta_rates = np.zeros(len(K))
+    nu_mass_shape = ((Q - K[index])**2 -mnu**2)**0.5
+    beta_rates[index] = GF**2.*Vud**2*Mnuc2/(2.*np.pi**3)*ephasespace(K[index], Q)*(Q - K[index])*nu_mass_shape
+    return beta_rates
+
+
+
 # Unsmeared eta spectrum without a lower energy bound
-def spectral_rate(K, Q, mnu):
-    if Q-mnu > K > 0:
-        return GF**2.*Vud**2*Mnuc2/(2.*np.pi**3)*ephasespace(K, Q)*(Q - K)*np.sqrt((Q - K)**2 - (mnu)**2)
+def spectral_rate(K, Q, mnu, final_states_file):
+
+    # load final state energies and probabilities
+    if final_states_file is None:
+        final_state_array = [0, 1]
     else:
-            return 0.
-    #np.heaviside(Q-mnu-K, 0.5)*np.heaviside(K-V0, 0.5)
+        with open(final_states_file, 'r') as infile:
+            a = json.load(infile)
+            index = np.where(np.array(a['Probability'])[:-1]>0)
+            final_state_array = [np.array(a['Binding energy'])[index], np.array(a['Probability'])[index]]
+
+
+
+    if isinstance(K, list) or isinstance(K, np.ndarray):
+        N_states = len(final_state_array[0])
+        beta_rates_array = np.zeros([N_states, len(K)])
+        Q_states = Q+final_state_array[0]
+
+
+        #K2d = [K for i in range(N_states)]
+        index = [np.where(K < Q_states[i]-mnu) for i in range(N_states)]
+
+
+        beta_rates_array = [beta_rates(K, Q_states[i], mnu, index[i])*final_state_array[1][i] for i in range(N_states)]
+
+        to_return = np.nansum(beta_rates_array, axis=0)/np.nansum(final_state_array[1])
+
+
+        #comparison = beta_rates(K, Q_states[0], mnu, index[0])
+
+        #plt.figure()
+        #plt.plot(K[comparison>0], 1-to_return[comparison>0]/comparison[comparison>0])
+        #plt.ylabel('Final states spectrum relative difference')
+        #plt.xlabel('Energy')
+        #plt.tight_layout()
+        #plt.yscale('log')
+        #plt.savefig('final_state_included_spectrum.png')
+        return to_return
+
+    else:
+
+        return_value = 0.
+
+        for i, e_binding in enumerate(final_state_array[0]):
+            # binding energies are negative
+            Q_state = Q+e_binding
+            if Q_state-mnu > K > 0:
+                return_value += final_state_array[1][i] *(GF**2.*Vud**2*Mnuc2/(2.*np.pi**3)*ephasespace(K, Q_state)*(Q_state - K)*np.sqrt((Q_state - K)**2 - (mnu)**2))
+
+        return return_value/np.sum(final_state_array[1])
 
 
 #Flat background with lower and upper bounds Kmin and Kmax
@@ -244,7 +300,8 @@ def convolved_bkgd_rate(K, Kmin, Kmax, lineshape, ls_params, min_energy, max_ene
 
 #Convolution of signal and lineshape using scipy.signal.convolve
 def convolved_spectral_rate_arrays(K, Q, mnu, Kmin,
-                                   lineshape, ls_params, min_energy, max_energy, complexLineShape):
+                                   lineshape, ls_params, min_energy, max_energy,
+                                   complexLineShape, final_states_file):
     """K is an array-like object
     """
     logger.info('Using scipy convolve')
@@ -265,9 +322,9 @@ def convolved_spectral_rate_arrays(K, Q, mnu, Kmin,
 
         lineshape_rates = complexLineShape.spectrum_func_1(K_lineshape/1000., ls_params[0], 0, 1, ls_params[1])
 
-    beta_rates = np.zeros(len(K))
-    for i,ke in enumerate(K):
-        beta_rates[i] = spectral_rate(ke, Q, mnu)
+    beta_rates = spectral_rate(K, Q, mnu, final_states_file) #np.zeros(len(K))
+    #for i,ke in enumerate(K):
+    #    beta_rates[i] = spectral_rate(ke, Q, mnu, final_states_file)
 
     #Convolving
     convolved = convolve(beta_rates, lineshape_rates, mode='same')
diff --git a/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py b/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
index c6c122f3..98ce0ba4 100644
--- a/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
+++ b/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
@@ -60,6 +60,7 @@ def InternalConfigure(self, params):
         - detailed_or_simplified_lineshape: If use lineshape, this string determines which lineshape model is used.
         - apply_efficiency (boolean): determines whether tritium spectrum is multiplied by efficiency
         - return_frequency: data is always generated as energies. If this parameter is true a list of frequencies is added to the dictionary
+        - molecular_final_states: if True molecular final states are included in tritium spectrum (using spectrum from Bodine et. al 2015)
         """
 
         # Read other parameters
@@ -103,12 +104,14 @@ def InternalConfigure(self, params):
         self.simplified_scattering_path = reader.read_param(params, 'simplified_scattering_path', '/host/input_data/simplified_scattering_params.txt')
         self.detailed_scatter_spectra_path = reader.read_param(params, 'path_to_detailed_scatter_spectra_dir', '.')
         self.efficiency_path = reader.read_param(params, 'efficiency_path', '')
+        self.final_states_file = reader.read_param(params, 'final_states_file', '')
 
         #options
         self.use_lineshape = reader.read_param(params, 'use_lineshape', True)
         self.detailed_or_simplified_lineshape = reader.read_param(params, 'detailed_or_simplified_lineshape', 'detailed')
         self.apply_efficiency = reader.read_param(params, 'apply_efficiency', False)
         self.return_frequency = reader.read_param(params, 'return_frequency', True)
+        self.molecular_final_states = reader.read_param(params, 'molecular_final_states', False)
 
         # will be replaced with complex lineshape object if detailed lineshape is used
         self.complexLineShape = None
@@ -164,12 +167,22 @@ def InternalConfigure(self, params):
                 logger.info('Checking existence of scatter spectra files')
                 self.complexLineShape.check_existence_of_scatter_file()
             else:
-                raise ValueError("'detailed_or_simplified' is neither 'detailed' nor 'simplified'")
+                logger.error("'detailed_or_simplified' is neither 'detailed' nor 'simplified'")
+                return False
 
         else:
             self.lineshape = 'gaussian'
             self.SimpParams = [self.broadening]
             logger.info('Lineshape is Gaussian')
+
+        # check final states file existence
+        if self.molecular_final_states:
+            if not os.path.exists(self.final_states_file):
+                logger.error('Final states file not found')
+                return False
+        else:
+            self.final_states_file = None
+
         return True
 
 
@@ -298,7 +311,9 @@ def generate_unbinned_data(self, Q_mean, mass, ROIbound, S, B_1kev, nsteps=10**4
         time0 = time.time()
 
         if array_method == True:
-            ratesS = convolved_spectral_rate_arrays(self.Koptions, Q_mean, mass, Kmin, lineshape, params, min_energy, max_energy, self.complexLineShape)
+            ratesS = convolved_spectral_rate_arrays(self.Koptions, Q_mean, mass, Kmin,
+                                                    lineshape, params, min_energy, max_energy,
+                                                    self.complexLineShape, self.final_states_file)
         else:
             ratesS = [convolved_spectral_rate(K, Q_mean, mass, Kmin, lineshape, params, min_energy, max_energy) for K in self.Koptions]
 
@@ -310,7 +325,9 @@ def generate_unbinned_data(self, Q_mean, mass, ROIbound, S, B_1kev, nsteps=10**4
 
         # background
         if array_method == True:
-            ratesB = convolved_bkgd_rate_arrays(self.Koptions, Kmin, Kmax, lineshape, params, min_energy, max_energy, self.complexLineShape)
+            ratesB = convolved_bkgd_rate_arrays(self.Koptions, Kmin, Kmax,
+                                                lineshape, params, min_energy, max_energy,
+                                                self.complexLineShape)
         else:
             ratesB = [convolved_bkgd_rate(K, Kmin, Kmax, lineshape, params, min_energy, max_energy) for K in self.Koptions]
 
diff --git a/tests/Fake_data_generator_test.py b/tests/Fake_data_generator_test.py
index fc809250..9333c7e8 100644
--- a/tests/Fake_data_generator_test.py
+++ b/tests/Fake_data_generator_test.py
@@ -31,7 +31,9 @@ def test_data_generation(self):
             "S": 4500, # number of tritium events
             "n_steps": 1000, # stepsize for pseudo continuous data is: (Kmax_eff-Kmin_eff)/nsteps
             "A_b": 1e-10, # background rate 1/eV/s
-            "poisson_stats": True
+            "poisson_stats": True,
+            "molecular_final_states": True,
+            "final_states_file": "../mermithid/misc/saenz_mfs.json"
         }
 
         specGen = FakeDataGenerator("specGen")

From 81279a9a4385eeeb377ab2ba4cfbd391f97765cb Mon Sep 17 00:00:00 2001
From: cclaessens <claessens.christine@hotmail.de>
Date: Sun, 16 Aug 2020 14:49:01 +0200
Subject: [PATCH 27/49] added final state spectrum

---
 mermithid/misc/saenz_mfs.json | 1 +
 1 file changed, 1 insertion(+)
 create mode 100644 mermithid/misc/saenz_mfs.json

diff --git a/mermithid/misc/saenz_mfs.json b/mermithid/misc/saenz_mfs.json
new file mode 100644
index 00000000..3d3b9d7d
--- /dev/null
+++ b/mermithid/misc/saenz_mfs.json
@@ -0,0 +1 @@
+{"Binding energy": [1.897, 1.844, 1.773, 1.65, 1.546, 1.455, 1.341, 1.232, 1.138, 1.047, 0.96, 0.849, 0.754, 0.647999999999999, 0.538, 0.446, 0.345, 0.24, 0.151999999999999, 0.0629999999999999, -0.0429999999999999, -0.147, -0.247, -0.347, -0.446999999999999, -0.613, -0.865, -1.112, -1.36, -1.61, -1.86, -2.186, -2.68199999999999, -3.23499999999999, -3.75, -16.603, -17.603, -18.799, -19.761, -20.73, -21.701, -22.676, -23.653, -24.632, -25.613, -26.596, -27.581, -28.567, -29.558, -30.593, -31.66, -32.637, -33.595, -34.562, -35.548, -36.566, -37.602, -38.609, -39.601, -40.601, -41.607, -42.614, -43.597, -44.584, -45.586, -46.616, -47.601, -48.565, -49.604, -50.599, -51.594, -52.605, -53.611, -54.629, -55.621, -56.632, -57.621, -58.608, -59.608, -60.604, -61.133, -62.615, -63.607, -64.613, -65.604, -66.595, -67.592, -68.589, -69.5759999999999, -70.5809999999999, -71.589, -72.5459999999999, -73.5489999999999, -74.568, -75.533, -76.6149999999999, -77.5669999999999, -78.607, -79.613, -80.6259999999999, -81.6079999999999, -82.6019999999999, -83.5929999999999, -84.5939999999999, -85.601, -86.601, -87.598, -89.0699999999999, -91.086, -93.0849999999999, -95.0849999999999, -97.084, -99.084, -101.086, -103.087, -105.088, -107.089, -109.089, -111.09, -113.09, -115.09, -117.091, -119.091, -121.091, -123.092, -125.092, -127.092, -129.092, -131.093, -133.093, -135.093, -137.093, -140.065, -144.067, -148.068, -152.069, -156.07, -160.072, -164.073, -168.074, -172.075, -176.076, -180.076, -184.077, -188.078, -192.079, -196.079, -200.08, -204.08, -208.081, -212.082, -216.082, -220.082, -224.083, -228.083, -232.084, -236.084], "Probability": [9.99999999999999e-08, 0.0006900000000000001, 0.00046, 0.00233, 0.00553, 0.00457, 0.02033, 0.01649, 0.03877, 0.038079999999999996, 0.06809, 0.11214, 0.10112, 0.24406, 0.32337000000000005, 0.40864, 0.68745, 0.66279, 0.51412, 0.6556100000000001, 0.54588, 0.37231000000000003, 0.25473, 0.16959, 0.11369, 0.16946999999999998, 0.10094, 0.05732, 0.02806, 0.013160000000000002, 0.00623, 0.0042, 0.0008, 0.00015, 0.0, 0.0, 0.0, 1.1999999999999999e-05, 0.000113, 0.0006560000000000001, 0.002567, 0.007149, 0.014804, 0.023583, 0.029715, 0.030307, 0.025527, 0.018080000000000002, 0.01107, 0.007377000000000001, 0.010637, 0.019095, 0.022178, 0.016434, 0.009037, 0.004989, 0.003978, 0.004124, 0.004152, 0.0039250000000000005, 0.003457, 0.003186, 0.0027010000000000003, 0.0027129999999999997, 0.002481, 0.002412, 0.001907, 0.001938, 0.0017599999999999998, 0.001575, 0.0015409999999999998, 0.001485, 0.001557, 0.001895, 0.002427, 0.003357, 0.004095, 0.004714, 0.005033999999999999, 0.005152, 0.005442000000000001, 0.005859, 0.006617, 0.0070940000000000005, 0.007404, 0.007164, 0.006563, 0.005620000000000001, 0.004691, 0.00368, 0.003049, 0.00221, 0.001928, 0.0017610000000000002, 0.0015300000000000001, 0.001215, 0.0013900000000000002, 0.001216, 0.0014219999999999999, 0.001384, 0.001368, 0.001316, 0.001153, 0.0010760000000000001, 0.000921, 0.0007570000000000001, 0.000696, 0.0006180000000000001, 0.00054, 0.00048300000000000003, 0.00043200000000000004, 0.000388, 0.00035150000000000003, 0.00031800000000000003, 0.000289, 0.000264, 0.00024150000000000002, 0.000222, 0.0002045, 0.000189, 0.00017500000000000003, 0.00016250000000000002, 0.000151, 0.000141, 0.0001315, 0.000123, 0.000115, 0.00010800000000000001, 0.0001015, 9.549999999999999e-05, 8.999999999999999e-05, 8.45e-05, 7.775e-05, 6.95e-05, 6.25e-05, 5.625e-05, 5.1000000000000006e-05, 4.625e-05, 4.225e-05, 3.85e-05, 3.525e-05, 3.25e-05, 3e-05, 2.775e-05, 2.5750000000000002e-05, 2.375e-05, 2.225e-05, 2.075e-05, 1.925e-05, 1.8e-05, 1.7e-05, 1.6000000000000003e-05, 1.5e-05, 1.4e-05, 1.325e-05, 1.1750000000000001e-05, NaN]}
\ No newline at end of file

From a2c146f725ba2224b6fe23d0e23c72ff7bbcf9b4 Mon Sep 17 00:00:00 2001
From: cclaessens <claessens.christine@hotmail.de>
Date: Tue, 18 Aug 2020 10:43:23 +0200
Subject: [PATCH 28/49] final state spectrum shifted so largest binding energy
 is 0.

also final state file is now read in in internal config and final state
array is passed on to spectral_rates function.
---
 mermithid/misc/FakeTritiumDataFunctions.py    | 46 ++++++++-----------
 .../TritiumSpectrum/FakeDataGenerator.py      | 12 +++--
 tests/Fake_data_generator_test.py             |  2 +-
 3 files changed, 27 insertions(+), 33 deletions(-)

diff --git a/mermithid/misc/FakeTritiumDataFunctions.py b/mermithid/misc/FakeTritiumDataFunctions.py
index 705748c7..feb95359 100644
--- a/mermithid/misc/FakeTritiumDataFunctions.py
+++ b/mermithid/misc/FakeTritiumDataFunctions.py
@@ -159,43 +159,35 @@ def beta_rates(K, Q, mnu, index):
 
 
 # Unsmeared eta spectrum without a lower energy bound
-def spectral_rate(K, Q, mnu, final_states_file):
-
-    # load final state energies and probabilities
-    if final_states_file is None:
-        final_state_array = [0, 1]
-    else:
-        with open(final_states_file, 'r') as infile:
-            a = json.load(infile)
-            index = np.where(np.array(a['Probability'])[:-1]>0)
-            final_state_array = [np.array(a['Binding energy'])[index], np.array(a['Probability'])[index]]
-
-
+def spectral_rate(K, Q, mnu, final_state_array):
 
     if isinstance(K, list) or isinstance(K, np.ndarray):
         N_states = len(final_state_array[0])
         beta_rates_array = np.zeros([N_states, len(K)])
-        Q_states = Q+final_state_array[0]
 
+        Q_states = Q+final_state_array[0]-np.max(final_state_array[0])
 
-        #K2d = [K for i in range(N_states)]
-        index = [np.where(K < Q_states[i]-mnu) for i in range(N_states)]
+        # plt.figure()
+        # plt.plot(Q_states, final_state_array[1])
+        # plt.plot(Q+final_state_array[0], final_state_array[1])
+        # plt.savefig('q_states.png')
 
+        index = [np.where(K < Q_states[i]-mnu) for i in range(N_states)]
 
         beta_rates_array = [beta_rates(K, Q_states[i], mnu, index[i])*final_state_array[1][i] for i in range(N_states)]
-
         to_return = np.nansum(beta_rates_array, axis=0)/np.nansum(final_state_array[1])
 
 
-        #comparison = beta_rates(K, Q_states[0], mnu, index[0])
+        # comparison = beta_rates(K, Q, mnu, np.where(K < Q - mnu))
 
-        #plt.figure()
-        #plt.plot(K[comparison>0], 1-to_return[comparison>0]/comparison[comparison>0])
-        #plt.ylabel('Final states spectrum relative difference')
-        #plt.xlabel('Energy')
-        #plt.tight_layout()
-        #plt.yscale('log')
-        #plt.savefig('final_state_included_spectrum.png')
+        # plt.figure()
+        # plt.plot(K[comparison>0], 1-to_return[comparison>0]/comparison[comparison>0], marker='', linestyle='-', markersize=1)
+        # plt.ylabel('Final states spectrum relative difference')
+        # plt.xlabel('Energy')
+        # plt.xlim(Q-100, Q+10)
+        # plt.tight_layout()
+        # #plt.yscale('log')
+        # plt.savefig('final_state_included_spectrum.png')
         return to_return
 
     else:
@@ -301,7 +293,7 @@ def convolved_bkgd_rate(K, Kmin, Kmax, lineshape, ls_params, min_energy, max_ene
 #Convolution of signal and lineshape using scipy.signal.convolve
 def convolved_spectral_rate_arrays(K, Q, mnu, Kmin,
                                    lineshape, ls_params, min_energy, max_energy,
-                                   complexLineShape, final_states_file):
+                                   complexLineShape, final_state_array):
     """K is an array-like object
     """
     logger.info('Using scipy convolve')
@@ -322,9 +314,9 @@ def convolved_spectral_rate_arrays(K, Q, mnu, Kmin,
 
         lineshape_rates = complexLineShape.spectrum_func_1(K_lineshape/1000., ls_params[0], 0, 1, ls_params[1])
 
-    beta_rates = spectral_rate(K, Q, mnu, final_states_file) #np.zeros(len(K))
+    beta_rates = spectral_rate(K, Q, mnu, final_state_array) #np.zeros(len(K))
     #for i,ke in enumerate(K):
-    #    beta_rates[i] = spectral_rate(ke, Q, mnu, final_states_file)
+    #    beta_rates[i] = spectral_rate(ke, Q, mnu, final_state_array)
 
     #Convolving
     convolved = convolve(beta_rates, lineshape_rates, mode='same')
diff --git a/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py b/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
index 98ce0ba4..c43bb1f1 100644
--- a/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
+++ b/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
@@ -177,11 +177,13 @@ def InternalConfigure(self, params):
 
         # check final states file existence
         if self.molecular_final_states:
-            if not os.path.exists(self.final_states_file):
-                logger.error('Final states file not found')
-                return False
+            logger.info('Going to use molecular final states from Bodine et al 2015')
+            with open(self.final_states_file, 'r') as infile:
+                a = json.load(infile)
+                index = np.where(np.array(a['Probability'])[:-1]>0)
+                self.final_state_array = [np.array(a['Binding energy'])[index], np.array(a['Probability'])[index]]
         else:
-            self.final_states_file = None
+            self.final_state_array = [0, 1]
 
         return True
 
@@ -313,7 +315,7 @@ def generate_unbinned_data(self, Q_mean, mass, ROIbound, S, B_1kev, nsteps=10**4
         if array_method == True:
             ratesS = convolved_spectral_rate_arrays(self.Koptions, Q_mean, mass, Kmin,
                                                     lineshape, params, min_energy, max_energy,
-                                                    self.complexLineShape, self.final_states_file)
+                                                    self.complexLineShape, self.final_state_array)
         else:
             ratesS = [convolved_spectral_rate(K, Q_mean, mass, Kmin, lineshape, params, min_energy, max_energy) for K in self.Koptions]
 
diff --git a/tests/Fake_data_generator_test.py b/tests/Fake_data_generator_test.py
index 9333c7e8..0b5207e8 100644
--- a/tests/Fake_data_generator_test.py
+++ b/tests/Fake_data_generator_test.py
@@ -29,7 +29,7 @@ def test_data_generation(self):
             "scatter_proportion": 0.8, # only used if use_lineshape = True and lineshape = detailed
             "B_field": 0.9578186017836624,
             "S": 4500, # number of tritium events
-            "n_steps": 1000, # stepsize for pseudo continuous data is: (Kmax_eff-Kmin_eff)/nsteps
+            "n_steps": 10000, # stepsize for pseudo continuous data is: (Kmax_eff-Kmin_eff)/nsteps
             "A_b": 1e-10, # background rate 1/eV/s
             "poisson_stats": True,
             "molecular_final_states": True,

From 36e8da5c9b72dad3bf4b27226044793a1924b55e Mon Sep 17 00:00:00 2001
From: cclaessens <claessens.christine@hotmail.de>
Date: Tue, 18 Aug 2020 10:45:49 +0200
Subject: [PATCH 29/49] back to get_param_states() because of backwards
 compatibility

---
 mermithid/processors/Fitters/BinnedDataFitter.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/mermithid/processors/Fitters/BinnedDataFitter.py b/mermithid/processors/Fitters/BinnedDataFitter.py
index e5c1a433..19359f43 100644
--- a/mermithid/processors/Fitters/BinnedDataFitter.py
+++ b/mermithid/processors/Fitters/BinnedDataFitter.py
@@ -118,7 +118,7 @@ def fit(self):
 
         # minimze
         m_binned.migrad(resume=False)
-        self.param_states = m_binned.params
+        self.param_states = m_binned.get_param_states()
         self.m_binned = m_binned
 
         # results

From a6d2ec51d67f4b89c7f1a8b3496d8080351ac353 Mon Sep 17 00:00:00 2001
From: cclaessens <claessens.christine@hotmail.de>
Date: Tue, 18 Aug 2020 10:57:43 +0200
Subject: [PATCH 30/49] Changed name of PDF to model because it isn't actually
 a pdf. Extended description.

---
 .../processors/Fitters/BinnedDataFitter.py     | 18 ++++++++++--------
 tests/Fitters_test.py                          |  4 ++--
 2 files changed, 12 insertions(+), 10 deletions(-)

diff --git a/mermithid/processors/Fitters/BinnedDataFitter.py b/mermithid/processors/Fitters/BinnedDataFitter.py
index 19359f43..82ac57bd 100644
--- a/mermithid/processors/Fitters/BinnedDataFitter.py
+++ b/mermithid/processors/Fitters/BinnedDataFitter.py
@@ -1,8 +1,10 @@
 '''
 Author: C. Claessens
 Date:6/12/2020
-Description: Configure with arbitrary histogram and pdf.
-Minimizes poisson loglikelihood using iMinuit.
+Description:
+    Minimizes poisson loglikelihood using iMinuit.
+    Configure with arbitrary histogram and model.
+    The model should consist of a pdf multiplied with a free amplitude parameter representing the total number of events.
 
 '''
 
@@ -85,9 +87,9 @@ def InternalRun(self):
 
 
 
-    def PDF(self, x, A, mu, sigma):
+    def model(self, x, A, mu, sigma):
         """
-        Overwrite by whatever PDF
+        Overwrite by whatever Model
         """
         f = A*(1/(sigma*np.sqrt(2*np.pi)))*np.exp(-(((x-mu)/sigma)**2.)/2.)
         return f
@@ -135,11 +137,11 @@ def PoissonLogLikelihood(self, params):
 
 
         # expectation
-        pdf_return = self.PDF(self.bin_centers, *params)
-        if np.shape(pdf_return)[0] == 2:
-            expectation, expectation_error = pdf_return
+        model_return = self.model(self.bin_centers, *params)
+        if np.shape(model_return)[0] == 2:
+            expectation, expectation_error = model_return
         else:
-            expectation = pdf_return
+            expectation = model_return
 
         if np.min(expectation) < 0:
             logger.error('Expectation contains negative numbers. They will be excluded but something could be horribly wrong.')
diff --git a/tests/Fitters_test.py b/tests/Fitters_test.py
index 6fa2849f..a8c00240 100644
--- a/tests/Fitters_test.py
+++ b/tests/Fitters_test.py
@@ -48,14 +48,14 @@ def test_BinnedDataFitter(self):
         error_list = results['param_errors']
         x = negll_fitter.bin_centers
         hist = negll_fitter.hist
-        hist_fit = negll_fitter.PDF(x, *result_list)
+        hist_fit = negll_fitter.model(x, *result_list)
         residuals = (hist-hist_fit)/np.sqrt(hist_fit)
 
         plt.rcParams.update({'font.size': 20})
         plt.figure(figsize=(10,10))
         plt.subplot(211)
         plt.errorbar(x, hist, np.sqrt(hist), drawstyle='steps-mid', label='Binned data')
-        plt.plot(x, negll_fitter.PDF(x, *result_list), label='Fit')
+        plt.plot(x, negll_fitter.model(x, *result_list), label='Fit')
         plt.xlabel(negll_fitter.namedata)
         plt.ylabel('Counts')
         plt.legend()

From 3d91729e7b6a98182e4b29a004d712a7287478ff Mon Sep 17 00:00:00 2001
From: cclaessens <claessens.christine@hotmail.de>
Date: Wed, 19 Aug 2020 14:34:13 +0200
Subject: [PATCH 31/49] Added option to pass already binned data

---
 .../processors/Fitters/BinnedDataFitter.py    | 34 ++++++++++++-------
 tests/Fitters_test.py                         | 32 ++++++++++++++---
 2 files changed, 49 insertions(+), 17 deletions(-)

diff --git a/mermithid/processors/Fitters/BinnedDataFitter.py b/mermithid/processors/Fitters/BinnedDataFitter.py
index 82ac57bd..dbaccc0f 100644
--- a/mermithid/processors/Fitters/BinnedDataFitter.py
+++ b/mermithid/processors/Fitters/BinnedDataFitter.py
@@ -29,21 +29,22 @@ class BinnedDataFitter(BaseProcessor):
     '''
     Processor that
     Args:
-        variables
-        parameter_names
-        initial_values
-        limits
-        fixed
-        bins
-        print_level
-        constrained_parameter_indices
-        constrained_parameter_means
-        constrained_parameter_widths
+        variables: dictionary key under which data is stored
+        parameter_names: names of model parameters
+        initial_values: list of parameter initial values
+        limits: parameter limits given as [lower, upper] for each parameter
+        fixed: boolean list of same length as parameters. Only un-fixed parameters will be fitted
+        bins: bins for data histogram
+        binned_data: if True data is assumed to already be histogrammed
+        print_level: if 1 fit plan and result summaries are printed
+        constrained_parameter_indices: list of indices indicating which parameters are contraint. Contstraints will be Gaussian.
+        constrained_parameter_means: Mean of Gaussian constraint
+        constrained_parameter_widths: Standard deviation of Gaussian constrained
 
     Inputs:
         data:
     Output:
-        result: dictionary containing
+        result: dictionary containing fit results and uncertainties
     '''
     def InternalConfigure(self, params):
         '''
@@ -57,6 +58,7 @@ def InternalConfigure(self, params):
         self.limits = reader.read_param(params, 'limits', [[None, None], [None, None], [None, None]])
         self.fixes = reader.read_param(params, 'fixed', [False, False, False])
         self.bins = reader.read_param(params, 'bins', np.linspace(-2, 2, 100))
+        self.binned_data = reader.read_param(params, 'binned_data', False)
         self.print_level = reader.read_param(params, 'print_level', 1)
 
         self.constrained_parameters = reader.read_param(params, 'constrained_parameter_indices', [])
@@ -72,7 +74,15 @@ def InternalConfigure(self, params):
     def InternalRun(self):
         logger.info('namedata: {}'.format(self.namedata))
 
-        self.hist, _ = np.histogram(self.data[self.namedata], self.bins)
+        if self.binned_data:
+            logger.info('Data is already binned. Getting binned data')
+            self.hist = self.data[self.namedata]
+            if len(self.hist) != len(self.bin_centers):
+                logger.error('Number of bins and histogram entries do not match')
+                raise ValueError('Number of bins and histogram entries do not match')
+        else:
+            logger.info('Data is unbinned. Will histogram before fitting.')
+            self.hist, _ = np.histogram(self.data[self.namedata], self.bins)
 
         result_array, error_array = self.fit()
 
diff --git a/tests/Fitters_test.py b/tests/Fitters_test.py
index a8c00240..debc2bef 100644
--- a/tests/Fitters_test.py
+++ b/tests/Fitters_test.py
@@ -21,34 +21,56 @@ def test_BinnedDataFitter(self):
 
         logger.info('iMinuit fit test')
         config_dict = {
-            'variables': 'K',
+            'variables': 'N',
             'bins': np.linspace(-3, 3, 100),
             'parameter_names': ['A', 'mu', 'sigma'],
             'initial_values': [100, 0, 1],
             'limits': [[0, None], [None, None], [0, None]],
             'constrained_parameter_indices': [],
             'constrained_parameter_means': [0.5],
-            'constrained_parameter_widths': [1]
+            'constrained_parameter_widths': [1],
+            'binned_data': True
 
             }
 
         random_data = {'K': np.random.randn(10000)*0.5+1}
+
+        # histogramming could be done by processor
+        binned_data, _ = np.histogram(random_data['K'], config_dict['bins'])
+        binned_data_dict = {'N': binned_data}
         negll_fitter = BinnedDataFitter('iminuit_processor')
         negll_fitter.Configure(config_dict)
-        negll_fitter.data = random_data
+        negll_fitter.data = binned_data_dict #random_data
+
+        # before running: overwrite model
+        def gaussian(x, A, mu, sigma):
+            """
+            This is the same function that is implemented as default model
+            """
+            f = A*(1/(sigma*np.sqrt(2*np.pi)))*np.exp(-(((x-mu)/sigma)**2.)/2.)
+            return f
+
+        negll_fitter.model = gaussian
+
+        # now run
         negll_fitter.Run()
 
+        # collect fit results
         results = negll_fitter.results
+        result_list = results['param_values']
+        error_list = results['param_errors']
 
         for k in results.keys():
             logger.info('{}: {}'.format(k, results[k]))
 
 
-        result_list = results['param_values']
-        error_list = results['param_errors']
+
+        # get bins histogram and fit curve from processor for plotting
         x = negll_fitter.bin_centers
         hist = negll_fitter.hist
         hist_fit = negll_fitter.model(x, *result_list)
+
+        # calculate normalized residuals
         residuals = (hist-hist_fit)/np.sqrt(hist_fit)
 
         plt.rcParams.update({'font.size': 20})

From fa42e795c5e6803cfac6fc196e4468e478546638 Mon Sep 17 00:00:00 2001
From: cclaessens <claessens.christine@hotmail.de>
Date: Wed, 19 Aug 2020 14:43:52 +0200
Subject: [PATCH 32/49] minor rearrangements in test script

---
 tests/Fitters_test.py | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/tests/Fitters_test.py b/tests/Fitters_test.py
index debc2bef..c54e78d5 100644
--- a/tests/Fitters_test.py
+++ b/tests/Fitters_test.py
@@ -38,11 +38,12 @@ def test_BinnedDataFitter(self):
         # histogramming could be done by processor
         binned_data, _ = np.histogram(random_data['K'], config_dict['bins'])
         binned_data_dict = {'N': binned_data}
+
+        # setup processor
         negll_fitter = BinnedDataFitter('iminuit_processor')
         negll_fitter.Configure(config_dict)
-        negll_fitter.data = binned_data_dict #random_data
 
-        # before running: overwrite model
+        # define new model and overwrite processor's model
         def gaussian(x, A, mu, sigma):
             """
             This is the same function that is implemented as default model
@@ -52,7 +53,8 @@ def gaussian(x, A, mu, sigma):
 
         negll_fitter.model = gaussian
 
-        # now run
+        # hand over data and run
+        negll_fitter.data = binned_data_dict #random_data
         negll_fitter.Run()
 
         # collect fit results
@@ -90,7 +92,7 @@ def gaussian(x, A, mu, sigma):
         plt.legend()
 
 
-        plt.savefig('iminit_fit.png')
+        plt.savefig('iminuit_fit.png')
 
 
 

From 22f365e552c7ec05610035377d0e621bce304036 Mon Sep 17 00:00:00 2001
From: cclaessens <claessens.christine@hotmail.de>
Date: Tue, 20 Oct 2020 20:25:39 +0200
Subject: [PATCH 33/49] small fix in fake data generator

---
 mermithid/processors/TritiumSpectrum/FakeDataGenerator.py | 3 ++-
 tests/Fake_data_generator_test.py                         | 2 +-
 2 files changed, 3 insertions(+), 2 deletions(-)

diff --git a/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py b/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
index c43bb1f1..9c5e80f9 100644
--- a/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
+++ b/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
@@ -183,7 +183,8 @@ def InternalConfigure(self, params):
                 index = np.where(np.array(a['Probability'])[:-1]>0)
                 self.final_state_array = [np.array(a['Binding energy'])[index], np.array(a['Probability'])[index]]
         else:
-            self.final_state_array = [0, 1]
+            logger.info('Not using molecular final state spectrum')
+            self.final_state_array = np.array([[0.], [1.]])
 
         return True
 
diff --git a/tests/Fake_data_generator_test.py b/tests/Fake_data_generator_test.py
index 0b5207e8..6e729297 100644
--- a/tests/Fake_data_generator_test.py
+++ b/tests/Fake_data_generator_test.py
@@ -32,7 +32,7 @@ def test_data_generation(self):
             "n_steps": 10000, # stepsize for pseudo continuous data is: (Kmax_eff-Kmin_eff)/nsteps
             "A_b": 1e-10, # background rate 1/eV/s
             "poisson_stats": True,
-            "molecular_final_states": True,
+            "molecular_final_states": False,
             "final_states_file": "../mermithid/misc/saenz_mfs.json"
         }
 

From 5e9b212511f079b1d6cf63d3dad8503f4203f7e8 Mon Sep 17 00:00:00 2001
From: cclaessens <claessens.christine@hotmail.de>
Date: Mon, 26 Oct 2020 15:27:24 +0100
Subject: [PATCH 34/49] I believe we are not using this file anymore

---
 mermithid/misc/KrLineshapeFunctions.py | 318 -------------------------
 1 file changed, 318 deletions(-)
 delete mode 100644 mermithid/misc/KrLineshapeFunctions.py

diff --git a/mermithid/misc/KrLineshapeFunctions.py b/mermithid/misc/KrLineshapeFunctions.py
deleted file mode 100644
index 4fa791dd..00000000
--- a/mermithid/misc/KrLineshapeFunctions.py
+++ /dev/null
@@ -1,318 +0,0 @@
-'''
-Various functions for the Krypton lineshape computation
-Adapted from https://github.com/project8/scripts/blob/feature/KrScatterFit/machado/fitting/data_fitting_rebuild.py#L429 - by E. M. Machado
-Changes from that version:
-- Kr specific functions (shake-up/shake-off) removed
-- Spectrum is normalized
-
-Author: T. Weiss, C. Claessens, Y. Sun
-Date:4/6/2020
-'''
-
-from __future__ import absolute_import
-
-import numpy as np
-
-from scipy import integrate
-from scipy.signal import convolve
-import os
-
-from morpho.utilities import morphologging
-logger = morphologging.getLogger(__name__)
-
-from mermithid.misc.Constants import *
-from mermithid.misc.ConversionFunctions import *
-
-
-# Natural constants
-kr_line = kr_k_line_e()*1e-3 #17.8260 # keV
-kr_line_width = kr_k_line_width() #2.83 # eV
-e_charge = e() #1.60217662*10**(-19) # Coulombs , charge of electron
-m_e = m_electron()/e()*c()**2 #9.10938356*10**(-31) # Kilograms , mass of electron
-mass_energy_electron = m_electron()*1e-3 #510.9989461 # keV
-max_scatters = 20
-gases = ["H2"]
-
-# This is an important parameter which determines how finely resolved
-# the scatter calculations are. 10000 seems to produce a stable fit, with minimal slowdown
-num_points_in_std_array = 10000
-
-# Establishes a standard energy loss array (SELA) from -1000 eV to 1000 eV
-# with number of points equal to num_points_in_std_array. All convolutions
-# will be carried out on this particular discretization
-def std_eV_array():
-    emin = -1000
-    emax = 1000
-    array = np.linspace(emin,emax,num_points_in_std_array)
-    return array
-
-# A lorentzian function
-def lorentzian(x_array,x0,FWHM):
-    HWHM = FWHM/2.
-    func = (1./np.pi)*HWHM/((x_array-x0)**2.+HWHM**2.)
-    return func
-
-# A lorentzian line centered at 0 eV, with 2.83 eV width on the SELA
-def std_lorentzian_17keV():
-    x_array = std_eV_array()
-    ans = lorentzian(x_array,0,kr_line_width)
-    return ans
-
-
-# A gaussian function
-def gaussian(x_array,A,sigma=0,mu=0):
-    if isinstance(A, list):
-        a = A
-        x = x_array
-        return 1/((2.*np.pi)**0.5*a[0])*(np.exp(-0.5*((x-a[1])/a[0])**2))
-    else:
-        f = A*(1./(sigma*np.sqrt(2*np.pi)))*np.exp(-(((x_array-mu)/sigma)**2.)/2.)
-        return f
-
-# A gaussian centered at 0 eV with variable width, on the SELA
-def std_gaussian(sigma):
-    x_array = std_eV_array()
-    ans = gaussian(x_array,1,sigma,0)
-    return ans
-
-def std_dirac():
-    x_array = std_eV_array()
-    ans = np.zeros(len(x_array))
-    min_x = np.min(np.abs(x_array))
-    ans[np.abs(x_array)==min_x] = 1.
-    logger.warning('Spectrum will be shifted by lineshape by {} eV'.format(min_x))
-    if min_x > 0.1:
-        logger.warning('Lineshape will shift spectrum by > 0.1 eV')
-    if min_x > 1.:
-        logger.warning('Lineshape will shift spectrum by > 1 eV')
-        raise ValueError('problem with std_eV_array()')
-    return ans
-
-# Converts a gaussian's FWHM to sigma
-def gaussian_FWHM_to_sigma(FWHM):
-    sigma = FWHM/(2*np.sqrt(2*np.log(2)))
-    return sigma
-
-# Converts a gaussian's sigma to FWHM
-def gaussian_sigma_to_FWHM(sigma):
-    FWHM = sigma*(2*np.sqrt(2*np.log(2)))
-    return FWHM
-
-# normalizes a function, but depends on binning.
-# Only to be used for functions evaluated on the SELA
-def normalize(f):
-    x_arr = std_eV_array()
-    f_norm = integrate.simps(f,x=x_arr)
-    # if  f_norm < 0.99 or f_norm > 1.01:
-    #     print(f_norm)
-    f_normed = f/f_norm
-    return f_normed
-
-#returns array with energy loss/ oscillator strength data
-#def read_oscillator_str_file(filename):
-#    f = open(filename, "r")
-#    lines = f.readlines()
-#    energyOsc = [[],[]] #2d array of energy losses, oscillator strengths
-#
-#    for line in lines:
-#        if line != "" and line[0]!="#":
-#            raw_data = [float(i) for i in line.split("\t")]
-#            energyOsc[0].append(raw_data[0])
-#            energyOsc[1].append(raw_data[1])
-#
-#    energyOsc = np.array(energyOsc)
-#    ### take data and sort by energy
-#    sorted_indices = energyOsc[0].argsort()
-#    energyOsc[0] = energyOsc[0][sorted_indices]
-#    energyOsc[1] = energyOsc[1][sorted_indices]
-#    return energyOsc
-
-# A sub function for the scatter function. Found in
-# "Energy loss of 18 keV electrons in gaseous T and quench condensed D films"
-# by V.N. Aseev et al. 2000
-#def aseev_func_tail(energy_loss_array, gas_type):
-#    if gas_type=="H2":
-#        A2, omeg2, eps2 = 0.195, 14.13, 10.60
-#    elif gas_type=="Kr":
-#        A2, omeg2, eps2 = 0.4019, 22.31, 16.725
-#    return A2*omeg2**2./(omeg2**2.+4*(energy_loss_array-eps2)**2.)
-
-#convert oscillator strength into energy loss spectrum
-#def get_eloss_spec(e_loss, oscillator_strength, kinetic_en = kr_line * 1000): #energies in eV
-#    e_rydberg = 13.605693009 #rydberg energy (eV)
-#    a0 = 5.291772e-11 #bohr radius
-#    return np.where(e_loss>0 , 4.*np.pi*a0**2 * e_rydberg / (kinetic_en * e_loss) * oscillator_strength * np.log(4. * kinetic_en * e_loss / (e_rydberg**3.) ), 0)
-
-# Function for energy loss from a single scatter of electrons by
-# V.N. Aseev et al. 2000
-# This function does the work of combining fit_func1 and fit_func2 by
-# finding the point where they intersect.
-# Evaluated on the SELA
-#def single_scatter_f(gas_type):
-#    energy_loss_array = std_eV_array()
-#    f = 0 * energy_loss_array
-#
-#    input_filename = "../data/" + gas_type + "OscillatorStrength.txt"
-#    energy_fOsc = read_oscillator_str_file(input_filename)
-#    fData = interpolate.interp1d(energy_fOsc[0], energy_fOsc[1], kind='linear')
-#    for i in range(len(energy_loss_array)):
-#        if energy_loss_array[i] < energy_fOsc[0][0]:
-#            f[i] = 0
-#        elif energy_loss_array[i] <= energy_fOsc[0][-1]:
-#            f[i] = fData(energy_loss_array[i])
-#        else:
-#            f[i] = aseev_func_tail(energy_loss_array[i], gas_type)
-#
-#    f_e_loss = get_eloss_spec(energy_loss_array, f)
-#    f_normed = normalize(f_e_loss)
-#    #plt.plot(energy_loss_array, f_e_loss)
-#    #plt.show()
-#    return f_normed
-
-# Convolves a function with the single scatter function, on the SELA
-#def another_scatter(input_spectrum, gas_type):
-#    single = single_scatter_f(gas_type)
-#    f = convolve(single,input_spectrum,mode='same')
-#    f_normed = normalize(f)
-#    return f_normed
-
-# Convolves the scatter function with itself over and over again and saves
-# the results to .npy files.
-#def generate_scatter_convolution_files(gas_type):
-#    t = time.time()
-#    first_scatter = single_scatter_f(gas_type)
-#    scatter_num_array = range(2,max_scatters+1)
-#    current_scatter = first_scatter
-#    np.save('scatter_spectra_files/scatter'+gas_type+"_"+str(1).zfill(2),current_scatter)
-#    # x = std_eV_array() # diagnostic
-#    for i in scatter_num_array:
-#        current_scatter = another_scatter(current_scatter, gas_type)
-#        np.save('scatter_spectra_files/scatter'+gas_type+"_"+str(i).zfill(2),current_scatter)
-#    # plt.plot(x,current_scatter) # diagnostic
-#    # plt.show() # diagnostic
-#    elapsed = time.time() - t
-#    logger.info('Files generated in '+str(elapsed)+'s')
-#    return
-
-# Returns the name of the current path
-def get_current_path():
-    path = os.path.abspath(os.getcwd())
-    return path
-
-# Prints a list of the contents of a directory
-def list_files(path):
-    directory = os.popen("ls "+path).readlines()
-    strippeddirs = [s.strip('\n') for s in directory]
-    return strippeddirs
-
-# Returns the name of the current directory
-def get_current_dir():
-    current_path = os.popen("pwd").readlines()
-    stripped_path = [s.strip('\n') for s in current_path]
-    stripped_path = stripped_path[0].split('/')
-    current_dir = stripped_path[len(stripped_path)-1]
-    return current_dir
-
-# Checks for the existence of a directory called 'scatter_spectra_files'
-# and checks that this directory contains the scatter spectra files.
-# If not, this function calls generate_scatter_convolution_files.
-# This function also checks to make sure that the scatter files have the correct
-# number of points in the SELA, and if not, it generates fresh files
-#def check_existence_of_scatter_files(gas_type):
-#    current_path = get_current_path()
-#    current_dir = get_current_dir()
-#    stuff_in_dir = list_files(current_path)
-#    if 'scatter_spectra_files' not in stuff_in_dir and current_dir != 'scatter_spectra_files':
-#        logger.warning('Scatter files not found, generating')
-#        os.popen("mkdir scatter_spectra_files")
-#        time.sleep(2)
-#        generate_scatter_convolution_files(gas_type)
-#    else:
-#        directory = os.popen("ls scatter_spectra_files").readlines()
-#        strippeddirs = [s.strip('\n') for s in directory]
-#        if len(directory) != len(gases) * max_scatters:
-#            generate_scatter_convolution_files(gas_type)
-#        test_file = 'scatter_spectra_files/scatter'+gas_type+'_01.npy'
-#        test_arr = np.load(test_file)
-#        if len(test_arr) != num_points_in_std_array:
-#            logger.warning('Scatter files do not match standard array binning, generating fresh files')
-#            generate_scatter_convolution_files(gas_type)
-#    return
-
-# Given a function evaluated on the SELA, convolves it with a gaussian
-def convolve_gaussian(func_to_convolve,gauss_FWHM_eV):
-    sigma = gaussian_FWHM_to_sigma(gauss_FWHM_eV)
-    resolution_f = std_gaussian(sigma)
-    ans = convolve(resolution_f,func_to_convolve,mode='same')
-    ans_normed = normalize(ans)
-    return ans_normed
-
-# Produces a full spectral shape on the SELA, given a gaussian resolution
-# and a scatter probability
-def make_spectrum(gauss_FWHM_eV,scatter_prob,gas_type, emitted_peak='dirac'):
-    current_path = get_current_path()
-    # check_existence_of_scatter_files()
-    #filenames = list_files('scatter_spectra_files')
-    scatter_num_array = range(1,max_scatters+1)
-    en_array = std_eV_array()
-    current_full_spectrum = np.zeros(len(en_array))
-    if emitted_peak == 'lorentzian':
-        current_working_spectrum = std_lorentzian_17keV() #Normalized
-    elif emitted_peak == 'shake':
-        current_working_spectrum = shake_spectrum()
-    else:
-        current_working_spectrum = std_dirac()
-    current_working_spectrum = convolve_gaussian(current_working_spectrum,gauss_FWHM_eV) #Still normalized
-    zeroth_order_peak = current_working_spectrum
-    current_full_spectrum += current_working_spectrum #I believe, still normalized
-    norm = 1
-    for i in scatter_num_array:
-        current_working_spectrum = np.load(os.path.join(current_path, 'scatter_spectra_files/scatter'+gas_type+'_'+str(i).zfill(2)+'.npy'))
-        current_working_spectrum = normalize(convolve(zeroth_order_peak,current_working_spectrum,mode='same'))
-        current_full_spectrum += current_working_spectrum*scatter_prob**scatter_num_array[i-1]
-        norm += scatter_prob**scatter_num_array[i-1]
-    # plt.plot(en_array,current_full_spectrum) # diagnostic
-    # plt.show() # diagnostic
-    return current_full_spectrum/norm
-
-# Takes only the nonzero bins of a histogram
-def get_only_nonzero_bins(bins,hist):
-    nonzero_idx = np.where(hist!=0)
-    hist_nonzero = hist[nonzero_idx]
-    hist_err = np.sqrt(hist_nonzero)
-    bins_nonzero = bins[nonzero_idx]
-    return bins_nonzero , hist_nonzero , hist_err
-
-# Flips an array left-to-right. Useful for converting between energy and frequency
-def flip_array(array):
-    flipped = np.fliplr([array]).flatten()
-    return flipped
-
-
-def spectrum_func(x_keV, *p0):
-    logger.info('Using detailed lineshape')
-    x_eV = x_keV*1000.
-    en_loss_array = std_eV_array()
-    en_loss_array_min = en_loss_array[0]
-    en_loss_array_max = en_loss_array[len(en_loss_array)-1]
-    en_array_rev = flip_array(-1*en_loss_array)
-    f = np.zeros(len(x_keV))
-
-    FWHM_G_eV = p0[0]
-    line_pos_keV = p0[1]
-
-    line_pos_eV = line_pos_keV*1000.
-    x_eV_minus_line = x_eV - line_pos_eV
-    zero_idx = np.r_[np.where(x_eV_minus_line<-1*en_loss_array_max)[0],np.where(x_eV_minus_line>-1*en_loss_array_min)[0]]
-    nonzero_idx = [i for i,_ in enumerate(x_keV) if i not in zero_idx]
-
-    for gas_index in range(len(gases)):
-        gas_type = gases[gas_index]
-        scatter_prob = p0[2*gas_index+2]
-        amplitude    = p0[2*gas_index+3]
-
-        full_spectrum = make_spectrum(FWHM_G_eV,scatter_prob, gas_type)
-        full_spectrum_rev = flip_array(full_spectrum)
-        f[nonzero_idx] += amplitude*np.interp(x_eV_minus_line[nonzero_idx],en_array_rev,full_spectrum_rev)
-    return f
-

From 407ba351ca9f90be0604da4fc149f0d95f182575 Mon Sep 17 00:00:00 2001
From: Mathieu Guigue <mathieu.guigue@sorbonne-universite.fr>
Date: Mon, 26 Oct 2020 15:34:29 +0100
Subject: [PATCH 35/49] Update setup.py

---
 setup.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/setup.py b/setup.py
index 56aa86ef..170d5474 100644
--- a/setup.py
+++ b/setup.py
@@ -20,7 +20,7 @@
 
 requirements = []
 extras_require = {
-    'core':['colorlog'],
+    'core':['colorlog', 'iminuit'],
     'doc': ['sphinx','sphinx_rtd_theme','sphinxcontrib-programoutput', 'six', 'colorlog']
 }
 

From c78f987000e8595c4bb25cc39db67307f14edc7b Mon Sep 17 00:00:00 2001
From: cclaessens <claessens.christine@hotmail.de>
Date: Mon, 26 Oct 2020 15:56:10 +0100
Subject: [PATCH 36/49] fix type error in IO_test script

---
 tests/IO_test.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/tests/IO_test.py b/tests/IO_test.py
index 7a1050b4..0ea2cd39 100644
--- a/tests/IO_test.py
+++ b/tests/IO_test.py
@@ -75,9 +75,9 @@ def test_MutliChannelReader(self):
         plt.ylabel('N')
 
         plt.subplot(212)
-        plt.hist(data['a']['TrueStartFrequenciesCut']-24.5e9, bins=bins-24.5e9, label='channel a: {} counts'.format(len(data['a']['TrueStartFrequenciesCut'])))
-        plt.hist(data['b']['TrueStartFrequenciesCut']-24.5e9, bins=bins-24.5e9, label='channel b: {} counts'.format(len(data['b']['TrueStartFrequenciesCut'])))
-        plt.hist(data['c']['TrueStartFrequenciesCut']-24.5e9, bins=bins-24.5e9, label='channel c: {} counts'.format(len(data['c']['TrueStartFrequenciesCut'])))
+        plt.hist(np.array(data['a']['TrueStartFrequenciesCut'])-24.5e9, bins=bins-24.5e9, label='channel a: {} counts'.format(len(data['a']['TrueStartFrequenciesCut'])))
+        plt.hist(np.array(data['b']['TrueStartFrequenciesCut'])-24.5e9, bins=bins-24.5e9, label='channel b: {} counts'.format(len(data['b']['TrueStartFrequenciesCut'])))
+        plt.hist(np.array(data['c']['TrueStartFrequenciesCut'])-24.5e9, bins=bins-24.5e9, label='channel c: {} counts'.format(len(data['c']['TrueStartFrequenciesCut'])))
         plt.xlabel('Start frequencies')
         plt.ylabel('N')
         plt.legend()

From 8a2079a0558799f1796d48c0cd44e63575c791fc Mon Sep 17 00:00:00 2001
From: Mathieu GUIGUE <mgguigue@gmail.com>
Date: Mon, 26 Oct 2020 19:02:49 +0100
Subject: [PATCH 37/49] Adding a iminiut super simple tester which should pass
 tests while the processor test doesn't...

---
 tests/Fitters_test.py | 16 ++++++++++++++++
 1 file changed, 16 insertions(+)

diff --git a/tests/Fitters_test.py b/tests/Fitters_test.py
index c54e78d5..6b2bde6b 100644
--- a/tests/Fitters_test.py
+++ b/tests/Fitters_test.py
@@ -15,6 +15,22 @@
 
 class FittersTest(unittest.TestCase):
 
+    def test_iminuit(self):
+        logger.info('iMinuit test')
+        from iminuit import Minuit
+
+        def f(x, y, z):
+            return (x - 2) ** 2 + (y - 3) ** 2 + (z - 4) ** 2
+
+        m = Minuit(f)
+
+        m.migrad()  # run optimiser
+        print(m.values)  # {'x': 2,'y': 3,'z': 4}
+
+        m.hesse()   # run covariance estimator
+        print(m.errors)  # {'x': 1,'y': 1,'z': 1}
+
+        logger.info('iMinuit test done')
 
     def test_BinnedDataFitter(self):
         from mermithid.processors.Fitters import BinnedDataFitter

From efc315ddf9d0db4e6d31cece41e6b9e9c26fdb58 Mon Sep 17 00:00:00 2001
From: cclaessens <claessens.christine@hotmail.de>
Date: Tue, 27 Oct 2020 13:45:39 +0100
Subject: [PATCH 38/49] added iminuit test more similar to BinnedDataFitter

---
 mermithid/processors/Fitters/BinnedDataFitter.py | 1 +
 tests/Fitters_test.py                            | 7 +++++++
 2 files changed, 8 insertions(+)

diff --git a/mermithid/processors/Fitters/BinnedDataFitter.py b/mermithid/processors/Fitters/BinnedDataFitter.py
index dbaccc0f..21165381 100644
--- a/mermithid/processors/Fitters/BinnedDataFitter.py
+++ b/mermithid/processors/Fitters/BinnedDataFitter.py
@@ -114,6 +114,7 @@ def fit(self):
             logger.info('Fit parameters: {}'.format(self.parameter_names))
             logger.info('Initial values: {}'.format(self.initial_values))
             logger.info('Initial error: {}'.format(self.parameter_errors))
+            logger.info('Limits: {}'.format(self.limits))
             logger.info('Fixed in fit: {}'.format(self.fixes))
             logger.info('Constrained parameters: {}'.format([self.parameter_names[i] for i in self.constrained_parameters]))
 
diff --git a/tests/Fitters_test.py b/tests/Fitters_test.py
index 6b2bde6b..b20973f6 100644
--- a/tests/Fitters_test.py
+++ b/tests/Fitters_test.py
@@ -22,6 +22,9 @@ def test_iminuit(self):
         def f(x, y, z):
             return (x - 2) ** 2 + (y - 3) ** 2 + (z - 4) ** 2
 
+        def g(params):
+            return f(*params)
+
         m = Minuit(f)
 
         m.migrad()  # run optimiser
@@ -30,6 +33,10 @@ def f(x, y, z):
         m.hesse()   # run covariance estimator
         print(m.errors)  # {'x': 1,'y': 1,'z': 1}
 
+        m2 = Minuit.from_array_func(g, [0, 0, 0], error=[0.1, 0.1, 0.1], name=['a', 'b', 'c'], errordef=1)
+        m2.migrad()
+        print(m2.values)
+
         logger.info('iMinuit test done')
 
     def test_BinnedDataFitter(self):

From 56b51617c7bb1feb51f2c01f3369da698111f34d Mon Sep 17 00:00:00 2001
From: cclaessens <claessens.christine@hotmail.de>
Date: Tue, 27 Oct 2020 19:55:46 +0100
Subject: [PATCH 39/49] print_level=1 was causing the problem

---
 mermithid/processors/Fitters/BinnedDataFitter.py |  2 +-
 tests/Fitters_test.py                            | 10 +++++++---
 2 files changed, 8 insertions(+), 4 deletions(-)

diff --git a/mermithid/processors/Fitters/BinnedDataFitter.py b/mermithid/processors/Fitters/BinnedDataFitter.py
index 21165381..13bd65ea 100644
--- a/mermithid/processors/Fitters/BinnedDataFitter.py
+++ b/mermithid/processors/Fitters/BinnedDataFitter.py
@@ -131,7 +131,7 @@ def fit(self):
 
         # minimze
         m_binned.migrad(resume=False)
-        self.param_states = m_binned.get_param_states()
+        #self.param_states = m_binned.get_param_states()
         self.m_binned = m_binned
 
         # results
diff --git a/tests/Fitters_test.py b/tests/Fitters_test.py
index b20973f6..db781c6d 100644
--- a/tests/Fitters_test.py
+++ b/tests/Fitters_test.py
@@ -33,8 +33,11 @@ def g(params):
         m.hesse()   # run covariance estimator
         print(m.errors)  # {'x': 1,'y': 1,'z': 1}
 
-        m2 = Minuit.from_array_func(g, [0, 0, 0], error=[0.1, 0.1, 0.1], name=['a', 'b', 'c'], errordef=1)
-        m2.migrad()
+        # repeat using from_array_func
+        m2 = Minuit.from_array_func(g, [0, 0, 0], error=[0.1, 0.1, 0.1], name=['a1', 'b1', 'c1'], errordef=1,
+                                    limit=[[None, None], [None, None], [None, None]],
+                                    print_level=0)
+        m2.migrad(resume=False)
         print(m2.values)
 
         logger.info('iMinuit test done')
@@ -52,7 +55,8 @@ def test_BinnedDataFitter(self):
             'constrained_parameter_indices': [],
             'constrained_parameter_means': [0.5],
             'constrained_parameter_widths': [1],
-            'binned_data': True
+            'binned_data': True,
+            'print_level': 0
 
             }
 

From 5d06520b6a97f41a27a5ab42314d748fab3fc29c Mon Sep 17 00:00:00 2001
From: cclaessens <claessens.christine@hotmail.de>
Date: Mon, 16 Nov 2020 11:55:55 +0100
Subject: [PATCH 40/49] added logger error output in BinnedDataFitter and made
 min_energy in FakeDAtaGenerator configurable.

---
 mermithid/processors/Fitters/BinnedDataFitter.py          | 5 +++--
 mermithid/processors/TritiumSpectrum/FakeDataGenerator.py | 5 +++--
 2 files changed, 6 insertions(+), 4 deletions(-)

diff --git a/mermithid/processors/Fitters/BinnedDataFitter.py b/mermithid/processors/Fitters/BinnedDataFitter.py
index 13bd65ea..6c5d03dd 100644
--- a/mermithid/processors/Fitters/BinnedDataFitter.py
+++ b/mermithid/processors/Fitters/BinnedDataFitter.py
@@ -108,7 +108,7 @@ def model(self, x, A, mu, sigma):
 
     def fit(self):
         # Now minimize neg log likelihood using iMinuit
-        if self.print_level == 1:
+        if self.print_level > 0:
             logger.info('This is the plan:')
             logger.info('Fitting data consisting of {} elements'.format(np.sum(self.hist)))
             logger.info('Fit parameters: {}'.format(self.parameter_names))
@@ -155,7 +155,8 @@ def PoissonLogLikelihood(self, params):
             expectation = model_return
 
         if np.min(expectation) < 0:
-            logger.error('Expectation contains negative numbers. They will be excluded but something could be horribly wrong.')
+            logger.error('Expectation contains negative numbers: Minimum {} -> {}. They will be excluded but something could be horribly wrong.'.format(np.argmin(expectation), np.min(expectation)))
+            logger.error('FYI, the parameters are: {}'.format(params))
 
         # exclude bins where expectation is <= zero or nan
         index = np.where(expectation>0)
diff --git a/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py b/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
index 9c5e80f9..7ba04d48 100644
--- a/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
+++ b/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py
@@ -99,6 +99,7 @@ def InternalConfigure(self, params):
         self.scattering_sigma = reader.read_param(params, 'scattering_sigma', 18.6)
         self.NScatters = reader.read_param(params, 'NScatters', 20)
         self.scatter_proportion = reader.read_param(params, 'scatter_proportion', 1.0)
+        self.min_energy = reader.read_param(params,'min_lineshape_energy', -1000)
 
         #paths
         self.simplified_scattering_path = reader.read_param(params, 'simplified_scattering_path', '/host/input_data/simplified_scattering_params.txt')
@@ -283,10 +284,10 @@ def generate_unbinned_data(self, Q_mean, mass, ROIbound, S, B_1kev, nsteps=10**4
         FWHM_convert = 2*math.sqrt(2*math.log(2))
         if lineshape=='gaussian':
             max_energy = nstdevs*params[0]
-            min_energy = -1000
+            min_energy = self.min_energy
         elif lineshape=='simplified_scattering' or lineshape=='simplified' or lineshape=='detailed_scattering' or lineshape=='detailed':
             max_energy = nstdevs/FWHM_convert*params[0]
-            min_energy = -1000
+            min_energy = self.min_energy
 
         Kmax_eff = Kmax+max_energy #Maximum energy for data is slightly above Kmax>Q-m
         Kmin_eff = Kmin+min_energy #Minimum is slightly below Kmin<Q-m

From f105e4192ac68eb1fda5ae75cbdd95e90bc216ca Mon Sep 17 00:00:00 2001
From: cclaessens <claessens.christine@hotmail.de>
Date: Tue, 5 Jan 2021 18:31:54 +0100
Subject: [PATCH 41/49] removed commented plots

---
 mermithid/misc/FakeTritiumDataFunctions.py | 18 ------------------
 1 file changed, 18 deletions(-)

diff --git a/mermithid/misc/FakeTritiumDataFunctions.py b/mermithid/misc/FakeTritiumDataFunctions.py
index feb95359..5e36fe6e 100644
--- a/mermithid/misc/FakeTritiumDataFunctions.py
+++ b/mermithid/misc/FakeTritiumDataFunctions.py
@@ -14,8 +14,6 @@
 from scipy.interpolate import interp1d
 from scipy.signal import convolve
 
-import matplotlib.pyplot as plt
-
 from morpho.utilities import morphologging
 logger = morphologging.getLogger(__name__)
 
@@ -167,27 +165,11 @@ def spectral_rate(K, Q, mnu, final_state_array):
 
         Q_states = Q+final_state_array[0]-np.max(final_state_array[0])
 
-        # plt.figure()
-        # plt.plot(Q_states, final_state_array[1])
-        # plt.plot(Q+final_state_array[0], final_state_array[1])
-        # plt.savefig('q_states.png')
-
         index = [np.where(K < Q_states[i]-mnu) for i in range(N_states)]
 
         beta_rates_array = [beta_rates(K, Q_states[i], mnu, index[i])*final_state_array[1][i] for i in range(N_states)]
         to_return = np.nansum(beta_rates_array, axis=0)/np.nansum(final_state_array[1])
 
-
-        # comparison = beta_rates(K, Q, mnu, np.where(K < Q - mnu))
-
-        # plt.figure()
-        # plt.plot(K[comparison>0], 1-to_return[comparison>0]/comparison[comparison>0], marker='', linestyle='-', markersize=1)
-        # plt.ylabel('Final states spectrum relative difference')
-        # plt.xlabel('Energy')
-        # plt.xlim(Q-100, Q+10)
-        # plt.tight_layout()
-        # #plt.yscale('log')
-        # plt.savefig('final_state_included_spectrum.png')
         return to_return
 
     else:

From 888cba466182c92b272de72bd1ffd4a659543769 Mon Sep 17 00:00:00 2001
From: cclaessens <claessens.christine@hotmail.de>
Date: Tue, 5 Jan 2021 20:18:46 +0100
Subject: [PATCH 42/49] iminuit.from_array_func and others were removed in
 iminuit >= 2

---
 .../processors/Fitters/BinnedDataFitter.py    | 26 ++++++++++++-------
 tests/Fitters_test.py                         |  9 ++++---
 2 files changed, 22 insertions(+), 13 deletions(-)

diff --git a/mermithid/processors/Fitters/BinnedDataFitter.py b/mermithid/processors/Fitters/BinnedDataFitter.py
index 6c5d03dd..e57c2f5a 100644
--- a/mermithid/processors/Fitters/BinnedDataFitter.py
+++ b/mermithid/processors/Fitters/BinnedDataFitter.py
@@ -118,25 +118,33 @@ def fit(self):
             logger.info('Fixed in fit: {}'.format(self.fixes))
             logger.info('Constrained parameters: {}'.format([self.parameter_names[i] for i in self.constrained_parameters]))
 
-        m_binned = Minuit.from_array_func(self.negPoissonLogLikelihood,
+        m_binned = Minuit(self.negPoissonLogLikelihood,
                                       self.initial_values,
-                                      error=self.parameter_errors,
-                                      errordef = 0.5, limit = self.limits,
+        #                              error=self.parameter_errors,
+        #                              errordef = 0.5, limit = self.limits,
                                       name=self.parameter_names,
-                                      fix=self.fixes,
-                                      print_level=self.print_level,
-                                      throw_nan=True
+        #                              fix=self.fixes,
+        #                              print_level=self.print_level,
+        #                              throw_nan=True
                                       )
 
+        m_binned.errordef = 0.5
+        m_binned.errors = self.parameter_errors
+        m_binned.throw_nan = True
+        m_binned.print_level = self.print_level
+
+        for i, name in enumerate(self.parameter_names):
+            m_binned.fixed[name] = self.fixes[i]
+            m_binned.limits[name] = self.limits[i]
 
         # minimze
-        m_binned.migrad(resume=False)
+        m_binned.migrad()
         #self.param_states = m_binned.get_param_states()
         self.m_binned = m_binned
 
         # results
-        result_array = m_binned.np_values()
-        error_array = m_binned.np_errors()
+        result_array = np.array(m_binned.values)
+        error_array = np.array(m_binned.errors)
 
         if self.print_level == 1:
             logger.info('Fit results: {}'.format(result_array))
diff --git a/tests/Fitters_test.py b/tests/Fitters_test.py
index db781c6d..77a08332 100644
--- a/tests/Fitters_test.py
+++ b/tests/Fitters_test.py
@@ -34,10 +34,11 @@ def g(params):
         print(m.errors)  # {'x': 1,'y': 1,'z': 1}
 
         # repeat using from_array_func
-        m2 = Minuit.from_array_func(g, [0, 0, 0], error=[0.1, 0.1, 0.1], name=['a1', 'b1', 'c1'], errordef=1,
-                                    limit=[[None, None], [None, None], [None, None]],
-                                    print_level=0)
-        m2.migrad(resume=False)
+        m2 = Minuit(g, [0, 0, 0], name=['a1', 'b1', 'c1'])
+        m2.errors = [0.1, 0.1, 0.1]
+        m2.errordef = 1
+        m2.print_level = 0
+        m2.migrad()
         print(m2.values)
 
         logger.info('iMinuit test done')

From 0d43ce8427ba1c8a78d5586ff819f269239b896e Mon Sep 17 00:00:00 2001
From: cclaessens <claessens.christine@hotmail.de>
Date: Tue, 5 Jan 2021 20:27:41 +0100
Subject: [PATCH 43/49] start values required...

---
 tests/Fitters_test.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/Fitters_test.py b/tests/Fitters_test.py
index 77a08332..7029d18b 100644
--- a/tests/Fitters_test.py
+++ b/tests/Fitters_test.py
@@ -25,7 +25,7 @@ def f(x, y, z):
         def g(params):
             return f(*params)
 
-        m = Minuit(f)
+        m = Minuit(f, x=0, y=0, z=0)
 
         m.migrad()  # run optimiser
         print(m.values)  # {'x': 2,'y': 3,'z': 4}

From 1afb3dac0580c7dcbeacfd1fe14d21fc06dc96b1 Mon Sep 17 00:00:00 2001
From: cclaessens <claessens.christine@hotmail.de>
Date: Thu, 7 Jan 2021 10:12:31 +0100
Subject: [PATCH 44/49] remove KrLineShapeFunctions from misc/__init__.py

---
 mermithid/misc/__init__.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/mermithid/misc/__init__.py b/mermithid/misc/__init__.py
index 0dc77499..afb4ebbe 100644
--- a/mermithid/misc/__init__.py
+++ b/mermithid/misc/__init__.py
@@ -8,4 +8,3 @@
 from . import TritiumFormFactor
 from . import FakeTritiumDataFunctions
 from . import ConversionFunctions
-from . import KrLineshapeFunctions

From ec131351f83de96a74c69a698025e5856a7ad499 Mon Sep 17 00:00:00 2001
From: cclaessens <claessens.christine@hotmail.de>
Date: Fri, 14 May 2021 16:11:44 -0700
Subject: [PATCH 45/49] Default settings were a problem if fit function had not
 exactly 3 fit parameters.

I also added the correlation matrix to the return and the call of
simplex() before migrad().
---
 mermithid/processors/Fitters/BinnedDataFitter.py | 12 +++++++-----
 1 file changed, 7 insertions(+), 5 deletions(-)

diff --git a/mermithid/processors/Fitters/BinnedDataFitter.py b/mermithid/processors/Fitters/BinnedDataFitter.py
index e57c2f5a..edd2faac 100644
--- a/mermithid/processors/Fitters/BinnedDataFitter.py
+++ b/mermithid/processors/Fitters/BinnedDataFitter.py
@@ -54,13 +54,12 @@ def InternalConfigure(self, params):
         # Read other parameters
         self.namedata = reader.read_param(params, 'variables', "required")
         self.parameter_names = reader.read_param(params, 'parameter_names', ['A', 'mu', 'sigma'])
-        self.initial_values = reader.read_param(params, 'initial_values', [1, 0, 1])
-        self.limits = reader.read_param(params, 'limits', [[None, None], [None, None], [None, None]])
-        self.fixes = reader.read_param(params, 'fixed', [False, False, False])
+        self.initial_values = reader.read_param(params, 'initial_values', [1]*len(self.parameter_names))
+        self.limits = reader.read_param(params, 'limits', [[None, None]]*len(self.parameter_names))
+        self.fixes = reader.read_param(params, 'fixed', [False]*len(self.parameter_names))
         self.bins = reader.read_param(params, 'bins', np.linspace(-2, 2, 100))
         self.binned_data = reader.read_param(params, 'binned_data', False)
         self.print_level = reader.read_param(params, 'print_level', 1)
-
         self.constrained_parameters = reader.read_param(params, 'constrained_parameter_indices', [])
         self.constrained_means = reader.read_param(params, 'constrained_parameter_means', [])
         self.constrained_widths = reader.read_param(params, 'constrained_parameter_widths', [])
@@ -90,6 +89,7 @@ def InternalRun(self):
         self.results = {}
         self.results['param_values'] = result_array
         self.results['param_errors'] = error_array
+        self.results['correlation_matrix'] = np.array(self.m_binned.covariance.correlation())
         for i, k in enumerate(self.parameter_names):
             self.results[k] = {'value': result_array[i], 'error': error_array[i]}
 
@@ -138,7 +138,8 @@ def fit(self):
             m_binned.limits[name] = self.limits[i]
 
         # minimze
-        m_binned.migrad()
+        m_binned.simplex().migrad()
+        m_binned.hesse()
         #self.param_states = m_binned.get_param_states()
         self.m_binned = m_binned
 
@@ -149,6 +150,7 @@ def fit(self):
         if self.print_level == 1:
             logger.info('Fit results: {}'.format(result_array))
             logger.info('Errors: {}'.format(error_array))
+
         return result_array, error_array
 
 

From ee8e58b3ac9740f0782cb9b5fe74fd6bffe75dba Mon Sep 17 00:00:00 2001
From: Florian Thomas <flthomas@students.uni-mainz.de>
Date: Wed, 30 Jun 2021 00:07:09 +0200
Subject: [PATCH 46/49] Add workflow for publishing the docker image

---
 .github/workflows/publish.yaml | 155 +++++++++++++++++++++++++++++++++
 Dockerfile                     |   9 +-
 2 files changed, 162 insertions(+), 2 deletions(-)
 create mode 100644 .github/workflows/publish.yaml

diff --git a/.github/workflows/publish.yaml b/.github/workflows/publish.yaml
new file mode 100644
index 00000000..d4feced7
--- /dev/null
+++ b/.github/workflows/publish.yaml
@@ -0,0 +1,155 @@
+
+name: Publish
+
+on:
+
+  push:
+    branches: [ master, main, develop, feature/github-actions]
+    tags: 'v*.*.*'
+  
+  release:
+    types: [published]
+
+  workflow_dispatch:
+
+jobs:
+
+  init:
+    runs-on: ubuntu-latest
+    outputs:
+      build_type: ${{ steps.setvars.outputs.build_type }}
+      build_tests: ${{ steps.setvars.outputs.build_tests }}
+
+    steps:
+      - name: Cancel previous workflow
+        uses: styfle/cancel-workflow-action@0.4.0
+        with:
+          access_token: ${{ github.token }}
+
+      - name: Set variables
+        id: setvars
+        #Only build the tests in a debug build, tests currently don't build in release mode
+        run: |
+          echo ${{github.base_ref}}
+          echo ${{github.ref}}
+          if [[ "${{github.base_ref}}" == "develop" || "${{github.ref}}" == "refs/heads/develop" ]]; then
+            echo "::set-output name=build_type::DEBUG"
+          else
+            echo "::set-output name=build_type::RELEASE"
+          fi
+
+  #~ test_docker:
+    #~ runs-on: ubuntu-latest
+    #~ steps:
+      #~ - name: Checkout the repo 
+        #~ uses: actions/checkout@v2
+        #~ with:
+          #~ submodules: recursive
+      #~ - name: Build and run
+        #~ run: |
+          #~ cd ${{github.workspace}}
+          #~ docker build \
+            #~ --build-arg IMG_USER=project8 \
+            #~ --build-arg IMG_REPO=p8compute_dependencies \
+            #~ --build-arg IMG_TAG=v1.0.0.beta \
+            #~ --build-arg PSYLLID_BASE_PREFIX=/usr/local/p8/psyllid \
+            #~ --build-arg PSYLLID_PREFIX_TAG=test \
+            #~ --tag project8/psyllid:test \
+            #~ .
+          
+          
+  #~ test_linux_gcc:
+    #~ runs-on: ubuntu-latest
+
+    #~ steps:
+      #~ - name: Checkout the repo 
+        #~ uses: actions/checkout@v2
+        #~ with:
+          #~ submodules: recursive
+      #~ - name: Install dependencies
+        #~ run: |
+          #~ sudo apt-get update
+          #~ sudo apt-get clean \
+          #~ sudo apt-get --fix-missing -y install \
+            #~ build-essential \
+            #~ cmake \
+            #~ libfftw3-3 \
+            #~ libfftw3-dev \
+            #~ gdb \
+            #~ libboost-all-dev \
+            #~ libhdf5-dev \
+            #~ librabbitmq-dev \
+            #~ rapidjson-dev \
+            #~ libyaml-cpp-dev \
+            #~ pybind11-dev \
+            #~ wget \
+      #~ - name: Configure CMake
+      #~ # Configure CMake in a 'build' subdirectory. `CMAKE_BUILD_TYPE` is only required if you are using a single-configuration generator such as make.
+      #~ # See https://cmake.org/cmake/help/latest/variable/CMAKE_BUILD_TYPE.html?highlight=cmake_build_type
+        #~ run: cmake -B ${{github.workspace}}/build -DCMAKE_BUILD_TYPE=${{env.BUILD_TYPE}}
+      #~ - name: Build
+        #~ run: |
+          #~ cd ${{github.workspace}}/build
+          #~ make -j3 install
+          
+      #TODO add tests for Psyllid
+      #~ - name: Run tests
+        #~ run: |
+          #~ cd build/testing
+          #~ ./run_tests
+
+  publish:
+    name: Push Docker image to Docker Hub
+    runs-on: ubuntu-latest
+   # needs: [test_docker, test_linux_gcc]
+    needs: init
+    steps:
+    
+      - name: Check out the repo
+        uses: actions/checkout@v2
+        with:
+          submodules: recursive
+          
+      - name: Docker meta
+        id: docker_meta
+        uses: crazy-max/ghaction-docker-meta@v1
+        with:
+          images: project8/locust_mc
+          tag-sha: false
+          tag-semver: |
+            {{raw}}
+          
+      - name: Set up QEMU
+        uses: docker/setup-qemu-action@v1
+        
+      - name: Set up Docker Buildx
+        id: setup_buildx
+        uses: docker/setup-buildx-action@v1
+        with:
+          buildkitd-flags: --debug
+          
+      - name: List platforms
+        run: echo ${{ steps.setup_buildx.outputs.platforms }}
+        
+      - name: Log in to Docker Hub
+        uses: docker/login-action@v1
+        with:
+          username: ${{ secrets.DOCKER_USERNAME }}
+          password: ${{ secrets.DOCKER_PASSWORD }}
+          
+      - name: Set env
+        run: echo "RELEASE_VERSION=${GITHUB_REF#refs/*/}" >> $GITHUB_ENV
+        
+      - name: Push to Docker Hub
+        id: build_push
+        uses: docker/build-push-action@v2
+        with:
+          context: .
+          push: true
+          build-args: |
+            IMG_USER=project8
+            IMG_REPO=p8compute_dependencies
+            IMG_TAG=v1.0.0.beta
+            MERMITHID_TAG=${{ env.RELEASE_VERSION }}
+            build_type=${{needs.init.outputs.build_type}}
+          tags: ${{ steps.docker_meta.outputs.tags }}
diff --git a/Dockerfile b/Dockerfile
index dbde3af5..1e3bde52 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -1,9 +1,14 @@
-FROM project8/p8compute_dependencies:v0.9.0 as mermithid_common
+ARG IMG_USER=project8
+ARG IMG_REPO=p8compute_dependencies
+ARG IMG_TAG=v1.0.0.beta
+
+FROM ${IMG_USER}/${IMG_REPO}:${IMG_TAG} as mermithid_common
 
 ARG build_type=Release
 ENV MERMITHID_BUILD_TYPE=$build_type
 
-ENV MERMITHID_TAG=v1.2.2
+ARG MERMITHID_TAG=beta
+ENV MERMITHID_TAG=${MERMITHID_TAG}
 ENV MERMITHID_BUILD_PREFIX=/usr/local/p8/mermithid/$MERMITHID_TAG
 
 RUN mkdir -p $MERMITHID_BUILD_PREFIX &&\

From 242fcf11890d542e12d7b10d247fd429e15bb203 Mon Sep 17 00:00:00 2001
From: Florian Thomas <flthomas@students.uni-mainz.de>
Date: Fri, 2 Jul 2021 15:46:23 +0200
Subject: [PATCH 47/49] Update to full workflow

---
 .github/workflows/publish.yaml | 170 ++++++++++++++++++++-------------
 Dockerfile                     |   5 +
 2 files changed, 106 insertions(+), 69 deletions(-)

diff --git a/.github/workflows/publish.yaml b/.github/workflows/publish.yaml
index d4feced7..ed621308 100644
--- a/.github/workflows/publish.yaml
+++ b/.github/workflows/publish.yaml
@@ -4,13 +4,20 @@ name: Publish
 on:
 
   push:
-    branches: [ master, main, develop, feature/github-actions]
+    branches: [ main, develop, feature/github-actions]
     tags: 'v*.*.*'
+    
+  pull_request:
+    branches: '*'
+    types: opened
   
   release:
     types: [published]
 
   workflow_dispatch:
+  
+  repository_dispatch:
+    types: [release-event]
 
 jobs:
 
@@ -18,7 +25,7 @@ jobs:
     runs-on: ubuntu-latest
     outputs:
       build_type: ${{ steps.setvars.outputs.build_type }}
-      build_tests: ${{ steps.setvars.outputs.build_tests }}
+      luna_version: ${{ steps.setversion.outputs.luna_version }}
 
     steps:
       - name: Cancel previous workflow
@@ -34,75 +41,67 @@ jobs:
           echo ${{github.ref}}
           if [[ "${{github.base_ref}}" == "develop" || "${{github.ref}}" == "refs/heads/develop" ]]; then
             echo "::set-output name=build_type::DEBUG"
+            echo Debug build
           else
             echo "::set-output name=build_type::RELEASE"
           fi
-
-  #~ test_docker:
-    #~ runs-on: ubuntu-latest
-    #~ steps:
-      #~ - name: Checkout the repo 
-        #~ uses: actions/checkout@v2
-        #~ with:
-          #~ submodules: recursive
-      #~ - name: Build and run
-        #~ run: |
-          #~ cd ${{github.workspace}}
-          #~ docker build \
-            #~ --build-arg IMG_USER=project8 \
-            #~ --build-arg IMG_REPO=p8compute_dependencies \
-            #~ --build-arg IMG_TAG=v1.0.0.beta \
-            #~ --build-arg PSYLLID_BASE_PREFIX=/usr/local/p8/psyllid \
-            #~ --build-arg PSYLLID_PREFIX_TAG=test \
-            #~ --tag project8/psyllid:test \
-            #~ .
           
-          
-  #~ test_linux_gcc:
-    #~ runs-on: ubuntu-latest
+      - name: Set Luna version
+        id: setversion
+        run: |
+          git clone https://github.com/project8/luna_base.git
+          cd luna_base
+          echo "::set-output name=luna_version::$(git describe --abbrev=0 --tags)"
+
+  test_linux:
+    runs-on: ubuntu-latest
+    needs: [init]
+    
+    strategy:
+      matrix:
+        build_type: [DEBUG, RELEASE]
+        compiler: [gcc]
+      fail-fast: false
+        
+    steps:
+      - name: Checkout the repo 
+        uses: actions/checkout@v2
+        with:
+          submodules: recursive
 
-    #~ steps:
-      #~ - name: Checkout the repo 
-        #~ uses: actions/checkout@v2
-        #~ with:
-          #~ submodules: recursive
-      #~ - name: Install dependencies
-        #~ run: |
-          #~ sudo apt-get update
-          #~ sudo apt-get clean \
-          #~ sudo apt-get --fix-missing -y install \
-            #~ build-essential \
-            #~ cmake \
-            #~ libfftw3-3 \
-            #~ libfftw3-dev \
-            #~ gdb \
-            #~ libboost-all-dev \
-            #~ libhdf5-dev \
-            #~ librabbitmq-dev \
-            #~ rapidjson-dev \
-            #~ libyaml-cpp-dev \
-            #~ pybind11-dev \
-            #~ wget \
-      #~ - name: Configure CMake
-      #~ # Configure CMake in a 'build' subdirectory. `CMAKE_BUILD_TYPE` is only required if you are using a single-configuration generator such as make.
-      #~ # See https://cmake.org/cmake/help/latest/variable/CMAKE_BUILD_TYPE.html?highlight=cmake_build_type
-        #~ run: cmake -B ${{github.workspace}}/build -DCMAKE_BUILD_TYPE=${{env.BUILD_TYPE}}
-      #~ - name: Build
-        #~ run: |
-          #~ cd ${{github.workspace}}/build
-          #~ make -j3 install
+      - name: Set env
+        run: |
+          if [[ "${{matrix.compiler}}" == "gcc" ]]; then
+            echo "CXX_VAL=g++" >> $GITHUB_ENV
+          else
+            echo "CXX_VAL=clang++" >> $GITHUB_ENV
+          fi
           
-      #TODO add tests for Psyllid
-      #~ - name: Run tests
-        #~ run: |
-          #~ cd build/testing
-          #~ ./run_tests
+      - name: Print config
+        run: |
+            echo Build type: ${{ matrix.build_type }}
+            echo Compiler: ${{ matrix.compiler }} ${{ env.CXX_VAL }} 
+            echo Luna version: ${{ needs.init.outputs.luna_version }}
+          
+      - name: Build with docker
+        run: |
+          cd ${{github.workspace}}
+          docker build \
+            --build-arg IMG_USER=project8 \
+            --build-arg IMG_REPO=p8compute_dependencies \
+            --build-arg IMG_TAG=${{ needs.init.outputs.luna_version }} \
+            --build-arg MERMITHID_TAG=test \
+            --build-arg build_type=${{matrix.build_type}} \
+            --build-arg CC_VAL=${{matrix.compiler}} \
+            --build-arg CXX_VAL=${{env.CXX_VAL}} \
+            --tag project8/mermithid:test \
+            .
 
   publish:
+    if: github.event_name != 'pull_request'
     name: Push Docker image to Docker Hub
     runs-on: ubuntu-latest
-   # needs: [test_docker, test_linux_gcc]
-    needs: init
+    needs: [init, test_linux]
     steps:
     
       - name: Check out the repo
@@ -112,12 +111,9 @@ jobs:
           
       - name: Docker meta
         id: docker_meta
-        uses: crazy-max/ghaction-docker-meta@v1
+        uses: docker/metadata-action@v3
         with:
-          images: project8/locust_mc
-          tag-sha: false
-          tag-semver: |
-            {{raw}}
+          images: project8/mermithid
           
       - name: Set up QEMU
         uses: docker/setup-qemu-action@v1
@@ -129,7 +125,8 @@ jobs:
           buildkitd-flags: --debug
           
       - name: List platforms
-        run: echo ${{ steps.setup_buildx.outputs.platforms }}
+        run: | 
+            echo Platforms: ${{ steps.setup_buildx.outputs.platforms }}
         
       - name: Log in to Docker Hub
         uses: docker/login-action@v1
@@ -138,7 +135,13 @@ jobs:
           password: ${{ secrets.DOCKER_PASSWORD }}
           
       - name: Set env
-        run: echo "RELEASE_VERSION=${GITHUB_REF#refs/*/}" >> $GITHUB_ENV
+        run: | 
+            echo "RELEASE_VERSION=${GITHUB_REF#refs/*/}" >> $GITHUB_ENV
+            
+      - name: Print config
+        run: |
+            echo Release version: ${{ env.RELEASE_VERSION }}
+            echo Luna version: ${{ needs.init.outputs.luna_version }}
         
       - name: Push to Docker Hub
         id: build_push
@@ -149,7 +152,36 @@ jobs:
           build-args: |
             IMG_USER=project8
             IMG_REPO=p8compute_dependencies
-            IMG_TAG=v1.0.0.beta
+            IMG_TAG=${{ needs.init.outputs.luna_version }}
             MERMITHID_TAG=${{ env.RELEASE_VERSION }}
             build_type=${{needs.init.outputs.build_type}}
           tags: ${{ steps.docker_meta.outputs.tags }}
+
+# disabled due to error
+#  doc:
+#    name: Build documentation
+#    runs-on: ubuntu-latest
+#    #needs: [init, test_linux]
+#    steps:
+#    
+#      - name: Check out the repo
+#        uses: actions/checkout@v2
+#        with:
+#          submodules: recursive
+#          fetch-depth: 0
+#      
+#      - name: Install dependencies
+#        run: |
+#          sudo apt-get update
+#          DEBIAN_FRONTEND=noninteractive sudo apt-get install -y \
+#            tree \
+#            doxygen \
+#            python3-sphinx \
+#            graphviz
+#          pip install sphinxcontrib-blockdiag
+#          pip install sphinxcontrib-contentui
+#      
+#      - name: Build docs
+#        run: |
+#          cd Documentation
+#          make
diff --git a/Dockerfile b/Dockerfile
index c99ffbde..4922e950 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -11,6 +11,11 @@ ARG MERMITHID_TAG=beta
 ENV MERMITHID_TAG=${MERMITHID_TAG}
 ENV MERMITHID_BUILD_PREFIX=/usr/local/p8/mermithid/$MERMITHID_TAG
 
+ARG CC_VAL=gcc
+ENV CC=${CC_VAL}
+ARG CXX_VAL=g++
+ENV CXX=${CXX_VAL}
+
 RUN mkdir -p $MERMITHID_BUILD_PREFIX &&\
     chmod -R 777 $MERMITHID_BUILD_PREFIX/.. &&\
     cd $MERMITHID_BUILD_PREFIX &&\

From a6e189cd4a36cd6e311983124ad42af651eb3cf6 Mon Sep 17 00:00:00 2001
From: Florian Thomas <flthomas@students.uni-mainz.de>
Date: Fri, 2 Jul 2021 18:27:34 +0200
Subject: [PATCH 48/49] Add trigger for p8compute build

---
 .github/workflows/publish.yaml | 19 +++++++++++++++++++
 1 file changed, 19 insertions(+)

diff --git a/.github/workflows/publish.yaml b/.github/workflows/publish.yaml
index ed621308..653f2988 100644
--- a/.github/workflows/publish.yaml
+++ b/.github/workflows/publish.yaml
@@ -185,3 +185,22 @@ jobs:
 #        run: |
 #          cd Documentation
 #          make
+
+  notify-packages:
+  
+      if: (github.event_name == 'push' && contains(github.ref, 'refs/tags/') ) 
+            || github.event_name == 'release' 
+            || github.event_name == 'workflow_dispatch'
+            || github.event_name == 'repository_dispatch'
+      name: Trigger new docker P8Compute image
+      runs-on: ubuntu-latest
+      needs: [publish]
+
+      steps:
+      
+        - name: Repository dispatch
+          uses: peter-evans/repository-dispatch@v1
+          with:
+            token: ${{ secrets.REPO_ACCESS_TOKEN }}
+            repository: project8/luna
+            event-type: release-event

From 726eca9511161fb57e4a2b2229dd74b0688504a1 Mon Sep 17 00:00:00 2001
From: Florian Thomas <flthomas@students.uni-mainz.de>
Date: Fri, 2 Jul 2021 18:43:33 +0200
Subject: [PATCH 49/49] Add master branch to triggers of publish workflow

---
 .github/workflows/publish.yaml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/publish.yaml b/.github/workflows/publish.yaml
index 653f2988..007b0cd5 100644
--- a/.github/workflows/publish.yaml
+++ b/.github/workflows/publish.yaml
@@ -4,7 +4,7 @@ name: Publish
 on:
 
   push:
-    branches: [ main, develop, feature/github-actions]
+    branches: [ master, main, develop, feature/github-actions]
     tags: 'v*.*.*'
     
   pull_request: