Feature/sensitivity curve link rf background and detection efficiency to snr

mermithid/misc/CRESFunctions_numericalunits.py

@@ -28,6 +28,10 @@ def frequency(kin_energy, magnetic_field):
     # cyclotron frequency
     return e/(2*np.pi*me)/gamma(kin_energy)*magnetic_field
 
+def cyclotron_radius(freq, kin_energy):
+    v = beta(kin_energy)*c0
+    return v/freq/(2*np.pi)
+
 def wavelength(kin_energy, magnetic_field):
     return c0/frequency(kin_energy, magnetic_field)
 

mermithid/sensitivity/SensitivityCavityFormulas.py

@@ -1,12 +1,15 @@
 '''
 Class calculating neutrino mass sensitivities based on analytic formulas from CDR.
-Author: R. Reimann, C. Claessens, T. Weiss, W. Van De Pontseele
-Date:06/07/2023
+Author: R. Reimann, C. Claessens, T. E. Weiss, W. Van De Pontseele
+Date: 06/07/2023
+Updated: December 2024
 
 The statistical method and formulas are described in
 CDR (CRES design report, Section 1.3) https://www.overleaf.com/project/5b9314afc673d862fa923d53.
 '''
 import numpy as np
+from scipy.stats import ncx2, chi2
+from scipy.integrate import quad
 
 from mermithid.misc.Constants_numericalunits import *
 from mermithid.misc.CRESFunctions_numericalunits import *
