Atomic effects

README.md

@@ -15,6 +15,8 @@ If you are not using a container with pre-installed dependencies, you will need
 
 ## Installation
 
+After cloning mermithid, run `git submodule update --init --recursive`.
+
 These are two possible ways of installing and working with mermithid.
 
 ### Virtual environment installation

docker-compose.yaml

@@ -1,6 +1,19 @@
+
 #version: "3.9"
-#Check running containers docker ps (detached) -d
-#Attach to container docker attach [name]
+
+# Check running containers docker ps (detached) -d
+#docker compose up -d 
+
+# Attach to container docker attach [name]
+#docker attach mermithid_sensitivity
+
+# Mermithid is a bind mount, this guarantees: live code edits on host, Python imports always resolve to your branch & no accidental use of the installed site-packages version
+#export PYTHONPATH=/mermithid/:$PYTHONPATH &&
+
+# /mermithid is on sys.path. Python should import from it but it does not, because mermithid is already imported earlier by ROOT’s import hook. Once a module is loaded, Python never re-resolves it.
+# environment: is already working. It just can’t solve this class of problem. PYTHONPATH is correct; ROOT imports mermithid first and locks it in — only an editable install can override that.
+#environment: #- PYTHONPATH=/mermithid:${PYTHONPATH}
+# Overriding sys for PYTHONPATH
 services:
   mermithid:
     image: mermithid:sensitivity_branch
@@ -12,9 +25,15 @@ services:
       - ~/repos/mermithid_share:/host_data
       - ~/repos/termite:/termite
     command: >
-      bash -c "
+      bash -lic "
         source /usr/local/p8/mermithid/v1.2.3/setup.sh &&
         source /usr/local/p8/mermithid/v1.2.3/bin/this_phylloxera.sh &&
         source /usr/local/p8/mermithid/v1.2.3/bin/this_cicada.sh &&
-        exec bash
-      "
+        pip install -e /mermithid &&
+        bash
+#bash -c "
+#source /usr/local/p8/mermithid/v1.2.3/setup.sh &&
+#source /usr/local/p8/mermithid/v1.2.3/bin/this_phylloxera.sh &&
+#source /usr/local/p8/mermithid/v1.2.3/bin/this_cicada.sh &&
+#exec bash
+#"

documentation/sensitivity_configurable_parameters.rst

@@ -7,6 +7,7 @@ The sensitivity calculation is configured using a configuration file. The config
 Our main goal is the calcualtion of sensitivity in a cavity experiment. The configurations below are to be used for the CavitySensitivity class in https://github.com/project8/mermithid/blob/feature/sensitivity_curve/mermithid/misc/SensitivityCavityFormulas.py
 This class is used by the CavitySensitivityCurveProcessor and the SensitivityParameterScanProcessor. 
 
+More information regarding parameter choices can be found in the Pre-CDR
 Structure of a config file
 --------------------------
 
@@ -15,6 +16,7 @@ Configuration files have several sections:
 
 * Experiment
 * Efficiency
+* Threshold
 * FrequencyExtraction
 * DopplerBroadening
 * MagneticField
@@ -39,14 +41,17 @@ Below is a list of the parameters with a short description of what role they pla
 * ``number_density``: This is the mean tritium gas density between the electron trapping coils. The density make be different at axial positions past the trap coils. The number_density together with the total volume, livetime, and the efficiency determines the statistical power of the experiment. Gas density also determines the track length and therefore the frequency resolution. The sensitivity curve processor can optimize this parameter to maximize the sensitivity. In that case this number is overwritten in the calculation. 
 * ``sri_factor``: The statistical rate increase factor articifially increases the number of observed events (it multiplies the total efficiency). It is highly recommended to set it to 1.
 * ``atomic``: If true, the calculation is done for atomic tritium. If false, moecular tritium is assumed. This affects the number of decays per gas molecule/atom (2 for molecular 1 for atomic), the track length in a given gas density (via electron scattering cross section), and the width of the final ground state.
-
+* ``active gas fraction``: The fraction of the gas within the CRES volume that's tritium. Very high (~1) in at atomic experiment. Lower in a molecular experiment due to presence of HT, H2, DT, and 3He, and possibly other gasses, too.
 
 **Efficiency**
 
