Combining multi gas scattering with simulated ins resolution in fake data generator

docker-compose.yaml

@@ -7,3 +7,5 @@ services:
     volumes:
       # share a subdirectory from the host to /host in the docker (can be edited)
       - ~/mermithid_share:/host
+      - ~/mermithid_FDG2/mermithid:/host-mermithid
+      - ~:/home_directory

mermithid/misc/ComplexLineShapeUtilities.py

@@ -39,9 +39,12 @@ def read_oscillator_str_file(filename):
 
     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])
+            try:
+                raw_data = [float(i) for i in line.split("\t")]
+                energyOsc[0].append(raw_data[0])
+                energyOsc[1].append(raw_data[1])
+            except:
+                continue
 
     energyOsc = np.array(energyOsc)
     ### take data and sort by energy
@@ -62,14 +65,21 @@ def aseev_func_tail(energy_loss_array, gas_type):
         A2, omeg2, eps2 = 0.1187, 33.40, 10.43
     elif gas_type=="Ar":
         A2, omeg2, eps2 = 0.3344, 21.91, 21.14
+    elif gas_type=="N2":
+        A2, omeg2, eps2 = 0.21754816, 44.99897054, 20.43916114
+    elif gas_type=="CO":
+        A2, omeg2, eps2 = 0.19583454, 55.21888452, 16.44972596
+    elif gas_type=="C2H4":
+        A2, omeg2, eps2 = 0.57492182, 23.77501391, 14.33107345
     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, kr_17keV_line): #energies in eV
-    kinetic_en = kr_17keV_line * 1000
+    kinetic_en = kr_17keV_line
     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)
+    argument_of_log = np.where(e_loss > 0, 4. * kinetic_en * e_rydberg / (e_loss**2.) , 1e-5)
+    return np.where(e_loss>0 , 1./(e_loss) * oscillator_strength* np.log(argument_of_log), 0)
 
 # Takes only the nonzero bins of a histogram
 def get_only_nonzero_bins(bins,hist):
@@ -109,11 +119,10 @@ def energy_guess_to_frequency(energy_guess, energy_guess_err, B_field_guess):
     return frequency , frequency_err
 
 # Given a frequency and error, converts those to B field values assuming the line is the 17.8 keV line
-def central_frequency_to_B_field(central_freq,central_freq_err):
+def central_frequency_to_B_field(central_freq):
     const = (2.*np.pi*m_e)*(1+kr_17keV_line/mass_energy_electron)/e_charge
     B_field = const*central_freq
-    B_field_err = const*central_freq_err
-    return B_field , B_field_err
+    return B_field
 
 # given a FWHM for the lorentian component and the FWHM for the gaussian component,
 # this function estimates the FWHM of the resulting voigt distribution

mermithid/misc/FakeTritiumDataFunctions.py

@@ -274,7 +274,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, scatter_peak_ratio_b, scatter_peak_ratio_c, scatter_fraction, min_energy, max_energy,
                                    complexLineShape, final_state_array):
     """K is an array-like object
     """
@@ -293,23 +293,22 @@ 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 = complexLineShape.make_spectrum_simulated_resolution_scaled_fit_scatter_peak_ratio(1, ls_params[1], scatter_peak_ratio_b, scatter_peak_ratio_c, scatter_fraction, gauss_FWHM_eV=ls_params[0], emitted_peak='dirac')
+        lineshape_rates = np.flipud(lineshape_rates)
 
-        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_state_array) #np.zeros(len(K))
-    #for i,ke in enumerate(K):
-    #    beta_rates[i] = spectral_rate(ke, Q, mnu, final_state_array)
+    beta_rates = spectral_rate(K, Q, mnu, final_state_array)
 
     #Convolving
     convolved = convolve(beta_rates, lineshape_rates, mode='same')
+
     below_Kmin = np.where(K < Kmin)
     np.put(convolved, below_Kmin, np.zeros(len(below_Kmin)))
     return convolved
 
 
 
 #Convolution of background and lineshape using scipy.signal.convolve
