Belle II Software  release-06-01-15
SVDValidationTTreeCluster.py
1 #!/usr/bin/env python3
2 # -*- coding: utf-8 -*-
3 
4 
11 
12 """
13 <header>
14  <contact> SVD Software Group, svd-software@belle2.org </contact>
15  <description>
16  This module is used for the SVD validation.
17  It gets information about clusters and truehits, saving in a ttree in a ROOT file.
18  </description>
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  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_ecl;\
64 };')
65 
66 
67 from ROOT import EventDataCluster # noqa
68 
69 
70 class SVDValidationTTreeCluster(b2.Module):
71  '''class to produce the ttree for cluster validation'''
72 
73  def __init__(self):
74  """Initialize the module"""
75 
76  super(SVDValidationTTreeCluster, self).__init__()
77 
78 
79  self.filefile = ROOT.TFile('../SVDValidationTTreeCluster.root', 'recreate')
80 
81  self.treetree = ROOT.TTree('tree', 'Event data of SVD validation events')
82 
83  self.datadata = EventDataCluster()
84 
85  # Declare tree branches
86  for key in EventDataCluster.__dict__:
87  if '__' not in key:
88  formstring = '/F'
89  if isinstance(self.datadata.__getattribute__(key), int):
90  formstring = '/I'
91  self.treetree.Branch(key, addressof(self.datadata, 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.datadata.eventt0_all = -1
100  self.datadata.eventt0_top = -1
101  self.datadata.eventt0_cdc = -1
102  self.datadata.eventt0_ecl = -1
103  top = Belle2.Const.DetectorSet(Belle2.Const.TOP)
104  cdc = Belle2.Const.DetectorSet(Belle2.Const.CDC)
105  ecl = Belle2.Const.DetectorSet(Belle2.Const.ECL)
106  if eventt0.hasEventT0():
107  self.datadata.eventt0_all = eventt0.getEventT0()
108  if eventt0.hasTemporaryEventT0(cdc):
109  tmp = eventt0.getTemporaryEventT0s(Belle2.Const.CDC)
110  self.datadata.eventt0_cdc = tmp.back().eventT0
111  if eventt0.hasTemporaryEventT0(top):
112  tmp = eventt0.getTemporaryEventT0s(Belle2.Const.TOP)
113  self.datadata.eventt0_top = tmp.back().eventT0
114  if eventt0.hasTemporaryEventT0(ecl):
115  evtT0List_ECL = eventt0.getTemporaryEventT0s(Belle2.Const.ECL)
116  # Select the event t0 value from the ECL as the one with the smallest chi squared value (defined as ".quality")
117  smallest_ECL_t0_minChi2 = evtT0List_ECL[0].quality
118  self.datadata.eventt0_ecl = evtT0List_ECL[0].eventT0
119  for tmp in evtT0List_ECL:
120  if tmp.quality < smallest_ECL_t0_minChi2:
121  smallest_ECL_t0_minChi2 = tmp.quality
122  self.datadata.eventt0_ecl = tmp.eventT0
123 
124  for cluster in clusters:
125  cluster_truehits = cluster.getRelationsTo('SVDTrueHits') # SVDClustersToSVDTrueHits
126  cluster_TrueHit_Length = len(cluster_truehits)
127  if (cluster_TrueHit_Length == 0) or (cluster_TrueHit_Length != 1):
128  # Sensor identification
129  sensorID = cluster.getSensorID()
130  self.datadata.sensor_id = int(sensorID)
131  sensorNum = sensorID.getSensorNumber()
132  self.datadata.sensor = sensorNum
133  layerNum = sensorID.getLayerNumber()
134  self.datadata.layer = layerNum
135  if (layerNum == 3):
136  sensorType = 1
137  else:
138  if (sensorNum == 1):
139  sensorType = 0
140  else:
141  sensorType = 1
142  self.datadata.sensor_type = sensorType
143  ladderNum = sensorID.getLadderNumber()
144  self.datadata.ladder = ladderNum
145  #
146  self.datadata.strip_dir = -1
147  self.datadata.cluster_truehits_number = cluster_TrueHit_Length
148  # Fill tree
149  self.filefile.cd()
150  self.treetree.Fill()
151  else:
152  # We want only clusters with exactly one associated TrueHit
153  for truehit in cluster_truehits:
154  self.datadata.cluster_truehits_number = cluster_TrueHit_Length
155  sensorInfo = Belle2.VXD.GeoCache.get(cluster.getSensorID())
156  # Sensor identification
157  sensorID = cluster.getSensorID()
158  self.datadata.sensor_id = int(sensorID)
159  sensorNum = sensorID.getSensorNumber()
160  self.datadata.sensor = sensorNum
161  layerNum = sensorID.getLayerNumber()
162  self.datadata.layer = layerNum
163  if (layerNum == 3):
164  sensorType = 1
165  else:
166  if (sensorNum == 1):
167  sensorType = 0
168  else:
169  sensorType = 1
170  self.datadata.sensor_type = sensorType
171  ladderNum = sensorID.getLadderNumber()
172  self.datadata.ladder = ladderNum
173  # Cluster information
174  self.datadata.cluster_clsTime = cluster.getClsTime()
175  self.datadata.cluster_clsTimeSigma = cluster.getClsTimeSigma()
176  self.datadata.cluster_charge = cluster.getCharge()
177  self.datadata.cluster_seedCharge = cluster.getSeedCharge()
178  self.datadata.cluster_size = cluster.getSize()
179  self.datadata.cluster_snr = cluster.getSNR()
180  cluster_position = cluster.getPosition()
181  if cluster.isUCluster():
182  cluster_position = cluster.getPosition(truehit.getV())
183  # Interstrip position calculations
184  if cluster.isUCluster():
185  strip_dir = 0
186  strip_pitch = sensorInfo.getUPitch(truehit.getV())
187  else:
188  strip_dir = 1
189  strip_pitch = sensorInfo.getVPitch(truehit.getU())
190  self.datadata.strip_dir = strip_dir
191  self.datadata.strip_pitch = strip_pitch
192  cluster_interstripPosition = cluster_position % strip_pitch / strip_pitch
193  self.datadata.cluster_interstripPosition = cluster_interstripPosition
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 = TVector3(uPos, vPos, 0) # sensor center at (0, 0, 0)
202  globalPosition = sensorInfo.pointToGlobal(localPosition, True)
203  x = globalPosition[0]
204  y = globalPosition[1]
205  z = globalPosition[2]
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.datadata.cluster_theta = theta
220  self.datadata.cluster_phi = phi
221  # Pull calculations
222  clusterPos = cluster_position
223  clusterPosSigma = cluster.getPositionSigma()
224  if cluster.isUCluster():
225  truehitPos = truehit.getU()
226  else:
227  truehitPos = truehit.getV()
228  cluster_residual = clusterPos - truehitPos
229  cluster_pull = cluster_residual / clusterPosSigma
230  self.datadata.cluster_position = clusterPos
231  self.datadata.cluster_positionSigma = clusterPosSigma
232  self.datadata.cluster_residual = cluster_residual
233  self.datadata.cluster_pull = cluster_pull
234  # Truehit information
235  self.datadata.truehit_position = truehitPos
236  truehit_interstripPosition = truehitPos % strip_pitch / strip_pitch
237  self.datadata.truehit_interstripPosition = truehit_interstripPosition
238  self.datadata.truehit_deposEnergy = truehit.getEnergyDep()
239  self.datadata.truehit_lossmomentum = truehit.getEntryMomentum().Mag() - truehit.getExitMomentum().Mag()
240  self.datadata.truehit_time = truehit.getGlobalTime()
241 
242  # Fill tree
243  self.filefile.cd()
244  self.treetree.Fill()
245 
246  def terminate(self):
247  """Close the output file. """
248  self.filefile.cd()
249  self.filefile.Write()
250  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:56
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