Getting started with WFSim

Hello friend. Welcome to the basic tutorial on how to simulate waveforms with the latest wfsim version in strax. Here we’ll just demonstrate the basic functionality. For more in depth analysis stuff, checkout the straxen tutorials for more detailed thing.

import numpy as np
import strax
import straxen
import wfsim

import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
from multihist import Histdd, Hist1d
from scipy import stats

import json

*** Detector definition message ***
You are currently using the default XENON10 template detector.

straxen.print_versions(modules=('strax', 'straxen', 'cutax', 'wfsim', 'ntauxfiles'))

ntauxfiles is not installed

	module	version	path	git
0	python	3.9.17	/opt/XENONnT/anaconda/envs/XENONnT_2023.07.1/b...	None
1	strax	1.5.2	/opt/XENONnT/anaconda/envs/XENONnT_2023.07.1/l...	None
2	straxen	2.1.1	/opt/XENONnT/anaconda/envs/XENONnT_2023.07.1/l...	None
3	cutax	1.15.1	/home/shenyangshi/cutax/cutax	branch:bdt_ms | 277f174
4	wfsim	1.0.2	/opt/XENONnT/anaconda/envs/XENONnT_2023.07.1/l...	None

Different simulators

Because we’d like to make simulations we have multiple different ways of making data. They are: * RawRecordsFromFaxNT(1T) * RawRecordsFromFaxEpix * RawRecordsFromFaxOptical * RawRecordsFromFaxNveto

The main difference is the way you specify instructions. The NT(1T) simulator depends on either a csv file or a default/user specified function to generate instructions. Epix takes input from G4 processed with epix. Optical takes optical input from G4 for photon channels and timings. And Nveto is a small modification to optical to provide nveto data

Setting everything up

First we need to define the right context. For this we made two contexts to help you out. One for 1T and one for nT. By default in your context RawRecordsFromFaxNT will be registered. If you want another one be sure to register it

For 1T do:

st = straxen.contexts.xenon1t_simulation()

Or nT do:

st = straxen.contexts.xenonnt_simulation(cmt_run_id_sim='00001')

