diff --git a/fuse/plugins/micro_physics/lineage_cluster.py b/fuse/plugins/micro_physics/lineage_cluster.py index b89dd5f5..3ff77f41 100644 --- a/fuse/plugins/micro_physics/lineage_cluster.py +++ b/fuse/plugins/micro_physics/lineage_cluster.py @@ -1,7 +1,6 @@ import numpy as np import strax import straxen -from numba import njit import re import periodictable as pt @@ -48,9 +47,18 @@ class LineageClustering(FuseBasePlugin): # Config options gamma_distance_threshold = straxen.URLConfig( - default=0.9, + default=0.0, type=(int, float), - help="Distance threshold to break lineage for gamma rays [cm]. Default from NEST code", + help="Distance threshold to break lineage for gamma rays [cm]. \ + Do not break if distance is smaller than threshold. \ + Default at 0 means we always break the lineage.", + ) + + brem_distance_threshold = straxen.URLConfig( + default=0.1, + type=(int, float), + help="Distance threshold to break lineage for bremsstrahlung [cm]. \ + Do not break if distance is smaller than threshold.", ) time_threshold = straxen.URLConfig( @@ -59,6 +67,19 @@ class LineageClustering(FuseBasePlugin): help="Time threshold to break the lineage [ns]", ) + classify_ic_as_gamma = straxen.URLConfig( + default=True, + type=bool, + help="Classify internal conversion electrons as gamma particles", + ) + + classify_phot_as_beta = straxen.URLConfig( + default=True, + type=bool, + help="Classify photoabsorption electrons as beta particles \ + (if False, classify as gamma particles)", + ) + def compute(self, geant4_interactions): """ Args: @@ -170,14 +191,21 @@ def build_lineage_for_event(self, event): parent, parent_lineage, self.gamma_distance_threshold, + self.brem_distance_threshold, self.time_threshold, ) if broken_lineage: # The lineage is broken. We can start a new one! running_lineage_index += 1 + tmp_result = start_new_lineage( - particle, tmp_result, i, running_lineage_index + particle, + tmp_result, + i, + running_lineage_index, + self.classify_ic_as_gamma, + self.classify_phot_as_beta, ) else: @@ -187,7 +215,15 @@ def build_lineage_for_event(self, event): else: # Particle without parent. Start a new lineage running_lineage_index += 1 - tmp_result = start_new_lineage(particle, tmp_result, i, running_lineage_index) + + tmp_result = start_new_lineage( + particle, + tmp_result, + i, + running_lineage_index, + self.classify_ic_as_gamma, + self.classify_phot_as_beta, + ) else: # We have seen this particle before. Now evaluate if we have to break the lineage @@ -202,13 +238,20 @@ def build_lineage_for_event(self, event): last_particle_interaction, last_particle_lineage, self.gamma_distance_threshold, + self.brem_distance_threshold, self.time_threshold, ) if broken_lineage: # New lineage! running_lineage_index += 1 + tmp_result = start_new_lineage( - particle, tmp_result, i, running_lineage_index + particle, + tmp_result, + i, + running_lineage_index, + self.classify_ic_as_gamma, + self.classify_phot_as_beta, ) else: @@ -295,10 +338,15 @@ def num_there(s): return any(i.isdigit() for i in s) -def classify_lineage(particle_interaction): +def classify_lineage(particle_interaction, classify_ic_as_gamma, classify_phot_as_beta): """Function to classify a new lineage based on the particle and its parent information.""" + # Excited states of nuclei, decaying electromagnetically + # this will become the lineage of internal conversion electrons + if "[" in particle_interaction["type"]: + return NEST_GAMMA if classify_ic_as_gamma else NEST_BETA + # NR interactions if (particle_interaction["parenttype"] == "neutron") & ( num_there(particle_interaction["type"]) @@ -317,7 +365,7 @@ def classify_lineage(particle_interaction): elif particle_interaction["creaproc"] == "conv": return NEST_BETA elif particle_interaction["creaproc"] == "phot": - return NEST_GAMMA + return NEST_BETA if classify_phot_as_beta else NEST_GAMMA else: # This case should not happen or? Classify it as nontype return NEST_NONE @@ -333,15 +381,8 @@ def classify_lineage(particle_interaction): elif particle_interaction["edproc"] == "conv": return NEST_BETA elif particle_interaction["edproc"] == "phot": - # Not sure about this, but was looking for a way to give beta yields - # to gammas that are coming directly from a radioactive decay - if particle_interaction["creaproc"] == "RadioactiveDecayBase": - return NEST_BETA - # Need this case for custom geant4 inputs... - elif particle_interaction["creaproc"] == "Null": - return NEST_BETA - else: - return NEST_GAMMA + # This is gamma photoabsorption. Return gamma + return NEST_BETA if classify_phot_as_beta else NEST_GAMMA else: # could be rayleigh scattering or something else. Classify it as gamma... return NEST_BETA @@ -351,6 +392,11 @@ def classify_lineage(particle_interaction): particle_interaction["parenttype"] == "none" ): + # If [ in type, it is a nucleus excitation + # we give it a beta for the possible conversion electrons + if "[" in particle_interaction["type"]: + return NEST_BETA + # Alpha particles if particle_interaction["type"] == "alpha": return NEST_ALPHA @@ -369,28 +415,37 @@ def classify_lineage(particle_interaction): return NEST_NONE -@njit() def is_lineage_broken( particle, parent, parent_lineage, gamma_distance_threshold, + brem_distance_threshold, time_threshold, ): """Function to check if the lineage is broken.""" - # In the nest code: Lineage is always broken if the parent is a ion - # But if it's an alpha particle, we want to keep the lineage - break_for_ion = parent_lineage["lineage_type"] == 6 - break_for_ion &= parent["type"] != "alpha" - break_for_ion &= particle["creaproc"] != "eIoni" + if ( + particle["creaproc"] == "RadioactiveDecayBase" + and particle["edproc"] == "RadioactiveDecayBase" + ): + # second step of a decay. We want to split the lineage + return True - if break_for_ion: + # In the nest code: Lineage is always broken if the parent is a ion + # this breaks the lineage for all ions, also for alpha decays (we need it) + # but if it's via an excited nuclear state, we want to keep the lineage + if (num_there(parent["type"])) and ("[" not in parent["type"]): return True # For gamma rays, check the distance between the parent and the particle if particle["type"] == "gamma": + if particle["creaproc"] == "phot" and particle["edproc"] == "phot": + # We do not want to split a photo absorption into two clusters + # The second photo absorption (that we see) could be x rays + return False + # Break the lineage for these transportation gammas # Transportations is a special case. They are not real gammas. # They are just used to transport the energy @@ -400,9 +455,14 @@ def is_lineage_broken( particle_position = np.array([particle["x"], particle["y"], particle["z"]]) parent_position = np.array([parent["x"], parent["y"], parent["z"]]) - distance = np.sqrt(np.sum((parent_position - particle_position) ** 2, axis=0)) + if particle["creaproc"] == "eBrem": + # we do not want to split a bremsstrahlung into two clusters + # if the distance is really small, it is most likely the same interaction + if distance < brem_distance_threshold: + return False + if distance > gamma_distance_threshold: return True @@ -497,9 +557,13 @@ def propagate_mega_type(mega_type): return main_cluster_types -def start_new_lineage(particle, tmp_result, i, running_lineage_index): +def start_new_lineage( + particle, tmp_result, i, running_lineage_index, classify_ic_as_gamma, classify_phot_as_beta +): - lineage_class, lineage_A, lineage_Z = classify_lineage(particle) + lineage_class, lineage_A, lineage_Z = classify_lineage( + particle, classify_ic_as_gamma, classify_phot_as_beta + ) tmp_result[i]["lineage_index"] = running_lineage_index tmp_result[i]["lineage_type"] = lineage_class tmp_result[i]["lineage_A"] = lineage_A diff --git a/fuse/plugins/micro_physics/yields.py b/fuse/plugins/micro_physics/yields.py index c4158a79..991e0857 100644 --- a/fuse/plugins/micro_physics/yields.py +++ b/fuse/plugins/micro_physics/yields.py @@ -2,7 +2,6 @@ import nestpy import strax import straxen -import pickle from ...dtypes import quanta_fields from ...plugin import FuseBasePlugin @@ -19,7 +18,7 @@ class NestYields(FuseBasePlugin): """Plugin that calculates the number of photons, electrons and excitons produced by energy deposit using nestpy.""" - __version__ = "0.2.2" + __version__ = "0.2.3" depends_on = ("interactions_in_roi", "electric_field_values") provides = "quanta" @@ -29,6 +28,37 @@ class NestYields(FuseBasePlugin): save_when = strax.SaveWhen.TARGET + return_yields_only = straxen.URLConfig( + default=False, + type=bool, + help="Set to True to return the yields model output directly instead of the \ + calculated actual quanta with NEST getQuanta function. Only for testing purposes.", + ) + + nest_width_parameters = straxen.URLConfig( + default={}, + type=dict, + help="Set to modify default NEST NRERWidthParameters to match recombination fluctuations. \ + From NEST code https://github.com/NESTCollaboration/nest/blob/v2.4.0/src/NEST.cpp \ + and NEST paper https://arxiv.org/abs/2211.10726 \ + See self.get_nest_width_parameters() for the options and default values. \ + Example use: {'fano_ER': -0.0015, 'A_ER': 0.096452}", + ) + + nest_er_yields_parameters = straxen.URLConfig( + default=[-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0], + type=list, + help="Set to modify default NEST ER yields parameters. Use -1 to keep default value. \ + From NEST code https://github.com/NESTCollaboration/nest/blob/v2.4.0/src/NEST.cpp \ + Used in the calcuations of BetaYieldsGR.", + ) + + fix_gamma_yield_field = straxen.URLConfig( + default=-1.0, + help="Field in V/cm to use for NEST gamma yield calculation. Only used if set to > 0.", + type=float, + ) + def setup(self): super().setup() @@ -39,7 +69,43 @@ def setup(self): else: self.log.debug("Generating random numbers with seed pulled from OS") - self.quanta_from_NEST = np.vectorize(self._quanta_from_NEST) + self.nc = nestpy.NESTcalc(nestpy.VDetector()) + self.vectorized_get_quanta = np.vectorize(self.get_quanta) + self.updated_nest_width_parameters = self.update_nest_width_parameters() + + def update_nest_width_parameters(self): + + # Get the default NEST NRERWidthsParam + free_parameters = self.nc.default_NRERWidthsParam + + # Map the parameters names to the index in the free_parameters list + parameters_key_map = { + "fano_ions_NR": 0, # Fano factor for NR Ions (default 0.4) + "fano_excitons_NR": 1, # Fano factor for NR Excitons (default 0.4) + "A_NR": 2, # A' - Amplitude for Recombinnation NR (default 0.04) + "xi_NR": 3, # ξ - Center for Recombination NR (default 0.50) + "omega_NR": 4, # ω - Width for Recombination NR (default 0.19) + "skewness_NR": 5, # Skewness for NR (default 2.25) + "fano_ER": 6, # Multiplier for Fano ER, if<0 field dep, else constant. (default 1) + "A_ER": 7, # A - Amplitude for Recombination ER, field dependent (default 0.096452) + "omega_ER": 8, # ω - Width for Recombination ER (default 0.205) + "xi_ER": 9, # ξ - Center for Recombination ER (default 0.45) + "alpha_skewness_ER": 10, # Skewness for ER (default -0.2) + } + + if self.nest_width_parameters is not None: + for key, value in self.nest_width_parameters.items(): + if key not in parameters_key_map: + raise ValueError( + f"Unknown NEST width parameter {key}.\ + Available parameters: {parameters_key_map.keys()}" + ) + self.log.debug(f"Updating NEST width parameter {key} to {value}") + free_parameters[parameters_key_map[key]] = value + + self.log.debug(f"Using NEST width parameters: {free_parameters}") + + return free_parameters def compute(self, interactions_in_roi): @@ -60,7 +126,15 @@ def compute(self, interactions_in_roi): # Generate quanta: if len(interactions_in_roi) > 0: - photons, electrons, excitons = self.get_quanta(interactions_in_roi) + photons, electrons, excitons = self.vectorized_get_quanta( + interactions_in_roi["ed"], + interactions_in_roi["nestid"], + interactions_in_roi["e_field"], + interactions_in_roi["A"], + interactions_in_roi["Z"], + interactions_in_roi["create_S2"], + interactions_in_roi["xe_density"], + ) result["photons"] = photons result["electrons"] = electrons result["excitons"] = excitons @@ -74,50 +148,18 @@ def compute(self, interactions_in_roi): return result - def get_quanta(self, interactions_in_roi): - - photons, electrons, excitons = self.quanta_from_NEST( - interactions_in_roi["ed"], - interactions_in_roi["nestid"], - interactions_in_roi["e_field"], - interactions_in_roi["A"], - interactions_in_roi["Z"], - interactions_in_roi["create_S2"], - log=self.log, - density=interactions_in_roi["xe_density"], - ) + def get_quanta(self, en, model, e_field, A, Z, create_s2, density): + """Function to get quanta for given parameters using NEST.""" - return photons, electrons, excitons + yields_result = self.get_yields_from_NEST(en, model, e_field, A, Z, density) - @staticmethod - def _quanta_from_NEST(en, model, e_field, A, Z, create_s2, log, **kwargs): + return self.process_yields(yields_result, create_s2) + + def get_yields_from_NEST(self, en, model, e_field, A, Z, density): """Function which uses NEST to yield photons and electrons for a given - set of parameters. - - Note: - In case the energy deposit is outside of the range of NEST a -1 - is returned. - Args: - en (numpy.array): Energy deposit of the interaction [keV] - model (numpy.array): Nest Id for qunata generation (integers) - e_field (numpy.array): Field value in the interaction site [V/cm] - A (numpy.array): Atomic mass number - Z (numpy.array): Atomic number - create_s2 (bool): Specifies if S2 can be produced by interaction, - in this case electrons are generated. - kwargs: Additional keyword arguments which can be taken by - GetYields e.g. density. - Returns: - photons (numpy.array): Number of generated photons - electrons (numpy.array): Number of generated electrons - excitons (numpy.array): Number of generated excitons - """ - nc = nestpy.NESTcalc(nestpy.VDetector()) - - # Fix for Kr83m events. - # Energies have to be very close to 32.1 keV or 9.4 keV - # See: https://github.com/NESTCollaboration/nest/blob/master/src/NEST.cpp#L567 - # and: https://github.com/NESTCollaboration/nest/blob/master/src/NEST.cpp#L585 + set of parameters.""" + + # Fix for Kr83m events max_allowed_energy_difference = 1 # keV if model == 11: if abs(en - 32.1) < max_allowed_energy_difference: @@ -125,123 +167,64 @@ def _quanta_from_NEST(en, model, e_field, A, Z, create_s2, log, **kwargs): if abs(en - 9.4) < max_allowed_energy_difference: en = 9.4 - # Some addition taken from - # https://github.com/NESTCollaboration/nestpy/blob/e82c71f864d7362fee87989ed642cd875845ae3e/src/nestpy/helpers.py#L94-L100 + # Some additions taken from NEST code if model == 0 and en > 2e2: - log.warning( - f"Energy deposition of {en} keV beyond NEST validity " "for NR model of 200 keV" + self.log.warning( + f"Energy deposition of {en} keV beyond NEST validity for NR model of 200 keV" ) - if model == 7 and en > 3e3: - log.warning( - f"Energy deposition of {en} keV beyond NEST validity " "for gamma model of 3 MeV" + self.log.warning( + f"Energy deposition of {en} keV beyond NEST validity for gamma model of 3 MeV" ) - if model == 8 and en > 3e3: - log.warning( - f"Energy deposition of {en} keV beyond NEST validity " "for beta model of 3 MeV" + self.log.warning( + f"Energy deposition of {en} keV beyond NEST validity for beta model of 3 MeV" + ) + + if model == 7 and self.fix_gamma_yield_field > 0: + e_field = self.fix_gamma_yield_field + + if e_field < 0: + raise ValueError( + f"Negative electric field {e_field} V/cm not allowed. \ + (no error will be raised by NEST)." ) - y = nc.GetYields( + yields_result = self.nc.GetYields( interaction=nestpy.INTERACTION_TYPE(model), energy=en, drift_field=e_field, A=A, Z=Z, - **kwargs, + density=density, + ERYieldsParam=self.nest_er_yields_parameters, ) - event_quanta = nc.GetQuanta(y) # Density argument is not use in function... - - photons = event_quanta.photons - excitons = event_quanta.excitons - electrons = 0 - if create_s2: - electrons = event_quanta.electrons - - return photons, electrons, excitons + return yields_result + def process_yields(self, yields_result, create_s2): + """Process the yields with NEST to get actual quanta.""" -@export -class BetaYields(NestYields): - """Plugin that calculates the number of photons, electrons and excitons - produced by energy deposit using nestpy.""" - - depends_on = ("interactions_in_roi", "electric_field_values") - provides = "quanta" - data_kind = "interactions_in_roi" - - # Config options - - use_recombination_fluctuation = straxen.URLConfig( - default=True, - type=bool, - help="Turn on or off the recombination fluctuation for beta interactions.", - ) - - recombination_fluctuation_std_factor = straxen.URLConfig( - default=3, - type=(int, float), - help="A factor that is defined to guess the recombination fluctuation", - ) - - beta_quanta_spline = straxen.URLConfig( - default=None, - help="Path to function that gives n_ph and n_e for a given energy, \ - calculated from beta spectrum. The function should be a pickle file.", - ) - - def get_quanta(self, interactions_in_roi): - - # for the non beta interactions we use the nestpy yields - photons, electrons, excitons = self.quanta_from_NEST( - interactions_in_roi["ed"], - interactions_in_roi["nestid"], - interactions_in_roi["e_field"], - interactions_in_roi["A"], - interactions_in_roi["Z"], - interactions_in_roi["create_S2"], - log=self.log, - density=interactions_in_roi["xe_density"], + # Density argument is not used in function... + event_quanta = self.nc.GetQuanta( + yields_result, free_parameters=self.updated_nest_width_parameters ) - mask_beta = interactions_in_roi["nestid"] == 8 + excitons = event_quanta.excitons + photons = event_quanta.photons + electrons = event_quanta.electrons - # now for the beta interactions we use the beta yields - photons_beta, electrons_beta = self.quanta_from_spline( - interactions_in_roi["ed"][mask_beta], - ) + # Only for testing purposes, return the yields directly + if self.return_yields_only: + photons = yields_result.PhotonYield + electrons = yields_result.ElectronYield - photons[mask_beta] = photons_beta - electrons[mask_beta] = electrons_beta + # If we don't want to create S2, set electrons to 0 + if not create_s2: + electrons = 0 return photons, electrons, excitons - def quanta_from_spline(self, energy): - - with open(self.beta_quanta_spline, "rb") as f: - cs1_poly, cs2_poly = pickle.load(f) - - beta_photons = cs1_poly(energy) - beta_electrons = cs2_poly(energy) - - if self.use_recombination_fluctuation: - - factor = self.recombination_fluctuation_std_factor - rf = self.rng.normal(0, energy * factor, len(energy)) - beta_photons += rf - beta_electrons -= rf - - # make integer quanta - beta_photons = np.round(beta_photons).astype(int) - beta_electrons = np.round(beta_electrons).astype(int) - - # make sure we don't have negative quanta, so clip at 0 - beta_photons = np.clip(beta_photons, 0, np.inf) - beta_electrons = np.clip(beta_electrons, 0, np.inf) - - return beta_photons, beta_electrons - @export class BBFYields(FuseBasePlugin): diff --git a/tests/test_MicroPhysics_alt_Yields.py b/tests/test_MicroPhysics_alt_Yields.py deleted file mode 100644 index 6e346b4f..00000000 --- a/tests/test_MicroPhysics_alt_Yields.py +++ /dev/null @@ -1,76 +0,0 @@ -import os -import shutil -import unittest -import tempfile -import timeout_decorator -import fuse -import straxen -from _utils import test_root_file_name -import pickle - -TIMEOUT = 60 - - -def yields_dummy_func(x): - """Dummy function that returns two values for the n_photon and n_electron. - - To be used as a dummy function for the BetaYields plugin. Needs to - be defined outside the test class to be picklable. - """ - return 20000 - - -class TestAlternativeYields(unittest.TestCase): - @classmethod - def setUpClass(cls): - cls.temp_dir = tempfile.TemporaryDirectory() - - cls.test_context = fuse.context.full_chain_context( - output_folder=cls.temp_dir.name, run_without_proper_corrections=True - ) - - cls.test_context.set_config( - { - "path": cls.temp_dir.name, - "file_name": test_root_file_name, - "entry_stop": 25, - } - ) - - cls.run_number = "TestRun_00000" - - @classmethod - def tearDownClass(cls): - cls.temp_dir.cleanup() - - def setUp(self): - downloader = straxen.MongoDownloader(store_files_at=(self.temp_dir.name,)) - downloader.download_single(test_root_file_name, human_readable_file_name=True) - - assert os.path.exists(os.path.join(self.temp_dir.name, test_root_file_name)) - - def tearDown(self): - shutil.rmtree(self.temp_dir.name) - os.makedirs(self.temp_dir.name) - - @timeout_decorator.timeout(TIMEOUT, exception_message="BetaYields timed out") - def test_BetaYields(self): - self.test_context.register(fuse.BetaYields) - - # Make a dummy pkl file, a function that returns two values - # one for the photon yield and one for the electron yield - # as a function of energy - spline_func_name = os.path.join(self.temp_dir.name, "beta_quanta_spline.pkl") - with open(spline_func_name, "wb") as f: - pickle.dump((yields_dummy_func, yields_dummy_func), f) - - self.test_context.set_config({"beta_quanta_spline": spline_func_name}) - - # Make the plugin - self.test_context.make(self.run_number, "quanta") - - # Add tests for the BBF yields later - - -if __name__ == "__main__": - unittest.main()