xenonnt / epix Goto Github PK

Electron and Photon Instructions generator for XENON

License: BSD 3-Clause "New" or "Revised" License

Python 100.00%

epix's Issues

Long Events

In some simulations delayed interactions (seconds, hours, days) are present. These interactions are cerated by e.g. neutron activations or long living isomers.

At the moment these interactions are cut to prevent possible issues with the waveform simulator. In GenerateGeant4.py these delayed interactions are postponed in new events.

We can think about adding this functionality to EPIX.

Keep track of removed interactions/events

We remove several interactions and events from the initial Geant4 simulation along the way in epix. One example would be the removal of interactions outside the specified detector volume. We could add the option to save all removed interactions/events in a separate file to help with debugging.

NRs beyond NEST validity

We currently get a lot of these while processing our files:

.../epix/quanta_generation.py:51: UserWarning: Energy deposition of 330.34259855747223 keV beyond
NEST validity for NR model of 200 keV - Remove Interaction
   warnings.warn(f"Energy deposition of {en} keV beyond NEST validity for
   NR model of 200 keV - Remove Interaction")

which could indicate that our NR clusters are probably too big, or our NR selection is unprecise. As expected, these messages pop up more often when we enlarge the MicroSeparation size, but the current default is already quite small (50 microns). I would like to revisit these two key points and hopefully figure out the cluster size that NEST "expects".

To reproduce the warning, you can just run:

bin/run_epix --InputFile /dali/lgrandi/xenonnt/simulations/testing_epix_wfsim/tpc_cryoneutrons_200.root\
--Detector XENONnT --OutputPath . --Debug

with the current epix release.

Space and time clustering

In EPIX we want to cluster interactions that are close to each other in space but prevent interactions to be merged when they are separated in time. At the moment we are using a 4D DBSCAN clustering where the time is scaled to a length.

We seem to agree that first a 1D clustering in time followed by a 3D clustering in time makes more sense than the 4D clustering. The 1D clustering can be performed with a simple algorithm and DBSCAN can be used for the 3D clustering.

Cluster Classification

The charge and light yield and thus the number of produced electrons and photons differs for interactions with the same energy deposition but different involved particles and interaction processes. As a result, NEST contains several models to calculate the yields for different particles.

In EPIX we have to specify which NEST model we want to use for each cluster. At the moment the clusters are classified based on the option energy or time.

Energy

The interaction with the highes energy deposition in a cluster determines the NEST model for the whole cluster.
This approach was used in GenerateGeant4.py.
Problem: In most cases electrons (from e.g. ionization processes) contribute in a large fraction to the total energy deposition. Since these individual energy depositions can not be resolved in an experiment, the NEST beta model is meant for electrons from e.g. beta decays so that we will (most likely) get a wrong estimate of the number of photons and electrons.

Time

The first interaction determines the NEST model for the whole cluster. In this case we assume that all following energy depositions (in this cluster) are somehow causally connected to this first interaction.
Problem: Interactions from independent processes can also be merged e.g. interactions from a beta decay with subsequent gamma emission. In this case we should use the gamma and beta model of NEST.

As both options have their advantages and drawbacks we need to discuss which option will give us better simulations so it can be used as default case. We could also think about other options like the reconstruction of causal relations of the individual energy depositions.

Add excitons to instructions

If we want to use scintillation delay times from NEST, we need to store number of excitons. It not important if we use epix with current WFsim hit time generation,so can be ignored in the code. However, it can be important if we want to model more realistic S1 times for G4 photon propagator or pulse shape discrimination.

`culstering._write_result` fails in numba `v0.56.0`

In numba==0.56.0, this line fails.

See e.g. https://github.com/XENONnT/WFSim/runs/7605518604?check_suite_focus=true.

Can we remove the numba decorator or fix it to work in the latest numba? I'm not very familiar with akward arrays and numba myself unfortunately

To reproduce:

source /cvmfs/xenon.opensciencegrid.org/releases/nT/2022.06.3/setup.sh
pip install numba==0.56.0 --user
cd software/wfsim/ # in case you don't have some wfsim version somewhere, do git clone https://github.com/XENONnT/WFSim first
pytest tests/test_wfsim.py -k test_sim_mc_chain

Traceback

(XENONnT_2022.06.3) [angevaare@dali-login1 wfsim]$ pytest tests/test_wfsim.py -k test_sim_mc_chain
============================================================================================================ test session starts ============================================================================================================
platform linux -- Python 3.8.13, pytest-7.1.2, pluggy-1.0.0
rootdir: /home/angevaare/software/wfsim, configfile: pytest.ini
plugins: anyio-3.6.1, hypothesis-6.47.0
collected 9 items / 8 deselected / 1 selected

tests/test_wfsim.py F                                                                                                                                                                                                                 [100%]

================================================================================================================= FAILURES ==================================================================================================================
_____________________________________________________________________________________________________________ test_sim_mc_chain _____________________________________________________________________________________________________________

    @skipIf(not straxen.utilix_is_configured(), 'utilix is not configured')
    def test_sim_mc_chain():
        """Test the nT simulator by Geant4 output file"""

        with tempfile.TemporaryDirectory() as tempdir:
            log.debug(f'Working in {tempdir}')

            # Download test file on the test directory
            import requests
            test_g4 = 'https://raw.githubusercontent.com/XENONnT/WFSim/master/tests/geant_test_data_small.root'
            url_data = requests.get(test_g4).content
            with open('test.root', mode='wb') as f:
                f.write(url_data)
            st = straxen.contexts.xenonnt_simulation(cmt_run_id_sim='010000',
                                                     cmt_version='global_ONLINE',
                                                     _config_overlap={},)
            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.,
                           'max_delay': 1e7, 'detector_config_override': None, 'micro_separation': 0.05,
                           'tag_cluster_by': 'time', 'nr_only': False}

            st.register(wfsim.RawRecordsFromMcChain)
            st.set_config(dict(
                detector='XENONnT',
                fax_file='./test.root',
                event_rate=100.,
                chunk_size=5.,
                entry_start=0,
                fax_config='fax_config_nt_design.json',
                fax_config_override=dict(
                    s1_model_type='nest',
                    s2_time_model="s2_time_spread around zero",
                    enable_electron_afterpulses=False),
                epix_config=epix_config,
                fax_config_nveto='fax_config_nt_nveto.json',
                fax_config_override_nveto=dict(enable_noise=False,
                                               enable_pmt_afterpulses=False,
                                               enable_electron_afterpulses=False,),
                targets=('tpc', 'nveto'),
                baseline_samples_nv=("nv_baseline_constant", 26, True),
            ))

            log.debug(f'Getting raw-records')
