wireCellGenericToy/matrixGeneration.py at master · OrghoN/wireCellGenericToy · 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
138
139
140
141
142
143
144
145
146
# External Dependencies
import numpy as np
import itertools

# Internal Dependencies
from dataTypes import *
import utilities

def blobInCell(blob, cell):
    """Check if a True Blob is in a cell

    Parameters
    ----------
    blob : Blob
        The true Blob
    cell : Cell
        Reconstructed Cell

    Returns
    -------
    Boolean
        True if blob in cell, false otherwise

    """
    for planeNo in range(len(cell.wires)):
        if blob.wires[planeNo][0] < cell.wires[planeNo][0] or blob.wires[planeNo][1] > cell.wires[planeNo][1]:
            return False
    return True


def generateTrueCellMatrix(blobs, cells):
    """Generate the True charge matrix for cells

    Parameters
    ----------
    blobs : list of Blob
        List of True Blobs
    cells : list of Cell
        List of cells reconstructed with geometry

    Returns
    -------
    np.matrix
        Matrix denoting charge of each cell

    """
    matrix = list(np.zeros((len(cells), 1)))

    for blob in blobs:
        for cellNo, cell in enumerate(cells):
            if(blobInCell(blob, cell)):
                matrix[cellNo][0] += blob.charge
                break

    return np.matrix(matrix)

def constructGeometryMatrix(planes, cells):
    """Make a matrix that associates merged wires and cells based on detector Geometry

    Parameters
    ----------
    planes : list of PlaneInfo
        A list containing information for all the planes in the detector
    cells : list of Cell
        List of cells Reconstructed from Geometric information

    Returns
    -------
    list of tuple[2] int, np.matrix
        list of merged channels, Matrix that associates merged wire with merged cells

    """
    splittingList = []

    # createSplittingList
    for cellNo, cell in enumerate(cells):
        for planeNo, wire in enumerate(cell.wires):
            splittingList.append(utilities.getChannelNo(planes, wire[0], planeNo))
            splittingList.append(utilities.getChannelNo(planes, wire[1], planeNo)+1)

    #Sort and make list unique
    splittingList = list(set(splittingList))
    splittingList.sort()

    channelList = []
    matrix = []

    for cellNo, cell in enumerate(cells):
        for planeNo, wire in enumerate(cell.wires):
            channel0 = utilities.getChannelNo(planes, wire[0], planeNo)
            channel1 = utilities.getChannelNo(planes, wire[1], planeNo)

            for i in range(splittingList.index(channel0),splittingList.index(channel1+1)):
                mergedChannel = (splittingList[i],splittingList[i+1]-1)
                fractionalAssociation = (mergedChannel[1]-mergedChannel[0]+1)/(channel1-channel0+1)

                #Check if wire is already in list
                if mergedChannel not in channelList:
                    channelList.append(mergedChannel)
                    matrix.append(np.zeros(len(cells)))

                matrix[channelList.index(mergedChannel)][cellNo] = fractionalAssociation

    return channelList, np.matrix(matrix)

def constructChargeList(planes,blobs):
    chargeList = []

    for plane in planes:
        chargeList.append(list(np.zeros(plane.noOfWires,dtype=int)))

    for blob in blobs:
        wires = utilities.fireWires(planes,blob.points)
        for planeNo, plane in enumerate(wires):
            charge = blob.charge/(len(plane))
            for wireNo in plane:
                chargeList[planeNo][wireNo-1] += charge

    return chargeList

def measureCharge(wireList,chargeList):
    charges = []
    chargeList = list(itertools.chain(*chargeList))

    for wire in wireList:
        charge = 0
        test = []
        for i in range(wire[0],wire[1]+1):
            try:
                charge += chargeList[i-1]
            except:
                print(len(chargeList),i)

        charges.append(charge)

    return np.reshape(np.matrix(charges),(len(wireList),1))

def addUncertainity(geometryMatrix,wireChargeMatrix,covarianceMatrix):
    invertedcovarianceMatrix = np.linalg.inv(covarianceMatrix)
    decomposedMatrix = np.linalg.cholesky(invertedcovarianceMatrix)

    #Adding Uncertaininty through Covariance Matrix
    wireChargeMatrixU = decomposedMatrix * wireChargeMatrix
    geometryMatrixU = decomposedMatrix * geometryMatrix

    return geometryMatrixU, wireChargeMatrixU