Belle II Software development
SVDValidationTTreeCluster.py
1#!/usr/bin/env python3
2
3
10
11"""
12<header>
13 <description>
14 This module is used for the SVD validation.
15 It gets information about clusters and truehits, saving in a ttree in a ROOT file.
16 </description>
17 <noexecute>SVD validation helper class</noexecute>
18</header>
19"""
20
21import math
22
23import basf2 as b2
24
25# Some ROOT tools
26import ROOT
27from ROOT import Belle2
28from ROOT import gROOT, addressof
29from ROOT import Math
30
31# Define a ROOT struct to hold output data in the TTree
32gROOT.ProcessLine('struct EventDataCluster {\
33 int sensor_id;\
34 int layer;\
35 int ladder;\
36 int sensor;\
37 int sensor_type;\
38 int strip_dir;\
39 int cluster_truehits_number;\
40 int matched;\
41 float strip_pitch;\
42 float cluster_theta;\
43 float cluster_phi;\
44 float cluster_position;\
45 float cluster_positionSigma;\
46 float cluster_clsTime;\
47 float cluster_clsTimeSigma;\
48 float cluster_charge;\
49 float cluster_seedCharge;\
50 float cluster_size;\
51 float cluster_snr;\
52 float cluster_interstripPosition;\
53 float cluster_pull;\
54 float cluster_residual;\
55 float truehit_position;\
56 float truehit_interstripPosition;\
57 float truehit_deposEnergy;\
58 float truehit_lossmomentum;\
59 float truehit_time;\
60 float eventt0_all;\
61 float eventt0_top;\
62 float eventt0_cdc;\
63 float eventt0_svd;\
64 float eventt0_ecl;\
65};')
66
67from ROOT import EventDataCluster # noqa
68
69
70class SVDValidationTTreeCluster(b2.Module):
71 '''class to produce the ttree for cluster reconstruction validation'''
72
73 def __init__(self):
74 """Initialize the module"""
75
76 super().__init__()
77
78
79 self.file = ROOT.TFile('../SVDValidationTTreeCluster.root', 'recreate')
80
81 self.tree = ROOT.TTree('tree', 'Event data of SVD validation events')
82
83 self.data = EventDataCluster()
84
85 # Declare tree branches
86 for key in EventDataCluster.__dict__:
87 if '__' not in key:
88 formstring = '/F'
89 if isinstance(self.data.__getattribute__(key), int):
90 formstring = '/I'
91 self.tree.Branch(key, addressof(self.data, key), key + formstring)
92
93 def event(self):
94 """ Find clusters with a truehit and save needed information """
95 clusters = Belle2.PyStoreArray('SVDClusters')
96 eventt0 = Belle2.PyStoreObj('EventT0')
97
98 # get EventT0: combined, TOP, CDC, ECL
99 self.data.eventt0_all = -1
100 self.data.eventt0_top = -1
101 self.data.eventt0_cdc = -1
102 self.data.eventt0_ecl = -1
103 self.data.eventt0_svd = -1
104 top = Belle2.Const.DetectorSet(Belle2.Const.TOP)
105 cdc = Belle2.Const.DetectorSet(Belle2.Const.CDC)
106 ecl = Belle2.Const.DetectorSet(Belle2.Const.ECL)
107 svd = Belle2.Const.DetectorSet(Belle2.Const.SVD)
108 if eventt0.hasEventT0():
109 self.data.eventt0_all = eventt0.getEventT0()
110 if eventt0.hasTemporaryEventT0(svd):
111 tmp = eventt0.getTemporaryEventT0s(Belle2.Const.SVD)
112 self.data.eventt0_svd = tmp.back().eventT0
113 if eventt0.hasTemporaryEventT0(cdc):
114 tmp = eventt0.getTemporaryEventT0s(Belle2.Const.CDC)
115 self.data.eventt0_cdc = tmp.back().eventT0
116 if eventt0.hasTemporaryEventT0(top):
117 tmp = eventt0.getTemporaryEventT0s(Belle2.Const.TOP)
118 self.data.eventt0_top = tmp.back().eventT0
119 if eventt0.hasTemporaryEventT0(ecl):
120 evtT0List_ECL = eventt0.getTemporaryEventT0s(Belle2.Const.ECL)
121 # Select the event t0 value from the ECL as the one with the smallest chi squared value (defined as ".quality")
122 smallest_ECL_t0_minChi2 = evtT0List_ECL[0].quality
123 self.data.eventt0_ecl = evtT0List_ECL[0].eventT0
124 for tmp in evtT0List_ECL:
125 if tmp.quality < smallest_ECL_t0_minChi2:
126 smallest_ECL_t0_minChi2 = tmp.quality
127 self.data.eventt0_ecl = tmp.eventT0
128
129 for cluster in clusters:
130
131 cluster_truehits = cluster.getRelationsTo('SVDTrueHits')
132
133 # Sensor identification
134 sensorInfo = Belle2.VXD.GeoCache.getInstance().getSensorInfo(cluster.getSensorID())
135 sensorID = cluster.getSensorID()
136 self.data.sensor_id = int(sensorID)
137 sensorNum = sensorID.getSensorNumber()
138 self.data.sensor = sensorNum
139 layerNum = sensorID.getLayerNumber()
140 self.data.layer = layerNum
141 if (layerNum == 3):
142 sensorType = 1
143 else:
144 if (sensorNum == 1):
145 sensorType = 0
146 else:
147 sensorType = 1
148 self.data.sensor_type = sensorType
149 ladderNum = sensorID.getLadderNumber()
150 self.data.ladder = ladderNum
151 if cluster.isUCluster():
152 strip_dir = 0
153 else:
154 strip_dir = 1
155
156 self.data.matched = 1
157 # We want only clusters with at least one associated TrueHit
158 # but, to compute purity, we need to store also some information
159 # for clusters not related to true hits
160 if len(cluster_truehits) == 0:
161 self.data.matched = 0
162 # Fill tree
163 self.file.cd()
164 self.tree.Fill()
165 continue
166
167 # find the trueHit with highest energy deposit (the "best" match)
168 energy = 0
169 bestTrueHitIndex = 0
170 for i, trueHit in enumerate(cluster_truehits):
171 if trueHit.getEnergyDep() > energy:
172 energy = trueHit.getEnergyDep()
173 bestTrueHitIndex = i
174 bestTrueHit = cluster_truehits[bestTrueHitIndex]
175
176 # Cluster information
177 self.data.cluster_clsTime = cluster.getClsTime()
178 self.data.cluster_clsTimeSigma = cluster.getClsTimeSigma()
179 self.data.cluster_charge = cluster.getCharge()
180 self.data.cluster_seedCharge = cluster.getSeedCharge()
181 self.data.cluster_size = cluster.getSize()
182 self.data.cluster_snr = cluster.getSNR()
183 cluster_position = cluster.getPosition()
184 if cluster.isUCluster():
185 cluster_position = cluster.getPosition(bestTrueHit.getV())
186 # Interstrip position calculations
187 strip_pitch = sensorInfo.getUPitch(bestTrueHit.getV())
188 else:
189 strip_pitch = sensorInfo.getVPitch(bestTrueHit.getU())
190 self.data.strip_dir = strip_dir
191 self.data.strip_pitch = strip_pitch
192 self.data.cluster_interstripPosition = cluster_position % strip_pitch / strip_pitch
193 # theta and phi definitions
194 if cluster.isUCluster():
195 uPos = cluster_position
196 vPos = 0
197 else:
198 uPos = 0
199 vPos = cluster_position
200 localPosition = Math.XYZVector(uPos, vPos, 0) # sensor center at (0, 0, 0)
201 globalPosition = sensorInfo.pointToGlobal(localPosition, True)
202 x = globalPosition.X()
203 y = globalPosition.Y()
204 z = globalPosition.Z()
205 # see https://d2comp.kek.jp/record/242?ln=en for the Belle II
206 # coordinate system and related variables
207 rho = math.sqrt(x * x + y * y)
208 r = math.sqrt(x * x + y * y + z * z)
209 # get theta as arccosine(z/r)
210 thetaRadians = math.acos(z / r)
211 theta = (thetaRadians * 180) / math.pi
212 # get phi as arccosine(x/rho)
213 phiRadians = math.acos(x / rho)
214 if y < 0:
215 phi = 360 - (phiRadians * 180) / math.pi
216 else:
217 phi = (phiRadians * 180) / math.pi
218 self.data.cluster_theta = theta
219 self.data.cluster_phi = phi
220 # Pull calculations
221 clusterPos = cluster_position
222 clusterPosSigma = cluster.getPositionSigma()
223 if cluster.isUCluster():
224 truehitPos = bestTrueHit.getU()
225 else:
226 truehitPos = bestTrueHit.getV()
227 cluster_residual = clusterPos - truehitPos
228 cluster_pull = cluster_residual / clusterPosSigma
229 self.data.cluster_position = clusterPos
230 self.data.cluster_positionSigma = clusterPosSigma
231 self.data.cluster_residual = cluster_residual
232 self.data.cluster_pull = cluster_pull
233 # Truehit information
234 self.data.truehit_position = truehitPos
235 truehit_interstripPosition = truehitPos % strip_pitch / strip_pitch
236 self.data.truehit_interstripPosition = truehit_interstripPosition
237 self.data.truehit_deposEnergy = bestTrueHit.getEnergyDep()
238 self.data.truehit_lossmomentum = bestTrueHit.getEntryMomentum().R() - bestTrueHit.getExitMomentum().R()
239 self.data.truehit_time = bestTrueHit.getGlobalTime()
240
241 # DEBUG options, commented out
242 # print("\n new cluster with at least one true hit ("+str(len(cluster_truehits))+")")
243 # print("best True hit has energy = "+str(bestTrueHit.getEnergyDep()*1e6)+" eV,\
244 # global time = "+str(bestTrueHit.getGlobalTime()))
245 # print("cluster isU = "+str(cluster.isUCluster())+", size = "+str(self.data.cluster_size)+"\
246 # position = "+str(clusterPos))
247 # print("cluster charge = "+str(self.data.cluster_charge)+", time = "+str(self.data.cluster_clsTime))
248
249 # Fill tree
250 self.file.cd()
251 self.tree.Fill()
252
253 def terminate(self):
254 """Close the output file. """
255 self.file.cd()
256 self.file.Write()
257 self.file.Close()
258
Belle2::Const::DetectorSet
The DetectorSet class for sets of detector IDs in the form of EDetector values.
Definition: Const.h:80
Belle2::PyStoreArray
A (simplified) python wrapper for StoreArray.
Definition: PyStoreArray.h:72
Belle2::PyStoreObj
a (simplified) python wrapper for StoreObjPtr.
Definition: PyStoreObj.h:67
Belle2::VXD::GeoCache::getInstance
static GeoCache & getInstance()
Return a reference to the singleton instance.
Definition: GeoCache.cc:214