Skip to content

Commit

Permalink
remove pulse_cleaning code
Browse files Browse the repository at this point in the history
  • Loading branch information
tianluyuan committed Oct 28, 2024
1 parent ac5dd4e commit 60e5548
Showing 1 changed file with 0 additions and 68 deletions.
68 changes: 0 additions & 68 deletions snowflake/library.py
Original file line number Diff line number Diff line change
Expand Up @@ -396,71 +396,3 @@ def refine_vertex_time(vertex, time, direction, pulses, omgeo):
if np.isinf(min_d):
return time
return min_t + adj_d/CCC


################## pulse cleaning
def _weighted_quantile_arg(values, weights, q=0.5):
indices = np.argsort(values)
sorted_indices = np.arange(len(values))[indices]
medianidx = (weights[indices].cumsum()/weights[indices].sum()).searchsorted(q)
if (0 <= medianidx) and (medianidx < len(values)):
return sorted_indices[medianidx]
else:
return np.nan


def weighted_quantile(values, weights, q=0.5):
if len(values) != len(weights):
raise ValueError("shape of `values` and `weights` don't match!")
index = _weighted_quantile_arg(values, weights, q=q)
if not np.isnan(index):
return values[index]
else:
return np.nan


def weighted_median(values, weights):
return weighted_quantile(values, weights, q=0.5)


def pulse_cleaning(frame, Pulses, Residual=1.5e3*I3Units.ns):
pulses = dataclasses.I3RecoPulseSeriesMap.from_frame(frame, Pulses)
mask = dataclasses.I3RecoPulseSeriesMapMask(frame, Pulses)
counter, charge = 0, 0
qtot = 0
times = dataclasses.I3TimeWindowSeriesMap()
all_ts = []
all_cs = []
for omkey, ps in pulses.items():
ts = np.asarray([p.time for p in ps])
all_ts.extend(ts)
cs = np.asarray([p.charge for p in ps])
all_cs.extend(cs)
tw_start = max(weighted_quantile(np.asarray(all_ts), np.asarray(all_cs), 0.1) - 1000, min(all_ts))
tw_stop = min(weighted_quantile(np.asarray(all_ts), np.asarray(all_cs), 0.95) + 1000, max(all_ts))
for omkey, ps in pulses.items():
ts = np.asarray([p.time for p in ps])
cs = np.asarray([p.charge for p in ps])
median = weighted_median(ts, cs)
dts = np.ediff1d(ts)
median_dts = np.median(dts)
qtot += cs.sum()
for p in ps:
if median_dts > 1200 and len(dts) > 1:
# attempt to mask out correlated noise
mask.set(omkey, p, False)
elif p.time >= min(median+Residual, tw_stop) or p.time < tw_start:
latest_time = min(median+Residual, tw_stop)
if omkey not in times:
tws = dataclasses.I3TimeWindowSeries()
tws.append(
dataclasses.I3TimeWindow(latest_time, np.inf)
) # this defines the **excluded** time window
times[omkey] = tws
mask.set(omkey, p, False)
counter += 1
charge += p.charge

frame[Pulses+"PulseCleaned"] = mask
frame[Pulses+"PulseCleanedTimeWindows"] = times
frame[Pulses+"PulseCleanedTimeRange"] = dataclasses.I3TimeWindow(tw_start, tw_stop)

0 comments on commit 60e5548

Please sign in to comment.