-def convolved_bkgd_rate_arrays(K, Kmin, Kmax, lineshape, ls_params, min_energy, max_energy, complexLineShape):
+def convolved_bkgd_rate_arrays(K, Kmin, Kmax, lineshape, ls_params, scatter_peak_ratio_b, scatter_peak_ratio_c, scatter_fraction, min_energy, max_energy, complexLineShape):
     """K is an array-like object
     """
     energy_half_range = max(max_energy, abs(min_energy))
@@ -325,8 +324,9 @@ 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, 1, ls_params[1])
-
+        lineshape_rates = complexLineShape.make_spectrum_simulated_resolution_scaled_fit_scatter_peak_ratio(1, ls_params[1], scatter_peak_ratio_b, scatter_peak_ratio_c, scatter_fraction, gauss_FWHM_eV=ls_params[0], emitted_peak='dirac')
+        lineshape_rates = np.flipud(lineshape_rates)
+
     bkgd_rates = np.full(len(K), bkgd_rate())
     if len(K) < len(K_lineshape):
         raise Exception("lineshape array is longer than Koptions")
@@ -335,30 +335,40 @@ def convolved_bkgd_rate_arrays(K, Kmin, Kmax, lineshape, ls_params, min_energy,
     convolved = convolve(bkgd_rates, lineshape_rates, mode='same')
     below_Kmin = np.where(K < Kmin)
     np.put(convolved, below_Kmin, np.zeros(len(below_Kmin)))
+
     return convolved
 
 
 
-##Fraction of events near the endpoint
-##Currently, this only holds for the last 13.6 eV of the spectrum
-#def frac_near_endpt(Kmin, Q, mass, atom_or_mol='atom'):
-#    A = integrate.quad(spectral_rate, Kmin, Q-mass, args=(Q,mass))
-#    B = integrate.quad(spectral_rate, V0, Q-mass, args=(Q,mass)) #Minimum at V0 because electrons with energy below screening barrier do not escape
-#    f = (A[0])/(B[0])
-#    if atom_or_mol=='atom':
-#        return 0.7006*f
-#    elif atom_or_mol=='mol' or atom_or_mol=='molecule':
-#        return 0.57412*f
-#    else:
-#        print("Choose 'atom' or 'mol'.")
+#Fraction of events near the endpoint
+def frac_near_endpt(Kmin, Q, mass, final_state_array, atom_or_mol='mol', range='wide'):
+    """
+    Options for range:
+        - 'narrow': Only extends ~18 eV (or less) below the endpoint, so that all decays are to the ground state
+        - 'wide': Wide enough that the probability of decay to a 3He electronic energy level that would shift Q below the ROI is very low
+    """
+    A = integrate.quad(spectral_rate, Kmin, Q-mass, args=(Q, mass, final_state_array))
+    B = integrate.quad(spectral_rate, V0, Q-mass, args=(Q, mass, final_state_array)) #Minimum at V0 because electrons with energy below screening barrier do not escape
+    f = (A[0])/(B[0])
+    if range=='narrow':
+        if atom_or_mol=='atom':
+            return 0.7006*f
+        elif atom_or_mol=='mol' or atom_or_mol=='molecule':
+            return 0.57412*f
+        else:
+            logger.warn("Choose 'atom' or 'mol'.")
+    elif range=='wide':
+        return f
+    else:
+        logger.warn("Choose range 'narrow' or 'wide'")
 
 
 #Convert [number of particles]=(density*volume*efficiency) to a signal activity A_s, measured in events/second.
 def find_signal_activity(Nparticles, m, Q, Kmin, atom_or_mol='atom', nTperMolecule=2):
     """
     Functions to calculate number of events to generate
     """
-    br = frac_near_endpt(Kmin, Q, m, atom_or_mol)
+    br = frac_near_endpt(Kmin, Q, m, final_state_array, atom_or_mol)
     Thalflife = 3.8789*10**8
     A_s = Nparticles*np.log(2)/(Thalflife)*br
     if atom_or_mol=='atom':