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 <pkavrigin@dali-login1.rcc.local>
Co-authored-by: Shenyang Shi <ss6109@columbia.edu>
Co-authored-by: ramirezdiego <diego.ramirez@physik.uni-freiburg.de>

README.md

@@ -39,6 +39,8 @@ The other keyword arguments are:
 | `--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` |
 | `--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` |
+| `--CutNeutron` | Add if you want to filter only nuclear recoil events (maximum ER energy deposit 10 keV) | `false` |
 | `--JobNumber`  | Job number in full chain simulation. Offset is computed as `JobNumber * n_simulated_events/SourceRate`, where `n_simulated_events` is read from file. | `0` |
 | `--OutputPath`  | Output file path | Same directory as input file |
 | `--Debug`  | Tell epix if you want timing outputs | `false` |
+

bin/run_epix

@@ -19,6 +19,9 @@ def pars_args():
     parser.add_argument('--CutOnEventid', dest='cut_by_eventid',
                         action='store_true', default=False,
                         help='If true eventy start/stop acts on eventid instead of rows.'),
+    parser.add_argument('--CutNeutron', dest='nr_only',
+                        action='store_true', default=False,
+                        help='Keep NR included events with ER components smaller than 10 keV'),
     parser.add_argument('--EntryStart', dest='entry_start', type=int,
                         action='store',
                         help='First event to be read. Defaulted is zero.'),
@@ -55,8 +58,7 @@ def pars_args():
                              'in the same dir as the input file.'))
     parser.add_argument('--Debug', dest='debug',
                        action='store_true', default=False,
-                       help=('If specifed additional information is printed to the console.')
-                        )
+                       help='If specifed additional information is printed to the console.')
 
     args = parser.parse_args(sys.argv[1:])
     return args

configs/example_xenonnt.ini

@@ -19,4 +19,4 @@ xe_density = 2.862
 [GasPhase]
 to_be_stored = False
 electric_field = 200
-xe_density = 0.0176
+xe_density = 0.0176

epix/io.py

@@ -90,6 +90,7 @@ def __init__(self,
                 outer_cylinder=None,
                 kwargs={},
                 cut_by_eventid=False,
+                cut_nr_only=False,
                 ):
 
         self.directory = directory
@@ -98,6 +99,7 @@ def __init__(self,
         self.outer_cylinder = outer_cylinder
         self.kwargs = kwargs
         self.cut_by_eventid = cut_by_eventid
+        self.cut_nr_only = cut_nr_only
 
         self.file = os.path.join(self.directory, self.file_name)
 
@@ -109,10 +111,10 @@ def __init__(self,
         #Prepare cut for root and csv case
         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}))')
+                               f' & ((zp >= {self.outer_cylinder["min_z"] * 10}) & (zp < {self.outer_cylinder["max_z"] * 10}))')            
         else:
             self.cut_string = None
-    
+
     def load_file(self):
         """ 
         Function which reads a root or csv file and removes 
@@ -140,6 +142,12 @@ def load_file(self):
             m = m & m2
         interactions = interactions[m]
 
+        if self.cut_nr_only:
+            m = ((interactions['type'] == "neutron")&(interactions['edproc'] == "hadElastic")) | (interactions['edproc'] == "ionIoni")
+            e_dep_er = ak.sum(interactions[~m]['ed'], axis=1)
+            e_dep_nr = ak.sum(interactions[m]['ed'], axis=1)
+            interactions = interactions[(e_dep_er<10) & (e_dep_nr>0)]
+
         # Removing all events with no interactions:
         m = ak.num(interactions['ed']) > 0
         interactions = interactions[m]

epix/run_epix.py

@@ -28,6 +28,7 @@ def main(args, return_df=False, return_wfsim_instructions=False, strax=False):
                                         kwargs={'entry_start': args['entry_start'],
                                                 'entry_stop': args['entry_stop']},
                                         cut_by_eventid=args['cut_by_eventid'],
+                                        cut_nr_only=args['nr_only'],
                                         )
     inter, n_simulated_events = epix_file_loader.load_file()