479 Allow the saving of events with 0 S2 or 0 S2 in dorothea

#479

[author: gonzaponte]

Allow the saving of events with 0 S2 or 0 S2 in dorothea.
Useful to keep track of detection efficiencies.

[reviewer: jjgomezcadenas]

All straight forward and badly needed for immediate analysis of data
popping out for the machine.

invisible_cities/cities/base_cities.py

@@ -649,8 +649,6 @@ def __init__(self, **kwds):
 
         self.cnt.init(n_events_tot                 = 0,
                       n_empty_pmaps                = 0,
-                      n_events_not_s1              = 0,
-                      n_events_not_s2              = 0,
                       n_events_not_s2si            = 0,
                       n_events_not_s1s2_filter     = 0,
                       n_events_not_s2si_filter     = 0,
@@ -746,12 +744,6 @@ def filter_event(self, pmap):
         """Filter the event in terms of s1, s2, s2si"""
         # loop event away if any signal (s1, s2 or s2si) not present
         empty = S12SelectorOutput(False, [], [])
-        if not pmap.s1s:
-            self.cnt.n_events_not_s1 += 1
-            return empty
-        elif not pmap.s2s:
-            self.cnt.n_events_not_s2 += 1
-            return empty
 
         # filters in s12 and s2si
         f = pmap_filter(self.s1s2_selector, pmap)
@@ -784,15 +776,15 @@ def compute_xy_position(self, sr):
                       msipm          =  self.conf.msipm,
                       masked_sipm    =  self.pos_sipm_masked)
 
-    def compute_z_and_dt(self, ts2, ts1):
+    def compute_z_and_dt(self, ts2, ts1s):
         """
         Computes dt & z
         dt = ts2 - ts1 (in mus)
         z = dt * v_drift (i natural units)
 
         """
-        dt  = ts2 - ts1
-        z = dt * self.drift_v
+        dt   = ts2 - np.array(ts1s)
+        z    = dt * self.drift_v
         dt  *= units.ns / units.mus  #in mus
         return z, dt
 
@@ -848,7 +840,7 @@ def create_kr_event(self, pmapVectors, filter_output):
                     c = clusters[loc_elem_1d(cQ, c_closest)]
                     print('c_chosen = {}'.format(c))
 
-                Z, DT = self.compute_z_and_dt(evt.S2t[-1], evt.S1t[0])
+                Z, DT = self.compute_z_and_dt(evt.S2t[-1], evt.S1t)
                 Zrms  = peak.rms / units.mus
 
                 evt.Nsipm.append(c.nsipm)

invisible_cities/evm/event_model.py

@@ -396,20 +396,43 @@ def __init__(self, event_number, event_time):
         self.S2t   = [] # time
 
         self.Nsipm = [] # number of SiPMs in S2Si
-        self.DT    = [] # drift time (wrt S1[0])
+        self.DT    = [] # drift time
         self.Z     = [] # Position (x,y,z,R,phi)
-        self.Zrms  = []
         self.X     = []
         self.Y     = []
         self.R     = []
         self.Phi   = []
         self.Xrms  = [] # error in position
         self.Yrms  = []
+        self.Zrms  = []
+
+    def fill_defaults(self):
+        if self.nS1 == 0:
+            for attribute in ["w", "h", "e", "t"]:
+                setattr(self, "S1" + attribute, [np.nan])
+
+        if self.nS2 == 0:
+            for attribute in ["w", "h", "e", "t", "q"]:
+                setattr(self, "S2" + attribute, [np.nan])
+
+            self.Nsipm = 0
+            for attribute in ["X", "Y", "R", "Phi", "Xrms", "Yrms", "Zrms"]:
+                setattr(self, attribute, [np.nan])
+
+        if not self.nS1 or not self.nS2:
+            for attribute in ["Z", "DT"]:
+                setattr(self, attribute, [[np.nan] * max(self.nS1, 1)] * max(self.nS2, 1))
 
     def store(self, table):
         row = table.row
-        for i in range(int(self.nS1)):
-            for j in range(int(self.nS2)):
+
+        s1_peaks = range(int(self.nS1)) if self.nS1 else [-1]
+        s2_peaks = range(int(self.nS2)) if self.nS2 else [-1]
+        self.fill_defaults()
+
+        print(self)
+        for i in s1_peaks:
+            for j in s2_peaks:
                 row["event"  ] = self.event
                 row["time"   ] = self.time
                 row["s1_peak"] = i
@@ -429,8 +452,8 @@ def store(self, table):
                 row["S2t"    ] = self.S2t  [j]
 
                 row["Nsipm"  ] = self.Nsipm[j]
-                row["DT"     ] = self.DT   [j]
-                row["Z"      ] = self.Z    [j]
+                row["DT"     ] = self.DT   [j][i]
+                row["Z"      ] = self.Z    [j][i]
                 row["Zrms"   ] = self.Zrms [j]
                 row["X"      ] = self.X    [j]
                 row["Y"      ] = self.Y    [j]

invisible_cities/io/kdst_io_test.py

@@ -29,6 +29,16 @@ def dump_df(write, df):
                         getattr(evt, col).append(value)
                     except AttributeError:
                         setattr(evt, col, value)
+
+            evt.nS1 = int(evt.nS1)
+            evt.nS2 = int(evt.nS2)
+            for col in ("Z", "DT"):
+                column = list(getattr(evt, col))
+                setattr(evt, col, [])
+                for i in range(evt.nS2):
+                    s1_data = column[:evt.nS1 ]
+                    column  = column[ evt.nS1:]
+                    getattr(evt, col).append(s1_data)
             write(evt)
 
     with tb.open_file(filename, 'w') as h5out: