updating scripts (e9ee8d22) · Commits · sns-hfir-scse / diffraction / pg3_calibration

README.md

+49 −0

Original line number	Diff line number	Diff line
		@@ -12,6 +12,23 @@ How-To
		To run the calibration routine locally, one first needs to fill in necessary
		configurations in the `/SNS/PG3/shared/CALIBRATION/cal_config.json` file.

		> Some of the entries in the configuration file will be used as the part of
		generated calibration file name -- refer to the comments in section-3 for
		detailed information about the file naming.

		> Manual mask file can be specified to be merged with the automatically
		generated masks during the calibration process. To obtain the manual mask
		is not part of the current routine and one is expected to load in the
		diamond data and load in the calibration file initially generated without
		the manual mask, then inspect the diamond data to identify bad pixels, and
		then save the manual mask to an XML file.
		N.B. The specified file path to the manual mask should be relative to the
		`/SNS/PG3/shared/CALIBRATION` directory.

		> The manual mask can be specified as either a single string or a list of
		strings. For the latter case, all those manual mask files included in the
		list will be merged into the final mask.

		2. Run

		Open a terminal and change directory to `/SNS/PG3/shared/CALIBRATION`, and then
		@@ -25,6 +42,38 @@ How-To
		autoreduction. Meanwhile, diagnostics plots will be generated and saved into
		locations parallel to each of the two saved calibration files.

		> Information of the diamond run will be collected from the central
		configuration file, which will be used to generate the name of the
		calibration file. For example, if a config file looks like the following
		example, the calibration file will be named as
		`PG3_MAG-He3_d55064_2023-02-14_ZYP.h5`.

		```json
		{
		"Diamond": "/SNS/PG3/IPTS-2767/nexus/PG3_55064.nxs.h5",
		"RunCycle": "2023-1_11A_CAL",
		"Instrument": "PG3",
		"Date": "2023-02-14",
		"SampleEnv": "MAG-He3",
		"Notes": "ZYP",
		"OutputDir": "/SNS/PG3/shared/CALIBRATION/autoreduce",
		"PDCalibration": {
		"TofBinningFirst": "300,-0.0008,16667",
		"TofBinningSecond": "300,-0.0005,16667",
		"TofBinningThird": "300,-0.0003,16667",
		"ConstrainPeakPositions": false,
		"MaxPeakWindow": 0.075,
		"MinimumPeakHeight": 3.0,
		"PeakWidthPercent": 0.008,
		"TZEROrange": "0,10"
		},
		"ManualMaskFile": "2018_2_11A_CAL/mask/PG3_40481_manual_mask.xml"
		}
		```

		> TODO-1: Diagnostics for multiple banks => monitor & local
		> TODO-2: Lining up diagnostics for sample run

		---

		Yuanpeng Zhang @ 2023-10-19 Thu 12:06 PM EDT

pg3_cal.py

+37 −10