>           rr = st.get_array(run_id, 'raw_records')

tests/test_wfsim.py:236:
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
/cvmfs/xenon.opensciencegrid.org/releases/nT/2022.06.3/anaconda/envs/XENONnT_2022.06.3/lib/python3.8/site-packages/strax/context.py:1440: in get_array
    results = [x.data for x in source]
/cvmfs/xenon.opensciencegrid.org/releases/nT/2022.06.3/anaconda/envs/XENONnT_2022.06.3/lib/python3.8/site-packages/strax/context.py:1440: in <listcomp>
    results = [x.data for x in source]
/cvmfs/xenon.opensciencegrid.org/releases/nT/2022.06.3/anaconda/envs/XENONnT_2022.06.3/lib/python3.8/site-packages/strax/context.py:1263: in get_iter
    components = self.get_components(run_id,
/cvmfs/xenon.opensciencegrid.org/releases/nT/2022.06.3/anaconda/envs/XENONnT_2022.06.3/lib/python3.8/site-packages/strax/context.py:1024: in get_components
    d.setup()
wfsim/strax_interface.py:556: in setup
    self.get_instructions()
wfsim/strax_interface.py:792: in get_instructions
    self.instructions_epix = epix.run_epix.main(
/cvmfs/xenon.opensciencegrid.org/releases/nT/2022.06.3/anaconda/envs/XENONnT_2022.06.3/lib/python3.8/site-packages/epix/run_epix.py:45: in main
    result = epix.cluster(inter, args['tag_cluster_by'] == 'energy')
/cvmfs/xenon.opensciencegrid.org/releases/nT/2022.06.3/anaconda/envs/XENONnT_2022.06.3/lib/python3.8/site-packages/epix/clustering.py:164: in cluster
    _cluster(x, y, z, ed, time, ci,
../../.local/lib/python3.8/site-packages/numba/core/dispatcher.py:468: in _compile_for_args
    error_rewrite(e, 'typing')
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

e = TypingError('Failed in nopython mode pipeline (step: nopython frontend)\nInternal error at <numba.core.typeinfer.CallC...simple but extensive with awkward...\n                _write_result(res, x_mean, y_mean, z_mean,\n                ^\n')
issue_type = 'typing'

    def error_rewrite(e, issue_type):
        """
        Rewrite and raise Exception `e` with help supplied based on the
        specified issue_type.
        """
        if config.SHOW_HELP:
            help_msg = errors.error_extras[issue_type]
            e.patch_message('\n'.join((str(e).rstrip(), help_msg)))
        if config.FULL_TRACEBACKS:
            raise e
        else:
>           raise e.with_traceback(None)
E           numba.core.errors.TypingError: Failed in nopython mode pipeline (step: nopython frontend)
E           Internal error at <numba.core.typeinfer.CallConstraint object at 0x7f7d4e54b1f0>.
E           Failed in nopython mode pipeline (step: native lowering)
E           module 'llvmlite' has no attribute 'llvmpy'
E
E           File "../../../../cvmfs/xenon.opensciencegrid.org/releases/nT/2022.06.3/anaconda/envs/XENONnT_2022.06.3/lib/python3.8/site-packages/epix/clustering.py", line 302:
E           def _write_result(res, x_mean, y_mean, z_mean,
E               <source elided>
E               res.begin_record()
E               res.field('x')
E               ^
E
E           During: lowering "$16call_method.6 = call $12load_method.4($const14.5, func=$12load_method.4, args=[Var($const14.5, clustering.py:302)], kws=(), vararg=None, varkwarg=None, target=None)" at /cvmfs/xenon.opensciencegrid.org/releases/nT/2022.06.3/anaconda/envs/XENONnT_2022.06.3/lib/python3.8/site-packages/epix/clustering.py (302)
E           During: resolving callee type: type(CPUDispatcher(<function _write_result at 0x7f7d52a2d310>))
E           During: typing of call at /cvmfs/xenon.opensciencegrid.org/releases/nT/2022.06.3/anaconda/envs/XENONnT_2022.06.3/lib/python3.8/site-packages/epix/clustering.py (211)
E
E           Enable logging at debug level for details.
E
E           File "../../../../cvmfs/xenon.opensciencegrid.org/releases/nT/2022.06.3/anaconda/envs/XENONnT_2022.06.3/lib/python3.8/site-packages/epix/clustering.py", line 211:
E           def _cluster(x, y, z, ed, time, ci,
E               <source elided>
E                           # Write result, simple but extensive with awkward...
E                           _write_result(res, x_mean, y_mean, z_mean,
E                           ^

../../.local/lib/python3.8/site-packages/numba/core/dispatcher.py:409: TypingError
----------------------------------------------------------------------------------------------------------- Captured stdout call ------------------------------------------------------------------------------------------------------------
Removing old incomplete data in ./strax_data/010000-raw_records_nv-4b5zg2gxrb
Removing old incomplete data in ./strax_data/010000-raw_records_he-4b5zg2gxrb
Removing old incomplete data in ./strax_data/010000-raw_records-4b5zg2gxrb
Removing old incomplete data in ./strax_data/010000-raw_records_aqmon-4b5zg2gxrb
Removing old incomplete data in ./strax_data/010000-truth-4b5zg2gxrb
Removing old incomplete data in ./strax_data/010000-truth_nv-4b5zg2gxrb
Using nestpy version 1.5.4
epix configuration:  {'cut_by_eventid': True, 'debug': True, 'source_rate': 0, 'micro_separation_time': 10.0, 'max_delay': 10000000.0, 'detector_config_override': None, 'micro_separation': 0.05, 'tag_cluster_by': 'time', 'nr_only': False, 'detector': 'XENONnT', 'entry_start': 0, 'entry_stop': None, 'input_file': './test.root', 'path': '.', 'file_name': 'test.root', 'detector_config': [<epix.detector_volumes.SensitiveVolume object at 0x7f7d52a34850>, <epix.detector_volumes.SensitiveVolume object at 0x7f7d52a34550>], 'outer_cylinder': {'max_z': 7.3936, 'min_z': -154.6555, 'max_r': 66.4}}
Total entries in input file = 10
Starting to read from g4 eventid 0
It took 0.1848 sec to load data.
Finding clusters of interactions with a dr = 0.05 mm and dt = 10.0 ns
It took 3.3882 sec to find clusters.
------------------------------------------------------------------------------------------------------------- Captured log call -------------------------------------------------------------------------------------------------------------
WARNING  strax:context.py:645 Option gain_model_nv not taken by any registered plugin
WARNING  strax:context.py:645 Option gain_model_mv not taken by any registered plugin
WARNING  strax:context.py:645 Option gain_model_nv not taken by any registered plugin
WARNING  strax:context.py:645 Option gain_model_mv not taken by any registered plugin
WARNING  strax:context.py:645 Option gain_model_nv not taken by any registered plugin
WARNING  strax:context.py:645 Option gain_model_mv not taken by any registered plugin
WARNING  strax:context.py:645 Option gain_model_nv not taken by any registered plugin
WARNING  strax:context.py:645 Option gain_model_mv not taken by any registered plugin
WARNING  strax:context.py:645 Option gain_model_nv not taken by any registered plugin
WARNING  strax:context.py:645 Option gain_model_mv not taken by any registered plugin
WARNING  strax:context.py:645 Option gain_model_nv not taken by any registered plugin
WARNING  strax:context.py:645 Option gain_model_mv not taken by any registered plugin
WARNING  strax:context.py:645 Option gain_model_mv not taken by any registered plugin
WARNING  strax:context.py:645 Option baseline_samples_nv not taken by any registered plugin
============================================================================================================= warnings summary ==============================================================================================================
../../../../cvmfs/xenon.opensciencegrid.org/releases/nT/2022.06.3/anaconda/envs/XENONnT_2022.06.3/lib/python3.8/site-packages/rucio_clients-1.23.14-py3.8.egg/rucio/client/baseclient.py:45
  /cvmfs/xenon.opensciencegrid.org/releases/nT/2022.06.3/anaconda/envs/XENONnT_2022.06.3/lib/python3.8/site-packages/rucio_clients-1.23.14-py3.8.egg/rucio/client/baseclient.py:45: DeprecationWarning: the imp module is deprecated in favour of importlib; see the module's documentation for alternative uses

../../../../cvmfs/xenon.opensciencegrid.org/releases/nT/2022.06.3/anaconda/envs/XENONnT_2022.06.3/lib/python3.8/site-packages/rucio_clients-1.23.14-py3.8.egg/rucio/common/utils.py:668
../../../../cvmfs/xenon.opensciencegrid.org/releases/nT/2022.06.3/anaconda/envs/XENONnT_2022.06.3/lib/python3.8/site-packages/rucio_clients-1.23.14-py3.8.egg/rucio/common/utils.py:668
  /cvmfs/xenon.opensciencegrid.org/releases/nT/2022.06.3/anaconda/envs/XENONnT_2022.06.3/lib/python3.8/site-packages/rucio_clients-1.23.14-py3.8.egg/rucio/common/utils.py:668: DeprecationWarning: invalid escape sequence \?

../../../../cvmfs/xenon.opensciencegrid.org/releases/nT/2022.06.3/anaconda/envs/XENONnT_2022.06.3/lib/python3.8/site-packages/rucio_clients-1.23.14-py3.8.egg/rucio/common/utils.py:669
../../../../cvmfs/xenon.opensciencegrid.org/releases/nT/2022.06.3/anaconda/envs/XENONnT_2022.06.3/lib/python3.8/site-packages/rucio_clients-1.23.14-py3.8.egg/rucio/common/utils.py:669
  /cvmfs/xenon.opensciencegrid.org/releases/nT/2022.06.3/anaconda/envs/XENONnT_2022.06.3/lib/python3.8/site-packages/rucio_clients-1.23.14-py3.8.egg/rucio/common/utils.py:669: DeprecationWarning: invalid escape sequence \?

../../../../cvmfs/xenon.opensciencegrid.org/releases/nT/2022.06.3/anaconda/envs/XENONnT_2022.06.3/lib/python3.8/site-packages/rucio_clients-1.23.14-py3.8.egg/rucio/common/utils.py:670
../../../../cvmfs/xenon.opensciencegrid.org/releases/nT/2022.06.3/anaconda/envs/XENONnT_2022.06.3/lib/python3.8/site-packages/rucio_clients-1.23.14-py3.8.egg/rucio/common/utils.py:670
  /cvmfs/xenon.opensciencegrid.org/releases/nT/2022.06.3/anaconda/envs/XENONnT_2022.06.3/lib/python3.8/site-packages/rucio_clients-1.23.14-py3.8.egg/rucio/common/utils.py:670: DeprecationWarning: invalid escape sequence \?

../../../../cvmfs/xenon.opensciencegrid.org/releases/nT/2022.06.3/anaconda/envs/XENONnT_2022.06.3/lib/python3.8/site-packages/rucio_clients-1.23.14-py3.8.egg/rucio/client/downloadclient.py:1070
../../../../cvmfs/xenon.opensciencegrid.org/releases/nT/2022.06.3/anaconda/envs/XENONnT_2022.06.3/lib/python3.8/site-packages/rucio_clients-1.23.14-py3.8.egg/rucio/client/downloadclient.py:1070
../../../../cvmfs/xenon.opensciencegrid.org/releases/nT/2022.06.3/anaconda/envs/XENONnT_2022.06.3/lib/python3.8/site-packages/rucio_clients-1.23.14-py3.8.egg/rucio/client/downloadclient.py:1070
../../../../cvmfs/xenon.opensciencegrid.org/releases/nT/2022.06.3/anaconda/envs/XENONnT_2022.06.3/lib/python3.8/site-packages/rucio_clients-1.23.14-py3.8.egg/rucio/client/downloadclient.py:1070
  /cvmfs/xenon.opensciencegrid.org/releases/nT/2022.06.3/anaconda/envs/XENONnT_2022.06.3/lib/python3.8/site-packages/rucio_clients-1.23.14-py3.8.egg/rucio/client/downloadclient.py:1070: DeprecationWarning: invalid escape sequence \i

../../../../cvmfs/xenon.opensciencegrid.org/releases/nT/2022.06.3/anaconda/envs/XENONnT_2022.06.3/lib/python3.8/site-packages/panel/util.py:33
  /cvmfs/xenon.opensciencegrid.org/releases/nT/2022.06.3/anaconda/envs/XENONnT_2022.06.3/lib/python3.8/site-packages/panel/util.py:33: DeprecationWarning: distutils Version classes are deprecated. Use packaging.version instead.
    bokeh_version = LooseVersion(bokeh.__version__)

tests/test_wfsim.py::test_sim_mc_chain
  /cvmfs/xenon.opensciencegrid.org/releases/nT/2022.06.3/anaconda/envs/XENONnT_2022.06.3/lib/python3.8/site-packages/epix/clustering.py:40: 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
    df['cluster_id'] = np.zeros(len(df.index), dtype=np.int)

tests/test_wfsim.py::test_sim_mc_chain
  /cvmfs/xenon.opensciencegrid.org/releases/nT/2022.06.3/anaconda/envs/XENONnT_2022.06.3/lib/python3.8/site-packages/epix/common.py:20: 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
    if (array.dtype == np.int) or (array.dtype == np.float64) or (array.dtype == np.float32):

-- Docs: https://docs.pytest.org/en/stable/how-to/capture-warnings.html
========================================================================================================== short test summary info ==========================================================================================================
FAILED tests/test_wfsim.py::test_sim_mc_chain - numba.core.errors.TypingError: Failed in nopython mode pipeline (step: nopython frontend)
=============================================================================================== 1 failed, 8 deselected, 14 warnings in 8.36s ================================================================================================

Unify quanta generation classes

Following #57. In words of @HenningSE: keep the code a bit cleaner if we set up a single quanta generation class containing both the NEST and BBF option. In run_epix.py we would then just call this one class.

Epix is untested

There are no unit-tests for epix. Certainly worth adding some as it saves a lot of time in the long run.

Calculation of average time per cluster can lead to numerical errors

After the DBSCAN clustering is applied, average values of coordinates and time are computed for each cluster. If the time values are large (in radioactive decay events), the average time calculation can lead to numerical errors.

The example of this issue is shown here:

/dali/lgrandi/pkavrigin/2022-12-29_epix_issue_kr85/2022-12-29_epix_issue_kr85.ipynb

Alpha events are always misclassified as NR

The classifier function defined in here always gives nestid=0 for alpha events instead of 6.

This is because, with how the classifier is defined, the type with None, None, None, Ion gets selected.

As a quick fix, just setting

classifier = classifier[::-1]
solves the problem, as the classifier with three None is checked last.

(All the alpha events were then removed by epix because of the too high energy for a NR event)

Incorrect input parameters for NEST

In 'clustering.py' some incorrect parameters for quanta generation with NEST are set, namely:

https://github.com/XENONnT/epix/blob/master/epix/clustering.py#L254-L255

A and Z should be the atomic mass and the atomic number of Xenon:

https://github.com/NESTCollaboration/nest/blob/cfecfdb8dc474eeff08163a361cfad58173ba93e/src/NEST.cpp#L746-L747

So the lines with these parameters should be changed to:

Xe_A = 131.293
Xe_Z = 54.0
classifier['A'] = [Xe_A]*7
classifier['Z'] = [Xe_Z]*7

As an illustration, we can use NEST to produce quanta using the NR spectrum which was generated from one of the AmBe calibration simulations (only elastic NRs were selected using a lineage clustering). Let's run it with A=infinity, Z=0, and then with A=131.3, Z=54.0, using the NR spectrum as an input, and apply g1=0.17, g2=14.1 (estimations from one of Henning's notes from June). Then we can compare it to the actual calibration measurement with base cuts applied to the data. This is the result:

EPIX doesn't make events except initial one after v0.1.2

Issue
EPIX returns only initial event, does not make instructions after that.
What I did to get it
This issue have observed in WFSim with Geant4 file simulation. In v0.1.1, this issue hasn't been observed.
What I investigated to find the cause
With debugger, I found the cause comes from

epix/epix/run_epix.py

Line 115 in 935ba6f

photons, electrons, excitons = epix.quanta_from_NEST(epix.awkward_to_flat_numpy(result['ed']),

After that, a several ten-s of photons, electrons, excitons has actual values, and the rest is 0.

The miminum example is here. (epix version v0.1.2)

$ ./run_epix --InputFile /dali/lgrandi/xenonnt/simulations/testing_epix_wfsim/tpc_and_nveto_cryoneutrons_200.root  --Debug --SourceRate 100
*** Detector definition message ***
You are currently using the default XENON10 template detector.

Using nestpy version 1.4.8
epix configuration:  {'input_file': '/dali/lgrandi/xenonnt/simulations/testing_epix_wfsim/tpc_and_nveto_cryoneutrons_200.root', 'detector': 'XENONnT', 'detector_config_override': '', 'cut_by_eventid': False, 'entry_start': None, 'entry_stop': None, 'micro_separation': 0.05, 'micro_separation_time': 10, 'tag_cluster_by': 'time', 'max_delay': 10000000.0, 'source_rate': 100.0, 'job_number': 0, 'output_path': '', 'debug': True, 'path': '/dali/lgrandi/xenonnt/simulations/testing_epix_wfsim', 'file_name': 'tpc_and_nveto_cryoneutrons_200.root', 'detector_config': [<epix.detector_volumes.SensitiveVolume object at 0x7ff4a65bea60>, <epix.detector_volumes.SensitiveVolume object at 0x7ff472db6fd0>], 'outer_cylinder': {'max_z': 7.3936, 'min_z': -154.6555, 'max_r': 66.4}}
Total entries in input file = 198
It took 1.7392 sec to load data.
Finding clusters of interactions with a dr = 0.05 mm and dt = 10 ns
It took 3.3007 sec to find clusters.
It took 3.2324 sec to merge clusters.
Removing clusters not in volumes: TPC BelowCathode
Number of clusters before: 2713
Number of clusters after: 2696
Assigning electric field to clusters
Generating photons and electrons for events
/home/mzks/official/epix/epix/quanta_generation.py:49: UserWarning: Energy deposition of 263.45294284820557 keV beyond NEST validity for NR model of 200 keV - Remove Interaction
  warnings.warn(f"Energy deposition of {en} keV beyond NEST validity for NR model of 200 keV - Remove Interaction")
/opt/XENONnT/anaconda/envs/XENONnT_development/lib/python3.8/site-packages/numpy/lib/function_base.py:2192: RuntimeWarning: invalid value encountered in ? (vectorized)
  outputs = ufunc(*inputs)
It took 0.8618 sec to get quanta.
Fixed event rate of 100.0 Hz
Done
Instructions saved to  tpc_and_nveto_cryoneutrons_200_wfsim_instructions.csv
It took 9.3694 sec to run epix.

$cat tpc_and_nveto_cryoneutrons_200_wfsim_instructions.csv
event_number,type,time,x,y,z,amp,recoil,e_dep,g4id,vol_id,local_field,n_excitons
0,1,64181824,19.575397,-62.51091,-38.670742,55,0,70.69686,11,1,200.0,54
0,1,64181824,21.687685,-58.404762,-38.932274,2834,8,58.00263,11,1,200.0,656
0,1,64181824,21.694613,-58.412357,-38.910706,3072,8,86.333595,11,1,200.0,897
0,1,64181824,21.70312,-58.441692,-38.907326,3908,8,94.53536,11,1,200.0,1113
0,1,64181824,21.708265,-58.3961,-38.94692,3390,8,81.34599,11,1,200.0,968
0,1,64181824,21.71108,-58.475006,-38.907944,6115,8,192.17854,11,1,200.0,2224
0,1,64181824,21.726461,-58.388557,-38.95812,3082,8,68.82258,11,1,200.0,798
0,1,64181824,21.746391,-58.384266,-38.96926,3888,8,98.20123,11,1,200.0,1151
0,1,64181824,21.76272,-58.39527,-38.949295,1648,7,29.4521,11,1,200.0,331
0,1,64181824,21.771029,-58.371216,-38.96195,1344,8,29.61329,11,1,200.0,309
0,1,64181824,21.771358,-58.37988,-38.969658,6643,7,159.28423,11,1,200.0,1745
0,1,64181824,21.782637,-58.37206,-38.96725,1249,8,28.531408,11,1,200.0,311
0,1,64181824,21.792103,-58.369766,-38.973232,1210,8,24.382124,11,1,200.0,255
0,1,64181824,21.802944,-58.372524,-38.975697,1870,8,35.91909,11,1,200.0,379
0,1,64181824,21.810009,-58.36582,-38.977016,1669,8,32.962933,11,1,200.0,404
0,1,64181824,21.816223,-58.35902,-38.979813,1437,8,31.141449,11,1,200.0,352
0,1,64181824,21.821163,-58.37088,-38.980762,1770,7,33.305866,11,1,200.0,388
0,1,64181824,21.82376,-58.356403,-38.982723,1343,8,26.706757,11,1,200.0,294
0,1,64181824,21.827198,-58.352818,-38.98878,870,8,17.632639,11,1,200.0,163
0,1,64181824,21.830738,-58.34713,-38.988533,1071,8,22.040703,11,1,200.0,221
0,1,64181824,21.835209,-58.3467,-38.98356,1148,8,21.878912,11,1,200.0,218
0,1,64181824,21.838873,-58.352108,-38.9808,1185,8,24.81744,11,1,200.0,247
0,1,64181824,21.842493,-58.356762,-38.979088,776,8,15.675897,11,1,200.0,135
0,1,64181824,21.846333,-58.365196,-38.967766,791,8,15.274351,11,1,200.0,130
0,1,64181824,21.846453,-58.377373,-38.95984,3765,8,83.938934,11,1,200.0,1007
0,1,64181824,21.847052,-58.358234,-38.975304,767,8,17.622437,11,1,200.0,170
0,1,64181824,21.848091,-58.36068,-38.97047,678,8,15.706513,11,1,200.0,140
0,1,64181824,21.855383,-58.36839,-38.94614,9332,8,323.62848,11,1,200.0,3613
0,1,64181824,27.006512,-48.68472,-39.216682,2404,7,46.001434,11,1,200.0,508
0,1,64181824,27.009596,-48.65983,-39.18694,6157,7,114.265724,11,1,200.0,1303
0,1,64181824,27.329262,-50.365635,-39.742687,101,7,2.6750844,11,1,200.0,4
0,1,64181824,27.72667,-49.465767,-39.161858,2285,7,42.793633,11,1,200.0,478
0,1,64181824,27.745577,-48.79939,-38.96093,4884,7,91.44903,11,1,200.0,1048
0,1,64181824,27.758389,-49.84528,-39.81129,8603,7,213.81401,11,1,200.0,2446
0,2,64181824,19.575397,-62.51091,-38.670742,8,0,70.69686,11,1,200.0,0
0,2,64181824,21.687685,-58.404762,-38.932274,1369,8,58.00263,11,1,200.0,0
0,2,64181824,21.694613,-58.412357,-38.910706,3280,8,86.333595,11,1,200.0,0
0,2,64181824,21.70312,-58.441692,-38.907326,3200,8,94.53536,11,1,200.0,0
0,2,64181824,21.708265,-58.3961,-38.94692,2588,8,81.34599,11,1,200.0,0
0,2,64181824,21.71108,-58.475006,-38.907944,8227,8,192.17854,11,1,200.0,0
0,2,64181824,21.726461,-58.388557,-38.95812,2005,8,68.82258,11,1,200.0,0
0,2,64181824,21.746391,-58.384266,-38.96926,3484,8,98.20123,11,1,200.0,0
0,2,64181824,21.76272,-58.39527,-38.949295,555,7,29.4521,11,1,200.0,0
0,2,64181824,21.771029,-58.371216,-38.96195,821,8,29.61329,11,1,200.0,0
0,2,64181824,21.771358,-58.37988,-38.969658,5156,7,159.28423,11,1,200.0,0
0,2,64181824,21.782637,-58.37206,-38.96725,776,8,28.531408,11,1,200.0,0
0,2,64181824,21.792103,-58.369766,-38.973232,486,8,24.382124,11,1,200.0,0
0,2,64181824,21.802944,-58.372524,-38.975697,752,8,35.91909,11,1,200.0,0
0,2,64181824,21.810009,-58.36582,-38.977016,785,8,32.962933,11,1,200.0,0
0,2,64181824,21.816223,-58.35902,-38.979813,868,8,31.141449,11,1,200.0,0
0,2,64181824,21.821163,-58.37088,-38.980762,661,7,33.305866,11,1,200.0,0
0,2,64181824,21.82376,-58.356403,-38.982723,679,8,26.706757,11,1,200.0,0
0,2,64181824,21.827198,-58.352818,-38.98878,398,8,17.632639,11,1,200.0,0
0,2,64181824,21.830738,-58.34713,-38.988533,504,8,22.040703,11,1,200.0,0
0,2,64181824,21.835209,-58.3467,-38.98356,483,8,21.878912,11,1,200.0,0
0,2,64181824,21.838873,-58.352108,-38.9808,620,8,24.81744,11,1,200.0,0
0,2,64181824,21.842493,-58.356762,-38.979088,341,8,15.675897,11,1,200.0,0
0,2,64181824,21.846333,-58.365196,-38.967766,329,8,15.274351,11,1,200.0,0
0,2,64181824,21.846453,-58.377373,-38.95984,2647,8,83.938934,11,1,200.0,0
0,2,64181824,21.847052,-58.358234,-38.975304,532,8,17.622437,11,1,200.0,0
0,2,64181824,21.848091,-58.36068,-38.97047,510,8,15.706513,11,1,200.0,0
0,2,64181824,21.855383,-58.36839,-38.94614,14530,8,323.62848,11,1,200.0,0
0,2,64181824,27.006512,-48.68472,-39.216682,1026,7,46.001434,11,1,200.0,0
0,2,64181824,27.009596,-48.65983,-39.18694,2269,7,114.265724,11,1,200.0,0
0,2,64181824,27.329262,-50.365635,-39.742687,101,7,2.6750844,11,1,200.0,0
0,2,64181824,27.72667,-49.465767,-39.161858,878,7,42.793633,11,1,200.0,0
0,2,64181824,27.745577,-48.79939,-38.96093,1831,7,91.44903,11,1,200.0,0
0,2,64181824,27.758389,-49.84528,-39.81129,7304,7,213.81401,11,1,200.0,0

No lines are omitted. It takes care about the only first event.

Conclusions / guess
I haven't go into the function but I guess it may come from nestpy or its versions. I couldn't find the cause.
Could someone provide an example which WFSim works well with the latest epix?

Slower code after fixing clustering index issue (#45)

Let's keep in mind to try to optimize this at some point.

The issue was spotted in #41 and resolved in #45, renouncing to use the (faster) groupby method for the index sorting after clustering.

Classification of gamma interactions

As for now, clusters originating from gamma interactions are classified as type 'gamma' by the type of the particle in the first interaction in time in the cluster. This type in NEST ("gammaRay") does not include Compton scattering:

https://github.com/NESTCollaboration/nestpy/blob/c6cadad63254f1d34d10c053d105e9207e81f111/src/nestpy/NEST.hh#L118

recoil types

Since #10 we are using the NEST ID as recoil type identifier. In WFsim theses IDs are mapped to NR, ER or Alpha in https://github.com/XENONnT/WFSim/blob/master/wfsim/core.py#L22. I ran some neutron simulations and got in about 5% of the generated instructions the NEST ID 12. This case is used in epix if the interaction does not match any implemented classification.

module 'awkward' has no attribute 'to_pandas' since awkward version 2.0

The aw.to_pandas function was annoyingly replaced with aw.to_dataframe in awkward version 2.0. This is used in clustering.py.

Kr83m deltaT_ns

When/If this issue (NESTCollaboration/nest#139) in NEST is resolved we may need to change the epix code to calculate deltaT_ns and pass it to NEST.

Numba - Awkward Array Version Problem

Looks like @cfuselli found a problem with numba and/or awkward:

(XENONnT_development) ../epix/bin/run_epix --InputFile /home/cfuselli/software/mc/events.root --Debug
*** Detector definition message ***
You are currently using the default XENON10 template detector.

Using nestpy version 2.0.0
epix configuration:  {'input_file': '/home/cfuselli/software/mc/events.root', 'detector': 'XENONnT', 'detector_config_override': '', 'cut_by_eventid': False, 'nr_only': False, 'entry_start': None, 'entry_stop': None, 'micro_separation': 0.005, 'micro_separation_time': 10, 'tag_cluster_by': 'energy', 'max_delay': 10000000.0, 'source_rate': 0, 'yield': 'nest', 'job_number': 0, 'output_path': '', 'debug': True, 'path': '/home/cfuselli/software/mc', 'file_name': 'events.root', 'detector_config': [<epix.detector_volumes.SensitiveVolume object at 0x7f14ae0e9ac0>, <epix.detector_volumes.SensitiveVolume object at 0x7f14ae0e9d60>], 'outer_cylinder': {'max_z': 7.3936, 'min_z': -154.6555, 'max_r': 66.4}}
Total entries in input file = 500
It took 48.4048 sec to load data.
Finding clusters of interactions with a dr = 0.005 mm and dt = 10 ns
/opt/XENONnT/anaconda/envs/XENONnT_development/lib/python3.8/site-packages/awkward/_connect/_numba/__init__.py:28: DeprecationWarning: distutils Version classes are deprecated. Use packaging.version instead.
  if not checked_version and distutils.version.LooseVersion(
/opt/XENONnT/anaconda/envs/XENONnT_development/lib/python3.8/site-packages/awkward/_connect/_numba/__init__.py:30: DeprecationWarning: distutils Version classes are deprecated. Use packaging.version instead.
  ) < distutils.version.LooseVersion("0.50"):
It took 10.0623 sec to find clusters.
Traceback (most recent call last):
  File "../epix/bin/run_epix", line 71, in <module>
    epix.run_epix.main(args, return_df=True)
  File "/home/cfuselli/software/epix/epix/run_epix.py", line 45, in main
    result = epix.cluster(inter, args['tag_cluster_by'] == 'energy')
  File "/home/cfuselli/software/epix/epix/clustering.py", line 164, in cluster
    _cluster(x, y, z, ed, time, ci,
  File "/opt/XENONnT/anaconda/envs/XENONnT_development/lib/python3.8/site-packages/numba/core/dispatcher.py", line 468, in _compile_for_args
    error_rewrite(e, 'typing')
  File "/opt/XENONnT/anaconda/envs/XENONnT_development/lib/python3.8/site-packages/numba/core/dispatcher.py", line 409, in error_rewrite
    raise e.with_traceback(None)
numba.core.errors.TypingError: Failed in nopython mode pipeline (step: nopython frontend)
Internal error at <numba.core.typeinfer.CallConstraint object at 0x7f1475f41460>.
Failed in nopython mode pipeline (step: native lowering)
module 'llvmlite' has no attribute 'llvmpy'

File "../epix/epix/clustering.py", line 302:
def _write_result(res, x_mean, y_mean, z_mean,
    <source elided>
    res.begin_record()
    res.field('x')
    ^

During: lowering "$16call_method.6 = call $12load_method.4($const14.5, func=$12load_method.4, args=[Var($const14.5, clustering.py:302)], kws=(), vararg=None, varkwarg=None, target=None)" at /home/cfuselli/software/epix/epix/clustering.py (302)
During: resolving callee type: type(CPUDispatcher(<function _write_result at 0x7f14ae0d05e0>))
During: typing of call at /home/cfuselli/software/epix/epix/clustering.py (211)

Enable logging at debug level for details.

File "../epix/epix/clustering.py", line 211:
def _cluster(x, y, z, ed, time, ci,
    <source elided>
                # Write result, simple but extensive with awkward...
                _write_result(res, x_mean, y_mean, z_mean,
                ^

Problem with the cluster indexing after DBSCAN

After DBSCAN is applied in epix, the clusters are indexed for further merging:

epix/epix/clustering.py

Lines 33 to 41 in 53d48c8

 # Splitting into individual events and time cluster and apply space clustering space: 

 groups = df.groupby(['entry', 'time_cluster']) 

 groups = groups.apply(lambda x: _find_cluster(x, cluster_size_space)) 

 for i in np.unique(groups.index.get_level_values(0)): 

 add_to_cluster = 0 

 for j in range(len(groups[i])): 

 groups[i][j]+=add_to_cluster 

 add_to_cluster = np.max(groups[i][j])+1

An order of cluster indices in the column which is created here does not match the order of rows in the dataframe with interactions, which creates a problem in this line:

epix/epix/clustering.py

Line 43 in 53d48c8

df['cluster_id'] = np.concatenate(groups.values)

Therefore we have a lot of cases where a cluster contains interactions which do not belong to this cluster. Since a weighted averaging is applied during cluster merging, this issue is not that noticeable in the resulting output. A detailed description of the issue and a possible solution is here:

/dali/lgrandi/pkavrigin/2021-08-25_FixEpixDemo/2021-08-25_FixForEpix_Demo.ipynb

NB: The notebook requires 'dbg_out' branch of epix if you want to re-run it.

Propagate LXe level from Geant4 file

Modify this value according to information stored in the ROOT file, if different than default.

Significant fraction of time is just spent in loading text objects

Just documenting something that is well-known since the early development stage, but is worth reporting for future upgrades: epix is fast and clean, but still uproot takes a long while loading the Geant4 branches that are stored in text format (type of interaction, energy depositing process, etc.). This was discussed with the uproot and awkward developer Jim Pivarski: https://gitter.im/Scikit-HEP/uproot?at=5fad4a7bc10273610aff87d1

And attaching a very recent output from @HenningSE for 100k, in which the file loading time is about 40% of the total:

Total entries in input file = 500000
Starting to read from output file entry 0
Ending read in at output file entry 100000
It took 182.5625 sec to load data.
Finding clusters of interactions with a dr = 0.05 mm and dt = 10.0 ns
It took 293.2970 sec to find clusters.
It took 6.8375 sec to merge clusters.
Removing clusters not in volumes: TPC BelowCathode
Number of clusters before: 127425
Number of clusters after: 126940
Assigning electric field to clusters
Generating photons and electrons for events
It took 5.7376 sec to get quanta.
Clean event separation
Min. S2 amp BEFORE MACRO-CLUSTERING: 1
Min. S2 amp AFTER MACRO-CLUSTERING: 1
Source finished!

An idea floating around at some point was to load data in uproot chunks and aim at some parallelization, although never explored in detail. We might want to reconsider if running into memory issues.

Typo in the Kr83m energy (9.4 keV)

I believe the following energy should be 9.4 keV, not 9.1 keV. I will fix this once I can talk to the experts since this is my first time looking at this.

Also, the comparison for nest to make equal SHOULD be fixed in the upcoming patch, so let's keep this issue open for a couple days in the hopes that NEST versions up soon and we can fix all these issues together.

epix/epix/quanta_generation.py

Line 45 in e3cc08f

if abs(en - 9.1) > max_allowed_energy_difference:

maxTime in NEST Kr83m model

In epix we set A and Z to zero in case they are not explicitly needed by e.g. the NEST Ion model. Inside NEST, A and Z are the variables massNum and atomNum. These variables are only explicitly used in the Ion model: https://github.com/NESTCollaboration/nest/blob/master/src/NEST.cpp#L718 or the NR model: https://github.com/NESTCollaboration/nest/blob/master/src/NEST.cpp#L713.

For some reason massNum is used as placeholder for a variable maxTimeSeparation when using the Kr83m model: https://github.com/NESTCollaboration/nest/blob/master/src/NEST.cpp#L723-L725. When A is set to 0 in epix the NEST Kr83m call freezes.

Can we just set A to an value other than 0?
Do we need to calculate a proper time separation?

Event Separation

At the moment EPIX separates all events by a fixed time. This allows "clean" simulations since effects like pile-up can not occur. I would propose to give the user of EPIX more flexibility in the event separation with these three options:

"Clean" simulations: Separation of events with a fixed time. This is basically the way it is implemented at the moment.
Event Rate: The user gives an event rate as argument to EPIX and the events are distributed randomly (uniform distribution) according to the rate.
Variable Rate: Maybe it is interesting for some users to simulate sources with not constant event rates (e.g. source opening, source closing ..). The user could give an csv file to EPIX containing informations about the event rate over time.

fast_simulator - Macro-clustering issue

This issue is related to 'fast_simulator' branch only. Apparently the macro-clustering function ('helpers.macro_cluster_events') can produce negative number of quanta after the clustering. You can see it in the 'fast_sim_test' branch, where I've added a notebook for tests with AmBe data. In this notebook you can look for "Min. S2 amp AFTER MACRO-CLUSTERING" in the output, which is the minimum value of ['amp'] for ['type'==2] in the epix instructions, and it is a negative value which naturally leads to an error. I guess the problem is this line:

epix/epix/simulator/helpers.py

Line 77 in 67b14c8

instructions[ix1]['amp']=-1 #flag to throw this instruction away later

This value is probably used later in the inner loop.

Chunk for large files

As @petergaemers pointed out EPIX has hard times to handle large files. We have to add some sort of chunking.

	# Splitting into individual events and time cluster and apply space clustering space:
	groups = df.groupby(['entry', 'time_cluster'])
	groups = groups.apply(lambda x: _find_cluster(x, cluster_size_space))

	for i in np.unique(groups.index.get_level_values(0)):
	add_to_cluster = 0
	for j in range(len(groups[i])):
	groups[i][j]+=add_to_cluster
	add_to_cluster = np.max(groups[i][j])+1

xenonnt / epix Goto Github PK

epix's Issues

To reproduce:

Traceback

Recommend Projects

Recommend Topics

Recommend Org