CSV Input

You have already learned that fuse can read Geant4 root files and simulate the microphysics effects as well as the detector response. However, it is also possible to run fuse simulations from csv input files. This can be useful if you want to simulate a specific event topology or in cases where a full Geant4 simulation is not necessary. In this notebook, we will show you how to run fuse simulations from csv input files.

Imports and config preparation

import fuse

Microphysics Simulation

At the moment fuse has two options for csv input. The first option is to give a csv file as input to the microphysics simulation similar as you would give a root file as input. The second option is to give a csv file as input to the detector simulation. In this case, the microphysics simulation is skipped and the csv file is directly given to the detector simulation.

Building some instructions

First we need to generate some instructions. This notebook provides two examples: a mono-energetic gamma source and Kr83m.

import pandas as pd
import numpy as np


def monoenergetic_source(n, energy):  # number of photons, energy of photons
    df = pd.DataFrame()

    # radius and depth are in mm for microphysics instructions
    r = np.sqrt(np.random.uniform(0, 435600, n))  # 435600 = radius of TPC in mm^2
    theta = np.random.uniform(-np.pi, np.pi, n)
    df["xp"] = r * np.cos(theta)
    df["yp"] = r * np.sin(theta)
    df["zp"] = np.random.uniform(-1500, 0, n)  # 0 is the position of the gate electrode

    df["xp_pri"] = df["xp"]
    df["yp_pri"] = df["yp"]
    df["zp_pri"] = df["zp"]

    df["ed"] = np.array([energy] * n)
    # if the times are all zero they are distributed according to source_rate by fuse
    # custom interaction times in nanoseconds can also be used e.g np.arange(n)*1e9
    df["time"] = np.zeros(n)
    df["eventid"] = np.arange(n)

    df["type"] = np.repeat("gamma", n)

    df["trackid"] = np.ones(n)
    df["parentid"] = np.zeros(n, dtype=np.int32)
    df["creaproc"] = np.repeat("None", n)
    df["parenttype"] = np.repeat("None", n)
    df["edproc"] = np.repeat("None", n)

    return df


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

    df = pd.DataFrame()

    r = np.sqrt(np.random.uniform(0, 435600, n))
    t = np.random.uniform(-np.pi, np.pi, n)
    df["xp"] = np.repeat(r * np.cos(t), 2)
    df["yp"] = np.repeat(r * np.sin(t), 2)
    df["zp"] = np.repeat(np.random.uniform(-1500, 0, n), 2)

    df["xp_pri"] = df["xp"]
    df["yp_pri"] = df["yp"]
    df["zp_pri"] = df["zp"]

    df["ed"] = np.tile(decay_energies, n)

    dt = np.random.exponential(half_life / np.log(2), n)
    df["time"] = np.array(list(zip(np.zeros(n), dt))).flatten() * 1e9  # in ns

    df["eventid"] = np.repeat(np.arange(n), 2)

    df["parenttype"] = np.tile(["Kr83[41.557]", "Kr83[9.405]"], n)

    # Not used:
    # a) since Kr83m is classified using only the parenttype
    # b) trackid, parentid are not used right now.
    # Please keep in mind that other "simulations" may require properly set
    # edproc, type and creaproc. Future epix updates using e.g. track reconstructions
    # may also need proper track- and parent-ids.
    df["trackid"] = np.tile([1, 2], n)
    df["parentid"] = np.zeros(2 * n, dtype=np.int32)
    df["creaproc"] = np.repeat("None", 2 * n)
    df["edproc"] = np.repeat("None", 2 * n)
    df["type"] = np.repeat("None", 2 * n)

    return df

We will use the mono-energetic souce here, but you can give Kr83m a try if you like. To pass the instructions to fuse we need to save them to a csv file first. The following code will generate 1000 gamma events with with an energy of 200 keV randomly distributed in the TPC.

microphysics_instructions = monoenergetic_source(1000, 200)
microphysics_instructions.to_csv("monoenergetic_200keV.csv")

Running the simulation

