Belle II Software  release-05-01-25
13_trackingEfficiency_createPlots.py
1 #!/usr/bin/env python3
2 # -*- coding: utf-8 -*-
3 
4 
13 
14 """
15 <header>
16  <input>trackingEfficiency_pt_0.05GeV.root,
17  trackingEfficiency_pt_0.10GeV.root,trackingEfficiency_pt_0.25GeV.root,
18  trackingEfficiency_pt_0.40GeV.root,trackingEfficiency_pt_0.60GeV.root,
19  trackingEfficiency_pt_1.00GeV.root,trackingEfficiency_pt_1.50GeV.root,
20  trackingEfficiency_pt_2.00GeV.root,trackingEfficiency_pt_3.00GeV.root,
21  trackingEfficiency_pt_4.00GeV.root</input>
22  <output>TrackingValidation.root</output>
23  <contact>software-tracking@belle2.org</contact>
24  <description>Create momentum resolution, impact parameter resolution and efficiency plots.</description>
25 </header>
26 """
27 
28 
29 import ROOT
30 
31 ROOT.PyConfig.IgnoreCommandLineOptions = True
32 from ROOT import TFile, TChain, TTree, TH1F, TCanvas, TGraphErrors, TGraph, \
33  gStyle, TNamed, TF1, TProfile
34 import sys
35 import math
36 import numpy as np
37 from optparse import OptionParser
39 
40 DELTA_PT = 0.0001
41 
42 PT_VALUES = get_generated_pt_values()
43 # PT_VALUES = [0.05, 0.1, 0.15, 0.25, 0.5, 1.0, 2.0, 3.0]
44 
45 # contact person information
46 # is added to the plot descriptions
47 CONTACT_PERSON = {'Name': 'tracking software mailing list',
48  'Email': 'software-tracking@belle2.org'}
49 
50 
51 def make_expert_plot(hist):
52  hist.GetListOfFunctions().Add(TNamed('MetaOptions', 'expert'))
53 
54 
55 def main():
56  """Function which is executed"""
57 
58  print('Tracking validation plots.')
59 
60  option_parser = OptionParser()
61  option_parser.add_option('-i', '--input-file', dest='input_file',
62  default='../trackingEfficiency_pt_*.root',
63  help='Root file with StandardTrackingPerformance output.'
64  )
65  option_parser.add_option('-o', '--output-file', dest='output_file',
66  default='TrackingValidation.root',
67  help='Root file with all histograms.')
68 
69  options = option_parser.parse_args()[0]
70 
71  gStyle.SetOptStat(0)
72 
73  # load data tree
74  tree_name = 'data'
75  data_tree = TChain(tree_name) # input_root_file.Get(tree_name)
76  data_tree.Add(options.input_file)
77 
78  number_entries = 0
79  try:
80  number_entries = data_tree.GetEntries()
81  except AttributeError:
82  print('Could not load tree with name %s.' % tree_name)
83 
84  if number_entries == 0:
85  print('Data tree \'%s\'is empty or does not exist. Exit.' % tree_name)
86  sys.exit(0)
87 
88  # output root file
89  output_file_name = options.output_file
90  output_root_file = TFile(output_file_name, 'recreate')
91 
92  # create efficiency in bins of pt plot
93  calculate_efficiency_in_pt(data_tree)
94 
95  pt_of_interest = [0.05, 0.25, 1.]
96  # pt_of_interest = PT_VALUES
97 
98  for pt_value in pt_of_interest:
99  # create plots of efficiency in bins of cos Theta for different pt
100  generate_cos_theta_plot(data_tree, pt_value)
101 
102  # create momentum resolution plot
103  create_momentum_resolution_plot(data_tree)
104 
105  draw_impact_parameter(data_tree)
106 
107  draw_pvalue(data_tree)
108 
109  draw_hit_counts(data_tree, pt_of_interest)
110 
111  draw_residua(data_tree, 'x', 'x_gen', pt_of_interest=pt_of_interest)
112  draw_residua(data_tree, 'y', 'y_gen', pt_of_interest=pt_of_interest)
113  draw_residua(data_tree, 'z', 'z_gen', pt_of_interest=pt_of_interest)
114 
115  draw_residua(
116  data_tree,
117  'pt',
118  'pt_gen',
119  bins=200,
120  ledge=-0.1,
121  uedge=0.1,
122  normalize=1,
123  pt_of_interest=pt_of_interest,
124  )
125 
126  # write additional plots in sub dir
127  sub_dir = 'additionalPlots'
128  output_root_file.mkdir(sub_dir)
129  output_root_file.cd(sub_dir)
130 
131  for pt_value in PT_VALUES:
132  # create plots of efficiency in bins of cos Theta for different pt
133  generate_cos_theta_plot(data_tree, pt_value)
134 
135  draw_residua(data_tree, 'x', 'x_gen')
136  draw_residua(data_tree, 'y', 'y_gen')
137  draw_residua(data_tree, 'z', 'z_gen')
138 
139  draw_residua(
140  data_tree,
141  'pt',
142  'pt_gen',
143  bins=200,
144  ledge=-0.1,
145  uedge=0.1,
146  normalize=1,
147  )
148 
149  # close output file
150  output_root_file.Write()
151  output_root_file.Close()
152 
153 
154 def draw_hit_counts(data_tree, pt_values):
155  """
156  Draw the hit count distribution.
157  """
158 
159  number_entries = data_tree.GetEntries()
160 
161  pt_lower = -0.025
162  pt_upper = 4.025
163  # 50 MeV wide bins. + 0.5 serves for correct rounding.
164  number_bins = int((pt_upper - pt_lower) / 0.05 + 0.5)
165 
166  hists = {}
167  for det in ['PXD', 'SVD', 'CDC']:
168  hists[det] = TProfile('hHitCounts%s' % det,
169  'Hit count profile for the %s;pT;nHits' % det,
170  number_bins, pt_lower, pt_upper)
171 
172  description = 'Distribution of Hit Counts in %s (Contact: %s).' \
173  % (det, CONTACT_PERSON['Email'])
174  check = ''
175 
176  hists[det].GetListOfFunctions().Add(TNamed('Description', description))
177  hists[det].GetListOfFunctions().Add(TNamed('Check', check))
178  make_expert_plot(hists[det])
179 
180  hNweights = TH1F('hNweights', 'number of weights stored', 201, 0, 201)
181  hWeightsProfile = TProfile('hWeightsProfile', 'profile of weights', 201,
182  0, 201)
183  hWeightsProfile.SetMinimum(0)
184  hWeightsProfile.SetMaximum(1.1)
185  hWeights = TH1F('hWeights', 'DAF weights', 50, 0, 1)
186 
187  for track in data_tree:
188  if not track.pValue == -999:
189  hists['PXD'].Fill(track.pt_gen, track.nPXDhits)
190  hists['SVD'].Fill(track.pt_gen, track.nSVDhits)
191  hists['CDC'].Fill(track.pt_gen, track.nCDChits)
192  hNweights.Fill(track.nWeights)
193  for i in range(track.nWeights):
194  hWeights.Fill(track.weights[i])
195  hWeightsProfile.Fill(i, track.weights[i])
196 
197  for key in hists:
198  hists[key].Write()
199  # hNweights.Write()
200  # hWeights.Write()
201  make_expert_plot(hWeightsProfile)
202  hWeightsProfile.Write()
203 
204 
205 def draw_pvalue(data_tree):
206  """
207  Create a histogram of the pValue of the track fits
208  """
209 
210  number_entries = data_tree.GetEntries()
211  # Normalize number of bins to number of events to keep plots comparable.
212  numBins = math.trunc(number_entries / 200)
213  hist_pvalue = TH1F('hpValue', 'p-Value Distribution', numBins, 0., 1.)
214  # Make axis range independent of first bin contents for more pleasant look
215  hist_pvalue.SetMinimum(0)
216  # hist_pvalue.SetMaximum(number_entries / 50.)
217 
218  for entry in range(number_entries):
219  data_tree.GetEntry(entry)
220 
221  try:
222  pvalue = data_tree.GetLeaf('pValue').GetValue()
223  except ReferenceError:
224  print('The variable "pValue" doesn\'t exist in the tree "%s".\nLeave this function without plotting the variable.'
225  % data_tree.GetName())
226  return
227 
228  if pvalue is -999:
229  continue
230  else:
231  hist_pvalue.Fill(pvalue)
232 
233  hist_pvalue.SetTitle('Distribution of pValue')
234  hist_pvalue.SetXTitle('pValue')
235  hist_pvalue.SetYTitle('number of entries')
236 
237  description = 'Distribution of pValues of the tracks (Contact: %s).' \
238  % CONTACT_PERSON['Email']
239  check = 'Should be a flat distribution.'
240 
241  hist_pvalue.GetListOfFunctions().Add(TNamed('Description', description))
242  hist_pvalue.GetListOfFunctions().Add(TNamed('Check', check))
243  make_expert_plot(hist_pvalue)
244 
245  hist_pvalue.Write()
246 
247 
248 def calculate_efficiency_in_pt(data_tree):
249  """Calculate single track reconstruction efficiency in bins of pt"""
250 
251  pt_lower = -0.025
252  pt_upper = 4.025
253  # 50 MeV wide bins. + 0.5 serves for correct rounding.
254  number_bins = int((pt_upper - pt_lower) / 0.05 + 0.5)
255 
256  hist_pt_gen = TH1F('hPtGen', 'hPtGen', number_bins, pt_lower, pt_upper)
257  hist_pt_rec = TH1F('hPtRec', 'hPtRec', number_bins, pt_lower, pt_upper)
258 
259  # draw data in defined histos
260  data_tree.Draw('pt_gen>>hPtGen', '', 'goff')
261  data_tree.Draw('pt_gen>>hPtRec', 'pt != -999', 'goff')
262 
263  # final hist
264  efficiency_hist = TH1F('hEfficiency', 'hEfficiency', number_bins,
265  pt_lower, pt_upper)
266  efficiency_hist.SetMinimum(0)
267  efficiency_hist.SetMaximum(1.05)
268 
269  description = ('Events with 10 muon tracks with fixed transverse '
270  'momentum are generated using the ParticleGun(500 '
271  'events for each pt value). The events are reconstructed '
272  'with VXDTF + CDC Track Finder + VXDCDCMerger.'
273  'plot shows the single track reconstruction efficiency '
274  'over the transverse momentum.')
275 
276  check = 'The efficiency should be stable for higher pt values.'
277  efficiency_hist.GetListOfFunctions().Add(TNamed('Description',
278  description))
279  efficiency_hist.GetListOfFunctions().Add(TNamed('Check', check))
280  efficiency_hist.GetListOfFunctions().Add(TNamed('Contact',
281  CONTACT_PERSON['Email']))
282 
283  # loop over bins and calculate efficiency and error of efficiency
284  for ibin in range(1, number_bins + 1):
285  number_generated = hist_pt_gen.GetBinContent(ibin)
286  efficiency = 0
287  error = 0
288 
289  if number_generated > 0:
290  number_reconstructed = hist_pt_rec.GetBinContent(ibin)
291  efficiency = number_reconstructed / number_generated
292  error = math.sqrt(number_reconstructed * (number_generated - number_reconstructed) / pow(number_generated, 3))
293  # set according bin to the efficiency value and add errorbar
294  efficiency_hist.SetBinContent(ibin, efficiency)
295  efficiency_hist.SetBinError(ibin, error)
296 
297  efficiency_hist.SetTitle('Tracking efficiency in bins of transverse momentum pt.'
298  )
299  efficiency_hist.SetXTitle('pt in GeV')
300  efficiency_hist.SetYTitle('efficiency')
301 
302  # write hist to the output file
303  efficiency_hist.Write()
304 
305 
306 def generate_cos_theta_plot(data_tree, pt_value):
307  """Creates a efficiency histo in bins of cos theta"""
308 
309  number_bins = 100
310  cos_lower = -1
311  cos_upper = 1
312 
313  hist_cos_gen = TH1F('hCosGen', 'hCosGen', number_bins, cos_lower,
314  cos_upper)
315  hist_cos_rec = TH1F('hCosRec', 'hCosRec', number_bins, cos_lower,
316  cos_upper)
317 
318  data_tree.Draw('cosTheta_gen>>hCosGen',
319  'pt_gen>(%.2f - %f) &&pt_gen<(%.2f + %f)' % (pt_value,
320  DELTA_PT, pt_value, DELTA_PT), 'goff')
321  data_tree.Draw('cosTheta_gen>>hCosRec',
322  'pt_gen>(%.2f - %f) &&pt_gen<(%.2f + %f) && pt != -999'
323  % (pt_value, DELTA_PT, pt_value, DELTA_PT), 'goff')
324 
325  description = ('Events with 10 muon tracks with fixed transverse '
326  'momentum are generated using the ParticleGun(500 '
327  'events for each pt value). The events are reconstructed '
328  'with VXDTF + CDCTF + MCTrackCandCombiner. The plot '
329  'shows the single track reconstruction efficiency in '
330  'bins of the polar angle for the fixed transverse '
331  'momentum pt = %.2f GeV.')
332  check = 'Stable efficiency over the hole range of the polar angle.'
333 
334  efficiency_hist = TH1F('hEfficiencyPt%.2fGeV' % pt_value,
335  'hEfficiencyPt%.2fGeV' % pt_value, number_bins,
336  cos_lower, cos_upper)
337  efficiency_hist.GetListOfFunctions().Add(TNamed('Description',
338  description))
339  efficiency_hist.GetListOfFunctions().Add(TNamed('Check', check))
340  efficiency_hist.GetListOfFunctions().Add(TNamed('Contact',
341  CONTACT_PERSON['Email']))
342 
343  for ibin in range(1, number_bins + 1):
344  efficiency = 0
345  error = 0
346  number_generated = hist_cos_gen.GetBinContent(ibin)
347 
348  if number_generated > 0:
349  number_reconstructed = hist_cos_rec.GetBinContent(ibin)
350  efficiency = number_reconstructed / number_generated
351  error = math.sqrt(number_reconstructed * (number_generated - number_reconstructed) / pow(number_generated, 3))
352 
353  efficiency_hist.SetBinContent(ibin, efficiency)
354  efficiency_hist.SetBinError(ibin, error)
355 
356  efficiency_hist.SetTitle('Tracks with pt = %.2f GeV' % pt_value)
357  efficiency_hist.SetXTitle('cos #Theta')
358  efficiency_hist.SetYTitle('efficiency')
359  make_expert_plot(efficiency_hist)
360  efficiency_hist.Write()
361 
362 
363 def create_momentum_resolution_plot(data_tree):
364  """Create momentum resolution plot"""
365 
366  pt_lower = -0.025
367  pt_upper = 4.025
368  # 50 MeV wide bins. + 0.5 serves for correct rounding.
369  number_bins = int((pt_upper - pt_lower) / 0.05 + 0.5)
370 
371  sigma_pt_values = calculate_momentum_resolution2(data_tree)
372 
373  hist_resolution = TH1F('hPtResolution', 'hPtResolution', number_bins,
374  pt_lower, pt_upper)
375 
376  fit_pt_values = []
377  fit_pt_res_values = []
378  fit_pt_res_errors = []
379 
380  for (pt_value, sigma_pt_value) in sigma_pt_values.items():
381  ibin = hist_resolution.FindBin(pt_value)
382  bin_center = hist_resolution.GetBinCenter(ibin)
383 
384  sigma_pt_over_pt = 0
385  if sigma_pt_value[0] != -999:
386  sigma_pt_over_pt = sigma_pt_value[0] / pt_value
387  sigma_pt_error_over_pt = sigma_pt_value[1] / pt_value
388 
389  fit_pt_res_values.append(sigma_pt_over_pt)
390  fit_pt_res_errors.append(sigma_pt_error_over_pt)
391  fit_pt_values.append(bin_center)
392 
393  hist_resolution.SetBinContent(ibin, sigma_pt_over_pt)
394  hist_resolution.SetBinError(ibin, sigma_pt_error_over_pt)
395 
396  hist_resolution.SetTitle('Transverse Momentum resolution')
397  hist_resolution.SetXTitle('pt in GeV/c')
398  hist_resolution.SetYTitle('#sigma_{pt}/pt')
399 
400  description = ('Events with 10 muon tracks with fixed transverse '
401  'momentum are generated using the ParticleGun(200 '
402  'events for each pt value). The events are reconstructed '
403  'with VXDTF + CDCTF + MCTrackCandCombiner. The plot '
404  'shows the relative momentum resolution of the '
405  'transverse momentum over transverse momentum.')
406 
407  check = ''
408  hist_resolution.GetListOfFunctions().Add(TNamed('Description',
409  description))
410  hist_resolution.GetListOfFunctions().Add(TNamed('Check', check))
411  hist_resolution.GetListOfFunctions().Add(TNamed('Contact',
412  CONTACT_PERSON['Email']))
413  hist_resolution.GetListOfFunctions().Add(TNamed('MetaOptions', 'logy'))
414  hist_resolution.Write()
415 
416  # fit_parameters = fit_resolution(fit_pt_values, fit_pt_res_values, 1, fit_pt_res_errors)
417 
418 
419 def get_scaled_rms_90(values, scale_factor=0.79):
420  """
421  Calculate the RMS90
422 
423  param values: list of numbers
424  return scaled RMS90 of the distribution in values
425  """
426 
427  number_entries = len(values)
428 
429  (vMin, vMax) = (np.nanpercentile(values, 5), np.nanpercentile(values, 95))
430  final_list = np.array(values)[np.logical_and(values > vMin, values < vMax)]
431 
432  rms_90 = np.sqrt(np.mean(np.square(final_list))) / scale_factor
433 
434  rms_error = rms_90 / (scale_factor * np.sqrt(2 * len(final_list)))
435 
436  return (rms_90, rms_error)
437 
438 
439 def value_in_list(test_value, value_list):
440  """
441  Function checks, if a given value is in the list values
442  """
443 
444  for value in value_list:
445  if value - DELTA_PT < test_value < value + DELTA_PT:
446  return (True, value)
447 
448  return (False, -999)
449 
450 
451 def calculate_momentum_resolution2(data_tree):
452  """
453  Calculate momentum resolution
454  """
455 
456  delta_pt_dict = {}
457 
458  for key in PT_VALUES:
459  delta_pt_dict[key] = []
460 
461  for entry_index in range(data_tree.GetEntries()):
462  # load entry
463  data_tree.GetEntry(entry_index)
464 
465  pt_gen = data_tree.GetLeaf('pt_gen').GetValue()
466 
467  test_object = value_in_list(pt_gen, PT_VALUES)
468 
469  if test_object[0]:
470  pt = data_tree.GetLeaf('pt').GetValue()
471  if pt != -999:
472  try:
473  delta_pt_dict[test_object[1]].append(pt - pt_gen)
474  except KeyError:
475  print('pt = %0.2f is not generated and not used as key. Abort!'
476  % test_object[1])
477  sys.exit(1)
478  pt_resolutions = {}
479 
480  print('********************************')
481  print('Momentum resolution:')
482  for key in sorted(delta_pt_dict):
483  (rms90, rms_error) = get_scaled_rms_90(delta_pt_dict[key])
484  pt_resolutions[key] = [rms90, rms_error]
485  print('pt = %0.2f: sigma/pt = %0.4f' % (key, rms90 / key))
486 
487  return pt_resolutions
488 
489 
490 def fit_resolution(
491  x_value_list,
492  y_value_list,
493  degree,
494  y_value_errors=None,
495 ):
496  """
497  Fit a polynomial to data
498 
499  @param x_value_list list of x values
500 
501  @param y_value_list list of y values
502 
503  @param degree degree of the fitting polynomial
504 
505  @param y_value_errors list of corresponding errors on y x_values
506  """
507 
508  try:
509  import scipy.optimize
510  except ImportError:
511  print('Module scipy.optimize is not installed. Fit cannot be done.')
512  return []
513 
514  # weights = None
515  # if y_value_errors is not None:
516  # weights = [1 / x ** 2 for x in y_value_errors]
517  #
518  # parameters = np.polyfit(x_value_list, y_value_list, degree, w=weights,
519  # cov=True)
520 
521  x_value_list = np.array(x_value_list)
522  y_value_list = np.array(y_value_list)
523  y_value_errors = np.array(y_value_errors)
524 
525  (par_opt, par_cov) = \
526  scipy.optimize.curve_fit(fit_function_sigma_pt_over_pt, x_value_list,
527  y_value_list, [0.002, 0.003],
528  sigma=y_value_errors)
529 
530  print(par_opt, par_cov)
531 
532  return par_opt
533 
534 
535 def draw_impact_parameter(data_tree):
536  """
537  Create a histogram with the impact parameter resolution
538 
539  @param data_tree ROOT tree with the input data
540  """
541 
542  print('********************************')
543  print('Impact parameter:')
544 
545  impact_param_d = {}
546  impact_param_z = {}
547  for key in PT_VALUES:
548  impact_param_d[key] = []
549  impact_param_z[key] = []
550 
551  try:
552  number_of_entries = data_tree.GetEntries()
553  except AttributeError:
554  print('ERROR: data_tree in function "draw_impact_parameter" is no TTree.')
555  sys.exit(2)
556 
557  for index_entry in range(number_of_entries):
558  data_tree.GetEntry(index_entry)
559 
560  pt_gen = data_tree.GetLeaf('pt_gen').GetValue()
561 
562  test_object = value_in_list(pt_gen, PT_VALUES)
563  if test_object[0] is True:
564  poca_x = data_tree.GetLeaf('x').GetValue()
565 
566  if poca_x == -999:
567  continue
568 
569  poca_y = data_tree.GetLeaf('y').GetValue()
570  poca_z = data_tree.GetLeaf('z').GetValue()
571 
572  z_gen = data_tree.GetLeaf('z_gen').GetValue()
573 
574  px = data_tree.GetLeaf('px').GetValue()
575  py = data_tree.GetLeaf('py').GetValue()
576 
577  d0 = d0_sign([poca_x, poca_y], [px, py]) * np.hypot(poca_x, poca_y)
578 
579  impact_param_d[test_object[1]].append(d0)
580  impact_param_z[test_object[1]].append(poca_z - z_gen)
581  # estimate sign of d0
582 
583  d0_resolutions = {}
584  z_resolutions = {}
585  for key in sorted(impact_param_d):
586  (d0_resolution, d0_error) = get_scaled_rms_90(impact_param_d[key])
587  d0_resolutions[key] = [d0_resolution, d0_error]
588  (z_resolution, z_error) = get_scaled_rms_90(impact_param_z[key])
589  z_resolutions[key] = [z_resolution, z_error]
590  print('pt = %0.2f: sigma_d0 = %0.4f, sigma_z = %0.4f' % (key,
591  d0_resolution, z_resolution))
592 
593  number_bins = 62
594  lower_edge = -0.025
595  upper_edge = 4.025
596 
597  hist_impact_parameter_d0 = TH1F('hd0resolution', 'hd0resolution',
598  number_bins, lower_edge, upper_edge)
599  hist_impact_parameter_d0.SetTitle('d0 resolution')
600  hist_impact_parameter_d0.SetXTitle('pt in GeV/c')
601  hist_impact_parameter_d0.SetYTitle('#sigma_{d0} [cm]')
602  hist_impact_parameter_d0.GetListOfFunctions().Add(TNamed('MetaOptions', 'logy'))
603  hist_impact_parameter_d0.GetListOfFunctions().Add(TNamed('Contact',
604  CONTACT_PERSON['Email']))
605 
606  hist_impact_parameter_z = TH1F('hzresolution', 'hzresolution',
607  number_bins, lower_edge, upper_edge)
608  hist_impact_parameter_z.SetTitle('z resolution')
609  hist_impact_parameter_z.SetXTitle('pt in GeV/c')
610  hist_impact_parameter_z.SetYTitle('#sigma_{z} [cm]')
611  hist_impact_parameter_z.GetListOfFunctions().Add(TNamed('MetaOptions', 'logy'))
612  hist_impact_parameter_z.GetListOfFunctions().Add(TNamed('Contact',
613  CONTACT_PERSON['Email']))
614 
615  for pt in PT_VALUES:
616  ibin = hist_impact_parameter_d0.FindBin(pt)
617  if d0_resolutions[pt][0] > 0 and d0_resolutions[pt][1] > 0:
618  hist_impact_parameter_d0.SetBinContent(ibin, d0_resolutions[pt][0])
619  hist_impact_parameter_d0.SetBinError(ibin, d0_resolutions[pt][1])
620  if z_resolutions[pt][0] > 0 and z_resolutions[pt][1] > 0:
621  hist_impact_parameter_z.SetBinContent(ibin, z_resolutions[pt][0])
622  hist_impact_parameter_z.SetBinError(ibin, z_resolutions[pt][1])
623 
624  hist_impact_parameter_d0.Write()
625  hist_impact_parameter_z.Write()
626 
627 
628 def scale_histogram(hist):
629  """
630  Scales a TH1F using the TH1::Scale function
631  """
632 
633  if isinstance(hist, TH1F):
634  if hist.GetEntries() > 0:
635  hist.Scale(1 / hist.GetEntries())
636  else:
637  print('Cannot scale. No entries in histogram ' + hist.GetName())
638  else:
639  print('Not a TH1F. Not able to scale.')
640 
641 
642 def d0_sign(vector_d, momentum):
643  """
644  Estimate the sign of the impact parameter d0_sign
645 
646  @return -1 or +1
647  """
648 
649  if len(vector_d) != 2 or len(momentum) != 2:
650  print('Need x and y values of d and p_d.')
651  sys.exit(1)
652 
653  phi_momentum = np.arctan2(momentum[1], momentum[0])
654  phi_d = np.arctan2(vector_d[1], vector_d[0])
655 
656  diff = phi_momentum - phi_d
657  # Catch corner cases. For clarity, don't return early.
658  if diff > np.pi:
659  diff = diff - 2 * np.pi
660  elif diff < -np.pi:
661  diff = 2 * np.pi + diff
662 
663  return np.sign(diff)
664 
665 
666 def draw_residua(
667  data_tree,
668  variable_name,
669  gen_variable_name,
670  bins=100,
671  ledge=-0.01,
672  uedge=0.01,
673  normalize=0,
674  pt_of_interest=None,
675 ):
676  """
677  Write histogram of difference of variable_name (reconstructed) and
678  gen_variable_name (generated) in gDirectory
679  """
680 
681  if pt_of_interest is None:
682  used_pts = PT_VALUES
683  else:
684  used_pts = [pt for pt in pt_of_interest if pt in PT_VALUES]
685 
686  number_entries = data_tree.GetEntries()
687 
688  histograms = {}
689  for pt_value in used_pts:
690  histograms[pt_value] = TH1F('h%sResiduum_%0.2fGeV' % (variable_name,
691  pt_value), 'h%sResiduum_%0.2fGeV'
692  % (variable_name, pt_value), bins, ledge,
693  uedge)
694 
695  for index in range(number_entries):
696  data_tree.GetEntry(index)
697 
698  variable_reconstructed = data_tree.GetLeaf(variable_name).GetValue()
699  if variable_reconstructed != -999:
700  variable_generated = \
701  data_tree.GetLeaf(gen_variable_name).GetValue()
702  pt_gen = round(data_tree.GetLeaf('pt_gen').GetValue(), 2)
703  if pt_gen in used_pts:
704  if normalize:
705  histograms[pt_gen].Fill((variable_reconstructed - variable_generated) / variable_generated)
706  else:
707  histograms[pt_gen].Fill(variable_reconstructed - variable_generated)
708 
709  for (pt, hist) in histograms.items():
710  scale_histogram(hist)
711  if normalize:
712  hist.SetXTitle('(%s - %s) / %s' % (variable_name,
713  gen_variable_name, gen_variable_name))
714  else:
715  hist.SetXTitle('(%s - %s)' % (variable_name, gen_variable_name))
716  hist.GetListOfFunctions().Add(TNamed('Contact', CONTACT_PERSON['Email'
717  ]))
718  make_expert_plot(hist)
719  hist.Write()
720 
721 
722 def tf1_fit_function(x, par):
723  return fit_function_sigma_pt_over_pt(x, par[0], par[1])
724 
725 
726 def fit_function_sigma_pt_over_pt(x, a, b):
727  return np.sqrt((a * x) ** 2 + b * b)
728 
729 
730 if __name__ == '__main__':
731  main()
tracking.validation.tracking_efficiency_helpers
Definition: tracking_efficiency_helpers.py:1
main
int main(int argc, char **argv)
Run all tests.
Definition: test_main.cc:77