-* ``usefixedvalue``: If true, fixed efficiency is used. If false, the efficiency is the product of radial, detection, and trapping efficienc. The trapping efficiency is calculated from the minimum pitch angle. 
+* ``usefixedvalue``: If true, fixed efficiency is used. If false, the efficiency is the product of radial, detection, and trapping efficiency. The trapping efficiency is calculated from the minimum pitch angle. Most of these are calculated in mermithid now. 
 * ``fixed_efficiency``: For example, set to roughly 2% for a 88deg minimum trapped pitch angle, assuming 100% detection efficiency of the trapped angles.
 * ``radial_efficiency``: Typically set to 0.67 from a calcualtion done for a 325MHz cavity with Halbach bite and radial cut on power of > 0.5 * maximum power.
 * ``detection_efficiency``: Fraction of events that is not detected.
+* ``unusable_dist_from_wall``: Cavity radius taken up by non-uniform field (The greater value between the Larmor radius and ioffe bite distance).
+
+**Threshold**
 
 **FrequencyExtraction**
 
@@ -75,10 +80,9 @@ We use the CRLB for calculating the frequency resolution. The CRLB is calculated
 * ``usefixedvalue``: If True ``default_systematic_smearing`` and ``default_systematic_uncertainty`` are used.
 * ``default_systematic_smearing``: Default systematic broadening for this category. Units must be eV.
 * ``default_systematic_uncertainty``: Default systematic uncertainty for this category. Units must be eV.
-* ``gas_temperature``: Temperature of the source gas. This should only be different from the cavity temperature if the gas is not in thermal equilibrium with the cavity. The gas temperature is used to calculate the Doppler broadening.
+* ``gas_temperature``: Temperature of the source gas. This should only be different from the cavity temperature if the gas is not in thermal equilibrium with the cavity. The gas temperature is used to calculate the Doppler broadening. In the molecular case, the molecules are in thermal equilibrium with the wall (85 K), but in the atomic case they are not, because they are prevented from contacting the wall by the Ioffe trap (4 mK - trapped gas temperature not wall temperature of 4 K)
 * ``gas_temperature_uncertainty``: Absolute uncertainty of the gas temperature.
-* ``fraction_uncertainty_on_doppler_broadening``: Fractional uncertainty on the Doppler broadening.
-
+* ``fraction_uncertainty_on_doppler_broadening``: Fractional uncertainty on the Doppler broadening. 
 
 **MagneticField**
 
@@ -93,6 +97,14 @@ We use the CRLB for calculating the frequency resolution. The CRLB is calculated
 * ``sigmae_theta``: Fixed input in eV. Energy broadening remaining after theta reconstruction, from electrons with lower pitch angles exploring high fields. Accounts for both the uncertainty on theta and uncertainties on the trap depth/boxiness.
 * ``sigmae_phi``: Fixed input in eV. Energy broadening from phi field inhomogeneity that remains after phi reconstruction.
 
+Each variable corresponds to the energy resolution contribution from field variation along a particular spatial direction / time, after accounting for position reconstruction in that direction / time. Then fraction_uncertainty_on_field_broadening accounts for how well that field broadening is known/calculated.
+Broken up between different position and time coordinates, for cases in which we want to study the effect of field variation along each coordinate individually, and then input the resulting resolution contributions into mermithid. Resolution contributions due to field variation along the different coordinates are added in quadrature.
+Currently using the variable sigmae_r to account for the total resolution contribution from field variation + reconstruction/calculation in all coordinates.
+
+In the atomic calculator, the input is a field inhomogeneity value deltaB/B in ppm, instead of being a resolution contribution value in eV. Field inhomogeneity is converted to a resolution value using a simple B-->E formula, based on the Larmor formula. Used to have something similar in mermithid but scrapped that.
+Don't love this approach because it implies that a certain physical field inhomogeneity corresponds to a certain resolution contribution via the Larmor formula and that's really not the case. Take, for example, radial variation. If there is no radial reconstruction, then to get the resolution contribution from the physical field variation, one needs to account for the fact that there are more electrons at higher radii, and that the detection efficiency depends on radius in some way. Those effects re-weight the impact of the field vs. radius profile. In addition, if there is radial reconstruction, that can further reduce the resolution relative to the physical radial field inhomogeneity, because we can use knowledge of electron radii to correct for the field vs. radius variation. So, there are several degrees of separation between physical field inhomogeneity and a resolution contribution.
+Since we don't have models of these various effects in mermithid or the atomic calculator, prefer to just have us directly input the resolution contributions. Yes, one could take the 0.085eV number and convert it to a deltaB/B requirement via the Larmor formula, but that number doesn't really mean anything.
+
 **FinalStates**
 
 * ``ground_state_width_uncertainty_fraction``: Uncertainty on the ground state width. Recommended to use 0.001.

