Add the energy range based on the total energy deposit in the sensitive volume

README.md

@@ -39,6 +39,8 @@ The other keyword arguments are:
 | `--MicroSeparationTime`  | Clustering time (ns) | `10` |
 | `--TagClusterBy`  | decide if you tag the cluster (particle type, energy depositing process) according to first interaction in it (`time`) or most energetic (`energy`) | `energy` |
 | `--MaxDelay`  | Time after which we cut the rest of the event (ns) | `1e7` |
+| `--MaxEnergy`  | Upper range for energy selection _in sensitive volume_ (keV) | `inf` |
+| `--MinEnergy`  | Lower range for energy selection _in sensitive volume_ (keV) | `-1`  |
 | `--YieldModel` | Model for yield/quanta generation (nest / bbf / beta) | `nest` |
 | `--ClusterMethod` | Microclustering method (dbscan / betadecay / brem) | `dbscan` |
 | `--SourceRate`  | Event rate for event separation<br /> - `0` for no time shift (G4 time remains)<br /> - `-1` for clean time shift between events<br /> - `>0` (Hz) for random spacing | `0` |

bin/run_epix

@@ -53,6 +53,12 @@ def pars_args():
     parser.add_argument('--ClusterMethod', dest='cluster_method', type=str,
                         default='dbscan', help = 'Switch between DBSCAN, BETA-DECAY, '
                                                  'and BREMSSTRAHLUNG microclustering, [ dbscan | betadecay | brem ]')
+    parser.add_argument('--MaxEnergy', dest='max_energy', type=float,
+                        action='store', default=float('inf'), #keV
+                        help='Maximam energy to be stored [keV].')
+    parser.add_argument('--MinEnergy', dest='min_energy', type=float,
+                        action='store', default=-1, #keV
+                        help='Minimam energy to be stored [keV].')
     parser.add_argument('--JobNumber', dest='job_number', type=int,
                         action='store', default=0,
                         help='Job number in full chain simulation. Offset is computed as '
@@ -67,6 +73,7 @@ def pars_args():
                        help='If specifed additional information is printed to the console.')
 
     args = parser.parse_args(sys.argv[1:])
+    print(args)
     return args
 
 if __name__ == "__main__":

epix/common.py

@@ -179,3 +179,38 @@ def apply_time_offset(result, dt):
     if len(result) == 0:
         return result
     return result['t'][:, :] + dt
+
+def calc_tot_e_dep(df):
+    """
+    Calculate the total energy deposit in the sensitive volume
+    :param df: dataframe
+    :return: total energy deposit in each g4 event
+    """
+    #merge all energy deposits in single G4 events
+    g4id = df['g4id']
+    u, s = np.unique(g4id, return_index=True)
+    #split with event_number
+    edeps= np.split(df['e_dep'],s[1:])
+    #calculate tot_edep
+    tot_edep = np.asarray(list(map(lambda x:np.sum(x),edeps))) #sum up the e_deps
+    num_edep = np.asarray(list(map(lambda x:len(x),edeps))) #get number of e_deps
+    #max_edep = list(map(lambda x:np.max(x),edeps)) #get max e_deps
+    return np.repeat(tot_edep,num_edep)
+
+
+def apply_energy_selection(instr,e_range):
+    """
+    Apply the energy selection with total energy deposit in the sensitive volume
+    :param instr: dataframe
+    :param e_range: (minimum E, maximum E)
+    :return: dataframe after the energy selection
+    """
+    minE,maxE = e_range[0],e_range[1]
+
+    #merge all energy deposits in a single G4 event
+    instr['tot_e'] = calc_tot_e_dep(instr)/2. #devide by 2 (S1, S2)
+    instr_aft_selection = instr[(instr['tot_e']>=minE) & (instr['tot_e']<=maxE) ]
+
+    return instr_aft_selection
+
+

epix/io.py

@@ -185,20 +185,17 @@ def _load_root_file(self):
 
         # If user specified entry start/stop we have to update number of
         # events for source rate computation:
-        if self.kwargs['entry_start'] is not None:
+        if self.kwargs['entry_start'] is not None and self.cut_by_eventid:
             start = self.kwargs['entry_start']
         else:
             start = 0
-        if self.kwargs['entry_stop'] is not None:
+        if self.kwargs['entry_stop'] is not None and self.cut_by_eventid:
             stop = self.kwargs['entry_stop']
         else:
             stop = n_simulated_events
 
-        if self.cut_by_eventid:
-            # Start/stop refers to eventid so drop start drop from kwargs
-            # dict if specified, otherwise we cut again on rows.
-            self.kwargs.pop('entry_start', None)
-            self.kwargs.pop('entry_stop', None)
+        self.kwargs.pop('entry_start', None)
+        self.kwargs.pop('entry_stop', None)
 
         # Conversions and parameters to be computed:
         alias = {'x': 'xp/10',  # converting "geant4" mm to "straxen" cm
@@ -357,3 +354,6 @@ def awkward_to_wfsim_row_style(interactions):
     # Remove entries with no quanta
     res = res[res['amp'] > 0]
     return res
+
+def empty_wfsim_instructions():
+    return np.empty(0, dtype=wfsim.instruction_dtype)

epix/run_epix.py

@@ -6,7 +6,7 @@
 import warnings
 import epix
 
-from .common import ak_num, calc_dt, apply_time_offset
+from .common import ak_num, calc_dt, apply_time_offset,apply_energy_selection
 
 
 def main(args, return_df=False, return_wfsim_instructions=False, strax=False):
@@ -189,6 +189,13 @@ def main(args, return_df=False, return_wfsim_instructions=False, strax=False):
     if args['debug'] & (len(result) == 0):
         warnings.warn('No interactions left, return empty DataFrame.')
     instructions = epix.awkward_to_wfsim_row_style(result)
+    # Apply energy range selection for instructions if required
+    if (args['min_energy'] > 0.0) or (args['max_energy'] < float('inf')):
+        instructions = apply_energy_selection(instructions,[args['min_energy'],args['max_energy']])
+        if len(instructions) == 0:
+            instructions = epix.empty_wfsim_instructions()
+            if args['debug']:
+                warnings.warn('No interactions left after the energy selection, return empty DataFrame.')
     if args['source_rate'] != 0:
         # Only sort by time again if source rates were applied, otherwise
         # things are already sorted within the events and should stay this way.
@@ -212,8 +219,6 @@ def main(args, return_df=False, return_wfsim_instructions=False, strax=False):
 
     if return_wfsim_instructions:
         return instructions
-
-
 def monitor_time(prev_time, task):
     t = time.time()
     print(f'It took {(t - prev_time):.4f} sec to {task}')