Update giovanni (#68)

* Add NR filter before quanta generation (#53) * New branch version. * NR-only cut string. * Fixed import of the NR cut config parameter * Fixed version label * nr cut with ER component smaller than 32kev * nr cut * field * change er threshold * add nr to bin run_epix * Add NR cut info to README * store_true for NR cut Co-authored-by: Pavel Kavrigin <[email protected]> Co-authored-by: Shenyang Shi <[email protected]> Co-authored-by: ramirezdiego <[email protected]> * Update HISTORY for release * Bump version: 0.2.2 → 0.3.0 * Set arg defaults (#55) * set nr_only cut default to false * add error for invalid input format * Save primary particle position (#56) * Added x_pri, y_pri, z_pri to the output * Debugged x_pri, y_pri, z_pri recording * WFSim requirements Co-authored-by: Diego Ramírez García <[email protected]> * Update HISTORY for release * Bump version: 0.3.0 → 0.3.1 * First implementation of BBF quanta generator (#57) * Added BBF yield quanta generation in epix * Fix white space * updated readme * New intermediate results as defaults * Nex/Ni ratio for NR is computed inside and does not need to be provided * update nr parameters * add back fixed parameters in BBF fitting * Correct syntax Co-authored-by: Zihao Xu <[email protected]> Co-authored-by: Diego Ramírez García <[email protected]> * Update HISTORY for release * Bump version: 0.3.1 → 0.3.2 * fix fast sim * remove prints * including Jaron ss-ms res * commit berfore pull * decision tree clustering working * Specify not working numba in req. file (#62) * fix loop for * updated * add prim position for materials and e dep from G4 * fixing bugs * updates for the meeting * a very bad example * Update requirements.txt (#63) * Bump version: 0.3.2 → 0.3.3 * Sync hystory * Update python-publish.yml * Add zenodo badge * Set S2 clustering algorithm in the config * New example notebook --------- Co-authored-by: Pavel Kavrigin <[email protected]> Co-authored-by: Pavel Kavrigin <[email protected]> Co-authored-by: Shenyang Shi <[email protected]> Co-authored-by: ramirezdiego <[email protected]> Co-authored-by: Andrii Terliuk <[email protected]> Co-authored-by: Zihao Xu <[email protected]> Co-authored-by: Jaron Grigat <[email protected]> Co-authored-by: Joran R. Angevaare <[email protected]> Co-authored-by: Joran Angevaare <[email protected]> Co-authored-by: PavelKavrigin <[email protected]> Co-authored-by: Pavel Kavrigin <[email protected]>
XENONnT · Mar 10, 2023 · f24fc4b · f24fc4b
1 parent 31cf81d
commit f24fc4b
Show file tree

Hide file tree

Showing 14 changed files with 2,366 additions and 1,113 deletions.
diff --git a/.bumpversion.cfg b/.bumpversion.cfg
@@ -1,6 +1,5 @@
 [bumpversion]
-current_version = 0.3.2
+current_version = 0.3.3
 files = setup.py epix/__init__.py
 commit = True
 tag = True
-
diff --git a/.github/workflows/python-publish.yml b/.github/workflows/python-publish.yml
@@ -19,7 +19,7 @@ jobs:
     - name: Set up Python
       uses: actions/setup-python@v2
       with:
-        python-version: '3.6'
+        python-version: '3.8'
     - name: Install dependencies
       run: |
         python -m pip install --upgrade pip

diff --git a/HISTORY.md b/HISTORY.md
@@ -1,3 +1,8 @@
+0.3.3 (2023-01-09)
+==================
+* Specify not working numba in req. file by @JoranAngevaare in https://github.com/XENONnT/epix/pull/62
+* Install scikit-learn by @JoranAngevaare in https://github.com/XENONnT/epix/pull/63
+
 0.3.2 (2022-08-15)
 ==================
  * First implementation of BBF quanta generator (#57)

diff --git a/README.md b/README.md
@@ -2,6 +2,7 @@
 
 [![PyPI version shields.io](https://img.shields.io/pypi/v/epix.svg)](https://pypi.python.org/pypi/epix/)
 [![CodeFactor](https://www.codefactor.io/repository/github/xenonnt/epix/badge)](https://www.codefactor.io/repository/github/xenonnt/epix)
+[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.7516941.svg)](https://doi.org/10.5281/zenodo.7516941)
 
 **E**lectron and **P**hoton **I**nstructions generator for **X**ENON
 

diff --git a/epix/__init__.py b/epix/__init__.py
@@ -1,4 +1,4 @@
-__version__ = '0.3.2'
+__version__ = '0.4.0'
 
 from .common import *
 from .io import *
@@ -10,4 +10,4 @@
 from .detectors import *
 from .run_epix import *
 from .quantgen_wrapper import BBF_quanta_generator
-from .simulator import *
+from .simulator import *
diff --git a/epix/io.py b/epix/io.py
@@ -110,11 +110,13 @@ def __init__(self,
                              "creaproc", "edproc"]
 
         #Prepare cut for root and csv case
+        # r and z are tranfsormted in the load_root
+        self.cut_string = '(ed > 0)'
         if self.outer_cylinder:
-            self.cut_string = (f'(r < {self.outer_cylinder["max_r"]})'
-                               f' & ((zp >= {self.outer_cylinder["min_z"] * 10}) & (zp < {self.outer_cylinder["max_z"] * 10}))')            
-        else:
-            self.cut_string = None
+            self.cut_string += (f' & (r < {self.outer_cylinder["max_r"]})'
+                                f' & ((z >= {self.outer_cylinder["min_z"]}) & (z < {self.outer_cylinder["max_z"]}))')            
+        # else:
+        #    self.cut_string = None
 
     def load_file(self):
         """ 

diff --git a/epix/simulator/GenerateEvents.py b/epix/simulator/GenerateEvents.py
@@ -32,18 +32,23 @@ def make_s1(instructions, config, resource):
         xyz = np.vstack([instructions['x'], instructions['y'], instructions['z']]).T
 
         num_ph = instructions['s1_area'].astype(np.int64)
-        n_photons = Helpers.get_s1_light_yield(n_photons=num_ph,
-                                               positions=xyz,
-                                               s1_lce_map=resource.s1_map,
-                                               config=config) * config['dpe_fraction']
+
+        # Here a WFsim function is called which remove the dpe
+        # we have to introduce it again in fast simulator
+        n_photons = Helpers.get_s1_light_yield(n_photons = num_ph,
+                                               positions = xyz,
+                                               s1_lce_map = resource.s1_map,
+                                               config = config) * (1 + config['p_double_pe_emision']) 
+
         instructions['s1_area'] = Helpers.get_s1_area_with_spe(resource.photon_area_distribution,
                                                                n_photons.astype(np.int64))
-
+                                                  
         num_ph = instructions['alt_s1_area'].astype(np.int64)
-        alt_n_photons = Helpers.get_s1_light_yield(n_photons=num_ph,
-                                                   positions=xyz,
-                                                   s1_lce_map=resource.s1_map,
-                                                   config=config) * config['dpe_fraction']
+        alt_n_photons = Helpers.get_s1_light_yield(n_photons = num_ph,
+                                                   positions = xyz,
+                                                   s1_lce_map = resource.s1_map,
+                                                   config = config) * (1 + config['p_double_pe_emision'])
+
         instructions['alt_s1_area'] = Helpers.get_s1_area_with_spe(resource.photon_area_distribution,
                                                                    alt_n_photons.astype(np.int64))
 
@@ -60,11 +65,12 @@ def make_s2(instructions, config, resource):
         alt_xy = np.array([instructions['alt_s2_x'], instructions['alt_s2_y']]).T
 
         n_el = instructions['s2_area'].astype(np.int64)
-        n_electron = Helpers.get_s2_charge_yield(n_electron=n_el,
-                                                 positions=xy,
-                                                 z_obs=instructions['z'],
-                                                 config=config,
-                                                 resource=resource)
+
+        n_electron = Helpers.get_s2_charge_yield(n_electron = n_el,
+                                                 positions = xy,
+                                                 z_obs = instructions['z'],
+                                                 config = config,
+                                                 resource = resource)
 
         n_el = instructions['alt_s2_area'].astype(np.int64)
         alt_n_electron = Helpers.get_s2_charge_yield(n_electron=n_el,
@@ -77,11 +83,13 @@ def make_s2(instructions, config, resource):
                                              config=config,
                                              resource=resource)
         sc_gain_sigma = np.sqrt(sc_gain)
-
-        instructions['s2_area'] = n_electron * np.random.normal(sc_gain, sc_gain_sigma) * config['dpe_fraction']
+
+        # Here a WFsim function is called which remove the dpe
+        # we have to introduce it again in fast simulator
+        instructions['s2_area'] = n_electron * np.random.normal(sc_gain, sc_gain_sigma) * (1 + config['p_double_pe_emision'])
         instructions['drift_time'] = -instructions['z'] / config['drift_velocity_liquid']
 
-        instructions['alt_s2_area'] = alt_n_electron * np.random.normal(sc_gain, sc_gain_sigma)
+        instructions['alt_s2_area'] = alt_n_electron * np.random.normal(sc_gain, sc_gain_sigma) * (1 + config['p_double_pe_emision'])
         instructions['alt_s2_drift_time'] = -instructions['alt_s2_z'] / config['drift_velocity_liquid']
 
     @staticmethod
@@ -96,8 +104,16 @@ def correction_s1(instructions, resource, **kwargs):
 
         event_positions = np.vstack([instructions['x'], instructions['y'], instructions['z']]).T
 
-        instructions['cs1'] = instructions['s1_area'] / resource.s1_map(event_positions)[:, 0]
-        instructions['alt_cs1'] = instructions['alt_s1_area'] / resource.s1_map(event_positions)[:, 0]
+        # Where does 0.15 come from ? Mean of s1_map in -130 < z < -20 and r < 50
+        # the_map = []
+        # for xyz in resource.s1_map.coordinate_system:
+        #     r = np.sqrt(xyz[0]**2 + xyz[1]**2)
+        #     if (-130 < xyz[2] < -20) & ( r < 50):
+        #         the_map.append(np.squeeze(resource.s1_map(xyz))[0])
+        # print(np.mean(the_map))
+
+        instructions['cs1'] = instructions['s1_area'] / (resource.s1_map(event_positions)[:, 0]/0.1581797073725071)
+        instructions['alt_cs1'] = instructions['alt_s1_area'] / (resource.s1_map(event_positions)[:, 0]/0.1581797073725071)
 
     @staticmethod
     @Helpers.assignOrder(4)
@@ -117,6 +133,7 @@ def correction_s2(instructions, config, resource):
         s2_positions = np.vstack([instructions['x'], instructions['y']]).T
         alt_s2_positions = np.vstack([instructions['alt_s2_x'], instructions['alt_s2_y']]).T
 
+        # Why S2 does not need the same threatment of S1 ?
         instructions['cs2'] = (instructions['s2_area'] * lifetime_corr / resource.s2_map(s2_positions))
         instructions['alt_cs2'] = (
                 instructions['alt_s2_area'] * alt_lifetime_corr / resource.s2_map(alt_s2_positions))

diff --git a/epix/simulator/fast_simulator.py b/epix/simulator/fast_simulator.py
@@ -5,12 +5,22 @@
 import epix
 import numpy as np
 import pandas as pd
+import os
 from copy import deepcopy
 from immutabledict import immutabledict
 from .GenerateEvents import GenerateEvents
 from .GenerateNveto import NVetoUtils
 from .helpers import Helpers
+import warnings
+import pickle
 
+# 2023-02-19: configuration_files:
+#   'nv_pmt_qe':'nveto_pmt_qe.json',
+#   'photon_area_distribution':'XENONnT_SR0_spe_distributions_20210713_no_noise_scaled.csv',
+#   's1_pattern_map': 'XENONnT_s1_xyz_patterns_LCE_MCvf051911_wires.pkl',
+#   's2_pattern_map': 'XENONnT_s2_xy_patterns_GXe_LCE_corrected_qes_MCv4.3.0_wires.pkl',
+#   's2_separation_bdt': 's2_separation_decision_tree_fast_sim.p'
+#   (FROM: /dali/lgrandi/jgrigat/s2_separation/s2_separation_decision_tree_fast_sim.p)
 
 class Simulator():
     '''Simulator class for epix to go from  epix instructions to fully processed data'''
@@ -29,6 +39,10 @@ def __init__(self, instructions_epix, config, resource):
         )
         self.instructions_epix = instructions_epix
 
+        #print(f"\n\nSimulator : Current directory [ {os.getcwd()} ]")
+        self.tree = pickle.load(open(self.config['configuration_files']['s2_separation_bdt'], 'rb+'))
+
+
     def cluster_events(self, ):
         """Events have more than 1 s1/s2. Here we throw away all of them except the largest 2
          Take the position to be that of the main s2
@@ -37,12 +51,13 @@ def cluster_events(self, ):
         First do the macro clustering. Clustered instructions will be flagged with amp=-1
         so we can safely through those out"""
 
-        Helpers.macro_cluster_events(self.instructions_epix)
-        self.instructions_epix = self.instructions_epix[self.instructions_epix['amp'] != -1]
-
-        Helpers.macro_cluster_events(self.instructions_epix)
+        Helpers.macro_cluster_events(self.tree, self.instructions_epix, self.config)
         self.instructions_epix = self.instructions_epix[self.instructions_epix['amp'] != -1]
 
+        # Why is it called for the second time??
+        #Helpers.macro_cluster_events(self.tree, self.instructions_epix, self.config)
+        #self.instructions_epix = self.instructions_epix[self.instructions_epix['amp'] != -1]
+
         event_numbers = np.unique(self.instructions_epix['event_number'])
         ins_size = len(event_numbers)
         instructions = np.zeros(ins_size, dtype=StraxSimulator.dtype['events_tpc'])
@@ -59,13 +74,19 @@ def cluster_events(self, ):
             if len(s1) < 1 or len(s2) < 1:
                 continue
 
-            i['s1_area'] = inst_s1[s1[-1]]['amp']
-            if len(s1) > 1:
-                i['alt_s1_area'] = inst_s1[s1[-2]]['amp']
+
+            # why we consider only the larger s1 ?
+            #i['s1_area'] = inst_s1[s1[-1]]['amp']
+            #if len(s1) > 1:
+                # I do not think it is correct
+            #    i['alt_s1_area'] = inst_s1[s1[-2]]['amp']
 
             i['s2_area'] = inst_s2[s2[-1]]['amp']
+            i['e_dep'] = inst_s2[s2[-1]]['e_dep']
+
             if len(s2) > 1:
                 i['alt_s2_area'] = inst_s2[s2[-2]]['amp']
+                i['alt_e_dep'] = inst_s2[s2[-2]]['e_dep']
                 i['alt_s2_x'] = inst_s2[s2[-2]]['x']
                 i['alt_s2_y'] = inst_s2[s2[-2]]['y']
                 i['alt_s2_z'] = inst_s2[s2[-2]]['z']
@@ -74,6 +95,12 @@ def cluster_events(self, ):
             i['y'] = inst_s2[s2[-1]]['y']
             i['z'] = inst_s2[s2[-1]]['z']
 
+            i['x_pri'] = inst_s2[s2[-1]]['x_pri']
+            i['y_pri'] = inst_s2[s2[-1]]['y_pri']
+            i['z_pri'] = inst_s2[s2[-1]]['z_pri']
+
+            i['g4id'] = inst_s2[s2[-1]]['g4id']
+
             # Strax wants begin and endtimes
             i['time'] = ix * 1000
             i['endtime'] = ix * 1000 + 1
@@ -89,14 +116,20 @@ def simulate(self, ):
 
     def run_simulator(self, ):
         self.cluster_events()
+
+        if isinstance(self.config['epix_config']['save_epix'], str):
+            epix_path = self.config['epix_config']['save_epix'] + self.config['epix_config']['file_name'][:-5] +'_epix_2'
+            print('Saving epix instruction: ', epix_path)
+            np.save(epix_path, self.instructions)
+
         self.simulate()
 
         return self.instructions
 
 
 # We should think a bit more to detector_config_override
 # and tell fast_sim to look into epix_args
-# also beacaus one entry of self.config is epix_config
+# also because one entry of self.config is epix_config
 @strax.takes_config(
     strax.Option('detector', default='XENONnT', help='Detector model'),
     strax.Option('detector_config_override', default='', help='For the electric field, otherwise 200 V/cm'),
@@ -105,14 +138,12 @@ def run_simulator(self, ):
     strax.Option('configuration_files', default=dict(), help='Files required for simulating'),
     strax.Option('fax_config', help='Fax configuration to load'),
     strax.Option('fax_config_overrides', help='Fax configuration to override', default=None),
-    strax.Option('xy_resolution', default=5, help='xy position resolution (cm)'),
-    strax.Option('z_resolution', default=1, help='xy position resolution (cm)'),
-    strax.Option('dpe_fraction', default=1, help="double photo electron emission probability. \
-                                                 Should be 1 since it's included in the LCE maps"),
+    strax.Option('xy_resolution', default=0, help='xy position resolution (cm)'),
+    strax.Option('z_resolution', default=0, help='xy position resolution (cm)'),
     strax.Option('nv_spe_resolution', default=0.40, help='nVeto SPE resolution'),
     strax.Option('nv_spe_res_threshold', default=0.50, help='nVeto SPE acceptance threshold'),
     strax.Option('nv_max_time_ns', default=1e7, help='nVeto maximum time for the acceptance of PMT hits in event'),
-    strax.Option('save_epix', default=False, help='save epix instruction'),
+    strax.Option('s2_clustering_algorithm', default='bdt', help='Macroclustering algorithm for S2, [ nsort | naive | bdt ]'),
 )
 class StraxSimulator(strax.Plugin):
     provides = ('events_tpc', 'events_nveto')
@@ -135,7 +166,13 @@ class StraxSimulator(strax.Plugin):
                              ('alt_s2_y', np.float),
                              ('alt_s2_z', np.float),
                              ('drift_time', np.float),
-                             ('alt_s2_drift_time', np.float)],
+                             ('alt_s2_drift_time', np.float),
+                             ('e_dep', np.float),
+                             ('alt_e_dep', np.float),
+                             ('x_pri', np.float),
+                             ('y_pri', np.float),
+                             ('z_pri', np.float),
+                             ('g4id', np.int)],
                  events_nveto=[('time', np.float),
                                ('endtime', np.float),
                                ('event_id', np.int),
@@ -148,15 +185,50 @@ def load_config(self):
         self.resource.s1_map = self.resource.s1_lce_correction_map
         self.resource.s2_map = self.resource.s2_correction_map
 
-        self.resource.nv_pmt_qe = straxen.get_resource(self.config['configuration_files']['nv_pmt_qe'], fmt='json')
-        self.resource.photon_area_distribution = epix.Helpers.average_spe_distribution(
-            get_resource(self.config['configuration_files']['photon_area_distribution'], fmt='csv'))
+        # Terrible way to handle the config
+        # we should have only one file 
+        # instead of several config
+        if 'nv_pmt_qe' in self.config['configuration_files'].keys():
+                    self.resource.nv_pmt_qe = straxen.get_resource(
+                        self.config['configuration_files']['nv_pmt_qe'], fmt='json')
+        else:
+            warnings.warn('The configuration_files should not exist!'
+                          'Everything should come by one config!'
+                          'Since nv_pmt_qe config is missing in configuration_files, '
+                          'the default one nveto_pmt_qe.json will be used')
+            self.resource.nv_pmt_qe = straxen.get_resource('nveto_pmt_qe.json', fmt='json')
+
+        if 'photon_area_distribution' in self.config['configuration_files'].keys():
+                self.resource.photon_area_distribution = epix.Helpers.average_spe_distribution(
+                                straxen.get_resource(self.config['configuration_files']['photon_area_distribution'], 
+                                fmt='csv'))
+        else:
+            warnings.warn('The configuration_files should not exist!'
+                          'Everything should come by one config!'
+                          'Since photon_area_distribution config is missing in configuration_files, '
+                          'the one in fax_config will be used')
+            self.resource.photon_area_distribution = epix.Helpers.average_spe_distribution(
+                                straxen.get_resource(self.config['photon_area_distribution'], fmt='csv'))
 
     def setup(self, ):
         overrides = self.config['fax_config_overrides']
         self.config.update(straxen.get_resource(self.config['fax_config'], fmt='json'))
         if overrides is not None:
             self.config.update(overrides)
+
+        if 'gains' not in self.config.keys():
+            warnings.warn('Why we have to provide gains here? '
+                          'No gains are passed in fax_config_overrides, default equal to 1')
+            self.config['gains'] =  [1 for i in range(494)]
+        if 'n_top_pmts' not in self.config.keys():
+            warnings.warn('This is a deault value, why we have to give it in fax_config_overrides? '
+                          'No n_top_pmts are passed in fax_config_overrides, default equal to 253')
+            self.config['n_top_pmts'] =  253
+        if 'n_tpc_pmts' not in self.config.keys():
+            warnings.warn('This is a deault value, why we have to give it in fax_config_overrides? '
+                          'No n_tpc_pmts are passed in fax_config_overrides, default equal to 494')
+            self.config['n_tpc_pmts'] =  494  
+
         self.load_config()
 
     def get_nveto_data(self, ):
@@ -185,6 +257,7 @@ def get_epix_instructions(self, ):
         epix_config['outer_cylinder'] = outer_cylinder
 
         epix_ins = epix.run_epix.main(epix_config, return_wfsim_instructions=True)
+
         if self.config['save_epix']:
             epix_path = self.config['epix_config']['path'] + self.config['epix_config']['file_name'][:-5] +'_epix'
             np.save(epix_path, epix_ins)
@@ -193,8 +266,14 @@ def get_epix_instructions(self, ):
 
     def compute(self):
         simulated_data_nveto = self.get_nveto_data()
-        epix_instructions = self.get_epix_instructions()
-        self.Simulator = Simulator(instructions_epix=epix_instructions,
+        self.epix_instructions = self.get_epix_instructions()
+
+        if isinstance(self.config['epix_config']['save_epix'], str):
+            epix_path = self.config['epix_config']['save_epix'] + self.config['epix_config']['file_name'][:-5] +'_epix_1'
+            print('Saving epix instruction: ', epix_path)
+            np.save(epix_path, self.epix_instructions)
+
+        self.Simulator = Simulator(instructions_epix=self.epix_instructions,
                                    config=self.config,
                                    resource=self.resource)
         simulated_data = self.Simulator.run_simulator()