Belle II Software prerelease-10-00-00a
module.py
1#!/usr/bin/env python3
2
3
10
11import math
12import collections
13import numpy as np
14
15from tracking.root_utils import root_save_name
16
17from tracking.validation.plot import ValidationPlot, compose_axis_label
18
19from tracking.validation.pull import PullAnalysis
20from tracking.validation.resolution import ResolutionAnalysis
21from tracking.validation.fom import ValidationFiguresOfMerit
22from tracking.validation.utilities import (
23 getHelixFromMCParticle,
24 getSeedTrackFitResult,
25 is_primary,
26 get_det_hit_ids,
27 calc_ndf_from_det_hit_ids
28)
29
30import basf2
31
32import ROOT
33from ROOT import Belle2
34import os
35
36ROOT.gSystem.Load("libtracking")
37
38
39class FilterProperties:
40 """
41 contains all information necessary for track filters to decide whether
42 track will be included into the processed list of tracks
43 This class is used for both providing information on pattern reco and
44 MC tracks
45 """
46
47 def __init__(
48 self,
49 trackCand=None,
50 mcParticle=None,
51 mcParticles=None,
52 wasFitted=False,
53 fitResult=None,
54 seedResult=None,
55 ):
56 """Constructor"""
57
58 ## cached value of the track candidate
59 self.trackCand = trackCand
60 ## cached value of the MC particle
61 self.mcParticle = mcParticle
62 ## cached value of the MCParticles StoreArray
63 self.mcParticles = mcParticles
64 ## cached value of the fitted flag
65 self.wasFitted = wasFitted
66 ## cached value of the fit result
67 self.fitResult = fitResult
68 ## cached value of the seed result
69 self.seedResult = seedResult
70
71
72# This class will accept all pattern reco and mctracks
73# It is used as the default option for the TrackingValidationModule
74#
75# It can be used as a starting point to implement your own particle filter class
76# and supply it to the TrackingValidationModule via the init method's
77# track_filter_object property.
78#
79# doesPrPass is called for all pattern reconstructed tracks
80# Objects available infilterProperties:
81# trackCand is guaranteed to be != None
82# If wasFitted == True, fitResult will be set
83# otherwise seedResult will be set
84#
85# doesMcPass is called for all MC tracks
86# Objects available infilterProperties:
87# mcParticle is guaranteed to be != None
88#
89
90class AlwaysPassFilter:
91 """Filter that always passes"""
92
93 def doesPrPass(self, filterProperties):
94 """Pattern-reconstructed track always passes"""
95 return True
96
97 def doesMcPass(self, filterProperties):
98 """MC track always passes"""
99 return True
100
101
102class TrackingValidationModule(basf2.Module):
103
104 """Module to collect matching information about the found particles and to
105 generate validation plots and figures of merit on the performance of track finding."""
106
107 def __init__(
108 self,
109 name,
110 contact,
111 fit=False,
112 pulls=False,
113 resolution=False,
114 output_file_name=None,
115 track_filter_object=AlwaysPassFilter(),
116 plot_name_postfix='',
117 plot_title_postfix='',
118 exclude_profile_mc_parameter='',
119 exclude_profile_pr_parameter='',
120 use_expert_folder=True,
121 trackCandidatesColumnName="RecoTracks",
122 mcTrackCandidatesColumName="MCRecoTracks",
123 non_expert_parameters=['p_{t}']
124 ):
125 """Constructor"""
126
127 super().__init__()
128
129 ## cached value of the tracking-validation name
130 self.validation_name = name
131 ## cached value of the contact person name
132 self.contact = contact
133 ## cached value of the track fit
134 self.fit = fit
135 ## cached values of the track-fit pulls
136 self.pulls = pulls
137 ## cached value of the resolution
138 self.resolution = resolution
139 ## cached value of the output ROOT TFile
140 self.output_file_name = output_file_name or self.validation_name \
141 + 'TrackingValidation.root'
142 ## cached value of the track-filter object
143 self.track_filter_object = track_filter_object
144 ## cached value of the suffix appended to the plot names
145 self.plot_name_postfix = plot_name_postfix
146 ## cached value of the suffix appended to the plot titles
147 self.plot_title_postfix = plot_title_postfix
148 ## cached list of perigee parameters excluded from PR side plots
149 self.exclude_profile_pr_parameter = exclude_profile_pr_parameter
150 ## cached list of perigee parameters excluded from MC side plots
151 self.exclude_profile_mc_parameter = exclude_profile_mc_parameter
152 ## cached flag to use the "expert" folder for the pull and residual plots
153 self.use_expert_folder = use_expert_folder
154 ## cached name of the RecoTracks StoreArray
155 self.trackCandidatesColumnName = trackCandidatesColumnName
156 ## cached name of the MCRecoTracks StoreArray
157 self.mcTrackCandidatesColumnName = mcTrackCandidatesColumName
158 ## list of parameters that determines which plots (all with corresponding x-axis) are marked as shifter plots
159 self.non_expert_parameters = non_expert_parameters
160
161 ## default binning used for resolution plots over pt
162 self.resolution_pt_binning = [0.05, 0.1, 0.25, 0.4, 0.6, 1., 1.5, 2., 3., 4.]
163
164 ## If this variable is set the code will open the file with same name as the file created here
165 ## and will read the binning from the TH1/TProfile with same name as the one created here. If you
166 ## do not want this feature either remove the corresponding root files from the validation
167 ## directory (this will trigger the default behaviour) or set the environmental variable DO_NOT_READ_BINNING
168 self.referenceFileName = None
169 if "DO_NOT_READ_BINNING" not in os.environ:
170 # the validity of the file will be checked later
171 self.referenceFileName = Belle2.FileSystem.findFile("tracking/validation/" + self.output_file_name, True)
172 basf2.B2INFO("Will read binning from: " + self.referenceFileName)
173 basf2.B2INFO("If this is not wanted set the environment variable DO_NOT_READ_BINNING or remove reference files.")
174 else:
175 basf2.B2INFO("Will not read binning from reference files.")
176
177 def initialize(self):
178 """Receive signal at the start of event processing"""
179
180 ## Track-match object that examines relation information from MCMatcherTracksModule
181 self.trackMatchLookUp = Belle2.TrackMatchLookUp(self.mcTrackCandidatesColumnName, self.trackCandidatesColumnName)
182
183 ## Use dequeues in favour of lists to prevent repeated memory allocation of cost O(n)
184
185 ## list of PR-track clones and matches
186 self.pr_clones_and_matches = collections.deque()
187 ## list of PR-track matches
188 self.pr_matches = collections.deque()
189 ## list of PR-track fakes
190 self.pr_fakes = collections.deque()
191
192 ## list of PR-track seed pt values
193 self.pr_seed_pt = collections.deque()
194 ## list of PR-track seed tan(lambda) values
195 self.pr_seed_tan_lambdas = collections.deque()
196 ## list of PR-track seed phi values
197 self.pr_seed_phi = collections.deque()
198
199 ## list of PR-track seed omega-truth values
200 self.pr_omega_truths = collections.deque()
201 ## list of PR-track seed omega-estimate values
202 self.pr_omega_estimates = collections.deque()
203 ## list of PR-track seed omega-variance values
204 self.pr_omega_variances = collections.deque()
205
206 ## list of PR-track seed tan(lambda)-truth values
207 self.pr_tan_lambda_truths = collections.deque()
208 ## list of PR-track seed tan(lambda)-estimate values
209 self.pr_tan_lambda_estimates = collections.deque()
210 ## list of PR-track seed tan(lambda)-variance values
211 self.pr_tan_lambda_variances = collections.deque()
212
213 ## list of PR-track seed d0-truth values
214 self.pr_d0_truths = collections.deque()
215 ## list of PR-track seed d0-estimate values
216 self.pr_d0_estimates = collections.deque()
217 ## list of PR-track seed d0-variance values
218 self.pr_d0_variances = collections.deque()
219
220 ## list of PR-track seed z0-truth values
221 self.pr_z0_truths = collections.deque()
222 ## list of PR-track seed z0-estimate values
223 self.pr_z0_estimates = collections.deque()
224
225 ## list of PR-track seed pt-truth values
226 self.pr_pt_truths = collections.deque()
227 ## list of PR-track seed pt-estimate values
228 self.pr_pt_estimates = collections.deque()
229
230 ## list of PR-track binning values
231 self.pr_bining_pt = collections.deque()
232
233 ## list of MC-track matches
234 self.mc_matches = collections.deque()
235 ## list of MC-track matches, including matched charge
236 self.mc_charge_matches = collections.deque()
237 ## list of MC-track matches charge asymmetry
238 self.mc_charge_asymmetry = collections.deque()
239 ## list of MC-track matches charge asymmetry weights
240 self.mc_charge_asymmetry_weights = collections.deque()
241 ## list of MC-track primaries
242 self.mc_primaries = collections.deque()
243 ## list of MC-track d0 values
244 self.mc_d0s = collections.deque()
245 ## list of MC-track tan(lambda) values
246 self.mc_tan_lambdas = collections.deque()
247 ## direction of the track in phi
248 self.mc_phi = collections.deque()
249 ## list of MC-track pt values
250 self.mc_pts = collections.deque()
251 ## list of MC-track hit efficiencies
252 self.mc_hit_efficiencies = collections.deque()
253 ## list of MC-track multiplicities
254 self.mc_multiplicities = collections.deque()
255 ## list of MC-track number of degrees of freedom
256 self.mc_ndf = collections.deque()
257
258 def event(self):
259 """Event method"""
260
261 self.examine_pr_tracks()
262 self.examine_mc_tracks()
263
264 def examine_pr_tracks(self):
265 """Looks at the individual pattern reconstructed tracks and store information about them"""
266
267 # Analyse from the pattern recognition side
268 trackMatchLookUp = self.trackMatchLookUp
269
270 trackCands = Belle2.PyStoreArray(self.trackCandidatesColumnName)
271 mcParticles = Belle2.PyStoreArray("MCParticles")
272 if not trackCands:
273 return
274
275 for trackCand in trackCands:
276 is_matched = trackMatchLookUp.isAnyChargeMatchedPRRecoTrack(trackCand)
277 is_clone = trackMatchLookUp.isAnyChargeClonePRRecoTrack(trackCand)
278
279 pt_truth = float('nan')
280 omega_truth = float('nan')
281 tan_lambda_truth = float('nan')
282 d0_truth = float('nan')
283 z0_truth = float('nan')
284
285 mcParticle = None
286 if is_matched or is_clone:
287 # Only matched and clone tracks have a related MCParticle
288 mcParticle = trackMatchLookUp.getRelatedMCParticle(trackCand)
289 mcHelix = getHelixFromMCParticle(mcParticle)
290 omega_truth = mcHelix.getOmega()
291 tan_lambda_truth = mcHelix.getTanLambda()
292 pt_truth = mcParticle.getMomentum().Rho()
293 d0_truth = mcHelix.getD0()
294 z0_truth = mcHelix.getZ0()
295
296 # fill the FilterProperties will all properties on this track
297 # gathered so far
298 filterProperties = FilterProperties(trackCand=trackCand,
299 mcParticle=mcParticle, mcParticles=mcParticles)
300
301 if self.fit:
302 prTrackFitResult = \
303 trackMatchLookUp.getRelatedTrackFitResult(trackCand)
304 filterProperties.wasFitted = prTrackFitResult is not None
305 filterProperties.fitResult = prTrackFitResult
306 else:
307 prTrackFitResult = getSeedTrackFitResult(trackCand)
308 filterProperties.seedResult = prTrackFitResult
309
310 # skip this track due to the filtering rules ?
311 if not self.track_filter_object.doesPrPass(filterProperties):
312 continue
313
314 omega_estimate = float('nan')
315 omega_variance = float('nan')
316 tan_lambda_estimate = float('nan')
317 tan_lambda_variance = float('nan')
318 d0_estimate = float('nan')
319 d0_variance = float('nan')
320 z0_estimate = float('nan')
321 pt_estimate = float('nan')
322
323 momentum = float('nan')
324
325 # store seed information, they are always available from the pattern reco
326 # even if the fit was no successful
327 # this information can we used when plotting fake tracks, for example
328 seed_momentum = trackCand.getMomentumSeed()
329 # Avoid zero division exception
330 seed_tan_lambda = np.divide(1.0, math.tan(seed_momentum.Theta()))
331 seed_phi = seed_momentum.Phi()
332 seed_pt = seed_momentum.Rho()
333
334 if prTrackFitResult:
335 omega_estimate = prTrackFitResult.getOmega()
336 omega_variance = prTrackFitResult.getCov()[9]
337
338 tan_lambda_estimate = prTrackFitResult.getCotTheta()
339 tan_lambda_variance = prTrackFitResult.getCov()[14]
340
341 d0_estimate = prTrackFitResult.getD0()
342 d0_variance = prTrackFitResult.getCov()[0]
343
344 z0_estimate = prTrackFitResult.getZ0()
345
346 momentum = prTrackFitResult.getMomentum()
347 pt_estimate = momentum.Rho()
348
349 # store properties of the seed
350 self.pr_seed_pt.append(seed_pt)
351 self.pr_seed_tan_lambdas.append(seed_tan_lambda)
352 self.pr_seed_phi.append(seed_phi)
353
354 self.pr_bining_pt.append(pt_truth)
355
356 # store properties resulting from this trackfit
357 isMatchedOrIsClone = is_matched or is_clone
358 self.pr_clones_and_matches.append(isMatchedOrIsClone)
359 self.pr_matches.append(is_matched)
360
361 self.pr_fakes.append(not isMatchedOrIsClone)
362
363 self.pr_omega_estimates.append(omega_estimate)
364 self.pr_omega_variances.append(omega_variance)
365 self.pr_omega_truths.append(omega_truth)
366
367 self.pr_tan_lambda_estimates.append(tan_lambda_estimate)
368 self.pr_tan_lambda_variances.append(tan_lambda_variance)
369 self.pr_tan_lambda_truths.append(tan_lambda_truth)
370
371 self.pr_d0_estimates.append(d0_estimate)
372 self.pr_d0_variances.append(d0_variance)
373 self.pr_d0_truths.append(d0_truth)
374
375 self.pr_z0_estimates.append(z0_estimate)
376 self.pr_z0_truths.append(z0_truth)
377
378 self.pr_pt_estimates.append(pt_estimate)
379 self.pr_pt_truths.append(pt_truth)
380
381 def examine_mc_tracks(self):
382 """Looks at the individual Monte Carlo tracks and store information about them"""
383
384 trackMatchLookUp = self.trackMatchLookUp
385
386 # Analyse from the Monte Carlo reference side
387 mcTrackCands = Belle2.PyStoreArray(self.mcTrackCandidatesColumnName)
388 mcParticles = Belle2.PyStoreArray('MCParticles')
389 if not mcTrackCands:
390 return
391
392 multiplicity = mcTrackCands.getEntries()
393 multiplicity_primaries = multiplicity
394
395 # measure the charge asymmetry
396 n_matched_plus = 0
397 n_matched_minus = 0
398
399 for mcTrackCand in mcTrackCands:
400 is_matched = trackMatchLookUp.isAnyChargeMatchedMCRecoTrack(mcTrackCand)
401
402 relatedPRtrackCand = trackMatchLookUp.getRelatedPRRecoTrack(mcTrackCand)
403 if relatedPRtrackCand:
404 is_chargeMatched = trackMatchLookUp.isChargeMatched(relatedPRtrackCand)
405 else:
406 is_chargeMatched = False
407
408 if is_chargeMatched:
409 if mcTrackCand.getChargeSeed() > 0:
410 n_matched_plus += 1
411 else:
412 n_matched_minus += 1
413
414 hit_efficiency = trackMatchLookUp.getRelatedEfficiency(mcTrackCand)
415 if math.isnan(hit_efficiency):
416 hit_efficiency = 0
417
418 mcParticle = trackMatchLookUp.getRelatedMCParticle(mcTrackCand)
419 mcHelix = getHelixFromMCParticle(mcParticle)
420
421 # fill the FilterProperties will all properties on this track
422 # gathered so far
423 filterProperties = FilterProperties(mcParticle=mcParticle,
424 mcParticles=mcParticles)
425
426 if not self.track_filter_object.doesMcPass(filterProperties):
427 continue
428
429 momentum = mcParticle.getMomentum()
430 pt = momentum.Rho()
431 tan_lambda = np.divide(1.0, math.tan(momentum.Theta())) # Avoid zero division exception
432 d0 = mcHelix.getD0()
433 det_hit_ids = get_det_hit_ids(mcTrackCand)
434 ndf = calc_ndf_from_det_hit_ids(det_hit_ids)
435
436 self.mc_matches.append(is_matched)
437 self.mc_charge_matches.append(is_chargeMatched and is_matched)
438 self.mc_primaries.append(is_primary(mcParticle))
439 self.mc_hit_efficiencies.append(hit_efficiency)
440 self.mc_pts.append(pt)
441 self.mc_d0s.append(d0)
442 self.mc_tan_lambdas.append(tan_lambda)
443 self.mc_multiplicities.append(multiplicity)
444 self.mc_phi.append(momentum.Phi())
445 self.mc_ndf.append(ndf)
446 if not is_primary(mcParticle):
447 multiplicity_primaries -= 1
448
449 charge_asymmetry = (n_matched_plus - n_matched_minus)/(n_matched_plus +
450 n_matched_minus) if (n_matched_plus + n_matched_minus) != 0 else 0
451 for mcTrackCand in mcTrackCands:
452 if is_primary(mcParticle):
453 self.mc_charge_asymmetry.append(charge_asymmetry)
454 self.mc_charge_asymmetry_weights.append(1./multiplicity_primaries)
455 else:
456 self.mc_charge_asymmetry.append(0)
457 self.mc_charge_asymmetry_weights.append(0)
458
459 def terminate(self):
460 """Receive signal at the end of event processing"""
461 name = self.validation_name
462 contact = self.contact
463
464 # Overall figures of merit #
465 ############################
466
467 mc_matched_primaries = np.logical_and(self.mc_primaries, self.mc_matches)
468
469 charge_asymmetry = np.average(self.mc_charge_asymmetry, weights=self.mc_charge_asymmetry_weights)
470 if len(mc_matched_primaries) > 0 and sum(mc_matched_primaries) > 0:
471 charge_efficiency = np.average(self.mc_charge_matches, weights=mc_matched_primaries)
472 hit_efficiency = np.average(self.mc_hit_efficiencies, weights=mc_matched_primaries)
473 else:
474 charge_efficiency = float('nan')
475 hit_efficiency = float('nan')
476 finding_charge_efficiency = np.average(self.mc_charge_matches, weights=self.mc_primaries)
477 finding_efficiency = np.average(self.mc_matches, weights=self.mc_primaries)
478 fake_rate = 1.0 - np.mean(self.pr_clones_and_matches)
479 # can only be computed if there are entries
480 if len(self.pr_clones_and_matches) > 0 and sum(self.pr_clones_and_matches) > 0:
481 clone_rate = 1.0 - np.average(self.pr_matches,
482 weights=self.pr_clones_and_matches)
483 else:
484 clone_rate = float('nan')
485
486 figures_of_merit = ValidationFiguresOfMerit(f'{name}_figures_of_merit')
487 figures_of_merit['finding_charge_efficiency'] = finding_charge_efficiency
488 figures_of_merit['finding_efficiency'] = finding_efficiency
489 figures_of_merit['charge_efficiency'] = charge_efficiency
490 figures_of_merit['charge_asymmetry'] = charge_asymmetry
491 figures_of_merit['fake_rate'] = fake_rate
492 figures_of_merit['clone_rate'] = clone_rate
493 figures_of_merit['hit_efficiency'] = hit_efficiency
494
495 figures_of_merit.description = \
496 """
497finding_efficiency - the ratio of matched Monte Carlo tracks to all primary Monte Carlo tracks <br/>
498charge_efficiency - the ratio of matched Monte Carlo tracks with correct charge to matched primary Monte Carlo tracks <br/>
499finding_charge_efficiency - the ratio of matched Monte Carlo tracks with correct charge to all primary Monte Carlo tracks <br/>
500fake_rate - ratio of pattern recognition tracks that are not related to a particle
501 (background, ghost) to all pattern recognition tracks <br/>
502clone_rate - ratio of clones divided the number of tracks that are related to a particle (clones and matches) <br/>
503
504"""
505 figures_of_merit.check = 'Compare for degradations with respect to the reference'
506 figures_of_merit.contact = contact
507 print(figures_of_merit)
508
509 # Validation plots #
510 ####################
511 validation_plots = []
512 pull_analyses = []
513
514 # Finding efficiency #
515 ######################
516 plots = self.profiles_by_mc_parameters(self.mc_matches,
517 'finding efficiency',
518 make_hist=False,
519 weights=self.mc_primaries)
520
521 validation_plots.extend(plots)
522
523 # Fake rate (all tracks not matched or clone #
524 # use TrackCand seeds for the fake track plotting #
525 # as the fit (if successful) is probably not meaningful #
526 #########################################################
527 print('fake list: ' + str(self.pr_fakes.count(1)))
528 plots = self.profiles_by_pr_parameters(self.pr_fakes, 'fake rate',
529 make_hist=False)
530
531 validation_plots.extend(plots)
532
533 # Charge efficiency of matched primary tracks #
534 #######################################
535 plots = self.profiles_by_mc_parameters(self.mc_charge_matches,
536 'charge efficiency for matched primary tracks',
537 weights=mc_matched_primaries)
538
539 validation_plots.extend(plots)
540
541 # Finding & Charge efficiency of primary tracks #
542 #######################################
543 plots = self.profiles_by_mc_parameters(self.mc_charge_matches,
544 'finding and charge efficiency for primary tracks',
545 weights=self.mc_primaries)
546
547 validation_plots.extend(plots)
548
549 # Charge asymmetry of primary tracks #
550 #######################################
551 plots = self.profiles_by_mc_parameters(self.mc_charge_asymmetry,
552 'charge asymmetry for primary tracks',
553 weights=self.mc_charge_asymmetry_weights,
554 is_asymmetry=True)
555
556 validation_plots.extend(plots)
557
558 # Hit efficiency #
559 ##################
560 plots = self.profiles_by_mc_parameters(self.mc_hit_efficiencies,
561 'hit efficiency with matched tracks',
562 weights=mc_matched_primaries)
563
564 validation_plots.extend(plots)
565
566 # Fit quality #
567 ###############
568 if self.pulls:
569 all_but_diagonal_plots = list(PullAnalysis.default_which_plots)
570 all_but_diagonal_plots.remove("diag_profile")
571 all_but_diagonal_plots.remove("diag_scatter")
572
573 plot_name_prefix = name + self.plot_name_postfix
574 if not self.fit:
575 plot_name_prefix += '_seed'
576
577 # Omega / curvature pull
578 pr_omega_truths = np.array(self.pr_omega_truths)
579 pr_omega_estimates = np.array(self.pr_omega_estimates)
580 pr_omega_variances = np.array(self.pr_omega_variances)
581
582 curvature_pull_analysis = PullAnalysis('#omega', unit='1/cm',
583 plot_name_prefix=plot_name_prefix + '_omega',
584 plot_title_postfix=self.plot_title_postfix,
585 referenceFileName=self.referenceFileName)
586
587 curvature_pull_analysis.analyse(pr_omega_truths,
588 pr_omega_estimates,
589 pr_omega_variances,
590 which_plots=all_but_diagonal_plots)
591
592 curvature_pull_analysis.contact = contact
593 pull_analyses.append(curvature_pull_analysis)
594
595 # Tan lambda pull
596 pr_tan_lambda_truths = np.array(self.pr_tan_lambda_truths)
597 pr_tan_lambda_estimates = np.array(self.pr_tan_lambda_estimates)
598 pr_tan_lambda_variances = np.array(self.pr_tan_lambda_variances)
599
600 curvature_pull_analysis = PullAnalysis('tan #lambda',
601 plot_name_prefix=plot_name_prefix + '_tan_lambda',
602 plot_title_postfix=self.plot_title_postfix,
603 referenceFileName=self.referenceFileName)
604
605 curvature_pull_analysis.analyse(pr_tan_lambda_truths,
606 pr_tan_lambda_estimates,
607 pr_tan_lambda_variances,
608 which_plots=all_but_diagonal_plots)
609
610 curvature_pull_analysis.contact = contact
611 pull_analyses.append(curvature_pull_analysis)
612
613 # d0 pull
614 curvature_pull_analysis = PullAnalysis('d0',
615 plot_name_prefix=plot_name_prefix + '_d0',
616 plot_title_postfix=self.plot_title_postfix,
617 referenceFileName=self.referenceFileName)
618
619 curvature_pull_analysis.analyse(np.array(self.pr_d0_truths),
620 np.array(self.pr_d0_estimates),
621 np.array(self.pr_d0_variances),
622 which_plots=all_but_diagonal_plots)
623
624 curvature_pull_analysis.contact = contact
625 pull_analyses.append(curvature_pull_analysis)
626
627 # Resolution plots
628 ##############################
629 if self.resolution:
630 # d0 impact parameter resolution plot
631 d0_resolution_analysis = ResolutionAnalysis('d0_res',
632 self.resolution_pt_binning,
633 'Pt',
634 plot_name_prefix=plot_name_prefix + '_d0_res',
635 plot_title_postfix=self.plot_title_postfix,
636 referenceFileName=self.referenceFileName)
637 d0_resolution_analysis.analyse(np.array(self.pr_bining_pt),
638 np.array(self.pr_d0_truths),
639 np.array(self.pr_d0_estimates))
640 d0_resolution_analysis.contact = contact
641 pull_analyses.append(d0_resolution_analysis)
642
643 # z0 impact parameter resolution plot
644 z0_resolution_analysis = ResolutionAnalysis('z0_res',
645 self.resolution_pt_binning,
646 "Pt",
647 plot_name_prefix=plot_name_prefix + '_z0_res',
648 plot_title_postfix=self.plot_title_postfix,
649 referenceFileName=self.referenceFileName)
650 z0_resolution_analysis.analyse(np.array(self.pr_bining_pt),
651 np.array(self.pr_z0_truths),
652 np.array(self.pr_z0_estimates))
653 z0_resolution_analysis.contact = contact
654 pull_analyses.append(z0_resolution_analysis)
655
656 # omega curvature parameter resolution plot
657 omega_resolution_analysis = ResolutionAnalysis('omega_res',
658 self.resolution_pt_binning,
659 "Pt",
660 plot_name_prefix=plot_name_prefix + '_omega_res',
661 plot_title_postfix=self.plot_title_postfix,
662 referenceFileName=self.referenceFileName)
663 omega_resolution_analysis.analyse(np.array(self.pr_bining_pt),
664 np.array(self.pr_omega_truths),
665 np.array(self.pr_omega_estimates))
666 omega_resolution_analysis.contact = contact
667 pull_analyses.append(omega_resolution_analysis)
668
669 # transverse momentum resolution plot
670 pt_resolution_analysis = ResolutionAnalysis('pt_res',
671 self.resolution_pt_binning,
672 "Pt",
673 plot_name_prefix=plot_name_prefix + '_pt_res',
674 plot_title_postfix=self.plot_title_postfix,
675 referenceFileName=self.referenceFileName)
676 pt_resolution_analysis.analyse(np.array(self.pr_bining_pt),
677 np.array(self.pr_pt_truths),
678 np.array(self.pr_pt_estimates))
679 pt_resolution_analysis.contact = contact
680 pull_analyses.append(pt_resolution_analysis)
681
682 # Saving #
683 ##########
684
685 # Save everything to a ROOT file
686 output_tfile = ROOT.TFile(self.output_file_name, 'recreate')
687
688 # Show all parameters and the fit result in the plots
689 # if viewed in the browser or the validation
690 opt_fit = 0o112
691 ROOT.gStyle.SetOptFit(opt_fit)
692
693 figures_of_merit.write()
694
695 for validation_plot in validation_plots:
696 validation_plot.write()
697
698 if self.use_expert_folder:
699 expert_tdirectory = output_tfile.mkdir('expert', 'Expert')
700 expert_tdirectory.cd()
701 ROOT.gStyle.SetOptFit(opt_fit)
702
703 for pull_analysis in pull_analyses:
704 pull_analysis.write()
705
706 output_tfile.Close()
707
708 def profiles_by_mc_parameters(
709 self,
710 xs,
711 quantity_name,
712 unit=None,
713 parameter_names=[
714 'd_0',
715 'p_t',
716 'tan_lambda',
717 'multiplicity',
718 'phi',
719 'ndf',
720 ],
721 make_hist=True,
722 weights=None,
723 is_asymmetry=False,
724 ):
725 """Create profile histograms by MC-track parameters"""
726
727 # apply exclusion list
728 new_parameter_names = [item for item in parameter_names if item
729 not in self.exclude_profile_mc_parameter]
730
731 # Profile versus the various parameters
732 profile_parameters = {
733 'd_{0}': self.mc_d0s,
734 'p_{t}': self.mc_pts,
735 'tan #lambda': self.mc_tan_lambdas,
736 '#phi': self.mc_phi,
737 'multiplicity': self.mc_multiplicities,
738 'ndf': self.mc_ndf,
739 }
740
741 return self.profiles_by_parameters_base(
742 xs,
743 quantity_name,
744 new_parameter_names,
745 profile_parameters,
746 unit,
747 make_hist,
748 weights=weights,
749 is_asymmetry=is_asymmetry
750 )
751
752 def profiles_by_pr_parameters(
753 self,
754 xs,
755 quantity_name,
756 unit=None,
757 parameter_names=['Seed_p_t', 'Seed tan #lambda', 'Seed #phi'],
758 make_hist=True,
759 ):
760 """Create profile histograms by PR-track parameters"""
761
762 # apply exclusion list
763 new_parameter_names = [item for item in parameter_names if item
764 not in self.exclude_profile_pr_parameter]
765
766 # Profile versus the various parameters
767 profile_parameters = {'Seed p_{t}': self.pr_seed_pt,
768 'Seed tan #lambda': self.pr_seed_tan_lambdas,
769 'Seed #phi': self.pr_seed_phi}
770
771 return self.profiles_by_parameters_base(
772 xs,
773 quantity_name,
774 new_parameter_names,
775 profile_parameters,
776 unit,
777 make_hist,
778 )
779
780 def profiles_by_parameters_base(
781 self,
782 xs,
783 quantity_name,
784 parameter_names,
785 profile_parameters,
786 unit,
787 make_hist,
788 weights=None,
789 is_asymmetry=False,
790 ):
791 """Create profile histograms for generic parameters"""
792
793 contact = self.contact
794
795 validation_plots = []
796 plot_name_prefix = self.validation_name + '_' + root_save_name(quantity_name) \
797 + self.plot_name_postfix
798
799 if make_hist:
800 # Histogram of the quantity
801 histogram = ValidationPlot(plot_name_prefix, self.referenceFileName)
802 histogram.hist(xs, weights=weights)
803
804 histogram.xlabel = quantity_name
805 histogram.description = 'Not a serious plot yet.'
806 histogram.check = ''
807 histogram.contact = contact
808 histogram.title = quantity_name + self.plot_title_postfix
809
810 validation_plots.append(histogram)
811
812 for (parameter_name, parameter_values) in list(profile_parameters.items()):
813 if parameter_name in parameter_names \
814 or root_save_name(parameter_name) in parameter_names:
815
816 is_expert = not (parameter_name in self.non_expert_parameters)
817
818 parameter_root_name = root_save_name(parameter_name)
819
820 # Apply some boundaries for the maximal tracking acceptance
821 # such that the plots look more instructive
822 if 'tan_lambda' in parameter_root_name:
823 lower_bound = -2.0
824 upper_bound = 5.0
825 # need different bounds for cosmics
826 if 'cosmics' in self.validation_name.lower() or \
827 'cosmics' in self.output_file_name.lower():
828 lower_bound = -1.5
829 upper_bound = 1.5
830 elif 'ndf' in parameter_root_name:
831 lower_bound = 0
832 upper_bound = min(200, np.max(parameter_values))
833 elif 'p_t' in parameter_root_name:
834 lower_bound = 0
835 upper_bound = 2.5
836 # need different upper_bound for cosmics
837 if 'cosmics' in self.validation_name.lower() or \
838 'cosmics' in self.output_file_name.lower():
839 upper_bound = 30
840 elif 'd_0' in parameter_root_name:
841 lower_bound = -0.06
842 upper_bound = 0.06
843 # need different bounds for cosmics
844 if 'cosmics' in self.validation_name.lower() or \
845 'cosmics' in self.output_file_name.lower():
846 lower_bound = -20
847 upper_bound = 20
848 else:
849 lower_bound = None
850 upper_bound = None
851
852 profile_plot_name = plot_name_prefix + '_by_' \
853 + root_save_name(parameter_name)
854 profile_plot = ValidationPlot(profile_plot_name, self.referenceFileName)
855 profile_plot.profile(parameter_values,
856 xs,
857 weights=weights,
858 outlier_z_score=10.0,
859 lower_bound=lower_bound,
860 upper_bound=upper_bound,
861 y_binary=True,
862 is_expert=is_expert,
863 is_asymmetry=is_asymmetry)
864
865 profile_plot.xlabel = compose_axis_label(parameter_name)
866 profile_plot.ylabel = compose_axis_label(quantity_name, unit)
867 profile_plot.title = quantity_name + ' by ' + parameter_name \
868 + ' profile' + self.plot_title_postfix
869
870 profile_plot.description = \
871 'Dependence of %s of the track over the true %s' \
872 % (quantity_name, parameter_name)
873 profile_plot.check = 'Variations should be low'
874 profile_plot.contact = contact
875 validation_plots.append(profile_plot)
876
877 return validation_plots
tracking.validation.module.FilterProperties.__init__
__init__(self, trackCand=None, mcParticle=None, mcParticles=None, wasFitted=False, fitResult=None, seedResult=None)
Definition module.py:55
tracking.root_utils
Definition root_utils.py:1