From 89c26ba7b801a06a3c6c4c3e56fedf13a8b04725 Mon Sep 17 00:00:00 2001 From: Jianyang Qi Date: Sun, 5 Jun 2022 08:49:47 -0700 Subject: [PATCH] XY-dependent SE gain and extraction efficiency (#363) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * Add extraction efficiency and garfield gas gap map * Update load_resource.py * got rid of the xy confinement * s2_correction_map fix * Update load_resource.py * Update load_resource.py * Update load_resource.py * draw se gain from map * make it work without an SE gain map * make compatible with previous wfsim * make compatible with previous wfsim * fix typo * fix conditional issue with se_gain_map * removed print * fix conditional for garfield stuff * remove unused config * change directory * change directory * add sr0 to filenames, change fname check * interpolate between excitation hists * fix bool * add print statements for testing * add print statements to test * add print statements for test * fix divide by zero error * get rid of print statements * update straxen * add param nph to docstring * add comments, copy things to the 1T config * Test without 1T additions * add g2_mean to 1T for test * comment * add _gg to s2_luminescence when using garfield_gas_gap * Add _gg to s2_luminescence when using garfield_gas_gap * allow for constant extraction efficiency * refactor * fix tests * tweak 1 * more debug * fix it more and more and more * fix this mess * one more test * add load resource * add gitignor * fix fail fast * what is the issue? * Update tests/test_wfsim.py * Update wfsim/core/s2.py * Fixing extraction efficiency calculation when no SE gain map is used * Correcting extraction efficiency when SE gain is used from fax config * add print statements to see why tests fail * add another print * flatten correction map from mc * get rid of print statements * change sampling index Co-authored-by: Melih Kara Co-authored-by: Diego Ramírez García Co-authored-by: Joran Angevaare Co-authored-by: Andrii Terliuk --- .github/workflows/pytest.yml | 7 +- .gitignore | 4 +- tests/test_contexts.py | 18 +---- tests/test_load_resource.py | 6 +- tests/test_wfsim.py | 113 ++++++++++++++++++---------- wfsim/core/rawdata.py | 4 +- wfsim/core/s2.py | 139 +++++++++++++++++++++++++++++------ wfsim/load_resource.py | 49 +++++++++--- wfsim/pax_interface.py | 4 +- 9 files changed, 247 insertions(+), 97 deletions(-) diff --git a/.github/workflows/pytest.yml b/.github/workflows/pytest.yml index 01547753..86d700e1 100644 --- a/.github/workflows/pytest.yml +++ b/.github/workflows/pytest.yml @@ -18,9 +18,6 @@ on: release: types: [created] pull_request: - branches: - - master - - stable push: branches: - master @@ -30,7 +27,7 @@ jobs: name: "${{ matrix.test }}_py${{ matrix.python-version }}" runs-on: ubuntu-latest strategy: - fail-fast: True + fail-fast: False matrix: python-version: [3.8, 3.9, "3.10"] test: ['coveralls', 'pytest', 'pytest_no_database'] @@ -84,7 +81,7 @@ jobs: env: NUMBA_DISABLE_JIT: 1 GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} - # We need to check if we had access to the secrets, otherwise coveralls + # We need to check if we have access to the secrets, otherwise coveralls # will yield a low coverage because of the lack of interfacing with the # database. HAVE_ACCESS_TO_SECTETS: ${{ secrets.RUNDB_API_URL }} diff --git a/.gitignore b/.gitignore index cce8bb6c..f6f7a8a1 100644 --- a/.gitignore +++ b/.gitignore @@ -21,6 +21,9 @@ custom_data test_input_data *.zip last_bootstrax_exception.txt +.xenon_config +test.root +pytestdebug.log # cProfile output *.prof @@ -58,4 +61,3 @@ docs/tutorials/* # coverage .coverage - diff --git a/tests/test_contexts.py b/tests/test_contexts.py index 9f344770..3a77e037 100644 --- a/tests/test_contexts.py +++ b/tests/test_contexts.py @@ -1,8 +1,9 @@ import strax import straxen import wfsim +from unittest import skipIf - +@skipIf(not straxen.utilix_is_configured(), 'utilix is not configured') def test_nt_context(register=None, context=None): """ Test a context if it is properly setup. To this end, we perform a @@ -13,9 +14,6 @@ def test_nt_context(register=None, context=None): :param context: Test with some other context than the nT simulation context """ - if not straxen.utilix_is_configured(): - return - if context is None: context = straxen.contexts.xenonnt_simulation(cmt_run_id_sim='010000', cmt_version='global_ONLINE') assert isinstance(context, strax.Context), f'{context} is not a context' @@ -25,15 +23,3 @@ def test_nt_context(register=None, context=None): # Search all plugins for the time field (each should have one) context.search_field('time') - - -# def test_mc_chain(): -# test_nt_context(wfsim.RawRecordsFromMcChain) - - -# def test_fax_nveto(): -# test_nt_context(wfsim.RawRecordsFromFaxnVeto) - - -# def test_1t_context(): -# test_nt_context(context=straxen.contexts.xenon1t_simulation()) diff --git a/tests/test_load_resource.py b/tests/test_load_resource.py index e5054c98..c619f62d 100644 --- a/tests/test_load_resource.py +++ b/tests/test_load_resource.py @@ -1,3 +1,4 @@ +import straxen from wfsim.load_resource import load_config @@ -11,6 +12,8 @@ def test_load_1t(): "enable_electron_afterpulses": True, "enable_noise": False, "field_distortion_on": True, + "g2_mean": 32.3, + 's2_time_model': 's2_time_spread around zero', } result = load_config(config) return result, config @@ -36,7 +39,8 @@ def test_load_nt(): "s2_pattern_map": ["constant dummy", 30e-5, [494,]], "s2_correction_map": ["constant dummy", 1, []], "field_dependencies_map": ["constant dummy", 1, []], - "gains": [1 for i in range(494)], + "gains": [1 for _ in range(straxen.n_tpc_pmts)], + "se_gain_map": ["constant dummy", 1, []], } result = load_config(config) return result, config diff --git a/tests/test_wfsim.py b/tests/test_wfsim.py index 84777cf6..5b8b6713 100644 --- a/tests/test_wfsim.py +++ b/tests/test_wfsim.py @@ -7,7 +7,8 @@ import straxen import wfsim -from .test_load_resource import test_load_nt +from .test_load_resource import test_load_nt, test_load_1t +from unittest import skipIf logging.basicConfig( level=logging.DEBUG, @@ -30,10 +31,12 @@ def test_sim_1T(): """Test the 1T simulator (should always work with the publicly available files)""" with tempfile.TemporaryDirectory() as tempdir: log.debug(f'Working in {tempdir}') - testing_config_1T = dict( + _, conf_1t = test_load_1t() + testing_config_1t = dict( hev_gain_model=("1T_to_pe_placeholder", True), gain_model=("1T_to_pe_placeholder", True), - gain_model_mc=("1T_to_pe_placeholder", True),) + gain_model_mc=("1T_to_pe_placeholder", True), + ) st = strax.Context( storage=tempdir, @@ -41,16 +44,18 @@ def test_sim_1T(): nchunk=2, event_rate=1, chunk_size=1, detector='XENON1T', fax_config=('https://raw.githubusercontent.com/XENONnT/strax_auxiliary_files' - '/36d352580b328ff057b1588b8af8c9a6ed8ae704/sim_files/fax_config_1t.json'), # noqa + '/a5b92102505d6d0bfcdb563b6117bd4040a93435/sim_files/fax_config_1t.json'), # noqa fax_config_override=dict( url_base=("https://raw.githubusercontent.com/XENONnT/strax_auxiliary_files" - "/36d352580b328ff057b1588b8af8c9a6ed8ae704/sim_files/"), ), + "/a5b92102505d6d0bfcdb563b6117bd4040a93435/sim_files/"), + **conf_1t, + ), **straxen.contexts.x1t_common_config), **straxen.contexts.x1t_context_config, ) st.register(wfsim.RawRecordsFromFax1T) - log.debug(f'Setting testing config {testing_config_1T}') - st.set_config(testing_config_1T) + log.debug(f'Setting testing config {testing_config_1t}') + st.set_config(testing_config_1t) log.debug(f'Getting raw-records') rr = st.get_array(run_id, 'raw_records') @@ -59,18 +64,19 @@ def test_sim_1T(): log.info(f'All done') -def test_sim_nT_basics(): +@skipIf(not straxen.utilix_is_configured(), 'utilix is not configured') +def test_sim_nt_basics(): """Test the nT simulator. Uses basic config so complicated steps are skipped. So this will test the simple s1 model and the simple s2 model""" with tempfile.TemporaryDirectory() as tempdir: log.debug(f'Working in {tempdir}') conf = copy.deepcopy(straxen.contexts.xnt_common_config) + resource, conf_override = test_load_nt() conf['gain_model'] = ("to_pe_placeholder", True) conf['gain_model_mc'] = ("to_pe_placeholder", True) conf['hev_gain_model'] = ("to_pe_placeholder", True) conf['hit_min_amplitude'] = 'pmt_commissioning_initial' - resource, conf_override = test_load_nt() # The SPE table in this package is for a single channel # We generate the full SPE file for testing here @@ -101,32 +107,27 @@ def test_sim_nT_basics(): log.info(f'All done') -def test_sim_nT_advanced(): - """Test the nT simulator. Works only if one has access to the XENONnT databases. - Clone the repo to dali and type 'pytest' to run. The first run will test simple s1, - garfield s2 and noise/afterpulses. The second run will test the s1 spline model""" - - if not straxen.utilix_is_configured(): - log.warning(f"Utilix is not configured, skipping database-requiring tests!") - return - +@skipIf(not straxen.utilix_is_configured(), 'utilix is not configured') +def test_sim_nt_advanced( + config = None +): + """ + Test the nT simulator. Works only if one has access to the XENONnT databases. + Clone the repo to dali and type 'pytest' to run. + """ with tempfile.TemporaryDirectory() as tempdir: log.debug(f'Working in {tempdir}') + st = straxen.contexts.xenonnt_simulation(cmt_run_id_sim='010000', cmt_version='global_ONLINE', - _config_overlap={},) - st.set_config(dict(gain_model_mc=("to_pe_placeholder", True), - gain_model=("to_pe_placeholder", True), - hit_min_amplitude='pmt_commissioning_initial' - )) + _config_overlap={}, + fax_config='fax_config_nt_sr0_v0.json' + ) st.set_config(dict(nchunk=1, event_rate=1, chunk_size=2,)) - st.set_config({'fax_config_override': dict(s2_luminescence_model='simple', - s2_time_model="s2_time_spread around zero", - s1_lce_correction_map='XENONnT_s1_xyz_LCE_corrected_qes_MCva43fa9b_wires.json.gz', - s1_time_spline='XENONnT_s1_proponly_va43fa9b_wires_20200625.json.gz', - s1_model_type='optical_propagation+simple',)}) - + if config is not None: + log.warning(f'Update config with {config}') + st.set_config(config) log.debug(f'Getting raw-records') rr = st.get_array(run_id, 'raw_records') log.debug(f'Getting peaks') @@ -135,13 +136,52 @@ def test_sim_nT_advanced(): log.info(f'All done') +@skipIf(not straxen.utilix_is_configured(), 'utilix is not configured') +def test_nt_advanced_alt_s2_model(): + config = dict( + fax_config_override=dict( + s2_luminescence_model='simple', + s2_time_model="s2_time_spread around zero", + s1_lce_correction_map='XENONnT_s1_xyz_LCE_corrected_qes_MCva43fa9b_wires.json.gz', + s1_time_spline='XENONnT_s1_proponly_va43fa9b_wires_20200625.json.gz', + s1_model_type='optical_propagation+simple', + ) + ) + test_sim_nt_advanced(config) + + +@skipIf(not straxen.utilix_is_configured(), 'utilix is not configured') +def test_nt_advanced_garfield(): + config = dict( + fax_config_override=dict( + s2_luminescence_model='garfield', + s2_correction_map=False, + s2_time_model="s2_time_spread around zero", + s1_lce_correction_map='XENONnT_s1_xyz_LCE_corrected_qes_MCva43fa9b_wires.json.gz', + s1_time_spline='XENONnT_s1_proponly_va43fa9b_wires_20200625.json.gz', + s1_model_type='optical_propagation+simple', + ) + ) + test_sim_nt_advanced(config) + +@skipIf(not straxen.utilix_is_configured(), 'utilix is not configured') +def test_nt_advanced_gas_gap_garfield(): + config = dict( + fax_config_override=dict( + s2_luminescence_model='garfield_gas_gap', + s2_correction_map="XENONnT_s2_xy_map_v4_210503_mlp_3_in_1_iterated.json", + s2_luminescence_gg= "garfield_timing_map_gas_gap_sr0.npy", + garfield_gas_gap_map= "garfield_gas_gap_map_sr0.json", + se_gain_map="XENONnT_se_xy_map_v1_mlp.json", + ) + ) + test_sim_nt_advanced(config) + + +@skipIf(not straxen.utilix_is_configured(), 'utilix is not configured') def test_sim_mc_chain(): """Test the nT simulator by Geant4 output file""" - if not straxen.utilix_is_configured(): - log.warning(f"Utilix is not configured, skipping database-requiring tests!") - return - with tempfile.TemporaryDirectory() as tempdir: log.debug(f'Working in {tempdir}') @@ -154,10 +194,7 @@ def test_sim_mc_chain(): st = straxen.contexts.xenonnt_simulation(cmt_run_id_sim='010000', cmt_version='global_ONLINE', _config_overlap={},) - st.set_config(dict(gain_model_mc=("to_pe_placeholder", True), - gain_model=("to_pe_placeholder", True), - gain_model_nv=("adc_nv", True), - hit_min_amplitude='pmt_commissioning_initial', + st.set_config(dict(gain_model_nv=("adc_nv", True), )) epix_config = {'cut_by_eventid': True, 'debug': True, 'source_rate': 0, 'micro_separation_time': 10., @@ -175,8 +212,6 @@ def test_sim_mc_chain(): fax_config_override=dict( s1_model_type='nest', s2_time_model="s2_time_spread around zero", - url_base='https://raw.githubusercontent.com/XENONnT/private_nt_aux_files/master/sim_files', - s1_lce_correction_map=["constant dummy", 1, []], enable_electron_afterpulses=False), epix_config=epix_config, fax_config_nveto='fax_config_nt_nveto.json', diff --git a/wfsim/core/rawdata.py b/wfsim/core/rawdata.py index 62ca73f4..e51aa8ed 100644 --- a/wfsim/core/rawdata.py +++ b/wfsim/core/rawdata.py @@ -377,6 +377,8 @@ def get_mean_xy_electron(self, peak_type, instruction, tb): _, xy_tmp = self.pulses['s2'].field_distortion_comsol(instruction['x'], instruction['y'], instruction['z'], self.resource) elif self.config.get('field_distortion_model', "none") == 'inverse_fdc': _, xy_tmp = self.pulses['s2'].inverse_field_distortion_correction(instruction['x'], instruction['y'], instruction['z'], self.resource) + else: + raise ValueError(f'{self.config.get("field_distortion_model", "none")} is invalid!') tb['x_mean_electron'] = np.mean(xy_tmp.T[0]) tb['y_mean_electron'] = np.mean(xy_tmp.T[1]) else: @@ -455,7 +457,7 @@ def digitizer_saturation(data, channel_mask): @export class RawDataOptical(RawData): - def __init__(self, config, channels=[], timings=[]): + def __init__(self, config, channels=tuple(), timings=tuple()): self.config = config self.pulses = dict( s1=Pulse(config), diff --git a/wfsim/core/s2.py b/wfsim/core/s2.py index b0824d16..e73538ea 100644 --- a/wfsim/core/s2.py +++ b/wfsim/core/s2.py @@ -2,7 +2,6 @@ from numba import njit import numpy as np from scipy.stats import skewnorm -from scipy import interpolate from strax import exporter from .pulse import Pulse from .. import units @@ -97,6 +96,7 @@ def __call__(self, instruction): sc_gain = self.get_s2_light_yield(positions=positions, config=self.config, resource=self.resource) + n_photons_per_xy, n_photons_per_ele, self._electron_timings = self.get_n_photons(t=t, @@ -138,21 +138,20 @@ def __call__(self, instruction): def get_avg_drift_velocity(z, xy, config, resource): """Calculate s2 drift time mean and spread - :param positions: 1d array of z (floats) + :param z: 1d array of z (floats) :param xy: 2d array of xy positions (floats) :param config: dict with wfsim config :param resource: instance of the resource class returns array of floats corresponding to average drift velocities from given point to the gate """ - drift_v_LXe=None if config['enable_field_dependencies']['drift_speed_map']: drift_v_LXe = resource.field_dependencies_map(z, xy, map_name='drift_speed_map') # mm/µs drift_v_LXe *= 1e-4 # cm/ns drift_v_LXe *= resource.drift_velocity_scaling else: - drift_v_LXe=config['drift_velocity_liquid'] - return(drift_v_LXe) + drift_v_LXe = config['drift_velocity_liquid'] + return drift_v_LXe @staticmethod def get_s2_drift_time_params(z_int, xy_int, config, resource): @@ -178,7 +177,7 @@ def get_s2_drift_time_params(z_int, xy_int, config, resource): drift_time_spread = np.sqrt(2 * diffusion_constant_longitudinal * drift_time_mean) drift_time_spread /= drift_velocity_liquid return drift_time_mean, drift_time_spread - + @staticmethod def get_s2_light_yield(positions, config, resource): """Calculate s2 light yield... @@ -189,19 +188,24 @@ def get_s2_light_yield(positions, config, resource): returns array of floats (mean expectation) """ - sc_gain = resource.s2_correction_map(positions) - # depending on if you use the data driven or mc pattern map for light yield for S2 - # the shape of n_photon_hits will change. Mc needs a squeeze + + if config.get('se_gain_from_map', False): + sc_gain = resource.se_gain_map(positions) + else: + # calculate it from MC pattern map directly if no "se_gain_map" is given + sc_gain = resource.s2_correction_map(positions) + sc_gain *= config['s2_secondary_sc_gain'] + + # depending on if you use the data driven or mc pattern map for light yield for S2 + # the shape of n_photon_hits will change. Mc needs a squeeze if len(sc_gain.shape) != 1: sc_gain=np.squeeze(sc_gain, axis=-1) # sc gain should has the unit of pe / electron, here we divide 1 + dpe to get nphoton / electron sc_gain /= 1 + config['p_double_pe_emision'] - sc_gain *= config['s2_secondary_sc_gain'] # data driven map contains nan, will be set to 0 here sc_gain[np.isnan(sc_gain)] = 0 - return sc_gain @staticmethod @@ -212,17 +216,32 @@ def get_electron_yield(n_electron, xy_int, z_int, config, resource): :param xy_int: 2d array of xy interaction positions (floats) :param z_int: 1d array of floats with the z interaction positions (floats) :param config: dict with wfsim config + :param resource: instance of the resource class returns 1d array ints with number of electrons """ # Average drift time of the electrons drift_time_mean, drift_time_spread = S2.get_s2_drift_time_params(z_int, xy_int, config, resource) - + # extraction efficiency in LXe/GXe interface + if config.get('ext_eff_from_map', False): + # Extraction efficiency is g2(x,y)/SE_gain(x,y) + rel_s2_cor=resource.s2_correction_map(xy_int) + #doesn't always need to be flattened, but if s2_correction_map = False, then map is made from MC + rel_s2_cor = rel_s2_cor.flatten() + + if config.get('se_gain_from_map', False): + se_gains=resource.se_gain_map(xy_int) + else: + # is in get_s2_light_yield map is scaled according to relative s2 correction + # we also need to do it here to have consistent g2 + se_gains=rel_s2_cor*config['s2_secondary_sc_gain'] + cy = config['g2_mean']*rel_s2_cor/se_gains + else: + cy = config['electron_extraction_yield'] # Absorb electrons during the drift electron_lifetime_correction = np.exp(- 1 * drift_time_mean / config['electron_lifetime_liquid']) - cy = config['electron_extraction_yield'] * electron_lifetime_correction - + cy*=electron_lifetime_correction # Remove electrons in insensitive volume if config['enable_field_dependencies']['survival_probability_map']: p_surv = resource.field_dependencies_map(z_int, xy_int, map_name='survival_probability_map') @@ -230,7 +249,9 @@ def get_electron_yield(n_electron, xy_int, z_int, config, resource): # FIXME: this is necessary due to map artefacts, such as negative or values >1 p_surv=np.clip(p_surv, a_min = 0, a_max = 1) cy *= p_surv + n_electron = np.random.binomial(n=n_electron, p=cy) + return n_electron @staticmethod @@ -271,7 +292,7 @@ def get_n_photons(t, n_electron, z_int, xy_int, sc_gain, config, resource): :param t: 1d int array time of s2 :param n_electron: 1d float array number of electrons to simulate :param z_int: float array. true Z interaction positions - :param positions: 2d float array, true xy interaction positions + :param xy_int: 2d float array, true xy interaction positions :param sc_gain: float, secondary s2 gain :param config: dict of the wfsim config :param resource: instance of the resource class """ @@ -354,7 +375,7 @@ def luminescence_timings_simple(xy, n_photons, config, resource): return S2._luminescence_timings_simple(len(xy), dG, E0, r, dr, rr, alpha, uE, pressure, n_photons) - + @staticmethod def luminescence_timings_garfield(xy, n_photons, config, resource, confine_position=None): """ @@ -364,7 +385,6 @@ def luminescence_timings_garfield(xy, n_photons, config, resource, confine_posit :param config: dict wfsim config :param resource: instance of wfsim resource :param confine_position: if float, confine extraction region +/- this position around anode wires - returns 2d array with ints for photon timings of input param 'shape' """ assert 's2_luminescence' in resource.__dict__, 's2_luminescence model not found' @@ -386,10 +406,84 @@ def luminescence_timings_garfield(xy, n_photons, config, resource, confine_posit avgt = np.average(resource.s2_luminescence['t']).astype(int) return resource.s2_luminescence['t'][index_row, index_col].astype(np.int64) - avgt + + @staticmethod + @njit + def draw_excitation_times(inv_cdf_list, hist_indices, nph, diff_nearest_gg, d_gas_gap): + + """ + Draws the excitation times from the GARFIELD electroluminescence map + + :param inv_cdf_list: List of inverse CDFs for the excitation time histograms + :param hist_indices: The index of the histogram which refers to the gas gap + :param nph: A 1-d array of the number of photons per electron + :param diff_nearest_gg: The difference between the gas gap from the map + (continuous value) and the nearest (discrete) value + of the gas gap corresponding to the excitation time + histograms + :param d_gas_gap: Spacing between two consecutive gas gap values + + returns time of each photon + """ + + inv_cdf_len = len(inv_cdf_list[0]) + timings = np.zeros(np.sum(nph)) + upper_hist_ind = np.clip(hist_indices+1, 0, len(inv_cdf_list)-1) + + count = 0 + for i, (hist_ind, u_hist_ind, n, dngg) in enumerate(zip(hist_indices, + upper_hist_ind, + nph, + diff_nearest_gg)): + + #There are only 10 values of gas gap separated by 0.1mm, so we interpolate + #between two histograms + + interp_cdf = ((inv_cdf_list[u_hist_ind]-inv_cdf_list[hist_ind])*(dngg/d_gas_gap) + +inv_cdf_list[hist_ind]) + + #Subtract 2 because this way we don't want to sample from this last strange tail + samples = np.random.uniform(0, inv_cdf_len-2, n) + t1 = interp_cdf[np.floor(samples).astype('int')] + t2 = interp_cdf[np.ceil(samples).astype('int')] + T = (t2-t1)*(samples - np.floor(samples))+t1 + if n!=0: + T = T-np.mean(T) + + #subtract mean to get proper drift time and z correlation + timings[count:count+n] = T + count+=n + return timings + + @staticmethod + def luminescence_timings_garfield_gasgap(xy, n_photons, resource): + """ + Luminescence time distribution computation according to garfield scintillation maps + which are ONLY drawn from below the anode, and at different gas gaps + :param xy: 1d array with positions + :param n_photons: 1d array with ints for number of xy positions + :param resource: instance of wfsim resource + returns 2d array with ints for photon timings of input param 'shape' + """ + assert 's2_luminescence_gg' in resource.__dict__, 's2_luminescence_gg model not found' + assert len(n_photons) == len(xy), 'Input number of n_electron should have same length as positions' + + d_gasgap = resource.s2_luminescence_gg['gas_gap'][1]-resource.s2_luminescence_gg['gas_gap'][0] + + cont_gas_gaps = resource.garfield_gas_gap_map(xy) + draw_index = np.digitize(cont_gas_gaps, resource.s2_luminescence_gg['gas_gap'])-1 + diff_nearest_gg = cont_gas_gaps - resource.s2_luminescence_gg['gas_gap'][draw_index] + + return S2.draw_excitation_times(resource.s2_luminescence_gg['timing_inv_cdf'], + draw_index, + n_photons, + diff_nearest_gg, + d_gasgap) + @staticmethod def optical_propagation(channels, config, spline): - """Function gettting times from s2 timing splines: + """Function getting times from s2 timing splines: :param channels: The channels of all s2 photon :param config: current configuration of wfsim :param spline: pointer to s2 optical propagation splines from resources @@ -425,7 +519,6 @@ def photon_timings(positions, n_photons_per_xy, _electron_timings, n_photons_per config=config, resource=resource) elif config['s2_luminescence_model']=='garfield': - # check to see if extraction region in Garfield needs to be confined confine_position=None if 's2_garfield_confine_position' in config: if config['s2_garfield_confine_position'] > 0.0: @@ -434,8 +527,12 @@ def photon_timings(positions, n_photons_per_xy, _electron_timings, n_photons_per config=config, resource=resource, confine_position=confine_position) + + elif config['s2_luminescence_model']=='garfield_gas_gap': + _photon_timings = S2.luminescence_timings_garfield_gasgap(positions, n_photons_per_xy, + resource=resource) else: - raise KeyError(f"{config['s2_luminescence_model']} is not valid! Use 'simple' or 'garfield'") + raise KeyError(f"{config['s2_luminescence_model']} is not valid! Use 'simple' or 'garfield' or 'garfield_gas_gap'") # Emission Delay _photon_timings += Pulse.singlet_triplet_delays(len(_photon_timings), config['singlet_fraction_gas'], config, phase) @@ -497,7 +594,7 @@ def s2_pattern_map_diffuse(n_electron, z, xy, config, resource): # Should we also output this xy position in truth? xy_multi = np.repeat(xy, n_electron, axis=0) + hdiff # One entry xy per electron # Remove points outside tpc, and the pattern will be the average inside tpc - # Should be done naturally with the s2 pattern map, however, there's some bug there so we apply this hard cut + # Should be done naturally with the s2 pattern map, however, there's some bug there, so we apply this hard cut mask = np.sum(xy_multi ** 2, axis=1) <= config['tpc_radius'] ** 2 if isinstance(resource.s2_pattern_map, DummyMap): diff --git a/wfsim/load_resource.py b/wfsim/load_resource.py index 4d5358ac..56aab8f7 100644 --- a/wfsim/load_resource.py +++ b/wfsim/load_resource.py @@ -89,10 +89,13 @@ def config_to_file(config): 's1_pattern_map': 'XENONnT_s1_xyz_patterns_LCE_corrected_qes_MCva43fa9b_wires.pkl', 's1_lce_correction_map': 'XENONnT_s1_xyz_LCE_corrected_qes_MCva43fa9b_wires.json.gz', 's2_pattern_map': 'XENONnT_s2_xy_patterns_LCE_corrected_qes_MCva43fa9b_wires.pkl', - 's2_correction_map': 'XENONnT_s2_xy_correction_corrected_qes_MCva43fa9b_wires.json.gz', + 's2_correction_map': 'XENONnT_s2_xy_map_v4_210503_mlp_3_in_1_iterated.json', + 'se_gain_map': 'XENONnT_se_xy_map_v1_mlp.json', 'photon_ap_cdfs': 'XENONnT_pmt_afterpulse_config_012605.json.gz', 's2_luminescence': 'XENONnT_GARFIELD_B1d5n_C30n_G1n_A6d5p_T1d5n_PMTs1d5n_FSR0d95n.npz', + "s2_luminescence_gg": "garfield_timing_map_gas_gap_sr0.npy", 'gas_gap_map': 'gas_gap_warping_map_January_2021.pkl', + 'garfield_gas_gap_map': 'garfield_gas_gap_map_sr0.json', 'ele_ap_pdfs': 'x1t_se_afterpulse_delaytime.pkl.gz', 'noise_file': 'x1t_noise_170203_0850_00_small.npz', 'fdc_3d': 'XnT_3D_FDC_xyt_dummy_all_zeros_v0.1.json.gz', @@ -180,7 +183,7 @@ def get_file_path(base, fname): # FileNotFoundError, ValueErrors can be raised if we # cannot load the requested config fpath = downloader.download_single(fname) - log.warning(f"Loading {fname} from mongo downloader") + log.warning(f"Loading {fname} from mongo downloader to {fpath}") return fname # Keep the name and let get_resource do its thing except (FileNotFoundError, ValueError, NameError, AttributeError): @@ -233,6 +236,8 @@ def __init__(self, config=None): pmt_mask = np.array(config['gains']) > 0 # Converted from to pe (from cmt by default) self.s1_pattern_map = make_patternmap(files['s1_pattern_map'], fmt='pkl', pmt_mask=pmt_mask) self.s2_pattern_map = make_patternmap(files['s2_pattern_map'], fmt='pkl', pmt_mask=pmt_mask) + self.se_gain_map = make_map(files['se_gain_map']) +# self.s2_correction_map = make_map(files['s2_correction_map'], fmt = 'json') # if there is a (data driven!) map, load it. If not make it from the pattern map if files['s1_lce_correction_map']: @@ -263,11 +268,11 @@ def __init__(self, config=None): # if there is a (data driven!) map, load it. If not make it from the pattern map if files['s2_correction_map']: - self.s2_correction_map = make_map(files['s2_correction_map']) + self.s2_correction_map = make_map(files['s2_correction_map'], fmt = 'json') else: s2cmap = deepcopy(self.s2_pattern_map) # Lower the LCE by removing contribution from dead PMTs - # AT: masking is a bit redundant due to PMT mask application in make_patternmap + # AT: masking is a bit redundant due to PMT mask application in make_patternmap s2cmap.data['map'] = np.sum(s2cmap.data['map'][:][:], axis=2, keepdims=True, where=pmt_mask) # Scale by median value s2cmap.data['map'] = s2cmap.data['map'] / np.median(s2cmap.data['map'][s2cmap.data['map'] > 0]) @@ -276,13 +281,32 @@ def __init__(self, config=None): # Garfield luminescence timing samples # if config.get('s2_luminescence_model', False) == 'garfield': - if 'garfield' in config.get('s2_luminescence_model', ''): - s2_luminescence_map = straxen.get_resource(files['s2_luminescence'], fmt='npy_pickle')['arr_0'] - # Get directly the map for the simulated level - liquid_level_available = np.unique(s2_luminescence_map['ll']) # available levels (cm) - liquid_level = config['gate_to_anode_distance'] - config['elr_gas_gap_length'] # cm - liquid_level = min(liquid_level_available, key=lambda x: abs(x - liquid_level)) - self.s2_luminescence = s2_luminescence_map[s2_luminescence_map['ll'] == liquid_level] + if 'garfield_gas_gap' in config.get('s2_luminescence_model', ''): + #garfield_gas_gap option is using (x,y) -> gas gap (from the map) -> s2 luminescence + #from garfield. This s2_luminescence_gg is indexed only by the gas gap, and + #corresponds to electrons drawn directly below the anode + + s2_luminescence_map = straxen.get_resource(files['s2_luminescence_gg'], fmt='npy') + self.s2_luminescence_gg = s2_luminescence_map + self.garfield_gas_gap_map = make_map(files['garfield_gas_gap_map'], fmt = 'json') + + elif 'garfield' in config.get('s2_luminescence_model', ''): + #This option indexes the luminescence times using the liquid level values + #as well as the position between the full pitch of two gate wires + gf_file_name = files['s2_luminescence'] + if gf_file_name.endswith('npy'): + s2_luminescence_map = straxen.get_resource(gf_file_name, fmt='npy') + self.s2_luminescence = s2_luminescence_map + elif gf_file_name.endswith('npz'): + # Backwards compatibility from before #363 / #370 + s2_luminescence_map = straxen.get_resource(gf_file_name, fmt='npy_pickle')['arr_0'] + # Get directly the map for the simulated level + liquid_level_available = np.unique(s2_luminescence_map['ll']) # available levels (cm) + liquid_level = config['gate_to_anode_distance'] - config['elr_gas_gap_length'] # cm + liquid_level = min(liquid_level_available, key=lambda x: abs(x - liquid_level)) + self.s2_luminescence = s2_luminescence_map[s2_luminescence_map['ll'] == liquid_level] + else: + raise ValueError(f'{gf_file_name} is of unknown format') if config.get('field_distortion_model', "none") == "inverse_fdc": self.fdc_3d = make_map(files['fdc_3d'], fmt='json.gz') @@ -368,7 +392,7 @@ def make_map(map_file, fmt=None, method='WeightedNearestNeighbors'): raise TypeError("Can't handle map_file except a string or a list") @export def make_patternmap(map_file, fmt=None, method='WeightedNearestNeighbors', pmt_mask=None): - """ This is special interpretation of the of previous make_map(), but designed + """ This is special interpretation of the of previous make_map(), but designed for pattern map loading with provided PMT mask. This way simplifies both S1 and S2 cases """ @@ -399,6 +423,7 @@ def make_patternmap(map_file, fmt=None, method='WeightedNearestNeighbors', pmt_m return straxen.InterpolatingMap(map_data, method=method) else: raise TypeError("Can't handle map_file except a string or a list") + @export class DummyMap: """Return constant results diff --git a/wfsim/pax_interface.py b/wfsim/pax_interface.py index 10b07413..ced04808 100644 --- a/wfsim/pax_interface.py +++ b/wfsim/pax_interface.py @@ -93,7 +93,9 @@ class PaxEventSimulator(object): Call compute to start the simulation process. """ - def __init__(self, config={}): + def __init__(self, config=None): + if config is None: + config = {} self.config = default_config self.config.update(get_resource(self.config['fax_config'], fmt='json')) self.config.update(config)