mermithid/misc/Constants_numericalunits.py

@@ -1,35 +1,94 @@
+# C033 - Bayesian Paper (unitless) 
 '''
 Some constants useful for various things...
 The constants here use the numericalunits package. For constants not using this package import form Constants.py
 '''
 
 import numpy as np
 
-from numericalunits import e, me, c0, eps0, kB, hbar
-from numericalunits import meV, eV, keV, MeV, mm, cm, m, ns, s, Hz, kHz, MHz, GHz, amu, nJ
-from numericalunits import nT, uT, mT, T, mK, K,  C, F, g, W
-from numericalunits import hour, year, day, s, ms
-from numericalunits import mu0, NA, kB, hbar, me, c0, e, eps0, hPlanck
-
+from numericalunits import e, me, c0, eps0, mu0, kB, hbar, Rgas, NA, hPlanck    # Constants
+from numericalunits import meV, eV, keV, MeV, nJ, J, mm, cm, m                  # Energy, Length
+from numericalunits import nT, uT, mT, T, mK, K, F, W                           # Magnetic field, Temp, Power
+from numericalunits import hour, year, day, s, ms, ns, Hz, kHz, MHz, GHz        # Time, Frequency
+from numericalunits import kg, g, amu, mol                                      # Mass
+from numericalunits import C, A, mA, uA, V, mV, nV, kV                          # Charge, Amps, Voltage
+from numericalunits import Pa, bar, mbar, atm, torr, mtorr, L, mL               # Pressure, Volume
 
 T0 = -273.15*K
-
 tritium_livetime = 5.605e8*s
 tritium_mass_atomic = 3.016* amu *c0**2
-tritium_electron_crosssection_atomic = 9.e-23*m**2 #Hamish extrapolated to 18.6keV using Shah et al. (1987): https://iopscience.iop.org/article/10.1088/0022-3700/20/14/022
+tritium_electron_crosssection_atomic = 9.e-23*m**2 #Inelastic cross-section Hamish extrapolated to 18.6keV using Shah et al. (1987): https://iopscience.iop.org/article/10.1088/0022-3700/20/14/022. Full Bethe formula must be exactly right (down to magnetic and QED corrections) for the hydrogen atom.
+#tritium_electron_crosssection_atomic = 1.32e-22*m**2 #Inelastic cross-section + Elastic cross-section for T-e
+tritium_tritium_crosssection_atomic = 4.40e-16*m**2 # T-T cross-section at T=0
 tritium_endpoint_atomic = 18563.251*eV
+# C033 - [eV^-3] Bayesian Paper
 last_1ev_fraction_atomic = 2.067914e-13/eV**3
 
 tritium_mass_molecular = 6.032099 * amu *c0**2
 tritium_electron_crosssection_molecular = 3.67*1e-22*m**2 #[Inelastic from Aseev (2000) for T2] + [Elastic from Liu (1987) for H2, extrapolated by Elise to 18.6keV]
 tritium_endpoint_molecular = 18574.01*eV
+# C034 - [eV^-3] Bayesian Paper
 last_1ev_fraction_molecular = 1.67364e-13/eV**3
 
 ground_state_width = 0.436 * eV
 ground_state_width_uncertainty = 0.001*0.436*eV
 
 gyro_mag_ratio_proton = 42.577*MHz/T
 
