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