Belle II Software  release-08-01-10
validationcomparison.py
1 #!/usr/bin/env python3
2 
3 
10 
11 """ Compare ROOT objects and perform e.g. chi2 tests.
12 A small command line interface for testing/debugging purposes is included.
13 Run `python3 validationcomparison.py --help` for more information. """
14 
15 # std
16 from abc import ABC, abstractmethod
17 import argparse
18 import numpy
19 import os.path
20 from typing import Optional, Any
21 
22 # 3rd
23 import ROOT
24 
25 # ours
26 from metaoptions import MetaOptionParser
27 
28 # Unfortunately doxygen has some trouble with inheritance of attributes, so
29 # we disable it.
30 # @cond SUPPRESS_DOXYGEN
31 
32 
33 # ==============================================================================
34 # Custom Exceptions
35 # ==============================================================================
36 
37 
38 class ComparisonFailed(Exception):
39  """
40  The comparison failed for some reason. For example
41  because ROOT was not able to compute the Chi^2 properly
42  """
43 
44 
45 class ObjectsNotSupported(Exception):
46  """
47  The type and/or combination of provided ROOT objects
48  is not supported for comparison
49  """
50 
51 
52 class DifferingBinCount(Exception):
53  """
54  The two ROOT objects provided have a different bin count
55  and therefor, cannot be compared using the Chi2 test
56  """
57 
58 
59 class TooFewBins(Exception):
60  """
61  Not sufficient bins to perform the Chi^2 test
62  """
63 
64 
65 # ==============================================================================
66 # Comparison class selector
67 # ==============================================================================
68 
69 
70 def get_comparison(
71  object_1, object_2, mop: Optional[MetaOptionParser] = None,
72 ) -> "ComparisonBase":
73  """ Uses the metaoptions to determine which comparison algorithm is used
74  and initializes the corresponding subclass of :class:`ComparisonBase` that
75  implements the actual comparison and holds the results.
76  @param object_1 ROOT TObject
77  @param object_2 ROOT TObject
78  @param mop Metaoption parser
79  """
80  if mop is None:
81  mop = MetaOptionParser()
82  if mop.has_option("kolmogorov"):
83  tester: Any = KolmogorovTest
84  elif mop.has_option("andersondarling"):
85  tester = AndersonDarlingTest
86  else:
87  tester = Chi2Test
88 
89  test = tester(object_1, object_2, mop=mop)
90 
91  return test
92 
93 
94 # ==============================================================================
95 # Comparison Base Class
96 # ==============================================================================
97 
98 
99 class ComparisonBase(ABC):
100  """
101  Base class for all comparison implementations.
102 
103  Follows 3 steps:
104 
105  1. Initialize the class together with two ROOT objects of different
106  revisions (that are to be compared) and the metaoptions (given in the
107  corresponding validation (steering) file), that determine how to compare
108  them.
109 
110  2. The Comparison class saves the ROOT objects and the metaoptions
111  internally, but does not compute anything yet
112 
113  3. If :meth:`ensure_compute` is called, or any property is accessed that
114  depends on computation, the internal implementation :meth:`_compute`
115  (to be implemented in the subclass) is called.
116 
117  4. :meth:`_compute` ensures that all values, like chi2, p-value etc. are
118  computed
119 
120  5. Two properties :meth:`comparison_result` (pass/warning/error) and
121  :meth:`comparison_result_long` (longer description of the comparison result)
122  allow to access the results.
123  """
124 
125  def __init__(
126  self,
127  object_a,
128  object_b,
129  mop: Optional[MetaOptionParser] = None,
130  debug=False,
131  ):
132  """
133  Initialize ComparisonBase class
134 
135  :param object_a:
136  :param object_b:
137  :param mop: MetaOptionParser
138  :param debug (bool): Debug mode enabled?
139  """
140 
141  self.object_a = object_a
142 
143 
144  self.object_b = object_b
145 
146 
147  if mop is None:
148  mop = MetaOptionParser()
149  self.mop = mop
150 
151 
152  self.debug = debug
153 
154 
155  self.computed = False
156 
157 
158  self._comparison_result = "not_compared"
159 
161  self._comparison_result_long = ""
162 
163  def ensure_compute(self):
164  """
165  Ensure all required quantities get computed and are cached inside the
166  class
167  """
168  if self.computed:
169  return
170 
171  if self.mop.has_option("nocompare"):
172  # is comparison disabled for this plot ?
173  self._comparison_result_long = "Testing is disabled for this plot"
174  return
175 
176  fail_message = "Comparison failed: "
177 
178  # Note: default for comparison_result is "not_compared"
179  try:
180  self._compute()
181  except ObjectsNotSupported as e:
182  self._comparison_result_long = fail_message + str(e)
183  except DifferingBinCount as e:
184  self._comparison_result = "error"
185  self._comparison_result_long = fail_message + str(e)
186  except TooFewBins as e:
187  self._comparison_result_long = fail_message + str(e)
188  except ComparisonFailed as e:
189  self._comparison_result = "error"
190  self._comparison_result_long = fail_message + str(e)
191  except Exception as e:
192  self._comparison_result = "error"
193  self._comparison_result_long = (
194  "Unknown error occurred. Please "
195  "submit a bug report. " + str(e)
196  )
197  else:
198  # Will be already set in case of errors above and we don't want
199  # to overwrite this.
200  self._comparison_result_long = self._get_comparison_result_long()
201  self._comparison_result = self._get_comparison_result()
202 
203  self.computed = True
204 
205  @abstractmethod
206  def _get_comparison_result(self) -> str:
207  """ Used to format the value of :attr:`_comparison_result`. """
208 
209  @abstractmethod
210  def _get_comparison_result_long(self) -> str:
211  """ Used to format the value of :attr:`_comparison_result_long`. """
212 
213  @property
214  def comparison_result(self):
215  """ Comparison result, i.e. pass/warning/error """
216  self.ensure_compute()
217  return self._comparison_result
218 
219  @property
220  def comparison_result_long(self):
221  """ Longer description of the comparison result """
222  self.ensure_compute()
223  return self._comparison_result_long
224 
225  @abstractmethod
226  def _compute(self):
227  """ This method performs the actual computations. """
228 
229  def can_compare(self):
230  """
231  @return: True if the two objects can be compared, False otherwise
232  """
233  return self._has_correct_types() and self._has_compatible_bins()
234 
235  def _has_correct_types(self) -> bool:
236  """
237  @return: True if the two objects have a) a type supported for
238  comparison and b) can be compared with each other
239  """
240  if self.object_a is None or self.object_b is None:
241  return False
242 
243  # check if the supplied object inherit from one of the supported types
244  # and if they are of the same type
245  supported_types = ["TProfile", "TH1D", "TH1F", "TEfficiency"]
246  if self.object_a.ClassName() != self.object_b.ClassName():
247  return False
248  if self.object_a.ClassName() not in supported_types:
249  return False
250 
251  if self.object_a.ClassName() == "TEfficiency":
252  # can only handle TEfficiencies with dimension one atm
253  if self.object_a.GetDimension() > 1:
254  return False
255 
256  return True
257 
258  def _raise_has_correct_types(self) -> None:
259  """
260  Raise Exception if not the two objects have a) a type supported for
261  comparison and b) can be compared with each other
262  @return: None
263  """
264  if not self._has_correct_types():
265  msg = (
266  "Comparison of {} (Type {}) with {} (Type {}) not "
267  "supported.\nPlease open a GitLab issue (validation "
268  "label) if you need this supported. "
269  )
270  raise ObjectsNotSupported(
271  msg.format(
272  self.object_a.GetName(),
273  self.object_a.ClassName(),
274  self.object_b.GetName(),
275  self.object_b.ClassName(),
276  )
277  )
278 
279  def _has_compatible_bins(self) -> bool:
280  """
281  Check if both ROOT obeject have the same amount of bins
282  @return: True if the bins are equal, otherwise False
283  """
284  if (
285  self.object_a.ClassName()
286  == "TEfficiency"
287  == self.object_b.ClassName()
288  ):
289  nbins_a = self.object_a.GetTotalHistogram().GetNbinsX()
290  nbins_b = self.object_b.GetTotalHistogram().GetNbinsX()
291  else:
292  nbins_a = self.object_a.GetNbinsX()
293  nbins_b = self.object_b.GetNbinsX()
294 
295  return nbins_a == nbins_b
296 
297  def _raise_has_compatible_bins(self) -> None:
298  """
299  Raise Exception if not both ROOT obeject have the same amount of bins
300  @return: None
301  """
302  if not self._has_compatible_bins():
303  msg = (
304  "The objects have differing x bin count: {} has {} vs. {} "
305  "has {}."
306  )
307  raise DifferingBinCount(
308  msg.format(
309  self.object_a.GetName(),
310  self.object_a.GetNbinsX(),
311  self.object_b.GetName(),
312  self.object_b.GetNbinsX(),
313  )
314  )
315 
316  @staticmethod
317  def _convert_teff_to_hist(teff_a):
318  """
319  Convert the content of a TEfficiency plot to a histogram and set
320  the bin content and errors
321  """
322  conv_hist = teff_a.GetTotalHistogram()
323  xbin_count = conv_hist.GetNbinsX()
324  xbin_low = conv_hist.GetXaxis().GetXmin()
325  xbin_max = conv_hist.GetXaxis().GetXmax()
326 
327  th1 = ROOT.TH1D(
328  teff_a.GetName() + "root_conversion",
329  teff_a.GetName(),
330  xbin_count,
331  xbin_low,
332  xbin_max,
333  )
334  # starting from the first to the last bin, ignoring the under/overflow
335  # bins
336  for i in range(1, xbin_count):
337  th1.SetBinContent(i, teff_a.GetEfficiency(i))
338  th1.SetBinError(i, teff_a.GetEfficiencyErrorLow(i))
339 
340  return th1
341 
342 
343 class PvalueTest(ComparisonBase):
344  """ Test with a pvalue """
345 
346 
348  _default_pvalue_warn = 0.99
349 
350 
352  _default_pvalue_error = 0.01
353 
354  def __init__(self, *args, **kwargs):
355  """ Initialize Pvalue test
356 
357  Args:
358  *args: Positional arguments to ComparisonBase
359  **kwargs: Keyword arguments to ComparisonBase
360  """
361  super().__init__(*args, **kwargs)
362 
363  self._pvalue = None
364 
365  self._pvalue_warn = self.mop.pvalue_warn()
366 
367  self._pvalue_error = self.mop.pvalue_error()
368 
369  if self._pvalue_warn is None:
370  self._pvalue_warn = self._default_pvalue_warn
371  if self._pvalue_error is None:
372  self._pvalue_error = self._default_pvalue_error
373 
374  def _get_comparison_result(self) -> str:
375  if self._pvalue is None:
376  return "error"
377 
378  if self._pvalue < self._pvalue_error:
379  return "error"
380  elif self._pvalue < self._pvalue_warn:
381  return "warning"
382  else:
383  return "equal"
384 
385  @abstractmethod
386  def _compute(self):
387  pass
388 
389  @abstractmethod
390  def _get_comparison_result_long(self):
391  pass
392 
393 
394 # ==============================================================================
395 # Implementation of specific comparison algorithms
396 # ==============================================================================
397 
398 # ------------------------------------------------------------------------------
399 # Chi2 Test
400 # ------------------------------------------------------------------------------
401 
402 
403 class Chi2Test(PvalueTest):
404 
405  """
406  Perform a Chi2Test for ROOT objects. The chi2 test method is e.g. described
407  in the documentation of TH1::Chi2Test. Basically this class wraps around
408  this Chi2Test function, and takes care that we can call perform these
409  tests for a wider selection of ROOT objects.
410  """
411 
412  def __init__(self, *args, **kwargs):
413  """
414  Initialize Chi2Test.
415  :param args: See arguments of :class:`ComparisonBase`
416  :param kwargs: See arguments of :class:`ComparisonBase`
417  """
418  super().__init__(*args, **kwargs)
419 
420  # The following attributes will be set in :meth:`_compute`
421 
422 
423  self._chi2 = None
424 
425  self._chi2ndf = None
426 
427  self._ndf = None
428 
429  def _ensure_zero_error_has_no_content(self, a, b):
430  """
431  Ensure there are no bins which have a content set, but 0 error
432  This bin content will be set to 0 to disable this bin completely during
433  the comparison
434  """
435  nbins = a.GetNbinsX()
436  for ibin in range(1, nbins + 1):
437  if a.GetBinError(ibin) <= 0.0 and b.GetBinError(ibin) <= 0.0:
438  # set the bin content of the profile plots to zero so ROOT
439  # will ignore this bin in its comparison
440  a.SetBinContent(ibin, 0.0)
441  b.SetBinContent(ibin, 0.0)
442  if self.debug:
443  print(
444  "DEBUG: Warning: Setting bin content of bin {} to "
445  "zero for both histograms, because both histograms "
446  "have vanishing errors there.".format(ibin)
447  )
448 
449  def _compute(self) -> None:
450  """
451  Performs the actual Chi^2 test
452  @return: None
453  """
454  self._raise_has_correct_types()
455  self._raise_has_compatible_bins()
456 
457  local_object_a = self.object_a
458  local_object_b = self.object_b
459 
460  # very special handling for TEfficiencies
461  if self.object_a.ClassName() == "TEfficiency":
462  local_object_a = self._convert_teff_to_hist(self.object_a)
463  local_object_b = self._convert_teff_to_hist(self.object_b)
464  if self.debug:
465  print("Converting TEfficiency objects to histograms.")
466 
467  nbins = local_object_a.GetNbinsX()
468 
469  if nbins < 2:
470  raise TooFewBins(
471  "{} bin(s) is too few to perform the Chi2 "
472  "test.".format(nbins)
473  )
474 
475  weighted_types = ["TProfile", "TH1D", "TH1F"]
476  comp_weight_a = local_object_a.ClassName() in weighted_types
477  comp_weight_b = local_object_b.ClassName() in weighted_types
478 
479  # clone, because possibly some content of profiles will
480  # be set to zero
481  first_obj = local_object_a.Clone()
482  second_obj = local_object_b.Clone()
483 
484  if comp_weight_a and not comp_weight_b:
485  # switch histograms, because ROOT can only have the first one
486  # to be unweighted
487  first_obj, second_obj = second_obj, first_obj
488  if self.debug:
489  print(
490  "Debug: Warning: Switching the two objects, because "
491  "ROOT can only have the first one to be unweighted"
492  )
493 
494  # Construct the option string for the Chi2Test call
495  comp_options = "P " # for debugging output
496  if comp_weight_a and comp_weight_b:
497  comp_options += "WW"
498  elif comp_weight_a or comp_weight_b:
499  comp_options += "UW"
500  else:
501  comp_options += "UU"
502 
503  if comp_weight_a and comp_weight_b:
504  self._ensure_zero_error_has_no_content(first_obj, second_obj)
505 
506  # use numpy arrays to support ROOT's pass-by-reference interface here
507  res_chi2 = numpy.array([1], numpy.float64)
508  res_igood = numpy.array([1], numpy.int32)
509  res_ndf = numpy.array([1], numpy.int32)
510 
511  res_pvalue = first_obj.Chi2TestX(
512  second_obj, res_chi2, res_ndf, res_igood, comp_options
513  )
514 
515  if self.debug:
516  print("Performing our own chi2 test, with bin-by-bin results: ")
517  print()
518  print_contents_and_errors(first_obj, second_obj)
519  print()
520  print(
521  "Here's what ROOT's Chi2Test gave us (comp_options: '{}'):"
522  " ".format(comp_options)
523  )
524 
525  tp = TablePrinter(3, width=(10, 10, 40))
526  print()
527  tp.print_divider()
528  tp.print(["Key", "Value", "Comment"])
529  tp.print_divider()
530  tp.print(
531  [
532  "chi2",
533  numpy.asscalar(res_chi2),
534  "Should roughly match above 'Total chi2'",
535  ]
536  )
537  tp.print(["ndf", numpy.asscalar(res_ndf), "#Non-empty bins - 1"])
538  tp.print(["chi2/ndf", numpy.asscalar(res_chi2 / res_ndf), ""])
539  tp.print(
540  [
541  "igood",
542  numpy.asscalar(res_igood),
543  "a debug indicator, 0 if all good",
544  ]
545  )
546  tp.print(["pvalue", res_pvalue, ""])
547  tp.print_divider()
548  print()
549  print(
550  "See https://root.cern.ch/doc/master/classTH1.html for more "
551  "information."
552  )
553  print()
554 
555  if res_ndf < 1:
556  msg = (
557  "Comparison failed, no Chi^2 could be computed. For "
558  "debugging, you can use the CLI of "
559  "'validation/scripts/validationcomparison.py' on your root "
560  "file and the reference. Run 'validationcomparison.py "
561  "--help' for info. If problem persists, please open "
562  "GitLab issue (validation label)."
563  )
564  raise ComparisonFailed(msg)
565 
566  res_chi2ndf = res_chi2 / res_ndf
567 
568  self._pvalue, self._chi2, self._chi2ndf, self._ndf = (
569  res_pvalue,
570  res_chi2[0],
571  res_chi2ndf[0],
572  res_ndf[0],
573  )
574 
575  def _get_comparison_result_long(self) -> str:
576  if self._pvalue is None or self._chi2ndf is None or self._chi2 is None:
577  return (
578  r"Could not perform $\chi^2$-Test between {{revision1}} "
579  r"and {{revision2}} due to an unknown error. Please "
580  r"submit a bug report."
581  )
582 
583  return (
584  r"Performed $\chi^2$-Test between {{revision1}} "
585  r"and {{revision2}} "
586  r"($\chi^2$ = {chi2:.4f}; NDF = {ndf}; "
587  r"$\chi^2/\text{{{{NDF}}}}$ = {chi2ndf:.4f})."
588  r" <b>p-value: {pvalue:.6f}</b> (p-value warn: {pvalue_warn}, "
589  r"p-value error: {pvalue_error})".format(
590  chi2=self._chi2,
591  ndf=self._ndf,
592  chi2ndf=self._chi2ndf,
593  pvalue=self._pvalue,
594  pvalue_warn=self._pvalue_warn,
595  pvalue_error=self._pvalue_error,
596  )
597  )
598 
599 
600 # ------------------------------------------------------------------------------
601 # Kolmogorov Test
602 # ------------------------------------------------------------------------------
603 
604 
605 class KolmogorovTest(PvalueTest):
606  """ Kolmogorov-Smirnov Test """
607 
608  def __init__(self, *args, **kwargs):
609  """
610  Initialize Kolmogorov test.
611  @param args: See arguments of :class:`ComparisonBase`
612  @param kwargs: See arguments of :class:`ComparisonBase`
613  """
614  super().__init__(*args, **kwargs)
615 
616  def _compute(self):
617  """
618  Perform the actual test
619  @return: None
620  """
621  self._raise_has_correct_types()
622  self._raise_has_compatible_bins()
623 
624  local_object_a = self.object_a
625  local_object_b = self.object_b
626 
627  # very special handling for TEfficiencies
628  if self.object_a.ClassName() == "TEfficiency":
629  local_object_a = self._convert_teff_to_hist(self.object_a)
630  local_object_b = self._convert_teff_to_hist(self.object_b)
631  if self.debug:
632  print("Converting TEfficiency objects to histograms.")
633 
634  option_str = "UON"
635  if self.debug:
636  option_str += "D"
637 
638  self._pvalue = local_object_a.KolmogorovTest(local_object_b, option_str)
639 
640  def _get_comparison_result_long(self) -> str:
641  if self._pvalue is None:
642  return (
643  r"Could not perform Kolmogorov test between {{revision1}} "
644  r"and {{revision2}} due to an unknown error. Please submit "
645  r"a bug report."
646  )
647 
648  return (
649  r"Performed Komlogorov test between {{revision1}} "
650  r"and {{revision2}} "
651  r" <b>p-value: {pvalue:.6f}</b> (p-value warn: {pvalue_warn}, "
652  r"p-value error: {pvalue_error})".format(
653  pvalue=self._pvalue,
654  pvalue_warn=self._pvalue_warn,
655  pvalue_error=self._pvalue_error,
656  )
657  )
658 
659 
660 # ------------------------------------------------------------------------------
661 # Anderson-Darling Test
662 # ------------------------------------------------------------------------------
663 
664 
665 class AndersonDarlingTest(PvalueTest):
666  """ Anderson-Darling test"""
667 
668  def __init__(self, *args, **kwargs):
669  """
670  Initialize Anderson-Darling test.
671  @param args: See arguments of :class:`ComparisonBase`
672  @param kwargs: See arguments of :class:`ComparisonBase`
673  """
674  super().__init__(*args, **kwargs)
675 
676  def _compute(self):
677  """
678  Perform the actual test
679  @return: None
680  """
681  self._raise_has_correct_types()
682  # description on
683  # https://root.cern.ch/doc/master/classTH1.html#aa6b386786876dc304d73ab6b2606d4f6
684  # sounds like we don't have to have the same bins
685 
686  local_object_a = self.object_a
687  local_object_b = self.object_b
688 
689  # very special handling for TEfficiencies
690  if self.object_a.ClassName() == "TEfficiency":
691  local_object_a = self._convert_teff_to_hist(self.object_a)
692  local_object_b = self._convert_teff_to_hist(self.object_b)
693  if self.debug:
694  print("Converting TEfficiency objects to histograms.")
695 
696  option_str = ""
697  if self.debug:
698  option_str += "D"
699 
700  self._pvalue = local_object_a.AndersonDarlingTest(
701  local_object_b, option_str
702  )
703 
704  def _get_comparison_result_long(self) -> str:
705  if self._pvalue is None:
706  return (
707  r"Could not perform-Anderson Darling test between "
708  r"{{revision1}} and {{revision2}} due to an unknown error."
709  r" Please support a bug report."
710  )
711 
712  return (
713  r"Performed Anderson-Darling test between {{revision1}} "
714  r"and {{revision2}} "
715  r" <b>p-value: {pvalue:.6f}</b> (p-value warn: {pvalue_warn}, "
716  r"p-value error: {pvalue_error})".format(
717  pvalue=self._pvalue,
718  pvalue_warn=self._pvalue_warn,
719  pvalue_error=self._pvalue_error,
720  )
721  )
722 
723 
724 # ==============================================================================
725 # Helpers
726 # ==============================================================================
727 
728 
729 class TablePrinter:
730  """ A tiny class to print columns of fixed width numbers. """
731 
732  def __init__(self, ncols, width=None):
733  """
734  Constructor.
735  @param ncols: Number of columns
736  @param width: Width of each column. Either int or list.
737  """
738 
739  self.ncols = ncols
740  if not width:
741  width = 10
742  if isinstance(width, int):
743 
744  self.widths = [width] * ncols
745  elif isinstance(width, list) or isinstance(width, tuple):
746  # let's hope this is a list then.
747  self.widths = width
748 
749  @property
750  def tot_width(self):
751  """ Total width of the table """
752  width = 0
753  # the widths of each column
754  width += sum(self.widths)
755  # three characters between each two columns
756  width += (self.ncols - 1) * 3
757  # 2 characters at the very left and right
758  width += 2 * 2
759  return width
760 
761  def print_divider(self, char="="):
762  """ Print a divider made up from repeated chars """
763  print(char * self.tot_width)
764 
765  def print(self, cols):
766  """ Print one row """
767  assert len(cols) == self.ncols
768  out = []
769  for icol, col in enumerate(cols):
770  width = self.widths[icol]
771  if isinstance(col, int):
772  form = f"{{:{width}d}}"
773  out.append(form.format(col))
774  elif isinstance(col, float):
775  form = "{{:{}.{}f}}".format(width, width // 2)
776  out.append(form.format(col))
777  else:
778  # convert everything else to a string if it isn't already
779  col = str(col)
780  col = col[:width].rjust(width)
781  out.append(col)
782  print("| " + " | ".join(out) + " |")
783 
784 
785 def print_contents_and_errors(obj_a, obj_b):
786  """
787  Print contents, errors and chi2 deviation for each bin as well as
788  some other information about two TH1-like objects.
789  @param obj_a: First TH1-like object
790  @param obj_b: Second TH1-like object
791  @return: None
792  """
793  nbins = obj_a.GetNbinsX()
794 
795  total_a = sum([obj_a.GetBinContent(ibin) for ibin in range(0, nbins + 2)])
796  total_b = sum([obj_b.GetBinContent(ibin) for ibin in range(0, nbins + 2)])
797 
798  print(f"Total events/summed weights in object 1: {total_a:10.5f}")
799  print(f"Total events/summed weights in object 2: {total_b:10.5f}")
800 
801  chi2_tot = 0
802 
803  cp = TablePrinter(6)
804  print()
805  cp.print_divider()
806  cp.print(["ibin", "a", "err a", "b", "err b", "chi2"])
807  cp.print_divider()
808  for ibin in range(1, nbins + 1):
809  content_a = obj_a.GetBinContent(ibin)
810  content_b = obj_b.GetBinContent(ibin)
811  error_a = obj_a.GetBinError(ibin)
812  error_b = obj_b.GetBinError(ibin)
813  # This is implemented according to
814  # https://root.cern.ch/doc/master/classTH1.html
815  try:
816  chi2 = (total_b * content_a - total_a * content_b) ** 2 / (
817  total_b ** 2 * error_a ** 2 + total_a ** 2 * error_b ** 2
818  )
819  chi2_tot += chi2
820  except ZeroDivisionError:
821  chi2 = "nan"
822  cp.print([ibin, content_a, error_a, content_b, error_b, chi2])
823  cp.print_divider()
824  print()
825 
826  print(f"Total chi2: {chi2_tot:10.5f}")
827 
828 
829 # ==============================================================================
830 # Command Line Interface
831 # ==============================================================================
832 
833 
834 def debug_cli():
835  """ A small command line interface for debugging purposes. """
836 
837  # 1. Get command line arguments
838  # =============================
839 
840  desc = (
841  "For testing purposes: Run the chi2 comparison with objects from "
842  "two root files."
843  )
844  parser = argparse.ArgumentParser(desc)
845 
846  _ = "Rootfile to read the first object from"
847  parser.add_argument("rootfile_a", help=_)
848 
849  _ = "Name of object inside first rootfile."
850  parser.add_argument("name_a", help=_)
851 
852  _ = "Rootfile to read the second object from"
853  parser.add_argument("rootfile_b", help=_)
854 
855  _ = "Name of object inside second rootfile."
856  parser.add_argument("name_b", help=_)
857 
858  args = parser.parse_args()
859 
860  # 2. Open rootfiles and get objects
861  # =================================
862 
863  if not os.path.exists(args.rootfile_a):
864  raise ValueError(f"Could not find '{args.rootfile_a}'.")
865 
866  if not os.path.exists(args.rootfile_b):
867  raise ValueError(f"Could not find '{args.rootfile_b}'.")
868 
869  rootfile_a = ROOT.TFile(args.rootfile_a)
870  obj_a = rootfile_a.Get(args.name_a)
871  if not obj_a:
872  raise ValueError(
873  f"Could not find object '{args.name_a}' "
874  f"in file '{args.rootfile_a}'."
875  )
876 
877  rootfile_b = ROOT.TFile(args.rootfile_b)
878  obj_b = rootfile_b.Get(args.name_b)
879  if not obj_b:
880  raise ValueError(
881  f"Could not find object '{args.name_b}' "
882  f"in file '{args.rootfile_b}'."
883  )
884 
885  # 3. Performe testing with debug option
886  # =====================================
887 
888  test = Chi2Test(obj_a, obj_b, debug=True)
889  test.ensure_compute()
890 
891  print("If you see this message, then no exception was thrown.")
892 
893  # 4. Close files
894  # ==============
895 
896  rootfile_a.Close()
897  rootfile_b.Close()
898 
899 
900 if __name__ == "__main__":
901  # Run command line interface for testing purposes.
902  debug_cli()
903 
904 # End suppression of doxygen checks
905 # @endcond