+# Atomic Calculator
+# C008 - [m^3] Volume of 1mol of ideal gas at 1atm. Higher pressure --> smaller volume; higher temp --> lower volume
+molar_volume = 0.0224*m**3
+# C009 - [m^2] H-He Cross-section at low temp but > 5K (Berlinsky)
+H_He_crosssection = 2e-19*m**2
+# C010 - [m^2] H-He Cross-section at low temps (Berlinsky)
+H_He_crosssection_low_temp = 3e-20*m**2
+# C012 - [m^2] Hard Spheres cross-section (289 pm kinetic diameter
+H_H2_crosssection = 2.62e-19*m**2
+# C019 - [eV/amu] Conversion
+eV_amu = 931494100*eV/amu
+# C020 - [kg/amu] Conversion
+kg_amu = 1.66e-27*kg/amu
+# C021 - Boltzmann Constant [eV/K] = kB * (1 eV / 1.61e-19 J)
+kB_eV = 8.617e-5*eV/K
+# C026 - [J/T] Bohr Magneton
+bohr_magneton = 9.27e-24*J/T
+# C027 - [eV/T] Bohr Magneton
+#bohr_magneton = 5.776e-5*eV/T
+# C028 - Curie Becquere Conversion: 1 Ci = 3.7e10 Bq. NOTE CURIE IS NOT A UNIT IN NUMERICAL UNITS SO DO THIS EXPLICTLY
+Ci_Bq = 3.7e10/s
+# C030 - [s^-1] Atomic tritium decay constant
+lambda_tritium = 1 / tritium_livetime
+# C031 - Ground state branch (atomic)
+ground_state_branch_atomic = 0.702
+# C032 - Ground state branch (molecular)
+ground_state_branch_molecular = 0.570
+# C035 - [eV] Difference between atomic and effective molecular endpoints, effective difference 10-2 eV
+endpoint_diff = tritium_endpoint_molecular - tritium_endpoint_atomic
+# C036 - [eV] Recoil energy of tritium atom (Bodine)
+atomic_tritium_recoil_energy = 3.409*eV
+# C038 - [eV] Binding energy of tritium molecule (Bodine)
+molecular_tritium_binding_energy = 4.59*eV
+# C039 - [eV] Molecular final-state g.s. manifold standard deviation
+molecular_final_state_manifold = 0.436*eV
+# C040 - Constant for saturated T2 vapor: A,  Souers et al.
+T2_vapor_A = 5.84605
+# C041 - Constant for saturated T2 vapor: B
+T2_vapor_B = -160.7
+# C042 - Constant for saturated T2 vapor: B'
+T2_vapor_B_prime = 2.3235
+# C043 - [m/s^2] gravitational constant
+gravity = 9.80 *m/s**2
+# C048 - Multiplier for Lagendijk G^d rates: Ben Jones, Morgan Elliott CM presentation 10/24
+LGd_rates = 50
+# C049 - per beta decay: see: https://www.overleaf.com/2817746228snnghrnzfthk
+molecules_desorbed_wall_beta = 1000
+# C050 - [eV] to [J] conversion
+eV_J = 1.605e-19 *J/eV
+# C108 - [m^2] Choose an area to determine injection speed and density
+area_atom_loading_aperture = 0.001 * m**2
+# C118 - [Ci] Total activity on wall at recyling. Ci (Bq = 1 decay per second)
+wall_activity = 100
+
 # units that do not show up in numericalunits
 # missing pre-factors
 fW = W*1e-15
@@ -43,4 +102,4 @@
 
 # radian and degree which are also not really units
 rad = 1
-deg = np.pi/180
+deg = np.pi/180

mermithid/misc/FakeTritiumDataFunctions.py

