7 from ROOT
import Belle2
22 import trackfindingcdc.harvest.cdc_peelers
as cdc_peelers
30 return logging.getLogger(__name__)
32 CONTACT =
"oliver.frost@desy.de"
36 """Harvester to read, postprocess and inspect MC events for track-segment fit validation"""
41 generator_module =
"low_gun"
43 root_input_file =
"low_gun.root"
50 flight_time_estimation =
"none"
52 flight_time_reestimation =
False
54 use_alpha_in_drift_length =
True
56 flight_time_mass_scale = float(
"nan")
59 fit_positions =
"recoPos"
61 fit_variance =
"proper"
65 """Get the output ROOT filename"""
66 file_name =
"karimaki" if self.
karimaki_fit else "riemann"
78 """Harvest and post-process the MC events"""
81 path.add_module(harvesting_module)
82 return harvesting_module
85 """Convert command-line arguments to basf2 argument list"""
88 argument_parser.add_argument(
91 choices=[
"no",
"medium",
"full"],
94 help=
'Amount of monte carlo information to be used in the segment generation.',
97 argument_parser.add_argument(
102 help=
'Use Karimaki fit instead of Riemann fit'
105 argument_parser.add_argument(
108 choices=[
"recoPos",
"rlDriftCircle",
"wirePos"],
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")
117 argument_parser.add_argument(
120 choices=[
"unit",
"driftLength",
"pseudo",
"proper"],
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")
131 argument_parser.add_argument(
133 "--flight-time-estimation",
134 choices=[
"none",
"outwards",
"downwards"],
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")
143 argument_parser.add_argument(
145 "--flight-time-reestimation",
147 dest=
"flight_time_reestimation",
148 help=(
"Switch to reestimate drift length before fitting.")
151 argument_parser.add_argument(
153 "--use-alpha-in-drift-length",
155 dest=
"use_alpha_in_drift_length",
156 help=(
"Switch to serve the alpha angle to the drift length translator.")
159 argument_parser.add_argument(
161 "--flight-time-mass-scale",
163 dest=
"flight_time_mass_scale",
164 help=(
"Mass parameter to estimate the velocity in the time of flight estimation")
167 return argument_parser
171 Sets up a path that plays back pregenerated events or generates events
172 based on the properties in the base class.
176 path.add_module(
"TFCDC_WireHitPreparer",
181 path.add_module(
"TFCDC_ClusterPreparer")
187 "TFCDC_SegmentFinderFacetAutomaton",
189 SegmentOrientation=
"none",
194 path.add_module(
"TFCDC_SegmentFinderFacetAutomaton",
196 FacetRelationFilter=
"truth",
197 SegmentOrientation=
"none",
203 path.add_module(
"TFCDC_SegmentCreatorMCTruth",
204 reconstructedDriftLength=
False,
205 reconstructedPositions=
True)
208 raise ValueError(
"Invalid degree of Monte Carlo information")
210 path.add_module(
"TFCDC_SegmentFitter",
211 inputSegments=
"CDCSegment2DVector",
224 """Module to collect information about the generated segments and
225 compose validation plots on terminate."""
229 super().
__init__(foreach=
'CDCSegment2DVector',
230 output_file_name=output_file_name)
236 """Receive signal at the start of event processing"""
245 """Initialize the MC-hit lookup method"""
249 """Select segments with 4 or more hits and a matching primary MC particle"""
251 mc_particle = mc_segment_lookup.getMCParticle(segment)
254 return mc_particle
and is_primary(mc_particle)
and segment.size() > 3
257 """Aggregate the track and MC information for track-segment analysis"""
264 mc_particle = mc_segment_lookup.getMCParticle(segment)
265 mc_particle_crops = peelers.peel_mc_particle(mc_particle)
267 fit3d_truth = mc_segment_lookup.getTrajectory3D(segment)
268 curvature_truth = fit3d_truth.getCurvatureXY()
270 reconstructed_backward = mc_segment_lookup.isForwardOrBackwardToMCTrack(segment)
271 is_fake = mc_segment_lookup.getMCTrackId(segment) < 0
273 if reconstructed_backward != 1:
274 curvature_truth = -curvature_truth
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,
288 n_total = segment.size()
289 first_i_incorrect = float(
"nan")
290 last_i_incorrect = float(
"nan")
291 first_l_incorrect = 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)
301 if rl_info != last_rl_info:
303 last_rl_info = rl_info
305 if true_rl_info == rl_info:
308 if first_i_incorrect != first_i_incorrect:
309 first_i_incorrect = i
310 first_l_incorrect = reco_hit2d.getRefDriftLength()
312 last_l_incorrect = reco_hit2d.getRefDriftLength()
315 elif rl_info == -1
or rl_info == 65535:
318 alias_score = segment.getAliasScore()
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
336 segment_crops.update(truth_crops)
337 segment_crops.update(mc_particle_crops)
338 segment_crops.update(rl_crops)
343 save_histograms = refiners.save_histograms(outlier_z_score=5.0, allow_discrete=
True)
345 save_tree = refiners.save_tree()
348 save_curvature_pull = refiners.save_pull_analysis(
349 part_name=
"curvature",
352 aux_names=[
"tan_lambda_truth",
"abs_curvature_truth"],
361 save_curvature_pull_rl_pure = refiners.save_pull_analysis(
362 part_name=
"curvature",
365 filter_on=
"rl_purity",
366 filter=
lambda rl_purity: rl_purity > 0.5,
367 aux_names=[
"tan_lambda_truth",
"abs_curvature_truth"],
374 folder_name=
"rl_pure/{groupby_addition}"
378 save_absolute_curvature_pull = refiners.save_pull_analysis(
379 part_name=
"curvature",
382 aux_names=[
"tan_lambda_truth",
"abs_curvature_truth"],
391 save_fit_quality_histograms = refiners.save_histograms(
396 "p_value":
"p-value",
397 "tan_lambda_truth":
"true tan #lambda",
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",
407 run.configure_and_execute_from_commandline()
409 if __name__ ==
"__main__":
410 logging.basicConfig(stream=sys.stdout, level=logging.INFO, format=
'%(levelname)s:%(message)s')