From 6ede17a90818744d49fd659e106b06b344d042f8 Mon Sep 17 00:00:00 2001 From: Lorenzo Principe <28869147+lorenzomag@users.noreply.github.com> Date: Thu, 8 Aug 2024 02:18:33 +0200 Subject: [PATCH 01/22] Initial draft of Migdal source --- flamedisx/nest/lxe_sources.py | 120 +++++++++++++++++++++++++++++++++- 1 file changed, 119 insertions(+), 1 deletion(-) diff --git a/flamedisx/nest/lxe_sources.py b/flamedisx/nest/lxe_sources.py index 628a4a4f..961ae395 100644 --- a/flamedisx/nest/lxe_sources.py +++ b/flamedisx/nest/lxe_sources.py @@ -672,4 +672,122 @@ def get_energy_hist( RATES, [time_bin_edges, energy_bin_edges], axis_names=("time", "energy") ) - return hist \ No newline at end of file + return hist + + +@export +class nestMigdalSource(nestERSource): + + model_blocks = (fd_nest.WIMPEnergySpectrum,) + nestERSource.model_blocks[1:] + + def __init__( + self, + *args, + wimp_mass=40, + sigma=1e-45, + fid_mass=1.0, + min_E=1e-2, + max_E=80.0, + n_energy_bins=800, + min_time="2019-09-01T08:28:00", + max_time="2020-09-01T08:28:00", + n_time_bins=25, + modulation=True, + migdal_model="Cox", + **kwargs + ): + if migdal_model not in ["Ibe", "Cox", "Cox_dipole"]: + raise ValueError("Invalid Migdal model. Choose from 'Ibe', 'Cox', 'Cox_dipole'") + + if "detector" not in kwargs: + kwargs["detector"] = "default" + + self.energy_hist = self.get_energy_hist( + wimp_mass=wimp_mass, + sigma=sigma, + min_E=min_E, + max_E=max_E, + n_energy_bins=n_energy_bins, + min_time=min_time, + max_time=max_time, + n_time_bins=n_time_bins, + modulation=modulation, + migdal_model=migdal_model, + ) + + livetime = pd.Timestamp(max_time) - pd.Timestamp(min_time) + livetime = livetime.days / 365.25 + scale = fid_mass * livetime + self.energy_hist *= scale + + self.n_time_bins = len(self.energy_hist.bin_edges[0]) + e_centers = fd_nest.WIMPEnergySpectrum.bin_centers( + self.energy_hist.bin_edges[1] + ) + self.energies = fd.np_to_tf(e_centers) + + self.array_columns = (("energy_spectrum", len(e_centers)),) + + super().__init__(*args, **kwargs) + + @staticmethod + def get_energy_hist( + wimp_mass=40.0, + sigma=1e-45, + min_E=1e-2, + max_E=40.0, + n_energy_bins=800, + min_time="2019-09-01T08:28:00", + max_time="2020-09-01T08:28:00", + n_time_bins=25, + modulation=True, + migdal_model="Cox", + ): + dipole = False + if migdal_model == "Cox_dipole": + migdal_model = "Cox" + dipole = True + + + energy_bin_edges = np.linspace(min_E, max_E, n_energy_bins + 1) + energy_values = (energy_bin_edges[1:] + energy_bin_edges[:-1]) / 2 + time_bin_edges = ( + pd.date_range(min_time, max_time, periods=n_time_bins + 1).to_julian_date() + - 2451545.0 # Convert to J2000 + ) + times = (time_bin_edges[1:] + time_bin_edges[:-1]) / 2 + + rates_list = [] + if modulation: + # This should probably be done in a more efficient way + for time in times: + rates_list.append( + wr.rate_wimp_std( + energy_values, + mw=wimp_mass, + sigma_nucleon=sigma, + t=time, + detection_mechanism="migdal", + migdal_model=migdal_model, + dipole=dipole, + ) + ) + else: + rates = wr.rate_wimp_std( + energy_values, + mw=wimp_mass, + sigma_nucleon=sigma, + detection_mechanism="migdal", + migdal_model=migdal_model, + dipole=dipole, + ) + for _ in times: + rates_list.append(rates) + + RATES = np.array(rates_list) + + hist = mh.Histdd.from_histogram( + RATES, [time_bin_edges, energy_bin_edges], axis_names=("time", "energy") + ) + + return hist From 56ca91ad5249343163b93d888a853b2188c8bec9 Mon Sep 17 00:00:00 2001 From: Lorenzo Principe <28869147+lorenzomag@users.noreply.github.com> Date: Sun, 11 Aug 2024 22:56:52 +0200 Subject: [PATCH 02/22] Use multiprocessing when modulation on --- flamedisx/nest/lxe_sources.py | 31 ++++++++++++++++++------------- 1 file changed, 18 insertions(+), 13 deletions(-) diff --git a/flamedisx/nest/lxe_sources.py b/flamedisx/nest/lxe_sources.py index 628a4a4f..033e7e14 100644 --- a/flamedisx/nest/lxe_sources.py +++ b/flamedisx/nest/lxe_sources.py @@ -1,3 +1,4 @@ +from concurrent.futures import ProcessPoolExecutor import configparser import math as m import os @@ -642,27 +643,31 @@ def get_energy_hist( energy_values = (energy_bin_edges[1:] + energy_bin_edges[:-1]) / 2 time_bin_edges = ( pd.date_range(min_time, max_time, periods=n_time_bins + 1).to_julian_date() - - 2451545.0 # Convert to J2000 + - 2451545.0 # Convert to J2000 ) times = (time_bin_edges[1:] + time_bin_edges[:-1]) / 2 rates_list = [] if modulation: - for time in times: - rates_list.append( - wr.rate_wimp_std( - energy_values, - mw=wimp_mass, - sigma_nucleon=sigma, - t=time, + with ProcessPoolExecutor(4) as executor: + for time in times: + rates_list.append( + executor.submit( + wr.rate_wimp_std, + energy_values, + mw=wimp_mass, + sigma_nucleon=sigma, + t=time, + ) ) - ) + + rates_list = [future.result() for future in rates_list] else: rates = wr.rate_wimp_std( - energy_values, - mw=wimp_mass, - sigma_nucleon=sigma, - ) + energy_values, + mw=wimp_mass, + sigma_nucleon=sigma, + ) for _ in times: rates_list.append(rates) From ad5ee93d590b1cc53a69062bd129d505b3ef58df Mon Sep 17 00:00:00 2001 From: Lorenzo Principe <28869147+lorenzomag@users.noreply.github.com> Date: Mon, 12 Aug 2024 14:32:43 +0200 Subject: [PATCH 03/22] Fix normalisation bug --- flamedisx/nest/lxe_blocks/energy_spectrum.py | 2 +- flamedisx/nest/lxe_sources.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/flamedisx/nest/lxe_blocks/energy_spectrum.py b/flamedisx/nest/lxe_blocks/energy_spectrum.py index e861a95f..91cb7c69 100644 --- a/flamedisx/nest/lxe_blocks/energy_spectrum.py +++ b/flamedisx/nest/lxe_blocks/energy_spectrum.py @@ -449,7 +449,7 @@ def clip_j2000_times(self, ts): return np.clip(ts, tbins[0], tbins[-1]) def mu_before_efficiencies(self, **params): - return self.energy_hist.n / self.n_time_bins + return self.energy_hist.n def random_truth(self, n_events, fix_truth=None, **params): """Draw n_events random energies and times from the energy/ diff --git a/flamedisx/nest/lxe_sources.py b/flamedisx/nest/lxe_sources.py index dbb3d61b..23a59ee0 100644 --- a/flamedisx/nest/lxe_sources.py +++ b/flamedisx/nest/lxe_sources.py @@ -597,8 +597,8 @@ def __init__( scale = fid_mass * livetime self.energy_hist *= scale - self.n_time_bins = len(self.energy_hist.bin_edges[0]) - 1 - e_centers = fd_nest.WIMPEnergySpectrum.bin_centers(self.energy_hist.bin_edges[1]) + self.n_time_bins = len(self.energy_hist.bin_centers()[0]) + e_centers = self.energy_hist.bin_centers()[1] self.energies = fd.np_to_tf(e_centers) self.array_columns = (("energy_spectrum", len(e_centers)),) From 6fc54877475d43229420d4ee03861e1d7e06964b Mon Sep 17 00:00:00 2001 From: Lorenzo Principe <28869147+lorenzomag@users.noreply.github.com> Date: Mon, 12 Aug 2024 14:32:43 +0200 Subject: [PATCH 04/22] Fix normalisation bug --- flamedisx/nest/lxe_blocks/energy_spectrum.py | 2 +- flamedisx/nest/lxe_sources.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/flamedisx/nest/lxe_blocks/energy_spectrum.py b/flamedisx/nest/lxe_blocks/energy_spectrum.py index e861a95f..91cb7c69 100644 --- a/flamedisx/nest/lxe_blocks/energy_spectrum.py +++ b/flamedisx/nest/lxe_blocks/energy_spectrum.py @@ -449,7 +449,7 @@ def clip_j2000_times(self, ts): return np.clip(ts, tbins[0], tbins[-1]) def mu_before_efficiencies(self, **params): - return self.energy_hist.n / self.n_time_bins + return self.energy_hist.n def random_truth(self, n_events, fix_truth=None, **params): """Draw n_events random energies and times from the energy/ diff --git a/flamedisx/nest/lxe_sources.py b/flamedisx/nest/lxe_sources.py index 6089480d..97887feb 100644 --- a/flamedisx/nest/lxe_sources.py +++ b/flamedisx/nest/lxe_sources.py @@ -565,8 +565,8 @@ def __init__(self, *args, wimp_mass=40, fid_mass=1., livetime=1., **kwargs): scale = fid_mass * livetime self.energy_hist *= scale - self.n_time_bins = len(self.energy_hist.bin_edges[0]) - 1 - e_centers = fd_nest.WIMPEnergySpectrum.bin_centers(self.energy_hist.bin_edges[1]) + self.n_time_bins = len(self.energy_hist.bin_centers()[0]) + e_centers = self.energy_hist.bin_centers()[1] self.energies = fd.np_to_tf(e_centers) self.array_columns = (('energy_spectrum', len(e_centers)),) From 31693d9de9aaf84142971ca8103a41ad9b9037e6 Mon Sep 17 00:00:00 2001 From: Lorenzo Principe <28869147+lorenzomag@users.noreply.github.com> Date: Mon, 12 Aug 2024 14:53:31 +0200 Subject: [PATCH 05/22] Revert "Revert "Merge branch 'RJ-XLZD_simple' of github.com:FlamTeam/flamedisx into RJ-XLZD_simple"" This reverts commit 65da548dd24e3396c9707322588f78ce6320eca9. --- flamedisx/nest/lxe_sources.py | 128 +++++++++++++++++++++++++++++++--- flamedisx/xlzd/xlzd.py | 40 +++++++++-- requirements_minimal.txt | 1 + 3 files changed, 153 insertions(+), 16 deletions(-) diff --git a/flamedisx/nest/lxe_sources.py b/flamedisx/nest/lxe_sources.py index 97887feb..23a59ee0 100644 --- a/flamedisx/nest/lxe_sources.py +++ b/flamedisx/nest/lxe_sources.py @@ -1,15 +1,17 @@ -import tensorflow as tf - import configparser +import math as m import os -import pickle as pkl +import numericalunits as nu +import numpy as np import pandas as pd +import tensorflow as tf import flamedisx as fd +import multihist as mh +import wimprates as wr from .. import nest as fd_nest -import math as m pi = tf.constant(m.pi) export, __all__ = fd.exporter() @@ -552,16 +554,46 @@ class nestWIMPSource(nestNRSource): """SI WIMP signal source. Reads in energy spectrum from .pkl file, generated with LZ's DMCalc. Normalised with an exposure of 1 tonne year and a cross section of 1e-45 cm^2. - Temporally varying between 2019-09-01T08:28:00 and 2020-09-01T08:28:00. + By default, temporally varying between 2019-09-01T08:28:00 and 2020-09-01T08:28:00. """ model_blocks = (fd_nest.WIMPEnergySpectrum,) + nestNRSource.model_blocks[1:] - def __init__(self, *args, wimp_mass=40, fid_mass=1., livetime=1., **kwargs): - if ('detector' not in kwargs): - kwargs['detector'] = 'default' + def __init__( + self, + *args, + wimp_mass=40, + sigma=1e-45, + fid_mass=1.0, + min_E=1e-2, + max_E=80.0, + n_energy_bins=800, + min_time="2019-09-01T08:28:00", + max_time="2020-09-01T08:28:00", + livetime=1.0, + n_time_bins=25, + modulation=True, + **kwargs + ): + + if (livetime < 0 or livetime > 1 )and not isinstance(livetime, float): + raise ValueError("Livetime should be a percentage between 0 and 1") + + if "detector" not in kwargs: + kwargs["detector"] = "default" + + self.energy_hist = self.get_energy_hist( + wimp_mass=wimp_mass, + sigma=sigma, + min_E=min_E, + max_E=max_E, + n_energy_bins=n_energy_bins, + min_time=min_time, + max_time=max_time, + n_time_bins=n_time_bins, + modulation=modulation, + ) - self.energy_hist = pkl.load(open(os.path.join(os.path.dirname(__file__), 'wimp_spectra/WIMP_spectra.pkl'), 'rb'))[wimp_mass] scale = fid_mass * livetime self.energy_hist *= scale @@ -569,6 +601,82 @@ def __init__(self, *args, wimp_mass=40, fid_mass=1., livetime=1., **kwargs): e_centers = self.energy_hist.bin_centers()[1] self.energies = fd.np_to_tf(e_centers) - self.array_columns = (('energy_spectrum', len(e_centers)),) + self.array_columns = (("energy_spectrum", len(e_centers)),) super().__init__(*args, **kwargs) + + @staticmethod + def get_energy_hist( + wimp_mass=40.0, + sigma=1e-45, + min_E=1e-2, + max_E=80.0, + n_energy_bins=800, + min_time="2019-09-01T08:28:00", + max_time="2020-09-01T08:28:00", + n_time_bins=25, + modulation=True, + ): + """ + Function to generate the energy histogram for a given WIMP mass + + Parameters + ---------- + wimp_mass : float + Mass of the WIMP in GeV/c^2 + min_E : float + Minimum energy of the histogram in keV + max_E : float + Maximum energy of the histogram in keV + sigma : float + Cross-section of the WIMP in pb + n_bins_energy : int + Number of bins for the energy histogram + n_bins_time : int + Number of bins for the time histogram. + To total timespan to bin is currently hardcoded to 1 year + modulation : bool + Flag to enable/disable modulation + """ + + energy_bin_edges = np.linspace(min_E, max_E, n_energy_bins + 1) + energy_bin_width = (energy_bin_edges[1] - energy_bin_edges[0]) * nu.keV + energy_values = (energy_bin_edges[1:] + energy_bin_edges[:-1]) / 2 + + time_bin_edges = ( + pd.date_range(min_time, max_time, periods=n_time_bins + 1).to_julian_date() + - 2451545.0 # Convert to J2000 + ) + time_bin_width = wr.j2000_to_datetime(time_bin_edges[1]) - wr.j2000_to_datetime(time_bin_edges[0]) + time_bin_width = time_bin_width.value * nu.ns + times = (time_bin_edges[:-1] + time_bin_edges[1:]) / 2 + + scale = time_bin_width / nu.year # Convert from [per year] to [per time_bin_width] + scale *= energy_bin_width / nu.keV # Convert from [per keV] to [per energy_bin_width] + + rates_list = [] + for i, time in enumerate(times): + if modulation: + rates = wr.rate_wimp_std( + energy_values, + mw=wimp_mass, + sigma_nucleon=sigma, + t=time, + ) + elif i == 0: + rates = wr.rate_wimp_std( + energy_values, + mw=wimp_mass, + sigma_nucleon=sigma, + ) + rates_list.append( + rates * scale + ) + + RATES = np.array(rates_list) + + hist = mh.Histdd.from_histogram( + RATES, [time_bin_edges, energy_bin_edges], axis_names=("time", "energy") + ) + + return hist \ No newline at end of file diff --git a/flamedisx/xlzd/xlzd.py b/flamedisx/xlzd/xlzd.py index cbef218d..b0f5e36f 100644 --- a/flamedisx/xlzd/xlzd.py +++ b/flamedisx/xlzd/xlzd.py @@ -133,12 +133,40 @@ def __init__(self, *args, **kwargs): @export class XLZDWIMPSource(XLZDSource, fd.nest.nestWIMPSource): - def __init__(self, *args, **kwargs): - if ('detector' not in kwargs): - kwargs['detector'] = 'xlzd' - if ('configuration' not in kwargs): - kwargs['configuration'] = '80t' - super().__init__(*args, **kwargs) + + def __init__( + self, + *args, + wimp_mass=40, + sigma=1e-45, + fid_mass=1.0, + min_E=1e-2, + max_E=80.0, + n_energy_bins=800, + min_time="2019-09-01T08:28:00", + max_time="2020-09-01T08:28:00", + n_time_bins=25, + modulation=True, + **kwargs + ): + if "detector" not in kwargs: + kwargs["detector"] = "xlzd" + if "configuration" not in kwargs: + kwargs["configuration"] = "80t" + super().__init__( + *args, + wimp_mass=wimp_mass, + sigma=sigma, + fid_mass=fid_mass, + min_E=min_E, + max_E=max_E, + n_energy_bins=n_energy_bins, + min_time=min_time, + max_time=max_time, + n_time_bins=n_time_bins, + modulation=modulation, + **kwargs + ) ## diff --git a/requirements_minimal.txt b/requirements_minimal.txt index 34362cc5..ee3922e4 100644 --- a/requirements_minimal.txt +++ b/requirements_minimal.txt @@ -1,3 +1,4 @@ +numericalunits numpy scipy pandas From f1395ac589c9f2d680471152dfd243e25af02b28 Mon Sep 17 00:00:00 2001 From: Lorenzo Principe <28869147+lorenzomag@users.noreply.github.com> Date: Thu, 8 Aug 2024 02:18:33 +0200 Subject: [PATCH 06/22] Initial draft of Migdal source --- flamedisx/nest/lxe_sources.py | 118 ++++++++++++++++++++++++++++++++++ 1 file changed, 118 insertions(+) diff --git a/flamedisx/nest/lxe_sources.py b/flamedisx/nest/lxe_sources.py index 3ad04587..6a075cc8 100644 --- a/flamedisx/nest/lxe_sources.py +++ b/flamedisx/nest/lxe_sources.py @@ -686,3 +686,121 @@ def get_energy_hist( ) return hist + + +@export +class nestMigdalSource(nestERSource): + + model_blocks = (fd_nest.WIMPEnergySpectrum,) + nestERSource.model_blocks[1:] + + def __init__( + self, + *args, + wimp_mass=40, + sigma=1e-45, + fid_mass=1.0, + min_E=1e-2, + max_E=80.0, + n_energy_bins=800, + min_time="2019-09-01T08:28:00", + max_time="2020-09-01T08:28:00", + n_time_bins=25, + modulation=True, + migdal_model="Cox", + **kwargs + ): + if migdal_model not in ["Ibe", "Cox", "Cox_dipole"]: + raise ValueError("Invalid Migdal model. Choose from 'Ibe', 'Cox', 'Cox_dipole'") + + if "detector" not in kwargs: + kwargs["detector"] = "default" + + self.energy_hist = self.get_energy_hist( + wimp_mass=wimp_mass, + sigma=sigma, + min_E=min_E, + max_E=max_E, + n_energy_bins=n_energy_bins, + min_time=min_time, + max_time=max_time, + n_time_bins=n_time_bins, + modulation=modulation, + migdal_model=migdal_model, + ) + + livetime = pd.Timestamp(max_time) - pd.Timestamp(min_time) + livetime = livetime.days / 365.25 + scale = fid_mass * livetime + self.energy_hist *= scale + + self.n_time_bins = len(self.energy_hist.bin_edges[0]) + e_centers = fd_nest.WIMPEnergySpectrum.bin_centers( + self.energy_hist.bin_edges[1] + ) + self.energies = fd.np_to_tf(e_centers) + + self.array_columns = (("energy_spectrum", len(e_centers)),) + + super().__init__(*args, **kwargs) + + @staticmethod + def get_energy_hist( + wimp_mass=40.0, + sigma=1e-45, + min_E=1e-2, + max_E=40.0, + n_energy_bins=800, + min_time="2019-09-01T08:28:00", + max_time="2020-09-01T08:28:00", + n_time_bins=25, + modulation=True, + migdal_model="Cox", + ): + dipole = False + if migdal_model == "Cox_dipole": + migdal_model = "Cox" + dipole = True + + + energy_bin_edges = np.linspace(min_E, max_E, n_energy_bins + 1) + energy_values = (energy_bin_edges[1:] + energy_bin_edges[:-1]) / 2 + time_bin_edges = ( + pd.date_range(min_time, max_time, periods=n_time_bins + 1).to_julian_date() + - 2451545.0 # Convert to J2000 + ) + times = (time_bin_edges[1:] + time_bin_edges[:-1]) / 2 + + rates_list = [] + if modulation: + # This should probably be done in a more efficient way + for time in times: + rates_list.append( + wr.rate_wimp_std( + energy_values, + mw=wimp_mass, + sigma_nucleon=sigma, + t=time, + detection_mechanism="migdal", + migdal_model=migdal_model, + dipole=dipole, + ) + ) + else: + rates = wr.rate_wimp_std( + energy_values, + mw=wimp_mass, + sigma_nucleon=sigma, + detection_mechanism="migdal", + migdal_model=migdal_model, + dipole=dipole, + ) + for _ in times: + rates_list.append(rates) + + RATES = np.array(rates_list) + + hist = mh.Histdd.from_histogram( + RATES, [time_bin_edges, energy_bin_edges], axis_names=("time", "energy") + ) + + return hist From 93b261bcac388f80e5c4b4d111303861f8627257 Mon Sep 17 00:00:00 2001 From: Lorenzo Principe <28869147+lorenzomag@users.noreply.github.com> Date: Tue, 13 Aug 2024 00:37:39 +0200 Subject: [PATCH 07/22] Updating Migdal source to match WIMP source --- flamedisx/nest/lxe_sources.py | 68 ++++++++++++++++++++--------------- 1 file changed, 39 insertions(+), 29 deletions(-) diff --git a/flamedisx/nest/lxe_sources.py b/flamedisx/nest/lxe_sources.py index 6a075cc8..5a838b2e 100644 --- a/flamedisx/nest/lxe_sources.py +++ b/flamedisx/nest/lxe_sources.py @@ -704,6 +704,7 @@ def __init__( n_energy_bins=800, min_time="2019-09-01T08:28:00", max_time="2020-09-01T08:28:00", + livetime=1.0, n_time_bins=25, modulation=True, migdal_model="Cox", @@ -711,7 +712,10 @@ def __init__( ): if migdal_model not in ["Ibe", "Cox", "Cox_dipole"]: raise ValueError("Invalid Migdal model. Choose from 'Ibe', 'Cox', 'Cox_dipole'") - + + if (livetime < 0 or livetime > 1) and not isinstance(livetime, float): + raise ValueError("Livetime should be a percentage between 0 and 1") + if "detector" not in kwargs: kwargs["detector"] = "default" @@ -728,15 +732,11 @@ def __init__( migdal_model=migdal_model, ) - livetime = pd.Timestamp(max_time) - pd.Timestamp(min_time) - livetime = livetime.days / 365.25 scale = fid_mass * livetime self.energy_hist *= scale - self.n_time_bins = len(self.energy_hist.bin_edges[0]) - e_centers = fd_nest.WIMPEnergySpectrum.bin_centers( - self.energy_hist.bin_edges[1] - ) + self.n_time_bins = len(self.energy_hist.bin_centers()[0]) + e_centers = self.energy_hist.bin_centers()[1] self.energies = fd.np_to_tf(e_centers) self.array_columns = (("energy_spectrum", len(e_centers)),) @@ -748,7 +748,7 @@ def get_energy_hist( wimp_mass=40.0, sigma=1e-45, min_E=1e-2, - max_E=40.0, + max_E=80.0, n_energy_bins=800, min_time="2019-09-01T08:28:00", max_time="2020-09-01T08:28:00", @@ -763,39 +763,49 @@ def get_energy_hist( energy_bin_edges = np.linspace(min_E, max_E, n_energy_bins + 1) + energy_bin_width = (energy_bin_edges[1] - energy_bin_edges[0]) * nu.keV energy_values = (energy_bin_edges[1:] + energy_bin_edges[:-1]) / 2 + time_bin_edges = ( pd.date_range(min_time, max_time, periods=n_time_bins + 1).to_julian_date() - 2451545.0 # Convert to J2000 ) - times = (time_bin_edges[1:] + time_bin_edges[:-1]) / 2 + time_bin_width = wr.j2000_to_datetime(time_bin_edges[1]) - wr.j2000_to_datetime(time_bin_edges[0]) + time_bin_width = time_bin_width.value * nu.ns + times = (time_bin_edges[:-1] + time_bin_edges[1:]) / 2 + + scale = time_bin_width / nu.year # Convert from [per year] to [per time_bin_width] + scale *= energy_bin_width / nu.keV # Convert from [per keV] to [per energy_bin_width] rates_list = [] if modulation: - # This should probably be done in a more efficient way - for time in times: - rates_list.append( - wr.rate_wimp_std( - energy_values, - mw=wimp_mass, - sigma_nucleon=sigma, - t=time, - detection_mechanism="migdal", - migdal_model=migdal_model, - dipole=dipole, + with ProcessPoolExecutor(4) as executor: + for time in times: + rates_list.append( + executor.submit( + wr.rate_wimp_std, + energy_values, + mw=wimp_mass, + sigma_nucleon=sigma, + t=time, + detection_mechanism="migdal", + migdal_model=migdal_model, + dipole=dipole, + ) ) - ) + + rates_list = [future.result() for future in rates_list] else: rates = wr.rate_wimp_std( - energy_values, - mw=wimp_mass, - sigma_nucleon=sigma, - detection_mechanism="migdal", - migdal_model=migdal_model, - dipole=dipole, - ) + energy_values, + mw=wimp_mass, + sigma_nucleon=sigma, + detection_mechanism="migdal", + migdal_model=migdal_model, + dipole=dipole, + ) for _ in times: - rates_list.append(rates) + rates_list.append(rates * scale) RATES = np.array(rates_list) From 30fb29baac0d153d8f6cf3663cd5dfaa198cc05d Mon Sep 17 00:00:00 2001 From: Lorenzo Principe <28869147+lorenzomag@users.noreply.github.com> Date: Tue, 13 Aug 2024 00:43:56 +0200 Subject: [PATCH 08/22] Remove duplicated class --- flamedisx/nest/lxe_sources.py | 120 +--------------------------------- 1 file changed, 1 insertion(+), 119 deletions(-) diff --git a/flamedisx/nest/lxe_sources.py b/flamedisx/nest/lxe_sources.py index 96f5ce36..49103d0f 100644 --- a/flamedisx/nest/lxe_sources.py +++ b/flamedisx/nest/lxe_sources.py @@ -813,122 +813,4 @@ def get_energy_hist( RATES, [time_bin_edges, energy_bin_edges], axis_names=("time", "energy") ) - return hist - - -@export -class nestMigdalSource(nestERSource): - - model_blocks = (fd_nest.WIMPEnergySpectrum,) + nestERSource.model_blocks[1:] - - def __init__( - self, - *args, - wimp_mass=40, - sigma=1e-45, - fid_mass=1.0, - min_E=1e-2, - max_E=80.0, - n_energy_bins=800, - min_time="2019-09-01T08:28:00", - max_time="2020-09-01T08:28:00", - n_time_bins=25, - modulation=True, - migdal_model="Cox", - **kwargs - ): - if migdal_model not in ["Ibe", "Cox", "Cox_dipole"]: - raise ValueError("Invalid Migdal model. Choose from 'Ibe', 'Cox', 'Cox_dipole'") - - if "detector" not in kwargs: - kwargs["detector"] = "default" - - self.energy_hist = self.get_energy_hist( - wimp_mass=wimp_mass, - sigma=sigma, - min_E=min_E, - max_E=max_E, - n_energy_bins=n_energy_bins, - min_time=min_time, - max_time=max_time, - n_time_bins=n_time_bins, - modulation=modulation, - migdal_model=migdal_model, - ) - - livetime = pd.Timestamp(max_time) - pd.Timestamp(min_time) - livetime = livetime.days / 365.25 - scale = fid_mass * livetime - self.energy_hist *= scale - - self.n_time_bins = len(self.energy_hist.bin_edges[0]) - e_centers = fd_nest.WIMPEnergySpectrum.bin_centers( - self.energy_hist.bin_edges[1] - ) - self.energies = fd.np_to_tf(e_centers) - - self.array_columns = (("energy_spectrum", len(e_centers)),) - - super().__init__(*args, **kwargs) - - @staticmethod - def get_energy_hist( - wimp_mass=40.0, - sigma=1e-45, - min_E=1e-2, - max_E=40.0, - n_energy_bins=800, - min_time="2019-09-01T08:28:00", - max_time="2020-09-01T08:28:00", - n_time_bins=25, - modulation=True, - migdal_model="Cox", - ): - dipole = False - if migdal_model == "Cox_dipole": - migdal_model = "Cox" - dipole = True - - - energy_bin_edges = np.linspace(min_E, max_E, n_energy_bins + 1) - energy_values = (energy_bin_edges[1:] + energy_bin_edges[:-1]) / 2 - time_bin_edges = ( - pd.date_range(min_time, max_time, periods=n_time_bins + 1).to_julian_date() - - 2451545.0 # Convert to J2000 - ) - times = (time_bin_edges[1:] + time_bin_edges[:-1]) / 2 - - rates_list = [] - if modulation: - # This should probably be done in a more efficient way - for time in times: - rates_list.append( - wr.rate_wimp_std( - energy_values, - mw=wimp_mass, - sigma_nucleon=sigma, - t=time, - detection_mechanism="migdal", - migdal_model=migdal_model, - dipole=dipole, - ) - ) - else: - rates = wr.rate_wimp_std( - energy_values, - mw=wimp_mass, - sigma_nucleon=sigma, - detection_mechanism="migdal", - migdal_model=migdal_model, - dipole=dipole, - ) - for _ in times: - rates_list.append(rates) - - RATES = np.array(rates_list) - - hist = mh.Histdd.from_histogram( - RATES, [time_bin_edges, energy_bin_edges], axis_names=("time", "energy") - ) - - return hist + return hist \ No newline at end of file From 256a99bf8cbc1a28604934f369247e3596f43223 Mon Sep 17 00:00:00 2001 From: Lorenzo Principe <28869147+lorenzomag@users.noreply.github.com> Date: Tue, 13 Aug 2024 14:44:02 +0200 Subject: [PATCH 09/22] Update XLZD sources --- flamedisx/xlzd/xlzd.py | 44 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 44 insertions(+) diff --git a/flamedisx/xlzd/xlzd.py b/flamedisx/xlzd/xlzd.py index b0f5e36f..334bab08 100644 --- a/flamedisx/xlzd/xlzd.py +++ b/flamedisx/xlzd/xlzd.py @@ -145,6 +145,7 @@ def __init__( n_energy_bins=800, min_time="2019-09-01T08:28:00", max_time="2020-09-01T08:28:00", + livetime=1.0, n_time_bins=25, modulation=True, **kwargs @@ -163,12 +164,55 @@ def __init__( n_energy_bins=n_energy_bins, min_time=min_time, max_time=max_time, + livetime=livetime, n_time_bins=n_time_bins, modulation=modulation, **kwargs ) +@export +class XLZDMigdalSource(XLZDSource, fd.nest.nestMigdalSource): + + def __init__( + self, + *args, + wimp_mass=40, + sigma=1e-45, + fid_mass=1.0, + min_E=1e-2, + max_E=80.0, + n_energy_bins=800, + min_time="2019-09-01T08:28:00", + max_time="2020-09-01T08:28:00", + livetime=1.0, + n_time_bins=25, + modulation=True, + migdal_model="Cox", + **kwargs + ): + if "detector" not in kwargs: + kwargs["detector"] = "xlzd" + if "configuration" not in kwargs: + kwargs["configuration"] = "80t" + super().__init__( + *args, + wimp_mass=wimp_mass, + sigma=sigma, + fid_mass=fid_mass, + min_E=min_E, + max_E=max_E, + n_energy_bins=n_energy_bins, + min_time=min_time, + max_time=max_time, + livetime=livetime, + n_time_bins=n_time_bins, + modulation=modulation, + migdal_model=migdal_model, + **kwargs + ) + + ## # Background sources ## From d7c29e4cb953df067fb66ec466d1677eb2eeb868 Mon Sep 17 00:00:00 2001 From: Lorenzo Principe <28869147+lorenzomag@users.noreply.github.com> Date: Tue, 13 Aug 2024 16:40:34 +0200 Subject: [PATCH 10/22] Fix scaling bug Scaling was only applied when modulation off --- flamedisx/nest/lxe_sources.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/flamedisx/nest/lxe_sources.py b/flamedisx/nest/lxe_sources.py index 49103d0f..dd320abc 100644 --- a/flamedisx/nest/lxe_sources.py +++ b/flamedisx/nest/lxe_sources.py @@ -657,7 +657,7 @@ def get_energy_hist( rates_list = [] if modulation: - with ProcessPoolExecutor(4) as executor: + with ProcessPoolExecutor() as executor: for time in times: rates_list.append( executor.submit( @@ -677,9 +677,9 @@ def get_energy_hist( sigma_nucleon=sigma, ) for _ in times: - rates_list.append(rates * scale) + rates_list.append(rates) - RATES = np.array(rates_list) + RATES = np.array(rates_list) * scale hist = mh.Histdd.from_histogram( RATES, [time_bin_edges, energy_bin_edges], axis_names=("time", "energy") @@ -779,7 +779,7 @@ def get_energy_hist( rates_list = [] if modulation: - with ProcessPoolExecutor(4) as executor: + with ProcessPoolExecutor() as executor: for time in times: rates_list.append( executor.submit( @@ -805,9 +805,9 @@ def get_energy_hist( dipole=dipole, ) for _ in times: - rates_list.append(rates * scale) + rates_list.append(rates) - RATES = np.array(rates_list) + RATES = np.array(rates_list) * scale hist = mh.Histdd.from_histogram( RATES, [time_bin_edges, energy_bin_edges], axis_names=("time", "energy") From b7771e65b06f8983b165e2cf81a7c2c6b4f4ab83 Mon Sep 17 00:00:00 2001 From: Lorenzo Principe <28869147+lorenzomag@users.noreply.github.com> Date: Wed, 14 Aug 2024 12:15:29 +0200 Subject: [PATCH 11/22] Code refactoring and normalisation fix The energy histogram is not normalised in the time dimension. Robert's expalation: I think all you need to do is remove scale = time_bin_width / nu.year # Convert from [per year] to [per time_bin_width] Whilst in principle keeping this and removing / self.n_time_bins would solve things at the level of source.mu_before_efficiencies(), it will mess up the tensor-based calculation of the differential rate, as the thing that is queried should be the spectrum corresponding to a given event time (such that the sum of the spectrum gives the expected counts within that energy range, after scaling by the exposure) --- flamedisx/nest/__init__.py | 1 + flamedisx/nest/get_energy_hist.py | 88 ++++++++++ flamedisx/nest/lxe_blocks/energy_spectrum.py | 2 +- flamedisx/nest/lxe_sources.py | 161 +------------------ 4 files changed, 94 insertions(+), 158 deletions(-) create mode 100644 flamedisx/nest/get_energy_hist.py diff --git a/flamedisx/nest/__init__.py b/flamedisx/nest/__init__.py index 15c530fc..b45a981a 100644 --- a/flamedisx/nest/__init__.py +++ b/flamedisx/nest/__init__.py @@ -1,4 +1,5 @@ from .parameter_calc import * +from .get_energy_hist import * from .lxe_blocks.energy_spectrum import * from .lxe_blocks.secondary_quanta_generation import * diff --git a/flamedisx/nest/get_energy_hist.py b/flamedisx/nest/get_energy_hist.py new file mode 100644 index 00000000..223f77e7 --- /dev/null +++ b/flamedisx/nest/get_energy_hist.py @@ -0,0 +1,88 @@ +from concurrent.futures import ProcessPoolExecutor + +import numericalunits as nu +import numpy as np +import pandas as pd + +import flamedisx as fd +import multihist as mh +import wimprates as wr + + +export, __all__ = fd.exporter() + + +@export +def get_energy_hist( + wimp_mass=40.0, + sigma=1e-45, + min_E=1e-2, + max_E=80.0, + n_energy_bins=800, + min_time="2019-09-01T08:28:00", + max_time="2020-09-01T08:28:00", + n_time_bins=25, + modulation=True, + migdal_model=None, +): + + wimprates_kwargs = dict( + mw=wimp_mass, + sigma_nucleon=sigma, + ) + + if migdal_model is not None: + if migdal_model not in ["Ibe", "Cox", "Cox_dipole"]: + raise ValueError( + "Invalid Migdal model. Choose from 'Ibe', 'Cox', 'Cox_dipole'" + ) + dipole = False + if migdal_model == "Cox_dipole": + migdal_model = "Cox" + dipole = True + + wimprates_kwargs.update( + dict( + detection_mechanism="migdal", + migdal_model=migdal_model, + dipole=dipole, + ) + ) + + energy_bin_edges = np.linspace(min_E, max_E, n_energy_bins + 1) + energy_bin_width = (energy_bin_edges[1] - energy_bin_edges[0]) * nu.keV + energy_values = (energy_bin_edges[1:] + energy_bin_edges[:-1]) / 2 + + time_bin_edges = ( + pd.date_range(min_time, max_time, periods=n_time_bins + 1).to_julian_date() + - 2451545.0 # Convert to J2000 + ) + times = (time_bin_edges[:-1] + time_bin_edges[1:]) / 2 + + scale = ( + energy_bin_width / nu.keV + ) # Convert from [per keV] to [per energy_bin_width] + + rates_list = [] + if modulation: + with ProcessPoolExecutor() as executor: + for time in times: + rates_list.append( + executor.submit( + wr.rate_wimp_std, energy_values, t=time, **wimprates_kwargs + ) + ) + + rates_list = [future.result() for future in rates_list] + else: + rates = wr.rate_wimp_std(energy_values, **wimprates_kwargs) + for _ in times: + rates_list.append(rates) + + RATES = np.array(rates_list) * scale + + hist = mh.Histdd.from_histogram( + RATES, [time_bin_edges, energy_bin_edges], axis_names=("time", "energy") + ) + + return hist diff --git a/flamedisx/nest/lxe_blocks/energy_spectrum.py b/flamedisx/nest/lxe_blocks/energy_spectrum.py index 91cb7c69..e861a95f 100644 --- a/flamedisx/nest/lxe_blocks/energy_spectrum.py +++ b/flamedisx/nest/lxe_blocks/energy_spectrum.py @@ -449,7 +449,7 @@ def clip_j2000_times(self, ts): return np.clip(ts, tbins[0], tbins[-1]) def mu_before_efficiencies(self, **params): - return self.energy_hist.n + return self.energy_hist.n / self.n_time_bins def random_truth(self, n_events, fix_truth=None, **params): """Draw n_events random energies and times from the energy/ diff --git a/flamedisx/nest/lxe_sources.py b/flamedisx/nest/lxe_sources.py index dd320abc..2020d986 100644 --- a/flamedisx/nest/lxe_sources.py +++ b/flamedisx/nest/lxe_sources.py @@ -573,7 +573,7 @@ def __init__( max_time="2020-09-01T08:28:00", livetime=1.0, n_time_bins=25, - modulation=True, + modulation=False, **kwargs ): @@ -583,7 +583,7 @@ def __init__( if "detector" not in kwargs: kwargs["detector"] = "default" - self.energy_hist = self.get_energy_hist( + self.energy_hist = fd_nest.get_energy_hist( wimp_mass=wimp_mass, sigma=sigma, min_E=min_E, @@ -606,87 +606,6 @@ def __init__( super().__init__(*args, **kwargs) - @staticmethod - def get_energy_hist( - wimp_mass=40.0, - sigma=1e-45, - min_E=1e-2, - max_E=80.0, - n_energy_bins=800, - min_time="2019-09-01T08:28:00", - max_time="2020-09-01T08:28:00", - n_time_bins=25, - modulation=True, - ): - """ - Function to generate the energy histogram for a given WIMP mass - - Parameters - ---------- - wimp_mass : float - Mass of the WIMP in GeV/c^2 - min_E : float - Minimum energy of the histogram in keV - max_E : float - Maximum energy of the histogram in keV - sigma : float - Cross-section of the WIMP in pb - n_bins_energy : int - Number of bins for the energy histogram - n_bins_time : int - Number of bins for the time histogram. - To total timespan to bin is currently hardcoded to 1 year - modulation : bool - Flag to enable/disable modulation - """ - - energy_bin_edges = np.linspace(min_E, max_E, n_energy_bins + 1) - energy_bin_width = (energy_bin_edges[1] - energy_bin_edges[0]) * nu.keV - energy_values = (energy_bin_edges[1:] + energy_bin_edges[:-1]) / 2 - - time_bin_edges = ( - pd.date_range(min_time, max_time, periods=n_time_bins + 1).to_julian_date() - - 2451545.0 # Convert to J2000 - ) - time_bin_width = wr.j2000_to_datetime(time_bin_edges[1]) - wr.j2000_to_datetime(time_bin_edges[0]) - time_bin_width = time_bin_width.value * nu.ns - times = (time_bin_edges[:-1] + time_bin_edges[1:]) / 2 - - scale = time_bin_width / nu.year # Convert from [per year] to [per time_bin_width] - scale *= energy_bin_width / nu.keV # Convert from [per keV] to [per energy_bin_width] - - rates_list = [] - if modulation: - with ProcessPoolExecutor() as executor: - for time in times: - rates_list.append( - executor.submit( - wr.rate_wimp_std, - energy_values, - mw=wimp_mass, - sigma_nucleon=sigma, - t=time, - ) - ) - - rates_list = [future.result() for future in rates_list] - else: - rates = wr.rate_wimp_std( - energy_values, - mw=wimp_mass, - sigma_nucleon=sigma, - ) - for _ in times: - rates_list.append(rates) - - RATES = np.array(rates_list) * scale - - hist = mh.Histdd.from_histogram( - RATES, [time_bin_edges, energy_bin_edges], axis_names=("time", "energy") - ) - - return hist - @export class nestMigdalSource(nestERSource): @@ -719,7 +638,7 @@ def __init__( if "detector" not in kwargs: kwargs["detector"] = "default" - self.energy_hist = self.get_energy_hist( + self.energy_hist = fd_nest.get_energy_hist( wimp_mass=wimp_mass, sigma=sigma, min_E=min_E, @@ -741,76 +660,4 @@ def __init__( self.array_columns = (("energy_spectrum", len(e_centers)),) - super().__init__(*args, **kwargs) - - @staticmethod - def get_energy_hist( - wimp_mass=40.0, - sigma=1e-45, - min_E=1e-2, - max_E=80.0, - n_energy_bins=800, - min_time="2019-09-01T08:28:00", - max_time="2020-09-01T08:28:00", - n_time_bins=25, - modulation=True, - migdal_model="Cox", - ): - dipole = False - if migdal_model == "Cox_dipole": - migdal_model = "Cox" - dipole = True - - - energy_bin_edges = np.linspace(min_E, max_E, n_energy_bins + 1) - energy_bin_width = (energy_bin_edges[1] - energy_bin_edges[0]) * nu.keV - energy_values = (energy_bin_edges[1:] + energy_bin_edges[:-1]) / 2 - - time_bin_edges = ( - pd.date_range(min_time, max_time, periods=n_time_bins + 1).to_julian_date() - - 2451545.0 # Convert to J2000 - ) - time_bin_width = wr.j2000_to_datetime(time_bin_edges[1]) - wr.j2000_to_datetime(time_bin_edges[0]) - time_bin_width = time_bin_width.value * nu.ns - times = (time_bin_edges[:-1] + time_bin_edges[1:]) / 2 - - scale = time_bin_width / nu.year # Convert from [per year] to [per time_bin_width] - scale *= energy_bin_width / nu.keV # Convert from [per keV] to [per energy_bin_width] - - rates_list = [] - if modulation: - with ProcessPoolExecutor() as executor: - for time in times: - rates_list.append( - executor.submit( - wr.rate_wimp_std, - energy_values, - mw=wimp_mass, - sigma_nucleon=sigma, - t=time, - detection_mechanism="migdal", - migdal_model=migdal_model, - dipole=dipole, - ) - ) - - rates_list = [future.result() for future in rates_list] - else: - rates = wr.rate_wimp_std( - energy_values, - mw=wimp_mass, - sigma_nucleon=sigma, - detection_mechanism="migdal", - migdal_model=migdal_model, - dipole=dipole, - ) - for _ in times: - rates_list.append(rates) - - RATES = np.array(rates_list) * scale - - hist = mh.Histdd.from_histogram( - RATES, [time_bin_edges, energy_bin_edges], axis_names=("time", "energy") - ) - - return hist \ No newline at end of file + super().__init__(*args, **kwargs) \ No newline at end of file From d4b2e3175a607f9ef300a5e157eafb4da70c8f3c Mon Sep 17 00:00:00 2001 From: Lorenzo Principe <28869147+lorenzomag@users.noreply.github.com> Date: Mon, 12 Aug 2024 14:32:43 +0200 Subject: [PATCH 12/22] Fix normalisation bug --- flamedisx/nest/lxe_blocks/energy_spectrum.py | 2 +- flamedisx/nest/lxe_sources.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/flamedisx/nest/lxe_blocks/energy_spectrum.py b/flamedisx/nest/lxe_blocks/energy_spectrum.py index e861a95f..91cb7c69 100644 --- a/flamedisx/nest/lxe_blocks/energy_spectrum.py +++ b/flamedisx/nest/lxe_blocks/energy_spectrum.py @@ -449,7 +449,7 @@ def clip_j2000_times(self, ts): return np.clip(ts, tbins[0], tbins[-1]) def mu_before_efficiencies(self, **params): - return self.energy_hist.n / self.n_time_bins + return self.energy_hist.n def random_truth(self, n_events, fix_truth=None, **params): """Draw n_events random energies and times from the energy/ diff --git a/flamedisx/nest/lxe_sources.py b/flamedisx/nest/lxe_sources.py index 6089480d..97887feb 100644 --- a/flamedisx/nest/lxe_sources.py +++ b/flamedisx/nest/lxe_sources.py @@ -565,8 +565,8 @@ def __init__(self, *args, wimp_mass=40, fid_mass=1., livetime=1., **kwargs): scale = fid_mass * livetime self.energy_hist *= scale - self.n_time_bins = len(self.energy_hist.bin_edges[0]) - 1 - e_centers = fd_nest.WIMPEnergySpectrum.bin_centers(self.energy_hist.bin_edges[1]) + self.n_time_bins = len(self.energy_hist.bin_centers()[0]) + e_centers = self.energy_hist.bin_centers()[1] self.energies = fd.np_to_tf(e_centers) self.array_columns = (('energy_spectrum', len(e_centers)),) From 783d393c23f35dc4a5de6e81c3da3811652d48f4 Mon Sep 17 00:00:00 2001 From: Lorenzo Principe <28869147+lorenzomag@users.noreply.github.com> Date: Mon, 12 Aug 2024 14:53:31 +0200 Subject: [PATCH 13/22] Revert "Revert "Merge branch 'RJ-XLZD_simple' of github.com:FlamTeam/flamedisx into RJ-XLZD_simple"" This reverts commit 65da548dd24e3396c9707322588f78ce6320eca9. --- flamedisx/nest/lxe_sources.py | 128 +++++++++++++++++++++++++++++++--- flamedisx/xlzd/xlzd.py | 40 +++++++++-- requirements_minimal.txt | 1 + 3 files changed, 153 insertions(+), 16 deletions(-) diff --git a/flamedisx/nest/lxe_sources.py b/flamedisx/nest/lxe_sources.py index 97887feb..23a59ee0 100644 --- a/flamedisx/nest/lxe_sources.py +++ b/flamedisx/nest/lxe_sources.py @@ -1,15 +1,17 @@ -import tensorflow as tf - import configparser +import math as m import os -import pickle as pkl +import numericalunits as nu +import numpy as np import pandas as pd +import tensorflow as tf import flamedisx as fd +import multihist as mh +import wimprates as wr from .. import nest as fd_nest -import math as m pi = tf.constant(m.pi) export, __all__ = fd.exporter() @@ -552,16 +554,46 @@ class nestWIMPSource(nestNRSource): """SI WIMP signal source. Reads in energy spectrum from .pkl file, generated with LZ's DMCalc. Normalised with an exposure of 1 tonne year and a cross section of 1e-45 cm^2. - Temporally varying between 2019-09-01T08:28:00 and 2020-09-01T08:28:00. + By default, temporally varying between 2019-09-01T08:28:00 and 2020-09-01T08:28:00. """ model_blocks = (fd_nest.WIMPEnergySpectrum,) + nestNRSource.model_blocks[1:] - def __init__(self, *args, wimp_mass=40, fid_mass=1., livetime=1., **kwargs): - if ('detector' not in kwargs): - kwargs['detector'] = 'default' + def __init__( + self, + *args, + wimp_mass=40, + sigma=1e-45, + fid_mass=1.0, + min_E=1e-2, + max_E=80.0, + n_energy_bins=800, + min_time="2019-09-01T08:28:00", + max_time="2020-09-01T08:28:00", + livetime=1.0, + n_time_bins=25, + modulation=True, + **kwargs + ): + + if (livetime < 0 or livetime > 1 )and not isinstance(livetime, float): + raise ValueError("Livetime should be a percentage between 0 and 1") + + if "detector" not in kwargs: + kwargs["detector"] = "default" + + self.energy_hist = self.get_energy_hist( + wimp_mass=wimp_mass, + sigma=sigma, + min_E=min_E, + max_E=max_E, + n_energy_bins=n_energy_bins, + min_time=min_time, + max_time=max_time, + n_time_bins=n_time_bins, + modulation=modulation, + ) - self.energy_hist = pkl.load(open(os.path.join(os.path.dirname(__file__), 'wimp_spectra/WIMP_spectra.pkl'), 'rb'))[wimp_mass] scale = fid_mass * livetime self.energy_hist *= scale @@ -569,6 +601,82 @@ def __init__(self, *args, wimp_mass=40, fid_mass=1., livetime=1., **kwargs): e_centers = self.energy_hist.bin_centers()[1] self.energies = fd.np_to_tf(e_centers) - self.array_columns = (('energy_spectrum', len(e_centers)),) + self.array_columns = (("energy_spectrum", len(e_centers)),) super().__init__(*args, **kwargs) + + @staticmethod + def get_energy_hist( + wimp_mass=40.0, + sigma=1e-45, + min_E=1e-2, + max_E=80.0, + n_energy_bins=800, + min_time="2019-09-01T08:28:00", + max_time="2020-09-01T08:28:00", + n_time_bins=25, + modulation=True, + ): + """ + Function to generate the energy histogram for a given WIMP mass + + Parameters + ---------- + wimp_mass : float + Mass of the WIMP in GeV/c^2 + min_E : float + Minimum energy of the histogram in keV + max_E : float + Maximum energy of the histogram in keV + sigma : float + Cross-section of the WIMP in pb + n_bins_energy : int + Number of bins for the energy histogram + n_bins_time : int + Number of bins for the time histogram. + To total timespan to bin is currently hardcoded to 1 year + modulation : bool + Flag to enable/disable modulation + """ + + energy_bin_edges = np.linspace(min_E, max_E, n_energy_bins + 1) + energy_bin_width = (energy_bin_edges[1] - energy_bin_edges[0]) * nu.keV + energy_values = (energy_bin_edges[1:] + energy_bin_edges[:-1]) / 2 + + time_bin_edges = ( + pd.date_range(min_time, max_time, periods=n_time_bins + 1).to_julian_date() + - 2451545.0 # Convert to J2000 + ) + time_bin_width = wr.j2000_to_datetime(time_bin_edges[1]) - wr.j2000_to_datetime(time_bin_edges[0]) + time_bin_width = time_bin_width.value * nu.ns + times = (time_bin_edges[:-1] + time_bin_edges[1:]) / 2 + + scale = time_bin_width / nu.year # Convert from [per year] to [per time_bin_width] + scale *= energy_bin_width / nu.keV # Convert from [per keV] to [per energy_bin_width] + + rates_list = [] + for i, time in enumerate(times): + if modulation: + rates = wr.rate_wimp_std( + energy_values, + mw=wimp_mass, + sigma_nucleon=sigma, + t=time, + ) + elif i == 0: + rates = wr.rate_wimp_std( + energy_values, + mw=wimp_mass, + sigma_nucleon=sigma, + ) + rates_list.append( + rates * scale + ) + + RATES = np.array(rates_list) + + hist = mh.Histdd.from_histogram( + RATES, [time_bin_edges, energy_bin_edges], axis_names=("time", "energy") + ) + + return hist \ No newline at end of file diff --git a/flamedisx/xlzd/xlzd.py b/flamedisx/xlzd/xlzd.py index d3a846d7..e8754b13 100644 --- a/flamedisx/xlzd/xlzd.py +++ b/flamedisx/xlzd/xlzd.py @@ -157,12 +157,40 @@ def __init__(self, *args, **kwargs): @export class XLZDWIMPSource(XLZDSource, fd.nest.nestWIMPSource): - def __init__(self, *args, **kwargs): - if ('detector' not in kwargs): - kwargs['detector'] = 'xlzd' - if ('configuration' not in kwargs): - kwargs['configuration'] = '80t' - super().__init__(*args, **kwargs) + + def __init__( + self, + *args, + wimp_mass=40, + sigma=1e-45, + fid_mass=1.0, + min_E=1e-2, + max_E=80.0, + n_energy_bins=800, + min_time="2019-09-01T08:28:00", + max_time="2020-09-01T08:28:00", + n_time_bins=25, + modulation=True, + **kwargs + ): + if "detector" not in kwargs: + kwargs["detector"] = "xlzd" + if "configuration" not in kwargs: + kwargs["configuration"] = "80t" + super().__init__( + *args, + wimp_mass=wimp_mass, + sigma=sigma, + fid_mass=fid_mass, + min_E=min_E, + max_E=max_E, + n_energy_bins=n_energy_bins, + min_time=min_time, + max_time=max_time, + n_time_bins=n_time_bins, + modulation=modulation, + **kwargs + ) ## diff --git a/requirements_minimal.txt b/requirements_minimal.txt index 34362cc5..ee3922e4 100644 --- a/requirements_minimal.txt +++ b/requirements_minimal.txt @@ -1,3 +1,4 @@ +numericalunits numpy scipy pandas From 96b76ee7958737a9aae5aa5a84bb8207addeab89 Mon Sep 17 00:00:00 2001 From: Lorenzo Principe <28869147+lorenzomag@users.noreply.github.com> Date: Sun, 11 Aug 2024 22:56:52 +0200 Subject: [PATCH 14/22] Use multiprocessing when modulation on --- flamedisx/nest/lxe_sources.py | 40 ++++++++++++++++++++--------------- 1 file changed, 23 insertions(+), 17 deletions(-) diff --git a/flamedisx/nest/lxe_sources.py b/flamedisx/nest/lxe_sources.py index 23a59ee0..87645ad0 100644 --- a/flamedisx/nest/lxe_sources.py +++ b/flamedisx/nest/lxe_sources.py @@ -1,3 +1,4 @@ +from concurrent.futures import ProcessPoolExecutor import configparser import math as m import os @@ -645,7 +646,7 @@ def get_energy_hist( time_bin_edges = ( pd.date_range(min_time, max_time, periods=n_time_bins + 1).to_julian_date() - - 2451545.0 # Convert to J2000 + - 2451545.0 # Convert to J2000 ) time_bin_width = wr.j2000_to_datetime(time_bin_edges[1]) - wr.j2000_to_datetime(time_bin_edges[0]) time_bin_width = time_bin_width.value * nu.ns @@ -655,23 +656,28 @@ def get_energy_hist( scale *= energy_bin_width / nu.keV # Convert from [per keV] to [per energy_bin_width] rates_list = [] - for i, time in enumerate(times): - if modulation: - rates = wr.rate_wimp_std( - energy_values, - mw=wimp_mass, - sigma_nucleon=sigma, - t=time, - ) - elif i == 0: - rates = wr.rate_wimp_std( - energy_values, - mw=wimp_mass, - sigma_nucleon=sigma, - ) - rates_list.append( - rates * scale + if modulation: + with ProcessPoolExecutor(4) as executor: + for time in times: + rates_list.append( + executor.submit( + wr.rate_wimp_std, + energy_values, + mw=wimp_mass, + sigma_nucleon=sigma, + t=time, + ) + ) + + rates_list = [future.result() for future in rates_list] + else: + rates = wr.rate_wimp_std( + energy_values, + mw=wimp_mass, + sigma_nucleon=sigma, ) + for _ in times: + rates_list.append(rates) RATES = np.array(rates_list) From 7eb577b6cdcb71f0862eb7e7668135f0fc9a2747 Mon Sep 17 00:00:00 2001 From: Lorenzo Principe <28869147+lorenzomag@users.noreply.github.com> Date: Thu, 8 Aug 2024 02:18:33 +0200 Subject: [PATCH 15/22] Initial draft of Migdal source --- flamedisx/nest/lxe_sources.py | 120 +++++++++++++++++++++++++++++++++- 1 file changed, 119 insertions(+), 1 deletion(-) diff --git a/flamedisx/nest/lxe_sources.py b/flamedisx/nest/lxe_sources.py index 87645ad0..37055e6b 100644 --- a/flamedisx/nest/lxe_sources.py +++ b/flamedisx/nest/lxe_sources.py @@ -685,4 +685,122 @@ def get_energy_hist( RATES, [time_bin_edges, energy_bin_edges], axis_names=("time", "energy") ) - return hist \ No newline at end of file + return hist + + +@export +class nestMigdalSource(nestERSource): + + model_blocks = (fd_nest.WIMPEnergySpectrum,) + nestERSource.model_blocks[1:] + + def __init__( + self, + *args, + wimp_mass=40, + sigma=1e-45, + fid_mass=1.0, + min_E=1e-2, + max_E=80.0, + n_energy_bins=800, + min_time="2019-09-01T08:28:00", + max_time="2020-09-01T08:28:00", + n_time_bins=25, + modulation=True, + migdal_model="Cox", + **kwargs + ): + if migdal_model not in ["Ibe", "Cox", "Cox_dipole"]: + raise ValueError("Invalid Migdal model. Choose from 'Ibe', 'Cox', 'Cox_dipole'") + + if "detector" not in kwargs: + kwargs["detector"] = "default" + + self.energy_hist = self.get_energy_hist( + wimp_mass=wimp_mass, + sigma=sigma, + min_E=min_E, + max_E=max_E, + n_energy_bins=n_energy_bins, + min_time=min_time, + max_time=max_time, + n_time_bins=n_time_bins, + modulation=modulation, + migdal_model=migdal_model, + ) + + livetime = pd.Timestamp(max_time) - pd.Timestamp(min_time) + livetime = livetime.days / 365.25 + scale = fid_mass * livetime + self.energy_hist *= scale + + self.n_time_bins = len(self.energy_hist.bin_edges[0]) + e_centers = fd_nest.WIMPEnergySpectrum.bin_centers( + self.energy_hist.bin_edges[1] + ) + self.energies = fd.np_to_tf(e_centers) + + self.array_columns = (("energy_spectrum", len(e_centers)),) + + super().__init__(*args, **kwargs) + + @staticmethod + def get_energy_hist( + wimp_mass=40.0, + sigma=1e-45, + min_E=1e-2, + max_E=40.0, + n_energy_bins=800, + min_time="2019-09-01T08:28:00", + max_time="2020-09-01T08:28:00", + n_time_bins=25, + modulation=True, + migdal_model="Cox", + ): + dipole = False + if migdal_model == "Cox_dipole": + migdal_model = "Cox" + dipole = True + + + energy_bin_edges = np.linspace(min_E, max_E, n_energy_bins + 1) + energy_values = (energy_bin_edges[1:] + energy_bin_edges[:-1]) / 2 + time_bin_edges = ( + pd.date_range(min_time, max_time, periods=n_time_bins + 1).to_julian_date() + - 2451545.0 # Convert to J2000 + ) + times = (time_bin_edges[1:] + time_bin_edges[:-1]) / 2 + + rates_list = [] + if modulation: + # This should probably be done in a more efficient way + for time in times: + rates_list.append( + wr.rate_wimp_std( + energy_values, + mw=wimp_mass, + sigma_nucleon=sigma, + t=time, + detection_mechanism="migdal", + migdal_model=migdal_model, + dipole=dipole, + ) + ) + else: + rates = wr.rate_wimp_std( + energy_values, + mw=wimp_mass, + sigma_nucleon=sigma, + detection_mechanism="migdal", + migdal_model=migdal_model, + dipole=dipole, + ) + for _ in times: + rates_list.append(rates) + + RATES = np.array(rates_list) + + hist = mh.Histdd.from_histogram( + RATES, [time_bin_edges, energy_bin_edges], axis_names=("time", "energy") + ) + + return hist From 26cd38ffb16964f4bc2162d4d25e2285e425881f Mon Sep 17 00:00:00 2001 From: Lorenzo Principe <28869147+lorenzomag@users.noreply.github.com> Date: Tue, 13 Aug 2024 00:37:39 +0200 Subject: [PATCH 16/22] Updating Migdal source to match WIMP source --- flamedisx/nest/lxe_sources.py | 68 ++++++++++++++++++++--------------- 1 file changed, 39 insertions(+), 29 deletions(-) diff --git a/flamedisx/nest/lxe_sources.py b/flamedisx/nest/lxe_sources.py index 37055e6b..8f55a5b6 100644 --- a/flamedisx/nest/lxe_sources.py +++ b/flamedisx/nest/lxe_sources.py @@ -704,6 +704,7 @@ def __init__( n_energy_bins=800, min_time="2019-09-01T08:28:00", max_time="2020-09-01T08:28:00", + livetime=1.0, n_time_bins=25, modulation=True, migdal_model="Cox", @@ -711,7 +712,10 @@ def __init__( ): if migdal_model not in ["Ibe", "Cox", "Cox_dipole"]: raise ValueError("Invalid Migdal model. Choose from 'Ibe', 'Cox', 'Cox_dipole'") - + + if (livetime < 0 or livetime > 1) and not isinstance(livetime, float): + raise ValueError("Livetime should be a percentage between 0 and 1") + if "detector" not in kwargs: kwargs["detector"] = "default" @@ -728,15 +732,11 @@ def __init__( migdal_model=migdal_model, ) - livetime = pd.Timestamp(max_time) - pd.Timestamp(min_time) - livetime = livetime.days / 365.25 scale = fid_mass * livetime self.energy_hist *= scale - self.n_time_bins = len(self.energy_hist.bin_edges[0]) - e_centers = fd_nest.WIMPEnergySpectrum.bin_centers( - self.energy_hist.bin_edges[1] - ) + self.n_time_bins = len(self.energy_hist.bin_centers()[0]) + e_centers = self.energy_hist.bin_centers()[1] self.energies = fd.np_to_tf(e_centers) self.array_columns = (("energy_spectrum", len(e_centers)),) @@ -748,7 +748,7 @@ def get_energy_hist( wimp_mass=40.0, sigma=1e-45, min_E=1e-2, - max_E=40.0, + max_E=80.0, n_energy_bins=800, min_time="2019-09-01T08:28:00", max_time="2020-09-01T08:28:00", @@ -763,39 +763,49 @@ def get_energy_hist( energy_bin_edges = np.linspace(min_E, max_E, n_energy_bins + 1) + energy_bin_width = (energy_bin_edges[1] - energy_bin_edges[0]) * nu.keV energy_values = (energy_bin_edges[1:] + energy_bin_edges[:-1]) / 2 + time_bin_edges = ( pd.date_range(min_time, max_time, periods=n_time_bins + 1).to_julian_date() - 2451545.0 # Convert to J2000 ) - times = (time_bin_edges[1:] + time_bin_edges[:-1]) / 2 + time_bin_width = wr.j2000_to_datetime(time_bin_edges[1]) - wr.j2000_to_datetime(time_bin_edges[0]) + time_bin_width = time_bin_width.value * nu.ns + times = (time_bin_edges[:-1] + time_bin_edges[1:]) / 2 + + scale = time_bin_width / nu.year # Convert from [per year] to [per time_bin_width] + scale *= energy_bin_width / nu.keV # Convert from [per keV] to [per energy_bin_width] rates_list = [] if modulation: - # This should probably be done in a more efficient way - for time in times: - rates_list.append( - wr.rate_wimp_std( - energy_values, - mw=wimp_mass, - sigma_nucleon=sigma, - t=time, - detection_mechanism="migdal", - migdal_model=migdal_model, - dipole=dipole, + with ProcessPoolExecutor(4) as executor: + for time in times: + rates_list.append( + executor.submit( + wr.rate_wimp_std, + energy_values, + mw=wimp_mass, + sigma_nucleon=sigma, + t=time, + detection_mechanism="migdal", + migdal_model=migdal_model, + dipole=dipole, + ) ) - ) + + rates_list = [future.result() for future in rates_list] else: rates = wr.rate_wimp_std( - energy_values, - mw=wimp_mass, - sigma_nucleon=sigma, - detection_mechanism="migdal", - migdal_model=migdal_model, - dipole=dipole, - ) + energy_values, + mw=wimp_mass, + sigma_nucleon=sigma, + detection_mechanism="migdal", + migdal_model=migdal_model, + dipole=dipole, + ) for _ in times: - rates_list.append(rates) + rates_list.append(rates * scale) RATES = np.array(rates_list) From 30753f467852b4383c0a529de21f537dc650df46 Mon Sep 17 00:00:00 2001 From: Lorenzo Principe <28869147+lorenzomag@users.noreply.github.com> Date: Tue, 13 Aug 2024 00:43:56 +0200 Subject: [PATCH 17/22] Remove duplicated class --- flamedisx/nest/lxe_sources.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/flamedisx/nest/lxe_sources.py b/flamedisx/nest/lxe_sources.py index 8f55a5b6..2404c7b5 100644 --- a/flamedisx/nest/lxe_sources.py +++ b/flamedisx/nest/lxe_sources.py @@ -813,4 +813,4 @@ def get_energy_hist( RATES, [time_bin_edges, energy_bin_edges], axis_names=("time", "energy") ) - return hist + return hist \ No newline at end of file From d85315284f4d30cd6ddd2ae2c794e5dd3ba1087f Mon Sep 17 00:00:00 2001 From: Lorenzo Principe <28869147+lorenzomag@users.noreply.github.com> Date: Tue, 13 Aug 2024 14:44:02 +0200 Subject: [PATCH 18/22] Update XLZD sources --- flamedisx/xlzd/xlzd.py | 44 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 44 insertions(+) diff --git a/flamedisx/xlzd/xlzd.py b/flamedisx/xlzd/xlzd.py index e8754b13..810935d8 100644 --- a/flamedisx/xlzd/xlzd.py +++ b/flamedisx/xlzd/xlzd.py @@ -169,6 +169,7 @@ def __init__( n_energy_bins=800, min_time="2019-09-01T08:28:00", max_time="2020-09-01T08:28:00", + livetime=1.0, n_time_bins=25, modulation=True, **kwargs @@ -187,12 +188,55 @@ def __init__( n_energy_bins=n_energy_bins, min_time=min_time, max_time=max_time, + livetime=livetime, n_time_bins=n_time_bins, modulation=modulation, **kwargs ) +@export +class XLZDMigdalSource(XLZDSource, fd.nest.nestMigdalSource): + + def __init__( + self, + *args, + wimp_mass=40, + sigma=1e-45, + fid_mass=1.0, + min_E=1e-2, + max_E=80.0, + n_energy_bins=800, + min_time="2019-09-01T08:28:00", + max_time="2020-09-01T08:28:00", + livetime=1.0, + n_time_bins=25, + modulation=True, + migdal_model="Cox", + **kwargs + ): + if "detector" not in kwargs: + kwargs["detector"] = "xlzd" + if "configuration" not in kwargs: + kwargs["configuration"] = "80t" + super().__init__( + *args, + wimp_mass=wimp_mass, + sigma=sigma, + fid_mass=fid_mass, + min_E=min_E, + max_E=max_E, + n_energy_bins=n_energy_bins, + min_time=min_time, + max_time=max_time, + livetime=livetime, + n_time_bins=n_time_bins, + modulation=modulation, + migdal_model=migdal_model, + **kwargs + ) + + ## # Background sources ## From 1fc573a25fc6235b857bf9284329b547e1d2d134 Mon Sep 17 00:00:00 2001 From: Lorenzo Principe <28869147+lorenzomag@users.noreply.github.com> Date: Tue, 13 Aug 2024 16:40:34 +0200 Subject: [PATCH 19/22] Fix scaling bug Scaling was only applied when modulation off --- flamedisx/nest/lxe_sources.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/flamedisx/nest/lxe_sources.py b/flamedisx/nest/lxe_sources.py index 2404c7b5..36831609 100644 --- a/flamedisx/nest/lxe_sources.py +++ b/flamedisx/nest/lxe_sources.py @@ -657,7 +657,7 @@ def get_energy_hist( rates_list = [] if modulation: - with ProcessPoolExecutor(4) as executor: + with ProcessPoolExecutor() as executor: for time in times: rates_list.append( executor.submit( @@ -679,7 +679,7 @@ def get_energy_hist( for _ in times: rates_list.append(rates) - RATES = np.array(rates_list) + RATES = np.array(rates_list) * scale hist = mh.Histdd.from_histogram( RATES, [time_bin_edges, energy_bin_edges], axis_names=("time", "energy") @@ -779,7 +779,7 @@ def get_energy_hist( rates_list = [] if modulation: - with ProcessPoolExecutor(4) as executor: + with ProcessPoolExecutor() as executor: for time in times: rates_list.append( executor.submit( @@ -805,9 +805,9 @@ def get_energy_hist( dipole=dipole, ) for _ in times: - rates_list.append(rates * scale) + rates_list.append(rates) - RATES = np.array(rates_list) + RATES = np.array(rates_list) * scale hist = mh.Histdd.from_histogram( RATES, [time_bin_edges, energy_bin_edges], axis_names=("time", "energy") From fb16bca0a7afb57f58ff2627f3d16edf7e174af7 Mon Sep 17 00:00:00 2001 From: Lorenzo Principe <28869147+lorenzomag@users.noreply.github.com> Date: Wed, 14 Aug 2024 12:15:29 +0200 Subject: [PATCH 20/22] Code refactoring and normalisation fix The energy histogram is not normalised in the time dimension. Robert's expalation: I think all you need to do is remove scale = time_bin_width / nu.year # Convert from [per year] to [per time_bin_width] Whilst in principle keeping this and removing / self.n_time_bins would solve things at the level of source.mu_before_efficiencies(), it will mess up the tensor-based calculation of the differential rate, as the thing that is queried should be the spectrum corresponding to a given event time (such that the sum of the spectrum gives the expected counts within that energy range, after scaling by the exposure) --- flamedisx/nest/__init__.py | 1 + flamedisx/nest/get_energy_hist.py | 88 ++++++++++ flamedisx/nest/lxe_blocks/energy_spectrum.py | 2 +- flamedisx/nest/lxe_sources.py | 161 +------------------ 4 files changed, 94 insertions(+), 158 deletions(-) create mode 100644 flamedisx/nest/get_energy_hist.py diff --git a/flamedisx/nest/__init__.py b/flamedisx/nest/__init__.py index 15c530fc..b45a981a 100644 --- a/flamedisx/nest/__init__.py +++ b/flamedisx/nest/__init__.py @@ -1,4 +1,5 @@ from .parameter_calc import * +from .get_energy_hist import * from .lxe_blocks.energy_spectrum import * from .lxe_blocks.secondary_quanta_generation import * diff --git a/flamedisx/nest/get_energy_hist.py b/flamedisx/nest/get_energy_hist.py new file mode 100644 index 00000000..223f77e7 --- /dev/null +++ b/flamedisx/nest/get_energy_hist.py @@ -0,0 +1,88 @@ +from concurrent.futures import ProcessPoolExecutor + +import numericalunits as nu +import numpy as np +import pandas as pd + +import flamedisx as fd +import multihist as mh +import wimprates as wr + + +export, __all__ = fd.exporter() + + +@export +def get_energy_hist( + wimp_mass=40.0, + sigma=1e-45, + min_E=1e-2, + max_E=80.0, + n_energy_bins=800, + min_time="2019-09-01T08:28:00", + max_time="2020-09-01T08:28:00", + n_time_bins=25, + modulation=True, + migdal_model=None, +): + + wimprates_kwargs = dict( + mw=wimp_mass, + sigma_nucleon=sigma, + ) + + if migdal_model is not None: + if migdal_model not in ["Ibe", "Cox", "Cox_dipole"]: + raise ValueError( + "Invalid Migdal model. Choose from 'Ibe', 'Cox', 'Cox_dipole'" + ) + dipole = False + if migdal_model == "Cox_dipole": + migdal_model = "Cox" + dipole = True + + wimprates_kwargs.update( + dict( + detection_mechanism="migdal", + migdal_model=migdal_model, + dipole=dipole, + ) + ) + + energy_bin_edges = np.linspace(min_E, max_E, n_energy_bins + 1) + energy_bin_width = (energy_bin_edges[1] - energy_bin_edges[0]) * nu.keV + energy_values = (energy_bin_edges[1:] + energy_bin_edges[:-1]) / 2 + + time_bin_edges = ( + pd.date_range(min_time, max_time, periods=n_time_bins + 1).to_julian_date() + - 2451545.0 # Convert to J2000 + ) + times = (time_bin_edges[:-1] + time_bin_edges[1:]) / 2 + + scale = ( + energy_bin_width / nu.keV + ) # Convert from [per keV] to [per energy_bin_width] + + rates_list = [] + if modulation: + with ProcessPoolExecutor() as executor: + for time in times: + rates_list.append( + executor.submit( + wr.rate_wimp_std, energy_values, t=time, **wimprates_kwargs + ) + ) + + rates_list = [future.result() for future in rates_list] + else: + rates = wr.rate_wimp_std(energy_values, **wimprates_kwargs) + for _ in times: + rates_list.append(rates) + + RATES = np.array(rates_list) * scale + + hist = mh.Histdd.from_histogram( + RATES, [time_bin_edges, energy_bin_edges], axis_names=("time", "energy") + ) + + return hist diff --git a/flamedisx/nest/lxe_blocks/energy_spectrum.py b/flamedisx/nest/lxe_blocks/energy_spectrum.py index 91cb7c69..e861a95f 100644 --- a/flamedisx/nest/lxe_blocks/energy_spectrum.py +++ b/flamedisx/nest/lxe_blocks/energy_spectrum.py @@ -449,7 +449,7 @@ def clip_j2000_times(self, ts): return np.clip(ts, tbins[0], tbins[-1]) def mu_before_efficiencies(self, **params): - return self.energy_hist.n + return self.energy_hist.n / self.n_time_bins def random_truth(self, n_events, fix_truth=None, **params): """Draw n_events random energies and times from the energy/ diff --git a/flamedisx/nest/lxe_sources.py b/flamedisx/nest/lxe_sources.py index 36831609..612c413b 100644 --- a/flamedisx/nest/lxe_sources.py +++ b/flamedisx/nest/lxe_sources.py @@ -573,7 +573,7 @@ def __init__( max_time="2020-09-01T08:28:00", livetime=1.0, n_time_bins=25, - modulation=True, + modulation=False, **kwargs ): @@ -583,7 +583,7 @@ def __init__( if "detector" not in kwargs: kwargs["detector"] = "default" - self.energy_hist = self.get_energy_hist( + self.energy_hist = fd_nest.get_energy_hist( wimp_mass=wimp_mass, sigma=sigma, min_E=min_E, @@ -606,87 +606,6 @@ def __init__( super().__init__(*args, **kwargs) - @staticmethod - def get_energy_hist( - wimp_mass=40.0, - sigma=1e-45, - min_E=1e-2, - max_E=80.0, - n_energy_bins=800, - min_time="2019-09-01T08:28:00", - max_time="2020-09-01T08:28:00", - n_time_bins=25, - modulation=True, - ): - """ - Function to generate the energy histogram for a given WIMP mass - - Parameters - ---------- - wimp_mass : float - Mass of the WIMP in GeV/c^2 - min_E : float - Minimum energy of the histogram in keV - max_E : float - Maximum energy of the histogram in keV - sigma : float - Cross-section of the WIMP in pb - n_bins_energy : int - Number of bins for the energy histogram - n_bins_time : int - Number of bins for the time histogram. - To total timespan to bin is currently hardcoded to 1 year - modulation : bool - Flag to enable/disable modulation - """ - - energy_bin_edges = np.linspace(min_E, max_E, n_energy_bins + 1) - energy_bin_width = (energy_bin_edges[1] - energy_bin_edges[0]) * nu.keV - energy_values = (energy_bin_edges[1:] + energy_bin_edges[:-1]) / 2 - - time_bin_edges = ( - pd.date_range(min_time, max_time, periods=n_time_bins + 1).to_julian_date() - - 2451545.0 # Convert to J2000 - ) - time_bin_width = wr.j2000_to_datetime(time_bin_edges[1]) - wr.j2000_to_datetime(time_bin_edges[0]) - time_bin_width = time_bin_width.value * nu.ns - times = (time_bin_edges[:-1] + time_bin_edges[1:]) / 2 - - scale = time_bin_width / nu.year # Convert from [per year] to [per time_bin_width] - scale *= energy_bin_width / nu.keV # Convert from [per keV] to [per energy_bin_width] - - rates_list = [] - if modulation: - with ProcessPoolExecutor() as executor: - for time in times: - rates_list.append( - executor.submit( - wr.rate_wimp_std, - energy_values, - mw=wimp_mass, - sigma_nucleon=sigma, - t=time, - ) - ) - - rates_list = [future.result() for future in rates_list] - else: - rates = wr.rate_wimp_std( - energy_values, - mw=wimp_mass, - sigma_nucleon=sigma, - ) - for _ in times: - rates_list.append(rates) - - RATES = np.array(rates_list) * scale - - hist = mh.Histdd.from_histogram( - RATES, [time_bin_edges, energy_bin_edges], axis_names=("time", "energy") - ) - - return hist - @export class nestMigdalSource(nestERSource): @@ -719,7 +638,7 @@ def __init__( if "detector" not in kwargs: kwargs["detector"] = "default" - self.energy_hist = self.get_energy_hist( + self.energy_hist = fd_nest.get_energy_hist( wimp_mass=wimp_mass, sigma=sigma, min_E=min_E, @@ -741,76 +660,4 @@ def __init__( self.array_columns = (("energy_spectrum", len(e_centers)),) - super().__init__(*args, **kwargs) - - @staticmethod - def get_energy_hist( - wimp_mass=40.0, - sigma=1e-45, - min_E=1e-2, - max_E=80.0, - n_energy_bins=800, - min_time="2019-09-01T08:28:00", - max_time="2020-09-01T08:28:00", - n_time_bins=25, - modulation=True, - migdal_model="Cox", - ): - dipole = False - if migdal_model == "Cox_dipole": - migdal_model = "Cox" - dipole = True - - - energy_bin_edges = np.linspace(min_E, max_E, n_energy_bins + 1) - energy_bin_width = (energy_bin_edges[1] - energy_bin_edges[0]) * nu.keV - energy_values = (energy_bin_edges[1:] + energy_bin_edges[:-1]) / 2 - - time_bin_edges = ( - pd.date_range(min_time, max_time, periods=n_time_bins + 1).to_julian_date() - - 2451545.0 # Convert to J2000 - ) - time_bin_width = wr.j2000_to_datetime(time_bin_edges[1]) - wr.j2000_to_datetime(time_bin_edges[0]) - time_bin_width = time_bin_width.value * nu.ns - times = (time_bin_edges[:-1] + time_bin_edges[1:]) / 2 - - scale = time_bin_width / nu.year # Convert from [per year] to [per time_bin_width] - scale *= energy_bin_width / nu.keV # Convert from [per keV] to [per energy_bin_width] - - rates_list = [] - if modulation: - with ProcessPoolExecutor() as executor: - for time in times: - rates_list.append( - executor.submit( - wr.rate_wimp_std, - energy_values, - mw=wimp_mass, - sigma_nucleon=sigma, - t=time, - detection_mechanism="migdal", - migdal_model=migdal_model, - dipole=dipole, - ) - ) - - rates_list = [future.result() for future in rates_list] - else: - rates = wr.rate_wimp_std( - energy_values, - mw=wimp_mass, - sigma_nucleon=sigma, - detection_mechanism="migdal", - migdal_model=migdal_model, - dipole=dipole, - ) - for _ in times: - rates_list.append(rates) - - RATES = np.array(rates_list) * scale - - hist = mh.Histdd.from_histogram( - RATES, [time_bin_edges, energy_bin_edges], axis_names=("time", "energy") - ) - - return hist \ No newline at end of file + super().__init__(*args, **kwargs) \ No newline at end of file From 3b280f5d12fd1ab4ba8f7c309bfa2ad1126bb8de Mon Sep 17 00:00:00 2001 From: Ananthu-Ravindran Date: Wed, 10 Sep 2025 16:10:02 +1000 Subject: [PATCH 21/22] 4 fold --- flamedisx/nest/config/xlzd.ini | 10 +++++----- flamedisx/nest/parameter_calc.py | 11 ++++++----- 2 files changed, 11 insertions(+), 10 deletions(-) diff --git a/flamedisx/nest/config/xlzd.ini b/flamedisx/nest/config/xlzd.ini index d8edac2f..8b219e09 100644 --- a/flamedisx/nest/config/xlzd.ini +++ b/flamedisx/nest/config/xlzd.ini @@ -11,12 +11,12 @@ elife_config = 13000000. ; DUMMY VALUE, DEFINED IN CONSTRUCTOR INSTEAD spe_res_config = 0.38 spe_thr_config = 0.375 -spe_eff_config = 0.999 +spe_eff_config = 1. num_pmts_config = 902 double_pe_fraction_config = 0.2 -coin_level_config = 3 +coin_level_config = 4 gas_gap_config = 10. ; mm g1_gas_config = 0.10 @@ -34,10 +34,10 @@ S1_max_config = 0. ; DUMMY VALUE, NOT USED CURRENTLY S2_min_config = 0. ; DUMMY VALUE, NOT USED CURRENTLY S2_max_config = 0. ; DUMMY VALUE, NOT USED CURRENTLY -cS1_min_config = 3.6 -cS1_max_config = 125. +cS1_min_config = 0. +cS1_max_config = 100. log10_cS2_min_config = 2.5 -log10_cS2_max_config = 5. +log10_cS2_max_config = 4.0 ; Dummy values - either redundant or overriden diff --git a/flamedisx/nest/parameter_calc.py b/flamedisx/nest/parameter_calc.py index 0bfb0168..ae22f944 100644 --- a/flamedisx/nest/parameter_calc.py +++ b/flamedisx/nest/parameter_calc.py @@ -181,16 +181,17 @@ def calculate_s1_mean_mult(spe_res): @export def get_coin_table(coin_level, num_pmts, spe_res, spe_thr, spe_eff, double_pe_fraction): - assert coin_level <= 3, 'This logic will not work well for coincidence levels higher than 3' + assert coin_level <= 4, 'This logic will not work well for coincidence levels higher than 3' + print(f'Running coing level: {coin_level}-fold') coin_dict = dict() coin_table = [] if (coin_level == 0): - for ph_det in np.arange(0, 6): + for ph_det in np.arange(0, 11): coin_table.append(1.) return coin_table - for spike in np.arange(1, 6): + for spike in np.arange(1, 11): numer = 0. denom = 0. i = spike @@ -212,7 +213,7 @@ def get_coin_table(coin_level, num_pmts, spe_res, spe_thr, spe_eff, double_pe_fr belowThresh_percentile = sPE_belowThresh_percentile * (1. - double_pe_fraction) + \ dPE_belowThresh_percentile * double_pe_fraction - for ph_det in np.arange(0, 6): + for ph_det in np.arange(0, 11): spe_eff_mod = spe_eff if (spe_eff_mod < 1.): spe_eff_mod += (1. - spe_eff_mod) / (2. * num_pmts) * ph_det @@ -221,7 +222,7 @@ def get_coin_table(coin_level, num_pmts, spe_res, spe_thr, spe_eff, double_pe_fr p = spe_eff_mod * (1. - belowThresh_percentile) binom_prob = 0 - for spike in np.arange(1, 6): + for spike in np.arange(1, 11): binom_prob += stats.binom.pmf(spike, ph_det, p) * coin_dict[spike] coin_table.append(binom_prob) From 9a658ab1aee5f4449e950efde1e0540feac88dc5 Mon Sep 17 00:00:00 2001 From: Ananthu-Ravindran Date: Tue, 4 Nov 2025 16:37:16 +1100 Subject: [PATCH 22/22] new lce from owen --- flamedisx/nest/config/xlzd.ini | 23 ++++++++-------------- flamedisx/xlzd/xlzd.py | 36 +++++++++++++++++----------------- 2 files changed, 26 insertions(+), 33 deletions(-) diff --git a/flamedisx/nest/config/xlzd.ini b/flamedisx/nest/config/xlzd.ini index a9509ea8..dcf23817 100644 --- a/flamedisx/nest/config/xlzd.ini +++ b/flamedisx/nest/config/xlzd.ini @@ -4,22 +4,21 @@ ; Detector parameters temperature_config = 174.1 ; K pressure_config = 1.79 ; bar -drift_field_config = 100. ; DUMMY VALUE, DEFINED IN CONSTRUCTOR INSTEAD -gas_field_config = 8. ; DUMMY VALUE, DEFINED IN CONSTRUCTOR INSTEAD +drift_field_config = 80. ; DUMMY VALUE, DEFINED IN CONSTRUCTOR INSTEAD +gas_field_config = 6. ; DUMMY VALUE, DEFINED IN CONSTRUCTOR INSTEAD -g1_config = 0.27 ; ; DUMMY VALUE, DEFINED IN CONSTRUCTOR INSTEAD +g1_config = 0.31 ; ; DUMMY VALUE, DEFINED IN CONSTRUCTOR INSTEAD elife_config = 13000000. ; DUMMY VALUE, DEFINED IN CONSTRUCTOR INSTEAD spe_res_config = 0.38 spe_thr_config = 0.375 spe_eff_config = 1. -num_pmts_config = 902 +num_pmts_config = 1155 double_pe_fraction_config = 0.2 coin_level_config = 4 -gas_gap_config = 10. ; mm g1_gas_config = 0.10 s2Fano_config = 2. @@ -30,15 +29,11 @@ S2_noise_config = 0. ; Selection parameters -S1_min_config = 0. ; DUMMY VALUE, NOT USED CURRENTLY -S1_max_config = 0. ; DUMMY VALUE, NOT USED CURRENTLY -S2_min_config = 0. ; DUMMY VALUE, NOT USED CURRENTLY -S2_max_config = 0. ; DUMMY VALUE, NOT USED CURRENTLY -cS1_min_config = 0. -cS1_max_config = 100. -log10_cS2_min_config = 2.5 -log10_cS2_max_config = 4.0 +cS1_min_config = 3. +cS1_max_config = 30. +log10_cS2_min_config = 3. +log10_cS2_max_config = 5.0 ; Dummy values - either redundant or overriden @@ -56,8 +51,6 @@ S1_min_config = 0. S1_max_config = 0. S2_min_config = 0. S2_max_config = 0. -s2_thr_config = 0. -coin_level_config = 3 [80t] diff --git a/flamedisx/xlzd/xlzd.py b/flamedisx/xlzd/xlzd.py index f9060675..836eeb5b 100644 --- a/flamedisx/xlzd/xlzd.py +++ b/flamedisx/xlzd/xlzd.py @@ -111,23 +111,23 @@ def s1_posDependence(self, z): Requires z to be in cm, and in the FV. """ if self.configuration == '80t': - a = 4.93526997e-01 - b = 5.58488459e-04 - c = -1.06051830e-07 - d = -1.02545243e-08 - e = -1.30897496e-11 + a = 5.01480764786202e-01 + b = 1.0987171117870357e-03 + c = 2.6949708579314157e-06 + d = -4.6066555019055335e-09 + e = -3.1658521366562203e-12 elif self.configuration == '60t': - a = 5.00828879e-01 - b = 4.25478465e-04 - c = 2.02171646e-07 - d = -1.31075129e-08 - e = -2.29858516e-11 + a = 5.550088308957059e-01 + b = 7.959774439474241e-04 + c = 1.9863250116287806e-06 + d = -1.4497010098307493e-08 + e = -2.028805722848711e-11 elif self.configuration == '40t': - a = 5.68299226e-01 - b = 7.29787915e-05 - c = -4.37731127e-06 - d = -6.39383595e-08 - e = -1.87387400e-10 + a = 6.324164580843357e-01 + b = 3.980052004436636e-04 + c = 8.151870713156558e-07 + d = -4.238165802951504e-08 + e = -1.0960362072784391e-10 LCE = a + b * z + c * z**2 + d * z**3 + e * z**4 @@ -137,11 +137,11 @@ def add_extra_columns(self, d): super().add_extra_columns(d) if self.configuration == '80t': - LCE_average = 0.471 + LCE_average = 0.4829 elif self.configuration == '60t': - LCE_average = 0.493 + LCE_average = 0.5601 elif self.configuration == '40t': - LCE_average = 0.570 + LCE_average = 0.6532 d['s1_pos_corr'] = self.s1_posDependence(d['z'].values) / LCE_average # normalise to volume-averaged LCE if 's1' in d.columns and 'cs1' not in d.columns: