13Validation of the Invariant Collision Energy calibration
17from prompt
import ValidationSettings
26import matplotlib.pyplot
as plt
27from matplotlib
import rcParams
33from datetime
import datetime
35from collections
import OrderedDict
36from shutil
import copyfile, rmtree
40settings = ValidationSettings(name=
'eCMS Calibrations',
42 download_files=[
'stdout'],
46def setPlotRange(df, tag):
48 Function which adjusts the y-axis range of the plot according to tag
52 df4S = df[df[
'Ecms'] > 10560][
'Ecms']
55 plt.ylim(bottom=yMin-3)
57 dfHigh = df[df[
'Ecms'] < 10600][
'Ecms']
63 dfOff = df[df[
'Ecms'] < 10560][
'Ecms']
71 Converts time from UTC to the JST
73 return np.vectorize(
lambda t: datetime.utcfromtimestamp((t + 9) * 3600))(times)
76def plotSplitLines(dfC):
78 Plot vertical lines in places where run energy type changes
80 indx = np.where(np.diff(dfC['pull'] == 0))[0]
82 tt1 = toJST(dfC[
't2'].loc[i]).item()
83 tt2 = toJST(dfC[
't1'].loc[i+1]).item()
84 tt = tt1 + (tt2 - tt1) / 2
85 plt.axvline(x=tt, color=
'g', linestyle=
'--')
88def get_Ecms_values(path):
90 Load the values of the Ecms properties from the payloads into the list
94 with open(path +
'/database.txt')
as fDB:
100 fN = ll[0].replace(
'/',
'_') +
'_rev_' + ll[1] +
'.root'
103 runDict[fN] = ((int(r[0]), int(r[1])), (int(r[2]), int(r[3])))
107 fList = glob(path +
'/*.root')
109 bName = os.path.basename(fName)
110 runExp = runDict[bName]
112 f = ROOT.TFile.Open(fName)
113 ecmsObj = f.Get(
"CollisionInvariantMass")
114 assert(ecmsObj.ClassName() ==
"Belle2::CollisionInvariantMass")
116 eCMS = ecmsObj.getMass()
117 eCMSe = ecmsObj.getMassError()
118 eCMSs = ecmsObj.getMassSpread()
120 arr.append((runExp, (eCMS, eCMSe, eCMSs)))
124 arr = sorted(arr, key=
lambda x: x[0])
130 def __init__(self, location):
132 Data are loaded from text files to pandas
136 self.dfB = pd.read_csv(f
'{location}/BonlyEcmsCalib.txt', delim_whitespace=
True)
139 self.dfC = pd.read_csv(f
'{location}/finalEcmsCalib.txt', delim_whitespace=
True)
142 self.dfM = pd.read_csv(f
'{location}/mumuEcalib.txt', delim_whitespace=
True)
150 for i
in range(len(dfM)):
151 if (dfM[
'id'][i]) == 0:
154 dfM.at[i,
'type'] = state
155 dfC.at[i,
'type'] = state
158 def plotLine(df, var, color, label):
160 Plot single line with error band
164 nan = np.full(len(df), np.nan)
165 avg = (df[
't1']+df[
't2']).to_numpy()/2
166 times = np.c_[df[[
't1',
't2']].to_numpy(), avg].ravel()
167 eCMS = np.c_[df[[var, var]].to_numpy(), nan].ravel()
168 eCMSu = np.c_[df[[varUnc, varUnc]].to_numpy(), nan].ravel()
170 timesg = np.c_[df[
't1'].to_numpy(), avg, df[
't2'].to_numpy()].ravel()
171 eCMSg = np.c_[df[var].to_numpy(), nan, df[var].to_numpy()].ravel()
174 timesg = toJST(timesg)
176 plt.fill_between(times, eCMS-eCMSu, eCMS+eCMSu, alpha=0.2, color=color)
177 plt.plot(times, eCMS, linewidth=2, color=color, label=label)
178 plt.plot(timesg, eCMSg, linewidth=2, color=color, alpha=0.35)
180 def plotEcmsComparison(self, limits=None, tag=''):
182 Plot Ecms obtained from combined, hadB
and mumu method
185 Plotter.plotLine(self.dfB, 'Ecms',
'green', label=
'B decay method')
186 Plotter.plotLine(self.dfC,
'Ecms',
'red', label=
'Combined method')
188 d = np.nanmedian(self.dfC[
'Ecms']-self.dfM[
'Ecms'])
190 dfMc = self.dfM.copy()
192 Plotter.plotLine(dfMc,
'Ecms',
'blue', label=f
'mumu method (+{round(d,1)} MeV)')
194 plotSplitLines(self.dfC)
195 setPlotRange(self.dfC, tag)
198 plt.ylabel(
'Ecms [MeV]')
202 loc =
'plots/allData'
203 if limits
is not None:
204 plt.xlim(datetime.strptime(limits[0],
'%Y-%m-%d'), datetime.strptime(limits[1],
'%Y-%m-%d'))
205 loc =
'plots/' + limits[0] +
'to' + limits[1]
207 os.makedirs(loc, exist_ok=
True)
208 plt.savefig(f
'{loc}/EcmsComparison{tag}.png')
211 def plotSpreadComparison(self, limits=None):
213 Plot Ecms spread estimated from the combined method
and B-decay method
216 Plotter.plotLine(self.dfB, 'spread',
'green', label=
'B decay method')
217 Plotter.plotLine(self.dfC,
'spread',
'red', label=
'Combined method')
219 plotSplitLines(self.dfC)
224 plt.ylabel(
'spread [MeV]')
226 loc =
'plots/allData'
227 if limits
is not None:
228 plt.xlim(datetime.strptime(limits[0],
'%Y-%m-%d'), datetime.strptime(limits[1],
'%Y-%m-%d'))
229 loc =
'plots/' + limits[0] +
'to' + limits[1]
231 os.makedirs(loc, exist_ok=
True)
232 plt.savefig(f
'{loc}/SpreadComparison.png')
236 def plotCurve(df, var, label, withBand=True, withCurve=True):
238 Plot curve with possible error band where large intervals are distinguished by color
241 def plotBands(df, var, color, withBand=True, withCurve=True):
243 nan = np.full(len(df), np.nan)
244 avg = (df[
't1']+df[
't2']).to_numpy()/2
245 times = np.c_[df[[
't1',
't2']].to_numpy(), avg].ravel()
246 eCMS = np.c_[df[[var, var]].to_numpy(), nan].ravel()
250 eCMSu = np.c_[df[[varUnc, varUnc]].to_numpy(), nan].ravel()
251 plt.fill_between(times, eCMS-eCMSu, eCMS+eCMSu, alpha=0.15, color=color)
253 plt.plot(times, eCMS, linewidth=2, color=color)
255 nan = np.full(len(df), np.nan)
256 avg = (df[
't1']+df[
't2']).to_numpy()/2
257 timesg = np.c_[df[
't1'].to_numpy(), avg, df[
't2'].to_numpy()].ravel()
258 eCMSg = np.c_[df[var].to_numpy(), nan, df[var].to_numpy()].ravel()
259 timesg = toJST(timesg)
262 plt.plot(timesg, eCMSg, linewidth=2, color=
'gray', alpha=0.35)
264 plotBands(df[df[
'type'] == 1], var,
'red', withBand, withCurve)
265 plotBands(df[df[
'type'] == -1], var,
'blue', withBand, withCurve)
267 def plotShift(self, limits=None):
269 Plot shift between Bhad and mumu method
272 Plotter.plotCurve(self.dfC, 'shift', label=
'Combined method')
274 plotSplitLines(self.dfC)
277 plt.ylabel(
'shift [MeV]')
279 loc =
'plots/allData'
280 if limits
is not None:
281 plt.xlim(datetime.strptime(limits[0],
'%Y-%m-%d'), datetime.strptime(limits[1],
'%Y-%m-%d'))
282 loc =
'plots/' + limits[0] +
'to' + limits[1]
284 os.makedirs(loc, exist_ok=
True)
285 plt.savefig(f
'{loc}/Shift.png')
289 def plotEcms(self, limits=None, tag=''):
291 Plot Ecms of the combined fit, 'tag' can be
'4S' means that y-axis
is zoomed to 4S mass region,
'Off' to off-resonance
294 Plotter.plotCurve(self.dfC, 'Ecms', label=
'Combined method')
296 plotSplitLines(self.dfC)
299 plt.ylabel(
'Ecms [MeV]')
301 setPlotRange(self.dfC, tag)
303 loc =
'plots/allData'
304 if limits
is not None:
305 plt.xlim(datetime.strptime(limits[0],
'%Y-%m-%d'), datetime.strptime(limits[1],
'%Y-%m-%d'))
306 loc =
'plots/' + limits[0] +
'to' + limits[1]
308 os.makedirs(loc, exist_ok=
True)
309 plt.savefig(f
'{loc}/Ecms{tag}.png')
312 def plotPulls(self, limits=None):
314 Plot pulls of the combined fit.
317 dfC = self.dfC.copy()
320 Plotter.plotCurve(dfC,
'pull', label=
'Combined method', withBand=
False)
321 Plotter.plotCurve(dfC,
'ref', label=
'Combined method', withCurve=
False)
323 plotSplitLines(self.dfC)
329 loc =
'plots/allData'
330 if limits
is not None:
331 plt.xlim(datetime.strptime(limits[0],
'%Y-%m-%d'), datetime.strptime(limits[1],
'%Y-%m-%d'))
332 loc =
'plots/' + limits[0] +
'to' + limits[1]
334 os.makedirs(loc, exist_ok=
True)
335 plt.savefig(f
'{loc}/Pull.png')
340def read_Combined_data(outputDir):
342 It reads the calibration table from the text file produced by the CAF calibration.
343 This text file includes the results
from the final calibration.
344 (combined
or mumu-only)
348 with open(outputDir +
'/finalEcmsCalib.txt')
as text_file:
349 for i, ll
in enumerate(text_file):
352 ll = ll.strip().split()
356 ll[1])), (int(ll[2]), int(ll[3]), int(ll[4]), int(ll[5])), int(
369def read_Bonly_data(outputDir):
371 It reads the calibration table from the text file produced by the CAF calibration.
372 This text file includes the results
from the B-only calibration.
376 with open(outputDir +
'/BonlyEcmsCalib.txt')
as text_file:
377 for i, ll
in enumerate(text_file):
380 ll = ll.strip().split()
384 ll[1])), (int(ll[2]), int(ll[3]), int(ll[4]), int(ll[5])),
385 (float(ll[2+4]), float(ll[3+4])), (float(ll[4+4]), float(ll[5+4]))))
389def read_mumu_data(outputDir):
391 It reads the calibration table from the text file produced by the CAF calibration.
392 This text file includes the results
from the mumu-based calibration.
396 with open(outputDir +
'/mumuEcalib.txt')
as text_file:
397 for i, ll
in enumerate(text_file):
400 ll = ll.strip().split()
404 ll[2+1])), (int(ll[2+2]), int(ll[2+3]), int(ll[2+4]), int(ll[2+5])),
405 (float(ll[2+6]), float(ll[2+7])), (int(ll[0]), int(ll[1]))))
410def create_hadB_fit_plots(outputDir, pdflatex):
412 Create multi-page pdf file with the fit plots
for hadronic B decays (
with mumu constraint).
413 The file
is created using pdflatex
416 arr = read_Combined_data(outputDir)
420 files = glob(outputDir+
'/'+dName +
'/*.pdf')
421 files = list(map(os.path.basename, files))
423 items = OrderedDict()
426 res = re.search(
'B[p0]_([0-9]*)_([0-9]*)\\.pdf', f)
427 t = int(res.group(1))
428 i = int(res.group(2))
432 items[t] = max(items[t], i + 1)
434 times = sorted(items)
436 header =
"""\\documentclass[aspectratio=169]{beamer}
437 \\usepackage{graphicx}
446 for i0Temp
in range(len(arr)):
447 if int(round(arr[i0Temp][0][0], 0)) == t:
450 assert(i0
is not None)
452 frac = 1. / (items[t] + 0.2)
454 frac = 1. / (items[t] + 0.3)
456 body +=
'\\begin{frame}[t]\n'
457 shift, shiftE = str(round(arr[i0][5][0], 1)), str(round(arr[i0][5][1], 1))
458 spread, spreadE = str(round(arr[i0][6][0], 1)), str(round(arr[i0][6][1], 1))
460 body +=
'$E_\\mathrm{shift} = (' + shift +
'\\pm' + shiftE +
')$~MeV \\hspace{2cm} \n'
461 body +=
'$E_\\mathrm{spread} = (' + spread +
'\\pm' + spreadE +
')$~MeV \\\\ \\vspace{0.5cm}\n'
462 for iShort
in range(items[t]):
464 tStart, tEnd = arr[i][0][0], arr[i][0][1]
466 exp1, run1, exp2, run2 = map(str, arr[i][1])
468 eCMS, eCMSe = str(round(arr[i][3][0], 1)), str(round(arr[i][3][1], 1))
469 pull = str(round(arr[i][4], 2))
471 t1 = datetime.utcfromtimestamp((tStart + 9) * 3600).strftime(
'%y-%m-%d %H:%M')
472 t2 = datetime.utcfromtimestamp((tEnd + 9) * 3600).strftime(
'%y-%m-%d %H:%M')
474 body +=
'\\begin{minipage}{' + str(frac) +
'\\textwidth}\n'
475 body +=
'\\begin{center}\n'
476 body +=
'\\scriptsize ' + exp1 +
' '+run1 +
' \\\\\n'
477 body +=
'\\scriptsize ' + exp2 +
' '+run2 +
' \\\\\n'
478 body +=
'\\scriptsize ' + t1 +
' \\\\\n'
479 body +=
'\\scriptsize ' + t2 +
' \\\\\n'
480 body +=
'\\scriptsize ' + str(round(tEnd - tStart, 1)) +
' hours \\\\ \\vspace{0.3cm}\n'
481 body +=
'$(' + eCMS +
'\\pm' + eCMSe +
')$~MeV \\\\\n'
482 body +=
'$\\mathrm{pull} = ' + pull +
'$ \\\\\n'
483 body +=
'\\includegraphics[width=1.0\\textwidth]{' + outputDir + \
484 '/' + dName +
'/B0_' + str(t) +
'_' + str(iShort) +
'.pdf}\\\\\n'
485 body +=
'\\includegraphics[width=1.0\\textwidth]{' + outputDir + \
486 '/' + dName +
'/Bp_' + str(t) +
'_' + str(iShort) +
'.pdf}\n'
487 body +=
'\\end{center}\n'
488 body +=
'\\end{minipage}\n'
490 body +=
'\\end{frame}\n\n'
492 tail =
'\n\\end{document}'
494 whole = header + body + tail
496 os.makedirs(
'tmp', exist_ok=
True)
498 with open(
"tmp/hadBfits.tex",
"w")
as text_file:
499 text_file.write(whole)
501 subprocess.call(f
'{pdflatex} tmp/hadBfits.tex', shell=
True)
503 os.makedirs(
'plots', exist_ok=
True)
504 copyfile(
'hadBfits.pdf',
'plots/hadBfits.pdf')
506 ext = [
'aux',
'log',
'nav',
'out',
'pdf',
'snm',
'toc']
508 if os.path.exists(f
'hadBfits.{e}'):
509 os.remove(f
'hadBfits.{e}')
515def create_hadBonly_fit_plots(outputDir, pdflatex):
517 Create multi-page pdf file with the fit plots
for hadronic B decays (without using mumu constraint).
518 The file
is created using pdflatex
521 arr = read_Bonly_data(outputDir)
523 dName = 'plotsHadBonly'
525 files = glob(outputDir+
'/'+dName +
'/*.pdf')
526 files = list(map(os.path.basename, files))
531 res = re.search(
'B[p0]Single_([0-9]*)\\.pdf', f)
532 t = int(res.group(1))
535 items = sorted(items)
537 header =
"""\\documentclass[aspectratio=169]{beamer}
538 \\usepackage{graphicx}
544 for i, t
in enumerate(items):
546 body +=
'\\begin{frame}[t]\n'
548 exp1, run1, exp2, run2 = map(str, arr[i][1])
549 eCMS, eCMSe = str(round(arr[i][2][0], 1)), str(round(arr[i][2][1], 1))
550 spread, spreadE = str(round(arr[i][3][0], 1)), str(round(arr[i][3][1], 1))
552 body +=
'$E_\\mathrm{cms} = (' + eCMS +
'\\pm' + eCMSe +
')$~MeV \\\\\n'
553 body +=
'$E_\\mathrm{spread} = (' + spread +
'\\pm' + spreadE +
')$~MeV \\\\ \\vspace{0.5cm}\n'
555 tStart, tEnd = arr[i][0][0], arr[i][0][1]
557 t1 = datetime.utcfromtimestamp((tStart + 9) * 3600).strftime(
'%y-%m-%d %H:%M')
558 t2 = datetime.utcfromtimestamp((tEnd + 9) * 3600).strftime(
'%y-%m-%d %H:%M')
560 body +=
'\\begin{minipage}{' + str(0.99) +
'\\textwidth}\n'
561 body +=
'\\begin{center}\n'
562 body +=
'\\scriptsize ' + exp1 +
' ' + run1 +
'\\hspace{1cm} ' + t1 +
' \\\\\n'
563 body +=
'\\scriptsize ' + exp2 +
' ' + run2 +
'\\hspace{1cm} ' + t2 +
' \\\\\n'
564 body +=
'\\scriptsize ' + str(round(tEnd - tStart, 1)) +
' hours \\\\ \\vspace{0.3cm}\n'
565 body +=
'\\includegraphics[width=0.48\\textwidth]{'+outputDir +
'/' + dName +
'/B0Single_' + str(t) +
'.pdf}\n'
566 body +=
'\\includegraphics[width=0.48\\textwidth]{'+outputDir +
'/' + dName +
'/BpSingle_' + str(t) +
'.pdf}\n'
567 body +=
'\\end{center}\n'
568 body +=
'\\end{minipage}\n'
570 body +=
'\\end{frame}\n\n'
572 tail =
'\n\\end{document}'
574 whole = header + body + tail
576 os.makedirs(
'tmp', exist_ok=
True)
578 with open(
"tmp/hadBonlyFits.tex",
"w")
as text_file:
579 text_file.write(whole)
581 subprocess.call(f
'{pdflatex} tmp/hadBonlyFits.tex', shell=
True)
583 os.makedirs(
'plots', exist_ok=
True)
584 copyfile(
'hadBonlyFits.pdf',
'plots/hadBonlyFits.pdf')
586 ext = [
'aux',
'log',
'nav',
'out',
'pdf',
'snm',
'toc']
588 if os.path.exists(f
'hadBonlyFits.{e}'):
589 os.remove(f
'hadBonlyFits.{e}')
591 if os.path.exists(
"tmp/hadBonlyFits.tex"):
592 os.remove(
"tmp/hadBonlyFits.tex")
596def create_mumu_fit_plots(outputDir, pdflatex):
598 Create multi-page pdf file with the fit plots
for mumu method.
599 The file
is created using pdflatex
602 arr = read_mumu_data(outputDir)
605 for i, a
in enumerate(arr):
606 if a[3][1] == a[3][0] - 1:
607 limits.append((i - a[3][1], a[3][0]))
611 files = glob(outputDir+
'/'+dName +
'/*.pdf')
612 files = list(map(os.path.basename, files))
617 res = re.search(
'mumu_([0-9]*)\\.pdf', f)
618 t = int(res.group(1))
621 items = sorted(items)
623 header =
"""\\documentclass[aspectratio=169]{beamer}
624 \\usepackage{graphicx}
648 body +=
'\\begin{frame}[t]\n'
650 for i
in range(k, k+n):
652 body +=
'\\begin{minipage}{' + str(frac) +
'\\textwidth}\n'
654 exp1, run1, exp2, run2 = map(str, arr[i][1])
655 eCMS, eCMSe = str(round(arr[i][2][0], 1)), str(round(arr[i][2][1], 1))
657 tStart, tEnd = arr[i][0][0], arr[i][0][1]
659 t1 = datetime.utcfromtimestamp((tStart + 9) * 3600).strftime(
'%y-%m-%d %H:%M')
660 t2 = datetime.utcfromtimestamp((tEnd + 9) * 3600).strftime(
'%y-%m-%d %H:%M')
662 body +=
'\\begin{center}\n'
663 body +=
'\\tiny $E_\\mathrm{cms} = (' + eCMS +
'\\pm' + eCMSe +
')$~MeV \\\\\n'
664 body +=
'\\tiny ' + exp1 +
' ' + run1 +
'\\hspace{0.3cm} ' + t1 +
' \\\\\n'
665 body +=
'\\tiny ' + exp2 +
' ' + run2 +
'\\hspace{0.3cm} ' + t2 +
' \\\\\n'
666 body +=
'\\tiny ' + str(round(tEnd - tStart, 1)) +
' hours \\vspace{0.3cm}\n'
667 body +=
'\\includegraphics[trim=0.3cm 0.0cm 1.3cm 0.7cm,clip=true,width=0.99\\textwidth]{' + \
668 outputDir +
'/' + dName +
'/mumu_' + str(items[i]) +
'.pdf}\n'
669 body +=
'\\end{center}\n'
670 body +=
'\\end{minipage}\n'
672 body +=
'\\end{frame}\n\n'
674 tail =
'\n\\end{document}'
676 whole = header + body + tail
678 os.makedirs(
'tmp', exist_ok=
True)
680 with open(
"tmp/mumuFits.tex",
"w")
as text_file:
681 text_file.write(whole)
683 subprocess.call(f
'{pdflatex} tmp/mumuFits.tex', shell=
True)
685 os.makedirs(
'plots', exist_ok=
True)
686 copyfile(
'mumuFits.pdf',
'plots/mumuFits.pdf')
688 ext = [
'aux',
'log',
'nav',
'out',
'pdf',
'snm',
'toc']
690 if os.path.exists(f
'mumuFits.{e}'):
691 os.remove(f
'mumuFits.{e}')
693 if os.path.exists(
"tmp/mumuFits.tex"):
694 os.remove(
"tmp/mumuFits.tex")
698def run_validation(job_path, input_data_path, requested_iov, expert_config):
700 Create validation plots related to the Ecms calibration
704 result = subprocess.run(
'locate pdflatex | grep "/pdflatex$"', stdout=subprocess.PIPE, shell=
True)
705 pdflatex = result.stdout.strip().decode(
'utf-8')
708 if expert_config !=
'':
709 expert_config = json.loads(expert_config)
712 if expert_config
is not None and 'plotsRanges' in expert_config:
713 allLimits += expert_config[
'plotsRanges']
715 plt.figure(figsize=(18, 9))
716 rcParams[
'axes.formatter.useoffset'] =
False
718 location = f
'{job_path}/eCMS/0/algorithm_output'
720 plotter = Plotter(location)
722 for limits
in allLimits:
723 plotter.plotEcmsComparison(limits, tag=
'4S')
724 plotter.plotEcmsComparison(limits, tag=
'Off')
725 plotter.plotEcmsComparison(limits)
726 plotter.plotSpreadComparison(limits)
728 plotter.plotEcms(limits, tag=
'4S')
729 plotter.plotEcms(limits, tag=
'Off')
730 plotter.plotEcms(limits, tag=
'')
732 plotter.plotShift(limits)
733 plotter.plotPulls(limits)
736 create_hadB_fit_plots(location, pdflatex)
737 create_hadBonly_fit_plots(location, pdflatex)
738 create_mumu_fit_plots(location, pdflatex)
741 copyfile(f
'{location}/BonlyEcmsCalib.txt',
'plots/BonlyEcmsCalib.txt')
742 copyfile(f
'{location}/finalEcmsCalib.txt',
'plots/finalEcmsCalib.txt')
743 copyfile(f
'{location}/mumuEcalib.txt',
'plots/mumuEcalib.txt')
746if __name__ ==
"__main__":