FASTPROD/utils.py at master · gondiaz/FASTPROD · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
import numpy  as np

# general IC imports
NN = -999999

# IRENE
from invisible_cities.reco.peak_functions import split_in_peaks
from invisible_cities.reco.peak_functions import select_peaks
from invisible_cities.types.ic_types      import minmax

def signals_selected_splits(s1_indices, s2_indices,
                            s1_stride , s2_stride ,
                            s1_tmin, s1_tmax, s1_lmin, s1_lmax,
                            s2_tmin, s2_tmax, s2_lmin, s2_lmax):

    indices_split   = split_in_peaks(s1_indices, s1_stride)
    s1_selected_splits = select_peaks  (indices_split,
                                        minmax(min = s1_tmin, max = s1_tmax),
                                        minmax(min = s1_lmin, max = s1_lmax))

    indices_split   = split_in_peaks(s2_indices, s2_stride)
    s2_selected_splits = select_peaks  (indices_split,
                                        minmax(min = s2_tmin, max = s2_tmax),
                                        minmax(min = s2_lmin, max = s2_lmax))

    return s1_selected_splits, s2_selected_splits


def _1s1_1s2(pmt_ccwfs, s2_selected_splits, s1_selected_splits,
             s1emin   , s1wmin):
    ######## 1S1 1S2 CUT #########
    # 1S1 cut
    if len(s1_selected_splits)==0:
        return None
    # 1S2 cut
    if len(s2_selected_splits)>1:
        return None

    # S1 energy and width cut
    s1es, s1ws = [], []
    for ss in s1_selected_splits:
        s1_pmt = np.sum( pmt_ccwfs[:, ss[0]: ss[-1]], axis=0)
        s1es.append( np.sum(s1_pmt)    )
        s1ws.append( (ss[-1]-ss[0])*25 )
    s1es, s1ws = np.array(s1es), np.array(s1ws)

    sel = (s1es>=s1emin) & (s1ws>=s1wmin)
    idxs = np.argwhere(sel).flatten()

    if len(idxs)==0:
        return None
    elif len(idxs)>1:
        return None
    else:
        idx = idxs[0]
        s1_pmt = np.sum( pmt_ccwfs[:, s1_selected_splits[idx][0]: s1_selected_splits[idx][-1]], axis=0)
        times  = np.arange(s1_selected_splits[idx][0], s1_selected_splits[idx][-1])*25

        S1_time = times[np.argmax(s1_pmt)]
        return S1_time

def create_penthesilea_hits(s2_pmts_penth, s2_sipms_penth,
                            sipm_xs      , sipm_ys       , sipm_ids,
                            times        , S1_time):
    ###### create penthesilea hits ########
    n_sipms = len(sipm_ids)
    X, Y = sipm_xs[sipm_ids], sipm_ys[sipm_ids]
    T = (times - S1_time)/1000

    E_per_slice = np.sum( s2_pmts_penth, axis=0)
    hits    = []
    nn_hits = []
    for t, e, q in zip(T, E_per_slice, s2_sipms_penth.T):
        if np.sum(q)==0:
            nn_hits.append( (0, 0, t, e, NN, -1) )
        else:
            E = e * q / np.sum(q)
            hits.append( (X, Y, np.full( n_sipms, t), E, q, np.full( n_sipms, -1) ) )
    hits = np.array( hits )
    hits = np.swapaxes(hits, axis1=1, axis2=2)
    hits = np.concatenate( hits )
    H = np.array(np.zeros(np.shape(hits)[0]),
                 dtype=[("X", int)  , ("Y", int)  , ("Z", float),
                        ("E", float), ("Q", float), ("Ec",float)])
    H["X"], H["Y"], H["Z"]  = hits[:, 0], hits[:, 1], hits[:, 2]
    H["E"], H["Q"], H["Ec"] = hits[:, 3], hits[:, 4], -1

    #### remove 0 charge hits and insert NN ####
    sel = ~(H["Q"]==0)
    H = np.insert( H[sel], 0, nn_hits)
    H = np.sort( H, order="Z")

    return H


def esmeralda_charge_cut(hits, qth_esmer):
    #### Charge cut ####
    sel = (hits["Q"]>=qth_esmer)
    hits["Q"][~sel] = 0

    slides = np.unique( hits["Z"] )
    for slide in slides:
        sel = (hits["Z"]==slide)
        slide_hits = hits[sel]
        q = slide_hits["Q"]
        e = slide_hits["E"]
        if np.sum( q ) == 0:
            idxs = np.argwhere(sel).flatten()
            hits = np.delete(hits, idxs)
            hits = np.insert(hits, 0, (0, 0, slide, np.sum(e), NN, -1))
        else:
            hits["E"][sel] = np.sum( e ) * q / np.sum(q)
    sel = (hits["Q"]==0)
    hits = np.delete( hits, np.argwhere(sel))
    hits = np.sort(hits, order="Z")
    return hits


def join_NN_hits(hits):

    ###### join NN hits ######
    sel = (hits["Q"]==NN)
    nn_hits = hits[ sel]
    hits    = hits[~sel]
    slides = np.unique( hits["Z"] )
    for nn_hit in nn_hits:
        #select slide to append
        d = np.abs( slides - nn_hit["Z"] )
        slide = slides[ np.argmin( d ) ]
        slide_hits = hits[hits["Z"]==slide]
        #new energy
        new_E = np.sum(slide_hits["E"]) + nn_hit["E"]
        q = hits[hits["Z"]==slide]["Q"]
        Q = np.sum( q )
        hits["E"][hits["Z"] == slide] = new_E * q / Q

    return hits