Original line number	Diff line number	Diff line
		@@ -23,13 +23,18 @@ from mantid.simpleapi import ApplyDiffCal, \
		import matplotlib.pyplot as plt
		import numpy as np
		import os
		import socket
		from Calibration.tofpd import diagnostics

		# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
		# \|\|\|\|\| Set-up session \|\|\|\|\|
		# vvvvv vvvvv
		# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
		host_name = socket.gethostname()
		if "autoreducer" in host_name:
		post_image = True
		else:
		post_image = False

		cal_w_dir = "/SNS/PG3/shared/CALIBRATION"
		cal_config_f = os.path.join(cal_w_dir, "cal_config.json")
		@@ -156,6 +161,7 @@ PDCalibration(InputWorkspace=wksp,
		DiagnosticWorkspaces='diag')

		# Merge in the manual mask
		if isinstance(manual_mask_f, str):
		manual_mask = os.path.join(working_dir, manual_mask_f)
		LoadMask(Instrument="POWGEN",
		InputFile=manual_mask,
		@@ -165,6 +171,27 @@ BinaryOperateMasks(InputWorkspace1="PDCalib_mask",
		InputWorkspace2="manual_mask",
		OperationType="OR",
		OutputWorkspace="Combined_Mask")
		elif isinstance(manual_mask_f, list):
		for count, mask_f in enumerate(manual_mask_f):
		manual_mask = os.path.join(working_dir, mask_f)
		LoadMask(Instrument="POWGEN",
		InputFile=manual_mask,
		RefWorkspace="PG3_group",
		OutputWorkspace="manual_mask")
		if count == 0:
		BinaryOperateMasks(InputWorkspace1="PDCalib_mask",
		InputWorkspace2="manual_mask",
		OperationType="OR",
		OutputWorkspace="Combined_Mask")
		else:
		BinaryOperateMasks(InputWorkspace1="Combined_Mask",
		InputWorkspace2="manual_mask",
		OperationType="OR",
		OutputWorkspace="Combined_Mask")
		else:
		err_msg = "The value for 'ManualMaskFile' should be either a single "
		err_msg += "path or a list of paths."
		raise TypeError(err_msg)

		# Save calibration file
		instr_name = instr_dict[cal_config["Instrument"]]

pg3_cal_local.py

deleted100644 → 0

+0 −327

Original line number	Diff line number	Diff line
		import base64
		from finddata.publish_plot import publish_plot
		from io import BytesIO
		import json
		from mantid.simpleapi import ApplyDiffCal, \
		BinaryOperateMasks, \
		ConvertToPointData, \
		ConvertUnits, \
		CreateDetectorTable, \
		CreateEmptyTableWorkspace, \
		CreateGroupingWorkspace, \
		DiffractionFocussing, \
		Load, \
		LoadDiffCal, \
		LoadEventNexus, \
		LoadMask, \
		MaskDetectors, \
		mtd, \
		PDCalibration, \
		Rebin, \
		SaveDiffCal, \
		RenameWorkspace
		import matplotlib.pyplot as plt
		import numpy as np
		import os
		from Calibration.tofpd import diagnostics

		# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
		# \|\|\|\|\| Set-up session \|\|\|\|\|
		# vvvvv vvvvv
		# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
		post_image = False

		cal_w_dir = "/SNS/PG3/shared/CALIBRATION"
		cal_config_f = os.path.join(cal_w_dir, "cal_config.json")
		with open(cal_config_f, "r") as f:
		cal_config = json.load(f)

		dia_file_name = os.path.basename(cal_config["Diamond"])
		if "." in dia_file_name:
		dia_run_num = dia_file_name.split(".")[0].split("_")[1]
		else:
		dia_run_num = dia_file_name.split("_")[1]
		wksp = dia_file_name.split(".")[0]

		# Diamond nominal peak positions and this SHOULD NOT be changed
		# unless we change to other type of calibrant than diamond.
		peakpositions = np.asarray((
		2.06,
		1.2615,
		1.0758,
		0.892,
		0.8186,
		0.7283,
		0.6867,
		0.6307,
		0.5642,
		0.5441,
		0.515,
		0.4996,
		0.4768,
		0.4645,
		0.4205,
		0.3916,
		0.3499,
		0.3257,
		0.3117))

		max_peak_window = cal_config["PDCalibration"]["MaxPeakWindow"]
		min_peak_height = cal_config["PDCalibration"]["MinimumPeakHeight"]
		peak_width_perct = cal_config["PDCalibration"]["PeakWidthPercent"]
		tzero_range = cal_config["PDCalibration"]["TZEROrange"]
		tof_bin_first = cal_config["PDCalibration"]["TofBinningFirst"]
		tof_bin_second = cal_config["PDCalibration"]["TofBinningSecond"]
		tof_bin_third = cal_config["PDCalibration"]["TofBinningThird"]
		constr_pp = cal_config["PDCalibration"]["ConstrainPeakPositions"]

		run_cycle = cal_config["RunCycle"]
		manual_mask_f = cal_config["ManualMaskFile"]
		working_dir = "/SNS/PG3/shared/CALIBRATION"

		instr_dict = {
		"POWGEN": "PG3",
		"PG3": "PG3"
		}
		# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
		# ^^^^^ Set-up session ^^^^^
		# \|\|\|\|\| \|\|\|\|\|
		# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+

		# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
		# \|\|\|\|\| Calib running session \|\|\|\|\|
		# vvvvv vvvvv
		# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
		condt1 = mtd.doesExist(wksp)
		if condt1:
		condt2 = mtd[wksp].getAxis(0).getUnit().name() == 'Time-of-flight'
		if not (condt1 and condt2):
		LoadEventNexus(Filename=wksp, OutputWorkspace=wksp)

		# Align and focus the diamond pattern across the whole instrument using the
		# engineering calibration values, given the surveyed detector locations.
		CreateGroupingWorkspace(InputWorkspace=wksp,
		OutputWorkspace='PG3_group',
		GroupDetectorsBy='All')
		ConvertUnits(InputWorkspace=wksp,
		OutputWorkspace=wksp + '_geom',
		Target='dSpacing',
		EMode='Elastic')
		DiffractionFocussing(InputWorkspace=wksp + '_geom',
		OutputWorkspace=wksp + '_geom',
		GroupingWorkspace='PG3_group')
		Rebin(InputWorkspace=wksp + '_geom',
		OutputWorkspace=wksp + '_geom',
		Params=(0.2, -.0004, 2.6),
		PreserveEvents=False)

		# This is the kernel part of the calibration, and fundamentally we are running
		# the `PDCalibration` algorithm for multiple times to get a reasonable
		# calibration outcome.
		PDCalibration(InputWorkspace=wksp,
		TofBinning=tof_bin_first,
		PeakPositions=peakpositions,
		ConstrainPeakPositions=constr_pp,
		PeakWindow=max_peak_window,
		MinimumPeakHeight=min_peak_height,
		PeakWidthPercent=peak_width_perct,
		TZEROrange=tzero_range,
		OutputCalibrationTable='PDCalib',
		DiagnosticWorkspaces='diag')
		RenameWorkspace(InputWorkspace='PDCalib_mask',
		OutputWorkspace='PG3_mask')

		peak_width_perct_sec = peak_width_perct / 2.
		PDCalibration(InputWorkspace=wksp,
		TofBinning=tof_bin_second,
		PreviousCalibrationTable='PDCalib',
		PeakPositions=peakpositions,
		ConstrainPeakPositions=constr_pp,
		PeakWindow=max_peak_window,
		MinimumPeakHeight=min_peak_height,
		PeakWidthPercent=peak_width_perct_sec,
		TZEROrange=tzero_range,
		OutputCalibrationTable='PDCalib',
		DiagnosticWorkspaces='diag')
		PDCalibration(InputWorkspace=wksp,
		TofBinning=tof_bin_third,
		PreviousCalibrationTable='PDCalib',
		PeakPositions=peakpositions,
		ConstrainPeakPositions=constr_pp,
		PeakWindow=max_peak_window,
		MinimumPeakHeight=min_peak_height,
		PeakWidthPercent=peak_width_perct_sec,
		TZEROrange=tzero_range,
		OutputCalibrationTable='PDCalib',
		DiagnosticWorkspaces='diag')

		# Merge in the manual mask
		manual_mask = os.path.join(working_dir, manual_mask_f)
		LoadMask(Instrument="POWGEN",
		InputFile=manual_mask,
		RefWorkspace="PG3_group",
		OutputWorkspace="manual_mask")
		BinaryOperateMasks(InputWorkspace1="PDCalib_mask",
		InputWorkspace2="manual_mask",
		OperationType="OR",
		OutputWorkspace="Combined_Mask")

		# Save calibration file
		instr_name = instr_dict[cal_config["Instrument"]]
		sam_env = cal_config["SampleEnv"]
		dia_name = f"d{dia_run_num}"
		date = cal_config["Date"]
		notes = cal_config["Notes"]
		if notes.strip() != "":
		cal_file_name = f"{instr_name}_{sam_env}_{dia_name}_{date}-{notes}.h5"
		else:
		cal_file_name = f"{instr_name}_{sam_env}_{dia_name}_{date}.h5"
		out_dir_1 = os.path.join(working_dir, cal_config["RunCycle"])
		out_dir_2 = cal_config["OutputDir"]
		if working_dir not in out_dir_2:
		out_dir_2 = os.path.join(working_dir, out_dir_2)
		if not os.path.exists(out_dir_1):
		os.makedirs(out_dir_1)
		if not os.path.exists(out_dir_2):
		os.makedirs(out_dir_2)
		cal_file_1 = os.path.join(out_dir_1, cal_file_name)
		cal_file_2 = os.path.join(out_dir_2, cal_file_name)

		SaveDiffCal(CalibrationWorkspace="PDCalib",
		GroupingWorkspace="PG3_group",
		MaskWorkspace="Combined_Mask",
		Filename=cal_file_1)
		SaveDiffCal(CalibrationWorkspace="PDCalib",
		GroupingWorkspace="PG3_group",
		MaskWorkspace="Combined_Mask",
		Filename=cal_file_2)
		# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
		# ^^^^^ Calib running session ^^^^^
		# \|\|\|\|\| \|\|\|\|\|
		# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+

		# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
		# Done with calibration
		# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+

		# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
		# Publish diagnostics plots \\|/
		# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
		#
		# Then we load in the saved calibration file, generate the diagnostics plots
		# and post them onto monitor.
		#
		# Here we first load in the raw diamond data and extract the engineering
		# calibration constants.
		#
		dia_check = Load(dia_file_name)
		CreateDetectorTable(InputWorkspace=dia_check,
		DetectorTableWorkspace="_det_table")
		num_det = mtd["dia_check"].getNumberHistograms()

		out_table = CreateEmptyTableWorkspace(OutputWorkspace="_cal_eng")
		out_table.addColumn("int", "detid")
		out_table.addColumn("double", "difc")
		out_table.addColumn("double", "difa")
		out_table.addColumn("double", "tzero")

		for wksp_index in range(num_det):
		det_id = int(mtd["_det_table"].row(wksp_index)["Detector ID(s)"])
		difc_eng = mtd["_det_table"].row(wksp_index)["DIFC"]
		difa = 0.
		tzero = 0.
		new_row = {
		"detid": det_id,
		"difc": difc_eng,
		"difa": difa,
		"tzero": tzero,
		}
		out_table.addRow(new_row)
		#
		# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
		# \|\|\|\|\| Diagnostics-1 \|\|\|\|\|
		# vvvvv Pixel-by-pixel line-up vvvvv
		# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
		Rebin(InputWorkspace=dia_check,
		OutputWorkspace=dia_check,
		Params="300,-0.0008,16667")
		LoadDiffCal(InstrumentName="PG3",
		Filename=cal_file_2,
		TofMin=300)
		ApplyDiffCal(InstrumentWorkspace=dia_check,
		CalibrationWorkspace="_cal")
		MaskDetectors(Workspace=dia_check,
		MaskedWorkspace="_mask")
		ConvertUnits(InputWorkspace=dia_check,
		OutputWorkspace="dia_check_d",
		Target="dSpacing",
		EMode="Elastic")
		Rebin(InputWorkspace="dia_check_d",
		OutputWorkspace="dia_check_d",
		Params="0.0,0.01,3.0")
		ConvertToPointData(InputWorkspace="dia_check_d",
		OutputWorkspace="dia_check_d")
		data = mtd["dia_check_d"].extractY()
		for det, mask in enumerate(mtd['_mask'].extractY()):
		if mask[0] == 1:
		data[det, :] = None

		max_all = np.nanmax(data)
		zmax_val = max_all / 10.

		fig, ax = plt.subplots()
		heatmap = ax.pcolormesh(data, cmap='viridis', vmax=zmax_val)
		colorbar = plt.colorbar(heatmap, ax=ax)
		ax.set_xticks(np.arange(0, 350, 50), np.arange(0., 3.5, 0.5))
		ax.set_xlabel('d-spacing')
		ax.set_ylabel('Detectors')

		diag_file_name = cal_file_name.replace(".h5", "_diag_1.png")
		diag_dir_1 = os.path.join(out_dir_1, "diagnostics")
		if not os.path.exists(diag_dir_1):
		os.makedirs(diag_dir_1)
		diag_dir_2 = os.path.join(out_dir_2, "diagnostics")
		if not os.path.exists(diag_dir_2):
		os.makedirs(diag_dir_2)
		diag_file_1 = os.path.join(diag_dir_1, diag_file_name)
		fig.savefig(diag_file_1, format='png', dpi=300)
		diag_file_2 = os.path.join(diag_dir_2, diag_file_name)
		fig.savefig(diag_file_2, format='png', dpi=300)

		if post_image:
		tmpfile = BytesIO()
		fig.savefig(tmpfile, format='png', dpi=200)
		encoded = base64.b64encode(tmpfile.getvalue()).decode('utf-8')

		html_tmp = '<h4>Pixel-by-pixel line-up for diamond</h4>'
		html_tmp += '<img src=\'data:image/png;base64,{}\'>'.format(encoded)
		html = "<div>{}</div>".format(html_tmp)

		# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
		# \|\|\|\|\| Diagnostics-2 \|\|\|\|\|
		# vvvvv Calib versus Engineer vvvvv
		# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+

		fig, ax = diagnostics.difc_plot2d("_cal", "_cal_eng", instr_ws="_group",
		mask="_mask", vrange=(0, 5))

		diag_file_name = cal_file_name.replace(".h5", "_diag_2.png")
		diag_file_1 = os.path.join(diag_dir_1, diag_file_name)
		fig.savefig(diag_file_1, format='png', dpi=300)
		diag_file_2 = os.path.join(diag_dir_2, diag_file_name)
		fig.savefig(diag_file_2, format='png', dpi=300)

		if post_image:
		tmpfile = BytesIO()
		fig.savefig(tmpfile, format='png', dpi=200)
		encoded = base64.b64encode(tmpfile.getvalue()).decode('utf-8')

		html_tmp = '<h4>Calibrated vs. Engineering calibration constant</h4>'
		html_tmp += '<img src=\'data:image/png;base64,{}\'>'.format(encoded)
		html += "<div>{}</div>".format(html_tmp)

		request = publish_plot(instr_name, dia_run_num, files={"file": html})

		# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
		# Done with calibration diagnostics plots publishing
		# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+