@@ -191,53 +194,78 @@ class CavitySensitivity(Sensitivity):
     """
     def __init__(self, config_path):
         Sensitivity.__init__(self, config_path)
-        self.Efficiency = NameSpace({opt: eval(self.cfg.get('Efficiency', opt)) for opt in self.cfg.options('Efficiency')})
 
+        ###
+        #Initialization related to the effective volume:
+        ###
         self.Jprime_0 = 3.8317
-        self.CRLB_constant = 12
-        #self.CRLB_constant = 90
-        if hasattr(self.FrequencyExtraction, "crlb_constant"):
-            self.CRLB_constant = self.FrequencyExtraction.crlb_constant
-            logger.info("Using configured CRLB constant")   
+        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'):
+            self.Experiment.trap_length = 2 * self.cavity_radius * self.Experiment.cavity_L_over_D
+
+        self.Efficiency = NameSpace({opt: eval(self.cfg.get('Efficiency', opt)) for opt in self.cfg.options('Efficiency')})
+        self.CavityVolume()
+        self.CavityPower()
 
         #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, 
                                                                     min_pitch_angle = self.FrequencyExtraction.minimum_angle_in_bandwidth, 
                                                                     trap_flat_fraction = self.MagneticField.trap_flat_fraction
-                                                                    )             
+                                                                    )          
+
+        #We may decide to remove the "Threshold" section and move the threshold-related parameters to the "Efficiency" section.
+        if not self.Efficiency.usefixedvalue:
+            self.Threshold = NameSpace({opt: eval(self.cfg.get('Threshold', opt)) for opt in self.cfg.options('Threshold')})
+
+        #Cyclotron radius is sometimes used in the effective volume calculation
+        self.cyc_rad = cyclotron_radius(self.cavity_freq, self.T_endpoint) 
+
+        #Calculate the effective volume and print out related quantities
+        self.EffectiveVolume()
+        logger.info("Trap radius: {} cm".format(round(self.cavity_radius/cm, 3), 2))
+        logger.info("Total trap volume: {} m^3".format(round(self.total_trap_volume/m**3), 2))
+        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")
+
+        ####
+        #Initialization related to the energy resolution:
+        ####
+        self.CRLB_constant = 12
+        #self.CRLB_constant = 90
+        if hasattr(self.FrequencyExtraction, "crlb_constant"):
+            self.CRLB_constant = self.FrequencyExtraction.crlb_constant
+            logger.info("Using configured CRLB constant")      
 
         #Numbr of steps in pitch angle between min_pitch and pi/2 for the frequency noise uncertainty calculation
         self.pitch_steps = 100
         if hasattr(self.FrequencyExtraction, "pitch_steps"):
             self.pitch_steps = self.FrequencyExtraction.pitch_steps
             logger.info("Using configured pitch_steps value")  
 
-        self.CavityRadius()
-        self.CavityVolume()
-        self.EffectiveVolume()
-        self.CavityPower()
-
-        #Get trap length from cavity length if not specified
-        if self.Experiment.trap_length == 0:
-            self.Experiment.trap_length = 2 * self.cavity_radius * self.Experiment.cavity_L_over_D
+        #Just calculated for comparison
+        self.larmor_power = rad_power(self.T_endpoint, np.pi/2, self.MagneticField.nominal_field) # currently not used
 
 
     # CAVITY
     def CavityRadius(self):
         axial_mode_index = 1
-        self.cavity_radius = c0/(2*np.pi*frequency(self.T_endpoint, self.MagneticField.nominal_field))*np.sqrt(self.Jprime_0**2+axial_mode_index**2*np.pi**2/(4*self.Experiment.cavity_L_over_D**2))
+        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
 
     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
 
-        logger.info("Frequency: {} MHz".format(round(frequency(self.T_endpoint, self.MagneticField.nominal_field)/MHz, 3)))
+        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("Total cavity volume {} m^3".format(round(self.total_cavity_volume/m**3)))\
+        logger.info("Total cavity volume: {} m^3".format(round(self.total_cavity_volume/m**3, 3)))\
 
         return self.total_cavity_volume
 
@@ -246,42 +274,50 @@ def CavityVolume(self):
     def TrapVolume(self):
         # Total volume of the electron traps in all cavities
         self.total_trap_volume = self.Experiment.trap_length*np.pi*(self.cavity_radius)**2*self.Experiment.n_cavities
-
-        logger.info("Trap radius: {} cm".format(round(self.cavity_radius/cm, 3)))
-        logger.info("Trap length: {} cm".format(round(self.Experiment.trap_length/cm, 3)))
-        logger.info("Total trap volume {} m^3".format(round(self.total_trap_volume/m**3)))
-
         return self.total_trap_volume
 
 
 
     def EffectiveVolume(self):
         self.total_trap_volume = self.TrapVolume()
+
         if self.Efficiency.usefixedvalue:
             self.effective_volume = self.total_trap_volume * self.Efficiency.fixed_efficiency
+            self.use_cyc_rad = False
         else:
-            # radial and detection efficiency are configured in the config file
-            #logger.info("Radial efficiency: {}".format(self.Efficiency.radial_efficiency))
-            #logger.info("Detection efficiency: {}".format(self.Efficiency.detection_efficiency))
-            #logger.info("Pitch angle efficiency: {}".format(self.PitchDependentTrappingEfficiency()))
-            #logger.info("SRI factor: {}".format(self.Experiment.sri_factor))
+            #Detection efficiency
+            if self.Threshold.use_detection_threshold:
+                #If the_threshold is True, calculate the detection efficiency given the SNR of data and the threshold.
+                self.assign_background_rate_from_threshold()
+                self.assign_detection_efficiency_from_threshold()
+            else:
+                #If use_threshold is False, use the inputted detection efficiency and RF background rate from the config file.
+                self.detection_efficiency = self.Efficiency.detection_efficiency
+                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
+                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
 
-            self.radial_efficiency = (self.cavity_radius - self.Efficiency.unusable_dist_from_wall)**2/self.cavity_radius**2
-            self.fa_cut_efficiency = trapping_efficiency( z_range = self.Experiment.trap_length /2,
+            #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, 
                                                                     min_pitch_angle = self.Efficiency.min_pitch_used_in_analysis, 
                                                                     trap_flat_fraction = self.MagneticField.trap_flat_fraction
                                                                     )/self.pos_dependent_trapping_efficiency 
-            self.effective_volume = self.total_trap_volume*self.radial_efficiency*self.Efficiency.detection_efficiency*self.fa_cut_efficiency*self.pos_dependent_trapping_efficiency   
 
-        #logger.info("Total efficiency: {}".format(self.effective_volume/self.total_volume))        
+            #The effective volume includes the three efficiency factors above, as well as the trapping efficiency
+            self.effective_volume = self.total_trap_volume*self.radial_efficiency*self.detection_efficiency*self.fa_cut_efficiency*self.pos_dependent_trapping_efficiency   
+
+        # The "signal rate improvement" factor can be toggled to test the increase in statistics required to reach some sensitivity
         self.effective_volume*=self.Experiment.sri_factor
-
-        # for parent SignalRate function
-        # self.Experiment.v_eff = self.effective_volume
-
         return self.effective_volume
 
+
     def BoxTrappingEfficiency(self):
         self.box_trapping_efficiency = np.cos(self.FrequencyExtraction.minimum_angle_in_bandwidth)
         return self.box_trapping_efficiency
@@ -300,7 +336,7 @@ def CavityPower(self):
                                                                             z_t, self.CavityLoadedQ(), 
                                                                             2*self.Experiment.cavity_L_over_D*self.cavity_radius, 
                                                                             self.cavity_radius, 
-                                                                            frequency(self.T_endpoint, self.MagneticField.nominal_field)))
+                                                                            self.cavity_freq))
         return self.signal_power
 
 
@@ -311,7 +347,7 @@ def CavityLoadedQ(self):
 
         #self.loaded_q =1/(0.22800*((90-self.FrequencyExtraction.minimum_angle_in_bandwidth)*np.pi/180)**2+2**2*0.01076**2/(4*0.22800))
 
-        endpoint_frequency = frequency(self.T_endpoint, self.MagneticField.nominal_field)
+        endpoint_frequency = self.cavity_freq
         #required_bw_axialfrequency = axial_frequency(self.Experiment.cavity_L_over_D*self.CavityRadius()*2, 
         #                                             self.T_endpoint, 
         #                                             self.FrequencyExtraction.minimum_angle_in_bandwidth/deg)
@@ -342,7 +378,7 @@ def calculate_tau_snr(self, time_window, power_fraction=1):
         power_fraction may be used as a carrier or a sideband power fraction,
         relative to the power of a 90 degree carrier.
         """
