Belle II Software  release-05-01-25
PrintTrueHits.py
1 #!/usr/bin/env python3
2 # -*- coding: utf-8 -*-
3 
4 from basf2 import *
5 from ROOT import Belle2
6 
7 logging.log_level = LogLevel.WARNING
8 
9 
10 class CheckTrueHits(Module):
11 
12  """
13  Lists TrueHits with MCParticles and SimHits that generated them.
14  """
15 
16  def __init__(self):
17  """Initialize the module"""
18 
19  super(CheckTrueHits, self).__init__()
20 
21  self.truehit_stats_pxd = {
22  'through': {'total': 0, 'secondary': 0, 'remapped': 0},
23  'into': {'total': 0, 'secondary': 0, 'remapped': 0},
24  'out': {'total': 0, 'secondary': 0, 'remapped': 0},
25  'inside': {'total': 0, 'secondary': 0, 'remapped': 0}
26  }
27 
28  self.truehit_stats_svd = {
29  'through': {'total': 0, 'secondary': 0, 'remapped': 0},
30  'into': {'total': 0, 'secondary': 0, 'remapped': 0},
31  'out': {'total': 0, 'secondary': 0, 'remapped': 0},
32  'inside': {'total': 0, 'secondary': 0, 'remapped': 0}
33  }
34 
35  def initialize(self):
36  """ Does nothing """
37 
38  def beginRun(self):
39  """ Does nothing """
40 
41  def event(self):
42  """
43  List VXD TrueHits, the MCParticles that generated them, related
44  SimHits, and check the reconstruction of the mid-point parameters.
45  """
46 
47  pxd_truehits = Belle2.PyStoreArray('PXDTrueHits')
48  svd_truehits = Belle2.PyStoreArray('SVDTrueHits')
49  geocache = Belle2.VXD.GeoCache.getInstance()
50 
51  for truehit in pxd_truehits:
52  # Get the VXD and sensor thickness
53  id = truehit.getSensorID()
54  layer = id.getLayerNumber()
55  ladder = id.getLadderNumber()
56  sensor = id.getSensorNumber()
57  sensor_info = geocache.get(id)
58  thickness = sensor_info.getThickness()
59  base_info = \
60  '\nPXDTrueHit {index}: layer:{ly} ladder:{ld} sensor:{sr}'\
61  .format(index=truehit.getArrayIndex(),
62  ly=layer, ld=ladder, sr=sensor)
63  # Classify the TrueHit
64  into_type = \
65  abs(abs(truehit.getEntryW()) - thickness / 2.0) < 1.0e-6
66  out_type = \
67  abs(abs(truehit.getExitW()) - thickness / 2.0) < 1.0e-6
68  through_type = into_type and out_type
69  truehit_type_text = 'inside' # should not create TrueHit
70  if through_type:
71  truehit_type_text = 'through'
72  elif into_type:
73  truehit_type_text = 'into'
74  elif out_type:
75  truehit_type_text = 'out'
76  base_info = '\n' + 'TrueHit type: ' + truehit_type_text
77  # Get the generating MCParticle and check sign of relation weight
78  mcparticle_relations = truehit.getRelationsFrom('MCParticles')
79  weight = mcparticle_relations.weight(0)
80  particle = mcparticle_relations[0]
81  particle_type_text = 'secondary'
82  if particle.hasStatus(Belle2.MCParticle.c_PrimaryParticle):
83  particle_type_text = 'primary'
84  mcparticle_message = 'MCParticle {index}: {text}'.format(
85  index=particle.getArrayIndex(),
86  text=particle_type_text + ' ')
87  self.truehit_stats_pxd[truehit_type_text]['total'] += 1
88  if particle_type_text == 'secondary':
89  self.truehit_stats_pxd[truehit_type_text]['secondary'] += 1
90  if weight < 0:
91  self.truehit_stats_pxd[truehit_type_text]['remapped'] += 1
92  mcparticle_message += 'remapped ' + str(weight)
93  else:
94  continue
95  B2INFO(base_info + '\n' + mcparticle_message)
96  # Now get and print out the SimHits
97  simhits = truehit.getRelationsTo('PXDSimHits')
98  B2INFO('SimHits:')
99  for simhit in simhits:
100  particle_type_text = 'secondary'
101  simhit_weight = 0.0
102  particle_index = -1
103 
104  simhit_relations = simhit.getRelationsFrom('MCParticles')
105  if simhit_relations.size() == 0:
106  particle_type_text = 'none'
107  else:
108  simhit_weight = simhit_relations.weight(0)
109  particle = simhit_relations[0]
110  particle_index = particle.getArrayIndex()
111  if particle.hasStatus(Belle2.MCParticle.c_PrimaryParticle):
112  particle_type_text = 'primary'
113  simhit_text = \
114  '{index} {ptext} {pindex} {wt:9.6f}'.format(
115  index=simhit.getArrayIndex(),
116  ptext=particle_type_text,
117  pindex=particle_index,
118  wt=simhit_weight)
119  simhit_text += \
120  ' {u_in:8.4f} {v_in:8.4f} {w_in:8.4f}'.format(
121  u_in=simhit.getPosIn().X(),
122  v_in=simhit.getPosIn().Y(),
123  w_in=simhit.getPosIn().Z())
124  simhit_text += \
125  '{time:8.5f}'.format(time=simhit.getGlobalTime())
126  B2INFO(simhit_text)
127  # Print the truhit position
128  midpoint_text = 'TrueHit position:\n'
129  midpoint_text += \
130  'IN: ({u_in:8.4f},{v_in:8.4f},{w_in:8.4f})\n'.format(
131  u_in=truehit.getEntryU(),
132  v_in=truehit.getEntryV(),
133  w_in=truehit.getEntryW())
134  midpoint_text += \
135  'MID: ({u_mid:8.4f},{v_mid:8.4f},{w_mid:8.4f})\n'.format(
136  u_mid=truehit.getU(),
137  v_mid=truehit.getV(),
138  w_mid=truehit.getW())
139  midpoint_text += \
140  'OUT: ({u_out:8.4f},{v_out:8.4f},{w_out:8.4f})\n'.format(
141  u_out=truehit.getExitU(),
142  v_out=truehit.getExitV(),
143  w_out=truehit.getExitW())
144  # Print the truhit momenta
145  midpoint_text += '\nTrueHit momentum:\n'
146  midpoint_text += \
147  'IN: ({pu_in:8.4f},{pv_in:8.4f},{pw_in:8.4f})\n'.format(
148  pu_in=truehit.getEntryMomentum().X(),
149  pv_in=truehit.getEntryMomentum().Y(),
150  pw_in=truehit.getEntryMomentum().Z())
151  midpoint_text += \
152  'MID: ({pu_mid:8.4f},{pv_mid:8.4f},{pw_mid:8.4f})\n'.format(
153  pu_mid=truehit.getMomentum().X(),
154  pv_mid=truehit.getMomentum().Y(),
155  pw_mid=truehit.getMomentum().Z())
156  midpoint_text += \
157  'OUT: ({pu_out:8.4f},{pv_out:8.4f},{pw_out:8.4f})\n'.format(
158  pu_out=truehit.getExitMomentum().X(),
159  pv_out=truehit.getExitMomentum().Y(),
160  pw_out=truehit.getExitMomentum().Z())
161  # Print the truhit time
162  midpoint_text += 'Time: {t:8.5f}\n'.format(
163  t=truehit.getGlobalTime())
164  B2INFO(midpoint_text)
165 
166  for truehit in svd_truehits:
167  # Get the VXD and sensor thickness
168  id = truehit.getSensorID()
169  layer = id.getLayerNumber()
170  ladder = id.getLadderNumber()
171  sensor = id.getSensorNumber()
172  sensor_info = geocache.get(id)
173  thickness = sensor_info.getThickness()
174  base_info = \
175  '\nSVDTrueHit {index}: layer:{ly} ladder:{ld} sensor:{sr}'\
176  .format(index=truehit.getArrayIndex(),
177  ly=layer, ld=ladder, sr=sensor)
178  # Classify the TrueHit
179  into_type = \
180  abs(abs(truehit.getEntryW()) - thickness / 2.0) < 1.0e-6
181  out_type = \
182  abs(abs(truehit.getExitW()) - thickness / 2.0) < 1.0e-6
183  through_type = into_type and out_type
184  truehit_type_text = 'inside' # should not create TrueHit
185  if through_type:
186  truehit_type_text = 'through'
187  elif into_type:
188  truehit_type_text = 'into'
189  elif out_type:
190  truehit_type_text = 'out'
191  base_info += '\n' + 'TrueHit type: ' + truehit_type_text
192  # Get the generating MCParticle and check sign of relation weight
193  mcparticle_relations = truehit.getRelationsFrom('MCParticles')
194  weight = mcparticle_relations.weight(0)
195  particle = mcparticle_relations[0]
196  particle_type_text = 'secondary'
197  if particle.hasStatus(Belle2.MCParticle.c_PrimaryParticle):
198  particle_type_text = 'primary'
199  mcparticle_message = 'MCParticle {index}: {text}'.format(
200  index=particle.getArrayIndex(),
201  text=particle_type_text + ' ')
202  self.truehit_stats_svd[truehit_type_text]['total'] += 1
203  if particle_type_text == 'secondary':
204  self.truehit_stats_svd[truehit_type_text]['secondary'] += 1
205  if weight < 0:
206  self.truehit_stats_svd[truehit_type_text]['remapped'] += 1
207  mcparticle_message += 'remapped ' + str(weight)
208  else:
209  continue
210  B2INFO(base_info + '\n' + mcparticle_message)
211  # Now get and print out the SimHits
212  simhits = truehit.getRelationsTo('SVDSimHits')
213  B2INFO('SimHits:')
214  for simhit in simhits:
215  particle_type_text = 'secondary'
216  simhit_weight = 0.0
217  particle_index = -1
218 
219  simhit_relations = simhit.getRelationsFrom('MCParticles')
220  if simhit_relations.size() == 0:
221  particle_type_text = 'none'
222  else:
223  simhit_weight = simhit_relations.weight(0)
224  particle = simhit_relations[0]
225  particle_index = particle.getArrayIndex()
226  if particle.hasStatus(Belle2.MCParticle.c_PrimaryParticle):
227  particle_type_text = 'primary'
228  simhit_text = \
229  '{index} {ptext} {pindex} {wt:9.6f}'.format(
230  index=simhit.getArrayIndex(),
231  ptext=particle_type_text,
232  pindex=particle_index,
233  wt=simhit_weight)
234  simhit_text += \
235  ' {u_in:8.4f} {v_in:8.4f} {w_in:8.4f}'.format(
236  u_in=simhit.getPosIn().X(),
237  v_in=simhit.getPosIn().Y(),
238  w_in=simhit.getPosIn().Z())
239  simhit_text += \
240  '{time:8.5f} '.format(time=simhit.getGlobalTime())
241  B2INFO(simhit_text)
242  # Print the truhit position
243  midpoint_text = 'TrueHit position:\n'
244  midpoint_text += \
245  'IN: ({u_in:8.4f},{v_in:8.4f},{w_in:8.4f})\n'.format(
246  u_in=truehit.getEntryU(),
247  v_in=truehit.getEntryV(),
248  w_in=truehit.getEntryW())
249  midpoint_text += \
250  'MID: ({u_mid:8.4f},{v_mid:8.4f},{w_mid:8.4f})\n'.format(
251  u_mid=truehit.getU(),
252  v_mid=truehit.getV(),
253  w_mid=truehit.getW())
254  midpoint_text += \
255  'OUT: ({u_out:8.4f},{v_out:8.4f},{w_out:8.4f})\n'.format(
256  u_out=truehit.getExitU(),
257  v_out=truehit.getExitV(),
258  w_out=truehit.getExitW())
259  # Print the truhit momenta
260  midpoint_text += '\nTrueHit momentum:\n'
261  midpoint_text += \
262  'IN: ({pu_in:8.4f},{pv_in:8.4f},{pw_in:8.4f})\n'.format(
263  pu_in=truehit.getEntryMomentum().X(),
264  pv_in=truehit.getEntryMomentum().Y(),
265  pw_in=truehit.getEntryMomentum().Z())
266  midpoint_text += \
267  'MID: ({pu_mid:8.4f},{pv_mid:8.4f},{pw_mid:8.4f})\n'.format(
268  pu_mid=truehit.getMomentum().X(),
269  pv_mid=truehit.getMomentum().Y(),
270  pw_mid=truehit.getMomentum().Z())
271  midpoint_text += \
272  'OUT: ({pu_out:8.4f},{pv_out:8.4f},{pw_out:8.4f})\n'.format(
273  pu_out=truehit.getExitMomentum().X(),
274  pv_out=truehit.getExitMomentum().Y(),
275  pw_out=truehit.getExitMomentum().Z())
276  # Print the truhit time
277  midpoint_text += 'Time: {t:8.5f}\n'.format(
278  t=truehit.getGlobalTime())
279  B2INFO(midpoint_text)
280 
281  def terminate(self):
282  """ Write results """
283  B2INFO(
284  '\nStatistics for PXD: {pxd};\nStatistics for SVD: {svd}\n'.format(
285  pxd=str(self.truehit_stats_pxd),
286  svd=str(self.truehit_stats_svd)
287  )
288  )
289 
290 # Particle gun module
291 particlegun = register_module('ParticleGun')
292 # Create Event information
293 eventinfosetter = register_module('EventInfoSetter')
294 # Show progress of processing
295 progress = register_module('Progress')
296 # Load parameters
297 gearbox = register_module('Gearbox')
298 # Create geometry
299 geometry = register_module('Geometry')
300 # Run simulation
301 simulation = register_module('FullSim')
302 # simulation.param('StoreAllSecondaries', True)
303 # PXD digitization module
304 printParticles = CheckTrueHits()
305 printParticles.set_log_level(LogLevel.INFO)
306 
307 # Specify number of events to generate
308 eventinfosetter.param({'evtNumList': [200], 'runList': [1]})
309 
310 # Set parameters for particlegun
311 particlegun.param({
312  'nTracks': 1,
313  'varyNTracks': True,
314  'pdgCodes': [211, -211, 11, -11],
315  'momentumGeneration': 'normalPt',
316  'momentumParams': [2, 1],
317  'phiGeneration': 'normal',
318  'phiParams': [0, 360],
319  'thetaGeneration': 'uniformCos',
320  'thetaParams': [17, 150],
321  'vertexGeneration': 'normal',
322  'xVertexParams': [0, 1],
323  'yVertexParams': [0, 1],
324  'zVertexParams': [0, 1],
325  'independentVertices': False,
326 })
327 
328 # Select subdetectors to be built
329 geometry.param('components', ['MagneticField', 'PXD', 'SVD'])
330 
331 # create processing path
332 main = create_path()
333 main.add_module(eventinfosetter)
334 main.add_module(progress)
335 main.add_module(particlegun)
336 main.add_module(gearbox)
337 main.add_module(geometry)
338 main.add_module(simulation)
339 main.add_module(printParticles)
340 
341 # generate events
342 process(main)
343 
344 # show call statistics
345 print(statistics)
PrintTrueHits.CheckTrueHits.terminate
def terminate(self)
Definition: PrintTrueHits.py:281
PrintTrueHits.CheckTrueHits.truehit_stats_pxd
truehit_stats_pxd
TrueHit statistics for the PXD.
Definition: PrintTrueHits.py:21
PrintTrueHits.CheckTrueHits.__init__
def __init__(self)
Definition: PrintTrueHits.py:16
Belle2::VXD::GeoCache::getInstance
static GeoCache & getInstance()
Return a reference to the singleton instance.
Definition: GeoCache.cc:215
PrintTrueHits.CheckTrueHits
Definition: PrintTrueHits.py:10
PrintTrueHits.CheckTrueHits.beginRun
def beginRun(self)
Definition: PrintTrueHits.py:38
PrintTrueHits.CheckTrueHits.truehit_stats_svd
truehit_stats_svd
Truehit statistics for the SVD.
Definition: PrintTrueHits.py:28
PrintTrueHits.CheckTrueHits.initialize
def initialize(self)
Definition: PrintTrueHits.py:35
Belle2::PyStoreArray
a (simplified) python wrapper for StoreArray.
Definition: PyStoreArray.h:58
PrintTrueHits.CheckTrueHits.event
def event(self)
Definition: PrintTrueHits.py:41