diff --git a/flamedisx/__init__.py b/flamedisx/__init__.py index a26bc3a00..bacbb748f 100644 --- a/flamedisx/__init__.py +++ b/flamedisx/__init__.py @@ -33,10 +33,6 @@ # Access through fd.nest.xxx from . import nest -# LUX models -# Access through fd.lux.xxx -from . import lux - # LZ models # Access through fd.lz.xxx from . import lz diff --git a/flamedisx/lux/__init__.py b/flamedisx/lux/__init__.py deleted file mode 100644 index 17ddcf019..000000000 --- a/flamedisx/lux/__init__.py +++ /dev/null @@ -1 +0,0 @@ -from .lux import * diff --git a/flamedisx/lux/lux.py b/flamedisx/lux/lux.py deleted file mode 100644 index 005c98aba..000000000 --- a/flamedisx/lux/lux.py +++ /dev/null @@ -1,132 +0,0 @@ -"""lux detector implementation - -""" -import tensorflow as tf - -import configparser -import os - -import flamedisx as fd -from .. import nest as fd_nest - -export, __all__ = fd.exporter() - - -## -# Flamedisx sources -## - -## -# Common to all LUX sources -## - - -class LUXSource: - def __init__(self, *args, **kwargs): - super().__init__(*args, **kwargs) - - assert kwargs['detector'] in ('default',) - - assert os.path.exists(os.path.join( - os.path.dirname(__file__), '../nest/config/', kwargs['detector'] + '.ini')) - - config = configparser.ConfigParser(inline_comment_prefixes=';') - config.read(os.path.join(os.path.dirname(__file__), '../nest/config/', - kwargs['detector'] + '.ini')) - - self.z_topDrift = config.getfloat('NEST', 'z_topDrift_config') - self.dt_cntr = config.getfloat('NEST', 'dt_cntr_config') - - self.density = fd_nest.calculate_density( - config.getfloat('NEST', 'temperature_config'), - config.getfloat('NEST', 'pressure_config')).item() - self.drift_velocity = fd_nest.calculate_drift_velocity( - config.getfloat('NEST', 'drift_field_config'), - self.density, - config.getfloat('NEST', 'temperature_config')).item() - - def s1_posDependence(self, r, z): - """ - Returns position-dependent S1 scale factor. - Requires r/z to be in cm, and in the FV. - """ - r_mm = r*10 - z_mm = z*10 - - amplitude = 307.9 - 0.3071 * z_mm + 0.0002257 * pow(z_mm, 2) - shape = 1.1525e-7 * tf.sqrt(abs(z_mm - 318.84)) - finalCorr = (-shape * pow(r_mm, 3) + amplitude) / 307.9 - - # We want to normalise to the detector centre - z_cntr = self.z_topDrift - self.drift_velocity * self.dt_cntr - z_cntr_mm = z_cntr * 10 - - amplitude_0 = 307.9 - 0.3071 * z_cntr_mm + 0.0002257 * pow(z_cntr_mm, 2) - finalCorr_0 = amplitude_0 / 307.9 - - return finalCorr / finalCorr_0 - - @staticmethod - def s2_posDependence(r): - """ - Returns position-dependent S2 scale factor. - Requires r to be in cm, and in the FV. - """ - r_mm = r*10 - - finalCorr = 9156.3 + 6.22750 * pow(r_mm, 1) + 0.38126 * pow(r_mm, 2) \ - - 0.017144 * pow(r_mm, 3) + 0.0002474 * pow(r_mm, 4) \ - - 1.6953e-6 * pow(r_mm, 5) + 5.6513e-9 * pow(r_mm, 6) \ - - 7.3989e-12 * pow(r_mm, 7) - - return finalCorr / 9156.3 - - -## -# Different interaction types: flat spectra -## - - -@export -class LUXERSource(LUXSource, fd.nest.nestERSource): - def __init__(self, *args, **kwargs): - if ('detector' not in kwargs): - kwargs['detector'] = 'default' - super().__init__(*args, **kwargs) - - -@export -class LUXGammaSource(LUXSource, fd.nest.nestGammaSource): - def __init__(self, *args, **kwargs): - if ('detector' not in kwargs): - kwargs['detector'] = 'default' - super().__init__(*args, **kwargs) - - -@export -class LUXERGammaWeightedSource(LUXSource, fd.nest.nestERGammaWeightedSource): - def __init__(self, *args, **kwargs): - if ('detector' not in kwargs): - kwargs['detector'] = 'default' - super().__init__(*args, **kwargs) - - -@export -class LUXNRSource(LUXSource, fd.nest.nestNRSource): - def __init__(self, *args, **kwargs): - if ('detector' not in kwargs): - kwargs['detector'] = 'default' - super().__init__(*args, **kwargs) - - -## -# Signal sources -## - - -@export -class LUXWIMPSource(LUXSource, fd.nest.nestWIMPSource): - def __init__(self, *args, **kwargs): - if ('detector' not in kwargs): - kwargs['detector'] = 'default' - super().__init__(*args, **kwargs) diff --git a/flamedisx/lz/WS2024_cuts_and_acceptances.py b/flamedisx/lz/WS2024_cuts_and_acceptances.py new file mode 100644 index 000000000..e672856c1 --- /dev/null +++ b/flamedisx/lz/WS2024_cuts_and_acceptances.py @@ -0,0 +1,144 @@ +from scipy.interpolate import RegularGridInterpolator +import numpy as np +##=============Cuts and acceptances=========== +def WS2024_S2splitting_reconstruction_efficiency(S2c, driftTime_us, hist): + """ + Returns the reconstruction efficiency based on S2 splitting + S2c: cs2 from data [phd] + dritTime_us: drift time from data [us] + hist: input histogram of S2c[Ne-] vs Drift Time [us] + adapted from BGSkimmer + """ + ## Make values more python friendly + weights = np.reshape(hist[0], hist[0].size) + zmax = np.max(weights) + zmin = np.min(weights[weights>0.]) + ## Read into interpolator + xentries_vals = np.array(hist[1][:-1]) + yentries_vals = np.array(hist[2][:-1]) + Interp = RegularGridInterpolator((xentries_vals,yentries_vals), hist[0]) + + ## Convert S2c to mean N_e + mean_SE_area = 44.5 # phd/e- + mean_Ne = S2c/mean_SE_area + + ## Initialize acceptance values + acceptance = np.ones_like(mean_Ne) + ## curves not defined for above 100 e- + ## Assume 100% eff. (probably ok...) + acceptance[mean_Ne>np.max(xentries_vals)] = 1. + ## Also not defined for < 10e- + ## ok with 15e- ROI threshold + acceptance[mean_Nenp.max(yentries_vals)] = np.max(yentries_vals) + + mask = (mean_Ne>=np.min(xentries_vals)) & (mean_Ne<=np.max(xentries_vals)) + ## Acceptances are provided in percent - divide by 100. + acceptance[mask] = Interp(np.vstack([mean_Ne[mask],temp_drift[mask]]).T)/100. + + + return acceptance + +##=====================CUTS===================== +##Build S2raw acceptance +def WS2024_trigger_acceptance(S2raw): + """ + Returns the Trigger efficiency of S2s + S2rw: s2 from data [phd] + adapted from BGSkimmer + """ + # 50% threshold = 106.89 +/- 0.43 phd + # 95% threshold = 160.7 +/- 3.5 phd + # k = 0.0547 +/- 0.0035 phd-1 + + x0 = 160.7 + k = 0.0547 + + accValues = ( 1./( 1 + np.exp( -k*(S2raw-x0) ) ) ) + + ## make sure acceptances can't go above 1. or below 0. + accValues[accValues>1.] = 1. + accValues[accValues<0.] = 0. + + return accValues + + +# Define polynomial coefficients associated with each azimuthal wall position; +# these start at the 7 o'clock position and progress anti-clockwise +phi_coeffs = [[-1.78880746e-13, 4.91268301e-10, -4.96134607e-07, 2.26430932e-04, -4.71792008e-02, 7.33811298e+01], + [-1.72264463e-13, 4.59149636e-10, -4.59325165e-07, 2.14612376e-04, -4.85599108e-02, 7.35290867e+01], + [-3.17099156e-14, 7.26336129e-11, -6.99495385e-08, 3.85531008e-05, -1.33386004e-02, 7.18002889e+01], + [-6.12280314e-14, 1.67968911e-10, -1.83625538e-07, 1.00457608e-04, -2.86728022e-02, 7.22754350e+01], + [-1.89897962e-14, 1.52777215e-11, -2.79681508e-09, 1.25689887e-05, -1.33093804e-02, 7.17662251e+01], + [-2.32118621e-14, 7.30043322e-11, -9.40606298e-08, 6.29728588e-05, -2.28150175e-02, 7.22661091e+01], + [-8.29749194e-14, 2.31096069e-10, -2.47867121e-07, 1.27576029e-04, -3.24702414e-02, 7.26357609e+01], + [-2.00718008e-13, 5.44135757e-10, -5.59484466e-07, 2.73028553e-04, -6.46879791e-02, 7.45264998e+01], + [-7.77420021e-14, 1.97357045e-10, -1.90016273e-07, 8.99659454e-05, -2.30169916e-02, 7.25038258e+01], + [-5.27296334e-14, 1.49415580e-10, -1.58205132e-07, 8.00275441e-05, -2.13559394e-02, 7.23995451e+01], + [-6.00198219e-14, 1.55333004e-10, -1.60367908e-07, 7.97754165e-05, -1.94435594e-02, 7.22714399e+01], + [-8.89919309e-14, 2.40830027e-10, -2.57060475e-07, 1.33002951e-04, -3.32969110e-02, 7.28696020e+01]] + + +# Use the above set of coefficients to define azimuthal wall position contours +phi_walls = [np.poly1d(phi_coeffs[i]) for i in range(len(phi_coeffs))] +FV_poly = np.poly1d([-4.44147071e-14, 1.43684777e-10, -1.82739476e-07, + 1.02160174e-04, -2.31617857e-02, -2.05932471e+00]) +def WS2024_fiducial_volume_cut(x,y,dt): + """ + Fiducial volume cute for WS2024 + x: x position [cm] + y: y position [cm] + dt: drift time [us] + adapted from BGSkimmer + """ + # Define radial contour for N_tot = 0.01 expected counts (drift-∆R_phi space) + contour=FV_poly(dt) + #===CALCULATE dR_DPHI + # Define azimuthal slices + n_phi_slices = 12 + phi_slices = np.linspace(-np.pi, np.pi, n_phi_slices + 1) + np.pi/4 + phi_slices[phi_slices > np.pi] -= 2*np.pi + + # Calculate event radii and angles, then mask them according to each slice + R = np.sqrt(x**2 + y**2) + phi = np.arctan2(y, x) + phi_cuts = [((phi >= phi_slices[i]) & (phi < phi_slices[i + 1])) if not (phi_slices[i] > phi_slices[i + 1]) else \ + ~((phi <= phi_slices[i]) & (phi > phi_slices[i + 1])) for i in range(n_phi_slices)] + + # Calculate dR_phi by replacing relevant points in loops over phi slices, then return + dR_phi = np.zeros(len(R)) + for i, p in enumerate(phi_cuts): + dR_phi[p] = R[p] - phi_walls[i](dt[p]) + + # Segment events by whether or not they're in the expandable part + expandable = (dt > 71) & (dt < 900) + + # Get the radial cut as a mask between two parts + expansion = 0 + mask = ((dR_phi < (contour + expansion)) & expandable) | ((dR_phi < contour) & ~expandable) + + #cut the drift time + dt_cut = (dt > 71) & (dt < 1030) + + return dt_cut&mask&(dR_phi<=0) + +def WS2024_resistor_XY_cut(x,y): + """ + Cut around hot (radioacitvity not like josh) resistors + x: x position [cm] + y: y position [cm] + """ + res1X = -69.8 + res1Y = 3.5 + res1R = 6 + res2X = -67.5 + res2Y = -14.3 + res2R = 6 + + not_inside_res1 = np.where(np.sqrt( (x-res1X)*(x-res1X) + (y-res1Y)*(y-res1Y) ) < res1R, 0., 1.) + not_inside_res2 = np.where(np.sqrt( (x-res2X)*(x-res2X) + (y-res2Y)*(y-res2Y) ) < res2R, 0., 1.) + + return not_inside_res1 * not_inside_res2 \ No newline at end of file diff --git a/flamedisx/lz/__init__.py b/flamedisx/lz/__init__.py index 96a373ee2..34c34fda0 100644 --- a/flamedisx/lz/__init__.py +++ b/flamedisx/lz/__init__.py @@ -1 +1,2 @@ from .lz import * +from .lz_WS2024 import * \ No newline at end of file diff --git a/flamedisx/lz/data.py b/flamedisx/lz/data.py index 84bfce2cb..cc65c7b8f 100644 --- a/flamedisx/lz/data.py +++ b/flamedisx/lz/data.py @@ -3,43 +3,50 @@ import getpass import os import os.path +import random +import string import flamedisx as fd export, __all__ = fd.exporter() # Path to the root folder of the private LZ data repository, if you have already cloned it -PATH = './lz_private_data' +PATH = 'lz_private_data' def ensure_token(token=None): """Requests for token if token is not already available.""" if token is None: - print(" - Create a token with full 'repo' permissions at " - "https://github.com/settings/tokens\n" + print(" - Create a token with 'read repository' permissions at " + "https://gitlab.com/-/profile/personal_access_tokens\n" " - Save or write it down somewhere \n" " - Type it in the prompt above\n\n" "'Repository not found' means you didn't give the token" - " full 'repo' permissions.\n" - "'Authentication failed' means you mistyped the token.\n") + " the correct permissions.\n" + "'Authentication failed' means you mistyped your username or the" + " token.\n") # We could prompt for username/password instead, but then the password # would be printed in plaintext if there is any problem during cloning. - token = getpass.getpass('Github OAuth token: ') - return token + user = input('Gitlab username:') + token = getpass.getpass('GitLab token: ') + return user, token -def ensure_repo(repo_name, repo_path, token=None): +def ensure_repo(repo_name, repo_path, user=None, token=None): """Clones private repository (prompting for credentials) if we do not have it""" if not os.path.exists(repo_path): print("Private data requested, we must clone repository folder.") - token = ensure_token() - fd.run_command(f'git clone https://{token}:x-oauth-basic' - f'@github.com/{repo_name}') + user, token = ensure_token() + temp_folder = ''.join(random.choices(string.ascii_lowercase, k=8)) + fd.run_command(f'git clone https://{user}:{token}' + f'@gitlab.com/{repo_name} {temp_folder}') + fd.run_command(f'mv {temp_folder}/{repo_path} .') + fd.run_command(f'rm -r -f {temp_folder}') @export def get_lz_file(data_file_name): """Return information from file in lz_private_data/... """ - ensure_repo('robertsjames/lz_private_data.git', PATH) + ensure_repo('luxzeplin/stats/LZFlameFit.git', PATH) return fd.get_resource(f'{PATH}/{data_file_name}') diff --git a/flamedisx/lz/lz.py b/flamedisx/lz/lz.py index 851719684..8caecf148 100644 --- a/flamedisx/lz/lz.py +++ b/flamedisx/lz/lz.py @@ -133,12 +133,12 @@ def __init__(self, *args, ignore_maps_acc=False, cap_upper_cs1=False, **kwargs): self.log10_cs2_acc_domain = None @staticmethod - def photon_detection_eff(z, *, g1=0.113569): - return g1 * tf.ones_like(z) + def photon_detection_eff(drift_time, *, g1=0.113569): + return g1 * tf.ones_like(drift_time) @staticmethod - def s2_photon_detection_eff(z, *, g1_gas=0.092103545): - return g1_gas * tf.ones_like(z) + def s2_photon_detection_eff(drift_time, *, g1_gas=0.092103545): + return g1_gas * tf.ones_like(drift_time) @staticmethod def get_elife(event_time): diff --git a/flamedisx/lz/lz_WS2024.py b/flamedisx/lz/lz_WS2024.py new file mode 100644 index 000000000..651b83a5a --- /dev/null +++ b/flamedisx/lz/lz_WS2024.py @@ -0,0 +1,799 @@ +"""lz detector implementation + +""" +import numpy as np +import math as m +import tensorflow as tf + +import configparser +import os +import pandas as pd + +import flamedisx as fd +from .. import nest as fd_nest +from .WS2024_cuts_and_acceptances import * +import pickle as pkl + +from multihist import Histdd +from scipy import interpolate + +export, __all__ = fd.exporter() +pi = tf.constant(m.pi) +GAS_CONSTANT = 8.314 +N_AVAGADRO = 6.0221409e23 +A_XENON = 131.293 +XENON_REF_DENSITY = 2.90 + +## +# Useful functions +## + + +def interpolate_acceptance(arg, domain, acceptances): + """ Function to interpolate signal acceptance curves + :param arg: argument values for domain interpolation + :param domain: domain values from interpolation map + :param acceptances: acceptance values from interpolation map + :return: Tensor of interpolated map values (same shape as x) + """ + return np.interp(x=arg, xp=domain, fp=acceptances) + +def build_position_map_from_data(map_file, axis_names, bins): + """ + """ + map_df= fd.get_lz_file(map_file) + assert isinstance(map_df, pd.DataFrame), 'Must pass in a dataframe to build position map hisotgram' + + mh = Histdd(bins=bins, axis_names=axis_names) + + add_args = [] + for axis_name in axis_names: + add_args.append(map_df[axis_name].values) + + try: + weights = map_df['weight'].values + except Exception: + weights = None + + mh.add(*add_args, weights=weights) + + return mh + + +## +# Flamedisx sources +## + +## +# Common to all LZ sources +## + + +class LZWS2024Source: + path_s1_corr_latest = 'WS2024/s1Area_Correction_TPC_WS2024_radon_31Jan2024.json' + path_s2_corr_latest = 'WS2024/s2Area_Correction_TPC_WS2024_radon_31Jan2024.json' + + path_s1_acc_curve = 'WS2024/cS1_tritium_acceptance_curve.json' + path_s2_splitting_curve='WS2024/WS2024_S2splittingReconEff_mean.pkl' + path_drift_map_dt='WS2024/drift_map_dt_WS2024.json' + path_drift_map_x='WS2024/drift_map_x_WS2024.json' + + def __init__(self, *args, ignore_LCE_maps=False, ignore_acc_maps=False,ignore_all_cuts=False, ignore_drift_map=False, cap_upper_cs1=False, **kwargs): + super().__init__(*args, **kwargs) + + self.cap_upper_cs1 = cap_upper_cs1 + self.ignore_all_cuts = ignore_all_cuts + assert kwargs['detector'] in ('lz_WS2024',) + + assert os.path.exists(os.path.join( + os.path.dirname(__file__), '../nest/config/', kwargs['detector'] + '.ini')) + + config = configparser.ConfigParser(inline_comment_prefixes=';') + config.read(os.path.join(os.path.dirname(__file__), '../nest/config/', + kwargs['detector'] + '.ini')) + + self.cS1_min = config.getfloat('NEST', 'cS1_min_config') * (1 + self.double_pe_fraction) # phd to phe + self.cS1_max = config.getfloat('NEST', 'cS1_max_config') * (1 + self.double_pe_fraction) # phd to phe + self.S2_min = config.getfloat('NEST', 'S2_min_config') * (1 + self.double_pe_fraction) # phd to phe + self.cS2_max = config.getfloat('NEST', 'cS2_max_config') * (1 + self.double_pe_fraction) # phd to phe + + if not ignore_drift_map: + try: + drift_map=fd.get_lz_file(self.path_drift_map_dt) + self.drift_map_dt = interpolate.LinearNDInterpolator(drift_map['coordinate_system'],drift_map['map'],fill_value=0) + drift_map_xy=fd.get_lz_file(self.path_drift_map_x) + self.drift_map_x = interpolate.LinearNDInterpolator(drift_map_xy['coordinate_system'],drift_map_xy['map'],fill_value=0) + except: + self.drift_map_dt = None + self.drift_map_x = None + + print('Could not load drift maps \n !Using default NEST Calculation!') + else: + print("Ignoring drift map") + + + + self.ignore_acceptances_maps=False + if ignore_acc_maps: + print("ignoring acceptances") + self.ignore_acceptances_maps = True + + self.cs1_acc_domain = None + self.cS2_drift_acceptance_hist = None + else: + try: + df_S1_acc = fd.get_lz_file(self.path_s1_acc_curve) + self.cs1_acc_domain = np.array(df_S1_acc['cS1_phd']) * (1 + self.double_pe_fraction) # phd to phe + self.cs1_acc_curve = np.array(df_S1_acc['cS1_acceptance']) + #TO-DO: Adapt to json for get_lz_file + self.cS2_drift_acceptance_hist= fd.get_lz_file(self.path_s2_splitting_curve) + except Exception: + print("Could not load acceptance curves; setting to 1") + + self.cs1_acc_domain = None + self.cS2_drift_acceptance_hist = None + + if ignore_LCE_maps: + print("ingoring LCE maps") + self.s1_map_latest = None + self.s2_map_latest = None + else: + try: + self.s1_map_latest = fd.InterpolatingMap(fd.get_lz_file(self.path_s1_corr_latest)) + self.s2_map_latest = fd.InterpolatingMap(fd.get_lz_file(self.path_s2_corr_latest)) + except Exception: + print("Could not load maps; setting position corrections to 1") + self.s1_map_latest = None + self.s2_map_latest = None + + @staticmethod + def photon_detection_eff(drift_time, *, g1=0.1122): + """ + g1_gas: Floatable (defined in function argument) + ensure this default is correct + """ + return g1 * tf.ones_like(drift_time) + + @staticmethod + def s2_photon_detection_eff(drift_time, *, g1_gas=0.076404): + """ + g1_gas: Floatable (defined in function argument) + ensure this default is correct + """ + return g1_gas * tf.ones_like(drift_time) + + @staticmethod + def get_elife(event_time): + """ + Static method so just need to know + Rob: Any reason I can't have self? + """ + return 9e6*np.ones_like(event_time) + + def electron_detection_eff(self, drift_time, electron_lifetime): + return self.extraction_eff * tf.exp(-drift_time / electron_lifetime) + + @staticmethod + def s1_posDependence(s1_pos_corr_latest): + return s1_pos_corr_latest + + @staticmethod + def s2_posDependence(s2_pos_corr_latest): + return s2_pos_corr_latest + + def s1_acceptance(self, s1, cs1, cs1_acc_curve): + if self.ignore_all_cuts: + return tf.ones_like(s1, dtype=fd.float_type()) + acceptance = tf.where((s1 >= self.spe_thr) & + (cs1 >= self.cS1_min) & (cs1 <= self.cS1_max), + tf.ones_like(s1, dtype=fd.float_type()), # if condition non-zero + tf.zeros_like(s1, dtype=fd.float_type())) # if false + + # multiplying by efficiency curve + if not self.ignore_acceptances_maps: + acceptance *= cs1_acc_curve + + return acceptance + + def s2_acceptance(self, s2, cs2, cs2_acc_curve, + fv_acceptance, resistor_acceptance, timestamp_acceptance): + if self.ignore_all_cuts: + return tf.ones_like(s2, dtype=fd.float_type()) + acceptance = tf.where((s2 >= self.s2_thr) & + (s2 >= self.S2_min) & (cs2 <= self.cS2_max), + tf.ones_like(s2, dtype=fd.float_type()), # if condition non-zero + tf.zeros_like(s2, dtype=fd.float_type())) # if false + + # multiplying by efficiency curve + if not self.ignore_acceptances_maps: + acceptance *= cs2_acc_curve + + # We will insert the FV acceptance here + acceptance *= fv_acceptance + # We will insert the resistor acceptance here + acceptance *= resistor_acceptance + # We will insert the timestamp acceptance here + acceptance *= timestamp_acceptance + + return acceptance + + def add_extra_columns(self, d): + super().add_extra_columns(d) + + + if 'x_obs' not in d.columns: + x_obs,y_obs=self.derive_observed_xy(d) + d['x_obs'] = x_obs + d['y_obs'] = y_obs + + if (self.s1_map_latest is not None) and (self.s2_map_latest is not None): + #LZLAMA uses correctedX and Y + #I think this is meant to represent cluster (and therfore True position) + d['s1_pos_corr_latest'] = self.s1_map_latest( + np.transpose([d['x'].values, + d['y'].values, + d['drift_time'].values * 1e-9 / 1e-6])) + d['s2_pos_corr_latest'] = self.s2_map_latest( + np.transpose([d['x'].values, + d['y'].values])) + else: + d['s1_pos_corr_latest'] = np.ones_like(d['x'].values) + d['s2_pos_corr_latest'] = np.ones_like(d['x'].values) + + if 'event_time' in d.columns and 'electron_lifetime' not in d.columns: + d['electron_lifetime'] = self.get_elife(d['event_time'].values) + + if 's1' in d.columns and 'cs1' not in d.columns: + d['cs1'] = d['s1'] / d['s1_pos_corr_latest'] + d['cs1_phd'] = d['cs1'] / (1 + self.double_pe_fraction) + if self.cap_upper_cs1 == True: + d['cs1'] = np.where(d['cs1'].values <= self.cS1_max, d['cs1'].values, self.cS1_max) + if 's2' in d.columns and 'cs2' not in d.columns: + d['cs2'] = ( + d['s2'] + / d['s2_pos_corr_latest'] + * np.exp(d['drift_time'] / d['electron_lifetime'])) + d['log10_cs2_phd'] = np.log10(d['cs2'] / (1 + self.double_pe_fraction)) + + if 'cs1' in d.columns and 's1' not in d.columns: + d['s1'] = d['cs1'] * d['s1_pos_corr_latest'] + if 'cs2' in d.columns and 's2' not in d.columns: + d['s2'] = ( + d['cs2'] + * d['s2_pos_corr_latest'] + / np.exp(d['drift_time'] / d['electron_lifetime'])) + + if 'cs1' in d.columns and 'cs2' in d.columns and 'ces_er_equivalent' not in d.columns: + g1 = self.photon_detection_eff(0.) + g1_gas = self.s2_photon_detection_eff(0.) + g2 = fd_nest.calculate_g2(self.gas_field, self.density_gas, self.gas_gap, + g1_gas, self.extraction_eff) + d['ces_er_equivalent'] = (d['cs1'] / fd.tf_to_np(g1) + d['cs2'] / fd.tf_to_np(g2)) * self.Wq_keV + + if 'cs1' in d.columns and 'cs1_acc_curve' not in d.columns: + if self.cs1_acc_domain is not None: + d['cs1_acc_curve'] = interpolate_acceptance( + d['cs1'].values, + self.cs1_acc_domain, + self.cs1_acc_curve) + else: + d['cs1_acc_curve'] = np.ones_like(d['cs1'].values) + if 'cs2' in d.columns and 'cs2_acc_curve' not in d.columns: + if self.cS2_drift_acceptance_hist is not None: + d['cs2_acc_curve'] = WS2024_S2splitting_reconstruction_efficiency( + d['cs2'].values/(1+self.double_pe_fraction),#phe->phd + d['drift_time'].values/1e3,#ns->us + self.cS2_drift_acceptance_hist) + d['cs2_acc_curve'] *=WS2024_trigger_acceptance(d['s2'].values/(1+self.double_pe_fraction)) + + else: + d['cs2_acc_curve'] = np.ones_like(d['cs2'].values) + + if 'fv_acceptance' not in d.columns: + x = d['x_obs'].values + y = d['y_obs'].values + dt=d['drift_time'].values/1e3 + d['fv_acceptance']=WS2024_fiducial_volume_cut(x,y,dt) + + if 'resistor_acceptance' not in d.columns: + x = d['x_obs'].values + y = d['y_obs'].values + d['resistor_acceptance'] = WS2024_resistor_XY_cut(x,y) + + if 'timestamp_acceptance' not in d.columns: + d['timestamp_acceptance'] = np.ones_like(d['event_time'],dtype=bool) + +## +# Different interaction types: flat spectra +## + +GAS_CONSTANT = 8.314 +N_AVAGADRO = 6.0221409e23 +A_XENON = 131.293 +XENON_REF_DENSITY = 2.90 + +@export +class LZ24ERSource(LZWS2024Source, fd.nest.nestERSource): + def __init__(self, *args, **kwargs): + if ('detector' not in kwargs): + kwargs['detector'] = 'lz_WS2024' + super().__init__(*args, **kwargs) + def mean_yield_electron(self, energy): + """ + Update the mean yields to WS2024 LZLAMA (!397) + ERYieldParams from NEST/LZLAMA + Constants are direct over-rides of eqn 6 in Arxiv: 2211.10726 + Energy: energy in keV + """ + m1=12.4886 + m2=85.0 + m3=0.6050 + m4= 2.14687 + m5=25.721 + m6=0. + m7=59.651 + m8=3.6869 + m9=0.2872 + m10=0.1121 + + Nq = energy / self.Wq_keV #equation is in keV + + Qy = m1 + (m2 - m1) / pow((1. + pow(energy /m3,m4)),m9) + \ + m5 + (m6 - m5) / pow((1. + pow(energy /m7, m8)), m10) + + coeff_TI = tf.cast(pow(1. / XENON_REF_DENSITY, 0.3), fd.float_type()) + coeff_Ni = tf.cast(pow(1. / XENON_REF_DENSITY, 1.4), fd.float_type()) + coeff_OL = tf.cast(pow(1. / XENON_REF_DENSITY, -1.7) / + fd.tf_log10(1. + coeff_TI * coeff_Ni * pow(XENON_REF_DENSITY, 1.7)), fd.float_type()) + + Qy *= coeff_OL * fd.tf_log10(1. + coeff_TI * coeff_Ni * pow(self.density, 1.7)) * pow(self.density, -1.7) + + nel_temp = Qy * energy + # Don't let number of electrons go negative + nel = tf.where(nel_temp < 0, + 0 * nel_temp, + nel_temp) + + return nel + def fano_factor(self, nq_mean): + """ + Update fano factor + ERNRWidthParams from NEST/LZLAMA + WS2024 directly over-ride and lose functional form of Arix:2211.10726 + nq_mean: mean number of quanta (self.mean_yield_quanta) + """ + er_free_a = 0.3 + return tf.constant(er_free_a,tf.float32) + + def variance(self, *args): + """ + Update variance + ERNRWidthParams from NEST/LZLAMA + Variance of skew guessian. I don't understand the args thing. + """ + nel_mean = args[0] + nq_mean = args[1] + recomb_p = args[2] + ni = args[3] + + if self.detector in ['lz','lz_WS2024']: + er_free_b = 0.04311 + else: + er_free_b = 0.0553 + er_free_c = 0.15505 + er_free_d = 0.46894 + er_free_e = -0.26564 + + elec_frac = nel_mean / nq_mean + ampl = er_free_b + + ampl = tf.cast(0.086036 + (er_free_b - 0.086036) / + pow((1. + pow(self.drift_field / 295.2, 251.6)), 0.0069114), + fd.float_type()) + wide = er_free_c + cntr = er_free_d + skew = er_free_e + + mode = cntr + 2. / (tf.sqrt(2. * pi)) * skew * wide / tf.sqrt(1. + skew * skew) + norm = 1. / (tf.exp(-0.5 * pow(mode - cntr, 2.) / (wide * wide)) * + (1. + tf.math.erf(skew * (mode - cntr) / (wide * tf.sqrt(2.))))) + + omega = norm * ampl * tf.exp(-0.5 * pow(elec_frac - cntr, 2.) / (wide * wide)) * \ + (1. + tf.math.erf(skew * (elec_frac - cntr) / (wide * tf.sqrt(2.)))) + omega = tf.where(nq_mean == 0, + tf.zeros_like(omega, dtype=fd.float_type()), + omega) + + return recomb_p * (1. - recomb_p) * ni + omega * omega * ni * ni + + +@export +class LZ24GammaSource(LZWS2024Source, fd.nest.nestGammaSource): + def __init__(self, *args, **kwargs): + if ('detector' not in kwargs): + kwargs['detector'] = 'lz_WS2024' + super().__init__(*args, **kwargs) + + +@export +class LZ24ERGammaWeightedSource(LZWS2024Source, fd.nest.nestERGammaWeightedSource): + def __init__(self, *args, **kwargs): + if ('detector' not in kwargs): + kwargs['detector'] = 'lz_WS2024' + super().__init__(*args, **kwargs) + + +@export +class LZ24NRSource(LZWS2024Source, fd.nest.nestNRSource): + def __init__(self, *args, **kwargs): + if ('detector' not in kwargs): + kwargs['detector'] = 'lz_WS2024' + super().__init__(*args, **kwargs) + def mean_yields(self, energy): + """ + Update the mean yields to WS2024 LZLAMA (!397) + NRYieldParams from NEST/LZLAMA + See section C. in Arxiv: 2211.10726 + Energy: energy in keV + """ + nr_nuis_alpha = 10.19 + nr_nuis_beta = 1.11 + nr_nuis_gamma = 0.0498 + nr_nuis_delta = -0.0533 + nr_nuis_epsilon = 12.46 + nr_nuis_zeta = 0.2942 + nr_nuis_eta = 1.899 + nr_nuis_theta = 0.3197 + nr_nuis_l = 2.066 + nr_nuis_p = 0.509 + nr_new_nuis_a = 0.996 + nr_new_nuis_b = 0.999 + + TIB = nr_nuis_gamma * tf.math.pow(self.drift_field, nr_nuis_delta) * pow(self.density / XENON_REF_DENSITY, 0.3) + Qy = 1. / (TIB * tf.math.pow(energy + nr_nuis_epsilon, nr_nuis_p)) + Qy *= (1. - (1. / tf.math.pow(1. + tf.math.pow(tf.math.divide_no_nan(energy , nr_nuis_zeta), nr_nuis_eta), nr_new_nuis_a))) + + nel_temp = Qy * energy + # Don't let number of electrons go negative + nel = tf.where(nel_temp < 0, + 0 * nel_temp, + nel_temp) + + nq_temp = nr_nuis_alpha * pow(energy, nr_nuis_beta) + + nph_temp = (nq_temp - nel) * (1. - (1. / tf.math.pow(1. + tf.math.pow(tf.math.divide_no_nan(energy , nr_nuis_theta), nr_nuis_l), nr_new_nuis_b))) + # Don't let number of photons go negative + nph = tf.where(nph_temp < 0, + 0 * nph_temp, + nph_temp) + + nq = nel + nph + + ni = (4. / TIB) * (tf.exp(nel * TIB / 4.) - 1.) + + nex = nq - ni + + ex_ratio = tf.cast(tf.math.divide_no_nan(nex , ni), fd.float_type()) + + ex_ratio = tf.where(tf.logical_and(ex_ratio < self.alpha, energy > 100.), + self.alpha * tf.ones_like(ex_ratio, dtype=fd.float_type()), + ex_ratio) + ex_ratio = tf.where(tf.logical_and(ex_ratio > 1., energy < 1.), + tf.ones_like(ex_ratio, dtype=fd.float_type()), + ex_ratio) + ex_ratio = tf.where(tf.math.is_nan(ex_ratio), + tf.zeros_like(ex_ratio, dtype=fd.float_type()), + ex_ratio) + + return nel, nq, ex_ratio + + def yield_fano(self, nq_mean): + """ + Update fano factor + ERNRWidthParams from NEST/LZLAMA + nq_mean: mean number of quanta (self.mean_yield_quanta) + """ + if self.detector in ['lz','lz_WS2024']: + nr_free_a = 0.404 + nr_free_b = 0.393 + else: + nr_free_a = 1. + nr_free_b = 1. + + ni_fano = tf.ones_like(nq_mean, dtype=fd.float_type()) * nr_free_a + nex_fano = tf.ones_like(nq_mean, dtype=fd.float_type()) * nr_free_b + + return ni_fano, nex_fano + + @staticmethod + def skewness(nq_mean): + """ + Update skewness + ERNRWidthParams from NEST/LZLAMA + nq_mean: mean number of quanta (self.mean_yield_quanta) + """ + nr_free_f = 2.220 + + mask = tf.less(nq_mean, 1e4 * tf.ones_like(nq_mean)) + skewness = tf.ones_like(nq_mean, dtype=fd.float_type()) * nr_free_f + skewness_masked = tf.multiply(skewness, tf.cast(mask, fd.float_type())) + + return skewness_masked + + def variance(self, *args): + """ + Update Variance + ERNRWidthParams from NEST/LZLAMA + I don't understand args + """ + nel_mean = args[0] + nq_mean = args[1] + recomb_p = args[2] + ni = args[3] + + if self.detector in ['lz','lz_WS2024']: + nr_free_c = 0.0383 + else: + nr_free_c = 0.1 + nr_free_d = 0.497 + nr_free_e = 0.1906 + + elec_frac = nel_mean / nq_mean + + omega = nr_free_c * tf.exp(-0.5 * pow(elec_frac - nr_free_d, 2.) / (nr_free_e * nr_free_e)) + omega = tf.where(nq_mean == 0, + tf.zeros_like(omega, dtype=fd.float_type()), + tf.cast(omega, dtype=fd.float_type())) + + return recomb_p * (1. - recomb_p) * ni + omega * omega * ni * ni + +## +# This is an architectural problem! +# can't specific changes like this, would need to make a whole new NEST! +# Maybe I can play with the inheritance +## +## +# Calibration sources +## + + +@export +class LZ24CH3TSource(LZ24ERSource, fd.nest.CH3TSource): + t_start = pd.to_datetime('2022-04-19T00:00:00') + t_start = t_start.tz_localize(tz='America/Denver') + + t_stop = pd.to_datetime('2022-04-19T00:00:00') + t_stop = t_stop.tz_localize(tz='America/Denver') + + def __init__(self, *args, **kwargs): + if ('detector' not in kwargs): + kwargs['detector'] = 'lz_WS2024' + super().__init__(*args, **kwargs) + + +@export +class LZ24DDSource(LZ24NRSource, fd.nest.DDSource): + t_start = pd.to_datetime('2022-04-19T00:00:00') + t_start = t_start.tz_localize(tz='America/Denver') + + t_stop = pd.to_datetime('2022-04-19T00:00:00') + t_stop = t_stop.tz_localize(tz='America/Denver') + + def __init__(self, *args, **kwargs): + if ('detector' not in kwargs): + kwargs['detector'] = 'lz_WS2024' + super().__init__(*args, **kwargs) + + +## +# Signal sources +## + + +@export +class LZWIMPSource(LZ24NRSource, fd.nest.nestWIMPSource): + def __init__(self, *args, **kwargs): + if ('detector' not in kwargs): + kwargs['detector'] = 'lz_WS2024' + super().__init__(*args, **kwargs) + + +@export +class LZ24FermionicDMSource(LZ24ERSource, fd.nest.FermionicDMSource): + def __init__(self, *args, **kwargs): + if ('detector' not in kwargs): + kwargs['detector'] = 'lz_WS2024' + super().__init__(*args, **kwargs) + + +## +# Background sources +## + + +@export +class LZ24Pb214Source(LZ24ERSource, fd.nest.Pb214Source, fd.nest.nestSpatialRateERSource): + def __init__(self, *args, bins=None, **kwargs): + if ('detector' not in kwargs): + kwargs['detector'] = 'lz_WS2024' + + if bins is None: + bins=(np.sqrt(np.linspace(0.**2, 67.8**2, num=21)), + np.linspace(86000., 936500., num=21)) + + mh = build_position_map_from_data('Pb214_spatial_map_data.pkl', ['r', 'drift_time'], bins) + self.spatial_hist = mh + + super().__init__(*args, **kwargs) + + +@export +class LZ24DetERSource(LZ24ERSource, fd.nest.DetERSource, fd.nest.nestSpatialRateERSource): + def __init__(self, *args, **kwargs): + if ('detector' not in kwargs): + kwargs['detector'] = 'lz_WS2024' + + mh = fd.get_lz_file('DetER_spatial_map_hist.pkl') + self.spatial_hist = mh + + super().__init__(*args, **kwargs) + + +@export +class LZ24BetaSource(LZ24ERSource, fd.nest.BetaSource): + def __init__(self, *args, **kwargs): + if ('detector' not in kwargs): + kwargs['detector'] = 'lz_WS2024' + super().__init__(*args, **kwargs) + + +@export +class LZ24Xe136Source(LZ24ERSource, fd.nest.Xe136Source): + def __init__(self, *args, **kwargs): + if ('detector' not in kwargs): + kwargs['detector'] = 'lz_WS2024' + super().__init__(*args, **kwargs) + + +@export +class LZ24vERSource(LZ24ERSource, fd.nest.vERSource, fd.nest.nestTemporalRateOscillationERSource): + def __init__(self, *args, amplitude=None, phase_ns=None, period_ns=None, **kwargs): + if ('detector' not in kwargs): + kwargs['detector'] = 'lz_WS2024' + + if amplitude is None: + self.amplitude = 2. * 0.01671 + else: + self.amplitude = amplitude + + if phase_ns is None: + self.phase_ns = pd.to_datetime('2022-01-04T00:00:00').value + else: + self.phase_ns = phase_ns + + if period_ns is None: + self.period_ns = 1. * 3600. * 24. * 365.25 * 1e9 + else: + self.period_ns = period_ns + + super().__init__(*args, **kwargs) + + +@export +class LZ24Ar37Source(LZ24ERSource, fd.nest.Ar37Source, fd.nest.nestTemporalRateDecayERSource): + def __init__(self, *args, time_constant_ns=None, **kwargs): + if ('detector' not in kwargs): + kwargs['detector'] = 'lz_WS2024' + + if time_constant_ns is None: + self.time_constant_ns = (35.0 / np.log(2)) * 1e9 * 3600. * 24. + else: + self.time_constant_ns = time_constant_ns + + super().__init__(*args, **kwargs) + + +@export +class LZXe124Source(LZWS2024Source, fd.nest.Xe124Source): + def __init__(self, *args, **kwargs): + if ('detector' not in kwargs): + kwargs['detector'] = 'lz_WS2024' + super().__init__(*args, **kwargs) + + +@export +class LZXe127Source(LZWS2024Source, fd.nest.Xe127Source, fd.nest.nestSpatialTemporalRateDecayERSource): + def __init__(self, *args, bins=None, time_constant_ns=None, **kwargs): + if ('detector' not in kwargs): + kwargs['detector'] = 'lz_WS2024' + + if bins is None: + bins=(np.sqrt(np.linspace(0.**2, 67.8**2, num=51)), + np.linspace(fd.lz.LZERSource().z_bottom, fd.lz.LZERSource().z_top, num=51)) + + mh = fd.get_lz_file('Xe127_spatial_map_hist.pkl') + self.spatial_hist = mh + + if time_constant_ns is None: + self.time_constant_ns = (36.4 / np.log(2)) * 1e9 * 3600. * 24. + else: + self.time_constant_ns = time_constant_ns + + super().__init__(*args, **kwargs) + + +@export +class LZ24B8Source(LZ24NRSource, fd.nest.B8Source, fd.nest.nestTemporalRateOscillationNRSource): + def __init__(self, *args, amplitude=None, phase_ns=None, period_ns=None, **kwargs): + if ('detector' not in kwargs): + kwargs['detector'] = 'lz_WS2024' + + if amplitude is None: + self.amplitude = 2. * 0.01671 + else: + self.amplitude = amplitude + + if phase_ns is None: + self.phase_ns = pd.to_datetime('2022-01-04T00:00:00').value + else: + self.phase_ns = phase_ns + + if period_ns is None: + self.period_ns = 1. * 3600. * 24. * 365.25 * 1e9 + else: + self.period_ns = period_ns + + super().__init__(*args, **kwargs) + + +@export +class LZ24DetNRSource(LZ24NRSource, fd.nest.nestSpatialRateNRSource): + """ + """ + + def __init__(self, *args, **kwargs): + if ('detector' not in kwargs): + kwargs['detector'] = 'lz_WS2024' + + df_DetNR = fd.get_lz_file('DetNR_spectrum.pkl') + + self.energies = tf.convert_to_tensor(df_DetNR['energy_keV'].values, dtype=fd.float_type()) + self.rates_vs_energy = tf.convert_to_tensor(df_DetNR['spectrum_value_norm'].values, dtype=fd.float_type()) + + mh = fd.get_lz_file('DetNR_spatial_map_hist.pkl') + self.spatial_hist = mh + + super().__init__(*args, **kwargs) + + +## TO DO: +## ADD in accidentals source + + +## +# Source groups +## + + +@export +class LZERSourceGroup(LZWS2024Source, fd.nest.nestERSourceGroup): + def __init__(self, *args, **kwargs): + if ('detector' not in kwargs): + kwargs['detector'] = 'lz_WS2024' + super().__init__(*args, **kwargs) + + +@export +class LZERGammaWeightedSourceGroup(LZWS2024Source, fd.nest.nestERGammaWeightedSourceGroup): + def __init__(self, *args, **kwargs): + if ('detector' not in kwargs): + kwargs['detector'] = 'lz_WS2024' + super().__init__(*args, **kwargs) + + +@export +class LZNRSourceGroup(LZWS2024Source, fd.nest.nestNRSourceGroup): + def __init__(self, *args, **kwargs): + if ('detector' not in kwargs): + kwargs['detector'] = 'lz_WS2024' + super().__init__(*args, **kwargs) diff --git a/flamedisx/nest/config/lz_WS2024.ini b/flamedisx/nest/config/lz_WS2024.ini new file mode 100644 index 000000000..038563d34 --- /dev/null +++ b/flamedisx/nest/config/lz_WS2024.ini @@ -0,0 +1,47 @@ +[NEST] + +; Detector parameters +temperature_config = 175.9 ; K +pressure_config = 1.86 ; bar +drift_field_config = 96.5 ; V/cm +gas_field_config = 8.301178 ; kV/cm + +radius_config = 72.8 ; cm not sure what to do here + +z_topDrift_config = 146.1 ; cm +dt_cntr_config = 517927 ; ns old: 462500. +gate_config = 145.6 ; cm +cathode_config = 0. ; cm + +g1_config = 0.1122 ; MADE FLOATABLE INSTEAD +elife_config = 9000000. ; ns + +spe_res_config = 0.338 +spe_thr_config = 0.1 +spe_eff_config = 0.999 + +num_pmts_config = 481 +double_pe_fraction_config = 0.21392 + +coin_level_config = 3 + +gas_gap_config = 8. ; mm +g1_gas_config = 0.076404 ; MADE FLOATABLE INSTEAD +s2Fano_config = 4.0 +;not sure about this +s2_thr_config = 570.58 ; + +S1_noise_config = 0.083 ; see the note +S2_noise_config = 0. + +; Selection parameters +S1_min_config = 0. ; REDUNDANT FOR THIS SOURCE +S1_max_config = 0. ; REDUNDANT FOR THIS SOURCE +;14.5e-=14.5*SE_size=14.5*48.5=703.25 +S2_min_config = 645.25 ; phd, not phe +S2_max_config = 0. ; REDUNDANT FOR THIS SOURCE + +cS1_min_config = 3. ; phd, not phe +cS1_max_config = 80. ; phd, not phe + +cS2_max_config = 100000. ; phd, not phe diff --git a/flamedisx/nest/lxe_blocks/detection.py b/flamedisx/nest/lxe_blocks/detection.py index 3f7170231..771fe1ba6 100644 --- a/flamedisx/nest/lxe_blocks/detection.py +++ b/flamedisx/nest/lxe_blocks/detection.py @@ -141,7 +141,7 @@ class DetectPhotons(DetectPhotonsOrElectrons): 's1_posDependence') + special_model_functions @staticmethod - def s1_posDependence(r, z): + def s1_posDependence(r): """ Override for specific detector. """ diff --git a/flamedisx/nest/lxe_blocks/energy_spectrum.py b/flamedisx/nest/lxe_blocks/energy_spectrum.py index ffcc186b6..480883394 100644 --- a/flamedisx/nest/lxe_blocks/energy_spectrum.py +++ b/flamedisx/nest/lxe_blocks/energy_spectrum.py @@ -10,7 +10,6 @@ export, __all__ = fd.exporter() o = tf.newaxis - @export class EnergySpectrum(fd.FirstBlock): dimensions = ('energy',) @@ -18,7 +17,8 @@ class EnergySpectrum(fd.FirstBlock): 'energies','max_dim_size', 'radius', 'z_top', 'z_bottom', 'z_topDrift', 'drift_velocity', - 't_start', 't_stop') + 't_start', 't_stop', + 'drift_map_dt', 'drift_map_x') # The default boundaries are at points where the WIMP wind is at its # average speed. @@ -32,13 +32,41 @@ class EnergySpectrum(fd.FirstBlock): # Just a dummy 0-10 keV spectrum energies = tf.cast(tf.linspace(0., 10., 1000), dtype=fd.float_type()) - default_size=100 max_dim_size : dict def setup(self): + """ + Setups up a managable max dim size + """ assert isinstance(self.max_dim_size,dict) self.max_dim_size={'energy':self.max_dim_size['energy']} - + + drift_map_dt = None + def derive_drift_time(self, data): + """ + Helper function for getting drift time from true z + data: Data DataFrame + """ + if self.drift_map_dt is None: + return (self.z_topDrift - data['z']) / self.drift_velocity + + return self.drift_map_dt(np.array([data['r'], data['z']]).T) + + drift_map_x = None + def derive_observed_xy(self, data): + """ + """ + if self.drift_map_x is None: + return data['x'], data['y'] + + cosTheta = data['x'] / data['r'] + sinTheta = data['y'] / data['r'] + + x_obs = self.drift_map_x(np.array([data['r'], data['z']]).T) * cosTheta + y_obs = self.drift_map_x(np.array([data['r'], data['z']]).T) * sinTheta + + return x_obs, y_obs + def domain(self, data_tensor): assert isinstance(self.energies, tf.Tensor) # see WIMPsource for why @@ -116,11 +144,16 @@ def draw_positions(self, n_events, **params): data = dict() data['r'] = (np.random.rand(n_events) * self.radius**2)**0.5 data['theta'] = np.random.uniform(0, 2*np.pi, size=n_events) + data['x'], data['y'] = fd.pol_to_cart(data['r'], data['theta']) data['z'] = np.random.uniform(self.z_bottom, self.z_top, size=n_events) - data['x'], data['y'] = fd.pol_to_cart(data['r'], data['theta']) - data['drift_time'] = (self.z_topDrift-data['z']) / self.drift_velocity + data['x_obs'], data['y_obs'] = self.derive_observed_xy(data) + data['r_obs'], data['theta_obs'] = fd.cart_to_pol(data['x_obs'], data['y_obs']) + + data['drift_time'] = self.derive_drift_time(data) + data.pop('z') + return data def draw_time(self, n_events, **params): @@ -167,7 +200,7 @@ def validate_fix_truth(self, d): """ # When passing in an event as DataFrame we select and set # only these columns: - cols = ['x', 'y', 'z', 'r', 'theta', 'event_time', 'drift_time'] + cols = ['x', 'y', 'r', 'theta', 'event_time', 'drift_time'] if d is None: return dict() elif isinstance(d, pd.DataFrame): @@ -186,21 +219,20 @@ def validate_fix_truth(self, d): assert isinstance(d, dict), \ "fix_truth needs to be a DataFrame, dict, or None" - if 'z' in d: + if 'drift_time' in d: # Position is fixed. Ensure both Cartesian and polar coordinates - # are available, and compute drift_time from z. + # are available. if 'x' in d and 'y' in d: d['r'], d['theta'] = fd.cart_to_pol(d['x'], d['y']) elif 'r' in d and 'theta' in d: d['x'], d['y'] = fd.pol_to_cart(d['r'], d['theta']) else: - raise ValueError("When fixing position, give (x, y, z), " - "or (r, theta, z).") - d['drift_time'] = (self.z_topDrift-d['z']) / self.drift_velocity + raise ValueError("When fixing position, give (x, y, drift_time), " + "or (r, theta, drift_time).") elif 'event_time' not in d and 'energy' not in d: # Neither position, time, nor energy given - raise ValueError(f"Dict should contain at least ['x', 'y', 'z'] " - "and/or ['r', 'theta', 'z'] and/or 'event_time' " + raise ValueError(f"Dict should contain at least ['x', 'y', 'drift_time'] " + "and/or ['r', 'theta', 'drift_time'] and/or 'event_time' " f"and/or 'energy', but it contains: {d.keys()}") return d @@ -287,24 +319,22 @@ class SpatialRateEnergySpectrum(FixedShapeEnergySpectrum): def setup(self): assert isinstance(self.spatial_hist, Histdd) - # Are we Cartesian, modified Cartesian (x, y, dt), polar, 2D (r, z or r, dt), + # Are we Cartesian (x, y, dt), polar (r, theta, dt), 2D (r, dt), # or in trouble? axes = tuple(self.spatial_hist.axis_names) - self.mod_cart = (axes == ('x', 'y', 'drift_time')) - self.polar = (axes == ('r', 'theta', 'z')) - self.r_z = (axes == ('r', 'z')) + self.polar = (axes == ('r', 'theta', 'drift_time')) self.r_dt = (axes == ('r', 'drift_time')) self.bin_volumes = self.spatial_hist.bin_volumes() # Volume element in cylindrical coords = r * (dr dq dz) if self.polar: self.bin_volumes *= self.spatial_hist.bin_centers('r')[:, None, None] - elif (self.r_z or self.r_dt): + elif (self.r_dt): self.bin_volumes *= self.spatial_hist.bin_centers('r')[:, None] else: - assert (axes == ('x', 'y', 'z')) or self.mod_cart, \ + assert (axes == ('x', 'y', 'drift_time')), \ ("axis_names of spatial_rate_hist must be " - "['r', 'theta', 'z'], ['r', 'z'], ['r', 'drift_time'], ['x', 'y', 'z'] " + "['r', 'theta', 'drift_time'], ['r', 'drift_time'] " "or ['x', 'y', 'drift_time']") # Normalize the histogram @@ -318,19 +348,13 @@ def setup(self): (self.spatial_hist.histogram / self.bin_volumes) * self.bin_volumes.sum()) - def energy_spectrum_rate_multiplier(self, x, y, z): + def energy_spectrum_rate_multiplier(self, x, y, drift_time): if self.polar: - positions = list(fd.cart_to_pol(x, y)) + [z] - elif self.r_z: - positions = [fd.cart_to_pol(x, y)[0]] + [z] + positions = list(fd.cart_to_pol(x, y)) + [drift_time] elif self.r_dt: - dt = (self.z_topDrift - z) / self.drift_velocity - positions = [fd.cart_to_pol(x, y)[0]] + [dt] - elif self.mod_cart: - dt = (self.z_topDrift - z) / self.drift_velocity - positions = [x, y, dt] + positions = [fd.cart_to_pol(x, y)[0]] + [drift_time] else: - positions = [x, y, z] + positions = [x, y, drift_time] return self.local_rate_multiplier.lookup(*positions) def draw_positions(self, n_events, **params): @@ -341,21 +365,14 @@ def draw_positions(self, n_events, **params): positions = self.spatial_hist.get_random(size=n_events) for idx, col in enumerate(self.spatial_hist.axis_names): data[col] = positions[:, idx] - if self.mod_cart: - data['z'] = self.z_topDrift - data['drift_time'] * self.drift_velocity if self.polar: data['x'], data['y'] = fd.pol_to_cart(data['r'], data['theta']) - elif self.r_z: - theta = np.random.uniform(0, 2*np.pi, size=n_events) - data['x'], data['y'] = fd.pol_to_cart(data['r'], theta) elif self.r_dt: theta = np.random.uniform(0, 2*np.pi, size=n_events) data['x'], data['y'] = fd.pol_to_cart(data['r'], theta) - data['z'] = self.z_topDrift - data['drift_time'] * self.drift_velocity else: data['r'], data['theta'] = fd.cart_to_pol(data['x'], data['y']) - data['drift_time'] = (self.z_topDrift-data['z']) / self.drift_velocity return data diff --git a/flamedisx/nest/lxe_sources.py b/flamedisx/nest/lxe_sources.py index a8ace77f2..ad78e60e5 100644 --- a/flamedisx/nest/lxe_sources.py +++ b/flamedisx/nest/lxe_sources.py @@ -22,7 +22,9 @@ class nestSource(fd.BlockModelSource): def __init__(self, *args, detector='default', **kwargs): - assert detector in ('default', 'lz') + + assert detector in ('default', 'lz','lz_WS2024') + self.detector = detector assert os.path.exists(os.path.join( @@ -49,10 +51,8 @@ def __init__(self, *args, detector='default', **kwargs): # energy_spectrum.py self.radius = config.getfloat('NEST', 'radius_config') self.z_topDrift = config.getfloat('NEST', 'z_topDrift_config') - self.z_top = self.z_topDrift - self.drift_velocity * \ - config.getfloat('NEST', 'dt_min_config') - self.z_bottom = self.z_topDrift - self.drift_velocity * \ - config.getfloat('NEST', 'dt_max_config') + self.z_top = config.getfloat('NEST', 'gate_config') + self.z_bottom = config.getfloat('NEST', 'cathode_config') # detection.py / pe_detection.py / double_pe.py / final_signals.py self.g1 = config.getfloat('NEST', 'g1_config') @@ -126,32 +126,32 @@ def mu_correction(*args): # detection.py - def photon_detection_eff(self, z): - return self.g1 * tf.ones_like(z) + def photon_detection_eff(self, drift_time): + return self.g1 * tf.ones_like(drift_time) def electron_detection_eff(self, drift_time): return self.extraction_eff * tf.exp(-drift_time / self.elife) - def s2_photon_detection_eff(self, z): - return self.g1_gas * tf.ones_like(z) + def s2_photon_detection_eff(self, drift_time): + return self.g1_gas * tf.ones_like(drift_time) # secondary_quanta_generation.py - def electron_gain_mean(self, z): + def electron_gain_mean(self, drift_time): elYield = ( 0.137 * self.gas_field * 1e3 - 4.70e-18 * (N_AVAGADRO * self.density_gas / A_XENON)) \ * self.gas_gap * 0.1 - return tf.cast(elYield, fd.float_type()) * tf.ones_like(z) + return tf.cast(elYield, fd.float_type()) * tf.ones_like(drift_time) - def electron_gain_std(self, z): + def electron_gain_std(self, drift_time): elYield = ( 0.137 * self.gas_field * 1e3 - 4.70e-18 * (N_AVAGADRO * self.density_gas / A_XENON)) \ * self.gas_gap * 0.1 - return tf.sqrt(self.s2Fano * elYield) * tf.ones_like(z) + return tf.sqrt(self.s2Fano * elYield) * tf.ones_like(drift_time) # pe_detection.py @@ -330,7 +330,7 @@ def skewness(self, nq_mean): skewness = tf.ones_like(nq_mean, dtype=fd.float_type()) * skew skewness_masked = tf.multiply(skewness, tf.cast(mask_product, fd.float_type())) - if self.detector == 'lz': + if self.detector in ['lz','lz_WS2024']: skewness_masked = tf.zeros_like(nq_mean, dtype=fd.float_type()) return skewness_masked