diff --git a/README.md b/README.md index d4af8611..172d422c 100644 --- a/README.md +++ b/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 diff --git a/docker-compose.yaml b/docker-compose.yaml index cb13bac2..7390e5b4 100644 --- a/docker-compose.yaml +++ b/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 +#" diff --git a/documentation/sensitivity_configurable_parameters.rst b/documentation/sensitivity_configurable_parameters.rst index ab42affd..2e270bc7 100644 --- a/documentation/sensitivity_configurable_parameters.rst +++ b/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. diff --git a/mermithid/misc/Constants_numericalunits.py b/mermithid/misc/Constants_numericalunits.py index e1577c41..7e3d97c3 100644 --- a/mermithid/misc/Constants_numericalunits.py +++ b/mermithid/misc/Constants_numericalunits.py @@ -1,3 +1,4 @@ +# 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 @@ -5,24 +6,28 @@ 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 @@ -30,6 +35,60 @@ 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 \ No newline at end of file +deg = np.pi/180 diff --git a/mermithid/misc/FakeTritiumDataFunctions.py b/mermithid/misc/FakeTritiumDataFunctions.py index 4bd26fe2..cbaa5cdc 100644 --- a/mermithid/misc/FakeTritiumDataFunctions.py +++ b/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'): diff --git a/mermithid/processors/Sensitivity/CavitySensitivityCurveProcessor.py b/mermithid/processors/Sensitivity/CavitySensitivityCurveProcessor.py index fe76187a..12ad06dc 100644 --- a/mermithid/processors/Sensitivity/CavitySensitivityCurveProcessor.py +++ b/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): diff --git a/mermithid/processors/Sensitivity/SensitivityCurveProcessor.py b/mermithid/processors/Sensitivity/SensitivityCurveProcessor.py index ddb2f129..06de5da5 100644 --- a/mermithid/processors/Sensitivity/SensitivityCurveProcessor.py +++ b/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() diff --git a/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py b/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py index 7ba04d48..4c67e5d5 100644 --- a/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py +++ b/mermithid/processors/TritiumSpectrum/FakeDataGenerator.py @@ -68,7 +68,7 @@ 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', 25813125000.0) #Minimum frequency + self.minf = reader.read_param(params, 'minf', 25.8e+9) #Minimum frequency self.maxf = reader.read_param(params, 'maxf', None) if self.Kmax <= self.Kmin: logger.error("Kmax <= Kmin!") @@ -120,9 +120,11 @@ def InternalConfigure(self, params): # get file content if needed # get efficiency dictionary if self.apply_efficiency: + print("00000") self.efficiency_dict = self.load_efficiency_curve() np.random.seed() else: + print("11111") self.efficiency_dict = None # generate data with lineshape @@ -201,6 +203,7 @@ def InternalRun(self): else: ROIbound = [self.Kmin, self.Kmax] + print(self.efficiency_dict) Kgen = self.generate_unbinned_data(self.Q, self.m, ROIbound, self.S, self.B_1kev, diff --git a/mermithid/sensitivity/AtomicCalculator.py b/mermithid/sensitivity/AtomicCalculator.py new file mode 100644 index 00000000..e1f71e46 --- /dev/null +++ b/mermithid/sensitivity/AtomicCalculator.py @@ -0,0 +1,279 @@ +import numpy as np +from scipy.special import roots_laguerre +from mermithid.misc.Constants_numericalunits import * +''' +Ben Jones calculation disagrees with evaporation losses (specifically out the top > 5 mK, agrees otherwise) in atomic calc and must be integrated: https://3.basecamp.com/3700981/buckets/3107037/uploads/8902667883 + +''' + +# Jins functions - Atomic Calculator +# Python decides whether a name is local or global at compile time (not runtime). The name must already exist in the module’s global namespace by the time the function runs, or Python has nowhere to bind it. + +# C157 - [m] Fiducial distance from the cavity wall between max(ioffe_bite, larmor_radius) +# Ioffe field at wall is now the quadrature of the central field and Ioffe field by itself, everywhere it is used. +def calculate_ioffe_bite(nominal_field, magnetic_inhomogenity, ioffe_field, ioffe_multipolarity, cavity_radius): + return ((2 * nominal_field**2 * magnetic_inhomogenity / (2 * nominal_field**2 * magnetic_inhomogenity + ioffe_field**2))**(1/(ioffe_multipolarity-2)) * cavity_radius * -1) + cavity_radius + +# C271 - [m/s] Atomic thermal speed of 3T or 3He (C270) +def calculate_average_velocity(temp): + average_velocity = c0 * np.sqrt((8 * 0.025*eV * temp / np.absolute(T0)) / (np.pi * tritium_mass_atomic)) + #average_velocity = c0 * np.sqrt((8 * kB * temp) / (np.pi * tritium_mass_atomic * eV_J)) + return average_velocity + +# C334 - [m^2] Cylinder wall; no endcaps included (2*pi*R*L) +def calculate_trap_wall_area(cavity_radius, cavity_length): + # (Cavity length is cavity top plate - top of cone / 2) + trap_wall_area = 2 * np.pi * cavity_radius * cavity_length + return trap_wall_area + +# C201 - [m^-1] Add the gravity-temperature scale parameter “b” for magnetogravitational trap. For pure magnetic trap, make gravity weaker by ratio of L/D. +def calculate_gravity_temperature_scale(trapped_gas_temp, pure_magnetic_flag, cavity_L_over_D): + gravity_temperature_scale = ((tritium_mass_atomic/c0**2 * gravity / (kB * trapped_gas_temp)) * np.absolute(pure_magnetic_flag * (1 + 1 / cavity_L_over_D) - 1)) + return gravity_temperature_scale + +# C202 - [m^-2] Add the surface density scale parameter “a/b” for magnetogravitational trap +def calculate_surface_density_scale(design_density, gravity_temperature_scale, top_plate_cavity): + surface_density_scale = design_density * top_plate_cavity / (1 - np.exp(-gravity_temperature_scale * top_plate_cavity)) + return surface_density_scale + +# C204 - [atom] Inventory in physical volume +def calculate_inventory(design_density, volume): + return design_density * volume + +# C203 - [m^-3] Mean Trap Density within coils. Density in the electron-trapping region as a function of the density in the full volume. Differ depending on whether the experiment is horizontal or vertical. (mean density) +# Full volume (everywhere where atoms are trapped) includes some length of the cavity past the ends of the trap, as well as the volume inside the Ioffe cones. +# In a horizontal experiment, the average atom density in the electron trapping region would be about the same as in the full volume, because the electron trapping region and the non-electron trapping region have the same range/distribution of heights in the them. +# In a vertical experiment, the atom density will be different in the full volume compared with the electron-trapping volume, since atoms pool toward the bottom Ioffe cone. +def calculate_mean_trap_density(design_density, top_plate_cavity, trapped_gas_temp, pure_magnetic_flag, cavity_L_over_D, trap_coil_1, trap_coil_2): + gravity_temp_scale = calculate_gravity_temperature_scale(trapped_gas_temp, pure_magnetic_flag, cavity_L_over_D) + surface_density_scale = calculate_surface_density_scale(design_density, gravity_temp_scale, top_plate_cavity) + mean_trap_density = design_density * top_plate_cavity * (np.exp(-gravity_temp_scale * trap_coil_1) - np.exp(-gravity_temp_scale * trap_coil_2)) / ((1 - np.exp(-gravity_temp_scale * top_plate_cavity)) * (trap_coil_2 - trap_coil_1)) + return mean_trap_density + +# C208 - [decays/s] Radioactivity in trap per cavity +def calculate_activity_in_trap(number_density, cavity_radius, trap_coil_1, trap_coil_2): + # C206 - [atom] Use exponential density vertical gradient, calculate mean density between trap coils. Calculate atoms and activity between the trap coils. + atoms_between_trap_coils = number_density * np.pi * cavity_radius**2 * np.absolute(trap_coil_2 - trap_coil_1) + return atoms_between_trap_coils / tritium_livetime + +# Background T2 in atomic trap: +def calculate_T2_background_atomic_trap(cavity_radius, cavity_length, cavity_wall_temp, max_ratio_T2_T, number_density): + # C306 - [mbar] Vapor Pressure of T2 w/ constants for saturated T2 vapor from Souers et al. + T2_vapor_pressure = mbar * np.exp(T2_vapor_A + (T2_vapor_B*K/cavity_wall_temp) + T2_vapor_B_prime * np.log(cavity_wall_temp/K)) / 0.76 + # C307 - [m^-3] Density of saturated vapor + sat_vapor_density = NA * (T2_vapor_pressure/mbar) * np.absolute(T0/cavity_wall_temp) / (1000 * molar_volume) + # C308 - [m^-3] = s^-1 * sqrt(kg/eV) / m^2 T2 density form desorption at end of a cycle + T2_density_desorp = 4 * molecules_desorbed_wall_beta * (wall_activity * Ci_Bq) * np.sqrt(molecules_desorbed_wall_beta * (tritium_mass_atomic / c0**2) / (2 * atomic_tritium_recoil_energy * eV_J )) / calculate_trap_wall_area(cavity_radius, cavity_length) + # C309 [m^-3] Total T2 density + T2_total_density = sat_vapor_density + T2_density_desorp + # C310 - T2/T number ratio. Activity ratio is 1.64 times bigger + T2_T_ratio = T2_total_density / number_density + #return T2_vapor_pressure, sat_vapor_density, T2_density_desorp, T2_total_density, T2_T_ratio + return T2_total_density, T2_T_ratio + +# Aperture Heat Leak: +def calculate_aperture_heat_leak(trapped_gas_temp, design_density, volume): + # C250 - [m/s] Average velocity of trapped gas + trapped_gas_velocity = calculate_average_velocity(trapped_gas_temp) + # C251 - [atom/s] Inward atom current required through aperture to balance losses (j_aperture) + current_aperture_leak = area_atom_loading_aperture * design_density * trapped_gas_velocity / 4 + # C252 - [s] Time constant for loss through aperture + time_constant_aperture = calculate_inventory(design_density, volume) / current_aperture_leak + return time_constant_aperture, current_aperture_leak + +# Radioactivity Heat Leak: +# Radioactivity loss is just a single number from Ben Clark’s thesis, needs to be done better. Dependent on a cross-section that varies with temp. +def calculate_rad_heat_leak(cavity_radius, number_density, trap_coil_1, trap_coil_2, design_density, volume, net_efficiency): + # C256 - [s] Input from Ben Clark's thesis with a mirror ratio of 0.5 + half_life = 11 * 24 * 3600 * s + trap_activity = calculate_activity_in_trap(number_density, cavity_radius, trap_coil_1, trap_coil_2) + # C257 - [atom/s] Atom current required to balance loss due to radioactivity + current_rad_leak = trap_activity * np.log(2) * tritium_livetime / (net_efficiency * half_life) + # C258 - [s] Time constant for loss due to radioactivty (tau_rad) + time_constant_rad = calculate_inventory(design_density, volume) / current_rad_leak + return time_constant_rad, current_rad_leak + +# T2 Desorption from the wall: T2 desorption from walls as a background and a trap-heating loss. +# Burst of molecules is emitted from the wall with each decay and can knock out atoms from the trap. Description in CDR 4.5.5. New parameter to enter is the number of molecules (choose 1000 for now). +def calculate_T2_desorption_from_wall(cavity_radius, cavity_length, design_density, volume): + # C261 - [s] Mean lifetime of atom in trap from desorption + time_constant_desorp = calculate_trap_wall_area(cavity_radius, cavity_length) / (2 * molecules_desorbed_wall_beta * (wall_activity * Ci_Bq) * H_H2_crosssection) + # C262 - [atom/s] Atom current required to keep up with desorption losses + current_desorp = calculate_inventory(design_density, volume) / time_constant_desorp + return time_constant_desorp, current_desorp + +# He Heat Leak: Notable difference comes from pumping speed, physical volume, and design density +def calculate_He_heat_leak(pumping_speed_theoretical, pumping_cavity_termination, turbopump_speed, cavity_temperature, design_density, volume): + He_velocity = calculate_average_velocity(cavity_temperature) + # C266 - [atom/s] Not including flow for He heat leak itself + He_production_rate = wall_activity * Ci_Bq + # C269 - [m^-3] Number density of Helium + He_density = He_production_rate * ((1/pumping_speed_theoretical) + (1/pumping_cavity_termination) + (1/turbopump_speed)) + # C272 - [s] Time constant for loss due to He-3 heat leak + time_constant_He = 1 / (He_density * He_velocity * H_He_crosssection) + # C273 - Fraction of total gas that is He + He_fraction = He_density / design_density + # C274 - [s^-1] Atom current required to keep up with He-3 heat leak + current_He_leak = calculate_inventory(design_density, volume) / time_constant_He + return time_constant_He, current_He_leak + +# Dipolar loss rate: Calculation with z-dependent density, cylinder & cone. Polynomial fit of nominal field. +def calculate_dipolar_loss(nominal_field, cavity_radius, design_density, volume, trapped_gas_temp, cavity_L_over_D, top_cone, top_plate_cavity, pure_magnetic_flag): + gravity_temperature_scale = calculate_gravity_temperature_scale(trapped_gas_temp, pure_magnetic_flag, cavity_L_over_D) + surface_density_scale = calculate_surface_density_scale(design_density, gravity_temperature_scale, top_plate_cavity) + # C298 - [m^3/s] Polynomial fit for Dipolar spin-flip rate (G_dd). Depends on field (Lagendijk et al). Polynomial-log fit used now. + dipolar_spin_flip_rate = (60.106 + 13.812 * np.log(nominal_field) - 4.7867 * np.log(nominal_field)**2 - 2.3192 * np.log(nominal_field)**3 \ + - 0.32663 * np.log(nominal_field)**4 - 0.015775 * np.log(nominal_field)**5) * 1e-22 * LGd_rates * m**3/s + # C299 - [atom/s] Flow to keep up with dipolar losses + current_dipolar = dipolar_spin_flip_rate * np.pi * cavity_radius**2 * surface_density_scale**2 * (gravity_temperature_scale * (np.exp(-gravity_temperature_scale * top_cone) \ + - np.exp(-gravity_temperature_scale * top_plate_cavity)) / 2 + (1 - np.exp(-2 * gravity_temperature_scale * top_cone) * (2 * gravity_temperature_scale**2 \ + * top_cone**2 + 2 * gravity_temperature_scale * top_cone + 1)) / (4 * top_cone**2 * gravity_temperature_scale)) + # C300 - [s] Time constant for loss due to dipolar spin-flip loss + time_constant_dipolar = calculate_inventory(design_density, volume) / current_dipolar + return time_constant_dipolar, current_dipolar + +# Evaporation loss rate: Does not take into account density of states with height. Magnetic potential limits evaporation. +# Cone and cylinder now separate because cone is weaker owing to azimuthal modulation of Ioffe field. Each now has its own density multiplier. The weaker cone field is handled in a separate Igor calculation Coneangle.pxp outside this SS and entered as a loss rate multiplier in C111. +def calculate_evaporation_loss(cavity_L_over_D, pure_magnetic_flag, trapped_gas_temp, cavity_radius, ioffe_field, nominal_field, cavity_length, top_cone, top_plate_cavity, design_density, relative_loss_rate_cone_wall, volume): + gravity_temperature_scale = calculate_gravity_temperature_scale(trapped_gas_temp, pure_magnetic_flag, cavity_L_over_D) + surface_density_scale = calculate_surface_density_scale(design_density, gravity_temperature_scale, top_plate_cavity) + + # C283 - [m] Mean free path at the base of cavity + mfp_cavity_base = 1 / (gravity_temperature_scale * surface_density_scale * tritium_tritium_crosssection_atomic) + # C284 - Pure Magnetic eta + eta = (bohr_magneton / eV_J) * (np.sqrt(ioffe_field**2 + nominal_field**2) - nominal_field) / (trapped_gas_temp * 0.025*eV / np.absolute(T0)) + # C285 - Magnetogravitational eta + eta_grav = (tritium_mass_atomic/c0**2) * kg_amu * gravity * cavity_length / (kB * trapped_gas_temp) + # C286 - [s^-1] Evaporation out the top; set to 1/1000 for pure magnetic trap + evaporation_top = (np.log(1 + np.log(2) * np.exp(gravity_temperature_scale * top_plate_cavity) / (surface_density_scale * tritium_tritium_crosssection_atomic)) \ + * surface_density_scale * gravity_temperature_scale * np.exp(-gravity_temperature_scale * top_plate_cavity) * np.sqrt(np.pi * kB * trapped_gas_temp / (2 * tritium_mass_atomic / c0**2)) * cavity_radius**2) \ + * np.absolute(pure_magnetic_flag * (1 + 0.001) -1) + # C287 - [s] Time constant for evaporation out of the top + time_constant_evap_top = calculate_inventory(design_density, volume) / evaporation_top + # C288 - [s^-1] Evaporation to the sides + evaporation_sides = (2 * np.pi * cavity_radius * cavity_length * design_density * np.sqrt(kB * trapped_gas_temp / (2 * np.pi * tritium_mass_atomic / c0**2)) \ + * np.exp(-eta) * top_plate_cavity * (np.exp(-gravity_temperature_scale * top_cone) - np.exp(-gravity_temperature_scale * top_plate_cavity)) \ + / ((1 - np.exp(-gravity_temperature_scale * top_plate_cavity)) * (top_plate_cavity - top_cone))) + # C290 - [s^-1] Evaporation in cone; for pure magnetic, double it to account for both ends + evaporation_cone = (relative_loss_rate_cone_wall * np.pi * cavity_radius * np.sqrt(cavity_radius**2 + top_cone**2) * design_density \ + * np.sqrt(kB * trapped_gas_temp / (2 * np.pi * tritium_mass_atomic/c0**2)) * np.exp(-eta) * top_plate_cavity * (1 - np.exp(-gravity_temperature_scale * top_cone)) \ + / ((1 - np.exp(-gravity_temperature_scale * top_plate_cavity)) * top_cone)) * np.absolute(pure_magnetic_flag + 1) + # C292 - [s] Time constant for evaporation loss of cone and sides + time_constant_evap_cone_sides = calculate_inventory(design_density, volume) / (evaporation_cone + evaporation_sides) + # C295 - [atom/s] Flow to keep up with evaporation losses aka Total evaporation_sum + current_evaporation = evaporation_top + evaporation_sides + evaporation_cone + # C294 - [s] Time constant for total evaporation loss + time_constant_total_evaporation = calculate_inventory(design_density, volume) / current_evaporation + return time_constant_total_evaporation, current_evaporation, evaporation_top, evaporation_sides, evaporation_cone + +# T2 Heat Leak: +def calculate_T2_heat_leak(aperture_current, rad_current, He_current, evap_current, dipolar_current, T2_total_density, cavity_wall_temp, design_density, volume): + # C277 - [atom/s] T2 production rate (all atoms entering, except for T2 heat leak itself) + T2_production_rate = aperture_current + rad_current + He_current + evap_current + dipolar_current + # C279 - [s] Time constant for loss due to T2 heat leak + time_constant_T2 = 1 / (T2_total_density * H_H2_crosssection * calculate_average_velocity(cavity_wall_temp) / np.sqrt(2)) + # C280 - [atom/s] + current_T2_leak = calculate_inventory(design_density, volume) / time_constant_T2 + return time_constant_T2, current_T2_leak + +# C303 - [s] Total lifetime of trap: +def calculate_trap_lifetime(time_rad, time_desorp, time_He, time_T2, time_evap, time_dipolar): + return 1 / ((1/time_rad) + (1/time_desorp) + (1/time_He) + (1/time_T2) + (1/time_evap) + (1/time_dipolar)) + +# Atom Supply into Trap: +# The inputted atom current determines the (average) density in the full volume, but the density for mass sensitivity is average atom density within the electron trap. Conversion is needed between these two densities +def calculate_trap_atom_supply(rad_current, evap_current, dipolar_current, He_current, T2_current, design_density, volume): + # The ‘atomic current required’ includes the c-state atoms as well as d-state, although those are lost almost immediately to spin exchange. + # C107 - C state flag for hyperfine states for total gas into trap. Trap is always d-state only. No c-states (1), include c-states (2) + c_states_flag = 1 + # C318 - [atom/s] Total current is atom current (d state only) + He heat leak + T2 heat leak + total_atom_current = rad_current + evap_current + dipolar_current + He_current + T2_current + # C319 - [s] Total time constant for atoms to remain in the trap + total_time_constant_trap = calculate_inventory(design_density, volume) / total_atom_current + # C320 - [atom/s] Atom current with c-states and d-states in calculation + total_atom_current_states = total_atom_current * c_states_flag + return total_time_constant_trap, total_atom_current_states +''' +There are two possible methods for pumping away the tritium, all of which is eventually in molecular form. One is to keep the tritium in circulation by using turbopumps and avoiding temperatures below 10K in the trap region, and the second is to cool parts of the trap region outside the magnetic wall to < 3K in order to cryopump the tritium. +The mechanical pumping method is impractical because of the pumping speed required, a fraction of a billion L/s, where an achievable upper limit is 4 orders of magnitude smaller. +The speed requirement is driven by the combination of the input atomic current and the need to maintain the molecular fraction below 10−4. +The cryopumping method ties up very large amounts of tritium, tens to hundreds of kCi, by the end of a day, and the need to warm up and recycle that tritium on such a short time scale would lead to low statistical precision and instabilities. +Still need turbos to handle He-3 and to do pumpout during recycling, but initial hope that we could run turbos only and keep the T2 pressure low enough was not realized. Must recycle. +''' +# Turbopump Calculations: Need to convert tritium mass amu to eV : 1 amu ~ 931 MeV/c^2 and kB = 8.6E-05 eV/K. Multiply by sqrt(2) if atomic; # Useful for Atomic and Helium-3 +# C233/C236 - [m^3/s] (molecular/atomic) theoretical pumping speed. The mean speed of the cylinder connected to a perfect pump (Dushman). Uniform source, outputs mean density in cavity. Cannot exceed (obstruction not included). +def turbopump_speed_limit(cavity_radius, cavity_temperature, cavity_L_over_D, atomic_flag): + return (np.pi * cavity_radius**2 * c0 * np.sqrt(kB_eV * cavity_temperature / (4 * np.pi * tritium_mass_atomic)) / (0.5 + cavity_L_over_D / 8)) * (1 + atomic_flag * (np.sqrt(2) - 1)) +# C235/C238 - [m^3/s] (molecular/atomic) Assumed ambient room air (28-29 amu) temperature (293 K, and a pumping speed of 0.5 L/s +def cavity_termination_speed(turbopumping_speed_air, cavity_top_plate_temp, atomic_flag): + return (turbopumping_speed_air * np.sqrt(cavity_top_plate_temp * 28 * amu * c0**2 / (293 * K * 2 * tritium_mass_atomic))) * (1 + atomic_flag * (np.sqrt(2) - 1)) +# C239 - [m^3/s] Turbopump in series speed (2.5 m^3/s * 2) +def turbopump_speed(number_turbos, turbopumping_speed_gas): + return number_turbos * turbopumping_speed_gas +# C327 - [m^-3] Molecular density allowed by molecular/atomic assuming total atom density in all of physical volume +def molecular_density_allowed(design_density, max_ratio_T2_T): + return design_density * max_ratio_T2_T * ground_state_branch_atomic / 2 +# The Pumping speed required (molecular) in C328 is kind of a dimensional calculation, not really a speed. As a result, the cryopumping speed limit is probably not real. Ideal cryopumping is the impingement rate: it all sticks. The actual molecular density is determined by the vapor pressure (Souers formula). +# C328 - [m^3/s] - Pumping speed required to keep molecular pressure from exceeding reference limit +def turbopump_speed_required(atom_current, molecular_density_limit): + return atom_current / 2 / molecular_density_limit +# C331 - If above 1000, critical +def ratio_turbopump_speed(pumping_speed_required, pumping_speed_limit, cavity_termination_speed, turbopump_speed): + return pumping_speed_required * ((1 / pumping_speed_limit) + (1 / cavity_termination_speed) + (1 / turbopump_speed)) + +# Cryopumping Calculations: +def calculate_cryopump_speed(cavity_wall_temp, cavity_radius, cavity_length, pumping_speed_required, atom_current, design_density, max_ratio_T2_T): + # C335 - [m/s] Cryopumping speed per area for mass-3 atoms at wall temperature stated above. sqrt(k T/(2 pi M)); see https://www.synsysco.com/wp-content/uploads/2016/02/Basics-of-Cryopumping-Booklet.pdf. Using M=3 because most gas is spin-flipped atoms. + cryopump_speed_per_area = np.sqrt(kB_eV * cavity_wall_temp / (2 * np.pi * tritium_mass_atomic)) * c0 + # C336 - [m^3/s] With given cylindrical dimensions at wall temperature stated above for mass-3 atoms + trap_interior_surface_speed = cryopump_speed_per_area * calculate_trap_wall_area(cavity_radius, cavity_length) + # C338 - Ratio of required molecular speed to cryopumping speed. Ideal cryopumping is the impingement rate (all sticks). + ratio_cryopump_speed = pumping_speed_required / trap_interior_surface_speed + # C339 - [W] Heat delivered to the cold surface + recombination_heat_load = atom_current * molecular_tritium_binding_energy * eV_J / 2 + # C340 - [Ci/day] Accumulated activity per day + accumulated_activity_day = atom_current * lambda_tritium * (86400*s/day) / Ci_Bq + # C342 - [day] Time between recycling cryosurface. For cryopumping atomic experiment; determines inventory on walls + time_recycling_cryosurface = wall_activity / accumulated_activity_day + # C344 - [m^3/s] Cryopumping speed for molecules + cryopump_speed = np.sqrt(kB_eV * cavity_wall_temp / (4 * np.pi * tritium_mass_atomic)) * c0 * calculate_trap_wall_area(cavity_radius, cavity_length) + # C345 - [molecules/s] max allowed injection rate of molecules + max_molecules_injection = cryopump_speed * molecular_density_allowed(design_density, max_ratio_T2_T) + return ratio_cryopump_speed, cryopump_speed, time_recycling_cryosurface + +# Injection Line Calculations: Polynomial fit of injection field +# Density in the injection line, mfp in injection line, dipolar loss in injection line. There are 2 new entries for this in the choices – beam temperature and beamline field. +def calculate_injection_line(atom_current, cavity_radius, design_density, trapped_gas_temp, nominal_field, injection_gas_temp, injection_field): + # C349 - [m/s] Vertical flow speed in trap at bottom and in the beamline + vertical_injection_speed = atom_current / (np.pi * cavity_radius**2 * design_density) + # C350 - [J/atom] Energy in trapped gas and in beam gas + trapped_gas_energy = (1.5 * kB * trapped_gas_temp) + (0.5 * (tritium_mass_atomic/c0**2) * vertical_injection_speed**2) + (bohr_magneton * nominal_field) + # C351 - [m^-3] Density in the injection line + injection_density = atom_current * np.sqrt((tritium_mass_atomic/c0**2) / (2 * trapped_gas_energy - 3 * kB * injection_gas_temp - 2 * bohr_magneton * injection_field)) / area_atom_loading_aperture + # C352 - [m] Mean free path in beamline + mfp_beamline = 1 / (injection_density * tritium_tritium_crosssection_atomic) + # C353 - [m/s] Vertical flow speed through aperture + vertical_aperture_speed = design_density * np.pi * cavity_radius**2 * vertical_injection_speed / (area_atom_loading_aperture * injection_density) + # C354 - [J/atom] Energy in injected gas + injection_gas_energy = (1.5 * kB * injection_gas_temp) + (0.5 * (tritium_mass_atomic/c0**2) * vertical_aperture_speed**2) + (bohr_magneton * injection_field) + # C355 - [m^3/s] 5th degree polynomial fit for injection line dipolar loss rate + injection_dipolar_loss_rate = (60.106 + 13.812 * np.log(injection_field/T) - 4.7867 * np.log(injection_field/T)**2 - 2.3192 * np.log(injection_field/T)**3 \ + - 0.32663 * np.log(injection_field/T)**4 - 0.015775 * np.log(injection_field/T)**5) * 1e-22 * LGd_rates * m**3/s + # C356 - [m^-1] Beamline dipolar loss + beamline_dipolar = injection_dipolar_loss_rate * injection_density / vertical_aperture_speed + return vertical_injection_speed, trapped_gas_energy, injection_density, vertical_aperture_speed, injection_gas_energy + + +def calculate_activity_last_1eV_spectrum(self, atomic_flag, number_density, cavity_radius, trap_coil_1, trap_coil_2, total_efficiency): + trap_activity = calculate_activity_in_trap(number_density, cavity_radius, trap_coil_1, trap_coil_2) + if atomic_flag: + # C209 - [decay/s] Activity in last 100 eV of spectrum * net efficiency + # Note last_1eV_fraction_atomic is the fraction of events in the last 1eV + activity_last_100eV_efficiency = trap_activity * last_1ev_fraction_atomic * total_efficiency * 1000000 / ground_state_branch_atomic + # C211 - [decay/s] + last_1eV_atomic = activity_last_100eV_efficiency * ground_state_branch_atomic / 1000000 + return last_1eV_atomic + else: + activity_last_100eV_efficiency = trap_activity * last_1ev_fraction_molecular * total_efficiency * 1000000 / ground_state_branch_molecular + last_1eV_molecular = activity_last_100eV_efficiency * ground_state_branch_molecular / 1000000 + return last_1eV_molecular diff --git a/mermithid/sensitivity/SensitivityCavityFormulas.py b/mermithid/sensitivity/SensitivityCavityFormulas.py index d2e461e2..56760faa 100755 --- a/mermithid/sensitivity/SensitivityCavityFormulas.py +++ b/mermithid/sensitivity/SensitivityCavityFormulas.py @@ -1,6 +1,6 @@ ''' Class calculating neutrino mass sensitivities based on analytic formulas from CDR. -Author: R. Reimann, C. Claessens, T. E. Weiss, W. Van De Pontseele +Author: R. Reimann, C. Claessens, T. E. Weiss, W. Van De Pontseele, M. Oueslati Date: 06/07/2023 Updated: December 2024 @@ -15,18 +15,16 @@ from mermithid.misc.Constants_numericalunits import * from mermithid.misc.CRESFunctions_numericalunits import * from mermithid.cavity.HannekeFunctions import * +from mermithid.sensitivity.AtomicCalculator import * from mermithid.sensitivity.SensitivityFormulas import * - try: from morpho.utilities import morphologging logger = morphologging.getLogger(__name__) except: print("Run without morpho!") - - # Wouters functinos def db_to_pwr_ratio(q_db): return 10**(q_db/10) @@ -260,16 +258,30 @@ def CalcDefaults(self, overwrite=False): self.Jprime_0 = 3.8317 self.cavity_freq = frequency(self.T_endpoint, self.MagneticField.nominal_field) self.CavityRadius() - - #Get trap length from cavity length if not specified - if ((not hasattr(self.Experiment, 'trap_length')) or overwrite): - self.Experiment.trap_length = 0.8 * 2 * self.cavity_radius * self.Experiment.cavity_L_over_D - logger.info("Calc'd trap length: {} m".format(round(self.Experiment.trap_length/m, 3), 2)) + self.trap_geometry() + self.CavityLength() self.Efficiency = NameSpace({opt: eval(self.cfg.get('Efficiency', opt)) for opt in self.cfg.options('Efficiency')}) self.CavityVolume() self.CavityPower() + #Get trap length from cavity length if not specified + if ((not hasattr(self.Experiment, 'trap_length')) or overwrite): + self.Experiment.trap_length = 0.8 * 2 * self.cavity_radius * self.Experiment.cavity_L_over_D + logger.info("Calc'd trap length: {:.3} m".format(self.Experiment.trap_length/m, 3)) + elif self.Experiment.trap_length_calc_flag: + # C129 - Trap L/D = Trap L / 2 / cavity radius + self.Experiment.trap_length = self.trap_coil_2 - self.trap_coil_1 + logger.info("Calc'd trap length from coils: {} m".format(round(self.Experiment.trap_length/m, 3), 2)) + + # C128 - Calculate Cavity L/D Ratio with actual dimensions: top of cone, first trap coil, second trap coil, top plate of cavity, top of vacuum containment. + if ((hasattr(self.Experiment, 'cavity_L_over_D')) and self.Experiment.cavity_L_over_D_calc_flag): + self.Experiment.cavity_L_over_D = self.cavity_length / (self.cavity_radius * 2) + + # C203 - Calculate mean density between trap coils if design density is specified + if ((hasattr(self.Experiment, 'design_density')) and self.Experiment.design_density_flag): + self.Experiment.number_density = calculate_mean_trap_density(self.Experiment.design_density, self.top_plate_cavity, self.DopplerBroadening.gas_temperature, self.Experiment.pure_magnetic_flag, self.Experiment.cavity_L_over_D, self.trap_coil_1, self.trap_coil_2) + #Calculate position dependent trapping efficiency self.pos_dependent_trapping_efficiency = trapping_efficiency( z_range = self.Experiment.trap_length /2, bg_magnetic_field = self.MagneticField.nominal_field, @@ -284,6 +296,13 @@ def CalcDefaults(self, overwrite=False): #Cyclotron radius is sometimes used in the effective volume calculation self.cyc_rad = cyclotron_radius(self.cavity_freq, self.T_endpoint) + #Ioffe bite used in radial efficiency and effective volume calculation + self.unusable_dist_from_wall = self.Efficiency.unusable_dist_from_wall + if self.Efficiency.calculate_ioffe_bite_flag and all(hasattr(self.MagneticField, attr) for attr in ("nominal_field", "magnetic_inhomogenity", "ioffe_field", "ioffe_multipolarity")): + self.unusable_dist_from_wall = calculate_ioffe_bite(self.MagneticField.nominal_field, self.MagneticField.magnetic_inhomogenity, self.MagneticField.ioffe_field, self.MagneticField.ioffe_multipolarity, self.cavity_radius) + elif self.Efficiency.calculate_ioffe_bite_flag: + logger.info("Error: Did not specify all attributes for Ioffe bite calculation in config file") + #Assigning the background constant if it's not in the config file if hasattr(self.Experiment, "bkgd_constant"): self.bkgd_constant = self.Experiment.bkgd_constant @@ -341,9 +360,12 @@ def CalcDefaults(self, overwrite=False): self.EffectiveVolume() logger.info("Trap radius: {} cm".format(round(self.cavity_radius/cm, 3), 2)) logger.info("Total trap volume: {} m^3".format(self.total_trap_volume/m**3)) - logger.info("Cyclotron radius: {}m".format(self.cyc_rad/m)) + logger.info("Ioffe bite: {} m".format(self.unusable_dist_from_wall/m)) + logger.info("Cyclotron radius: {} m".format(self.cyc_rad/m)) if self.use_cyc_rad: - logger.info("Using cyclotron radius as unusable distance from wall, for radial efficiency calculation") + logger.info("Using cyclotron (Larmor) radius as unusable distance from wall, for radial efficiency calculation") + else: + logger.info("Using ioffe bite as unusable distance from wall, for radial efficiency calculation") #### #Initialization related to the energy resolution: @@ -376,20 +398,43 @@ def CavityRadius(self): axial_mode_index = 1 self.cavity_radius = c0/(2*np.pi*self.cavity_freq)*np.sqrt(self.Jprime_0**2+axial_mode_index**2*np.pi**2/(4*self.Experiment.cavity_L_over_D**2)) return self.cavity_radius - + + # Geometry of trap + def trap_geometry(self): + # C089 - [m] origin is the virtual point of the cone height + self.top_cone = 0.60*m + # C090 - [m] height of first trap coil + self.trap_coil_1 = 0.75*m + # C091 - [m] height of second trap coil + self.trap_coil_2 = 4.80*m + # C092 - [m] height of top plate of the cavity + self.top_plate_cavity = 5.30*m + # C093 - [m] height of the top of the vacuum system + self.top_vacuum_system = 7.50*m + return self.top_cone, self.trap_coil_1, self.trap_coil_2, self.top_plate_cavity, self.top_vacuum_system + # C126 - [m] Cavity Length from true L and f (estimated as L-z1/2); effective cavity length enters first third of ioffe cone + def CavityLength(self): + self.cavity_length = (self.top_plate_cavity - self.top_cone/2) + return self.cavity_length + def CavityVolume(self): - #radius = 0.5*wavelength(self.T_endpoint, self.MagneticField.nominal_field) - self.total_cavity_volume = 2*self.cavity_radius*self.Experiment.cavity_L_over_D*np.pi*(self.cavity_radius)**2*self.Experiment.n_cavities - + # Calculate vacuum volume including cone and top service volume - important for vertical LFA due to higher density at the bottom + if self.Experiment.cavity_cone_flag: + # V = V_cyl + V_cone = pi * r^2 * (L + h/3), total height is height of trap coils + height of cone; + # See Robertson_H&V_2025-10-29 for explanation; nothing changes about physical volume for second cone in horizontal configuration. + self.total_cavity_volume = np.pi * self.cavity_radius**2 * (self.top_plate_cavity - (2/3) * self.top_cone) * self.Experiment.n_cavities + else: + #radius = 0.5*wavelength(self.T_endpoint, self.MagneticField.nominal_field) + self.total_cavity_volume = 2*self.cavity_radius*self.Experiment.cavity_L_over_D*np.pi*(self.cavity_radius)**2*self.Experiment.n_cavities logger.info("Frequency: {} MHz".format(round(self.cavity_freq/MHz, 3))) logger.info("Wavelength: {} cm".format(round(wavelength(self.T_endpoint, self.MagneticField.nominal_field)/cm, 3))) logger.info("Cavity radius: {} cm".format(round(self.cavity_radius/cm, 3))) - logger.info("Cavity length: {} cm".format(round(2*self.cavity_radius*self.Experiment.cavity_L_over_D/cm, 3))) + logger.info("Cavity length: {} cm".format(round(self.cavity_length/cm, 3))) + #logger.info("Cavity length: {} cm".format(round(2*self.cavity_radius*self.Experiment.cavity_L_over_D/cm, 3))) logger.info("Total cavity volume: {} m^3".format(round(self.total_cavity_volume/m**3, 3)))\ return self.total_cavity_volume - # ELECTRON TRAP def TrapVolume(self): # Total volume of the electron traps in all cavities @@ -418,14 +463,14 @@ def EffectiveVolume(self): self.RF_background_rate_per_eV = self.Experiment.RF_background_rate_per_eV - #Radial efficiency - if self.Efficiency.unusable_dist_from_wall >= self.cyc_rad: - self.radial_efficiency = (self.cavity_radius - self.Efficiency.unusable_dist_from_wall)**2/self.cavity_radius**2 + #Radial efficiency - efficiency hit from Ioffe/Larmor bite + if self.unusable_dist_from_wall >= self.cyc_rad: + self.radial_efficiency = (self.cavity_radius - self.unusable_dist_from_wall)**2/self.cavity_radius**2 self.use_cyc_rad = False else: self.radial_efficiency = (self.cavity_radius - self.cyc_rad)**2/self.cavity_radius**2 self.use_cyc_rad = True - + #Efficiency from a cut during analysis on the axial frequency self.fa_cut_efficiency = trapping_efficiency(z_range = self.Experiment.trap_length /2, bg_magnetic_field = self.MagneticField.nominal_field, @@ -445,10 +490,6 @@ def BoxTrappingEfficiency(self): self.box_trapping_efficiency = np.cos(self.FrequencyExtraction.minimum_angle_in_bandwidth) return self.box_trapping_efficiency - def TrapLength(self): - self.Experiment.trap_length = 0.8 * 2 * self.cavity_radius * self.Experiment.cavity_L_over_D - logger.info("Calc'd trap length: {} m".format(round(self.Experiment.trap_length/m, 3), 2)) - def CavityPower(self): #Jprime_0 = 3.8317 max_ax_freq, mean_field, z_t = axial_motion(self.MagneticField.nominal_field, @@ -739,6 +780,28 @@ def syst_magnetic_field(self): else: return 0, 0 + def syst_plasma_effects(self): + if self.PlasmaEffects.UseFixedValue: + sigma = self.PlasmaEffects.Default_Systematic_Smearing + delta = self.PlasmaEffects.Default_Systematic_Uncertainty + return sigma, delta + elif not self.PlasmaEffects.UseFixedValue and self.Experiment.design_density_flag: + # C212 [s] - Differs from atomic calculator due to In + Es Crosssection for T-e at 18.6 keV and trap density + mean_track_duration = track_length(self.Experiment.number_density, self.T_endpoint, molecular=(not self.Experiment.atomic)) + #mean_track_duration = self.time_window + # C061 - Only used to estimate plasma broadening and mean track number present. Rough guess. + tracks_per_event = 10 + # C131 - Box trap approximation + min_pitch_angle_acceptance = 0.089 + # C217 - [eV] Plasma Broadening Calculation. Conservative upper limit based on dominance of 1 charge + # 2.35 converts stddev to FWHM + sigma = 7.2e-10 * (calculate_inventory(self.Experiment.design_density, self.total_cavity_volume) * lambda_tritium / Ci_Bq) * mean_track_duration * tracks_per_event \ + * Ci_Bq * min_pitch_angle_acceptance * np.log(self.cavity_length / self.cavity_radius) / 2.35 * eV + delta = self.PlasmaEffects.Default_Systematic_Uncertainty + return sigma, delta + else: + return 0, 0 + def det_efficiency_track_duration(self): """ Detection efficiency implemented based on René's slides, with faster and stable implementation using Gauss-Laguerre quadrature (G-L method): @@ -829,7 +892,7 @@ def print_SNRs(self, rho=None): track_duration = self.time_window logger.info("SNR-related parameters are printed for pre-set number density.") else: - track_duration = track_length(rho, self.T_endpoint, molecular=(not self.Experiment.atomic)) + track_duration = track_length(rho, self.T_endpoint, molecular=(not self.Experiment.atomic)) tau_snr_90deg = self.calculate_tau_snr(track_duration, power_fraction=1) #For an example carrier: @@ -872,10 +935,8 @@ def print_SNRs(self, rho=None): def print_Efficiencies(self): - logger.info("Effective volume: {} mm^3".format(round(self.effective_volume/mm**3, 3))) logger.info("Total efficiency: {}".format(self.effective_volume/self.total_trap_volume)) - if not self.Efficiency.usefixedvalue: # radial and detection efficiency are configured in the config file logger.info("Radial efficiency: {}".format(self.radial_efficiency)) @@ -885,7 +946,89 @@ def print_Efficiencies(self): logger.info("Efficiency from axial frequency cut: {}".format(self.fa_cut_efficiency)) logger.info("SRI factor: {}".format(self.Experiment.sri_factor)) - + def print_T2_background_atomic_trap(self): + #logger.info("T2 background: {}".format(self.Efficiency.T2_background_atomic_trap)) + if self.Efficiency.T2_background_atomic_trap: + self.T2_total_density, self.T2_T_ratio = calculate_T2_background_atomic_trap(self.cavity_radius, self.cavity_length, self.FrequencyExtraction.cavity_temperature, self.Efficiency.max_ratio_T2_T, self.Experiment.number_density) + if self.Efficiency.usefixedratio: + self.T2_T_ratio = self.Efficiency.T2_T_ratio + logger.info("T2_total_density: {:.4e} m^-3".format(self.T2_total_density*m**3)) + logger.info("Ratio T2/T: {:.4e}".format(self.T2_T_ratio)) + + def print_pumping_requirements(self): + #logger.info("Pumping Calculation: {}".format(self.Efficiency.pumping_calculation)) + if self.Efficiency.pumping_calculation: + self.turbopump_speed = turbopump_speed(self.Efficiency.number_turbopumps, self.Efficiency.turbopump_speed_gas_T2) + logger.info("Turbopump Speed: {:.4} m^3/s".format(self.turbopump_speed * s / m**3)) + self.turbopump_speed_limit_atomic = turbopump_speed_limit(self.cavity_radius, self.FrequencyExtraction.cavity_temperature, self.Experiment.cavity_L_over_D, self.Experiment.atomic) + self.cavity_termination_speed_atomic = cavity_termination_speed(self.Efficiency.turbopump_speed_cavity_termination_air, self.Efficiency.cavity_top_plate_temperature, self.Experiment.atomic) + self.turbopump_speed_limit_molecular = turbopump_speed_limit(self.cavity_radius, self.FrequencyExtraction.cavity_temperature, self.Experiment.cavity_L_over_D, not self.Experiment.atomic) + self.cavity_termination_speed_molecular = cavity_termination_speed(self.Efficiency.turbopump_speed_cavity_termination_air, self.Efficiency.cavity_top_plate_temperature, not self.Experiment.atomic) + logger.info("Pumping Speed Limit (Atomic): {:.4} m^3/s".format(self.turbopump_speed_limit_atomic * s / m**3)) + logger.info("Cavity Termination Speed (Atomic): {:.4} m^3/s".format(self.cavity_termination_speed_atomic * s / m**3)) + logger.info("Pumping Speed Limit (Molecular): {:.4} m^3/s".format(self.turbopump_speed_limit_molecular * s / m**3)) + logger.info("Cavity Termination Speed (Molecular): {:.4} m^3/s".format(self.cavity_termination_speed_molecular * s / m**3)) + if self.Efficiency.T_atom_supply_trap: + if self.Efficiency.usefixedvalue: + self.total_efficiency = self.Efficiency.fixed_efficiency + else: + self.total_efficiency = self.effective_volume/self.total_trap_volume + self.time_constant_He, self.current_He_leak = calculate_He_heat_leak(self.turbopump_speed_limit_atomic, self.cavity_termination_speed_atomic, self.turbopump_speed, self.FrequencyExtraction.cavity_temperature, self.Experiment.design_density, self.total_cavity_volume) + logger.info("He Time Constant: {:.4} s".format(self.time_constant_He/s)) + logger.info("Atom current required for He-3 leak: {:.4e} atoms/s".format(self.current_He_leak*s)) + self.time_constant_aperture, self.current_aperture_leak = calculate_aperture_heat_leak(self.DopplerBroadening.gas_temperature, self.Experiment.design_density, self.total_cavity_volume) + logger.info("Aperture Time Constant: {:.4} s".format(self.time_constant_aperture/s)) + logger.info("Atom current required for aperture leak: {:.4e} atoms/s".format(self.current_aperture_leak*s)) + self.time_constant_rad, self.current_rad_leak = calculate_rad_heat_leak(self.cavity_radius, self.Experiment.number_density, self.trap_coil_1, self.trap_coil_2, self.Experiment.design_density, self.total_cavity_volume, self.total_efficiency) + logger.info("Radiation Time Constant: {:.4} s".format(self.time_constant_rad/s)) + logger.info("Atom current required for radiation leak: {:.4e} atoms/s".format(self.current_rad_leak*s)) + self.time_constant_desorp, self.current_desorp = calculate_T2_desorption_from_wall(self.cavity_radius, self.cavity_length, self.Experiment.design_density, self.total_cavity_volume) + logger.info("Desorption Time Constant: {:.4} s".format(self.time_constant_desorp/s)) + logger.info("Atom current required for desorption losses: {:.4e} atoms/s".format(self.current_desorp*s)) + if not self.Experiment.pure_magnetic_flag: + logger.info("***Magnetogravitional Trap***") + else: + logger.info("***Pure Magnetic Trap***") + self.time_constant_dipolar, self.current_dipolar = calculate_dipolar_loss(self.MagneticField.nominal_field/T, self.cavity_radius, self.Experiment.design_density, self.total_cavity_volume, self.DopplerBroadening.gas_temperature, self.Experiment.cavity_L_over_D, self.top_cone, self.top_plate_cavity, self.Experiment.pure_magnetic_flag) + logger.info("Dipolar Time Constant: {:.4} s".format(self.time_constant_dipolar/s)) + logger.info("Atom current required for dipolar losses: {:.4e} atoms/s".format(self.current_dipolar*s)) + self.time_constant_evaporation, self.current_evaporation, self.evaporation_top, self.evaporation_sides, self.evaporation_cone = calculate_evaporation_loss(self.Experiment.cavity_L_over_D, self.Experiment.pure_magnetic_flag, self.DopplerBroadening.gas_temperature, self.cavity_radius, self.MagneticField.ioffe_field, self.MagneticField.nominal_field, self.cavity_length, self.top_cone, self.top_plate_cavity, self.Experiment.design_density, self.Efficiency.relative_loss_rate_cone_wall, self.total_cavity_volume) + logger.info("Evaporation Time Constant: {:.4} s".format(self.time_constant_evaporation/s)) + logger.info("Atom current required for Evaporation losses: {:.4e} atoms/s".format(self.current_evaporation*s)) + self.time_constant_T2, self.current_T2_leak = calculate_T2_heat_leak(self.current_aperture_leak, self.current_rad_leak, self.current_He_leak, self.current_evaporation, self.current_dipolar, self.T2_total_density, self.FrequencyExtraction.cavity_temperature, self.Experiment.design_density, self.total_cavity_volume) + logger.info("T2 Time Constant: {:.4} s".format(self.time_constant_T2/s)) + logger.info("Atom current required for T2 leak: {:.4e} atoms/s".format(self.current_T2_leak*s)) + self.trap_lifetime = calculate_trap_lifetime(self.time_constant_rad, self.time_constant_desorp, self.time_constant_He, self.time_constant_T2, self.time_constant_evaporation, self.time_constant_dipolar) + self.total_time_constant_trap, self.total_atom_current = calculate_trap_atom_supply(self.current_rad_leak, self.current_evaporation, self.current_dipolar, self.current_He_leak, self.current_T2_leak, self.Experiment.design_density, self.total_cavity_volume) + logger.info("Total trap time constant: {:.4} s".format(self.total_time_constant_trap/s)) + logger.info("Total atom current: {:.4e} atoms/s".format(self.total_atom_current*s)) + logger.info("Trap lifetime: {:.4} s".format(self.trap_lifetime/s)) + + self.molecular_density_allowed = molecular_density_allowed(self.Experiment.design_density, self.Efficiency.max_ratio_T2_T) + self.turbopump_speed_required = turbopump_speed_required(self.total_atom_current, self.molecular_density_allowed) + logger.info("Turbopump Speed required: {} m^3/s".format(self.turbopump_speed_required * s / m**3)) + self.ratio_turbopump_speed = ratio_turbopump_speed(self.turbopump_speed_required, self.turbopump_speed_limit_molecular, self.cavity_termination_speed_molecular, self.turbopump_speed) + logger.info("Ratio of Turbopump Speed Limit: {:.4e}".format(self.ratio_turbopump_speed)) + if (self.ratio_turbopump_speed >= 1.0): logger.info("CRITICAL PUMPING REACHED! Experiment cannot proceed with turbopumps.") + self.ratio_cryopump_speed, self.cryopump_speed, self.time_recycling_cryosurface = calculate_cryopump_speed(self.FrequencyExtraction.cavity_temperature, self.cavity_radius, self.cavity_length, self.turbopump_speed_required, self.total_atom_current, self.Experiment.design_density, self.Efficiency.max_ratio_T2_T) + logger.info("Ratio of Cryopump Speed: {:.4}".format(self.ratio_cryopump_speed)) + if (self.ratio_cryopump_speed >= 1.0): logger.info("CRITICAL PUMPING REACHED! The molecular pumping speed required is a dimensional calculation, not really a speed. As a result, the cryopumping speed limit is not real. Ideal cryopumping is the impingement rate: it all sticks. The actual molecular density is determined by the vapor pressure.") + logger.info("Cryopump speed of molecules: {:.4} m^3/s".format(self.cryopump_speed * s/m**3)) + logger.info("Time Recycling Cryosurface: {:.4} days".format(self.time_recycling_cryosurface / day)) + # C289 [Bq - Decay/s] - Accumulated activity in the sides by end of cycle + self.accumulated_activity_sides = self.evaporation_sides * self.time_recycling_cryosurface * (86400*s/day) * lambda_tritium / Ci_Bq + # C291 [Bq - Decay/s] - Accumulated activity in the cone by end of cycle + self.accumulated_activity_cone = self.evaporation_cone * self.time_recycling_cryosurface * (86400*s/day) * lambda_tritium / Ci_Bq + logger.info("Accumulated activity on the sides by end of cycle: {:.4} Ci".format(self.accumulated_activity_sides)) + logger.info("Accumulated activity on the cone by end of cycle: {:.4} Ci".format(self.accumulated_activity_cone)) + + self.vertical_injection_speed, self.trapped_gas_energy, self.injection_density, self.vertical_aperture_speed, self.injection_gas_energy = calculate_injection_line(self.total_atom_current, self.cavity_radius, self.Experiment.design_density, self.DopplerBroadening.gas_temperature, self.MagneticField.nominal_field, self.DopplerBroadening.injection_temperature, self.MagneticField.injection_field) + logger.info("Injection_density: {:.4e} m^-3".format(self.injection_density*m**3)) + if (self.injection_density*m**3 > 1e20): logger.info("Too high! Lower temperature and/or magnetic field in the injection beamline.") + logger.info("Vertical injection speed in trap at bottom: {:.4e} m/s".format(self.vertical_injection_speed * s/m)) + logger.info("Trapped gas energy: {:.4e} J".format(self.trapped_gas_energy/J)) + logger.info("Vertical injection speed in aperture: {:.4e} m/s".format(self.vertical_aperture_speed * s/m)) + logger.info("Injection gas energy: {:.4e} J".format(self.injection_gas_energy/J)) """ # Cramer-Rao lower bound / how much worse are we than the lower bound diff --git a/mermithid/sensitivity/SensitivityFormulas.py b/mermithid/sensitivity/SensitivityFormulas.py index 46187b75..07e0499e 100755 --- a/mermithid/sensitivity/SensitivityFormulas.py +++ b/mermithid/sensitivity/SensitivityFormulas.py @@ -9,9 +9,64 @@ ''' import numpy as np import configparser - +''' +from numpy import pi + +# Numericalunits is a package to handle units and some natural constants +# natural constants +from numericalunits import e, me, c0, eps0, kB, hbar +from numericalunits import meV, eV, keV, MeV, cm, m, ns, s, Hz, kHz, MHz, GHz, amu +from numericalunits import nT, uT, mT, T, mK, K, C, F, g, W +from numericalunits import hour, year, day +from numericalunits import mu0, NA, kB, hbar, me, c0, e, eps0, hPlanck + +T0 = -273.15*K + +tritium_livetime = 5.605e8*s +tritium_mass_atomic = 3.016* amu *c0**2 +tritium_electron_crosssection_atomic = 1.1e-22*m**2 +tritium_endpoint_atomic = 18563.251*eV +last_1ev_fraction_atomic = 2.067914e-13/eV**3 + +tritium_mass_molecular = 6.032099 * amu *c0**2 +tritium_electron_crosssection_molecular = 3.487*1e-22*m**2 +tritium_endpoint_molecular = 18573.24*eV +last_1ev_fraction_molecular = 1.67364e-13/eV**3 + +ground_state_width = 0.436 * eV +ground_state_width_uncertainty = 0.01*0.436*eV + +gyro_mag_ratio_proton = 42.577*MHz/T + +#Atomic tritium-related constant +ideal_gas_molar_volume = 2.24e-02*m**3 #Volume of 1mol of ideal gas at 1atm. Higher pressure --> smaller volume; higher temp --> lower volume. +#Cross sections are from Berlinsky +H_He_cross_sec_high_temp = 2.00e-19*m**2 #"High temp" means >5K +H_He_cross_sec_low_temp = 3.00e-20*m**2 +#H_H_cross_sec = +#H-H2 (289 pm kinetic diameter: Wikipedia) +#18.6-keV electrons on atomic T +#18.6-keV electrons on molecular T23.00E-20 + + +# units that do not show up in numericalunits +# missing pre-factors +fW = W*1e-15 + +# unitless units, relative fractions +pc = 0.01 +ppm = 1e-6 +ppb = 1e-9 +ppt = 1e-12 +ppq = 1e-15 + +# radian and degree which are also not really units +rad = 1 +deg = np.pi/180 +''' from mermithid.misc.Constants_numericalunits import * from mermithid.misc.CRESFunctions_numericalunits import * +from mermithid.sensitivity.AtomicCalculator import * try: from morpho.utilities import morphologging @@ -55,12 +110,13 @@ def __init__(self, config_path, verbose=True): self.Experiment = NameSpace({opt: eval(self.cfg.get('Experiment', opt)) for opt in self.cfg.options('Experiment')}) - # seetings fro molecular or atomic tritium + # settings for molecular or atomic tritium self.tau_tritium = tritium_livetime if self.Experiment.atomic: self.T_mass = tritium_mass_atomic self.Te_crosssection = tritium_electron_crosssection_atomic self.T_endpoint = tritium_endpoint_atomic + # Initially defined as number of events divided by energy window of last eV. For info: https://3.basecamp.com/3700981/buckets/3107037/documents/10099116348 self.last_1ev_fraction = last_1ev_fraction_atomic else: self.T_mass = tritium_mass_molecular @@ -90,16 +146,16 @@ def __init__(self, config_path, verbose=True): def SignalRate(self): """signal events in the energy interval before the endpoint, scale with DeltaE**3""" self.EffectiveVolume() - signal_rate = self.Experiment.number_density*self.effective_volume*self.last_1ev_fraction/self.tau_tritium + self.signal_rate = self.Experiment.number_density*self.effective_volume*self.last_1ev_fraction/self.tau_tritium if not self.Experiment.atomic: if hasattr(self.Experiment, 'gas_fractions'): avg_n_T_atoms = self.AvgNumTAtomsPerParticle_MolecularExperiment(self.Experiment.gas_fractions, self.Experiment.H2_type_gas_fractions) - signal_rate *= avg_n_T_atoms + self.signal_rate *= avg_n_T_atoms else: - signal_rate *= 2 + self.signal_rate *= 2 if hasattr(self.Experiment, 'active_gas_fraction'): - signal_rate *= self.Experiment.active_gas_fraction - return signal_rate + self.signal_rate *= self.Experiment.active_gas_fraction + return self.signal_rate def BackgroundRate(self): """background rate, can be calculated from multiple components. @@ -117,19 +173,58 @@ def BackgroundEvents(self): """Number of background events.""" return self.BackgroundRate()*self.Experiment.LiveTime*self.DeltaEWidth() + def SignalRatio(self): + if self.Efficiency.T2_background_atomic_trap: + self.T2_total_density, self.T2_T_ratio = calculate_T2_background_atomic_trap(self.cavity_radius, self.cavity_length, self.FrequencyExtraction.cavity_temperature, self.Efficiency.max_ratio_T2_T, self.Experiment.number_density) + if self.Efficiency.usefixedratio: + self.T2_T_ratio = self.Efficiency.T2_T_ratio + # C240 - Activity ratio of the rate of T2/T events in the last eV of the T2/T spectrum. (Difference between the two spectra) + signal_ratio = self.T2_T_ratio * 2 / ground_state_branch_atomic + """ + # Atomic Calculator activity in last eV of spectrum + if self.Efficiency.usefixedvalue: + self.total_efficiency = self.Efficiency.fixed_efficiency + else: + self.total_efficiency = self.effective_volume/self.total_trap_volume + sig_rate = calculate_activity_last_1eV_spectrum(self.Experiment.atomic, self.Experiment.number_density, self.cavity_radius, self.trap_coil_1, self.trap_coil_2, self.total_efficiency): + """ + return signal_ratio + return None + def DeltaEWidth(self): """optimal energy bin width""" labels, sigmas, deltas = self.get_systematics() - return np.sqrt(self.BackgroundRate()/self.SignalRate() - + 8*np.log(2)*(np.sum(sigmas**2))) + sig_rate = self.SignalRate() + bkg_rate = self.BackgroundRate() + zeta = self.SignalRatio() + #logger.info(f"zeta value: {repr(zeta)}") + #logger.info(f"zeta type: {type(zeta)}") + #logger.info(f"zeta dtype: {getattr(zeta, 'dtype', None)}") + #logger.info(f"zeta finite check: {np.isfinite(zeta)}") + if self.Efficiency.T2_background_atomic_trap: + # zeta = molecular background / atomic signal (r_m / r_a); if r_m = 0 then zeta = 0 and equation reduces to else statement + # delta_E = sqrt((b + 3 * delta_endpoint * zeta * r_a) / ((1 + zeta) * r_a) + 8*ln(2) * Systematics) + return np.sqrt((bkg_rate + 3 * endpoint_diff**2 * sig_rate * zeta) / (sig_rate * (1 + zeta)) \ + + 8*np.log(2)*(np.sum(sigmas**2))) + return np.sqrt(bkg_rate/sig_rate + 8*np.log(2)*(np.sum(sigmas**2))) def StatSens(self): """Pure statistic sensitivity assuming Poisson count experiment in a single bin As defined, it needs to be squared before being added to the systematic component""" + bkg_rate = self.BackgroundRate() sig_rate = self.SignalRate() DeltaE = self.DeltaEWidth() - sens = 2/(3*sig_rate*self.Experiment.LiveTime)*np.sqrt(sig_rate*self.Experiment.LiveTime*DeltaE - +self.BackgroundRate()*self.Experiment.LiveTime/DeltaE) + zeta = self.SignalRatio() + n_cavities = self.Experiment.n_cavities + if self.Efficiency.T2_background_atomic_trap: + # Stat^2 = 4/(9*r_a*t*#_cav) * ((1 + zeta)*delta_E + b*#_cav/r_a/delta_E + 3*#_cav*zeta*delta_endpoint*(1 + delta_endpoint/delta_E)) + sens = 2/(3*sig_rate*self.Experiment.LiveTime*n_cavities)*np.sqrt(sig_rate*self.Experiment.LiveTime*n_cavities*(1+zeta)*DeltaE \ + + bkg_rate*(n_cavities**2*self.Experiment.LiveTime/DeltaE) + 3*n_cavities**2*zeta*endpoint_diff*sig_rate \ + * self.Experiment.LiveTime*(1 + endpoint_diff/DeltaE)) + else: + # Stat^2 = 4/(9*r_a*t) * (delta_E + b/r_a/delta_E) if zeta=0 + sens = 2/(3*sig_rate*self.Experiment.LiveTime)*np.sqrt(sig_rate*self.Experiment.LiveTime*DeltaE \ + + bkg_rate*self.Experiment.LiveTime/DeltaE) return sens def SystSens(self): @@ -263,7 +358,7 @@ def syst_doppler_broadening(self): delta = self.DopplerBroadening.Default_Systematic_Uncertainty return sigma, delta - # termal doppler broardening + # thermal doppler broadening gasTemp = self.DopplerBroadening.gas_temperature mass_T = self.T_mass endpoint = self.T_endpoint @@ -282,7 +377,7 @@ def syst_doppler_broadening(self): Ee = endpoint + me*c0**2 p_rec = np.sqrt( Emax**2-me**2*c0**4 + (Emax - Ee - E_rec)**2 - mbeta**2 + 2*Ee*(Emax - Ee - E_rec)*betae*betanu*cosThetaenu ) sigma_trans = np.sqrt(p_rec**2/(2*mass_T)*2*kB*gasTemp) - + # sigma_trans = 2 * Ke * np.sqrt(kB_eV * trapped_gas_temp / (mass (eV) * beta^2)) if self.Experiment.atomic == True: delta_trans = np.sqrt(p_rec**2/(2*mass_T)*kB/gasTemp*self.DopplerBroadening.gas_temperature_uncertainty**2) else: @@ -428,4 +523,4 @@ def syst_plasma_effects(self): delta = self.PlasmaEffects.Default_Systematic_Uncertainty return sigma, delta else: - raise NotImplementedError("Plasma effect sysstematic is not implemented.") + raise NotImplementedError("Plasma effect systematic is not implemented.") diff --git a/test_analysis/fake_data_stan_analysis.py b/test_analysis/fake_data_stan_analysis.py index 00ebdce0..fa8e92f4 100644 --- a/test_analysis/fake_data_stan_analysis.py +++ b/test_analysis/fake_data_stan_analysis.py @@ -40,6 +40,7 @@ aposterioriPlotter = APosterioriDistribution("posterioriDistrib") divPlotter = Histo2dDivergence("2dDivergence") +Bfield = 0.9574762 def DefineGeneratorInputs(root_file='./results/tritium_analysis.root'): """ @@ -52,6 +53,7 @@ def DefineGeneratorInputs(root_file='./results/tritium_analysis.root'): "fixed_inputs": { 'Nscatters': 20, 'minf': 1353.125e+06 - 40e+06 + 24.5e+09, #In Hz + 'maxf': 1353.125e+06 + 2*40e+06 + 24.5e+09, 'err_from_B': 0. }, "priors": [ @@ -61,7 +63,7 @@ def DefineGeneratorInputs(root_file='./results/tritium_analysis.root'): {'name': 'S', 'prior_dist': 'poisson', 'prior_params': [3300.]}, {'name': 'B_1kev', 'prior_dist': 'lognormal', 'prior_params': [-3.826151398234498, 2.1673316326073935]}, {'name': 'survival_prob', 'prior_dist': 'beta', 'prior_params': [55, 17]}, #Centered approximately around 0.77. To be replaced given complex lineshape result/systematic assessment - {'name': 'Bfield', 'prior_dist': 'normal', 'prior_params': [0.9574762, 1.51e-6]}, #From complex lineshape fit to calibration data. More sig figs on mean needed? + {'name': 'Bfield', 'prior_dist': 'normal', 'prior_params': [Bfield, 1.51e-6]}, #From complex lineshape fit to calibration data. More sig figs on mean needed? ] } @@ -72,6 +74,7 @@ def DefineGeneratorInputs(root_file='./results/tritium_analysis.root'): "variables": [ {"variable": "Nscatters", "type":"int"}, {"variable": "minf", "type": "float"}, + {"variable": "maxf", "type": "float"}, {"variable": "err_from_B", "type": "float"}, {"variable": "Q", "type": "float"}, {"variable": "mass", "type": "float"}, @@ -112,12 +115,13 @@ def GenerateFakeData(inputs_dict): """ specGen_config = { "apply_efficiency": True, - "efficiency_path": "../phase2_detection_efficiency_curve/combined_energy_corrected_count_rates/combined_energy_corrected_eff_at_quad_trap_frequencies.json", + "efficiency_path": "../tests/combined_energy_corrected_eff_at_quad_trap_frequencies.json", "detailed_or_simplified_lineshape": "detailed", "return_frequency": True, "Q": inputs_dict["Q"], "mass": inputs_dict["mass"], "minf": inputs_dict["minf"], + "maxf": inputs_dict["maxf"], "scattering_sigma": inputs_dict["sigma"], "S": inputs_dict["S"], "B_1kev": inputs_dict["B_1kev"], @@ -126,6 +130,7 @@ def GenerateFakeData(inputs_dict): "Nscatters": inputs_dict["Nscatters"], "B_field": inputs_dict["Bfield"], "n_steps": 100000, + "use_lineshape": False } histo_config = { @@ -154,7 +159,7 @@ def BinAndSaveData(tritium_data, nbins, root_file="./results/tritium_analysis.ro "energy_or_frequency": 'frequency', "variables": "F", "title": "corrected_spectrum", - "efficiency_filepath": "../phase2_detection_efficiency_curve/combined_energy_corrected_count_rates/combined_energy_corrected_eff_at_quad_trap_frequencies.json", + "efficiency_filepath": "../tests/combined_energy_corrected_eff_at_quad_trap_frequencies.json", 'bins': np.linspace(tritium_data['minf'], tritium_data['maxf'], nbins), 'fss_bins': False # If fss_bins is True, bins is ignored and overridden } @@ -210,7 +215,7 @@ def SaveUnbinnedData(tritium_data, root_file="./results/tritium_analysis.root"): -def StanTritiumAnalysis(tritium_data, fit_parameters=None, root_file='./results/tritium_analysis.root', stan_files_location='../../morpho_models/', model_code='tritium_model/models/tritium_phase_II_analyzer_binned.stan', scattering_params_R='simplified_scattering_params.R'): +def StanTritiumAnalysis(tritium_data, fit_parameters=None, root_file='./results/tritium_analysis.root', stan_files_location='./', model_code='main_model/tritium_phase_II_analyzer_binned.stan', scattering_params_R='simplified_scattering_params.R'): """ Analyzes frequency or kinetic energy data using a Stan model. Saves and plots posteriors. @@ -241,14 +246,14 @@ def StanTritiumAnalysis(tritium_data, fit_parameters=None, root_file='./results/ "cache_dir": stan_files_location+"tritium_model/cache", "warmup": 4000, #Increase for real run (to 3000-5000) "iter": 8000, #Increase for real run (to 6000-9000) - "chain": 3, #Increase for real run (to 3-4) + "chain": 1, #Increase for real run (to 3-4) "control": {'adapt_delta':0.97}, "init": { "sigma": 17.7084, "survival_prob": 0.77, "Q": 18573.24, "mass": 0.2, - "Bfield": 0.9574762, + "Bfield": Bfield, "S": 3300., "B_1kev": 0.0217933282798889, "B": 0.3, @@ -257,7 +262,7 @@ def StanTritiumAnalysis(tritium_data, fit_parameters=None, root_file='./results/ "sigma_ctr": 17.7084, #sigma params from Ali's complex lineshape fits. "sigma_std": 1.14658, #Final states not yet included. "err_from_B": 0.001, #Tiny smearing from f_c->K conversion - "Bfield_ctr": 0.9574762, #From complex lineshape fit to calibration + "Bfield_ctr": Bfield, #From complex lineshape fit to calibration "Bfield_std": 1.51e-06, #data. More sig figs on mean needed? "survival_prob_alpha": 55, #Centered around ~0.77. To be replaced "survival_prob_beta": 17, #given complex lineshape result+systematics @@ -270,7 +275,9 @@ def StanTritiumAnalysis(tritium_data, fit_parameters=None, root_file='./results/ "KEscale": 16323, #This enables the option of cmdstan running # "slope": 0.000390369173, #For efficiency modeling with unbinned data # "intercept": -6.00337656, - "Nscatters": 16 #Because peaks>16 in simplified linesahpe have means->inf as FWHM->0 + "Nscatters": 16, #Because peaks>16 in simplified linesahpe have means->inf as FWHM->0 + "slope": 1, + "intercept": 20380.5153 }, "interestParams": ['Q', 'mass', 'survival_prob', 'Bfield', 'sigma', 'S', 'B_1kev', 'KEmin', 'KE_sample', 'Nfit_signal', 'Nfit_bkgd'], } @@ -410,7 +417,7 @@ def CalibrateResults(root_filenames, vars_to_calibrate, cred_interval=[0.05, 0.9 calibrator.Run() -def FakeExperimentEnsemble(n_runs, root_basename, wait_before_runs=0, parallelize=True, n_processes=4, vars_to_calibrate=['Q']): +def FakeExperimentEnsemble(n_runs, root_basename, wait_before_runs=0, parallelize=True, n_processes=1, vars_to_calibrate=['Q']): """ To-do: add parallelization option for a Slurm environment. @@ -438,8 +445,8 @@ def FakeExperimentEnsemble(n_runs, root_basename, wait_before_runs=0, paralleliz with Pool(n_processes) as p: p.map(PerformFakeExperiment, root_filenames) - coverages = CalibrateResults(root_filenames, vars_to_calibrate) - + #coverages = CalibrateResults(root_filenames, vars_to_calibrate) + coverages = None return coverages