15 from ROOT
import Belle2
18 set_random_seed(10346)
23 Histogram the difference position residuals and pulls between clusters and truehits.
28 Create histograms for pulls and residuals
40 self.
rfile = ROOT.TFile(
"PXDPositionEstimation.root",
"RECREATE")
55 kind),
'Particle momentum kind={:d}'.format(kind), 5000, 0.0, 10.0)
57 kind),
'Particle angle thetaU kind={:d}'.format(kind), 180, -90, +90)
59 kind),
'Particle angle thetaV kind={:d}'.format(kind), 180, -90, +90)
61 kind),
'Cluster charge kind={:d}'.format(kind), 255, 0.0, 255.0)
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)
81 kind, mode),
'PXD residual pull U kind={:d} mode={:d}'.format(kind, mode), 200, -10, +10)
83 kind, mode),
'PXD residual pull V kind={:d} mode={:d}'.format(kind, mode), 200, -10, +10)
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(
98 """Fill the residual and pull histograms"""
103 for truehit
in truehits:
106 clusters = truehit.getRelationsFrom(
"PXDClusters")
109 for j, cls
in enumerate(clusters):
112 if clusters.weight(j) < 100:
115 mom = truehit.getMomentum()
118 thetaU = math.atan(tu) * 180 / math.pi
119 thetaV = math.atan(tv) * 180 / math.pi
122 for mcp
in truehit.getRelationsFrom(
"MCParticles"):
123 if not mcp.hasStatus(1):
128 clusterkind = cls.getKind()
133 if clusterkind <= 3
and mom.Mag() > 0.02:
151 pull_u = (truehit.getU() - cls.getU()) / cls.getUSigma()
152 pull_v = (truehit.getV() - cls.getV()) / cls.getVSigma()
166 shape_likelyhood = PositionEstimator.getShapeLikelyhood(cls, tu, tv)
167 if shape_likelyhood > 0:
170 offset = PositionEstimator.getClusterOffset(cls, tu, tv)
177 shiftU = sensor_info.getUCellPosition(cls.getUStart())
178 shiftV = sensor_info.getVCellPosition(cls.getVStart())
182 pull_u = (truehit.getU() - shiftU - offset.getU()) / (math.sqrt(offset.getUSigma2()))
183 pull_v = (truehit.getV() - shiftV - offset.getV()) / (math.sqrt(offset.getVSigma2()))
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())
192 truehit.getV() - shiftV - offset.getV())
195 truehit.getV() - shiftV - offset.getV())
198 truehit.getV() - shiftV - offset.getV())
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())
209 truehit.getV() - shiftV - offset.getV())
212 truehit.getV() - shiftV - offset.getV())
215 truehit.getV() - shiftV - offset.getV())
221 pull_u = (truehit.getU() - cls.getU()) / cls.getUSigma()
222 pull_v = (truehit.getV() - cls.getV()) / cls.getVSigma()
238 Format and write all histograms and plot them
241 for kind
in range(5):
258 for kind
in range(4):
259 for mode
in range(3):
281 hcoverage = ROOT.TH1F(
"hist_coverage",
'Coverage of corrections', 2, 1, 2)
284 hcoverage.SetLineWidth(2)
285 hcoverage.SetYTitle(
'coverage / %')
286 hcoverage.SetTitle(
'Coverage of cluster shape corrections')
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)))
295 if __name__ ==
"__main__":
301 parser = argparse.ArgumentParser(description=
"Test cluster shape corrections on generic BBbar events")
305 default=
'/home/benjamin/prerelease-01-00-00b-validation/samples',
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()
312 bkgfiles = glob.glob(args.bglocation +
'/*.root')
313 if len(bkgfiles) == 0:
314 print(
'No BG overlay files found')
319 main.add_module(
"EventInfoSetter", evtNumList=[10000])
320 main.add_module(
"Gearbox")
322 main.add_module(
"Geometry", components=[
'MagneticField',
'BeamPipe',
'PXD'], useDB=
False)
325 main.add_module(
"EvtGenInput")
330 bkginput = register_module(
'BGOverlayInput')
331 bkginput.param(
'inputFileNames', bkgfiles)
332 main.add_module(bkginput)
334 main.add_module(
"FullSim")
335 main.add_module(
"PXDDigitizer")
340 main.add_module(
'BGOverlayExecutor', PXDDigitsName=
'')
341 main.add_module(
"PXDDigitSorter", digits=
'')
343 main.add_module(
"ActivatePXDClusterPositionEstimator")
344 main.add_module(
"PXDClusterizer")
347 main.add_module(positionestimation)
348 main.add_module(
"Progress")