11 from .plot
import ValidationPlot, compose_axis_label
13 from .pull
import PullAnalysis
14 from .resolution
import ResolutionAnalysis
15 from .fom
import ValidationFiguresOfMerit
16 from .utilities
import (
17 getHelixFromMCParticle,
18 getSeedTrackFitResult,
21 calc_ndf_from_det_hit_ids
29 ROOT.gSystem.Load(
"libtracking")
30 from ROOT
import Belle2
37 contains all informations necessary for track filters to decide whether
38 track will be included into the processed list of tracks
39 This class is used for both providing information on pattern reco and
87 """Filter that always passes"""
90 """Pattern-reconstructed track always passes"""
94 """MC track always passes"""
100 """Module to collect matching information about the found particles and to
101 generate validation plots and figures of merit on the performance of track finding."""
110 output_file_name=None,
111 track_filter_object=AlwaysPassFilter(),
112 plot_name_postfix=
'',
113 plot_title_postfix=
'',
114 exclude_profile_mc_parameter=
'',
115 exclude_profile_pr_parameter=
'',
116 use_expert_folder=
True,
117 trackCandidatesColumnName=
"RecoTracks",
118 mcTrackCandidatesColumName=
"MCRecoTracks"
136 +
'TrackingValidation.root'
162 if "DO_NOT_READ_BINNING" not in os.environ:
166 basf2.B2INFO(
"If this is not wanted set the environment variable DO_NOT_READ_BINNING or remove reference files.")
168 basf2.B2INFO(
"Will not read binning from reference files.")
171 """Receive signal at the start of event processing"""
254 """Looks at the individual pattern reconstructed tracks and store information about them"""
264 for trackCand
in trackCands:
265 is_matched = trackMatchLookUp.isMatchedPRRecoTrack(trackCand)
266 is_clone = trackMatchLookUp.isClonePRRecoTrack(trackCand)
268 pt_truth = float(
'nan')
269 omega_truth = float(
'nan')
270 tan_lambda_truth = float(
'nan')
271 d0_truth = float(
'nan')
272 z0_truth = float(
'nan')
275 if is_matched
or is_clone:
277 mcParticle = trackMatchLookUp.getRelatedMCParticle(trackCand)
278 mcHelix = getHelixFromMCParticle(mcParticle)
279 omega_truth = mcHelix.getOmega()
280 tan_lambda_truth = mcHelix.getTanLambda()
281 pt_truth = mcParticle.getMomentum().Perp()
282 d0_truth = mcHelix.getD0()
283 z0_truth = mcHelix.getZ0()
288 mcParticle=mcParticle, mcParticles=mcParticles)
292 trackMatchLookUp.getRelatedTrackFitResult(trackCand)
293 filterProperties.wasFitted = prTrackFitResult
is not None
294 filterProperties.fitResult = prTrackFitResult
296 prTrackFitResult = getSeedTrackFitResult(trackCand)
297 filterProperties.seedResult = prTrackFitResult
303 omega_estimate = float(
'nan')
304 omega_variance = float(
'nan')
305 tan_lambda_estimate = float(
'nan')
306 tan_lambda_variance = float(
'nan')
307 d0_estimate = float(
'nan')
308 d0_variance = float(
'nan')
309 z0_estimate = float(
'nan')
310 pt_estimate = float(
'nan')
312 momentum_pt = float(
'nan')
313 momentum = float(
'nan')
318 seed_position = trackCand.getPositionSeed()
319 seed_momentum = trackCand.getMomentumSeed()
321 seed_tan_lambda = np.divide(1.0, math.tan(seed_momentum.Theta()))
322 seed_phi = seed_position.Phi()
323 seed_theta = seed_position.Theta()
326 omega_estimate = prTrackFitResult.getOmega()
327 omega_variance = prTrackFitResult.getCov()[9]
329 tan_lambda_estimate = prTrackFitResult.getCotTheta()
330 tan_lambda_variance = prTrackFitResult.getCov()[14]
332 d0_estimate = prTrackFitResult.getD0()
333 d0_variance = prTrackFitResult.getCov()[0]
335 z0_estimate = prTrackFitResult.getZ0()
337 momentum = prTrackFitResult.getMomentum()
338 pt_estimate = momentum.Perp()
348 isMatchedOrIsClone = is_matched
or is_clone
351 self.
pr_fakes.append(
not isMatchedOrIsClone)
372 """Looks at the individual Monte Carlo tracks and store information about them"""
382 multiplicity = mcTrackCands.getEntries()
384 for mcTrackCand
in mcTrackCands:
385 is_matched = trackMatchLookUp.isMatchedMCRecoTrack(mcTrackCand)
387 hit_efficiency = trackMatchLookUp.getRelatedEfficiency(mcTrackCand)
388 if math.isnan(hit_efficiency):
391 mcParticle = trackMatchLookUp.getRelatedMCParticle(mcTrackCand)
392 mcHelix = getHelixFromMCParticle(mcParticle)
397 mcParticles=mcParticles)
402 momentum = mcParticle.getMomentum()
404 tan_lambda = np.divide(1.0, math.tan(momentum.Theta()))
406 det_hit_ids = get_det_hit_ids(mcTrackCand)
407 ndf = calc_ndf_from_det_hit_ids(det_hit_ids)
416 self.
mc_theta.append(momentum.Theta())
417 self.
mc_phi.append(momentum.Phi())
421 """Receive signal at the end of event processing"""
431 clone_rate = 1.0 - np.average(self.
pr_matches,
434 clone_rate = float(
'nan')
441 figures_of_merit[
'finding_efficiency'] = finding_efficiency
442 figures_of_merit[
'fake_rate'] = fake_rate
443 figures_of_merit[
'clone_rate'] = clone_rate
444 figures_of_merit[
'hit_efficiency'] = hit_efficiency
446 figures_of_merit.description = \
448 finding_efficiency - the ratio of matched Monte Carlo tracks to all Monte Carlo tracks <br/>
449 fake_rate - ratio of pattern recognition tracks that are not related to a particle
450 (background, ghost) to all pattern recognition tracks <br/>
451 clone_rate - ratio of clones divided the number of tracks that are related to a particle (clones and matches) <br/>
453 figures_of_merit.check =
'Compare for degradations with respect to the reference'
454 figures_of_merit.contact = contact
455 print(figures_of_merit)
459 validation_plots = []
465 'finding efficiency',
469 validation_plots.extend(plots)
475 print(
'fake list: ' + str(len(self.
pr_fakes)))
478 validation_plots.extend(plots)
483 'hit efficiency with matched tracks',
484 weights=mc_matched_primaries)
485 validation_plots.extend(plots)
490 all_but_diagonal_plots = list(PullAnalysis.default_which_plots)
491 all_but_diagonal_plots.remove(
"diag_profile")
492 all_but_diagonal_plots.remove(
"diag_scatter")
496 plot_name_prefix +=
'_seed'
503 curvature_pull_analysis =
PullAnalysis(
'#omega', unit=
'1/cm',
504 plot_name_prefix=plot_name_prefix +
'_omega',
508 curvature_pull_analysis.analyse(pr_omega_truths,
511 which_plots=all_but_diagonal_plots)
513 curvature_pull_analysis.contact = contact
514 pull_analyses.append(curvature_pull_analysis)
522 plot_name_prefix=plot_name_prefix +
'_tan_lambda',
526 curvature_pull_analysis.analyse(pr_tan_lambda_truths,
527 pr_tan_lambda_estimates,
528 pr_tan_lambda_variances,
529 which_plots=all_but_diagonal_plots)
531 curvature_pull_analysis.contact = contact
532 pull_analyses.append(curvature_pull_analysis)
536 plot_name_prefix=plot_name_prefix +
'_d0',
540 curvature_pull_analysis.analyse(np.array(self.
pr_d0_truths),
543 which_plots=all_but_diagonal_plots)
545 curvature_pull_analysis.contact = contact
546 pull_analyses.append(curvature_pull_analysis)
555 plot_name_prefix=plot_name_prefix +
'_d0_res',
558 d0_resolution_analysis.analyse(np.array(self.
pr_bining_pt),
561 d0_resolution_analysis.contact = contact
562 pull_analyses.append(d0_resolution_analysis)
568 plot_name_prefix=plot_name_prefix +
'_z0_res',
571 z0_resolution_analysis.analyse(np.array(self.
pr_bining_pt),
574 z0_resolution_analysis.contact = contact
575 pull_analyses.append(z0_resolution_analysis)
581 plot_name_prefix=plot_name_prefix +
'_omega_res',
584 omega_resolution_analysis.analyse(np.array(self.
pr_bining_pt),
587 omega_resolution_analysis.contact = contact
588 pull_analyses.append(omega_resolution_analysis)
594 plot_name_prefix=plot_name_prefix +
'_pt_res',
597 pt_resolution_analysis.analyse(np.array(self.
pr_bining_pt),
600 pt_resolution_analysis.contact = contact
601 pull_analyses.append(pt_resolution_analysis)
612 ROOT.gStyle.SetOptFit(opt_fit)
614 figures_of_merit.write()
616 for validation_plot
in validation_plots:
617 validation_plot.write()
620 expert_tdirectory = output_tfile.mkdir(
'expert',
'Expert')
621 expert_tdirectory.cd()
622 ROOT.gStyle.SetOptFit(opt_fit)
624 for pull_analysis
in pull_analyses:
625 pull_analysis.write()
643 non_expert_parameters=['p_{t}'],
647 """Create profile histograms by MC-track parameters"""
650 new_parameter_names = [item
for item
in parameter_names
if item
654 profile_parameters = {
671 non_expert_parameters=non_expert_parameters,
680 parameter_names=['Seed tan
683 """Create profile histograms by PR-track parameters"""
686 new_parameter_names = [item
for item
in parameter_names
if item
711 non_expert_parameters=[],
714 """Create profile histograms for generic parameters"""
718 validation_plots = []
719 plot_name_prefix = self.
validation_name +
'_' + root_save_name(quantity_name) \
725 histogram.hist(xs, weights=weights)
727 histogram.xlabel = quantity_name
728 histogram.description =
'Not a serious plot yet.'
730 histogram.contact = contact
733 validation_plots.append(histogram)
735 for (parameter_name, parameter_values)
in list(profile_parameters.items()):
736 if parameter_name
in parameter_names \
737 or root_save_name(parameter_name)
in parameter_names:
739 is_expert = not(parameter_name
in non_expert_parameters)
743 if root_save_name(parameter_name) ==
'tan_lambda':
746 elif root_save_name(parameter_name) ==
'theta':
747 lower_bound = 17 * math.pi / 180
748 upper_bound = 150 * math.pi / 180
749 elif root_save_name(parameter_name) ==
'ndf':
751 upper_bound = min(200, np.max(parameter_values))
756 profile_plot_name = plot_name_prefix +
'_by_' \
757 + root_save_name(parameter_name)
759 profile_plot.profile(parameter_values,
762 outlier_z_score=10.0,
763 lower_bound=lower_bound,
764 upper_bound=upper_bound,
768 profile_plot.xlabel = compose_axis_label(parameter_name)
769 profile_plot.ylabel = compose_axis_label(quantity_name, unit)
770 profile_plot.title = quantity_name +
' by ' + parameter_name \
773 profile_plot.description = \
774 'Dependence of %s of the track over the true %s' \
775 % (quantity_name, parameter_name)
776 profile_plot.check =
'Variations should be low'
777 profile_plot.contact = contact
778 validation_plots.append(profile_plot)
780 return validation_plots
785 standard_reco_run = StandardReconstructionEventsRun()
786 standard_reco_run.configure_from_commandline()
788 validation_module = SeparatedTrackingValidationModule(name=
"test_run",
790 output_file_name=
"test_separated_module.root",
793 standard_reco_run.add_module(validation_module)
794 standard_reco_run.execute()
797 if __name__ ==
"__main__":
798 logging.basicConfig(level=logging.INFO)