To use the csv file as input we just need to specify the correct path and file_name in the config. The ChunkInput plugin can handle both csv and root files as input.

When dealing with csv input files we should reduce n_interactions_per_chunk so that our chunks will not get too large later on.

st = fuse.context.xenonnt_fuse_full_chain_simulation(
    output_folder="./fuse_data",
    simulation_config="sr0_dev",
)

st.set_config(
    {
        "path": ".",
        "file_name": "monoenergetic_200keV.csv",
        "n_interactions_per_chunk": 250,
    }
)

run_number = "00000"

st.make(run_number, "microphysics_summary")

microphysics_summary = st.get_df(run_number, "microphysics_summary")

Now that we have the data loaded we can plot the number of photons per interaction against the number of electrons.

import matplotlib.pyplot as plt
from scipy.stats import gaussian_kde

fig = plt.figure()
ax = fig.gca()
ax.set_xlim(11000, 14000)
ax.set_ylim(1000, 4000)

xdata = microphysics_summary["photons"].values
ydata = microphysics_summary["electrons"].values

xy = np.vstack([xdata, ydata])
z = gaussian_kde(xy)(xy)

ax.scatter(xdata, ydata, c=z, s=7.5, cmap="viridis")

ax.set_xlabel("Number of photons]")
ax.set_ylabel("Number of electrons")

plt.show()

Detector Simulation

As already mentioned we can also run the detector physics simulation from a csv file. In this case, the microphysics simulation is skipped. To do so, we need to register a new plugin that will do the csv file handling and chunking for us but first lets create some simulation instructions.

Building the instructions

To run the detector simulation we need to specify the following information for each event: - x: the x position of the interaction - y: the y position of the interaction - z: the z position of the interaction - t: the time of the interaction - photons: the number of photons produced in the interaction - electrons: the number of electrons produced in the interaction - e_field: the electric field at the interaction position - nest_id: the nest id of the interaction - eventid: the event id of the interaction - ed: the energy deposited in the interaction - cluster_id: The number of the specific cluster. Needs to start from 1 and increase by 1 for each new cluster.

Not all of the informations are necessary to run the detector simulation. I will set them to zero in the example below.

import pandas as pd


def build_random_instructions(n):
    df = pd.DataFrame()

    r = np.sqrt(np.random.uniform(0, 2500, n))
    t = np.random.uniform(-np.pi, np.pi, n)
    df["x"] = r * np.cos(t)
    df["y"] = r * np.sin(t)
    df["z"] = np.random.uniform(-150, -1, n)

    df["photons"] = np.random.uniform(100, 5000, n)
    df["electrons"] = np.random.uniform(100, 5000, n)
    df["excitons"] = np.zeros(n)

    df["e_field"] = np.array([23] * n)
    df["nestid"] = np.array([7] * n)
    df["ed"] = np.zeros(n)

    # just set the time with respect to the start of the event
    # The events will be distributed in time by fuse
    df["t"] = np.zeros(n)

    df["eventid"] = np.arange(n)
    df["cluster_id"] = np.arange(n) + 1

    return df

detectorphysics_instructions = build_random_instructions(100)
detectorphysics_instructions.to_csv("random_detectorphysics_instructions.csv", index=False)

Run the simulation

Now that the simulation instructions are prepared we can set up our simulation context. We will use the xenonnt_fuse_full_chain_simulation and register the ChunkCsvInput plugin.

st = fuse.context.xenonnt_fuse_full_chain_simulation(
    output_folder="./fuse_data",
    simulation_config="sr0_dev",
)

st.register(fuse.detector_physics.ChunkCsvInput)

st.set_config(
    {
        "input_file": "./random_detectorphysics_instructions.csv",
        "n_interactions_per_chunk": 50,
    }
)

run_number = "00042"

st.make(run_number, "raw_records", progress_bar=True)

After we finished simulating raw_records we can process the data up to event_info.

st.make(run_number, "event_info", progress_bar=True)

event_info_data = st.get_df(run_number, "event_info")

event_info_data.head()