-        endpoint_frequency = frequency(self.T_endpoint, self.MagneticField.nominal_field)
+        endpoint_frequency = self.cavity_freq
 
         # Cavity coupling
         self.CavityLoadedQ()
@@ -412,15 +448,14 @@ def syst_frequency_extraction(self):
             return sigma, delta
 
 
-        endpoint_frequency = frequency(self.T_endpoint, self.MagneticField.nominal_field)
+        endpoint_frequency = self.cavity_freq
         # using Pe and alpha (aka slope) from above
-        Pe = self.signal_power #/self.FrequencyExtraction.mode_coupling_efficiency 
-        self.larmor_power = rad_power(self.T_endpoint, np.pi/2, self.MagneticField.nominal_field) # currently not used
+        Pe = self.signal_power #/self.FrequencyExtraction.mode_coupling_efficiency
 
         self.slope = endpoint_frequency * 2 * np.pi * Pe/me/c0**2 # track slope
         self.time_window = track_length(self.Experiment.number_density, self.T_endpoint, molecular=(not self.Experiment.atomic))
 
-        self.time_window_slope_zero = abs(frequency(self.T_endpoint, self.MagneticField.nominal_field)-frequency(self.T_endpoint+20*meV, self.MagneticField.nominal_field))/self.slope
+        self.time_window_slope_zero = abs(self.cavity_freq-frequency(self.T_endpoint+20*meV, self.MagneticField.nominal_field))/self.slope
 
         tau_snr_full_length = self.calculate_tau_snr(self.time_window, self.FrequencyExtraction.carrier_power_fraction)
         tau_snr_part_length = self.calculate_tau_snr(self.time_window_slope_zero, self.FrequencyExtraction.carrier_power_fraction)
@@ -470,7 +505,7 @@ def syst_frequency_extraction(self):
                                     np.pi/2-phis, self.Experiment.trap_length,
                                     self.FrequencyExtraction.minimum_angle_in_bandwidth, 
                                     self.T_endpoint, flat_fraction=self.MagneticField.trap_flat_fraction)
-            fc0_endpoint = frequency(self.T_endpoint, self.MagneticField.nominal_field)
+            fc0_endpoint = self.cavity_freq
             p_array = ax_freq_array/fc0_endpoint/phis
             self.p = np.mean(p_array[1:]) #Cut out theta=pi/2 (ill defined there)
 
@@ -545,6 +580,35 @@ def syst_magnetic_field(self):
             return sigma, frac_uncertainty*sigma
         else:
             return 0, 0