@@ -341,7 +341,7 @@ def convolved_bkgd_rate_arrays(K, Kmin, Kmax, lineshape, ls_params, min_energy,
     return convolved
 
 
-
+# https://3.basecamp.com/3700981/buckets/3107037/documents/10099116348
 ##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'):

mermithid/processors/Sensitivity/CavitySensitivityCurveProcessor.py

@@ -71,7 +71,7 @@ def InternalConfigure(self, params):
         # options
         self.optimize_main_density = reader.read_param(params, 'optimize_main_density', True)
         self.optimize_comparison_density = reader.read_param(params, 'optimize_comparison_density', True)
-        self.verbose = reader.read_param(params, 'verbose', True)
+        self.verbose = reader.read_param(params, 'verbose', False)
         self.comparison_curve = reader.read_param(params, 'comparison_curve', False)
         self.B_error = reader.read_param(params, 'B_inhomogeneity', 7e-6)
         self.B_error_uncertainty = reader.read_param(params, 'B_inhom_uncertainty', 0.05)
@@ -239,9 +239,14 @@ def InternalRun(self):
         self.sens_main.BackgroundRate()
         logger.info('RF background: {}/eV/s'.format(self.sens_main.RF_background_rate_per_eV*eV*s))
         logger.info('Total background: {}/eV/s'.format(self.sens_main.background_rate*eV*s))
+        if self.sens_main.Efficiency.T2_background_atomic_trap:
+            logger.info("***T2 background in atomic trap:***")
+            self.sens_main.print_T2_background_atomic_trap()
+        if self.sens_main.Efficiency.pumping_calculation:
+            logger.info("***Pumping Requirements:***")
+            self.sens_main.print_pumping_requirements()
         logger.info("***Done printing pre-optimization***")
 
-
         #Optimizing the detection threshold for the comparison config files
         #Before the density optimization
         if self.comparison_curve:
@@ -590,7 +595,7 @@ def print_disclaimers(self):
                         scenarios, if the minimum allowed density is 1e-20 atoms/m^3, the optimization \
                         over density still works.")
         logger.info("Once you have read these disclaimers and are familiar with them, you can set \
-                    verbose==False in your config dictionary to stop seeing them.")
+                    verbose==False in your config dictionary (CavitySensitivtyCurveProcessor or test_anaylsis/) to stop seeing them.")
 
 
     def create_plot(self):

mermithid/processors/Sensitivity/SensitivityCurveProcessor.py

@@ -125,6 +125,11 @@ def InternalRun(self):
 
         self.create_plot()
 
+        if self.sens_main_is_atomic:
+            num_atoms_per_particle = 1
+        else:
+            num_atoms_per_particle = 2
+
         # add second and third x axis for track lengths
         if self.track_length_axis:
             self.add_track_length_axis()
@@ -170,15 +175,18 @@ def InternalRun(self):
 
 
         logger.info('Main curve (veff = {} cm**3, rho = {} /m**3):'.format(self.sens_main.Experiment.v_eff/(cm**3), rho_opt*(m**3)))
-        logger.info('CL90 limit: {}'.format(self.sens_main.CL90(Experiment={"number_density": rho_opt})/eV))
-        logger.info('T2 in Veff: {}'.format(rho_opt*self.sens_main.Experiment.v_eff))
-        logger.info('Total signal: {}'.format(rho_opt*self.sens_main.Experiment.v_eff*
+        #logger.info('CL90 limit: {}'.format(self.sens_main.CL90(Experiment={"number_density": rho_opt})/eV))
+        #logger.info('T2 in Veff: {}'.format(rho_opt*self.sens_main.Experiment.v_eff))
+        #logger.info('Total signal: {}'.format(rho_opt*self.sens_main.Experiment.v_eff*
+        logger.info('Neutrino mass 90% C.L. : {}eV'.format(self.sens_main.CL90(Experiment={"number_density": rho_opt})/eV))
+        logger.info('Number of tritium particles in Veff at optimal density: {}'.format(rho_opt*self.sens_main.Experiment.v_eff))
+        logger.info('Total number of signal events: {}'.format(rho_opt*self.sens_main.Experiment.v_eff*
                                                    self.sens_main.Experiment.LiveTime/
-                                                   self.sens_main.tau_tritium*2))
-        logger.info('Signal in last eV: {}'.format(self.sens_main.last_1ev_fraction*eV**3*
+                                                   self.sens_main.tau_tritium*num_atoms_per_particle))
+        logger.info('Total number of signal events in last eV: {}'.format(self.sens_main.last_1ev_fraction*eV**3*
                                                    rho_opt*self.sens_main.Experiment.v_eff*
                                                    self.sens_main.Experiment.LiveTime/
-                                                   self.sens_main.tau_tritium*2))
+                                                   self.sens_main.tau_tritium*num_atoms_per_particle))
 
         self.sens_main.print_statistics()
         self.sens_main.print_systematics()