Feature/post taup 2021 frequentist analysis

.gitignore

@@ -101,3 +101,6 @@ ENV/
 .mypy_cache/
 *.root
 *.png
+
+# CMake build directory
+build*/

CMakeLists.txt

@@ -9,13 +9,13 @@ include( PackageBuilder )
 
 pbuilder_prepare_project()
 
+set_option( Scarab_BUILD_CODEC_JSON OFF )
+set_option( Scarab_BUILD_CODEC_YAML OFF )
+set_option( Scarab_BUILD_AUTHENTICATION OFF )
+set_option( Scarab_BUILD_CLI OFF )
+set_option( Scarab_BUILD_PARAM OFF )
+set_option( Scarab_ENABLE_EXECUTABLES OFF )
+set_option( Cicada_ENABLE_KATYDID_NAMESPACE OFF )
+
 pbuilder_add_submodule( Cicada Cicada )
 pbuilder_add_submodule( Phylloxera Phylloxera )
-
-set( Scarab_BUILD_CODEC_JSON OFF CACHE BOOL "No coded" FORCE)
-set( Scarab_BUILD_CODEC_YAML OFF CACHE BOOL "No codec" FORCE)
-set( Scarab_BUILD_AUTHENTICATION OFF CACHE BOOL "No auth" FORCE)
-set( Scarab_BUILD_PARAM OFF CACHE BOOL "No param" FORCE)
-
-set( Cicada_ENABLE_KATYDID_NAMESPACE OFF CACHE BOOL "Use Cicada namespace" FORCE)
-

Dockerfile

@@ -22,10 +22,6 @@ RUN mkdir -p $MERMITHID_BUILD_PREFIX &&\
     echo 'export PYTHONPATH=$MERMITHID_BUILD_PREFIX/$(python3 -m site --user-site | sed "s%$(python3 -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
 
@@ -45,6 +41,7 @@ COPY tests $MERMITHID_BUILD_PREFIX/tests
 
 # repeat the cmake command to get the change of install prefix to set correctly (a package_builder known issue)
 RUN source $MERMITHID_BUILD_PREFIX/setup.sh &&\
+#    pip3 install --upgrade pip &&\
     cd /tmp_source &&\
     mkdir -p build &&\
     cd build &&\
@@ -64,3 +61,4 @@ RUN source $MERMITHID_BUILD_PREFIX/setup.sh &&\
 FROM mermithid_common
 
 COPY --from=mermithid_done $MERMITHID_BUILD_PREFIX $MERMITHID_BUILD_PREFIX
+RUN pip3 install seaborn statsmodels

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
@@ -50,9 +53,14 @@ def read_oscillator_str_file(filename):
     energyOsc[1] = energyOsc[1][sorted_indices]
     return energyOsc
 
-# A sub function for the scatter function. Found in
+# A sub function for the scatter function. The oscillator strength tails are from
+# this database: https://nl.lxcat.net/home/. 
+# To get the parameters below, the tails are then fitted using the function in:
 # "Energy loss of 18 keV electrons in gaseous T and quench condensed D films"
 # by V.N. Aseev et al. 2000
+# The tails in the LXCAT database sometimes only go out to a certain energy.
+# Sometimes, we need energies extending beyond that value. That is why we
+# need to use this parameterized form of the tail.
 def aseev_func_tail(energy_loss_array, gas_type):
     if gas_type=="H2":
         A2, omeg2, eps2 = 0.195, 14.13, 10.60
@@ -62,14 +70,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 +124,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
@@ -200,4 +214,4 @@ def shake_spectrum(self):
         x_array = flip_array(x_array)
         shake_spectrum = self.full_shake_spectrum(x_array, 0, 24)
         return shake_spectrum
-    ###############################################################################
+    ###############################################################################

mermithid/misc/Constants.py

@@ -32,8 +32,8 @@ def GF(): return 1.1663787*10**(-23) #Gf/(hc)^3, in eV^(-2)
 def Vud(): return 0.97425 #CKM element
 
 #Beta decay-specific physical constants
-def QT(): return 18563.251 #For atomic tritium (eV), from Bodine et al. (2015)
-def QT2(): return 18573.24 #For molecular tritium (eV), Bodine et al. (2015)
+def QT(): return 18563.251 #SHOULD BE DOUBLE-CHECKED. For atomic tritium (eV), from Bodine et al. (2015) 
+def QT2(): return 18574.01 #For molecular tritium (eV). Calculation here: https://projecteight.slack.com/archives/CG5TY2UE7/p1649963449399179, based on Bodine et al. (2015).
 def Rn(): return 2.8840*10**(-3) #Helium-3 nuclear radius in units of me, from Kleesiek et al. (2018): https://arxiv.org/pdf/1806.00369.pdf
 def M_3He_in_me(): return 5497.885 #Helium-3 mass in units of me, Kleesiek et al. (2018)
 def atomic_num(): return 2. #For helium-3