Event mixer

bin/event_mixer

@@ -0,0 +1,68 @@
+#!/usr/bin/env python
+
+import os
+import argparse
+import numpy  as np
+import pandas as pd
+
+from invisible_cities.core.system_of_units import kg, dalton, year
+from invisible_cities.core.configure       import read_config_file
+from invisible_cities.evm.mixer            import Event_Mixer, _check_enough_nevents
+from invisible_cities.evm.mixer            import get_mixer_nevents, get_reco_and_sim_nevents
+
+# parse the config filename
+parser = argparse.ArgumentParser(description="Run the event mixer by parsing the configuration file")
+parser.add_argument("config", type=str, help="config file")
+args = parser.parse_args()
+config_filename = os.path.expandvars(args.config)
+
+# 2nubb half-life taken from EXO-200 Phys. Rev. C 89, 015502
+T12_2nubb = 2.165e+21 * year
+
+if __name__ == "__main__":
+
+    conf = read_config_file(config_filename)
+
+    detector_db = conf.get("detector_db")
+    inpath      = conf.get("inpath")
+    outpath     = conf.get("outpath")
+    isotopes    = conf.get("isotopes")
+    xenon_mass  = conf.get("xenon_mass")
+    enrichment  = conf.get("enrichment")
+    exposure    = conf.get("exposure")
+    T12_0nubb   = conf.get("T12_0nubb")
+    verbosity   = conf.get("verbosity", default=0)
+    ic_efficiencies  = conf.get("ic_efficiencies")
+    nevents_per_file = conf.get("nevents_per_file")
+
+    # compute initial number of 136Xe isotopes
+    assert (0. <= enrichment) & (enrichment <= 1.)
+    N0 = enrichment*(xenon_mass/(136. * dalton))
+
+    # get event df with (g4volume, isotope, nevts)
+    nevent_df = get_mixer_nevents(exposure, detector_db, isotopes)
+
+    # get number of decays of signal-like events
+    if "0nubb" in isotopes:
+        if not T12_0nubb:
+            raise Exception("0nubb half-life (T12_0nubb) is not provided")
+        nevts = N0 * (np.log(2)/T12_0nubb) * exposure
+        nevent_df.loc[len(nevent_df)] = ("ACTIVE", "0nubb", nevts)
+
+    if "2nubb" in isotopes:
+        nevts = N0 * (np.log(2)/T12_2nubb) * exposure
+        nevent_df.loc[len(nevent_df)] = ("ACTIVE", "2nubb", nevts)
+
+    # load processing efficiencies
+    eff_df = pd.read_csv(os.path.expandvars(ic_efficiencies))
+
+    # add reconstruction efficiency and poissonize
+    nevent_df.nevts = nevent_df.nevts * (eff_df.nreco/eff_df.nsim)
+    nevent_df.nevts = np.random.poisson(nevent_df.nevts)
+
+    # check input events are enough to run the mixer
+    _check_enough_nevents(nevent_df, eff_df)
+
+    # run mixer
+    mixer = Event_Mixer(inpath, outpath, nevent_df, nevents_per_file, verbosity)
+    mixer.run()

invisible_cities/config/mixer.conf

@@ -0,0 +1,58 @@
+"""
+------------
+Event mixer:
+------------
+
+:verbosity: (optional) default 0
+    verbosity level, either 0 (no verbosity) or 1
+
+:detector_db:
+    detector database name
+
+:inpath:
+  address of the isaura input files for background/signal events.
+  - The files for each isotope must be (by convention) at different locations as a function of
+  the isotope and g4volume names. See the example path below.
+  - The files must be (by convention) named such that they contain a file-number placed after the
+  first "_" in the file. Namely, the file name must start by "name_{file_number}"
+
+:isotopes:
+  isotopes to include in the mixer. Any of 0nubb, 2nubb, 214Bi, 208Tl, 40K, 60Co.
+
+:ic_efficiencies:
+  csv file with number of reconstructed (after IC processing) and simulated events for each component.
+  Numbers can be computed from get_reco_and_sim_nevents function in evm/mixer.py, and saved with
+  pandas.DataFrame.to_csv function
+
+:xenon_mass:
+  total xenon mass of the detector
+
+:enrichment:
+  136Xe enrichment factor (value between 0 and 1)
+
+:exposure:
+  total exposure of the production.
+
+:T12_0nubb:
+    0nubb half-life. If 0nubb is not included in :isotopes:, this parameter is not required
+
+:nevents_per_file:
+  number of events per output file.
+
+:out_path:
+  output filename structure as a function of a file number variable, "file_number".
+  By convention "file_number" must be placed after a first "_". See example below.
+"""
+verbosity = 1
+detector_db = "next100"
+inpath = "$ICDIR/database/test_data/mixer/{g4volume}/{isotope}/*_test.h5"
+isotopes = ["0nubb", "2nubb"]
+ic_efficiencies = "/path/to/ic_efficiencies.csv"
+
+xenon_mass = 100. * kg
+enrichment = 0.9
+exposure   = 1. * year
+T12_0nubb  = 1.e+27 * year
+
+nevents_per_file = 10
+outpath = "/tmp/mixer_{file_number}_tag.h5"

invisible_cities/core/system_of_units.py

@@ -204,11 +204,13 @@
 milligram = 1.e-3 * gram
 ton = 1.e+3 * kilogram
 kiloton = 1.e+3 * ton
+dalton = 1.660539e-27 * kilogram
 
 # symbols
 kg = kilogram
 g = gram
 mg = milligram
+Da = dalton
 
 #
 # Power [E][T^-1]