Belle II Software  release-05-02-19
record.py
1 #!/usr/bin/env python3
2 # -*- coding: utf-8 -*-
3 
4 import basf2
5 
6 import ROOT
7 from ROOT import Belle2 # make Belle2 namespace available
8 from ROOT import std
9 
10 import os
11 import sys
12 
13 from tracking.run.utilities import NonstrictChoices
14 from tracking.validation.utilities import prob, is_primary
15 from tracking.validation.plot import ValidationPlot
16 
17 import tracking.harvest.peelers as peelers
18 import tracking.harvest.harvesting as harvesting
19 import tracking.harvest.refiners as refiners
20 from tracking.harvest.run import HarvestingRun
21 
22 import trackfindingcdc.harvest.cdc_peelers as cdc_peelers
23 
24 import argparse
25 
26 import logging
27 
28 
29 def get_logger():
30  return logging.getLogger(__name__)
31 
32 CONTACT = "oliver.frost@desy.de"
33 
34 
36  """Harvester to read, postprocess and inspect MC events for track-segment fit validation"""
37 
38 
39  n_events = 10000
40 
41  generator_module = "low_gun" # Rather high momentum tracks should make the tracks rather straight.
42 
43  root_input_file = "low_gun.root"
44 
45 
46  monte_carlo = "no"
47 
48  karimaki_fit = False
49 
50  flight_time_estimation = "none"
51 
52  flight_time_reestimation = False
53 
54  use_alpha_in_drift_length = True
55 
56  flight_time_mass_scale = float("nan")
57 
58 
59  fit_positions = "recoPos"
60 
61  fit_variance = "proper"
62 
63  @property
64  def output_file_name(self):
65  """Get the output ROOT filename"""
66  file_name = "karimaki" if self.karimaki_fit else "riemann"
67  file_name += "-mc-" + self.monte_carlo
69  file_name += "-re"
70  if self.root_input_file:
71  file_name += "-" + os.path.split(self.root_input_file)[1]
72  else:
73  file_name += ".root"
74 
75  return file_name
76 
77  def harvesting_module(self, path=None):
78  """Harvest and post-process the MC events"""
79  harvesting_module = SegmentFitValidationModule(self.output_file_name)
80  if path:
81  path.add_module(harvesting_module)
82  return harvesting_module
83 
84  def create_argument_parser(self, **kwds):
85  """Convert command-line arguments to basf2 argument list"""
86  argument_parser = super().create_argument_parser(**kwds)
87 
88  argument_parser.add_argument(
89  '-m',
90  '--monte-carlo',
91  choices=["no", "medium", "full"],
92  default=self.monte_carlo,
93  dest='monte_carlo',
94  help='Amount of monte carlo information to be used in the segment generation.',
95  )
96 
97  argument_parser.add_argument(
98  "-k",
99  "--karimaki",
100  dest="karimaki_fit",
101  action="store_true",
102  help='Use Karimaki fit instead of Riemann fit'
103  )
104 
105  argument_parser.add_argument(
106  "-fp",
107  "--fit-pos",
108  choices=["recoPos", "rlDriftCircle", "wirePos"],
109  default=self.fit_positions,
110  dest="fit_positions",
111  help=("Choose which positional information the segment fit should be used. \n"
112  "* 'wirePos' means only the wire position\n"
113  "* 'recoPos' means only the reconstructed position\n"
114  "* 'rlDriftCircle' means only the drift circle with the right left passage\n")
115  )
116 
117  argument_parser.add_argument(
118  "-fv",
119  "--fit-var",
120  choices=["unit", "driftLength", "pseudo", "proper"],
121  default=self.fit_variance,
122  dest="fit_variance",
123  help=("Choose which variance information the segment fit should be used. \n"
124  "* 'unit' means equal variance of 1\n"
125  "* 'driftLength' means inserting the drift length as variance, very improper because dimension mismatch\n"
126  "* 'pseudo' means the squared dirft length + plus the drift length variance "
127  "(correct dimension, proper lower bound)\n"
128  "* 'proper' means only the drift length variance\n")
129  )
130 
131  argument_parser.add_argument(
132  "-ft",
133  "--flight-time-estimation",
134  choices=["none", "outwards", "downwards"],
135  default=self.flight_time_estimation,
136  dest="flight_time_estimation",
137  help=("Choose which estimation method for the time of flight should be use. \n"
138  "* 'none' no time of flight corrections\n"
139  "* 'outwards' means the minimal time needed to travel to the wire from the interaction point \n"
140  "* 'downwards' means the minimal time needed to travel to the wire from the y = 0 plane downwards \n")
141  )
142 
143  argument_parser.add_argument(
144  "-fr",
145  "--flight-time-reestimation",
146  action="store_true",
147  dest="flight_time_reestimation",
148  help=("Switch to reestimate drift length before fitting.")
149  )
150 
151  argument_parser.add_argument(
152  "-fa",
153  "--use-alpha-in-drift-length",
154  action="store_true",
155  dest="use_alpha_in_drift_length",
156  help=("Switch to serve the alpha angle to the drift length translator.")
157  )
158 
159  argument_parser.add_argument(
160  "-fm",
161  "--flight-time-mass-scale",
162  type=float,
163  dest="flight_time_mass_scale",
164  help=("Mass parameter to estimate the velocity in the time of flight estimation")
165  )
166 
167  return argument_parser
168 
169  def create_path(self):
170  """
171  Sets up a path that plays back pregenerated events or generates events
172  based on the properties in the base class.
173  """
174  path = super().create_path()
175 
176  path.add_module("TFCDC_WireHitPreparer",
177  flightTimeEstimation=self.flight_time_estimation,
178  UseNLoops=1
179  )
180 
181  path.add_module("TFCDC_ClusterPreparer")
182 
183 
184  if self.monte_carlo == "no":
185  # MC free - default
186  path.add_module(
187  "TFCDC_SegmentFinderFacetAutomaton",
188  # investigate=[],
189  SegmentOrientation="none",
190  )
191 
192  elif self.monte_carlo == "medium":
193  # Medium MC - proper generation logic, but true facets and facet relations
194  path.add_module("TFCDC_SegmentFinderFacetAutomaton",
195  FacetFilter="truth",
196  FacetRelationFilter="truth",
197  SegmentOrientation="none",
198  # SegmentOrientation="outwards"
199  )
200 
201  elif self.monte_carlo == "full":
202  # Only true monte carlo segments, but make the positions realistic
203  path.add_module("TFCDC_SegmentCreatorMCTruth",
204  reconstructedDriftLength=False,
205  reconstructedPositions=True)
206 
207  else:
208  raise ValueError("Invalid degree of Monte Carlo information")
209 
210  path.add_module("TFCDC_SegmentFitter",
211  inputSegments="CDCSegment2DVector",
212  karimakiFit=self.karimaki_fit,
213  fitPos=self.fit_positions,
214  fitVariance=self.fit_variance,
215  updateDriftLength=self.flight_time_reestimation,
216  useAlphaInDriftLength=self.use_alpha_in_drift_length,
217  tofMassScale=self.flight_time_mass_scale)
218 
219  return path
220 
221 
222 class SegmentFitValidationModule(harvesting.HarvestingModule):
223 
224  """Module to collect information about the generated segments and
225  compose validation plots on terminate."""
226 
227  def __init__(self, output_file_name):
228  """Constructor"""
229  super().__init__(foreach='CDCSegment2DVector',
230  output_file_name=output_file_name)
231 
232 
233  self.mc_segment_lookup = None
234 
235  def initialize(self):
236  """Receive signal at the start of event processing"""
237  super().initialize()
238 
239 
241 
243 
244  def prepare(self):
245  """Initialize the MC-hit lookup method"""
247 
248  def pick(self, segment):
249  """Select segments with 4 or more hits and a matching primary MC particle"""
250  mc_segment_lookup = self.mc_segment_lookup
251  mc_particle = mc_segment_lookup.getMCParticle(segment)
252 
253  # Check that mc_particle is not a nullptr
254  return mc_particle and is_primary(mc_particle) and segment.size() > 3
255 
256  def peel(self, segment):
257  """Aggregate the track and MC information for track-segment analysis"""
258 
259  mc_segment_lookup = self.mc_segment_lookup
260  mc_hit_lookup = self.mc_hit_lookup
261 
262  segment_crops = cdc_peelers.peel_segment2d(segment)
263 
264  mc_particle = mc_segment_lookup.getMCParticle(segment)
265  mc_particle_crops = peelers.peel_mc_particle(mc_particle)
266 
267  fit3d_truth = mc_segment_lookup.getTrajectory3D(segment)
268  curvature_truth = fit3d_truth.getCurvatureXY()
269 
270  reconstructed_backward = mc_segment_lookup.isForwardOrBackwardToMCTrack(segment)
271  is_fake = mc_segment_lookup.getMCTrackId(segment) < 0
272 
273  if reconstructed_backward != 1:
274  curvature_truth = -curvature_truth
275 
276  truth_crops = dict(
277  tan_lambda_truth=fit3d_truth.getTanLambda(),
278  curvature_truth=curvature_truth,
279  abs_curvature_truth=abs(curvature_truth),
280  curvature_residual=segment_crops["curvature_estimate"] - curvature_truth,
281  curvature_pull=(segment_crops["curvature_estimate"] - curvature_truth) / segment_crops["curvature_variance"],
282  reconstructed_backward=reconstructed_backward,
283  is_fake=is_fake
284  )
285 
286  rl_sum = 0
287  n_correct = 0
288  n_total = segment.size()
289  first_i_incorrect = float("nan")
290  last_i_incorrect = float("nan")
291  first_l_incorrect = 0
292  last_l_incorrect = 0
293 
294  n_rl_switch = 0
295  last_rl_info = 0
296  for i, reco_hit2d in enumerate(segment):
297  rl_info = reco_hit2d.getRLInfo()
298  hit = reco_hit2d.getWireHit().getHit()
299  true_rl_info = mc_hit_lookup.getRLInfo(hit)
300 
301  if rl_info != last_rl_info:
302  n_rl_switch += 1
303  last_rl_info = rl_info
304 
305  if true_rl_info == rl_info:
306  n_correct += 1
307  else:
308  if first_i_incorrect != first_i_incorrect:
309  first_i_incorrect = i
310  first_l_incorrect = reco_hit2d.getRefDriftLength()
311  last_i_incorrect = i
312  last_l_incorrect = reco_hit2d.getRefDriftLength()
313  if rl_info == 1:
314  rl_sum += 1
315  elif rl_info == -1 or rl_info == 65535: # <- The root interface mistakes the signed enum value for an unsigned value
316  rl_sum -= 1
317 
318  alias_score = segment.getAliasScore()
319 
320  rl_crops = dict(
321  n_total=n_total,
322  n_correct=n_correct,
323  n_incorrect=n_total - n_correct,
324  rl_purity=n_correct / n_total,
325  first_incorrect_location=first_i_incorrect / (n_total - 1),
326  first_l_incorrect=first_l_incorrect,
327  last_incorrect_location=last_i_incorrect / (n_total - 1),
328  last_l_incorrect=last_l_incorrect,
329  n_rl_switch=n_rl_switch,
330  n_rl_asymmetry=rl_sum / n_total,
331  n_abs_rl_asymmetry=abs(rl_sum) / n_total,
332  may_alias=alias_score == 0,
333  alias_score=alias_score
334  )
335 
336  segment_crops.update(truth_crops)
337  segment_crops.update(mc_particle_crops)
338  segment_crops.update(rl_crops)
339  return segment_crops
340 
341  # Refiners to be executed at the end of the harvesting / termination of the module
342 
343  save_histograms = refiners.save_histograms(outlier_z_score=5.0, allow_discrete=True)
344 
345  save_tree = refiners.save_tree()
346 
347 
348  save_curvature_pull = refiners.save_pull_analysis(
349  part_name="curvature",
350  unit="1/cm",
351  absolute=False,
352  aux_names=["tan_lambda_truth", "abs_curvature_truth"],
353  groupby=[
354  "stereo_kind",
355  "superlayer_id"
356  ],
357  outlier_z_score=4.0,
358  title_postfix="")
359 
360 
361  save_curvature_pull_rl_pure = refiners.save_pull_analysis(
362  part_name="curvature",
363  unit="1/cm",
364  absolute=True,
365  filter_on="rl_purity",
366  filter=lambda rl_purity: rl_purity > 0.5,
367  aux_names=["tan_lambda_truth", "abs_curvature_truth"],
368  groupby=[
369  "stereo_kind",
370  "superlayer_id"
371  ],
372  outlier_z_score=4.0,
373  title_postfix="",
374  folder_name="rl_pure/{groupby_addition}"
375  )
376 
377 
378  save_absolute_curvature_pull = refiners.save_pull_analysis(
379  part_name="curvature",
380  unit="1/cm",
381  absolute=True,
382  aux_names=["tan_lambda_truth", "abs_curvature_truth"],
383  groupby=[
384  "stereo_kind",
385  "superlayer_id"
386  ],
387  outlier_z_score=4.0,
388  title_postfix="")
389 
390 
391  save_fit_quality_histograms = refiners.save_histograms(
392  outlier_z_score=5.0,
393  select={
394  "ndf": "ndf",
395  "chi2": "#chi2",
396  "p_value": "p-value",
397  "tan_lambda_truth": "true tan #lambda",
398  },
399  groupby=[None, "stereo_kind", "superlayer_id"],
400  title="Histogram of {part_name}{stacked_by_indication}{stackby}",
401  description="Distribution of {part_name} in the segment fits",
402  )
403 
404 
405 def main():
407  run.configure_and_execute_from_commandline()
408 
409 if __name__ == "__main__":
410  logging.basicConfig(stream=sys.stdout, level=logging.INFO, format='%(levelname)s:%(message)s')
411  main()
cdc_peelers.peel_segment2d
def peel_segment2d(segment, key="{part_name}")
Definition: cdc_peelers.py:8
record.SegmentFitValidationRun.flight_time_reestimation
bool flight_time_reestimation
do not re-estimate the flight time
Definition: record.py:52
record.SegmentFitValidationModule.__init__
def __init__(self, output_file_name)
Definition: record.py:227
tracking.harvest.refiners
Definition: refiners.py:1
record.SegmentFitValidationRun.harvesting_module
def harvesting_module(self, path=None)
Definition: record.py:77
record.SegmentFitValidationRun.create_path
def create_path(self)
Definition: record.py:169
tracking.run.utilities
Definition: utilities.py:1
record.SegmentFitValidationRun.fit_positions
string fit_positions
fit the reconstructed hit position for the segment fit
Definition: record.py:59
tracking.harvest.harvesting
Definition: harvesting.py:1
Belle2::TrackFindingCDC::CDCMCSegment2DLookUp::getInstance
static const CDCMCSegment2DLookUp & getInstance()
Getter for the singletone instance.
Definition: CDCMCSegment2DLookUp.cc:23
tracking.harvest.peelers
Definition: peelers.py:1
record.SegmentFitValidationRun
Definition: record.py:35
record.SegmentFitValidationRun.karimaki_fit
bool karimaki_fit
use the Riemann fit instead of the Karimaki fit
Definition: record.py:48
record.SegmentFitValidationRun.use_alpha_in_drift_length
bool use_alpha_in_drift_length
use the alpha parameter for drift length
Definition: record.py:54
record.SegmentFitValidationModule
Definition: record.py:222
tracking.harvest.run
Definition: run.py:1
main
int main(int argc, char **argv)
Run all tests.
Definition: test_main.cc:77
tracking.run.event_generation.ReadOrGenerateEventsRun.root_input_file
root_input_file
generating events, so there is no ROOT input file
Definition: event_generation.py:120
record.SegmentFitValidationRun.flight_time_estimation
string flight_time_estimation
do not estimate the flight time
Definition: record.py:50
tracking.harvest.run.HarvestingRunMixin.output_file_name
output_file_name
Disable the writing of an output ROOT file.
Definition: run.py:12
Belle2::TrackFindingCDC::CDCMCHitLookUp::getInstance
static const CDCMCHitLookUp & getInstance()
Getter for the singletone instance.
Definition: CDCMCHitLookUp.cc:32
tracking.validation.utilities
Definition: utilities.py:1
record.SegmentFitValidationRun.flight_time_mass_scale
flight_time_mass_scale
assume a very heavy particle for flight time calculation
Definition: record.py:56
record.SegmentFitValidationModule.prepare
def prepare(self)
Definition: record.py:244
tracking.harvest.run.HarvestingRun
Definition: run.py:69
record.SegmentFitValidationRun.monte_carlo
string monte_carlo
do not generate new events
Definition: record.py:46
record.SegmentFitValidationModule.peel
def peel(self, segment)
Definition: record.py:256
tracking.validation.plot
Definition: plot.py:1
record.SegmentFitValidationModule.initialize
def initialize(self)
Definition: record.py:235
record.SegmentFitValidationRun.fit_variance
string fit_variance
use the drift-length variance in the segment fit
Definition: record.py:61
record.SegmentFitValidationModule.mc_hit_lookup
mc_hit_lookup
Method to find matching MC hits.
Definition: record.py:242
record.SegmentFitValidationModule.pick
def pick(self, segment)
Definition: record.py:248
record.SegmentFitValidationModule.mc_segment_lookup
mc_segment_lookup
by default, there is no method to find matching MC track segments
Definition: record.py:233
record.SegmentFitValidationRun.create_argument_parser
def create_argument_parser(self, **kwds)
Definition: record.py:84