/opt/XENONnT/anaconda/envs/XENONnT_2023.07.1/lib/python3.9/site-packages/straxen/url_config.py:711: UserWarning: From straxen version 2.1.0 onward, URLConfig parameterswill be sorted alphabetically before being passed to the plugins, this will change the lineage hash for non-sorted URLs. To load data processed with non-sorted URLs, you will need to use an older version.
  warnings.warn("From straxen version 2.1.0 onward, URLConfig parameters"

Now we need to define a run id. What you give it doesn’t really matter, since strax will look for data and make new if it doesn’t find anything. And this is what you want. Strax will use the run id to get the electron lifetime and pmt gains from a database, and returns placeholders if the run doesn’t exist.

run_id = '1'

Defining instructions

For the instructions there are multiple different ways to do it. The simulator has this option called “fax_file”. If it has a value (None by default) the simulator will either read it as a csv or root file. If not it will use some predefined functions to make your events. The number of event you’ll simulate based on the product of the config values “nchunk”,”event_rate” and “chunk_size”. Which you can set as follows:

st.set_config(dict(nchunk=1, event_rate=5, chunk_size=10,))

Strax groups data together in chunks based on time (for low level data). nchunk is the number of chunks you want to simulate event_rate is the number of events per seconds, so this effects the amount of spacing between events. Finally chunk_size is the length in seconds of your chunk For the DAQ this is 5 seconds. For simulations you can do whatever you want. It is important to note that Strax will write out data per chunk So when your chunks are small you’ll, among other things, call Strax’ IO functions a lot, giving a substantial overhead. On the other hand, to large chunks will hog all your memory and your kernel might crash. Based on my experimentation setting chunk_size to 500 gives best performance

These are ways you can give instructions ### Random This is the default where simulator will generate some random instructions for you. ### Custom For this you will need to overwrite the instruction generator function ### CSV You can provide a csv file with the instruction (Like the output of nSort) ### Geant4 For this you’ll need to use epix to do some clustering for you.

Random

I guess this is pretty self explanatory. The simulator has this function called “rand_instructions” which will make something up for you.

wfsim.strax_interface.rand_instructions??

Signature: wfsim.strax_interface.rand_instructions(c)
Source:
def rand_instructions(c):
    """Random instruction generator function. This will be called by wfsim if you do not specify 
    specific instructions.
    :params c: wfsim configuration dict"""
    log.warning('rand_instructions is deprecated, please use wfsim.random_instructions')
    if 'drift_field' not in c:
        log.warning('drift field not specified!')
    # Do the conversion for now, don't break everything all at once
    return _rand_instructions(event_rate=c.get('event_rate', 10),
                              chunk_size=c.get('chunk_size', 5),
                              n_chunk=c.get('n_chunk', 2),
                              energy_range=[1, 100],  # keV
                              drift_field=c.get('drift_field', 100),  # V/cm
                              tpc_radius=c.get('tpc_radius', straxen.tpc_r),
                              tpc_length=c.get('tpc_length', straxen.tpc_z),
                              nest_inst_types=[7]
                              )
File:      /opt/XENONnT/anaconda/envs/XENONnT_2023.07.1/lib/python3.9/site-packages/wfsim/strax_interface.py
Type:      function

Custom

This is some more fun. To do this we’ll write a new function which returns a structured numpy array with the correct dtype.

wfsim.instruction_dtype

[(('Waveform simulator event number.', 'event_number'), numpy.int32),
 (('Quanta type (S1 photons or S2 electrons)', 'type'), numpy.int8),
 (('Time of the interaction [ns]', 'time'), numpy.int64),
 (('X position of the cluster [cm]', 'x'), numpy.float32),
 (('Y position of the cluster [cm]', 'y'), numpy.float32),
 (('Z position of the cluster [cm]', 'z'), numpy.float32),
 (('Number of quanta', 'amp'), numpy.int32),
 (('Recoil type of interaction.', 'recoil'), numpy.int8),
 (('Energy deposit of interaction', 'e_dep'), numpy.float32),
 (('Eventid like in geant4 output rootfile', 'g4id'), numpy.int32),
 (('Volume id giving the detector subvolume', 'vol_id'), numpy.int32),
 (('Local field [ V / cm ]', 'local_field'), numpy.float64),
 (('Number of excitons', 'n_excitons'), numpy.int32),
 (('X position of the primary particle [cm]', 'x_pri'), numpy.float32),
 (('Y position of the primary particle [cm]', 'y_pri'), numpy.float32),
 (('Z position of the primary particle [cm]', 'z_pri'), numpy.float32)]

Event number is just a lable which peaks are together. type is either 1(S1) or 2 (S2). In the truth higher numbers are also used to refer to different types of afterpulses. time,x,y,z are the time and positions of the signal. Amp is the number of photons or electrons generated, and recoil can be used for different types of recoil (but only Electronic recoil is supported). Recoil is a number according to the NEST convention which will be used to indicate which interaction. For the key look here: https://github.com/XENONnT/WFSim/blob/2c614b0f7b0d7c7adc516f6188e281857e8d7e22/wfsim/core.py#L22

Now, lets say we want some krypton peaks. For this we’ll need to change the default instruction function to include this double decay and use Nestpy to convert energy deposits into a number of photons and electrons.

In this case I’ll use 1 “event” per full decay, that where all the 4’s are coming from.

def super_awesome_custom_instruction(c):
    import nestpy
    half_life = 156.94e-9 #Kr intermediate state half-life in ns
    decay_energies = [32.2,9.4] # Decay energies in kev

    n = c['nevents'] = c['event_rate'] * c['chunk_size']
    c['total_time'] = c['chunk_size']

    instructions = np.zeros(4 * n, dtype=wfsim.instruction_dtype)
    instructions['event_number'] = np.repeat(np.arange(n), 4)

    instructions['type'] = np.tile([1, 2], 2 * n)
    instructions['recoil'] = [7 for i in range(4 * n)]

    r = np.sqrt(np.random.uniform(0, 2500, n))
    t = np.random.uniform(-np.pi, np.pi, n)
    instructions['x'] = np.repeat(r * np.cos(t), 4)
    instructions['y'] = np.repeat(r * np.sin(t), 4)
    instructions['z'] = np.repeat(np.random.uniform(-100, 0, n), 4)

    #To get the correct times we'll need to include the 156.94 ns half life of the intermediate state.

    uniform_times = c['total_time'] * (np.arange(n) + 0.5) / n
    delayed_times = uniform_times + np.random.exponential(half_life/np.log(2),len(uniform_times))
    instructions['time'] = np.repeat(list(zip(uniform_times,delayed_times)),2) * 1e9


    # Here we'll define our XENON-like detector
    nc = nestpy.NESTcalc(nestpy.VDetector())
    A = 131.293
    Z = 54.
    density = 2.862  # g/cm^3   #SR1 Value
    drift_field = 82  # V/cm    #SR1 Value
    interaction = nestpy.INTERACTION_TYPE(7)

    energy = np.tile(decay_energies,n)
    quanta = []
    exciton = []
    for en in energy:
        y = nc.GetYields(interaction,
                         en,
                         density,
                         drift_field,
                         A,
                         Z,
                         (1, 1))
        q = nc.GetQuanta(y, density)
        quanta.append(q.photons)
        quanta.append(q.electrons)
        exciton.append(q.excitons)
        exciton.append(0)

    instructions['amp'] = quanta
    instructions['local_field'] = drift_field
    instructions['n_excitons'] = exciton

    return instructions

Now here comes the magic line:

wfsim.strax_interface.rand_instructions = super_awesome_custom_instruction

This changes the default rand_instruction function in our own super awesome function. So when the simulator will call rand_instruction the code from super_awesome_custom_instruction will be executed

CSV

The format for csv files is the same as the instructions dtype. So on every line specify event_number,type,time ,x,y,z, amp and recoil in that order. Then tell the simulator it exists:

st.set_config(dict(fax_file='instructions.csv'))

Ofcourse if you do not have this file it will not work

Geant4

For starters you’ll need to register the Epix plugin. This requires an epix configuration (it’s on private_nt_aux_files), and if you want a start and stop event.

st.register(wfsim.strax_interface.RawRecordsFromMcChain)
st.set_config(dict(fax_file    = 'path_to_g4_file',
                   epix_config = 'path_to_config'))

Optical

Similar to above, register correct plugin and file:

st.register(wfsim.strax_interface.RawRecordsFromMcChain)
st.set_config(dict(fax_file = 'path_to_g4_file',))

nVeto

Because the nveto is a different system it has a slightly different configuration. This is not entirely finished. Once it is you can run like so:

st.register(wfsim.strax_interface.RawRecordsFromFaxnVeto)
st.set_config(dict(fax_config_nveto = 'path_to_nveto_config_file',
                  fax_file          = 'path_to_g4_file',))

Configuration customization

The simulator using a larger large amount of configuration settings to do it’s magic. Some of them are best left along, like pmt_circuit_load_resistor. Others on the other hand are things you might want change a bit to see how the data will change. Unfortunately currently the full list is spread out over two different places. One is the fax config json which is on github. The other is the option list in strax. Besides those things like pattern maps are hardcoded in load_resource.py. These files are overrideable with the option “fax_config_override”

The strax config is viewable like this and can be changed by st.set_config(dict(option you want=value you want))

st.show_config('raw_records')

	option	default	current	applies_to	help
0	gain_model_mc	cmt://to_pe_model?version=ONLINE&run_id=plugin...	cmt://to_pe_model?run_id=00001&version=ONLINE	(raw_records, raw_records_he, raw_records_aqmo...	PMT gain model. Specify as (model_type, model_...
1	detector	XENONnT	XENONnT	(raw_records, raw_records_he, raw_records_aqmo...
2	event_rate	1000	5	(raw_records, raw_records_he, raw_records_aqmo...	Average number of events per second
3	chunk_size	100	1	(raw_records, raw_records_he, raw_records_aqmo...	Duration of each chunk in seconds
4	n_chunk	10	<OMITTED>	(raw_records, raw_records_he, raw_records_aqmo...	Number of chunks to simulate
5	per_pmt_truth	False	<OMITTED>	(raw_records, raw_records_he, raw_records_aqmo...	Store the info per channel in the truth file
6	fax_file	None	path_to_g4_file	(raw_records, raw_records_he, raw_records_aqmo...	Directory with fax instructions
7	fax_config	fax_config_nt_design.json	fax_config_nt_design.json	(raw_records, raw_records_he, raw_records_aqmo...
8	fax_config_override	None	<OMITTED>	(raw_records, raw_records_he, raw_records_aqmo...	Dictionary with configuration option overrides
9	fax_config_override_from_cmt	None	{'drift_time_gate': ('cmt_run_id', '00001', 'e...	(raw_records, raw_records_he, raw_records_aqmo...	Dictionary of fax parameter names (key) mapped...
10	channel_map	<OMITTED>	(tpc, he, aqmon, aqmon_nv, tpc_blank, mv, aux_...	(raw_records, raw_records_he, raw_records_aqmo...	immutabledict mapping subdetector to (min, max...
11	n_tpc_pmts	<OMITTED>	494	(raw_records, raw_records_he, raw_records_aqmo...	Number of pmts in tpc. Provided by context
12	n_top_pmts	<OMITTED>	253	(raw_records, raw_records_he, raw_records_aqmo...	Number of pmts in top array. Provided by context
13	right_raw_extension	100000	<OMITTED>	(raw_records, raw_records_he, raw_records_aqmo...
14	seed	False	<OMITTED>	(raw_records, raw_records_he, raw_records_aqmo...	Option for setting the seed of the random numb...
15	gain_model_nv	<OMITTED>	cmt://to_pe_model_nv?run_id=00001&version=ONLINE	(raw_records, raw_records_he, raw_records_aqmo...	nveto gain model, provided by context
16	epix_config	{}	path_to_config	(raw_records, raw_records_he, raw_records_aqmo...	Dict with epix configuration
17	entry_start	0	<OMITTED>	(raw_records, raw_records_he, raw_records_aqmo...
18	entry_stop	None	<OMITTED>	(raw_records, raw_records_he, raw_records_aqmo...	G4 id event number to stop at. If -1 process t...
19	fax_config_nveto	None	path_to_nveto_config_file	(raw_records, raw_records_he, raw_records_aqmo...
20	fax_config_override_nveto	None	<OMITTED>	(raw_records, raw_records_he, raw_records_aqmo...	Dictionary with configuration option overrides
21	targets	(tpc,)	<OMITTED>	(raw_records, raw_records_he, raw_records_aqmo...	tuple with what data to simulate (tpc, nveto o...

With the following block you can visualized which fax config you are using. Remeber that the fax config files are stored in private_nt_aux_files github repository, in the sim_files folder.

st.config['fax_config']

'fax_config_nt_design.json'

The config from github can be loaded as:

straxen.get_resource('fax_config_nt_design.json',fmt='json').keys()

dict_keys(['s1_model_type', 's2_time_model', 's2_luminescence_model', 's2_garfield_confine_position', 'se_gain_from_map', 'ext_eff_from_map', 'enable_gas_gap_warping', 'enable_pmt_afterpulses', 'enable_electron_afterpulses', 'enable_noise', 'field_distortion_on', 'field_distortion_model', 'enable_field_dependencies', 'url_base', 'photon_area_distribution', 's1_pattern_map', 's2_pattern_map', 's1_lce_correction_map', 's2_correction_map', 'se_gain_map', 's1_time_spline', 's2_time_spline', 'field_distortion_comsol_map', 'photon_ap_cdfs', 'noise_file', 's2_luminescence_gg', 's2_luminescence', 'gas_gap_map', 'garfield_gas_gap_map', 'field_dependencies_map', 'temperature', 'pressure', 'lxe_dielectric_constant', 'tpc_length', 'tpc_radius', 'anode_wire_radius', 'anode_field_domination_distance', 'elr_gas_gap_length', 'gate_to_anode_distance', 'drift_field', 'anode_voltage', 'diffusion_constant_longitudinal', 'diffusion_constant_transverse', 'drift_time_gate', 'drift_velocity_liquid', 'singlet_fraction_gas', 'singlet_lifetime_gas', 'singlet_lifetime_liquid', 'triplet_lifetime_gas', 'triplet_lifetime_liquid', 's1_ER_alpha_singlet_fraction', 's1_ER_primary_singlet_fraction', 's1_ER_recombination_fraction', 's1_ER_recombination_time', 's1_ER_secondary_singlet_fraction', 's1_NR_singlet_fraction', 'maximum_recombination_time', 's1_decay_spread', 's1_decay_time', 's1_detection_efficiency', 's2_mean_area_fraction_top', 's2_aft_sigma', 's2_aft_skewness', 's2_secondary_sc_gain', 's2_time_spread', 'electron_extraction_yield', 'electron_lifetime_liquid', 'electron_trapping_time', 'gas_drift_velocity_slope', 'g2_mean', 'p_double_pe_emision', 'pe_pulse_ts', 'pe_pulse_ys', 'pmt_ap_modifier', 'pmt_ap_t_modifier', 'pmt_pulse_time_rounding', 'pmt_transit_time_mean', 'pmt_transit_time_spread', 'photoelectric_modifier', 'photoelectric_p', 'photoelectric_t_center', 'photoelectric_t_spread', 'photoionization_modifier', 'sample_duration', 'samples_after_pulse_center', 'samples_before_pulse_center', 'samples_to_store_after', 'samples_to_store_before', 'pmt_circuit_load_resistor', 'external_amplification', 'high_energy_deamplification_factor', 'trigger_window', 'digitizer_bits', 'digitizer_reference_baseline', 'digitizer_voltage_range', 'special_thresholds', 'zle_threshold'])

Changing things in this guy goes slightly different. In the strax option list there is the option called “fax_config_override”. This takes a dict which will be used to override any values in the json config. So changing the ‘s2_secondary_sc_gain’ is done as:

st.set_config(dict(fax_config_override = dict(s2_secondary_sc_gain=23)))

Things you might want to change

There are 4 main things you might be interested in changing. They are: * Noise (True/False) * PMTAfterpulses (True/False) * Electron Afterpulses (True/False) * S2 luminescence (Simple or garfield)

They are controled by config settings you can change as so:

st.set_config(dict(fax_config_override=dict(
            s2_luminescence_model='simple',
            enable_noise=False,
            enable_pmt_afterpulses=False,
            enable_electron_afterpulses=False,)))

What actually happens?

What happens behind the scenes is that the instructions are first grouped together in chunks. Then we loop over the instructions and the full chunk is returned before starting with the next one.

We use a S1 and S2 class to calculate the arrival times of the photons and the channels which have been hit. Then we’ll hand them over to the Pulse class to calculate the currents in the channels. Finally the currents go to a RawData class where we fake the digitizer response.

S1

For S1s we start with calculating the light yield based on the position of the interaction, and draw the number of photons seen from a Poisson distribution.

Second we calculate the arrival times of the photons. This is based on the scintillation of the xenon atoms. It is dependend on the recombination time, the singlet and triplet fractions.

Finally the channels are calculated. Based on the pattern map we use a interpolation map to get a probability distribution for channels to be hit for a S1 signal based on the position of the interaction, and then we draw from this distribution for every photon.

S2

S2s are slightly more complicated. First we need to drift the electrons up, and while doing so we’ll lose some of them. To get the photon timings, we first need to get the arrival times of the electrons at the gas interface based on a diffusion model. Then we can calculate the photon timings based on a luminescence model for every individual electron. And for the channels we do the same trick with the interpolating map.

Pulse

When we have our lists of channels and timing we can generate actual pulses. First we add a pmt transition time. Then we loop over all channels, calculate the double pe emission probabilities, and add a current in the pmt channel based on the arrival time. This is all stored in a big dictionary. Afterwards this is passed to our fake digitizers which then returns you with your very own pretty data

Getting down to bussiness

Now we have access to all the normal strax data types, and another one called ‘truth’ which holds the simulation instructions. Calling it follows the normal strax convention.

st = straxen.contexts.xenonnt_simulation(cmt_run_id_sim='00001')
st.set_config(dict(event_rate=5, chunk_size=1, fax_file=None))

# Remove any previously simulated data, if such exists
# !rm -r strax_data

records = st.get_array(run_id,'records')
# peaks = st.get_array(run_id, ['peaks','peak_basics'])
# data = st.get_df(run_id, 'event_info')

truth = st.get_df(run_id, 'truth')

Now it is time to make pretty plots and see if what we makes actually makes any sense

peak_basics = st.get_df(run_id, 'peak_basics')

peak_basics.head(10)

	time	endtime	center_time	area	n_hits	n_channels	max_pmt	max_pmt_area	n_saturated_channels	range_50p_area	range_90p_area	area_fraction_top	length	dt	rise_time	tight_coincidence	type	max_diff	min_diff
0	100000030	100000670	100000205	45.194511	14	13	345	5.892519	0	198.728180	327.420563	0.138998	64	10	94.678757	7	2	80	10
1	100622470	100626210	100624258	143.898376	41	17	50	45.115543	0	1164.394287	2533.276855	0.875330	187	20	930.115356	7	2	500	0
2	300000020	300000410	300000128	22.863621	7	6	483	5.560007	0	73.276711	153.896225	0.276379	39	10	53.755302	5	1	40	10
3	300724040	300724470	300724142	7.752567	3	3	461	3.555671	0	92.052193	296.905457	0.000000	43	10	25.530293	2	1	180	10
4	500000050	500000560	500000248	51.385868	11	11	408	20.320269	0	166.012741	310.153564	0.129670	51	10	142.900955	2	2	90	0
5	500589980	500593580	500591982	37.373817	15	10	241	8.749174	0	1940.112183	3275.072021	0.688202	180	20	2096.490967	3	2	-1	-1
6	700000040	700000610	700000139	15.011408	4	4	324	4.820670	0	57.526924	345.084930	0.000000	57	10	49.048229	3	1	270	20
7	700699840	700706240	700703462	102.122833	25	14	129	33.482544	0	3729.012451	5596.686035	0.880986	160	40	2050.706543	4	2	-1	-1
8	900000050	900000720	900000210	28.657455	8	8	346	7.821901	0	151.675720	443.599609	0.229542	67	10	58.658257	5	1	290	0
9	900703540	900707830	900705461	99.501129	24	16	64	27.903099	0	1144.565796	3354.698730	0.942021	143	30	838.284546	2	2	610	10

Matching

To do matching with the truth the easiest way is to write a new strax plugin where you loop over peaks and get the truth arrays where the mean arrival time of the photons are within the time window of the peak So that will look something like this:

class MatchedPeaks(strax.LoopPlugin):
    depends_on = ('peak_basics', 'truth')
    provides = 'matched_peaks'
    __version__ = '0.0.2'

    dtype = [('time', np.int),
             ('endtime', np.int),
             ('area', np.int),
             ('n_photon', np.int)]

    def compute(self, peaks, truth):
        result = np.zeros(len(truth), self.dtype)

        temp = strax.touching_windows(peaks, truth)
        idx_truth = []
        idx_peak = []
        for ix, (l, r) in enumerate(temp):
            res = result[ix]
            res['time'] = truth['time'][ix]
            res['endtime'] = strax.endtime(truth)[ix]
            res['n_photon'] = truth['n_photon'][ix]

            if r > l:
                i_peak = l+np.argmax(peaks['area'][l:r])
                res['area'] = peaks['area'][i_peak]

        return result

/tmp/jobs/29284557/ipykernel_277/3800134942.py:6: DeprecationWarning: `np.int` is a deprecated alias for the builtin `int`. To silence this warning, use `int` by itself. Doing this will not modify any behavior and is safe. When replacing `np.int`, you may wish to use e.g. `np.int64` or `np.int32` to specify the precision. If you wish to review your current use, check the release note link for additional information.
Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations
  dtype = [('time', np.int),
/tmp/jobs/29284557/ipykernel_277/3800134942.py:7: DeprecationWarning: `np.int` is a deprecated alias for the builtin `int`. To silence this warning, use `int` by itself. Doing this will not modify any behavior and is safe. When replacing `np.int`, you may wish to use e.g. `np.int64` or `np.int32` to specify the precision. If you wish to review your current use, check the release note link for additional information.
Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations
  ('endtime', np.int),
/tmp/jobs/29284557/ipykernel_277/3800134942.py:8: DeprecationWarning: `np.int` is a deprecated alias for the builtin `int`. To silence this warning, use `int` by itself. Doing this will not modify any behavior and is safe. When replacing `np.int`, you may wish to use e.g. `np.int64` or `np.int32` to specify the precision. If you wish to review your current use, check the release note link for additional information.
Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations
  ('area', np.int),
/tmp/jobs/29284557/ipykernel_277/3800134942.py:9: DeprecationWarning: `np.int` is a deprecated alias for the builtin `int`. To silence this warning, use `int` by itself. Doing this will not modify any behavior and is safe. When replacing `np.int`, you may wish to use e.g. `np.int64` or `np.int32` to specify the precision. If you wish to review your current use, check the release note link for additional information.
Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations
  ('n_photon', np.int)]

Of course this doesn’t actually work. An electron afterpulse can be very spread out leading it to be interpreted as multiple peaks while coming from 1 instruction falling outside of the specified range. It would be very much appreciated if someone wants make a more sturdy selection criteria :)

For example, checkout out peak matching algorithm for WFSim: https://github.com/XENONnT/pema

st.register(MatchedPeaks)

__main__.MatchedPeaks

st.get_df(run_id, 'matched_peaks').head(10)

	time	endtime	area	n_photon
0	100000029	100000630	45	251
1	100000333	100000791	45	69
2	100620431	100627474	143	856
3	100622911	100627768	143	211
4	300000006	300000565	22	255
5	300000067	300000556	22	69
6	300723015	300728784	7	186
7	300723660	300727606	7	128
8	500000003	500000619	51	235
9	500000111	500000642	51	49

Externally call functions

The s1/s2 functions used for generating photon timings and channels are now doocumented and callable from the outside. This you can use if you want to do something like compare photon arrival timings between wfsims results and g4

st = straxen.contexts.xenonnt_simulation(cmt_run_id_sim='00001')

Object `wfsim.core.S1.get_n_photons` not found.

config = straxen.get_resource('fax_config_nt_sr0_v4.json', fmt='json')
config.update({'detector':'XENONnT', 'right_raw_extension':50000})
st.config['gains'] = [1 for _ in range(straxen.n_tpc_pmts)]
config.update(st.config)

s2 = wfsim.core.s2.S2(config)

WARNING:wfsim.resource:Loading XENONnT_SR0_spe_distributions_20210713_no_noise_scaled.csv from mongo downloader to /dali/lgrandi/strax/resource_cache/99c2cbc580cfa8eeebe831456076e136
WARNING:wfsim.resource:Loading XENONnT_s1_xyz_patterns_LCE_MCvf051911_wires.pkl from mongo downloader to /dali/lgrandi/strax/resource_cache/02ed2634596793f6e7aa13783745089a
WARNING:wfsim.resource:A file has value False, assuming this is intentional.
WARNING:wfsim.resource:Loading XENONnT_s2_xy_patterns_GXe_LCE_corrected_qes_MCv4.3.0_wires.pkl from mongo downloader to /dali/lgrandi/strax/resource_cache/6639b01b73d4fdbb703174ad424801f5
WARNING:wfsim.resource:Loading XENONnT_s2_xy_map_v4_210503_mlp_3_in_1_iterated.json from mongo downloader to /dali/lgrandi/strax/resource_cache/ca348171778ae5cdb5de6c1f5814f156
WARNING:wfsim.resource:Loading XENONnT_se_xy_map_v1_mlp.json from mongo downloader to /dali/lgrandi/strax/resource_cache/4ed37f1d06853503899ac37a0a2237ba
WARNING:wfsim.resource:Loading XENONnT_pmt_afterpulse_config_018435.json.gz from mongo downloader to /dali/lgrandi/strax/resource_cache/a38b18cd61ca67a065a3a32d9b57c2b2
WARNING:wfsim.resource:Loading XENONnT_GARFIELD_SR0_B2d75n_C2d75n_G0d3p_A4d9p_T0d9n_PMTs1d3n_FSR0d65p.npz from mongo downloader to /dali/lgrandi/strax/resource_cache/42d87bc750ceb2f3c0491a586e016341
WARNING:wfsim.resource:Loading garfield_timing_map_gas_gap_sr0.npy from mongo downloader to /dali/lgrandi/strax/resource_cache/803c5f08cb89063a9671ac8da5cb7e75
WARNING:wfsim.resource:Loading gas_gap_warping_map_January_2021.pkl from mongo downloader to /dali/lgrandi/strax/resource_cache/d37b34dda3e865ef607b0bee2057772c
WARNING:wfsim.resource:Loading garfield_gas_gap_map_sr0.json from mongo downloader to /dali/lgrandi/strax/resource_cache/d39d7eda03373635c60e209067470dc8
WARNING:wfsim.resource:Loading x1t_se_afterpulse_delaytime.pkl.gz from mongo downloader to /dali/lgrandi/strax/resource_cache/014aa048758a51f154cdc3b32dac2773
WARNING:wfsim.resource:Loading XENONnT_noise_tpc_only_2ms_25118.npz from mongo downloader to /dali/lgrandi/strax/resource_cache/e560759e28f74d5c049425c25e5ae0d5
WARNING:wfsim.resource:Loading XnT_3D_FDC_xyt_dummy_all_zeros_v0.1.json.gz from mongo downloader to /dali/lgrandi/strax/resource_cache/55b30b2d26d151ae8e2aa55ffe6693ad
WARNING:wfsim.resource:Loading field_dependent_radius_depth_maps_B2d75n_C2d75n_G0d3p_A4d9p_T0d9n_PMTs1d3n_FSR0d65p_QPTFE_0d5n_0d4p.json.gz from mongo downloader to /dali/lgrandi/strax/resource_cache/487272e33b4c832ae715350b6275ba0b
WARNING:wfsim.resource:Loading data_driven_diffusion_map_XENONnTSR0V2.json.gz from mongo downloader to /dali/lgrandi/strax/resource_cache/8a6c6dc1bfe5380f48e96462161a3313
WARNING:wfsim.resource:Loading XENONnT_s1_proponly_pc_reflection_optPhot_perPMT_S1_local_20220510.json.gz from mongo downloader to /dali/lgrandi/strax/resource_cache/91f4a6162d5e335c4416165349222061
WARNING:wfsim.resource:Loading XENONnT_s2_opticalprop_time_v0.json.gz from mongo downloader to /dali/lgrandi/strax/resource_cache/1d257e215e00873e37ae49f16a04d717
WARNING:wfsim.resource:Using json for unspecified XENONnT_se_xy_map_v1_mlp.json
/opt/XENONnT/anaconda/envs/XENONnT_2023.07.1/lib/python3.9/site-packages/wfsim/load_resource.py:262: RuntimeWarning: invalid value encountered in true_divide
  orig_aft_=np.mean((s2map_topeff_/s2map_toteff_)[s2map_toteff_>0.0])
WARNING:wfsim.resource:Using json.gz for unspecified XENONnT_s2_opticalprop_time_v0.json.gz
WARNING:wfsim.resource:Using noise data XENONnT_noise_tpc_only_2ms_25118.npz with 494 channels for XENONnT

s2.electron_timings?

Signature:
s2.electron_timings(
    t,
    n_electron,
    drift_time_mean,
    drift_time_spread,
    sc_gain,
    timings,
    gains,
    electron_trapping_time,
)
Call signature:  s2.electron_timings(*args, **kwargs)
Type:            CPUDispatcher
String form:     CPUDispatcher(<function S2.electron_timings at 0x7f8d3dae25e0>)
File:            /opt/XENONnT/anaconda/envs/XENONnT_2023.07.1/lib/python3.9/site-packages/wfsim/core/s2.py
Docstring:
Calculate arrival times of the electrons. Data is written to the timings and gains arrays
:param t: 1d array of ints
:param n_electron:1 d array of ints
:param drift_time_mean: 1d array of floats
:param drift_time_spread: 1d array of floats
:param sc_gain: secondary scintillation gain
:param timings: empty array with length sum(n_electron)
:param gains: empty array with length sum(n_electron)
:param electron_trapping_time: configuration values
Class docstring:
Implementation of user-facing dispatcher objects (i.e. created using
the @jit decorator).
This is an abstract base class. Subclasses should define the targetdescr
class attribute.
Init docstring:
Parameters
----------
py_func: function object to be compiled
locals: dict, optional
    Mapping of local variable names to Numba types.  Used to override
    the types deduced by the type inference engine.
targetoptions: dict, optional
    Target-specific config options.
impl_kind: str
    Select the compiler mode for `@jit` and `@generated_jit`
pipeline_class: type numba.compiler.CompilerBase
    The compiler pipeline type.