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