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>
31 ROOT.PyConfig.IgnoreCommandLineOptions =
True
32 from ROOT
import TFile, TChain, TTree, TH1F, TCanvas, TGraphErrors, TGraph, \
33 gStyle, TNamed, TF1, TProfile
37 from optparse
import OptionParser
42 PT_VALUES = get_generated_pt_values()
47 CONTACT_PERSON = {
'Name':
'tracking software mailing list',
48 'Email':
'software-tracking@belle2.org'}
51 def make_expert_plot(hist):
52 hist.GetListOfFunctions().Add(TNamed(
'MetaOptions',
'expert'))
56 """Function which is executed"""
58 print(
'Tracking validation plots.')
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.'
65 option_parser.add_option(
'-o',
'--output-file', dest=
'output_file',
66 default=
'TrackingValidation.root',
67 help=
'Root file with all histograms.')
69 options = option_parser.parse_args()[0]
75 data_tree = TChain(tree_name)
76 data_tree.Add(options.input_file)
80 number_entries = data_tree.GetEntries()
81 except AttributeError:
82 print(
'Could not load tree with name %s.' % tree_name)
84 if number_entries == 0:
85 print(
'Data tree \'%s\'is empty or does not exist. Exit.' % tree_name)
89 output_file_name = options.output_file
90 output_root_file = TFile(output_file_name,
'recreate')
93 calculate_efficiency_in_pt(data_tree)
95 pt_of_interest = [0.05, 0.25, 1.]
98 for pt_value
in pt_of_interest:
100 generate_cos_theta_plot(data_tree, pt_value)
103 create_momentum_resolution_plot(data_tree)
105 draw_impact_parameter(data_tree)
107 draw_pvalue(data_tree)
109 draw_hit_counts(data_tree, pt_of_interest)
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)
123 pt_of_interest=pt_of_interest,
127 sub_dir =
'additionalPlots'
128 output_root_file.mkdir(sub_dir)
129 output_root_file.cd(sub_dir)
131 for pt_value
in PT_VALUES:
133 generate_cos_theta_plot(data_tree, pt_value)
135 draw_residua(data_tree,
'x',
'x_gen')
136 draw_residua(data_tree,
'y',
'y_gen')
137 draw_residua(data_tree,
'z',
'z_gen')
150 output_root_file.Write()
151 output_root_file.Close()
154 def draw_hit_counts(data_tree, pt_values):
156 Draw the hit count distribution.
159 number_entries = data_tree.GetEntries()
164 number_bins = int((pt_upper - pt_lower) / 0.05 + 0.5)
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)
172 description =
'Distribution of Hit Counts in %s (Contact: %s).' \
173 % (det, CONTACT_PERSON[
'Email'])
176 hists[det].GetListOfFunctions().Add(TNamed(
'Description', description))
177 hists[det].GetListOfFunctions().Add(TNamed(
'Check', check))
178 make_expert_plot(hists[det])
180 hNweights = TH1F(
'hNweights',
'number of weights stored', 201, 0, 201)
181 hWeightsProfile = TProfile(
'hWeightsProfile',
'profile of weights', 201,
183 hWeightsProfile.SetMinimum(0)
184 hWeightsProfile.SetMaximum(1.1)
185 hWeights = TH1F(
'hWeights',
'DAF weights', 50, 0, 1)
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])
201 make_expert_plot(hWeightsProfile)
202 hWeightsProfile.Write()
205 def draw_pvalue(data_tree):
207 Create a histogram of the pValue of the track fits
210 number_entries = data_tree.GetEntries()
212 numBins = math.trunc(number_entries / 200)
213 hist_pvalue = TH1F(
'hpValue',
'p-Value Distribution', numBins, 0., 1.)
215 hist_pvalue.SetMinimum(0)
218 for entry
in range(number_entries):
219 data_tree.GetEntry(entry)
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())
231 hist_pvalue.Fill(pvalue)
233 hist_pvalue.SetTitle(
'Distribution of pValue')
234 hist_pvalue.SetXTitle(
'pValue')
235 hist_pvalue.SetYTitle(
'number of entries')
237 description =
'Distribution of pValues of the tracks (Contact: %s).' \
238 % CONTACT_PERSON[
'Email']
239 check =
'Should be a flat distribution.'
241 hist_pvalue.GetListOfFunctions().Add(TNamed(
'Description', description))
242 hist_pvalue.GetListOfFunctions().Add(TNamed(
'Check', check))
243 make_expert_plot(hist_pvalue)
248 def calculate_efficiency_in_pt(data_tree):
249 """Calculate single track reconstruction efficiency in bins of pt"""
254 number_bins = int((pt_upper - pt_lower) / 0.05 + 0.5)
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)
260 data_tree.Draw(
'pt_gen>>hPtGen',
'',
'goff')
261 data_tree.Draw(
'pt_gen>>hPtRec',
'pt != -999',
'goff')
264 efficiency_hist = TH1F(
'hEfficiency',
'hEfficiency', number_bins,
266 efficiency_hist.SetMinimum(0)
267 efficiency_hist.SetMaximum(1.05)
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.')
276 check =
'The efficiency should be stable for higher pt values.'
277 efficiency_hist.GetListOfFunctions().Add(TNamed(
'Description',
279 efficiency_hist.GetListOfFunctions().Add(TNamed(
'Check', check))
280 efficiency_hist.GetListOfFunctions().Add(TNamed(
'Contact',
281 CONTACT_PERSON[
'Email']))
284 for ibin
in range(1, number_bins + 1):
285 number_generated = hist_pt_gen.GetBinContent(ibin)
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))
294 efficiency_hist.SetBinContent(ibin, efficiency)
295 efficiency_hist.SetBinError(ibin, error)
297 efficiency_hist.SetTitle(
'Tracking efficiency in bins of transverse momentum pt.'
299 efficiency_hist.SetXTitle(
'pt in GeV')
300 efficiency_hist.SetYTitle(
'efficiency')
303 efficiency_hist.Write()
306 def generate_cos_theta_plot(data_tree, pt_value):
307 """Creates a efficiency histo in bins of cos theta"""
313 hist_cos_gen = TH1F(
'hCosGen',
'hCosGen', number_bins, cos_lower,
315 hist_cos_rec = TH1F(
'hCosRec',
'hCosRec', number_bins, cos_lower,
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')
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.'
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',
339 efficiency_hist.GetListOfFunctions().Add(TNamed(
'Check', check))
340 efficiency_hist.GetListOfFunctions().Add(TNamed(
'Contact',
341 CONTACT_PERSON[
'Email']))
343 for ibin
in range(1, number_bins + 1):
346 number_generated = hist_cos_gen.GetBinContent(ibin)
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))
353 efficiency_hist.SetBinContent(ibin, efficiency)
354 efficiency_hist.SetBinError(ibin, error)
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()
363 def create_momentum_resolution_plot(data_tree):
364 """Create momentum resolution plot"""
369 number_bins = int((pt_upper - pt_lower) / 0.05 + 0.5)
371 sigma_pt_values = calculate_momentum_resolution2(data_tree)
373 hist_resolution = TH1F(
'hPtResolution',
'hPtResolution', number_bins,
377 fit_pt_res_values = []
378 fit_pt_res_errors = []
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)
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
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)
393 hist_resolution.SetBinContent(ibin, sigma_pt_over_pt)
394 hist_resolution.SetBinError(ibin, sigma_pt_error_over_pt)
396 hist_resolution.SetTitle(
'Transverse Momentum resolution')
397 hist_resolution.SetXTitle(
'pt in GeV/c')
398 hist_resolution.SetYTitle(
'#sigma_{pt}/pt')
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.')
408 hist_resolution.GetListOfFunctions().Add(TNamed(
'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()
419 def get_scaled_rms_90(values, scale_factor=0.79):
423 param values: list of numbers
424 return scaled RMS90 of the distribution in values
427 number_entries = len(values)
429 (vMin, vMax) = (np.nanpercentile(values, 5), np.nanpercentile(values, 95))
430 final_list = np.array(values)[np.logical_and(values > vMin, values < vMax)]
432 rms_90 = np.sqrt(np.mean(np.square(final_list))) / scale_factor
434 rms_error = rms_90 / (scale_factor * np.sqrt(2 * len(final_list)))
436 return (rms_90, rms_error)
439 def value_in_list(test_value, value_list):
441 Function checks, if a given value is in the list values
444 for value
in value_list:
445 if value - DELTA_PT < test_value < value + DELTA_PT:
451 def calculate_momentum_resolution2(data_tree):
453 Calculate momentum resolution
458 for key
in PT_VALUES:
459 delta_pt_dict[key] = []
461 for entry_index
in range(data_tree.GetEntries()):
463 data_tree.GetEntry(entry_index)
465 pt_gen = data_tree.GetLeaf(
'pt_gen').GetValue()
467 test_object = value_in_list(pt_gen, PT_VALUES)
470 pt = data_tree.GetLeaf(
'pt').GetValue()
473 delta_pt_dict[test_object[1]].append(pt - pt_gen)
475 print(
'pt = %0.2f is not generated and not used as key. Abort!'
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))
487 return pt_resolutions
497 Fit a polynomial to data
499 @param x_value_list list of x values
501 @param y_value_list list of y values
503 @param degree degree of the fitting polynomial
505 @param y_value_errors list of corresponding errors on y x_values
509 import scipy.optimize
511 print(
'Module scipy.optimize is not installed. Fit cannot be done.')
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)
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)
530 print(par_opt, par_cov)
535 def draw_impact_parameter(data_tree):
537 Create a histogram with the impact parameter resolution
539 @param data_tree ROOT tree with the input data
542 print(
'********************************')
543 print(
'Impact parameter:')
547 for key
in PT_VALUES:
548 impact_param_d[key] = []
549 impact_param_z[key] = []
552 number_of_entries = data_tree.GetEntries()
553 except AttributeError:
554 print(
'ERROR: data_tree in function "draw_impact_parameter" is no TTree.')
557 for index_entry
in range(number_of_entries):
558 data_tree.GetEntry(index_entry)
560 pt_gen = data_tree.GetLeaf(
'pt_gen').GetValue()
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()
569 poca_y = data_tree.GetLeaf(
'y').GetValue()
570 poca_z = data_tree.GetLeaf(
'z').GetValue()
572 z_gen = data_tree.GetLeaf(
'z_gen').GetValue()
574 px = data_tree.GetLeaf(
'px').GetValue()
575 py = data_tree.GetLeaf(
'py').GetValue()
577 d0 = d0_sign([poca_x, poca_y], [px, py]) * np.hypot(poca_x, poca_y)
579 impact_param_d[test_object[1]].append(d0)
580 impact_param_z[test_object[1]].append(poca_z - z_gen)
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))
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']))
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']))
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])
624 hist_impact_parameter_d0.Write()
625 hist_impact_parameter_z.Write()
628 def scale_histogram(hist):
630 Scales a TH1F using the TH1::Scale function
633 if isinstance(hist, TH1F):
634 if hist.GetEntries() > 0:
635 hist.Scale(1 / hist.GetEntries())
637 print(
'Cannot scale. No entries in histogram ' + hist.GetName())
639 print(
'Not a TH1F. Not able to scale.')
642 def d0_sign(vector_d, momentum):
644 Estimate the sign of the impact parameter d0_sign
649 if len(vector_d) != 2
or len(momentum) != 2:
650 print(
'Need x and y values of d and p_d.')
653 phi_momentum = np.arctan2(momentum[1], momentum[0])
654 phi_d = np.arctan2(vector_d[1], vector_d[0])
656 diff = phi_momentum - phi_d
659 diff = diff - 2 * np.pi
661 diff = 2 * np.pi + diff
677 Write histogram of difference of variable_name (reconstructed) and
678 gen_variable_name (generated) in gDirectory
681 if pt_of_interest
is None:
684 used_pts = [pt
for pt
in pt_of_interest
if pt
in PT_VALUES]
686 number_entries = data_tree.GetEntries()
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,
695 for index
in range(number_entries):
696 data_tree.GetEntry(index)
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:
705 histograms[pt_gen].Fill((variable_reconstructed - variable_generated) / variable_generated)
707 histograms[pt_gen].Fill(variable_reconstructed - variable_generated)
709 for (pt, hist)
in histograms.items():
710 scale_histogram(hist)
712 hist.SetXTitle(
'(%s - %s) / %s' % (variable_name,
713 gen_variable_name, gen_variable_name))
715 hist.SetXTitle(
'(%s - %s)' % (variable_name, gen_variable_name))
716 hist.GetListOfFunctions().Add(TNamed(
'Contact', CONTACT_PERSON[
'Email'
718 make_expert_plot(hist)
722 def tf1_fit_function(x, par):
723 return fit_function_sigma_pt_over_pt(x, par[0], par[1])
726 def fit_function_sigma_pt_over_pt(x, a, b):
727 return np.sqrt((a * x) ** 2 + b * b)
730 if __name__ ==
'__main__':