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