Belle II Software  release-05-01-25
test_cluster_position_estimator.py
1 #!/usr/bin/env python3
2 # -*- coding: utf-8 -*-
3 
4 
5 # This steering file steers fills cluster position residuals and pulls with truehits to
6 # test the cluster position estimator payloads from CAF.
7 #
8 # Execute as: basf2 test_cluster_position_estimator.py
9 #
10 # author: benjamin.schwenker@phys.uni-goettingen.de
11 
12 import math
13 from basf2 import *
14 import ROOT
15 from ROOT import Belle2
16 
17 # set some random seed
18 set_random_seed(10346)
19 
20 
21 class PXDPositionEstimation(Module):
22  """
23  Histogram the difference position residuals and pulls between clusters and truehits.
24  """
25 
26  def initialize(self):
27  """
28  Create histograms for pulls and residuals
29  """
30 
31 
32  self.nclusters = 0
33 
34  self.nfound_shapes = 0
35 
36  self.nfound_offset = 0
37 
38  # Let's create a root file to store all profiles
39 
40  self.rfile = ROOT.TFile("PXDPositionEstimation.root", "RECREATE")
41  self.rfile.cd()
42 
43 
45 
46  self.hist_map_theta_u = {}
47 
48  self.hist_map_theta_v = {}
49 
51 
52  # Loop over clusterkinds (=4 means 'all' kinds)
53  for kind in range(5):
54  self.hist_map_momentum[kind] = ROOT.TH1F("hist_momentum_kind_{:d}".format(
55  kind), 'Particle momentum kind={:d}'.format(kind), 5000, 0.0, 10.0)
56  self.hist_map_theta_u[kind] = ROOT.TH1F("hist_theta_u_kind_{:d}".format(
57  kind), 'Particle angle thetaU kind={:d}'.format(kind), 180, -90, +90)
58  self.hist_map_theta_v[kind] = ROOT.TH1F("hist_theta_v_kind_{:d}".format(
59  kind), 'Particle angle thetaV kind={:d}'.format(kind), 180, -90, +90)
60  self.hist_map_clustercharge[kind] = ROOT.TH1F("hist_clustercharge_kind_{:d}".format(
61  kind), 'Cluster charge kind={:d}'.format(kind), 255, 0.0, 255.0)
62 
63 
65 
67 
69 
71 
73 
74  for kind in range(4):
75  for mode in range(3):
76  self.hist_map_residual_u[(kind, mode)] = ROOT.TH1F('hist_map_residual_u_kind_{:d}_mode_{:d}'.format(
77  kind, mode), 'PXD residual U kind={:d} mode={:d}'.format(kind, mode), 400, -0.007, +0.007)
78  self.hist_map_residual_v[(kind, mode)] = ROOT.TH1F('hist_map_residual_v_kind_{:d}_mode_{:d}'.format(
79  kind, mode), 'PXD residual V kind={:d} mode={:d}'.format(kind, mode), 400, -0.007, +0.007)
80  self.hist_map_residual_pull_u[(kind, mode)] = ROOT.TH1F('hist_map_residual_pull_u_kind_{:d}_mode_{:d}'.format(
81  kind, mode), 'PXD residual pull U kind={:d} mode={:d}'.format(kind, mode), 200, -10, +10)
82  self.hist_map_residual_pull_v[(kind, mode)] = ROOT.TH1F('hist_map_residual_pull_v_kind_{:d}_mode_{:d}'.format(
83  kind, mode), 'PXD residual pull V kind={:d} mode={:d}'.format(kind, mode), 200, -10, +10)
84 
85 
86  self.binlimits = {}
87  self.binlimits[0] = (-90, -30)
88  self.binlimits[1] = (-30, +30)
89  self.binlimits[2] = (+30, +90)
90 
91  for bin in range(3):
92  name = 'hist_map_residual_v_kind_{:d}_mode_{:d}_special_{:d}'.format(kind, mode, bin)
93  title = 'PXD residual V kind={:d} mode={:d} {:.0f}<thetaV<{:.0f}'.format(
94  kind, mode, self.binlimits[bin][0], self.binlimits[bin][1])
95  self.hist_map_residual_v_special[(kind, mode, bin)] = ROOT.TH1F(name, title, 400, -0.007, +0.007)
96 
97  def event(self):
98  """Fill the residual and pull histograms"""
99 
100  # Get truehits
101  truehits = Belle2.PyStoreArray("PXDTrueHits")
102 
103  for truehit in truehits:
104  if isinstance(truehit, Belle2.PXDTrueHit):
105  sensor_info = Belle2.VXD.GeoCache.get(truehit.getSensorID())
106  clusters = truehit.getRelationsFrom("PXDClusters")
107 
108  # now check if we find a cluster
109  for j, cls in enumerate(clusters):
110  # we ignore all clusters where less then 100 electrons come from
111  # our truehit
112  if clusters.weight(j) < 100:
113  continue
114 
115  mom = truehit.getMomentum()
116  tu = mom[0] / mom[2]
117  tv = mom[1] / mom[2]
118  thetaU = math.atan(tu) * 180 / math.pi
119  thetaV = math.atan(tv) * 180 / math.pi
120 
121  # Only look at primary particles
122  for mcp in truehit.getRelationsFrom("MCParticles"):
123  if not mcp.hasStatus(1):
124  reject = True
125 
126  # Get instance of position estimator
128  clusterkind = cls.getKind()
129 
130  # Clusterkinds 0,1,2,3 refer to all cases which can currently
131  # be corrected. Cases where a cluster pixel touches a sensor
132  # edge or contains pixel with varying vPitch are excluded here.
133  if clusterkind <= 3 and mom.Mag() > 0.02:
134 
135  self.nclusters += 1
136 
137  # Fill momentum and angles for clusterkind
138  self.hist_map_momentum[clusterkind].Fill(mom.Mag())
139  self.hist_map_theta_u[clusterkind].Fill(thetaU)
140  self.hist_map_theta_v[clusterkind].Fill(thetaV)
141  self.hist_map_clustercharge[clusterkind].Fill(cls.getCharge())
142 
143  # Fill clusterkind=4 for all PXD sensors
144  self.hist_map_momentum[4].Fill(mom.Mag())
145  self.hist_map_theta_u[4].Fill(thetaU)
146  self.hist_map_theta_v[4].Fill(thetaV)
147  self.hist_map_clustercharge[4].Fill(cls.getCharge())
148 
149  # Fill the histograms (mode=2)
150  mode = 2
151  pull_u = (truehit.getU() - cls.getU()) / cls.getUSigma()
152  pull_v = (truehit.getV() - cls.getV()) / cls.getVSigma()
153 
154  self.hist_map_residual_u[(clusterkind, mode)].Fill(truehit.getU() - cls.getU())
155  self.hist_map_residual_v[(clusterkind, mode)].Fill(truehit.getV() - cls.getV())
156  self.hist_map_residual_pull_u[(clusterkind, mode)].Fill(pull_u)
157  self.hist_map_residual_pull_v[(clusterkind, mode)].Fill(pull_v)
158 
159  if thetaV >= self.binlimits[0][0] and thetaV < self.binlimits[0][1]:
160  self.hist_map_residual_v_special[(clusterkind, mode, 0)].Fill(truehit.getV() - cls.getV())
161  elif thetaV >= self.binlimits[1][0] and thetaV < self.binlimits[1][1]:
162  self.hist_map_residual_v_special[(clusterkind, mode, 1)].Fill(truehit.getV() - cls.getV())
163  else:
164  self.hist_map_residual_v_special[(clusterkind, mode, 2)].Fill(truehit.getV() - cls.getV())
165 
166  shape_likelyhood = PositionEstimator.getShapeLikelyhood(cls, tu, tv)
167  if shape_likelyhood > 0:
168  self.nfound_shapes += 1
169 
170  offset = PositionEstimator.getClusterOffset(cls, tu, tv)
171  if offset:
172  # Now, we can safely querry the correction
173  self.nfound_offset += 1
174 
175  # We need to explicitely add a shift to the offsets
176  # This is not needed when working with PXDRecoHits
177  shiftU = sensor_info.getUCellPosition(cls.getUStart())
178  shiftV = sensor_info.getVCellPosition(cls.getVStart())
179 
180  # Fill the histograms (mode=0)
181  mode = 0
182  pull_u = (truehit.getU() - shiftU - offset.getU()) / (math.sqrt(offset.getUSigma2()))
183  pull_v = (truehit.getV() - shiftV - offset.getV()) / (math.sqrt(offset.getVSigma2()))
184 
185  self.hist_map_residual_u[(clusterkind, mode)].Fill(truehit.getU() - shiftU - offset.getU())
186  self.hist_map_residual_v[(clusterkind, mode)].Fill(truehit.getV() - shiftV - offset.getV())
187  self.hist_map_residual_pull_u[(clusterkind, mode)].Fill(pull_u)
188  self.hist_map_residual_pull_v[(clusterkind, mode)].Fill(pull_v)
189 
190  if thetaV >= self.binlimits[0][0] and thetaV < self.binlimits[0][1]:
191  self.hist_map_residual_v_special[(clusterkind, mode, 0)].Fill(
192  truehit.getV() - shiftV - offset.getV())
193  elif thetaV >= self.binlimits[1][0] and thetaV < self.binlimits[1][1]:
194  self.hist_map_residual_v_special[(clusterkind, mode, 1)].Fill(
195  truehit.getV() - shiftV - offset.getV())
196  else:
197  self.hist_map_residual_v_special[(clusterkind, mode, 2)].Fill(
198  truehit.getV() - shiftV - offset.getV())
199 
200  # Fill the histograms (mode=1)
201  mode = 1
202  self.hist_map_residual_u[(clusterkind, mode)].Fill(truehit.getU() - shiftU - offset.getU())
203  self.hist_map_residual_v[(clusterkind, mode)].Fill(truehit.getV() - shiftV - offset.getV())
204  self.hist_map_residual_pull_u[(clusterkind, mode)].Fill(pull_u)
205  self.hist_map_residual_pull_v[(clusterkind, mode)].Fill(pull_v)
206 
207  if thetaV >= self.binlimits[0][0] and thetaV < self.binlimits[0][1]:
208  self.hist_map_residual_v_special[(clusterkind, mode, 0)].Fill(
209  truehit.getV() - shiftV - offset.getV())
210  elif thetaV >= self.binlimits[1][0] and thetaV < self.binlimits[1][1]:
211  self.hist_map_residual_v_special[(clusterkind, mode, 1)].Fill(
212  truehit.getV() - shiftV - offset.getV())
213  else:
214  self.hist_map_residual_v_special[(clusterkind, mode, 2)].Fill(
215  truehit.getV() - shiftV - offset.getV())
216 
217  else:
218 
219  # Fill the histograms (mode=1)
220  mode = 1
221  pull_u = (truehit.getU() - cls.getU()) / cls.getUSigma()
222  pull_v = (truehit.getV() - cls.getV()) / cls.getVSigma()
223 
224  self.hist_map_residual_u[(clusterkind, mode)].Fill(truehit.getU() - cls.getU())
225  self.hist_map_residual_v[(clusterkind, mode)].Fill(truehit.getV() - cls.getV())
226  self.hist_map_residual_pull_u[(clusterkind, mode)].Fill(pull_u)
227  self.hist_map_residual_pull_v[(clusterkind, mode)].Fill(pull_v)
228 
229  if thetaV >= self.binlimits[0][0] and thetaV < self.binlimits[0][1]:
230  self.hist_map_residual_v_special[(clusterkind, mode, 0)].Fill(truehit.getV() - cls.getV())
231  elif thetaV >= self.binlimits[1][0] and thetaV < self.binlimits[1][1]:
232  self.hist_map_residual_v_special[(clusterkind, mode, 1)].Fill(truehit.getV() - cls.getV())
233  else:
234  self.hist_map_residual_v_special[(clusterkind, mode, 2)].Fill(truehit.getV() - cls.getV())
235 
236  def terminate(self):
237  """
238  Format and write all histograms and plot them
239  """
240 
241  for kind in range(5):
242  self.hist_map_momentum[kind].SetLineWidth(2)
243  self.hist_map_momentum[kind].SetXTitle('momentum / GeV')
244  self.hist_map_momentum[kind].SetYTitle('number of particles')
245 
246  self.hist_map_theta_u[kind].SetLineWidth(2)
247  self.hist_map_theta_u[kind].SetXTitle('thetaU / degree')
248  self.hist_map_theta_u[kind].SetYTitle('number of particles')
249 
250  self.hist_map_theta_v[kind].SetLineWidth(2)
251  self.hist_map_theta_v[kind].SetXTitle('thetaV / degree')
252  self.hist_map_theta_v[kind].SetYTitle('number of particles')
253 
254  self.hist_map_clustercharge[kind].SetLineWidth(2)
255  self.hist_map_clustercharge[kind].SetXTitle('cluster charge / ADU')
256  self.hist_map_clustercharge[kind].SetYTitle('number of particles')
257 
258  for kind in range(4):
259  for mode in range(3):
260  self.hist_map_residual_pull_u[(kind, mode)].SetLineWidth(2)
261  self.hist_map_residual_pull_u[(kind, mode)].SetXTitle('pull u')
262  self.hist_map_residual_pull_u[(kind, mode)].SetYTitle('number of particles')
263 
264  self.hist_map_residual_pull_v[(kind, mode)].SetLineWidth(2)
265  self.hist_map_residual_pull_v[(kind, mode)].SetXTitle('pull v')
266  self.hist_map_residual_pull_v[(kind, mode)].SetYTitle('number of particles')
267 
268  self.hist_map_residual_u[(kind, mode)].SetLineWidth(2)
269  self.hist_map_residual_u[(kind, mode)].SetXTitle('residuals u / cm')
270  self.hist_map_residual_u[(kind, mode)].SetYTitle('number of particles')
271 
272  self.hist_map_residual_v[(kind, mode)].SetLineWidth(2)
273  self.hist_map_residual_v[(kind, mode)].SetXTitle('residuals v / cm')
274  self.hist_map_residual_v[(kind, mode)].SetYTitle('number of particles')
275 
276  for bin in range(3):
277  self.hist_map_residual_v_special[(kind, mode, bin)].SetLineWidth(2)
278  self.hist_map_residual_v_special[(kind, mode, bin)].SetXTitle('residuals v / cm')
279  self.hist_map_residual_v_special[(kind, mode, bin)].SetYTitle('number of particles')
280 
281  hcoverage = ROOT.TH1F("hist_coverage", 'Coverage of corrections', 2, 1, 2)
282  hcoverage.SetBinContent(1, 100.0 * float(self.nfound_offset / self.nclusters))
283  hcoverage.SetBinContent(2, 100.0 * float(self.nfound_shapes / self.nclusters))
284  hcoverage.SetLineWidth(2)
285  hcoverage.SetYTitle('coverage / %')
286  hcoverage.SetTitle('Coverage of cluster shape corrections')
287 
288  print("Coverage of cluster shape corrections is {:.2f}% ".format(100.0 * float(self.nfound_offset / self.nclusters)))
289  print("Coverage of cluster shape likelyhoods is {:.2f}% ".format(100.0 * float(self.nfound_shapes / self.nclusters)))
290 
291  self.rfile.Write()
292  self.rfile.Close()
293 
294 
295 if __name__ == "__main__":
296 
297  import argparse
298  import glob
299  import sys
300 
301  parser = argparse.ArgumentParser(description="Test cluster shape corrections on generic BBbar events")
302  parser.add_argument(
303  '--bglocation',
304  dest='bglocation',
305  default='/home/benjamin/prerelease-01-00-00b-validation/samples',
306  type=str,
307  help='Location of bg overlay files')
308  parser.add_argument('--bkgOverlay', dest='bkgOverlay', action="store_true", help='Perform background overlay')
309  args = parser.parse_args()
310 
311  # Find background overlay files
312  bkgfiles = glob.glob(args.bglocation + '/*.root')
313  if len(bkgfiles) == 0:
314  print('No BG overlay files found')
315  sys.exit()
316 
317  # Now let's create a path to simulate our events.
318  main = create_path()
319  main.add_module("EventInfoSetter", evtNumList=[10000])
320  main.add_module("Gearbox")
321  # We only need the pxd for this
322  main.add_module("Geometry", components=['MagneticField', 'BeamPipe', 'PXD'], useDB=False)
323 
324  # Generate BBbar events
325  main.add_module("EvtGenInput")
326 
327  # Background overlay input
328  if bkgfiles:
329  if args.bkgOverlay:
330  bkginput = register_module('BGOverlayInput')
331  bkginput.param('inputFileNames', bkgfiles)
332  main.add_module(bkginput)
333 
334  main.add_module("FullSim")
335  main.add_module("PXDDigitizer")
336 
337  # Background overlay executor - after digitizer
338  if bkgfiles:
339  if args.bkgOverlay:
340  main.add_module('BGOverlayExecutor', PXDDigitsName='')
341  main.add_module("PXDDigitSorter", digits='')
342 
343  main.add_module("ActivatePXDClusterPositionEstimator")
344  main.add_module("PXDClusterizer")
345 
346  positionestimation = PXDPositionEstimation()
347  main.add_module(positionestimation)
348  main.add_module("Progress")
349 
350  process(main)
351  print(statistics)
test_cluster_position_estimator.PXDPositionEstimation.hist_map_theta_u
hist_map_theta_u
Histograms for true particle angle thetaU.
Definition: test_cluster_position_estimator.py:46
test_cluster_position_estimator.PXDPositionEstimation.hist_map_residual_u
hist_map_residual_u
Histograms for u residuals.
Definition: test_cluster_position_estimator.py:64
Belle2::PXDTrueHit
Class PXDTrueHit - Records of tracks that either enter or leave the sensitive volume.
Definition: PXDTrueHit.h:35
Belle2::VXD::GeoCache::get
static const SensorInfoBase & get(Belle2::VxdID id)
Return a reference to the SensorInfo of a given SensorID.
Definition: GeoCache.h:141
test_cluster_position_estimator.PXDPositionEstimation.hist_map_residual_pull_v
hist_map_residual_pull_v
Histograms for v residual pulls.
Definition: test_cluster_position_estimator.py:72
test_cluster_position_estimator.PXDPositionEstimation.nfound_offset
nfound_offset
Counter for clusters where position correction was found in payload.
Definition: test_cluster_position_estimator.py:36
test_cluster_position_estimator.PXDPositionEstimation.hist_map_clustercharge
hist_map_clustercharge
Histograms for cluster charge related to particle.
Definition: test_cluster_position_estimator.py:50
test_cluster_position_estimator.PXDPositionEstimation.binlimits
binlimits
ThetaV angle ranges for v residuals.
Definition: test_cluster_position_estimator.py:86
test_cluster_position_estimator.PXDPositionEstimation.nclusters
nclusters
Counter for all clusters.
Definition: test_cluster_position_estimator.py:32
test_cluster_position_estimator.PXDPositionEstimation.rfile
rfile
Output file to store all plots.
Definition: test_cluster_position_estimator.py:40
test_cluster_position_estimator.PXDPositionEstimation.hist_map_theta_v
hist_map_theta_v
Histograms for true particle angle thetaV.
Definition: test_cluster_position_estimator.py:48
test_cluster_position_estimator.PXDPositionEstimation
Definition: test_cluster_position_estimator.py:21
test_cluster_position_estimator.PXDPositionEstimation.initialize
def initialize(self)
Definition: test_cluster_position_estimator.py:26
Belle2::PXD::PXDClusterPositionEstimator::getInstance
static PXDClusterPositionEstimator & getInstance()
Main (and only) way to access the PXDClusterPositionEstimator.
Definition: PXDClusterPositionEstimator.cc:51
test_cluster_position_estimator.PXDPositionEstimation.hist_map_residual_v
hist_map_residual_v
Histograms for v residuals.
Definition: test_cluster_position_estimator.py:66
test_cluster_position_estimator.PXDPositionEstimation.hist_map_momentum
hist_map_momentum
Histograms for true particle momenta.
Definition: test_cluster_position_estimator.py:44
test_cluster_position_estimator.PXDPositionEstimation.hist_map_residual_v_special
hist_map_residual_v_special
Histograms for v residuals for smaller thetaV range.
Definition: test_cluster_position_estimator.py:68
test_cluster_position_estimator.PXDPositionEstimation.hist_map_residual_pull_u
hist_map_residual_pull_u
Histograms for u residual pulls.
Definition: test_cluster_position_estimator.py:70
Belle2::PyStoreArray
a (simplified) python wrapper for StoreArray.
Definition: PyStoreArray.h:58
test_cluster_position_estimator.PXDPositionEstimation.event
def event(self)
Definition: test_cluster_position_estimator.py:97
test_cluster_position_estimator.PXDPositionEstimation.terminate
def terminate(self)
Definition: test_cluster_position_estimator.py:236
test_cluster_position_estimator.PXDPositionEstimation.nfound_shapes
nfound_shapes
Counter for cluster where shape likelyhood was found in payload.
Definition: test_cluster_position_estimator.py:34