diff --git a/Notebooks/Simulation_Refactor.ipynb b/Notebooks/Simulation_Refactor.ipynb index 117a8b57..ce70c2dc 100644 --- a/Notebooks/Simulation_Refactor.ipynb +++ b/Notebooks/Simulation_Refactor.ipynb @@ -10,6 +10,7 @@ ] }, { + "attachments": {}, "cell_type": "markdown", "id": "6e4fcd89-bd8d-4e5e-8fe3-ed133fa19874", "metadata": {}, @@ -25,8 +26,12 @@ "outputs": [], "source": [ "import strax\n", + "import straxen\n", "import cutax\n", - "import fuse" + "import fuse\n", + "\n", + "import numpy as np\n", + "from straxen import URLConfig" ] }, { @@ -39,14 +44,13 @@ "\n", "In fuse all simulation steps are handled in dedicated plugins just like straxen does for our data processing. The first part of our full chain simulations is the simulation of microphysics effects. The corresponding plugins are grouped in the `micro_physics` module.\n", "\n", - "To set up a simulation we first need to define a `Context`. At this stage we dont't need a lot of inputs so we can run the simulation inside a generic `strax.Context`. Make sure to register all needed plugins and define an `DataDirectory`. Afterwards we can change the configuration of the simulation using `st.set_config`. If you want to try your own `.root` or `.csv`file feel free to change the given example file.\n", + "To set up a simulation we first need to define a `Context`. At this stage we dont't need a lot of inputs so we can run the simulation inside a generic `strax.Context`. fuse has a predefined `microphyics_context` that registers all needed plugins for us. Afterwards we can change the configuration of the simulation using `st.set_config`. If you want to try your own `.root` or `.csv`file feel free to change the given example file.\n", "\n", "There are some more options we can change: \n", "- `debug` will print out some debug informations during the simulation (For now there are basically no debug statements.)\n", "- `source_rate`: You can specify a rate in Hz. fuse will distribute your events according to this rate. \n", "- `n_interactions_per_chunk`: fuse builds strax chunks based on an algorithm called 'dynamic chunking'. It searches for empty time intervalls and cuts the data into chunks there. With this parameter you can set an lower limit to the number of interactions per chunk. For Geant4 files something in the order of 10k - 100k can work. For csv files i would reccomend something in the order of 10 - 1000 depending on you specific data. \n", - "- `cut_delayed`: fuse will remove interactions that happen after the last sampled time + the cut_delayed time. We dont want to deal with interactions that happen way after our \"run\" finished (yet)\n", - "- `DetectorConfigOverride`: Part of the epix detector config - we will update this part at some point. (Here you can set your electric field for the microphysics simulation)" + "- `cut_delayed`: fuse will remove interactions that happen after the last sampled time + the cut_delayed time. We dont want to deal with interactions that happen way after our \"run\" finished (yet)" ] }, { @@ -56,23 +60,12 @@ "metadata": {}, "outputs": [], "source": [ - "st = strax.Context(register = [fuse.micro_physics.ChunkInput,\n", - " fuse.micro_physics.FindCluster,\n", - " fuse.micro_physics.MergeCluster,\n", - " fuse.micro_physics.ElectricField,\n", - " fuse.micro_physics.NestYields,\n", - " #fuse.micro_physics.BBFYields,\n", - " fuse.micro_physics.MicroPhysicsSummary],\n", - " storage = [strax.DataDirectory('path/to/output/folder')]\n", - " )\n", + "st = fuse.context.microphysics_context(\"/path/to/output/folder\")\n", "\n", "st.set_config({\"path\": \"/project2/lgrandi/xenonnt/simulations/testing\",\n", " \"file_name\": \"pmt_neutrons_100.root\",\n", " \"debug\": True,\n", " \"source_rate\": 1,\n", - " #\"n_interactions_per_chunk\": 100,\n", - " \"cut_delayed\": 4e14,\n", - " \"detector_config_override\" : \"sr0_epix_detectorconfig.ini\",\n", " })" ] }, @@ -82,7 +75,7 @@ "id": "f9c7149d-b2bb-4ced-987b-f4ab82f3e0c0", "metadata": {}, "source": [ - "Now that we have our context we can start to run fuse. Just use `st.make`. Probably there will be a lot of warnings..." + "Now that we have our context we can start to run fuse. Just use `st.make` and pick the data types you would like to make." ] }, { @@ -105,6 +98,7 @@ ] }, { + "attachments": {}, "cell_type": "markdown", "id": "655c40ff-a51c-44e6-8321-d21098b5094c", "metadata": {}, @@ -124,6 +118,7 @@ ] }, { + "attachments": {}, "cell_type": "markdown", "id": "5eaa1305-baf8-4675-bbdc-3810a5afb796", "metadata": {}, @@ -143,6 +138,7 @@ ] }, { + "attachments": {}, "cell_type": "markdown", "id": "8eb4b98d-2d8c-45f6-88b3-15dde93ca3f3", "metadata": {}, @@ -151,43 +147,32 @@ ] }, { + "attachments": {}, "cell_type": "markdown", "id": "b4902b2c-2ce6-445d-84b5-35e5fb7e68be", "metadata": {}, "source": [ - "Now that you know how to run microphysics simulations we can continue with the detector physics, PMT and DAQ simulations. The corresponding plugins can be found in `detector_physics` and `pmt_and_daq`. " + "Now that you know how to run microphysics simulations we can continue with the detector physics, PMT and DAQ simulations. The corresponding plugins can be found in the`detector_physics` and `pmt_and_daq` modules. " ] }, { - "cell_type": "code", - "execution_count": null, - "id": "2a987408-bbe2-413b-9f38-4727bb566157", + "attachments": {}, + "cell_type": "markdown", + "id": "a68764f4-72a2-47eb-be3f-36cf9f3761d2", "metadata": {}, - "outputs": [], "source": [ - "full_chain_modules = [fuse.micro_physics.ChunkInput,\n", - " fuse.micro_physics.FindCluster,\n", - " fuse.micro_physics.MergeCluster,\n", - " fuse.micro_physics.ElectricField,\n", - " fuse.micro_physics.NestYields,\n", - " fuse.micro_physics.MicroPhysicsSummary,\n", - " fuse.detector_physics.S1PhotonPropagation,\n", - " fuse.detector_physics.ElectronDrift,\n", - " fuse.detector_physics.ElectronExtraction,\n", - " fuse.detector_physics.ElectronTiming,\n", - " fuse.detector_physics.SecondaryScintillation,\n", - " fuse.detector_physics.S2PhotonPropagation,\n", - " fuse.pmt_and_daq.PMTAfterPulses,\n", - " fuse.pmt_and_daq.PMTResponseAndDAQ,\n", - " ]" + "We want to simulate now up to the raw_records level and then process the data with our usual straxen processing routines. To do so, we will use the `full_chain_context` that is prepared in fuse. To build the context you need to provide a config dictionary, we will load our XENONnT SR0 config here. Remember to select an output folder with sufficient free disk space as `propagated_s2_photons` and `raw_records` data can be a bit heavy. " ] }, { - "cell_type": "markdown", - "id": "a68764f4-72a2-47eb-be3f-36cf9f3761d2", + "cell_type": "code", + "execution_count": null, + "id": "ad9c4eff", "metadata": {}, + "outputs": [], "source": [ - "We want to simulate now up to the raw_records level and then process the data with our usual straxen processing routines. To do so, we will now hijack the simulation context from cutax and register our modules. In the future we aim for a dedicated context where we don't need to do this manualy. Remember to select an output folder with sufficient free disk space as raw_records data can be a bit heavy. " + "url_string = 'simple_load://resource://format://fax_config_nt_sr0_v4.json?&fmt=json'\n", + "config = URLConfig.evaluate_dry(url_string)" ] }, { @@ -197,28 +182,24 @@ "metadata": {}, "outputs": [], "source": [ - "st = cutax.contexts.xenonnt_sim_SR0v3_cmt_v9(output_folder = \"path/to/output/folder\")\n", - "\n", - "for module in full_chain_modules:\n", - " st.register(module)\n", + "st = fuse.context.full_chain_context(out_dir = \"/path/to/output/folder\",\n", + " config = config)\n", "\n", "st.set_config({\"path\": \"/project2/lgrandi/xenonnt/simulations/testing\",\n", " \"file_name\": \"pmt_neutrons_100.root\",\n", " \"debug\": True,\n", " \"source_rate\": 1,\n", - " #\"n_interactions_per_chunk\": 100,\n", - " \"cut_delayed\": 4e14,\n", - " \"detector_config_override\" : \"sr0_epix_detectorconfig.ini\",\n", " })" ] }, { + "attachments": {}, "cell_type": "markdown", "id": "c46efb1e-5bb9-45d6-9d2c-a16076b35423", "metadata": {}, "source": [ "### Microphysics\n", - "Just like before we first run the microphysics simulation. Lets make `microphysics_summary` first. " + "Just like before we first run the microphysics simulation. Lets make `microphysics_summary` first. Note that you dont need to make every data type manually but let strax take care of it." ] }, { @@ -241,7 +222,26 @@ "source": [ "### S1 Simulation\n", "\n", - "After the microphysics simulation is done we can continue and distribute our S1 photons to our PMTs. The output will be a long list where each row represents a single photon with the information of `time` and `channel`attached. Strax can show you a progress bar during simulation. I found it to be a bit buggy from times in combination with fuse. " + "After the microphysics simulation is done we can continue and distribute our S1 photons to our PMTs. First we need to determin how many photons are actually recorded at an PMT. This step is quite fast and can be used in fastsim in the future (I guess)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "358033c6", + "metadata": {}, + "outputs": [], + "source": [ + "st.make(run_number,\"s1_photons\" ,progress_bar = True)" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "25aec3a3", + "metadata": {}, + "source": [ + "In a second step we distribute the photons to our PMTs. The output will be a long list where each row represents a single photon with the information of `time` and `channel`attached. Strax can show you a progress bar during simulation. I found it to be a bit buggy misleading as it shows the fraction of simulated time and not number of interactions. " ] }, { @@ -255,13 +255,14 @@ ] }, { + "attachments": {}, "cell_type": "markdown", "id": "0623b54a-fa02-4b29-bf88-a1c4c8929f01", "metadata": {}, "source": [ "### S2 Simulations\n", "\n", - "The S2 simulations are a bit more complicated than the S1 simulation. We first need to drift our electrons to the liquid-gas interface, extract them and calculate the timing of the electrons when reaching the gas phase. Afterwards we can simulate how many photons each electron generates and then distribute the photons on the PMTs. The output of the S2 simulation has basically the same format as the S1 simulation output. " + "The S2 simulations are a bit more complicated than the S1 simulation and thus split into more plugins. We first need to drift our electrons to the liquid-gas interface, extract them and calculate the timing of the electrons when reaching the gas phase. Afterwards we can simulate how many photons each electron generates and then distribute the photons on the PMTs. The output of the S2 simulation has basically the same format as the S1 simulation output. " ] }, { @@ -280,6 +281,7 @@ ] }, { + "attachments": {}, "cell_type": "markdown", "id": "b33903bf-7ad2-46ea-9763-f6ab716c0e37", "metadata": {}, @@ -300,6 +302,7 @@ ] }, { + "attachments": {}, "cell_type": "markdown", "id": "a4dd9cbc-39df-40cf-8695-2251dd268b3a", "metadata": {}, @@ -320,6 +323,7 @@ ] }, { + "attachments": {}, "cell_type": "markdown", "id": "e1a534dc-6449-46f5-86e9-f477a280d1a1", "metadata": {}, @@ -328,6 +332,7 @@ ] }, { + "attachments": {}, "cell_type": "markdown", "id": "b2cad900-8ae6-4fe3-a0fe-479f529c8e43", "metadata": {}, @@ -365,6 +370,15 @@ "event_info.head(20)" ] }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "c46aebc0", + "metadata": {}, + "source": [ + "If you take a look at the readme of fuse you will find that `S1PhotonHits`, `ElectronDrift`, `ElectronExtraction` and `SecondaryScintillation` produce data with the `interactions_in_roi` data_kind just like e.g. the microphysics_summary. As a result we can load this data togeter. For each interaction cluster we thus have the number of S1 and S2 photons. This can be exploited when integrating fastsim. " + ] + }, { "attachments": {}, "cell_type": "markdown", @@ -485,7 +499,7 @@ "id": "2171b669", "metadata": {}, "source": [ - "We can now set upt the context like before. fuse will automatically detect that a csv file is given and treat it accordingly. As our csv file contains a single energy deposit per event, we should set `n_interactions_per_chunk` to something reasonable like 100. Now each chunk will consist of at least 100 events. " + "We can use the full_chain_context again as fuse will automatically detect that a csv file is given and treat it accordingly. As our csv file contains a single energy deposit per event, we should set `n_interactions_per_chunk` to something reasonable like 100. Now each chunk will consist of at least 100 events. " ] }, { @@ -495,10 +509,8 @@ "metadata": {}, "outputs": [], "source": [ - "st = cutax.contexts.xenonnt_sim_SR0v3_cmt_v9(output_folder = \"path/to/output/folder\")\n", - "\n", - "for module in full_chain_modules:\n", - " st.register(module)\n", + "st = fuse.context.full_chain_context(out_dir = \"/path/to/output/folder\",\n", + " config = config)\n", "\n", "st.set_config({\"path\": \".\",\n", " \"file_name\": \"monoenergetic_200keV.csv\",\n", @@ -506,7 +518,6 @@ " \"source_rate\": 1,\n", " \"n_interactions_per_chunk\": 100,\n", " \"cut_delayed\": 4e14,\n", - " \"detector_config_override\" : \"sr0_epix_detectorconfig.ini\",\n", " })\n", "\n", "run_number = \"00001\"" @@ -658,7 +669,7 @@ "id": "d91a353e", "metadata": {}, "source": [ - "We can now set up the simulation context. As we are just running the detector physics simulation we do not need to register the `micro_pyhsics` plugins but need to register the `ChunkCsvInput` plugin." + "We can now set up the simulation context. As i was to lazy to set up a context for this use case yet i just hijack the `full_chain_context` and register the `ChunkCsvInput` plugin. " ] }, { @@ -668,23 +679,11 @@ "metadata": {}, "outputs": [], "source": [ - "st = cutax.contexts.xenonnt_sim_SR0v3_cmt_v9(output_folder = \"path/to/output/folder\")\n", - "\n", - "detector_simulation_modules = [\n", - " fuse.detector_physics.ChunkCsvInput,\n", - " fuse.detector_physics.S1PhotonPropagation,\n", - " fuse.detector_physics.ElectronDrift,\n", - " fuse.detector_physics.ElectronExtraction,\n", - " fuse.detector_physics.ElectronTiming,\n", - " fuse.detector_physics.SecondaryScintillation,\n", - " fuse.detector_physics.S2PhotonPropagation,\n", - " fuse.pmt_and_daq.PMTAfterPulses,\n", - " fuse.pmt_and_daq.PMTResponseAndDAQ,\n", - " ]\n", + "st = fuse.context.full_chain_context(out_dir = \"/path/to/output/folder\",\n", + " config = config)\n", "\n", "\n", - "for module in detector_simulation_modules:\n", - " st.register(module)\n", + "st.register(fuse.detector_physics.ChunkCsvInput)\n", "\n", "st.set_config({\"input_file\": \"./random_detectorphysics_instructions.csv\",\n", " \"debug\": True,\n", diff --git a/Notebooks/sr0_epix_detectorconfig.ini b/Notebooks/sr0_epix_detectorconfig.ini deleted file mode 100644 index c59c4207..00000000 --- a/Notebooks/sr0_epix_detectorconfig.ini +++ /dev/null @@ -1,20 +0,0 @@ -# Settings for detector volumes -# If defaults should be used remove field from file. -[TPC] -# If interactions in this volume should be stored: -to_be_stored = True -# Electric field map to be used can be either a number (constant field -# in V/cm) or a path to a field map file. -electric_field = ../../private_nt_aux_files/sim_files/fieldmap_2D_B2d75n_C2d75n_G0d3p_A4d9p_T0d9n_PMTs1d3n_FSR0d65p_QPTFE_0d5n_0d4p.json.gz -# Electric field [V/cm] for interactions outside of the map if any: -efield_outside_map = 23 -# Xe density in given volume, needed for quanta generation: -xe_density = 2.862 - -[BelowCathode] -to_be_stored = True -electric_field = ../../private_nt_aux_files/sim_files/fieldmap_2D_B2d75n_C2d75n_G0d3p_A4d9p_T0d9n_PMTs1d3n_FSR0d65p_QPTFE_0d5n_0d4p.json.gz -xe_density = 2.862 - -[GasPhase] - diff --git a/README.md b/README.md index a3471544..5bfd796d 100644 --- a/README.md +++ b/README.md @@ -12,12 +12,11 @@ At the moment the intallation procedure is not very advanced. I would recommend 1. Clone the fuse repository. 2. Clone the private_nt_aux_files repository to the same directory as you cloned fuse. -3. Install fuse using `pip install -e .` in the fuse directory. +2. Install fuse using `pip install -e .` in the fuse directory. ## Plugin Structure The full simulation chain in split into multiple plugins. An overview of the simulation structure can be found below. -![Simulation_Refactor_Plugins](https://user-images.githubusercontent.com/27280678/235156990-7fa63aae-21c4-45b5-9b71-a42a4173f0da.jpg) - +![Simulation Refactor](https://github.com/XENONnT/fuse/assets/27280678/2604ff67-ae7b-4d5b-968a-206af6d3e34a) diff --git a/fuse/__init__.py b/fuse/__init__.py index a843e29d..1b8885b8 100644 --- a/fuse/__init__.py +++ b/fuse/__init__.py @@ -1,2 +1,7 @@ from . import plugins -from .plugins import * \ No newline at end of file +from .plugins import * + +from .context import * + +from .vertical_merger_plugin import * +from .volume_plugin import * \ No newline at end of file diff --git a/fuse/common.py b/fuse/common.py index 90df8b6e..bdc56c0b 100644 --- a/fuse/common.py +++ b/fuse/common.py @@ -1,9 +1,34 @@ import numpy as np import awkward as ak -import straxen -import strax import numba -from copy import deepcopy + + +@numba.njit() +def dynamic_chunking(data, scale, n_min): + + idx_sort = np.argsort(data) + idx_undo_sort = np.argsort(idx_sort) + + data_sorted = data[idx_sort] + + diff = data_sorted[1:] - data_sorted[:-1] + + clusters = np.array([0]) + c = 0 + for value in diff: + if value <= scale: + clusters = np.append(clusters, c) + + elif len(clusters[clusters == c]) < n_min: + clusters = np.append(clusters, c) + + elif value > scale: + c = c + 1 + clusters = np.append(clusters, c) + + clusters_undo_sort = clusters[idx_undo_sort] + + return clusters_undo_sort def full_array_to_numpy(array, dtype): @@ -17,63 +42,26 @@ def full_array_to_numpy(array, dtype): return numpy_data -#WFSim functions -def make_map(map_file, fmt=None, method='WeightedNearestNeighbors'): - """Fetch and make an instance of InterpolatingMap based on map_file - Alternatively map_file can be a list of ["constant dummy", constant: int, shape: list] - return an instance of DummyMap""" +#This is a modified version of the corresponding WFsim code.... +@numba.njit() +def uniform_to_pe_arr(p, channel, __uniform_to_pe_arr): + indices = np.int64(p * 2000) + 1 + return __uniform_to_pe_arr[channel, indices] - if isinstance(map_file, list): - assert map_file[0] == 'constant dummy', ('Alternative file input can only be ' - '("constant dummy", constant: int, shape: list') - return DummyMap(map_file[1], map_file[2]) - - elif isinstance(map_file, str): - if fmt is None: - fmt = parse_extension(map_file) - - #log.debug(f'Initialize map interpolator for file {map_file}') - map_data = straxen.get_resource(map_file, fmt=fmt) - return straxen.InterpolatingMap(map_data, method=method) - - else: - raise TypeError("Can't handle map_file except a string or a list") - -def make_patternmap(map_file, fmt=None, method='WeightedNearestNeighbors', pmt_mask=None): - """ 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 - """ - # making tests not failing, we can probably overwrite it completel - if isinstance(map_file, list): - #log.warning(f'Using dummy map with pattern mask! This has no effect here!') - assert map_file[0] == 'constant dummy', ('Alternative file input can only be ' - '("constant dummy", constant: int, shape: list') - return DummyMap(map_file[1], map_file[2]) - elif isinstance(map_file, str): - if fmt is None: - fmt = parse_extension(map_file) - map_data = deepcopy(straxen.get_resource(map_file, fmt=fmt)) - # XXX: straxed deals with pointers and caches resources, it means that resources are global - # what is bad, so we make own copy here and modify it locally - if 'compressed' in map_data: - compressor, dtype, shape = map_data['compressed'] - map_data['map'] = np.frombuffer( - strax.io.COMPRESSORS[compressor]['decompress'](map_data['map']), - dtype=dtype).reshape(*shape) - del map_data['compressed'] - if 'quantized' in map_data: - map_data['map'] = map_data['quantized']*map_data['map'].astype(np.float32) - del map_data['quantized'] - if not (pmt_mask is None): - assert (map_data['map'].shape[-1]==pmt_mask.shape[0]), "Error! Pattern map and PMT gains must have same dimensions!" - map_data['map'][..., ~pmt_mask]=0.0 - return straxen.InterpolatingMap(map_data, method=method) - else: - raise TypeError("Can't handle map_file except a string or a list") +#In WFSim uniform_to_pe_arr is called inside a loop over the channels +#I needed to change the code to run on all channels at once +@numba.njit() +def loop_uniform_to_pe_arr(p, channel, __uniform_to_pe_arr): + result = [] + for i in range(len(p)): + result.append(uniform_to_pe_arr(p[i], + channel=channel[i], + __uniform_to_pe_arr=__uniform_to_pe_arr) ) + return np.array(result) +#WFSim functions def parse_extension(name): """Get the extention from a file name. If zipped or tarred, can contain a dot""" split_name = name.split('.') diff --git a/fuse/context.py b/fuse/context.py new file mode 100644 index 00000000..1eeb1e1e --- /dev/null +++ b/fuse/context.py @@ -0,0 +1,241 @@ +# Define simulation cotexts here +# When fuse gets public this needs to be moved to e.g. cutax? + +import strax +import cutax +import straxen +import fuse +import os +import numpy as np + +from straxen import URLConfig + +#Microphysics context +def microphysics_context(out_dir): + st = strax.Context(register = [fuse.micro_physics.ChunkInput, + fuse.micro_physics.FindCluster, + fuse.micro_physics.MergeCluster, + fuse.micro_physics.XENONnT_TPC, + fuse.micro_physics.XENONnT_BelowCathode, + fuse.micro_physics.VolumesMerger, + fuse.micro_physics.ElectricField, + fuse.micro_physics.NestYields, + fuse.micro_physics.MicroPhysicsSummary], + storage = [strax.DataDirectory(out_dir)] + ) + + st.set_config({ + "efield_map": 'itp_map://resource://format://' + 'fieldmap_2D_B2d75n_C2d75n_G0d3p_A4d9p_T0d9n_PMTs1d3n_FSR0d65p_QPTFE_0d5n_0d4p.json.gz?' + '&fmt=json.gz' + '&method=RegularGridInterpolator', + #'g1_value' : 0.151, + #'g2_value' : 16.45, + #'cs1_spline_path': '/project2/lgrandi/pkavrigin/2023-04-24_epix_data_files/cs1_func_E_option2.pkl', + #'cs2_spline_path' : '/project2/lgrandi/pkavrigin/2023-04-24_epix_data_files/cs2_func_E_option2.pkl', + }) + + return st + + +#Full Chain Context +#For now lets just hijack the cutax simulations context + +full_chain_modules = [fuse.micro_physics.ChunkInput, + fuse.micro_physics.FindCluster, + fuse.micro_physics.MergeCluster, + fuse.micro_physics.XENONnT_TPC, + fuse.micro_physics.XENONnT_BelowCathode, + fuse.micro_physics.VolumesMerger, + fuse.micro_physics.ElectricField, + fuse.micro_physics.NestYields, + fuse.micro_physics.MicroPhysicsSummary, + fuse.detector_physics.S1PhotonHits, + fuse.detector_physics.S1PhotonPropagation, + fuse.detector_physics.ElectronDrift, + fuse.detector_physics.ElectronExtraction, + fuse.detector_physics.ElectronTiming, + fuse.detector_physics.SecondaryScintillation, + fuse.detector_physics.S2PhotonPropagation, + fuse.pmt_and_daq.PMTAfterPulses, + fuse.pmt_and_daq.PhotonSummary, + fuse.pmt_and_daq.PMTResponseAndDAQ, + ] + +def full_chain_context(out_dir, config): + + st = cutax.contexts.xenonnt_sim_SR0v3_cmt_v9(output_folder = out_dir) + + for module in full_chain_modules: + st.register(module) + + st.set_config({ + "detector": "XENONnT", + "efield_map": 'itp_map://resource://format://' + 'fieldmap_2D_B2d75n_C2d75n_G0d3p_A4d9p_T0d9n_PMTs1d3n_FSR0d65p_QPTFE_0d5n_0d4p.json.gz?' + '&fmt=json.gz' + '&method=RegularGridInterpolator', + "drift_velocity_liquid": config["drift_velocity_liquid"], + "drift_time_gate": config["drift_time_gate"], + "diffusion_constant_longitudinal": config["diffusion_constant_longitudinal"], + "electron_lifetime_liquid": config["electron_lifetime_liquid"], + "enable_field_dependencies": config["enable_field_dependencies"], + "tpc_length": config["tpc_length"], + "field_distortion_model": config["field_distortion_model"], + "field_dependencies_map_tmp": 'itp_map://resource://format://' + 'field_dependent_radius_depth_maps_B2d75n_C2d75n_G0d3p_A4d9p_T0d9n_PMTs1d3n_FSR0d65p_QPTFE_0d5n_0d4p.json.gz?' + '&fmt=json.gz' + '&method=RectBivariateSpline', + "diffusion_longitudinal_map_tmp":'itp_map://resource://format://' + 'data_driven_diffusion_map_XENONnTSR0V2.json.gz?' + '&fmt=json.gz' + '&method=WeightedNearestNeighbors', + "fdc_map_fuse": 'itp_map://resource://format://' + 'init_to_final_position_mapping_B2d75n_C2d75n_G0d3p_A4d9p_T0d9n_PMTs1d3n_FSR0d65p_QPTFE_0d5n_0d4p.json.gz?' + '&fmt=json.gz' + '&method=RectBivariateSpline', + "digitizer_voltage_range": config["digitizer_voltage_range"], + "digitizer_bits": config["digitizer_bits"], + "pmt_circuit_load_resistor": config["pmt_circuit_load_resistor"], + "s2_secondary_sc_gain": config["s2_secondary_sc_gain"], + "g2_mean": config["g2_mean"], + "electron_extraction_yield": config["electron_extraction_yield"], + "ext_eff_from_map": config["ext_eff_from_map"], + "se_gain_from_map": config["se_gain_from_map"], + "gains": 'pmt_gains://resource://format://' + 'to_pe_nt.npy?' + '&fmt=npy' + '&digitizer_voltage_range=plugin.digitizer_voltage_range' + '&digitizer_bits=plugin.digitizer_bits' + '&pmt_circuit_load_resistor=plugin.pmt_circuit_load_resistor', + "s2_correction_map": 'itp_map://resource://format://' + 'XENONnT_s2_xy_map_v4_210503_mlp_3_in_1_iterated.json?' + '&fmt=json', + "se_gain_map": 'itp_map://resource://format://' + 'XENONnT_se_xy_map_v1_mlp.json?' + '&fmt=json', + "s2_pattern_map": 'pattern_map://resource://format://' + 'XENONnT_s2_xy_patterns_GXe_LCE_corrected_qes_MCv4.3.0_wires.pkl?' + '&fmt=pkl' + '&pmt_mask=plugin.pmt_mask', + "electron_trapping_time": config["electron_trapping_time"], + "p_double_pe_emision": config["p_double_pe_emision"], + "pmt_transit_time_spread": config["pmt_transit_time_spread"], + "pmt_transit_time_mean": config["pmt_transit_time_mean"], + "maximum_recombination_time": config["maximum_recombination_time"], + "s1_decay_spread": config["s1_decay_spread"], + "s1_decay_time": config["s1_decay_time"], + "s1_model_type": config["s1_model_type"], + "s1_model_type": config["s1_model_type"], + "n_top_pmts": 253, + "n_tpc_pmts": 494, + "s1_detection_efficiency": 1, + "s1_lce_correction_map": 'itp_map://resource://format://' + 'XENONnT_s1_xyz_LCE_corrected_qes_MCva43fa9b_wires.json.gz?' + '&fmt=json.gz', + "s1_pattern_map":'pattern_map://resource://format://' + 'XENONnT_s1_xyz_patterns_corrected_qes_MCva43fa9b_wires.pkl?' + '&fmt=pkl' + '&pmt_mask=None', + "s1_optical_propagation_spline": 'itp_map://resource://format://' + 'XENONnT_s1_proponly_va43fa9b_wires_20200625.json.gz?' + '&fmt=json.gz' + '&method=RegularGridInterpolator', + "photon_area_distribution": 'simple_load://resource://format://' + f'{config["photon_area_distribution"]}?' + '&fmt=csv', + "s2_gain_spread": 0, + "triplet_lifetime_gas": config["triplet_lifetime_gas"], + "singlet_lifetime_gas": config["singlet_lifetime_gas"], + "triplet_lifetime_liquid": config["triplet_lifetime_liquid"], + "singlet_lifetime_liquid": config["singlet_lifetime_liquid"], + "s2_time_spread": config["s2_time_spread"], + "s2_time_model": config["s2_time_model"], + "singlet_fraction_gas": config["singlet_fraction_gas"], + "s2_luminescence_model": config["s2_luminescence_model"], + "s2_luminescence_model": config["s2_luminescence_model"], + "tpc_radius": config["tpc_radius"], + "diffusion_constant_transverse": config["diffusion_constant_transverse"], + "s2_aft_skewness": config["s2_aft_skewness"], + "s2_aft_sigma": config["s2_aft_sigma"], + "s2_optical_propagation_spline":'itp_map://resource://format://' + 'XENONnT_s2_opticalprop_time_v0.json.gz?' + '&fmt=json.gz', + "s2_luminescence_map":'simple_load://resource://format://' + 'garfield_timing_map_gas_gap_sr0.npy?' + '&fmt=npy', + "garfield_gas_gap_map": 'itp_map://resource://format://' + 'garfield_gas_gap_map_sr0.json?' + '&fmt=json', + "pmt_ap_t_modifier": config["pmt_ap_t_modifier"], + "pmt_ap_modifier": config["pmt_ap_modifier"], + "photon_ap_cdfs": 'simple_load://resource://format://' + f'{config["photon_ap_cdfs"]}?' + '&fmt=json.gz', + "zle_threshold": config["zle_threshold"], + "digitizer_reference_baseline": config["digitizer_reference_baseline"], + "enable_noise": config["enable_noise"], + "high_energy_deamplification_factor": config["high_energy_deamplification_factor"], + "trigger_window": config["trigger_window"], + "external_amplification": config["external_amplification"], + "pmt_pulse_time_rounding": config["pmt_pulse_time_rounding"], + "samples_after_pulse_center": config["samples_after_pulse_center"], + "samples_to_store_after": config["samples_to_store_after"], + "samples_before_pulse_center": config["samples_before_pulse_center"], + "samples_to_store_before": config["samples_to_store_before"], + "dt": config["sample_duration"], + "pe_pulse_ts": config["pe_pulse_ts"], + "pe_pulse_ys": config["pe_pulse_ys"], + "rext": 100000, + "special_thresholds": config["special_thresholds"], + "noise_data_tmp": 'simple_load://resource://format://' + f'{config["noise_file"]}?' + '&fmt=npy', + + }) + + return st + + + +@URLConfig.register('pmt_gains') +def pmt_gains(to_pe, digitizer_voltage_range, digitizer_bits, pmt_circuit_load_resistor): + """Build PMT Gains""" + + to_pe = to_pe[0][1] + + adc_2_current = (digitizer_voltage_range + / 2 ** (digitizer_bits) + / pmt_circuit_load_resistor) + + gains = np.divide(adc_2_current, + to_pe, + out=np.zeros_like(to_pe), + where=to_pe != 0, + ) + return gains + + +@URLConfig.register('pattern_map') +def pattern_map(map_data, pmt_mask, method='WeightedNearestNeighbors'): + """Pattern map handling""" + + if 'compressed' in map_data: + compressor, dtype, shape = map_data['compressed'] + map_data['map'] = np.frombuffer( + strax.io.COMPRESSORS[compressor]['decompress'](map_data['map']), + dtype=dtype).reshape(*shape) + del map_data['compressed'] + if 'quantized' in map_data: + map_data['map'] = map_data['quantized']*map_data['map'].astype(np.float32) + del map_data['quantized'] + if not (pmt_mask is None): + assert (map_data['map'].shape[-1]==pmt_mask.shape[0]), "Error! Pattern map and PMT gains must have same dimensions!" + map_data['map'][..., ~pmt_mask]=0.0 + return straxen.InterpolatingMap(map_data, method=method) + +#Probably not needed! +@URLConfig.register('simple_load') +def load(data): + """Some Documentation""" + return data \ No newline at end of file diff --git a/fuse/plugins/detector_physics/__init__.py b/fuse/plugins/detector_physics/__init__.py index ded59b57..b615ab52 100644 --- a/fuse/plugins/detector_physics/__init__.py +++ b/fuse/plugins/detector_physics/__init__.py @@ -7,11 +7,14 @@ from . import electron_timing from .electron_timing import * -from . import S1_Signal -from .S1_Signal import * +from . import s1_photon_hits +from .s1_photon_hits import * -from . import S2_Signal -from .S2_Signal import * +from . import s1_photon_propagation +from .s1_photon_propagation import * + +from . import s2_photon_propagation +from .s2_photon_propagation import * from . import secondary_scintillation from .secondary_scintillation import * diff --git a/fuse/plugins/detector_physics/csv_input.py b/fuse/plugins/detector_physics/csv_input.py index 193b4617..38bfc314 100644 --- a/fuse/plugins/detector_physics/csv_input.py +++ b/fuse/plugins/detector_physics/csv_input.py @@ -1,4 +1,5 @@ import strax +import straxen import os import numba import logging @@ -6,6 +7,8 @@ import pandas as pd import numpy as np +from ...common import dynamic_chunking + export, __all__ = strax.exporter() logging.basicConfig(handlers=[logging.StreamHandler()]) @@ -13,18 +16,6 @@ log.setLevel('WARNING') @export -@strax.takes_config( - strax.Option('input_file', track=False, infer_type=False, - help="CSV file to read"), - strax.Option('debug', default=False, track=False, infer_type=False, - help="Show debug informations"), - strax.Option('source_rate', default=1, track=False, infer_type=False, - help="source_rate"), - strax.Option('separation_scale', default=1e8, track=False, infer_type=False, - help="Add Description"), - strax.Option('n_interactions_per_chunk', default=25, track=False, infer_type=False, - help="n_interactions_per_chunk"), -) class ChunkCsvInput(strax.Plugin): """ Plugin which reads a CSV file containing instructions for the detector physics simulation @@ -57,6 +48,33 @@ class ChunkCsvInput(strax.Plugin): ] dtype = dtype + strax.time_fields + #Config options + debug = straxen.URLConfig( + default=False, type=bool, track=False, + help='Show debug informations', + ) + + input_file = straxen.URLConfig( + track=False, + infer_type=False, + help='CSV file to read', + ) + + separation_scale = straxen.URLConfig( + default=1e8, type=(int, float), + help='separation_scale', + ) + + source_rate = straxen.URLConfig( + default=1, type=(int, float), + help='source_rate', + ) + + n_interactions_per_chunk = straxen.URLConfig( + default=25, type=(int, float), + help='n_interactions_per_chunk', + ) + def setup(self): if self.debug: @@ -208,32 +226,4 @@ def __load_csv_file(self): for column in df.columns: instructions[column] = df[column] - return instructions, n_simulated_events - - -@numba.njit() -def dynamic_chunking(data, scale, n_min): - - idx_sort = np.argsort(data) - idx_undo_sort = np.argsort(idx_sort) - - data_sorted = data[idx_sort] - - diff = data_sorted[1:] - data_sorted[:-1] - - clusters = np.array([0]) - c = 0 - for value in diff: - if value <= scale: - clusters = np.append(clusters, c) - - elif len(clusters[clusters == c]) < n_min: - clusters = np.append(clusters, c) - - elif value > scale: - c = c + 1 - clusters = np.append(clusters, c) - - clusters_undo_sort = clusters[idx_undo_sort] - - return clusters_undo_sort \ No newline at end of file + return instructions, n_simulated_events \ No newline at end of file diff --git a/fuse/plugins/detector_physics/electron_drift.py b/fuse/plugins/detector_physics/electron_drift.py index 820c535e..2d47a32b 100644 --- a/fuse/plugins/detector_physics/electron_drift.py +++ b/fuse/plugins/detector_physics/electron_drift.py @@ -4,67 +4,20 @@ import os import logging -from ...common import make_map - export, __all__ = strax.exporter() logging.basicConfig(handlers=[logging.StreamHandler()]) log = logging.getLogger('fuse.detector_physics.electron_drift') log.setLevel('WARNING') -#private_files_path = "path/to/private/files" - -base_path = os.path.abspath(os.getcwd()) -private_files_path = os.path.join("/",*base_path.split("/")[:-2], "private_nt_aux_files") - -config = straxen.get_resource(os.path.join(private_files_path, 'sim_files/fax_config_nt_sr0_v4.json') , fmt='json') - @export -@strax.takes_config( - strax.Option('field_distortion_model', default=config["field_distortion_model"], track=False, infer_type=False, - help="field_distortion_model"), - strax.Option('drift_velocity_liquid', default=config["drift_velocity_liquid"], track=False, infer_type=False, - help="drift_velocity_liquid"), - strax.Option('field_distortion_model', default=config["field_distortion_model"], track=False, infer_type=False, - help="field_distortion_model"), - strax.Option('fdc_3d', - default=os.path.join(private_files_path,"sim_files/XnT_3D_FDC_xyz_24_Jun_2022_MC.json.gz"), - track=False, - infer_type=False, - help="fdc_3d map"), - strax.Option('field_distortion_comsol_map', - default=os.path.join(private_files_path,"sim_files/init_to_final_position_mapping_B2d75n_C2d75n_G0d3p_A4d9p_T0d9n_PMTs1d3n_FSR0d65p_QPTFE_0d5n_0d4p.json.gz"), - track=False, - infer_type=False, - help="field_distortion_comsol_map"), - strax.Option('tpc_length', default=config["tpc_length"], track=False, infer_type=False, - help="tpc_length"), - strax.Option('field_dependencies_map', - default=os.path.join(private_files_path,"sim_files/field_dependent_radius_depth_maps_B2d75n_C2d75n_G0d3p_A4d9p_T0d9n_PMTs1d3n_FSR0d65p_QPTFE_0d5n_0d4p.json.gz"), - track=False, - infer_type=False, - help="field_dependencies_map"), - strax.Option('electron_lifetime_liquid', default=config["electron_lifetime_liquid"], track=False, infer_type=False, - help="electron_lifetime_liquid"), - strax.Option('diffusion_longitudinal_map', - default=os.path.join(private_files_path,"sim_files/data_driven_diffusion_map_XENONnTSR0V2.json.gz"), - track=False, - infer_type=False, - help="diffusion_longitudinal_map"), - strax.Option('diffusion_constant_longitudinal', default=config["diffusion_constant_longitudinal"], track=False, infer_type=False, - help="diffusion_constant_longitudinal"), - strax.Option('drift_time_gate', default=config["drift_time_gate"], track=False, infer_type=False, - help="drift_time_gate"), - strax.Option('debug', default=False, track=False, infer_type=False, - help="Show debug informations"), -) class ElectronDrift(strax.Plugin): - __version__ = "0.0.0" + __version__ = "0.0.1" depends_on = ("microphysics_summary") provides = "drifted_electrons" - data_kind = 'electron_cloud' + data_kind = 'interactions_in_roi' dtype = [('n_electron_interface', np.int64), ('drift_time_mean', np.int64), @@ -79,61 +32,105 @@ class ElectronDrift(strax.Plugin): dtype = dtype + strax.time_fields + #Config options + debug = straxen.URLConfig( + default=False, type=bool,track=False, + help='Show debug informations', + ) + + drift_velocity_liquid = straxen.URLConfig( + type=(int, float), + help='drift_velocity_liquid', + ) + + drift_time_gate = straxen.URLConfig( + type=(int, float), + help='drift_time_gate', + ) + + diffusion_constant_longitudinal = straxen.URLConfig( + type=(int, float), + help='diffusion_constant_longitudinal', + ) + + electron_lifetime_liquid = straxen.URLConfig( + type=(int, float), + help='electron_lifetime_liquid', + ) + + enable_field_dependencies = straxen.URLConfig( + help='enable_field_dependencies', + ) + + tpc_length = straxen.URLConfig( + type=(int, float), + help='tpc_length', + ) + + field_distortion_model = straxen.URLConfig( + help='field_distortion_model', + ) + + field_dependencies_map_tmp = straxen.URLConfig( + help='field_dependencies_map', + ) + + diffusion_longitudinal_map_tmp = straxen.URLConfig( + help='diffusion_longitudinal_map', + ) + + fdc_map_fuse = straxen.URLConfig( + cache=True, + help='fdc_map', + ) + def setup(self): if self.debug: log.setLevel('DEBUG') log.debug("Running ElectronDrift in debug mode") + #Can i do this scaling in the url config? if self.field_distortion_model == "inverse_fdc": - self.fdc_3d = make_map(self.fdc_3d, fmt='json.gz') - self.fdc_3d.scale_coordinates([1., 1., - self.drift_velocity_liquid]) - - if self.field_distortion_model == "comsol": - self.fd_comsol = make_map(self.field_distortion_comsol_map, fmt='json.gz', method='RectBivariateSpline') - + self.fdc_map_fuse.scale_coordinates([1., 1., - self.drift_velocity_liquid]) + # Field dependencies - # This config entry a dictionary of 5 items - self.enable_field_dependencies = config['enable_field_dependencies'] #This is not so nice if any(self.enable_field_dependencies.values()): - field_dependencies_map_tmp = make_map(self.field_dependencies_map, fmt='json.gz', method='RectBivariateSpline') self.drift_velocity_scaling = 1.0 # calculating drift velocity scaling to match total drift time for R=0 between cathode and gate if "norm_drift_velocity" in self.enable_field_dependencies.keys(): if self.enable_field_dependencies['norm_drift_velocity']: - norm_dvel = field_dependencies_map_tmp(np.array([ [0], [- self.tpc_length]]).T, map_name='drift_speed_map')[0] + norm_dvel = self.field_dependencies_map_tmp(np.array([ [0], [- self.tpc_length]]).T, map_name='drift_speed_map')[0] norm_dvel*=1e-4 drift_velocity_scaling = self.drift_velocity_liquid/norm_dvel def rz_map(z, xy, **kwargs): r = np.sqrt(xy[:, 0]**2 + xy[:, 1]**2) - return field_dependencies_map_tmp(np.array([r, z]).T, **kwargs) + return self.field_dependencies_map_tmp(np.array([r, z]).T, **kwargs) self.field_dependencies_map = rz_map # Data-driven longitudinal diffusion map # TODO: Change to the best way to accommodate simulation/data-driven map if self.enable_field_dependencies["diffusion_longitudinal_map"]: - diffusion_longitudinal_map_tmp = make_map(self.diffusion_longitudinal_map, fmt='json.gz', - method='WeightedNearestNeighbors') def _rz_map(z, xy, **kwargs): r = np.sqrt(xy[:, 0]**2 + xy[:, 1]**2) - return diffusion_longitudinal_map_tmp(np.array([r, z]).T, **kwargs) + return self.diffusion_longitudinal_map_tmp(np.array([r, z]).T, **kwargs) self.diffusion_longitudinal_map = _rz_map - def compute(self, clustered_interactions): + def compute(self, interactions_in_roi): - #Just apply this to clusters with free electrons - instruction = clustered_interactions[clustered_interactions["electrons"] > 0] + #Just apply this to clusters with photons + mask = interactions_in_roi["electrons"] > 0 - if len(instruction) == 0: + if len(interactions_in_roi[mask]) == 0: return np.zeros(0, self.dtype) - t = instruction["time"] - x = instruction["x"] - y = instruction["y"] - z = instruction["z"] - n_electron = instruction["electrons"].astype(np.int64) - recoil_type = instruction["nestid"] + t = interactions_in_roi[mask]["time"] + x = interactions_in_roi[mask]["x"] + y = interactions_in_roi[mask]["y"] + z = interactions_in_roi[mask]["z"] + n_electron = interactions_in_roi[mask]["electrons"].astype(np.int64) + recoil_type = interactions_in_roi[mask]["nestid"] recoil_type = np.where(np.isin(recoil_type, [0, 6, 7, 8, 11]), recoil_type, 8) # Reverse engineering FDC @@ -160,17 +157,17 @@ def compute(self, clustered_interactions): n_electron = n_electron*electron_lifetime_correction - result = np.zeros(len(n_electron), dtype = self.dtype) - result["time"] = instruction["time"] - result["endtime"] = instruction["endtime"] - result["n_electron_interface"] = n_electron - result["drift_time_mean"] = drift_time_mean - result["drift_time_spread"] = drift_time_spread + result = np.zeros(len(interactions_in_roi), dtype = self.dtype) + result["time"] = interactions_in_roi["time"] + result["endtime"] = interactions_in_roi["endtime"] + result["n_electron_interface"][mask] = n_electron + result["drift_time_mean"][mask] = drift_time_mean + result["drift_time_spread"][mask] = drift_time_spread #These ones are needed later - result["x"] = positions.T[0] - result["y"] = positions.T[1] - result["z_obs"] = z_obs + result["x"][mask] = positions.T[0] + result["y"][mask] = positions.T[1] + result["z_obs"][mask] = z_obs return result @@ -187,7 +184,7 @@ def inverse_field_distortion_correction(self, x, y, z): """ positions = np.array([x, y, z]).T for i_iter in range(6): # 6 iterations seems to work - dr = self.fdc_3d(positions) + dr = self.fdc_map_fuse(positions) if i_iter > 0: dr = 0.5 * dr + 0.5 * dr_pre # Average between iter dr_pre = dr @@ -211,7 +208,7 @@ def field_distortion_comsol(self, x, y, z): """ positions = np.array([np.sqrt(x**2 + y**2), z]).T theta = np.arctan2(y, x) - r_obs = self.fd_comsol(positions, map_name='r_distortion_map') + r_obs = self.fdc_map_fuse(positions, map_name='r_distortion_map') x_obs = r_obs * np.cos(theta) y_obs = r_obs * np.sin(theta) @@ -247,7 +244,7 @@ def get_s2_drift_time_params(self, xy_int, z_int): drift_time_mean = - z_int / \ drift_velocity_liquid + self.drift_time_gate drift_time_mean = np.clip(drift_time_mean, 0, np.inf) - drift_time_spread = np.sqrt(2 * self.diffusion_constant_longitudinal * drift_time_mean) + 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 @@ -263,6 +260,4 @@ def get_avg_drift_velocity(self, z, xy): drift_v_LXe *= self.drift_velocity_scaling else: drift_v_LXe = self.drift_velocity_liquid - return drift_v_LXe - - + return drift_v_LXe \ No newline at end of file diff --git a/fuse/plugins/detector_physics/electron_extraction.py b/fuse/plugins/detector_physics/electron_extraction.py index e4614702..d32b4713 100644 --- a/fuse/plugins/detector_physics/electron_extraction.py +++ b/fuse/plugins/detector_physics/electron_extraction.py @@ -1,69 +1,23 @@ import strax import straxen import numpy as np -from copy import deepcopy import os import logging -from ...common import make_map, make_patternmap - export, __all__ = strax.exporter() logging.basicConfig(handlers=[logging.StreamHandler()]) log = logging.getLogger('fuse.detector_physics.electron_extraction') log.setLevel('WARNING') - -#private_files_path = "path/to/private/files" -base_path = os.path.abspath(os.getcwd()) -private_files_path = os.path.join("/",*base_path.split("/")[:-2], "private_nt_aux_files") -config = straxen.get_resource(os.path.join(private_files_path, 'sim_files/fax_config_nt_sr0_v4.json') , fmt='json') - @export -@strax.takes_config( - strax.Option('s2_correction_map_file', - default=os.path.join(private_files_path, "strax_files/XENONnT_s2_xy_map_v4_210503_mlp_3_in_1_iterated.json"), - track=False, - infer_type=False, - help="s2_correction_map"), - strax.Option('s2_pattern_map_file', - default=os.path.join(private_files_path, "sim_files/XENONnT_s2_xy_patterns_GXe_LCE_corrected_qes_MCv4.3.0_wires.pkl"), - track=False, - infer_type=False, - help="s2_pattern_map"), - strax.Option('to_pe_file', default=os.path.join(private_files_path, "sim_files/to_pe_nt.npy"), track=False, infer_type=False, - help="to_pe file"), - strax.Option('digitizer_voltage_range', default=config['digitizer_voltage_range'], track=False, infer_type=False, - help="digitizer_voltage_range"), - strax.Option('digitizer_bits', default=config['digitizer_bits'], track=False, infer_type=False, - help="digitizer_bits"), - strax.Option('pmt_circuit_load_resistor', default=config['pmt_circuit_load_resistor'], track=False, infer_type=False, - help="pmt_circuit_load_resistor"), - strax.Option('se_gain_from_map', default=config['se_gain_from_map'], track=False, infer_type=False, - help="se_gain_from_map"), - strax.Option('se_gain_map', - default=os.path.join(private_files_path, "strax_files/XENONnT_se_xy_map_v1_mlp.json"), - track=False, - infer_type=False, - help="se_gain_map"), - strax.Option('s2_secondary_sc_gain', default=config['s2_secondary_sc_gain'], track=False, infer_type=False, - help="s2_secondary_sc_gain"), - strax.Option('g2_mean', default=config['g2_mean'], track=False, infer_type=False, - help="g2_mean"), - strax.Option('electron_extraction_yield', default=config['electron_extraction_yield'], track=False, infer_type=False, - help="electron_extraction_yield"), - strax.Option('ext_eff_from_map', default=config['ext_eff_from_map'], track=False, infer_type=False, - help="ext_eff_from_map"), - strax.Option('debug', default=False, track=False, infer_type=False, - help="Show debug informations"), -) class ElectronExtraction(strax.Plugin): __version__ = "0.0.0" depends_on = ("microphysics_summary", "drifted_electrons") provides = "extracted_electrons" - data_kind = "electron_cloud" + data_kind = "interactions_in_roi" #Forbid rechunking rechunk_on_save = False @@ -73,6 +27,73 @@ class ElectronExtraction(strax.Plugin): ] dtype = dtype + strax.time_fields + + #Config options + debug = straxen.URLConfig( + default=False, type=bool,track=False, + help='Show debug informations', + ) + + digitizer_voltage_range = straxen.URLConfig( + type=(int, float), + help='digitizer_voltage_range', + ) + + digitizer_bits = straxen.URLConfig( + type=(int, float), + help='digitizer_bits', + ) + + pmt_circuit_load_resistor = straxen.URLConfig( + type=(int, float), + help='pmt_circuit_load_resistor', + ) + + s2_secondary_sc_gain = straxen.URLConfig( + type=(int, float), + help='s2_secondary_sc_gain', + ) + #Rename? -> g2_value in beta_yields model + g2_mean = straxen.URLConfig( + type=(int, float), + help='g2_mean', + ) + + electron_extraction_yield = straxen.URLConfig( + type=(int, float), + help='electron_extraction_yield', + ) + + ext_eff_from_map = straxen.URLConfig( + type=bool, + help='ext_eff_from_map', + ) + + se_gain_from_map = straxen.URLConfig( + type=bool, + help='se_gain_from_map', + ) + + gains = straxen.URLConfig( + cache=True, + help='pmt gains', + ) + + s2_correction_map = straxen.URLConfig( + cache=True, + help='s2_correction_map', + ) + + se_gain_map = straxen.URLConfig( + cache=True, + help='se_gain_map', + ) + + s2_pattern_map = straxen.URLConfig( + cache=True, + help='s2_pattern_map', + ) + def setup(self): @@ -80,46 +101,31 @@ def setup(self): log.setLevel('DEBUG') log.debug("Running ElectronExtraction in debug mode") - to_pe = straxen.get_resource(self.to_pe_file, fmt='npy') - self.to_pe = to_pe[0][1] + self.pmt_mask = np.array(self.gains) > 0 # Converted from to pe (from cmt by default) - adc_2_current = (self.digitizer_voltage_range - / 2 ** (self.digitizer_bits) - / self.pmt_circuit_load_resistor) - - gains = np.divide(adc_2_current, - self.to_pe, - out=np.zeros_like(self.to_pe), - where=self.to_pe != 0) - - self.pmt_mask = np.array(gains) > 0 # Converted from to pe (from cmt by default) - - self.s2_pattern_map = make_patternmap(self.s2_pattern_map_file, fmt='pkl', pmt_mask=self.pmt_mask) - - if self.s2_correction_map_file: - self.s2_correction_map = make_map(self.s2_correction_map_file, 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 - s2cmap.data['map'] = np.sum(s2cmap.data['map'][:][:], axis=2, keepdims=True, where=self.pmt_mask) - # Scale by median value - s2cmap.data['map'] = s2cmap.data['map'] / np.median(s2cmap.data['map'][s2cmap.data['map'] > 0]) - s2cmap.__init__(s2cmap.data) - self.s2_correction_map = s2cmap - - self.se_gain_map = make_map(self.se_gain_map, fmt = "json") + #Is this else case ever used? if no -> remove + #if self.s2_correction_map_file: + # self.s2_correction_map = make_map(self.s2_correction_map_file, 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 + # s2cmap.data['map'] = np.sum(s2cmap.data['map'][:][:], axis=2, keepdims=True, where=self.pmt_mask) + # # Scale by median value + # s2cmap.data['map'] = s2cmap.data['map'] / np.median(s2cmap.data['map'][s2cmap.data['map'] > 0]) + # s2cmap.__init__(s2cmap.data) + # self.s2_correction_map = s2cmap - def compute(self, clustered_interactions, electron_cloud): + def compute(self, interactions_in_roi): - #Just apply this to clusters with free electrons - instruction = clustered_interactions[clustered_interactions["electrons"] > 0] + #Just apply this to clusters with photons + mask = interactions_in_roi["electrons"] > 0 - if len(instruction) == 0: + if len(interactions_in_roi[mask]) == 0: return np.zeros(0, self.dtype) - x = instruction["x"] - y = instruction["y"] + x = interactions_in_roi[mask]["x"] + y = interactions_in_roi[mask]["y"] xy_int = np.array([x, y]).T # maps are in R_true, so orginal position should be here @@ -140,11 +146,11 @@ def compute(self, clustered_interactions, electron_cloud): else: cy = self.electron_extraction_yield - n_electron = np.random.binomial(n=electron_cloud["n_electron_interface"], p=cy) + n_electron = np.random.binomial(n=interactions_in_roi[mask]["n_electron_interface"], p=cy) - result = np.zeros(len(n_electron), dtype=self.dtype) - result["n_electron_extracted"] = n_electron - result["time"] = instruction["time"] - result["endtime"] = instruction["endtime"] + result = np.zeros(len(interactions_in_roi), dtype=self.dtype) + result["n_electron_extracted"][mask] = n_electron + result["time"] = interactions_in_roi["time"] + result["endtime"] = interactions_in_roi["endtime"] return result \ No newline at end of file diff --git a/fuse/plugins/detector_physics/electron_timing.py b/fuse/plugins/detector_physics/electron_timing.py index 6d033f2b..e5562363 100644 --- a/fuse/plugins/detector_physics/electron_timing.py +++ b/fuse/plugins/detector_physics/electron_timing.py @@ -10,24 +10,14 @@ log = logging.getLogger('fuse.detector_physics.electron_timing') log.setLevel('WARNING') -#private_files_path = "path/to/private/files" -base_path = os.path.abspath(os.getcwd()) -private_files_path = os.path.join("/",*base_path.split("/")[:-2], "private_nt_aux_files") -config = straxen.get_resource(os.path.join(private_files_path, 'sim_files/fax_config_nt_sr0_v4.json') , fmt='json') - @export -@strax.takes_config( - strax.Option('electron_trapping_time', default=config["electron_trapping_time"], track=False, infer_type=False, - help="electron_trapping_time"), - strax.Option('debug', default=False, track=False, infer_type=False, - help="Show debug informations"), -) class ElectronTiming(strax.Plugin): __version__ = "0.0.0" depends_on = ("drifted_electrons", "extracted_electrons") provides = "electron_time" + data_kind = "individual_electrons" #Forbid rechunking rechunk_on_save = False @@ -39,7 +29,17 @@ class ElectronTiming(strax.Plugin): ('y', np.float64), ] dtype = dtype + strax.time_fields - #dtype = strax.time_fields + + #Config options + debug = straxen.URLConfig( + default=False, type=bool,track=False, + help='Show debug informations', + ) + + electron_trapping_time = straxen.URLConfig( + type=(int, float), + help='electron_trapping_time', + ) def setup(self): @@ -47,20 +47,23 @@ def setup(self): log.setLevel('DEBUG') log.debug("Running ElectronTiming in debug mode") - def compute(self, electron_cloud): + def compute(self, interactions_in_roi): + + #Just apply this to clusters with photons + mask = interactions_in_roi["n_electron_extracted"] > 0 - if len(electron_cloud) == 0: + if len(interactions_in_roi[mask]) == 0: return np.zeros(0, dtype=self.dtype) - timing = self.electron_timing(electron_cloud["time"], - electron_cloud["n_electron_extracted"], - electron_cloud["drift_time_mean"], - electron_cloud["drift_time_spread"] + timing = self.electron_timing(interactions_in_roi[mask]["time"], + interactions_in_roi[mask]["n_electron_extracted"], + interactions_in_roi[mask]["drift_time_mean"], + interactions_in_roi[mask]["drift_time_spread"] ) - x = np.repeat(electron_cloud["x"], electron_cloud["n_electron_extracted"]) - y = np.repeat(electron_cloud["y"], electron_cloud["n_electron_extracted"]) + x = np.repeat(interactions_in_roi[mask]["x"], interactions_in_roi[mask]["n_electron_extracted"]) + y = np.repeat(interactions_in_roi[mask]["y"], interactions_in_roi[mask]["n_electron_extracted"]) result = np.zeros(len(timing), dtype = self.dtype) result["time"] = timing diff --git a/fuse/plugins/detector_physics/s1_photon_hits.py b/fuse/plugins/detector_physics/s1_photon_hits.py new file mode 100644 index 00000000..857649ef --- /dev/null +++ b/fuse/plugins/detector_physics/s1_photon_hits.py @@ -0,0 +1,104 @@ +import strax +import straxen +import logging + +import numpy as np + +export, __all__ = strax.exporter() + +logging.basicConfig(handlers=[logging.StreamHandler()]) +log = logging.getLogger('fuse.detector_physics.s1_photon_hits') +log.setLevel('WARNING') + +@export +class S1PhotonHits(strax.Plugin): + + __version__ = '0.0.0' + + depends_on = ("microphysics_summary") + provides = "s1_photons" + data_kind = "interactions_in_roi" + + #Forbid rechunking + rechunk_on_save = False + + dtype = [('n_s1_photon_hits', np.int64), + ] + dtype = dtype + strax.time_fields + + #Config options + debug = straxen.URLConfig( + default=False, type=bool,track=False, + help='Show debug informations', + ) + + s1_lce_correction_map = straxen.URLConfig( + cache=True, + help='s1_lce_correction_map', + ) + + p_double_pe_emision = straxen.URLConfig( + type=(int, float), + help='p_double_pe_emision', + ) + + s1_detection_efficiency = straxen.URLConfig( + type=(int, float), + help='Some placeholder for s1_detection_efficiency', + ) + + def setup(self): + + if self.debug: + log.setLevel('DEBUG') + log.debug("Running S1PhotonHits in debug mode") + + def compute(self, interactions_in_roi): + + #Just apply this to clusters with photons + mask = interactions_in_roi["photons"] > 0 + + if len(interactions_in_roi[mask]) == 0: + return np.zeros(0, self.dtype) + + x = interactions_in_roi[mask]['x'] + y = interactions_in_roi[mask]['y'] + z = interactions_in_roi[mask]['z'] + n_photons = interactions_in_roi[mask]['photons'].astype(np.int64) + + positions = np.array([x, y, z]).T + + n_photon_hits = self.get_n_photons(n_photons=n_photons, + positions=positions, + ) + + result = np.zeros(interactions_in_roi.shape[0], dtype = self.dtype) + + result["time"] = interactions_in_roi["time"] + result["endtime"] = interactions_in_roi["endtime"] + + result["n_s1_photon_hits"][mask] = n_photon_hits + + return result + + + def get_n_photons(self, n_photons, positions): + + """Calculates number of detected photons based on number of photons in total and the positions + :param n_photons: 1d array of ints with number of emitted S1 photons: + :param positions: 2d array with xyz positions of interactions + :param s1_lce_correction_map: interpolator instance of s1 light yield map + :param config: dict wfsim config + + return array with number photons""" + ly = self.s1_lce_correction_map(positions) + # depending on if you use the data driven or mc pattern map for light yield + #the shape of n_photon_hits will change. Mc needs a squeeze + if len(ly.shape) != 1: + ly = np.squeeze(ly, axis=-1) + ly /= 1 + self.p_double_pe_emision + ly *= self.s1_detection_efficiency + + n_photon_hits = np.random.binomial(n=n_photons, p=ly) + + return n_photon_hits \ No newline at end of file diff --git a/fuse/plugins/detector_physics/S1_Signal.py b/fuse/plugins/detector_physics/s1_photon_propagation.py similarity index 64% rename from fuse/plugins/detector_physics/S1_Signal.py rename to fuse/plugins/detector_physics/s1_photon_propagation.py index 7759f575..27dddfc8 100644 --- a/fuse/plugins/detector_physics/S1_Signal.py +++ b/fuse/plugins/detector_physics/s1_photon_propagation.py @@ -2,85 +2,25 @@ import strax import straxen import nestpy -import os import logging -from numba import njit from strax import deterministic_hash from scipy.interpolate import interp1d -export, __all__ = strax.exporter() +from ...common import loop_uniform_to_pe_arr -from ...common import make_map, make_patternmap +export, __all__ = strax.exporter() logging.basicConfig(handlers=[logging.StreamHandler()]) log = logging.getLogger('fuse.detector_physics.S1_Signal') log.setLevel('WARNING') -#private_files_path = "path/to/private/files" -base_path = os.path.abspath(os.getcwd()) -private_files_path = os.path.join("/",*base_path.split("/")[:-2], "private_nt_aux_files") -config = straxen.get_resource(os.path.join(private_files_path, 'sim_files/fax_config_nt_sr0_v4.json') , fmt='json') - @export -@strax.takes_config( - strax.Option('s1_detection_efficiency', default=1, track=False, infer_type=False, - help="Some placeholder for s1_detection_efficiency"), - strax.Option('p_double_pe_emision', default=config["p_double_pe_emision"], track=False, infer_type=False, - help="Some placeholder for p_double_pe_emision"), - strax.Option('s1_lce_correction_map', - default=os.path.join(private_files_path, "sim_files/XENONnT_s1_xyz_LCE_corrected_qes_MCva43fa9b_wires.json.gz"), - track=False, - infer_type=False, - help="S1 LCE correction map"), - strax.Option('s1_pattern_map', - default=os.path.join(private_files_path, "sim_files/XENONnT_s1_xyz_patterns_corrected_qes_MCva43fa9b_wires.pkl"), - track=False, - infer_type=False, - help="S1 pattern map"), - strax.Option('n_tpc_pmts', default=494, track=False, infer_type=False, - help="Number of PMTs in the TPC"), - strax.Option('n_top_pmts', default=253, track=False, infer_type=False, - help="Number of PMTs on top array"), - strax.Option('digitizer_voltage_range', default=config['digitizer_voltage_range'], track=False, infer_type=False, - help="digitizer_voltage_range"), - strax.Option('digitizer_bits', default=config['digitizer_bits'], track=False, infer_type=False, - help="digitizer_bits"), - strax.Option('pmt_circuit_load_resistor', default=config['pmt_circuit_load_resistor'], track=False, infer_type=False, - help="pmt_circuit_load_resistor"), - strax.Option('to_pe_file', default=os.path.join(private_files_path, "sim_files/to_pe_nt.npy"), track=False, infer_type=False, - help="to_pe file"), - strax.Option('s1_model_type', default=config['s1_model_type'], track=True, infer_type=False, - help="s1_model_type"), - strax.Option('s1_time_spline', - default=os.path.join(private_files_path, "sim_files/XENONnT_s1_proponly_va43fa9b_wires_20200625.json.gz"), - track=False, - infer_type=False, - help="S1 Time Spline"), - strax.Option('s1_decay_time', default=config['s1_decay_time'], track=False, infer_type=False, - help="s1_decay_time"), - strax.Option('s1_decay_spread', default=config['s1_decay_spread'], track=False, infer_type=False, - help="s1_decay_spread"), - strax.Option('phase', default="liquid", track=False, infer_type=False, - help="xenon phase"), - strax.Option('maximum_recombination_time', default=config["maximum_recombination_time"], track=False, infer_type=False, - help="maximum_recombination_time"), - strax.Option('pmt_transit_time_mean', default=config['pmt_transit_time_mean'], track=False, infer_type=False, - help="pmt_transit_time_mean"), - strax.Option('pmt_transit_time_spread', default=config['pmt_transit_time_spread'], track=False, infer_type=False, - help="pmt_transit_time_spread"), - strax.Option('p_double_pe_emision', default=config['p_double_pe_emision'], track=False, infer_type=False, - help="p_double_pe_emision"), - strax.Option('photon_area_distribution', default=config['photon_area_distribution'], track=False, infer_type=False, - help="photon_area_distribution"), - strax.Option('debug', default=False, track=False, infer_type=False, - help="Show debug informations"), -) class S1PhotonPropagation(strax.Plugin): __version__ = "0.0.0" - depends_on = ("microphysics_summary") + depends_on = ("s1_photons", "microphysics_summary") provides = "propagated_s1_photons" data_kind = "S1_photons" @@ -92,48 +32,116 @@ class S1PhotonPropagation(strax.Plugin): ('photon_gain', np.int64), ] dtype = dtype + strax.time_fields + + #Config options + debug = straxen.URLConfig( + default=False, type=bool,track=False, + help='Show debug informations', + ) + + p_double_pe_emision = straxen.URLConfig( + type=(int, float), + help='p_double_pe_emision', + ) + + pmt_transit_time_spread = straxen.URLConfig( + type=(int, float), + help='pmt_transit_time_spread', + ) + + pmt_transit_time_mean = straxen.URLConfig( + type=(int, float), + help='pmt_transit_time_mean', + ) + + maximum_recombination_time = straxen.URLConfig( + type=(int, float), + help='maximum_recombination_time', + ) + + phase = straxen.URLConfig( + default="liquid", + help='phase', + ) + + s1_decay_spread = straxen.URLConfig( + type=(int, float), + help='s1_decay_spread', + ) + + s1_decay_time = straxen.URLConfig( + type=(int, float), + help='s1_decay_time', + ) + + s1_model_type = straxen.URLConfig( + help='s1_model_type', + ) + + pmt_circuit_load_resistor = straxen.URLConfig( + help='pmt_circuit_load_resistor', type=(int, float), + ) + + digitizer_bits = straxen.URLConfig( + type=(int, float), + help='digitizer_bits', + ) + + digitizer_voltage_range = straxen.URLConfig( + type=(int, float), + help='digitizer_voltage_range', + ) + + n_top_pmts = straxen.URLConfig( + type=(int), + help='Number of PMTs on top array', + ) + + n_tpc_pmts = straxen.URLConfig( + type=(int), + help='Number of PMTs in the TPC', + ) + + s1_pattern_map = straxen.URLConfig( + cache=True, + help='s1_pattern_map', + ) + + gains = straxen.URLConfig( + cache=True, + help='pmt gains', + ) + + s1_optical_propagation_spline = straxen.URLConfig( + cache=True, + help='s1_optical_propagation_spline', + ) + + photon_area_distribution = straxen.URLConfig( + cache=True, + help='photon_area_distribution', + ) def setup(self): if self.debug: log.setLevel('DEBUG') log.debug("Running S1PhotonPropagation in debug mode") - - self.s1_lce_correction_map = make_map(self.s1_lce_correction_map, fmt='json.gz') - self.s1_pattern_map = make_patternmap(self.s1_pattern_map, fmt='pkl', pmt_mask=None) - - to_pe = straxen.get_resource(self.to_pe_file, fmt='npy') - self.to_pe = to_pe[0][1] - - adc_2_current = (self.digitizer_voltage_range - / 2 ** (self.digitizer_bits) - / self.pmt_circuit_load_resistor) - - self.gains = np.divide(adc_2_current, - self.to_pe, - out=np.zeros_like(self.to_pe), - where=self.to_pe != 0) self.turned_off_pmts = np.arange(len(self.gains))[np.array(self.gains) == 0] - - self.s1_optical_propagation_spline = make_map(self.s1_time_spline, - fmt='json.gz', - method='RegularGridInterpolator') - if 'nest' in self.s1_model_type: #and (self.nestpy_calc is None): log.info('Using NEST for scintillation time without set calculator\n' 'Creating new nestpy calculator') self.nestpy_calc = nestpy.NESTcalc(nestpy.DetectorExample_XENON10()) - self.photon_area_distribution = straxen.get_resource(self.photon_area_distribution, fmt='csv') self._cached_uniform_to_pe_arr = {} self.__uniform_to_pe_arr = self.init_spe_scaling_factor_distributions() - def compute(self, clustered_interactions): + def compute(self, interactions_in_roi): - #Just apply this to clusters with free electrons - instruction = clustered_interactions[clustered_interactions["photons"] > 0] + #Just apply this to clusters with photons hitting a PMT + instruction = interactions_in_roi[interactions_in_roi["n_s1_photon_hits"] > 0] if len(instruction) == 0: return np.zeros(0, self.dtype) @@ -146,13 +154,9 @@ def compute(self, clustered_interactions): recoil_type = instruction['nestid'] positions = np.array([x, y, z]).T # For map interpolation - n_photon_hits = self.get_n_photons(n_photons=n_photons, - positions=positions, - ) - # The new way interpolation is written always require a list _photon_channels = self.photon_channels(positions=positions, - n_photon_hits=n_photon_hits, + n_photon_hits=instruction["n_s1_photon_hits"], ) extra_targs = {} @@ -164,7 +168,7 @@ def compute(self, clustered_interactions): extra_targs['nestpy_calc'] = self.nestpy_calc _photon_timings = self.photon_timings(t=t, - n_photon_hits=n_photon_hits, + n_photon_hits=instruction["n_s1_photon_hits"], recoil_type=recoil_type, channels=_photon_channels, positions=positions, @@ -183,6 +187,7 @@ def compute(self, clustered_interactions): self.pmt_transit_time_spread / 2.35482, len(_photon_timings)).astype(np.int64) + #Why is this done here and additionally in the get_n_photons function of S1PhotonHits?? _photon_is_dpe = np.random.binomial(n=1, p=self.p_double_pe_emision, size=len(_photon_timings)).astype(np.bool_) @@ -206,27 +211,6 @@ def compute(self, clustered_interactions): return result - def get_n_photons(self, n_photons, positions): - - """Calculates number of detected photons based on number of photons in total and the positions - :param n_photons: 1d array of ints with number of emitted S1 photons: - :param positions: 2d array with xyz positions of interactions - :param s1_lce_correction_map: interpolator instance of s1 light yield map - :param config: dict wfsim config - - return array with number photons""" - ly = self.s1_lce_correction_map(positions) - # depending on if you use the data driven or mc pattern map for light yield - #the shape of n_photon_hits will change. Mc needs a squeeze - if len(ly.shape) != 1: - ly = np.squeeze(ly, axis=-1) - ly /= 1 + self.p_double_pe_emision - ly *= self.s1_detection_efficiency - - n_photon_hits = np.random.binomial(n=n_photons, p=ly) - - return n_photon_hits - def photon_channels(self, positions, n_photon_hits): """Calculate photon arrival channels @@ -361,7 +345,7 @@ def optical_propagation(self, channels, z_positions): def init_spe_scaling_factor_distributions(self): #This code will be duplicate with the corresponding S2 class # Improve!! - + config="PLACEHOLDER_FIX_THIS_PART" h = deterministic_hash(config) # What is this part doing? #if h in self._cached_uniform_to_pe_arr: # __uniform_to_pe_arr = self._cached_uniform_to_pe_arr[h] @@ -396,22 +380,4 @@ def init_spe_scaling_factor_distributions(self): self._cached_uniform_to_pe_arr[h] = __uniform_to_pe_arr log.debug('Spe scaling factors created, cached with key %s' % h) - return __uniform_to_pe_arr - - -#This is a modified version of the corresponding WFsim code.... -@njit() -def uniform_to_pe_arr(p, channel, __uniform_to_pe_arr): - indices = np.int64(p * 2000) + 1 - return __uniform_to_pe_arr[channel, indices] - -#In WFSim uniform_to_pe_arr is called inside a loop over the channels -#I needed to change the code to run on all channels at once -@njit() -def loop_uniform_to_pe_arr(p, channel, __uniform_to_pe_arr): - result = [] - for i in range(len(p)): - result.append(uniform_to_pe_arr(p[i], - channel=channel[i], - __uniform_to_pe_arr=__uniform_to_pe_arr) ) - return np.array(result) \ No newline at end of file + return __uniform_to_pe_arr \ No newline at end of file diff --git a/fuse/plugins/detector_physics/S2_Signal.py b/fuse/plugins/detector_physics/s2_photon_propagation.py similarity index 68% rename from fuse/plugins/detector_physics/S2_Signal.py rename to fuse/plugins/detector_physics/s2_photon_propagation.py index dfe79603..70d802c9 100644 --- a/fuse/plugins/detector_physics/S2_Signal.py +++ b/fuse/plugins/detector_physics/s2_photon_propagation.py @@ -11,84 +11,14 @@ from strax import deterministic_hash export, __all__ = strax.exporter() -from ...common import make_map, make_patternmap, DummyMap +from ...common import DummyMap, loop_uniform_to_pe_arr logging.basicConfig(handlers=[logging.StreamHandler()]) log = logging.getLogger('fuse.detector_physics.S2_Signal') log.setLevel('WARNING') -#private_files_path = "path/to/private/files" -base_path = os.path.abspath(os.getcwd()) -private_files_path = os.path.join("/",*base_path.split("/")[:-2], "private_nt_aux_files") -config = straxen.get_resource(os.path.join(private_files_path, 'sim_files/fax_config_nt_sr0_v4.json') , fmt='json') @export -@strax.takes_config( - strax.Option('s2_aft_sigma', default=config["s2_aft_sigma"], track=False, infer_type=False, - help="s2_aft_sigma"), - strax.Option('s2_aft_skewness', default=config["s2_aft_skewness"], track=False, infer_type=False, - help="s2_aft_skewness"), - strax.Option('diffusion_constant_transverse', default=config["diffusion_constant_transverse"], track=False, infer_type=False, - help="diffusion_constant_transverse"), - strax.Option('n_top_pmts', default=253, track=False, infer_type=False, - help="n_top_pmts"), - strax.Option('n_tpc_pmts', default=494, track=False, infer_type=False, - help="n_tpc_pmts"), - strax.Option('tpc_radius', default=config["tpc_radius"], track=False, infer_type=False, - help="tpc_radius"), - strax.Option('to_pe_file', default=os.path.join(private_files_path,"sim_files/to_pe_nt.npy"), track=False, infer_type=False, - help="to_pe file"), - strax.Option('digitizer_voltage_range', default=config['digitizer_voltage_range'], track=False, infer_type=False, - help="digitizer_voltage_range"), - strax.Option('digitizer_bits', default=config['digitizer_bits'], track=False, infer_type=False, - help="digitizer_bits"), - strax.Option('pmt_circuit_load_resistor', default=config['pmt_circuit_load_resistor'], track=False, infer_type=False, - help="pmt_circuit_load_resistor"), - strax.Option('s2_pattern_map_file', - default=os.path.join(private_files_path,"sim_files/XENONnT_s2_xy_patterns_GXe_LCE_corrected_qes_MCv4.3.0_wires.pkl"), - track=False, - infer_type=False, - help="s2_pattern_map"), - strax.Option('field_dependencies_map', - default=os.path.join(private_files_path,"sim_files/field_dependent_radius_depth_maps_B2d75n_C2d75n_G0d3p_A4d9p_T0d9n_PMTs1d3n_FSR0d65p_QPTFE_0d5n_0d4p.json.gz"), - track=False, - infer_type=False, - help="field_dependencies_map"), - strax.Option('tpc_length', default=config['tpc_length'], track=False, infer_type=False, - help="tpc_length"), - strax.Option('drift_velocity_liquid', default=config['drift_velocity_liquid'], track=False, infer_type=False, - help="drift_velocity_liquid"), - strax.Option('s2_luminescence_model', default=config['s2_luminescence_model'], track=False, infer_type=False, - help="s2_luminescence_model"), - strax.Option('phase_s2', default="gas", track=False, infer_type=False, - help="phase"), - strax.Option('singlet_fraction_gas', default=config['singlet_fraction_gas'], track=False, infer_type=False, - help="singlet_fraction_gas"), - strax.Option('s2_time_model', default=config['s2_time_model'], track=False, infer_type=False, - help="s2_time_model"), - strax.Option('s2_time_spline', default=config['s2_time_spline'], track=False, infer_type=False, - help="s2_time_spline"), - strax.Option('s2_time_spread', default=config['s2_time_spread'], track=False, infer_type=False, - help="s2_time_spread"), - strax.Option('singlet_lifetime_liquid', default=config['singlet_lifetime_liquid'], track=False, infer_type=False, - help="singlet_lifetime_liquid"), - strax.Option('triplet_lifetime_liquid', default=config['triplet_lifetime_liquid'], track=False, infer_type=False, - help="triplet_lifetime_liquid"), - strax.Option('singlet_lifetime_gas', default=config['singlet_lifetime_gas'], track=False, infer_type=False, - help="singlet_lifetime_gas"), - strax.Option('triplet_lifetime_gas', default=config['triplet_lifetime_gas'], track=False, infer_type=False, - help="triplet_lifetime_gas"), - strax.Option('pmt_transit_time_mean', default=config['pmt_transit_time_mean'], track=False, infer_type=False, - help="pmt_transit_time_mean"), - strax.Option('pmt_transit_time_spread', default=config['pmt_transit_time_spread'], track=False, infer_type=False, - help="pmt_transit_time_spread"), - strax.Option('p_double_pe_emision', default=config['p_double_pe_emision'], track=False, infer_type=False, - help="p_double_pe_emision"), - strax.Option('photon_area_distribution', default=config['photon_area_distribution'], track=False, infer_type=False, - help="photon_area_distribution"), - strax.Option('debug', default=False, track=False, infer_type=False, - help="Show debug informations"), -) class S2PhotonPropagation(strax.Plugin): __version__ = "0.0.0" @@ -105,84 +35,217 @@ class S2PhotonPropagation(strax.Plugin): ('photon_gain', np.int64), ] dtype = dtype + strax.time_fields + + #Config options + debug = straxen.URLConfig( + default=False, type=bool,track=False, + help='Show debug informations', + ) + + p_double_pe_emision = straxen.URLConfig( + type=(int, float), + help='p_double_pe_emision', + ) + + pmt_transit_time_spread = straxen.URLConfig( + type=(int, float), + help='pmt_transit_time_spread', + ) + + pmt_transit_time_mean = straxen.URLConfig( + type=(int, float), + help='pmt_transit_time_mean', + ) + + triplet_lifetime_gas = straxen.URLConfig( + type=(int, float), + help='triplet_lifetime_gas', + ) + + singlet_lifetime_gas = straxen.URLConfig( + type=(int, float), + help='singlet_lifetime_gas', + ) + + triplet_lifetime_liquid = straxen.URLConfig( + type=(int, float), + help='triplet_lifetime_liquid', + ) + + singlet_lifetime_liquid = straxen.URLConfig( + type=(int, float), + help='singlet_lifetime_liquid', + ) + + s2_time_spread = straxen.URLConfig( + type=(int, float), + help='s2_time_spread', + ) + + s2_time_model = straxen.URLConfig( + help='s2_time_model', + ) + + singlet_fraction_gas = straxen.URLConfig( + type=(int, float), + help='singlet_fraction_gas', + ) + #needed as config? + phase_s2 = straxen.URLConfig( + default="gas", + help='phase_s2', + ) + + s2_luminescence_model = straxen.URLConfig( + help='s2_luminescence_model', + ) + + drift_velocity_liquid = straxen.URLConfig( + type=(int, float), + help='drift_velocity_liquid', + ) + + tpc_length = straxen.URLConfig( + type=(int, float), + help='tpc_length', + ) + + pmt_circuit_load_resistor = straxen.URLConfig( + type=(int, float), + help='pmt_circuit_load_resistor', + ) + + digitizer_bits = straxen.URLConfig( + type=(int, float), + help='digitizer_bits', + ) + + digitizer_voltage_range = straxen.URLConfig( + type=(int, float), + help='digitizer_voltage_range', + ) + + tpc_radius = straxen.URLConfig( + type=(int, float), + help='tpc_radius', + ) + + n_top_pmts = straxen.URLConfig( + type=(int), + help='Number of PMTs on top array', + ) + + n_tpc_pmts = straxen.URLConfig( + type=(int), + help='Number of PMTs in the TPC', + ) + + diffusion_constant_transverse = straxen.URLConfig( + type=(int, float), + help='diffusion_constant_transverse', + ) + + s2_aft_skewness = straxen.URLConfig( + type=(int, float), + help='s2_aft_skewness', + ) + + s2_aft_sigma = straxen.URLConfig( + type=(int, float), + help='s2_aft_sigma', + ) + + enable_field_dependencies = straxen.URLConfig( + help='enable_field_dependencies', + ) + + gains = straxen.URLConfig( + cache=True, + help='pmt gains', + ) + + s2_pattern_map = straxen.URLConfig( + cache=True, + help='s2_pattern_map', + ) + photon_area_distribution = straxen.URLConfig( + cache=True, + help='photon_area_distribution', + ) + + s2_optical_propagation_spline = straxen.URLConfig( + cache=True, + help='s2_optical_propagation_spline', + ) + + s2_luminescence_map = straxen.URLConfig( + cache=True, + help='s2_luminescence_map', + ) + + garfield_gas_gap_map = straxen.URLConfig( + cache=True, + help='garfield_gas_gap_map', + ) + #stupid naming problem... + field_dependencies_map_tmp = straxen.URLConfig( + help='field_dependencies_map', + ) + def setup(self): if self.debug: log.setLevel('DEBUG') log.debug("Running S2PhotonPropagation in debug mode") - - to_pe = straxen.get_resource(self.to_pe_file, fmt='npy') - self.to_pe = to_pe[0][1] - - adc_2_current = (self.digitizer_voltage_range - / 2 ** (self.digitizer_bits) - / self.pmt_circuit_load_resistor) - - self.gains = np.divide(adc_2_current, - self.to_pe, - out=np.zeros_like(self.to_pe), - where=self.to_pe != 0) self.pmt_mask = np.array(self.gains) > 0 # Converted from to pe (from cmt by default) self.turned_off_pmts = np.arange(len(self.gains))[np.array(self.gains) == 0] - - self.s2_pattern_map = make_patternmap(self.s2_pattern_map_file, fmt='pkl', pmt_mask=self.pmt_mask) - + #Move this part into a nice URLConfig protocol? # Field dependencies - # This config entry a dictionary of 5 items - self.enable_field_dependencies = config['enable_field_dependencies'] #This is not so nice if any(self.enable_field_dependencies.values()): - field_dependencies_map_tmp = make_map(self.field_dependencies_map, fmt='json.gz', method='RectBivariateSpline') self.drift_velocity_scaling = 1.0 # calculating drift velocity scaling to match total drift time for R=0 between cathode and gate if "norm_drift_velocity" in self.enable_field_dependencies.keys(): if self.enable_field_dependencies['norm_drift_velocity']: - norm_dvel = field_dependencies_map_tmp(np.array([ [0], [- self.tpc_length]]).T, map_name='drift_speed_map')[0] + norm_dvel = self.field_dependencies_map_tmp(np.array([ [0], [- self.tpc_length]]).T, map_name='drift_speed_map')[0] norm_dvel*=1e-4 drift_velocity_scaling = self.drift_velocity_liquid/norm_dvel def rz_map(z, xy, **kwargs): r = np.sqrt(xy[:, 0]**2 + xy[:, 1]**2) - return field_dependencies_map_tmp(np.array([r, z]).T, **kwargs) + return self.field_dependencies_map_tmp(np.array([r, z]).T, **kwargs) self.field_dependencies_map = rz_map - - if 'garfield_gas_gap' in self.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 - self.s2_luminescence_map = straxen.get_resource(os.path.join(private_files_path,"sim_files/garfield_timing_map_gas_gap_sr0.npy"), fmt='npy') - self.garfield_gas_gap_map = make_map(os.path.join(private_files_path,"sim_files/garfield_gas_gap_map_sr0.json"), fmt = 'json') - - if self.s2_time_spline: - self.s2_optical_propagation_spline = make_map(os.path.join(private_files_path,"sim_files/XENONnT_s2_opticalprop_time_v0.json.gz"), fmt="json.gz") - self.photon_area_distribution = straxen.get_resource(self.photon_area_distribution, fmt='csv') + self._cached_uniform_to_pe_arr = {} self.__uniform_to_pe_arr = self.init_spe_scaling_factor_distributions() - def compute(self, individual_electrons, electron_cloud): + def compute(self, individual_electrons, interactions_in_roi): + + #Just apply this to clusters with photons + mask = interactions_in_roi["n_electron_extracted"] > 0 if len(individual_electrons) == 0: return np.zeros(0, dtype=self.dtype) - positions = np.array([electron_cloud["x"], electron_cloud["y"]]).T + positions = np.array([interactions_in_roi[mask]["x"], interactions_in_roi[mask]["y"]]).T - _photon_channels = self.photon_channels(electron_cloud["n_electron_extracted"], - electron_cloud["z_obs"], + _photon_channels = self.photon_channels(interactions_in_roi[mask]["n_electron_extracted"], + interactions_in_roi[mask]["z_obs"], positions, - electron_cloud["drift_time_mean"] , - electron_cloud["sum_photons"], + interactions_in_roi[mask]["drift_time_mean"] , + interactions_in_roi[mask]["sum_s2_photons"], ) #_photon_channels = _photon_channels.astype(np.int64) _photon_timings = self.photon_timings(positions, - electron_cloud["sum_photons"], + interactions_in_roi[mask]["sum_s2_photons"], _photon_channels, ) #repeat for n photons per electron # Should this be before adding delays? - _photon_timings += np.repeat(individual_electrons["time"], individual_electrons["n_photons"]) + _photon_timings += np.repeat(individual_electrons["time"], individual_electrons["n_s2_photons"]) #Do i want to save both -> timings with and without pmt transition time spread? # Correct for PMT Transition Time Spread (skip for pmt after-pulses) @@ -417,7 +480,7 @@ def luminescence_timings_garfield_gasgap(self, xy, n_photons): def init_spe_scaling_factor_distributions(self): #This code will be duplicate with the corresponding S2 class # Improve!! - + config="PLACEHOLDER_FIX_THIS_PART" h = deterministic_hash(config) # What is this part doing? #if h in self._cached_uniform_to_pe_arr: # __uniform_to_pe_arr = self._cached_uniform_to_pe_arr[h] @@ -501,21 +564,4 @@ def draw_excitation_times(inv_cdf_list, hist_indices, nph, diff_nearest_gg, d_ga #subtract mean to get proper drift time and z correlation timings[count:count+n] = T count+=n - return timings - -#This is a modified version of the corresponding WFsim code.... -@njit() -def uniform_to_pe_arr(p, channel, __uniform_to_pe_arr): - indices = np.int64(p * 2000) + 1 - return __uniform_to_pe_arr[channel, indices] - -#In WFSim uniform_to_pe_arr is called inside a loop over the channels -#I needed to change the code to run on all channels at once -@njit() -def loop_uniform_to_pe_arr(p, channel, __uniform_to_pe_arr): - result = [] - for i in range(len(p)): - result.append(uniform_to_pe_arr(p[i], - channel=channel[i], - __uniform_to_pe_arr=__uniform_to_pe_arr) ) - return np.array(result) \ No newline at end of file + return timings \ No newline at end of file diff --git a/fuse/plugins/detector_physics/secondary_scintillation.py b/fuse/plugins/detector_physics/secondary_scintillation.py index 8f1bc825..b6854380 100644 --- a/fuse/plugins/detector_physics/secondary_scintillation.py +++ b/fuse/plugins/detector_physics/secondary_scintillation.py @@ -7,53 +7,11 @@ export, __all__ = strax.exporter() -from ...common import make_map, make_patternmap - logging.basicConfig(handlers=[logging.StreamHandler()]) log = logging.getLogger('fuse.detector_physics.secondary_scintillation') log.setLevel('WARNING') -#private_files_path = "path/to/private/files" -base_path = os.path.abspath(os.getcwd()) -private_files_path = os.path.join("/",*base_path.split("/")[:-2], "private_nt_aux_files") -config = straxen.get_resource(os.path.join(private_files_path, 'sim_files/fax_config_nt_sr0_v4.json') , fmt='json') - @export -@strax.takes_config( - strax.Option('p_double_pe_emision', default=config["p_double_pe_emision"], track=False, infer_type=False, - help="p_double_pe_emision"), - strax.Option('se_gain_from_map', default=config["se_gain_from_map"], track=False, infer_type=False, - help="se_gain_from_map"), - strax.Option('se_gain_map', - default=os.path.join(private_files_path,"strax_files/XENONnT_se_xy_map_v1_mlp.json"), - track=False, - infer_type=False, - help="se_gain_map"), - strax.Option('s2_correction_map_file', - default=os.path.join(private_files_path,"strax_files/XENONnT_s2_xy_map_v4_210503_mlp_3_in_1_iterated.json"), - track=False, - infer_type=False, - help="s2_correction_map"), - strax.Option('to_pe_file', default=os.path.join(private_files_path,"sim_files/to_pe_nt.npy"), track=False, infer_type=False, - help="to_pe file"), - strax.Option('digitizer_voltage_range', default=config['digitizer_voltage_range'], track=False, infer_type=False, - help="digitizer_voltage_range"), - strax.Option('digitizer_bits', default=config['digitizer_bits'], track=False, infer_type=False, - help="digitizer_bits"), - strax.Option('pmt_circuit_load_resistor', default=config['pmt_circuit_load_resistor'], track=False, infer_type=False, - help="pmt_circuit_load_resistor"), - strax.Option('s2_pattern_map_file', - default=os.path.join(private_files_path,"sim_files/XENONnT_s2_xy_patterns_GXe_LCE_corrected_qes_MCv4.3.0_wires.pkl"), - track=False, - infer_type=False, - help="s2_pattern_map"), - strax.Option('s2_secondary_sc_gain', default=config['s2_secondary_sc_gain'], track=False, infer_type=False, - help="s2_secondary_sc_gain"), - strax.Option('s2_gain_spread', default=0, track=False, infer_type=False, - help="s2_gain_spread"), - strax.Option('debug', default=False, track=False, infer_type=False, - help="Show debug informations"), -) class SecondaryScintillation(strax.Plugin): __version__ = "0.0.0" @@ -61,11 +19,11 @@ class SecondaryScintillation(strax.Plugin): depends_on = ("drifted_electrons","extracted_electrons" ,"electron_time") provides = ("s2_photons", "s2_photons_sum") data_kind = {"s2_photons": "individual_electrons", - "s2_photons_sum" : "electron_cloud" + "s2_photons_sum" : "interactions_in_roi" } - dtype_photons = [('n_photons', np.int64),] + strax.time_fields - dtype_sum_photons = [('sum_photons', np.int64),] + strax.time_fields + dtype_photons = [('n_s2_photons', np.int64),] + strax.time_fields + dtype_sum_photons = [('sum_s2_photons', np.int64),] + strax.time_fields dtype = dict() dtype["s2_photons"] = dtype_photons @@ -74,75 +32,126 @@ class SecondaryScintillation(strax.Plugin): #Forbid rechunking rechunk_on_save = False + #Config options + debug = straxen.URLConfig( + default=False, type=bool,track=False, + help='Show debug informations', + ) + + s2_gain_spread = straxen.URLConfig( + type=(int, float), + help='s2_gain_spread', + ) + + s2_secondary_sc_gain = straxen.URLConfig( + type=(int, float), + help='s2_secondary_sc_gain', + ) + + pmt_circuit_load_resistor = straxen.URLConfig( + type=(int, float), + help='pmt_circuit_load_resistor', + ) + + digitizer_bits = straxen.URLConfig( + type=(int, float), + help='digitizer_bits', + ) + + digitizer_voltage_range = straxen.URLConfig( + type=(int, float), + help='digitizer_voltage_range', + ) + + se_gain_from_map = straxen.URLConfig( + help='se_gain_from_map', + ) + + p_double_pe_emision = straxen.URLConfig( + type=(int, float), + help='p_double_pe_emision', + ) + + se_gain_map = straxen.URLConfig( + cache=True, + help='se_gain_map', + ) + + s2_correction_map = straxen.URLConfig( + cache=True, + help='s2_correction_map', + ) + + gains = straxen.URLConfig( + cache=True, + help='pmt gains', + ) + + s2_pattern_map = straxen.URLConfig( + cache=True, + help='s2_pattern_map', + ) + def setup(self): if self.debug: log.setLevel('DEBUG') log.debug("Running SecondaryScintillation in debug mode") - - if self.se_gain_from_map: - self.se_gain_map = make_map(self.se_gain_map, fmt = "json") - else: - if self.s2_correction_map_file: - self.s2_correction_map = make_map(self.s2_correction_map_file, fmt = 'json') - else: - - to_pe = straxen.get_resource(self.to_pe_file, fmt='npy') - self.to_pe = to_pe[0][1] - - adc_2_current = (self.digitizer_voltage_range - / 2 ** (self.digitizer_bits) - / self.pmt_circuit_load_resistor) - - gains = np.divide(adc_2_current, - self.to_pe, - out=np.zeros_like(self.to_pe), - where=self.to_pe != 0) - - self.pmt_mask = np.array(gains) > 0 # Converted from to pe (from cmt by default) - - self.s2_pattern_map = make_patternmap(self.s2_pattern_map_file, fmt='pkl', pmt_mask=self.pmt_mask) - - - 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 - s2cmap.data['map'] = np.sum(s2cmap.data['map'][:][:], axis=2, keepdims=True, where=self.pmt_mask) - # Scale by median value - s2cmap.data['map'] = s2cmap.data['map'] / np.median(s2cmap.data['map'][s2cmap.data['map'] > 0]) - s2cmap.__init__(s2cmap.data) - self.s2_correction_map = s2cmap + + self.pmt_mask = np.array(self.gains) + + #Are these if cases needed?? -> If no remove, if yes, correct the code + #if self.s2_correction_map_file: + # self.s2_correction_map = make_map(self.s2_correction_map_file, fmt = 'json') + #else: + + # self.pmt_mask = np.array(self.gains) > 0 # Converted from to pe (from cmt by default) + + # self.s2_pattern_map = make_patternmap(self.s2_pattern_map_file, fmt='pkl', pmt_mask=self.pmt_mask) + + + # 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 + # s2cmap.data['map'] = np.sum(s2cmap.data['map'][:][:], axis=2, keepdims=True, where=self.pmt_mask) + # Scale by median value + # s2cmap.data['map'] = s2cmap.data['map'] / np.median(s2cmap.data['map'][s2cmap.data['map'] > 0]) + # s2cmap.__init__(s2cmap.data) + # self.s2_correction_map = s2cmap - def compute(self, electron_cloud, individual_electrons ): + def compute(self, interactions_in_roi, individual_electrons): + + #Just apply this to clusters with photons + mask = interactions_in_roi["n_electron_extracted"] > 0 - if len(electron_cloud) == 0: + if len(interactions_in_roi[mask]) == 0: return dict(s2_photons=np.zeros(0, self.dtype["s2_photons"]), s2_photons_sum=np.zeros(0, self.dtype["s2_photons_sum"])) - positions = np.array([electron_cloud["x"], electron_cloud["y"]]).T + positions = np.array([interactions_in_roi[mask]["x"], interactions_in_roi[mask]["y"]]).T sc_gain = self.get_s2_light_yield(positions=positions) - electron_gains = np.repeat(sc_gain, electron_cloud["n_electron_extracted"]) + electron_gains = np.repeat(sc_gain, interactions_in_roi[mask]["n_electron_extracted"]) n_photons_per_ele = np.random.poisson(electron_gains) if self.s2_gain_spread: n_photons_per_ele += np.random.normal(0, self.s2_gain_spread, len(n_photons_per_ele)).astype(np.int64) - sum_photons_per_interaction = [np.sum(x) for x in np.split(n_photons_per_ele, np.cumsum(electron_cloud["n_electron_extracted"]))[:-1]] + sum_photons_per_interaction = [np.sum(x) for x in np.split(n_photons_per_ele, np.cumsum(interactions_in_roi[mask]["n_electron_extracted"]))[:-1]] n_photons_per_ele[n_photons_per_ele < 0] = 0 result_photons = np.zeros(len(n_photons_per_ele), dtype = self.dtype["s2_photons"]) - result_photons["n_photons"] = n_photons_per_ele + result_photons["n_s2_photons"] = n_photons_per_ele result_photons["time"] = individual_electrons["time"] result_photons["endtime"] = individual_electrons["endtime"] - result_sum_photons = np.zeros(len(sum_photons_per_interaction), dtype = self.dtype["s2_photons_sum"]) - result_sum_photons["sum_photons"] = sum_photons_per_interaction - result_sum_photons["time"] = electron_cloud["time"] - result_sum_photons["endtime"] = electron_cloud["endtime"] + result_sum_photons = np.zeros(len(interactions_in_roi), dtype = self.dtype["s2_photons_sum"]) + result_sum_photons["sum_s2_photons"][mask] = sum_photons_per_interaction + result_sum_photons["time"][mask] = interactions_in_roi[mask]["time"] + result_sum_photons["endtime"][mask] = interactions_in_roi[mask]["endtime"] return dict(s2_photons=result_photons, s2_photons_sum=result_sum_photons) diff --git a/fuse/plugins/micro_physics/__init__.py b/fuse/plugins/micro_physics/__init__.py index d3b69c3e..ac12e784 100644 --- a/fuse/plugins/micro_physics/__init__.py +++ b/fuse/plugins/micro_physics/__init__.py @@ -17,4 +17,7 @@ from .yields import * from . import wfsim_connection -from .wfsim_connection import * \ No newline at end of file +from .wfsim_connection import * + +from . import detector_volumes +from .detector_volumes import * \ No newline at end of file diff --git a/fuse/plugins/micro_physics/detector_volumes.py b/fuse/plugins/micro_physics/detector_volumes.py new file mode 100644 index 00000000..f90dceef --- /dev/null +++ b/fuse/plugins/micro_physics/detector_volumes.py @@ -0,0 +1,314 @@ +import strax +import numpy as np +import straxen +import logging +from ...volume_plugin import * +from ...vertical_merger_plugin import * + +logging.basicConfig(handlers=[logging.StreamHandler()]) +log = logging.getLogger('fuse.micro_physics.detector_volumes') +log.setLevel('WARNING') + + +class VolumesMerger(VerticalMergerPlugin): + """ + Plugin that concatenates the interactions that are in the + XENONnT TPC or the volume below the cathode. + """ + + depends_on = ("tpc_interactions", "below_cathode_interactions") + + provides = "interactions_in_roi" + data_kind = "interactions_in_roi" + __version__ = "0.0.0" + + #Forbid rechunking + rechunk_on_save = False + + def compute(self, **kwargs): + return super().compute(**kwargs) + + +# Fixed detector dimensions of XENONnT: +# See also: https://xe1t-wiki.lngs.infn.it/doku.php?id=xenon:xenonnt:analysis:coordinate_system + +class XENONnT_TPC(VolumePlugin): + """ + Plugin that evaluates if interactions are in the XENONnT TPC. + Only these interactions are returned. The output of this plugin is + meant to go into a vertical merger plugin. + """ + + depends_on = ("clustered_interactions") + + provides = "tpc_interactions" + data_kind = "tpc_interactions" + __version__ = "0.0.0" + + #Forbid rechunking + rechunk_on_save = False + + #Can we import this from MergeCluster and just add the needed fields? + dtype = [('x', np.float32), + ('y', np.float32), + ('z', np.float32), + ('ed', np.float64), + ('nestid', np.int64), + ('A', np.int64), + ('Z', np.int64), + ('evtid', np.int64), + ('x_pri', np.float32), + ('y_pri', np.float32), + ('z_pri', np.float32), + ('xe_density', np.float32), + ('vol_id', np.int64), + ('create_S2', np.bool8), + ] + + dtype = dtype + strax.time_fields + + #Config options + debug = straxen.URLConfig( + default=False, type=bool, + help='Show debug informations', + ) + #Define the TPC volume + xenonnt_z_cathode = straxen.URLConfig( + default = -148.6515, # cm ... top of the cathode electrode + type=(int, float), + help='xenonnt_z_cathode', + ) + + xenonnt_z_gate_mesh = straxen.URLConfig( + default = 0., # bottom of the gate electrode + type=(int, float), + help='xenonnt_z_gate_mesh', + ) + + xenonnt_sensitive_volume_radius = straxen.URLConfig( + default = 66.4, # cm + type=(int, float), + help='xenonnt_sensitive_volume_radius', + ) + + xenon_density_tpc = straxen.URLConfig( + default = 2.862, + type=(int, float), + help='xenon_density', + ) + + create_S2_xenonnt_TPC = straxen.URLConfig( + default=True, type=bool, + help='Create S2s in the XENONnT TPC', + ) + + def setup(self): + + if self.debug: + log.setLevel('DEBUG') + log.debug("Running XENONnT_TPC in debug mode") + + def compute(self, clustered_interactions): + + #Call the ROI function: + mask = self.in_ROI(clustered_interactions, + self.xenonnt_z_cathode, + self.xenonnt_z_gate_mesh, + self.xenonnt_sensitive_volume_radius + ) + + tpc_interactions = clustered_interactions[mask] + tpc_interactions["xe_density"] = self.xenon_density_tpc + tpc_interactions["vol_id"] = 1 #Do we need this? -> Now just copied from epix + tpc_interactions["create_S2"] = self.create_S2_xenonnt_TPC + + return tpc_interactions + + + +class XENONnT_BelowCathode(VolumePlugin): + """ + Plugin that evaluates if interactions are below the Cathode of XENONnT. + Only these interactions are returned. The output of this plugin is + meant to go into a vertical merger plugin. + """ + + depends_on = ("clustered_interactions") + + provides = "below_cathode_interactions" + data_kind = "below_cathode_interactions" + __version__ = "0.0.0" + + #Forbid rechunking + rechunk_on_save = False + + #Can we import this from MergeCluster and just add the needed fields? + dtype = [('x', np.float32), + ('y', np.float32), + ('z', np.float32), + ('ed', np.float64), + ('nestid', np.int64), + ('A', np.int64), + ('Z', np.int64), + ('evtid', np.int64), + ('x_pri', np.float32), + ('y_pri', np.float32), + ('z_pri', np.float32), + ('xe_density', np.float32), + ('vol_id', np.int64), + ('create_S2', np.bool8), + ] + + dtype = dtype + strax.time_fields + + #Config options + debug = straxen.URLConfig( + default=False, type=bool, + help='Show debug informations', + ) + #Define the volume + xenonnt_z_cathode = straxen.URLConfig( + default = -148.6515, # cm ... top of the cathode electrode + type=(int, float), + help='xenonnt_z_cathode', + ) + + xenonnt_z_bottom_pmts = straxen.URLConfig( + default = -154.6555, # cm ... top surface of the bottom PMT window + type=(int, float), + help='xenonnt_z_bottom_pmts', + ) + + xenonnt_sensitive_volume_radius = straxen.URLConfig( + default = 66.4, # cm + type=(int, float), + help='xenonnt_sensitive_volume_radius', + ) + + xenon_density_below_cathode = straxen.URLConfig( + default = 2.862, + type=(int, float), + help='xenon_density', + ) + + create_S2_xenonnt_below_cathode = straxen.URLConfig( + default=False, type=bool, + help='No S2s below the Cathode!', + ) + + def setup(self): + + if self.debug: + log.setLevel('DEBUG') + log.debug("Running XENONnT_BelowCathode in debug mode") + + def compute(self, clustered_interactions): + + #Call the ROI function: + mask = self.in_ROI(clustered_interactions, + self.xenonnt_z_bottom_pmts, + self.xenonnt_z_cathode, + self.xenonnt_sensitive_volume_radius + ) + + tpc_interactions = clustered_interactions[mask] + tpc_interactions["xe_density"] = self.xenon_density_below_cathode + tpc_interactions["vol_id"] = 2 #Do we need this? -> Now just copied from epix + tpc_interactions["create_S2"] = self.create_S2_xenonnt_below_cathode + + return tpc_interactions + + + +class XENONnT_GasPhase(VolumePlugin): + """ + Plugin that evaluates if interactions are in the gas phase of XENONnT. + Only these interactions are returned. The output of this plugin is + meant to go into a vertical merger plugin. + """ + + depends_on = ("clustered_interactions") + + provides = "gas_phase_interactions" + data_kind = "gas_phase_interactions" + __version__ = "0.0.0" + + #Forbid rechunking + rechunk_on_save = False + + #Can we import this from MergeCluster and just add the needed fields? + dtype = [('x', np.float32), + ('y', np.float32), + ('z', np.float32), + ('ed', np.float64), + ('nestid', np.int64), + ('A', np.int64), + ('Z', np.int64), + ('evtid', np.int64), + ('x_pri', np.float32), + ('y_pri', np.float32), + ('z_pri', np.float32), + ('xe_density', np.float32), + ('vol_id', np.int64), + ('create_S2', np.bool8), + ] + + dtype = dtype + strax.time_fields + + #Config options + debug = straxen.URLConfig( + default=False, type=bool, + help='Show debug informations', + ) + #Define the volume + xenonnt_z_top_pmts = straxen.URLConfig( + default = 7.3936, # cm + type=(int, float), + help='xenonnt_z_top_pmts', + ) + + xenonnt_z_lxe = straxen.URLConfig( + default = 0.416, # cm ... liquid-gas interface + type=(int, float), + help='xenonnt_z_lxe', + ) + + xenonnt_sensitive_volume_radius = straxen.URLConfig( + default = 66.4, # cm + type=(int, float), + help='xenonnt_sensitive_volume_radius', + ) + + xenon_density_gas_phase = straxen.URLConfig( + default = 0.0177, + type=(int, float), + help='xenon_density', + ) + + create_S2_xenonnt_gas_phase = straxen.URLConfig( + default=False, type=bool, + help='No S2s in gas', + ) + + def setup(self): + + if self.debug: + log.setLevel('DEBUG') + log.debug("Running XENONnT_GasPhase in debug mode") + + def compute(self, clustered_interactions): + + #Call the ROI function: + mask = self.in_ROI(clustered_interactions, + self.xenonnt_z_lxe, + self.xenonnt_z_top_pmts, + self.xenonnt_sensitive_volume_radius + ) + + tpc_interactions = clustered_interactions[mask] + tpc_interactions["xe_density"] = self.xenon_density_gas_phase + tpc_interactions["vol_id"] = 3 #Do we need this? -> Now just copied from epix + tpc_interactions["create_S2"] = self.create_S2_xenonnt_gas_phase + + return tpc_interactions + diff --git a/fuse/plugins/micro_physics/electric_field.py b/fuse/plugins/micro_physics/electric_field.py index 3f17b78d..56f1fca1 100644 --- a/fuse/plugins/micro_physics/electric_field.py +++ b/fuse/plugins/micro_physics/electric_field.py @@ -1,7 +1,7 @@ import strax -import epix import numpy as np import logging +import straxen export, __all__ = strax.exporter() @@ -10,14 +10,6 @@ log.setLevel('WARNING') @export -@strax.takes_config( - strax.Option('debug', default=False, track=False, infer_type=False, - help="Show debug information"), - strax.Option('detector', default="XENONnT", track=False, infer_type=False, - help="Detector to be used. Has to be defined in epix.detectors"), - strax.Option('detector_config_override', default=None, track=False, infer_type=False, - help="Config file to overwrite default epix.detectors settings; see examples in the configs folder") -) class ElectricField(strax.Plugin): """ Plugin that calculates the electric field values for the detector. @@ -25,9 +17,9 @@ class ElectricField(strax.Plugin): __version__ = "0.0.0" - depends_on = ("clustered_interactions",) + depends_on = ("interactions_in_roi",) provides = "electric_field_values" - data_kind = "clustered_interactions" + data_kind = "interactions_in_roi" #Forbid rechunking rechunk_on_save = False @@ -37,52 +29,43 @@ class ElectricField(strax.Plugin): *strax.time_fields ] - def setup(self): - """ - Do the volume cuts here. + #Config options + debug = straxen.URLConfig( + default=False, type=bool,track=False, + help='Show debug informations', + ) - Initialize the detector config. - """ - self.detector_config = epix.init_detector( - self.detector.lower(), - self.detector_config_override - ) + efield_map = straxen.URLConfig( + cache=True, + help='electric field map', + ) + + def setup(self): if self.debug: log.setLevel('DEBUG') log.debug("Running ElectricField in debug mode") - def compute(self, clustered_interactions): + def compute(self, interactions_in_roi): """ Calculate the electric field values for the given clustered interactions. Args: - clustered_interactions (numpy.ndarray): array of clustered interactions. + interactions_in_roi (numpy.ndarray): array of clustered interactions. Returns: numpy.ndarray: array of electric field values. """ - if len(clustered_interactions) == 0: + if len(interactions_in_roi) == 0: return np.zeros(0, dtype=self.dtype) - electric_field_array = np.zeros(len(clustered_interactions), dtype=self.dtype) - electric_field_array['time'] = clustered_interactions['time'] - electric_field_array['endtime'] = clustered_interactions['endtime'] - - efields = electric_field_array['e_field'] - - for volume in self.detector_config: - if isinstance(volume.electric_field, (float, int)): - volume_id_mask = clustered_interactions['vol_id'] == volume.volume_id - efields[volume_id_mask] = volume.electric_field - else: - efields = volume.electric_field( - clustered_interactions['x'], - clustered_interactions['y'], - clustered_interactions['z'] - ) + electric_field_array = np.zeros(len(interactions_in_roi), dtype=self.dtype) + electric_field_array['time'] = interactions_in_roi['time'] + electric_field_array['endtime'] = interactions_in_roi['endtime'] - electric_field_array['e_field'] = efields + r = np.sqrt(interactions_in_roi['x'] ** 2 + interactions_in_roi['y'] ** 2) + positions = np.stack((r, interactions_in_roi['z']), axis=1) + electric_field_array['e_field'] = self.efield_map(positions) return electric_field_array diff --git a/fuse/plugins/micro_physics/find_cluster.py b/fuse/plugins/micro_physics/find_cluster.py index 2d3064d8..0ac500e3 100644 --- a/fuse/plugins/micro_physics/find_cluster.py +++ b/fuse/plugins/micro_physics/find_cluster.py @@ -4,6 +4,7 @@ import strax import awkward as ak import logging +import straxen from sklearn.cluster import DBSCAN export, __all__ = strax.exporter() @@ -15,14 +16,6 @@ log.setLevel('WARNING') @export -@strax.takes_config( - strax.Option('micro_separation', default=0.005, track=False, infer_type=False, - help="DBSCAN clustering distance (mm)"), - strax.Option('micro_separation_time', default = 10, track=False, infer_type=False, - help="Clustering time (ns)"), - strax.Option('debug', default=False, track=False, infer_type=False, - help="Show debug informations"), -) class FindCluster(strax.Plugin): __version__ = "0.0.0" @@ -38,6 +31,22 @@ class FindCluster(strax.Plugin): #Forbid rechunking rechunk_on_save = False + #Config options + debug = straxen.URLConfig( + default=False, type=bool,track=False, + help='Show debug informations', + ) + + micro_separation_time = straxen.URLConfig( + default=10, type=(int, float), + help='Clustering time (ns)', + ) + + micro_separation = straxen.URLConfig( + default=0.005, type=(int, float), + help='DBSCAN clustering distance (mm)', + ) + def setup(self): if self.debug: diff --git a/fuse/plugins/micro_physics/input.py b/fuse/plugins/micro_physics/input.py index 079cc107..e470abad 100644 --- a/fuse/plugins/micro_physics/input.py +++ b/fuse/plugins/micro_physics/input.py @@ -1,4 +1,5 @@ import strax +import straxen import uproot import os import numba @@ -10,43 +11,14 @@ export, __all__ = strax.exporter() -from ...common import full_array_to_numpy, reshape_awkward - -import epix +from ...common import full_array_to_numpy, reshape_awkward, dynamic_chunking logging.basicConfig(handlers=[logging.StreamHandler()]) log = logging.getLogger('fuse.micro_physics.input') log.setLevel('WARNING') +#Remove the path and file name option from the config and do this with the run_number?? @export -@strax.takes_config( - strax.Option('path', default=".", track=False, infer_type=False, - help="Path to search for data"), - strax.Option('file_name', track=False, infer_type=False, - help="File to open"), - strax.Option('separation_scale', default=1e8, track=False, infer_type=False, - help="Add Description"), - strax.Option('source_rate', default=1, track=False, infer_type=False, - help="source_rate"), - strax.Option('cut_delayed', default=4e12, track=False, infer_type=False, - help="cut_delayed"), - strax.Option('n_interactions_per_chunk', default=10000, track=False, infer_type=False, - help="Add n_interactions_per_chunk"), - strax.Option('debug', default=False, track=False, infer_type=False, - help="Show debug informations"), - strax.Option('entry_start', default=0, track=False, infer_type=False, - help="First event to be read"), - strax.Option('entry_stop', default=None, track=False, infer_type=False, - help="How many entries from the ROOT file you want to process. I think it is not working at the moment"), - strax.Option('cut_by_eventid', default=False, track=False, infer_type=False, - help="If selected, the next two arguments act on the G4 event id, and not the entry number (default)"), - strax.Option('nr_only', default=False, track=False, infer_type=False, - help="Add if you want to filter only nuclear recoil events (maximum ER energy deposit 10 keV)"), - strax.Option('detector', default="XENONnT", track=False, infer_type=False, - help="Detector to be used. Has to be defined in epix.detectors"), - strax.Option('detector_config_override', default=None, track=False, infer_type=False, - help="Config file to overwrite default epix.detectors settings; see examples in the configs folder"), -) class ChunkInput(strax.Plugin): __version__ = "0.0.0" @@ -80,16 +52,67 @@ class ChunkInput(strax.Plugin): source_done = False + #Config options + debug = straxen.URLConfig( + default=False, type=bool,track=False, + help='Show debug informations', + ) + + path = straxen.URLConfig( + track=False, + help='Path to search for data', + ) + + file_name = straxen.URLConfig( + track=False, + help='File to open', + ) + + separation_scale = straxen.URLConfig( + default=1e8, type=(int, float), + help='separation_scale', + ) + + source_rate = straxen.URLConfig( + default=1, type=(int, float), + help='source_rate', + ) + + cut_delayed = straxen.URLConfig( + default=4e14, type=(int, float), + help='cut_delayed', + ) + + n_interactions_per_chunk = straxen.URLConfig( + default=1e5, type=(int, float), + help='n_interactions_per_chunk', + ) + + entry_start = straxen.URLConfig( + default=0, type=(int, float), + help='entry_start', + ) + + entry_stop = straxen.URLConfig( + default=None, + help='How many entries from the ROOT file you want to process. I think it is not working at the moment', + ) + + cut_by_eventid = straxen.URLConfig( + default=False, type=bool, + help='If selected, the next two arguments act on the G4 event id, and not the entry number (default)', + ) + + nr_only = straxen.URLConfig( + default=False, type=bool, + help='Add if you want to filter only nuclear recoil events (maximum ER energy deposit 10 keV)', + ) + def setup(self): if self.debug: log.setLevel('DEBUG') log.debug("Running ChunkInput in debug mode") - - #Do the volume cuts here #Maybe we can move these lines somewhere else? - self.detector_config = epix.init_detector(self.detector.lower(), self.detector_config_override) - outer_cylinder = getattr(epix.detectors, self.detector.lower()) - outer_cylinder = outer_cylinder() self.file_reader = file_loader(self.path, self.file_name, @@ -106,8 +129,6 @@ def setup(self): ) self.file_reader_iterator = self.file_reader.output_chunk() - - def compute(self): try: @@ -470,33 +491,4 @@ def _awkwardify_df(df): "evtid": reshape_awkward(df["evtid"].values , evt_offsets), } - return ak.Array(dictionary) - - - -@numba.njit() -def dynamic_chunking(data, scale, n_min): - - idx_sort = np.argsort(data) - idx_undo_sort = np.argsort(idx_sort) - - data_sorted = data[idx_sort] - - diff = data_sorted[1:] - data_sorted[:-1] - - clusters = np.array([0]) - c = 0 - for value in diff: - if value <= scale: - clusters = np.append(clusters, c) - - elif len(clusters[clusters == c]) < n_min: - clusters = np.append(clusters, c) - - elif value > scale: - c = c + 1 - clusters = np.append(clusters, c) - - clusters_undo_sort = clusters[idx_undo_sort] - - return clusters_undo_sort \ No newline at end of file + return ak.Array(dictionary) \ No newline at end of file diff --git a/fuse/plugins/micro_physics/merge_cluster.py b/fuse/plugins/micro_physics/merge_cluster.py index ff439c2a..915fae3d 100644 --- a/fuse/plugins/micro_physics/merge_cluster.py +++ b/fuse/plugins/micro_physics/merge_cluster.py @@ -1,6 +1,6 @@ import strax +import straxen import numpy as np -import epix import awkward as ak import numba import logging @@ -14,21 +14,6 @@ log.setLevel('WARNING') @export -@strax.takes_config( - strax.Option('debug', default=False, track=False, infer_type=False, - help="Show debug informations"), - strax.Option('tag_cluster_by', default=False, track=False, infer_type=False, - help="decide if you tag the cluster (particle type, energy depositing process)\ - according to first interaction in it (time) or most energetic (energy)"), - strax.Option('detector', default="XENONnT", track=False, infer_type=False, - help="Detector to be used. Has to be defined in epix.detectors"), - strax.Option('detector_config_override', default=None, track=False, infer_type=False, - help="Config file to overwrite default epix.detectors settings; see examples in the configs folder"), - strax.Option('max_delay', default=1e7, track=False, infer_type=False, - help="Time after which we cut the rest of the event (ns)"), - strax.Option('debug', default=False, track=False, infer_type=False, - help="Show debug informations"), -) class MergeCluster(strax.Plugin): __version__ = "0.0.0" @@ -52,12 +37,30 @@ class MergeCluster(strax.Plugin): ('x_pri', np.float32), ('y_pri', np.float32), ('z_pri', np.float32), - ('xe_density', np.float32), - ('vol_id', np.int64), - ('create_S2', np.bool8), + ('xe_density', np.float32), #Will be set i a later plugin + ('vol_id', np.int64), #Will be set i a later plugin + ('create_S2', np.bool8), #Will be set i a later plugin ] dtype = dtype + strax.time_fields + + #Config options + debug = straxen.URLConfig( + default=False, type=bool,track=False, + help='Show debug informations', + ) + + tag_cluster_by = straxen.URLConfig( + default=False, type=bool, + help='decide if you tag the cluster (particle type, energy depositing process)\ + according to first interaction in it (time) or most energetic (energy))', + ) + + #Check this variable again in combination with cut_delayed + max_delay = straxen.URLConfig( + default=1e7, type=(int, float), + help='Time after which we cut the rest of the event (ns)', + ) def setup(self): @@ -65,9 +68,6 @@ def setup(self): log.setLevel('DEBUG') log.debug("Running MergeCluster in debug mode") - #Do the volume cuts here #Maybe we can move these lines somewhere else? - self.detector_config = epix.init_detector(self.detector.lower(), self.detector_config_override) - def compute(self, geant4_interactions): if len(geant4_interactions) == 0: @@ -87,17 +87,6 @@ def compute(self, geant4_interactions): result['y_pri'] = ak.broadcast_arrays(inter['y_pri'][:, 0], result['ed'])[0] result['z_pri'] = ak.broadcast_arrays(inter['z_pri'][:, 0], result['ed'])[0] - - res_det = epix.in_sensitive_volume(result, self.detector_config) - for field in res_det.fields: - result[field] = res_det[field] - m = result['vol_id'] > 0 # All volumes have an id larger zero - result = result[m] - - # Removing now empty events as a result of the selection above: - m = ak_num(result['ed']) > 0 - result = result[m] - # Sort entries (in an event) by in time, then chop all delayed # events which are too far away from the rest. # (This is a requirement of WFSim) diff --git a/fuse/plugins/micro_physics/microphysics_summary.py b/fuse/plugins/micro_physics/microphysics_summary.py index 4003779a..367c1b4f 100644 --- a/fuse/plugins/micro_physics/microphysics_summary.py +++ b/fuse/plugins/micro_physics/microphysics_summary.py @@ -8,7 +8,7 @@ class MicroPhysicsSummary(strax.MergeOnlyPlugin): Plugin which summarizes the MicroPhysics simulation into a single output """ - depends_on = ['clustered_interactions', + depends_on = ['interactions_in_roi', 'quanta', 'electric_field_values', ] @@ -17,10 +17,4 @@ class MicroPhysicsSummary(strax.MergeOnlyPlugin): __version__ = '0.0.0' #Forbid rechunking - rechunk_on_save = False - - def compute(self, **kwargs): - - microphysics_summary = super().compute(**kwargs) - - return microphysics_summary \ No newline at end of file + rechunk_on_save = False \ No newline at end of file diff --git a/fuse/plugins/micro_physics/wfsim_connection.py b/fuse/plugins/micro_physics/wfsim_connection.py index 6df92c59..d0769904 100644 --- a/fuse/plugins/micro_physics/wfsim_connection.py +++ b/fuse/plugins/micro_physics/wfsim_connection.py @@ -2,6 +2,7 @@ #We can keep this for validation of fuse but it can be removed later on import strax +import straxen import numpy as np import awkward as ak import logging @@ -15,10 +16,6 @@ log.setLevel('WARNING') @export -@strax.takes_config( - strax.Option('debug', default=False, track=False, infer_type=False, - help="Show debug informations"), -) class output_plugin(strax.Plugin): __version__ = "0.0.0" @@ -49,6 +46,12 @@ class output_plugin(strax.Plugin): (('Y position of the primary particle [cm]', 'y_pri'), np.float32), (('Z position of the primary particle [cm]', 'z_pri'), np.float32), ] + + #Config options + debug = straxen.URLConfig( + default=False, type=bool,track=False, + help='Show debug informations', + ) def setup(self): diff --git a/fuse/plugins/micro_physics/yields.py b/fuse/plugins/micro_physics/yields.py index f469f5ad..39c11c67 100644 --- a/fuse/plugins/micro_physics/yields.py +++ b/fuse/plugins/micro_physics/yields.py @@ -1,6 +1,7 @@ import numpy as np import nestpy import strax +import straxen import logging import pickle @@ -11,17 +12,13 @@ log.setLevel('WARNING') @export -@strax.takes_config( - strax.Option('debug', default=False, track=False, infer_type=False, - help="Show debug informations"), -) class NestYields(strax.Plugin): __version__ = "0.0.0" - depends_on = ["clustered_interactions", "electric_field_values"] + depends_on = ["interactions_in_roi", "electric_field_values"] provides = "quanta" - data_kind = "clustered_interactions" + data_kind = "interactions_in_roi" dtype = [('photons', np.float64), ('electrons', np.float64), @@ -32,6 +29,12 @@ class NestYields(strax.Plugin): #Forbid rechunking rechunk_on_save = False + + #Config options + debug = straxen.URLConfig( + default=False, type=bool,track=False, + help='Show debug informations', + ) def setup(self): if self.debug: @@ -41,33 +44,33 @@ def setup(self): self.quanta_from_NEST = np.vectorize(self._quanta_from_NEST) - def compute(self, clustered_interactions): + def compute(self, interactions_in_roi): """ Computes the charge and light quanta for a list of clustered interactions. Args: - clustered_interactions (numpy.ndarray): An array of clustered interactions. + interactions_in_roi (numpy.ndarray): An array of clustered interactions. Returns: numpy.ndarray: An array of quanta, with fields for time, endtime, photons, electrons, and excitons. """ - if len(clustered_interactions) == 0: + if len(interactions_in_roi) == 0: return np.zeros(0, dtype=self.dtype) - result = np.zeros(len(clustered_interactions), dtype=self.dtype) - result["time"] = clustered_interactions["time"] - result["endtime"] = clustered_interactions["endtime"] + result = np.zeros(len(interactions_in_roi), dtype=self.dtype) + result["time"] = interactions_in_roi["time"] + result["endtime"] = interactions_in_roi["endtime"] # Generate quanta: - if len(clustered_interactions) > 0: + if len(interactions_in_roi) > 0: photons, electrons, excitons = self.quanta_from_NEST( - clustered_interactions['ed'], - clustered_interactions['nestid'], - clustered_interactions['e_field'], - clustered_interactions['A'], - clustered_interactions['Z'], - clustered_interactions['create_S2'], - density=clustered_interactions['xe_density'] + interactions_in_roi['ed'], + interactions_in_roi['nestid'], + interactions_in_roi['e_field'], + interactions_in_roi['A'], + interactions_in_roi['Z'], + interactions_in_roi['create_S2'], + density=interactions_in_roi['xe_density'] ) result['photons'] = photons result['electrons'] = electrons @@ -102,7 +105,6 @@ def _quanta_from_NEST(en, model, e_field, A, Z, create_s2, **kwargs): excitons (numpy.array): Number of generated excitons """ nc = nestpy.NESTcalc(nestpy.VDetector()) - density = 2.862 # g/cm^3 # Fix for Kr83m events. # Energies have to be very close to 32.1 keV or 9.4 keV @@ -145,29 +147,13 @@ def _quanta_from_NEST(en, model, e_field, A, Z, create_s2, **kwargs): return photons, electrons, excitons - - -@strax.takes_config( - strax.Option('debug', default=False, track=False, infer_type=False, - help="Show debug informations"), - strax.Option('use_recombination_fluctuation', default = True, track=False, infer_type=False, - help="use_recombination_fluctuation"), - strax.Option('g1_value', default=0.151, track=False, infer_type=False, - help="g1_value"), - strax.Option('g2_value', default=16.45, track=False, infer_type=False, - help="g2_value"), - strax.Option('cs1_spline_path', default="/project2/lgrandi/pkavrigin/2023-04-24_epix_data_files/cs1_func_E_option2.pkl", track=False, infer_type=False, - help="cs1_spline_path"), - strax.Option('cs2_spline_path', default="/project2/lgrandi/pkavrigin/2023-04-24_epix_data_files/cs2_func_E_option2.pkl", track=False, infer_type=False, - help="cs2_spline_path"), -) class BetaYields(strax.Plugin): __version__ = "0.0.0" - depends_on = ["clustered_interactions", "electric_field_values"] + depends_on = ["interactions_in_roi", "electric_field_values"] provides = "quanta" - data_kind = "clustered_interactions" + data_kind = "interactions_in_roi" dtype = [('photons', np.float64), ('electrons', np.float64), @@ -178,15 +164,45 @@ class BetaYields(strax.Plugin): #Forbid rechunking rechunk_on_save = False + + #Config options + debug = straxen.URLConfig( + default=False, type=bool, + help='Show debug informations', + ) + + use_recombination_fluctuation = straxen.URLConfig( + default=True, type=bool, + help='use_recombination_fluctuation', + ) + + g1_value = straxen.URLConfig( + type=(int, float), + help='g1', + ) + + g2_value = straxen.URLConfig( + type=(int, float), + help='g2', + ) + + cs1_spline_path = straxen.URLConfig( + help='cs1_spline_path', + ) + + cs2_spline_path = straxen.URLConfig( + help='cs2_spline_path', + ) def setup(self): if self.debug: log.setLevel('DEBUG') - log.debug("Running NestYields in debug mode") + log.debug("Running BetaYields in debug mode") log.debug(f'Using nestpy version {nestpy.__version__}') self.get_quanta_vectorized = np.vectorize(self.get_quanta, excluded="self") + #This can be moved into an URLConfig protocol before merging the PR! with open(self.cs1_spline_path, 'rb') as f: self.cs1_spline = pickle.load(f) with open(self.cs2_spline_path, 'rb') as f: @@ -196,24 +212,24 @@ def setup(self): for i in range(np.random.randint(100)): self.nc.GetQuanta(self.nc.GetYields(energy=np.random.uniform(10, 100))) - def compute(self, clustered_interactions): + def compute(self, interactions_in_roi): """ Computes the charge and light quanta for a list of clustered interactions using custom yields. Args: - clustered_interactions (numpy.ndarray): An array of clustered interactions. + interactions_in_roi (numpy.ndarray): An array of clustered interactions. Returns: numpy.ndarray: An array of quanta, with fields for time, endtime, photons, electrons, and excitons. """ - if len(clustered_interactions) == 0: + if len(interactions_in_roi) == 0: return np.zeros(0, dtype=self.dtype) - result = np.zeros(len(clustered_interactions), dtype=self.dtype) - result["time"] = clustered_interactions["time"] - result["endtime"] = clustered_interactions["endtime"] + result = np.zeros(len(interactions_in_roi), dtype=self.dtype) + result["time"] = interactions_in_roi["time"] + result["endtime"] = interactions_in_roi["endtime"] - photons, electrons, excitons = self.get_quanta_vectorized(clustered_interactions["ed"], clustered_interactions["e_field"]) + photons, electrons, excitons = self.get_quanta_vectorized(interactions_in_roi["ed"], interactions_in_roi["e_field"]) result['photons'] = photons result['electrons'] = electrons result['excitons'] = excitons @@ -245,18 +261,11 @@ def get_quanta(self, energy, field): - - - -@strax.takes_config( - strax.Option('debug', default=False, track=False, infer_type=False, - help="Show debug informations"), -) class BBFYields(strax.Plugin): __version__ = "0.0.0" - depends_on = ["clustered_interactions", "electric_field_values"] + depends_on = ["interactions_in_roi", "electric_field_values"] provides = "quanta" dtype = [('photons', np.float64), @@ -266,6 +275,12 @@ class BBFYields(strax.Plugin): dtype = dtype + strax.time_fields + #Config options + debug = straxen.URLConfig( + default=False, type=bool, + help='Show debug informations', + ) + def setup(self): self.bbfyields = BBF_quanta_generator() @@ -273,19 +288,19 @@ def setup(self): log.setLevel("DEBUG") log.debug("Running BBFYields in debug mode") - def compute(self, geant4_interactions): + def compute(self, interactions_in_roi): - result = np.zeros(len(geant4_interactions), dtype=self.dtype) - result["time"] = geant4_interactions["time"] - result["endtime"] = geant4_interactions["endtime"] + result = np.zeros(len(interactions_in_roi), dtype=self.dtype) + result["time"] = interactions_in_roi["time"] + result["endtime"] = interactions_in_roi["endtime"] # Generate quanta: - if len(geant4_interactions) > 0: + if len(interactions_in_roi) > 0: photons, electrons, excitons = self.bbfyields.get_quanta_vectorized( - energy=geant4_interactions['ed'], - interaction=geant4_interactions['nestid'], - field=geant4_interactions['e_field'] + energy=interactions_in_roi['ed'], + interaction=interactions_in_roi['nestid'], + field=interactions_in_roi['e_field'] ) diff --git a/fuse/plugins/pmt_and_daq/__init__.py b/fuse/plugins/pmt_and_daq/__init__.py index 715e4658..61492f83 100644 --- a/fuse/plugins/pmt_and_daq/__init__.py +++ b/fuse/plugins/pmt_and_daq/__init__.py @@ -2,4 +2,7 @@ from .pmt_afterpulses import * from . import pmt_response_and_daq -from .pmt_response_and_daq import * \ No newline at end of file +from .pmt_response_and_daq import * + +from . import photon_summary +from .photon_summary import * \ No newline at end of file diff --git a/fuse/plugins/pmt_and_daq/photon_summary.py b/fuse/plugins/pmt_and_daq/photon_summary.py new file mode 100644 index 00000000..e99f073e --- /dev/null +++ b/fuse/plugins/pmt_and_daq/photon_summary.py @@ -0,0 +1,20 @@ +import strax + +from ...vertical_merger_plugin import VerticalMergerPlugin + +export, __all__ = strax.exporter() + +@export +class PhotonSummary(VerticalMergerPlugin): + """ + Plugin which concatenates the output of the S1, S2 and AP plugins + """ + + depends_on = ("propagated_s2_photons", "propagated_s1_photons", "pmt_afterpulses") + + provides = 'photon_summary' + data_kind = 'propagated_photons' + __version__ = '0.0.0' + + #Forbid rechunking + rechunk_on_save = False \ No newline at end of file diff --git a/fuse/plugins/pmt_and_daq/pmt_afterpulses.py b/fuse/plugins/pmt_and_daq/pmt_afterpulses.py index a4094cf6..bb09a31a 100644 --- a/fuse/plugins/pmt_and_daq/pmt_afterpulses.py +++ b/fuse/plugins/pmt_and_daq/pmt_afterpulses.py @@ -10,32 +10,7 @@ log = logging.getLogger('fuse.pmt_and_daq.pmt_afterpulses') log.setLevel('WARNING') -#private_files_path = "path/to/private/files" -base_path = os.path.abspath(os.getcwd()) -private_files_path = os.path.join("/",*base_path.split("/")[:-2], "private_nt_aux_files") -config = straxen.get_resource(os.path.join(private_files_path, 'sim_files/fax_config_nt_sr0_v4.json') , fmt='json') - @export -@strax.takes_config( - strax.Option('enable_pmt_afterpulses', default=config["enable_pmt_afterpulses"], track=False, infer_type=False, - help="enable_pmt_afterpulses"), - strax.Option('photon_ap_cdfs', default=config["photon_ap_cdfs"], track=False, infer_type=False, - help="photon_ap_cdfs"), - strax.Option('to_pe_file', default=os.path.join(private_files_path,"sim_files/to_pe_nt.npy"), track=False, infer_type=False, - help="to_pe file"), - strax.Option('digitizer_voltage_range', default=config['digitizer_voltage_range'], track=False, infer_type=False, - help="digitizer_voltage_range"), - strax.Option('digitizer_bits', default=config['digitizer_bits'], track=False, infer_type=False, - help="digitizer_bits"), - strax.Option('pmt_circuit_load_resistor', default=config['pmt_circuit_load_resistor'], track=False, infer_type=False, - help="pmt_circuit_load_resistor"), - strax.Option('pmt_ap_modifier', default=config['pmt_ap_modifier'], track=False, infer_type=False, - help="pmt_ap_modifier"), - strax.Option('pmt_ap_t_modifier', default=config['pmt_ap_t_modifier'], track=False, infer_type=False, - help="pmt_ap_t_modifier"), - strax.Option('debug', default=False, track=False, infer_type=False, - help="Show debug informations"), -) class PMTAfterPulses(strax.Plugin): __version__ = "0.0.0" @@ -54,32 +29,59 @@ class PMTAfterPulses(strax.Plugin): ] dtype = dtype + strax.time_fields + #Config options + debug = straxen.URLConfig( + default=False, type=bool,track=False, + help='Show debug informations', + ) + + pmt_ap_t_modifier = straxen.URLConfig( + type=(int, float), + help='pmt_ap_t_modifier', + ) + + pmt_ap_modifier = straxen.URLConfig( + type=(int, float), + help='pmt_ap_modifier', + ) + + pmt_circuit_load_resistor = straxen.URLConfig( + type=(int, float), + help='pmt_circuit_load_resistor', + ) + + digitizer_bits = straxen.URLConfig( + type=(int, float), + help='digitizer_bits', + ) + + digitizer_voltage_range = straxen.URLConfig( + type=(int, float), + help='digitizer_voltage_range', + ) + + gains = straxen.URLConfig( + cache=True, + help='pmt gains', + ) + photon_ap_cdfs = straxen.URLConfig( + cache=True, + help='photon_ap_cdfs', + ) + def setup(self): if self.debug: log.setLevel('DEBUG') log.debug("Running PMTAfterPulses in debug mode") - self.uniform_to_pmt_ap = straxen.get_resource(self.photon_ap_cdfs, fmt='json.gz') + self.uniform_to_pmt_ap = self.photon_ap_cdfs for k in self.uniform_to_pmt_ap.keys(): for q in self.uniform_to_pmt_ap[k].keys(): if isinstance(self.uniform_to_pmt_ap[k][q], list): self.uniform_to_pmt_ap[k][q] = np.array(self.uniform_to_pmt_ap[k][q]) - - - to_pe = straxen.get_resource(self.to_pe_file, fmt='npy') - self.to_pe = to_pe[0][1] - - adc_2_current = (self.digitizer_voltage_range - / 2 ** (self.digitizer_bits) - / self.pmt_circuit_load_resistor) - - self.gains = np.divide(adc_2_current, - self.to_pe, - out=np.zeros_like(self.to_pe), - where=self.to_pe != 0) def compute(self, S1_photons, S2_photons): @@ -184,12 +186,7 @@ def photon_afterpulse(self, merged_photon_timings, merged_photon_channels, merge _photon_amplitude = np.hstack(_photon_amplitude) _photon_gains = np.array(self.gains)[_photon_channels] * _photon_amplitude - #sorted_index = np.argsort(_photon_channels) - #return _photon_timings[sorted_index], _photon_channels[sorted_index], _photon_gains[sorted_index] return _photon_timings, _photon_channels, _photon_gains else: - return np.zeros(0, np.int64), np.zeros(0, np.int64), np.zeros(0) - - - \ No newline at end of file + return np.zeros(0, np.int64), np.zeros(0, np.int64), np.zeros(0) \ No newline at end of file diff --git a/fuse/plugins/pmt_and_daq/pmt_response_and_daq.py b/fuse/plugins/pmt_and_daq/pmt_response_and_daq.py index 4dc4221a..82611c7e 100644 --- a/fuse/plugins/pmt_and_daq/pmt_response_and_daq.py +++ b/fuse/plugins/pmt_and_daq/pmt_response_and_daq.py @@ -17,97 +17,161 @@ log = logging.getLogger('fuse.pmt_and_daq.pmt_response_and_daq') log.setLevel('WARNING') -#private_files_path = "path/to/private/files" -base_path = os.path.abspath(os.getcwd()) -private_files_path = os.path.join("/",*base_path.split("/")[:-2], "private_nt_aux_files") -config = straxen.get_resource(os.path.join(private_files_path, 'sim_files/fax_config_nt_sr0_v4.json') , fmt='json') - @export -@strax.takes_config( - strax.Option('rext', default=100000, track=False, infer_type=False, - help="right raw extension"), - strax.Option('pmt_transit_time_spread', default=config['pmt_transit_time_spread'], track=False, infer_type=False, - help="pmt_transit_time_spread"), - strax.Option('dt', default=config['sample_duration'], track=False, infer_type=False, - help="sample_duration"), - strax.Option('to_pe_file', default=os.path.join(private_files_path,"sim_files/to_pe_nt.npy"), track=False, infer_type=False, - help="to_pe file"), - strax.Option('digitizer_voltage_range', default=config['digitizer_voltage_range'], track=False, infer_type=False, - help="digitizer_voltage_range"), - strax.Option('digitizer_bits', default=config['digitizer_bits'], track=False, infer_type=False, - help="digitizer_bits"), - strax.Option('pmt_circuit_load_resistor', default=config['pmt_circuit_load_resistor'], track=False, infer_type=False, - help="pmt_circuit_load_resistor"), - strax.Option('photon_area_distribution', default=config['photon_area_distribution'], track=False, infer_type=False, - help="photon_area_distribution"), - strax.Option('samples_to_store_before', default=config['samples_to_store_before'], track=False, infer_type=False, - help="samples_to_store_before"), - strax.Option('samples_before_pulse_center', default=config['samples_before_pulse_center'], track=False, infer_type=False, - help="samples_before_pulse_center"), - strax.Option('samples_to_store_after', default=config['samples_to_store_after'], track=False, infer_type=False, - help="samples_to_store_after"), - strax.Option('samples_after_pulse_center', default=config['samples_after_pulse_center'], track=False, infer_type=False, - help="samples_after_pulse_center"), - strax.Option('pmt_pulse_time_rounding', default=config['pmt_pulse_time_rounding'], track=False, infer_type=False, - help="pmt_pulse_time_rounding"), - strax.Option('external_amplification', default=config['external_amplification'], track=False, infer_type=False, - help="external_amplification"), - strax.Option('trigger_window', default=config['trigger_window'], track=False, infer_type=False, - help="trigger_window"), - strax.Option('n_top_pmts', default=253, track=False, infer_type=False, - help="n_top_pmts"), - strax.Option('n_tpc_pmts', default=494, track=False, infer_type=False, - help="n_tpc_pmts"), - strax.Option('detector', default="XENONnT", track=False, infer_type=False, - help="detector"), - strax.Option('high_energy_deamplification_factor', default=config['high_energy_deamplification_factor'], track=False, infer_type=False, - help="high_energy_deamplification_factor"), - strax.Option('enable_noise', default=config['enable_noise'], track=False, infer_type=False, - help="enable_noise"), - strax.Option('digitizer_reference_baseline', default=config['digitizer_reference_baseline'], track=False, infer_type=False, - help="digitizer_reference_baseline"), - strax.Option('zle_threshold', default=config['zle_threshold'], track=False, infer_type=False, - help="zle_threshold"), - strax.Option('debug', default=False, track=False, infer_type=False, - help="Show debug informations"), -) class PMTResponseAndDAQ(strax.Plugin): __version__ = "0.0.0" - - depends_on = ("propagated_s2_photons", "propagated_s1_photons", "pmt_afterpulses") + + depends_on = ("photon_summary") provides = ('raw_records', 'raw_records_he', 'raw_records_aqmon')#, 'truth') data_kind = immutabledict(zip(provides, provides)) + #Config options + debug = straxen.URLConfig( + default=False, type=bool,track=False, + help='Show debug informations', + ) + + zle_threshold = straxen.URLConfig( + type=(int, float), + help='zle_threshold', + ) + + digitizer_reference_baseline = straxen.URLConfig( + type=(int, float), + help='digitizer_reference_baseline', + ) + + enable_noise = straxen.URLConfig( + type=bool, + help='enable_noise', + ) + + high_energy_deamplification_factor = straxen.URLConfig( + type=(int, float), + help='high_energy_deamplification_factor', + ) + + detector = straxen.URLConfig( + help='Detector to be simulated', + ) + + n_top_pmts = straxen.URLConfig( + type=(int), + help='Number of PMTs on top array', + ) + + n_tpc_pmts = straxen.URLConfig( + type=(int), + help='Number of PMTs in the TPC', + ) + + trigger_window = straxen.URLConfig( + type=(int, float), + help='trigger_window', + ) + + external_amplification = straxen.URLConfig( + type=(int, float), + help='external_amplification', + ) + + pmt_pulse_time_rounding = straxen.URLConfig( + type=(int, float), + help='pmt_pulse_time_rounding', + ) + + samples_after_pulse_center = straxen.URLConfig( + type=(int, float), + help='samples_after_pulse_center', + ) + + samples_to_store_after = straxen.URLConfig( + type=(int, float), + help='samples_to_store_after', + ) + + samples_before_pulse_center = straxen.URLConfig( + type=(int, float), + help='samples_before_pulse_center', + ) + + samples_to_store_before = straxen.URLConfig( + type=(int, float), + help='samples_to_store_before', + ) + + pmt_circuit_load_resistor = straxen.URLConfig( + help='pmt_circuit_load_resistor', type=(int, float), + ) + + digitizer_bits = straxen.URLConfig( + type=(int, float), + help='digitizer_bits', + ) + + digitizer_voltage_range = straxen.URLConfig( + type=(int, float), + help='digitizer_voltage_range', + ) + + dt = straxen.URLConfig( + type=(int), + help='sample_duration', + ) + + pmt_transit_time_spread = straxen.URLConfig( + type=(int, float), + help='pmt_transit_time_spread', + ) + + rext = straxen.URLConfig( + type=(int), + help='right raw extension', + ) + + pe_pulse_ts = straxen.URLConfig( + help='pe_pulse_ts', + ) + + pe_pulse_ys = straxen.URLConfig( + help='pe_pulse_ys', + ) + + gains = straxen.URLConfig( + cache=True, + help='pmt gains', + ) + + photon_area_distribution = straxen.URLConfig( + cache=True, + help='photon_area_distribution', + ) + + special_thresholds = straxen.URLConfig( + help='special_thresholds', + ) + + noise_data_tmp = straxen.URLConfig( + cache=True, + help='noise_data', + ) + + channel_map = straxen.URLConfig( + track=False, type=immutabledict, + help="immutabledict mapping subdetector to (min, max) " + "channel number.") + def setup(self): if self.debug: log.setLevel('DEBUG') log.debug("Running pmt_response_and_daq in debug mode") - - - to_pe = straxen.get_resource(self.to_pe_file, fmt='npy') - self.to_pe = to_pe[0][1] - - adc_2_current = (self.digitizer_voltage_range - / 2 ** (self.digitizer_bits) - / self.pmt_circuit_load_resistor) - - self.gains = np.divide(adc_2_current, - self.to_pe, - out=np.zeros_like(self.to_pe), - where=self.to_pe != 0) self.pmt_mask = np.array(self.gains) > 0 # Converted from to pe (from cmt by default) self.turned_off_pmts = np.arange(len(self.gains))[np.array(self.gains) == 0] - - self.photon_area_distribution = straxen.get_resource(self.photon_area_distribution, fmt='csv') - - - self.pe_pulse_ts = config["pe_pulse_ts"] - self.pe_pulse_ys = config["pe_pulse_ys"] - + #setup for Digitize Pulse Cache self.current_2_adc = self.pmt_circuit_load_resistor \ @@ -115,54 +179,35 @@ def setup(self): / (self.digitizer_voltage_range / 2 ** (self.digitizer_bits)) self.channels_bottom = np.arange(self.n_top_pmts, self.n_tpc_pmts) - #self.channel_map = dict(st.config['channel_map']) - self.channel_map = {'tpc': (0, 493), - 'he': (500, 752), - 'aqmon': (790, 807), - 'aqmon_nv': (808, 815), - 'tpc_blank': (999, 999), - 'mv': (1000, 1083), - 'aux_mv': (1084, 1087), - 'mv_blank': (1999, 1999), - 'nveto': (2000, 2119), - 'nveto_blank': (2999, 2999) - } #I can get this from the context...but how?? - + + #Ugly hack below + self.channel_map_mutable = dict(self.channel_map) + self.channel_map_mutable['sum_signal'] = 800 - self.channel_map['sum_signal'] = 800 self.channels_bottom = np.arange(self.n_top_pmts, self.n_tpc_pmts) + #We can most likely get rid of this hack and make it properly.... if self.enable_noise: - self.noise_data = straxen.get_resource(config['noise_file'], fmt='npy')['arr_0'] + self.noise_data = self.noise_data_tmp['arr_0'] #ZLE Part (Now building actual raw_records data!) self.samples_per_record = strax.DEFAULT_RECORD_LENGTH #self.blevel = 0 # buffer_filled_level self.record_buffer = np.zeros(5000000, dtype=strax.raw_record_dtype(samples_per_record=strax.DEFAULT_RECORD_LENGTH)) - - self.special_thresholds = config.get('special_thresholds', {}) #Im sure there is a better way to handle this part self._cached_pmt_current_templates = {} - def compute(self, S1_photons, S2_photons, AP_photons): + def compute(self, propagated_photons): - if len(S1_photons) == 0 and len(S2_photons) == 0 and len(AP_photons) == 0: + if len(propagated_photons) == 0: return dict(raw_records=np.zeros(0, dtype=strax.raw_record_dtype(samples_per_record=strax.DEFAULT_RECORD_LENGTH)), raw_records_he=np.zeros(0, dtype=strax.raw_record_dtype(samples_per_record=strax.DEFAULT_RECORD_LENGTH)), raw_records_aqmon=np.zeros(0, dtype=strax.raw_record_dtype(samples_per_record=strax.DEFAULT_RECORD_LENGTH))) - merged_photons = np.concatenate([S1_photons, S2_photons, AP_photons]) - S1_photons = None - S2_photons = None - AP_photons = None - - #Sort all photons by time - sortind = np.argsort(merged_photons["time"]) - merged_photons = merged_photons[sortind] #Split the photons into groups that will be simualted at once - split_photons = np.split(merged_photons, np.where(np.diff(merged_photons["time"]) > self.rext)[0]+1) + split_photons = np.split(propagated_photons, np.where(np.diff(propagated_photons["time"]) > self.rext)[0]+1) self.blevel = 0 # buffer_filled_level @@ -253,8 +298,8 @@ def compute(self, S1_photons, S2_photons, AP_photons): if self.detector == 'XENONnT': adc_wave_he = adc_wave * int(self.high_energy_deamplification_factor) if ch < self.n_top_pmts: - ch_he = np.arange(self.channel_map['he'][0], - self.channel_map['he'][1] + 1)[ch] + ch_he = np.arange(self.channel_map_mutable['he'][0], + self.channel_map_mutable['he'][1] + 1)[ch] _raw_data[ch_he, _slice] += adc_wave_he _channel_mask[ch_he] = True _channel_mask['left'][ch_he] = _channel_mask['left'][ch] @@ -263,10 +308,10 @@ def compute(self, S1_photons, S2_photons, AP_photons): sum_signal(adc_wave_he, _pulse['left'] - left, _pulse['right'] - left + 1, - _raw_data[self.channel_map['sum_signal']]) + _raw_data[self.channel_map_mutable['sum_signal']]) - _channel_mask['left'] -= left + config['trigger_window'] - _channel_mask['right'] -= left - config['trigger_window'] + _channel_mask['left'] -= left + self.trigger_window + _channel_mask['right'] -= left - self.trigger_window # Adding noise, baseline and digitizer saturation if self.enable_noise: @@ -347,20 +392,21 @@ def compute(self, S1_photons, S2_photons, AP_photons): records = self.record_buffer[:self.blevel] records = strax.sort_by_time(records) - return dict(raw_records=records[records['channel'] < self.channel_map['he'][0]], - raw_records_he=records[(records['channel'] >= self.channel_map['he'][0]) & - (records['channel'] <= self.channel_map['he'][-1])], + return dict(raw_records=records[records['channel'] < self.channel_map_mutable['he'][0]], + raw_records_he=records[(records['channel'] >= self.channel_map_mutable['he'][0]) & + (records['channel'] <= self.channel_map_mutable['he'][-1])], raw_records_aqmon=records[records['channel'] == 800], #truth=_truth ) + #Check if we can just move this function into the setup method!!! def init_pmt_current_templates(self): """ Create spe templates, for 10ns sample duration and 1ns rounding we have: _pmt_current_templates[i] : photon timing fall between [10*m+i, 10*m+i+1) (i, m are integers) """ - h = deterministic_hash(config) + h = deterministic_hash("Placeholder!") #if h in self._cached_pmt_current_templates: # _pmt_current_templates = self._cached_pmt_current_templates[h] # return _pmt_current_templates diff --git a/fuse/vertical_merger_plugin.py b/fuse/vertical_merger_plugin.py new file mode 100644 index 00000000..c354092d --- /dev/null +++ b/fuse/vertical_merger_plugin.py @@ -0,0 +1,37 @@ +import strax +from strax import Plugin, SaveWhen + +import numpy as np +from itertools import groupby + +class VerticalMergerPlugin(Plugin): + "Plugin that concatenates data from the dependencies along the fist axis" + + save_when = SaveWhen.EXPLICIT + + def infer_dtype(self): + incoming_dtypes = [self.deps[d].dtype_for(d) for d in sorted(self.depends_on)] + + eq = self.all_equal(incoming_dtypes) + if eq == False: + raise ValueError("VerticalMergerPlugin can only merge data " + "with the same dtype! " + ) + + return incoming_dtypes[0] + + + def compute(self, **kwargs): + + merged_data = np.concatenate([kwargs[x] for x in kwargs]) + + #Sort everything by time + sortind = np.argsort(merged_data["time"]) + merged_data = merged_data[sortind] + + return merged_data + + @staticmethod + def all_equal(iterable): + g = groupby(iterable) + return next(g, True) and not next(g, False) diff --git a/fuse/volume_plugin.py b/fuse/volume_plugin.py new file mode 100644 index 00000000..09d85b2f --- /dev/null +++ b/fuse/volume_plugin.py @@ -0,0 +1,40 @@ +from strax import Plugin, SaveWhen + +import numba +import numpy as np + +class VolumePlugin(Plugin): + """ + Plugin that evaluates if interactions are in a defined detector volume. + """ + + def in_ROI(self, interactions, min_z, max_z, max_r): + """ + Function that evaluates if an interaction is in the ROI. + """ + mask = in_cylinder(interactions['x'], + interactions['y'], + interactions['z'], + min_z, + max_z, + max_r + ) + return mask + +@numba.njit +def in_cylinder(x, y, z, min_z, max_z, max_r): + """ + Function which checks if a given set of coordinates is within the + boundaries of the specified cylinder. + + Args: + x,y,z: Coordinates of the interaction + min_z: Inclusive lower z boundary + max_z: Exclusive upper z boundary + max_r: Exclusive radial boundary + """ + r = np.sqrt(x ** 2 + y ** 2) + m = r < max_r + m = m & (z < max_z) + m = m & (z >= min_z) + return m \ No newline at end of file