+
+    def det_efficiency_track_duration(self):
+        # Detection efficiency implemented based on René's slides:
+        # https://3.basecamp.com/3700981/buckets/3107037/documents/8013439062
+        # Also check the antenna paper for more details. 
+        # Especially the section on the signal detection with matched filtering.
+        mean_track_duration = track_length(self.Experiment.number_density, self.T_endpoint, molecular=(not self.Experiment.atomic))
+        tau_snr_ex_carrier = self.calculate_tau_snr(mean_track_duration, self.FrequencyExtraction.carrier_power_fraction)
+        result, abs_err = quad(lambda track_duration: ncx2(df=2, nc=track_duration/tau_snr_ex_carrier).sf(self.Threshold.detection_threshold)*1/mean_track_duration*np.exp(-track_duration/mean_track_duration), 0, np.inf)
+        return result, abs_err
+
+    def assign_detection_efficiency_from_threshold(self):
+        self.detection_efficiency, self.abs_err = self.det_efficiency_track_duration()
+        return self.detection_efficiency
+
+    def rf_background_rate_cavity(self):
+        # Detection efficiency implemented based on René's slides
+        # https://3.basecamp.com/3700981/buckets/3107037/documents/8013439062
+        # Also check the antenna paper for more details, especially the section
+        # on the signal detection with matched filtering.
+        # Assuming background rate constant of 1/(eV*s) for now. This constant 
+        # will need to be determined from Monte Carlo simulations.
+        return chi2(df=2).sf(self.Threshold.detection_threshold)/(eV*s)
+
+    def assign_background_rate_from_threshold(self):
+        self.RF_background_rate_per_eV = self.rf_background_rate_cavity()
+        return self.RF_background_rate_per_eV
+
+
 
 
     # PRINTS
@@ -563,7 +627,7 @@ def print_SNRs(self, rho=None):
         tau_snr_ex_carrier = self.calculate_tau_snr(track_duration, self.FrequencyExtraction.carrier_power_fraction)
 
 
-        eV_bandwidth = np.abs(frequency(self.T_endpoint, self.MagneticField.nominal_field) - frequency(self.T_endpoint + 1*eV, self.MagneticField.nominal_field))
+        eV_bandwidth = np.abs(self.cavity_freq - frequency(self.T_endpoint + 1*eV, self.MagneticField.nominal_field))
         SNR_1eV_90deg = 1/eV_bandwidth/tau_snr_90deg
         SNR_track_duration_90deg = track_duration/tau_snr_90deg
         SNR_1ms_90deg = 0.001*s/tau_snr_90deg
@@ -601,7 +665,8 @@ def print_Efficiencies(self):
         if not self.Efficiency.usefixedvalue:
             # radial and detection efficiency are configured in the config file
             logger.info("Radial efficiency: {}".format(self.radial_efficiency))
-            logger.info("Detection efficiency: {}".format(self.Efficiency.detection_efficiency))
+            logger.info("Detection efficiency: {}".format(self.detection_efficiency))
+            logger.info("Detection efficiency integration error: {}".format(self.abs_err))
             logger.info("Trapping efficiency: {}".format(self.pos_dependent_trapping_efficiency))
             logger.info("Efficiency from axial frequency cut: {}".format(self.fa_cut_efficiency))
             logger.info("SRI factor: {}".format(self.Experiment.sri_factor))
@@ -650,4 +715,4 @@ def print_Efficiencies(self):
 
 """fc_endpoint_array = frequency(self.T_endpoint, mean_field_array)
 self.q_array = 1/phis**2*(fc_endpoint_array/fc0_endpoint - 1)
-self.q = np.mean(self.q_array[1:])"""
+self.q = np.mean(self.q_array[1:])"""

mermithid/sensitivity/SensitivityFormulas.py

@@ -89,6 +89,7 @@ def __init__(self, config_path):
     # SENSITIVITY
     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
         if not self.Experiment.atomic:
             if hasattr(self.Experiment, 'gas_fractions'):
@@ -105,7 +106,7 @@ def BackgroundRate(self):
         Currently, RF noise and cosmic ray backgrounds are included.
         Assumes that background rate is constant over considered energy / frequency range."""
         self.cosmic_ray_background = self.Experiment.cosmic_ray_bkgd_per_tritium_particle*self.Experiment.number_density*self.effective_volume
-        self.background_rate = self.Experiment.RF_background_rate_per_eV + self.cosmic_ray_background
+        self.background_rate = self.RF_background_rate_per_eV + self.cosmic_ray_background
         return self.background_rate
 
     def SignalEvents(self):