Commit 733530de authored by Zhang, Yuanpeng's avatar Zhang, Yuanpeng ⭐
Browse files

add in GENII upgrade test scripts

parent bc2f24d4

Scripts/GENII_Tests/cal_config.json

0 → 100644
+19 −0
Original line number Diff line number Diff line 
{

    "Diamond": "/SNS/NOWD/IPTS-36028/nexus/NOWD_60444.nxs.h5",

    "RunCycle": "2025-1_11A_CAL",

    "Instrument": "PG3",

    "Date": "2025-05-27",

    "SampleEnv": "PAC",

    "Notes": "",

    "OutputDir": "/SNS/PG3/shared/CALIBRATION/autoreduce",

    "PDCalibration": {

        "TofBinningFirst": "34500,-0.0008,51000",

        "TofBinningSecond": "34500,-0.0005,51000",

        "TofBinningThird": "34500,-0.0003,51000",

        "ConstrainPeakPositions": false,

        "MaxPeakWindow": 0.075,

        "MinimumPeakHeight": 3.0,

        "PeakWidthPercent": 0.008,

        "TZEROrange": "0,10"

    }

}

Scripts/GENII_Tests/pg3_cal_legacy_gen2_test.py

0 → 100644
+394 −0
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

import socket

from Calibration.tofpd import diagnostics

import mantid



# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+

#                        ||||| 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")

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.get("ManualMaskFile", None)

working_dir = "/SNS/PG3/shared/CALIBRATION"



instr_dict = {

    "POWGEN": "PG3",

    "PG3": "PG3"

}

# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+

#                        ^^^^^ Set-up session ^^^^^

#                        |||||                |||||

# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+



# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+

#                     ||||| Calib running session |||||

#                     vvvvv                       vvvvv

# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+

print("[Info] Loading calibrant data...")

condt1 = mtd.doesExist(wksp)

if condt1:

    condt2 = mtd[wksp].getAxis(0).getUnit().name() == 'Time-of-flight'

if not (condt1 and condt2):

    LoadEventNexus(Filename=cal_config["Diamond"], 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.

print("[Info] Calibrating...")

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')



print("[Info] Merging manual mask...")

# 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,

             RefWorkspace="PG3_group",

             OutputWorkspace="manual_mask")

    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:

    pass



# Save calibration file

print("[Info] Saving 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"],

    "auto_calib_to_check"

)

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)



# write out log information

mantid_version = mantid.__version__

out_log_dir_1 = os.path.join(

    out_dir_1,

    "logs"

)

if not os.path.exists(out_log_dir_1):

    os.makedirs(out_log_dir_1)



out_log_dir_2 = os.path.join(

    out_dir_2,

    "logs"

)

if not os.path.exists(out_log_dir_2):

    os.makedirs(out_log_dir_2)



if notes.strip() != "":

    log_file_name = f"{instr_name}_{sam_env}_{dia_name}_{date}-{notes}.log"

else:

    log_file_name = f"{instr_name}_{sam_env}_{dia_name}_{date}.log"



with open(os.path.join(out_log_dir_1, log_file_name), "w") as f:

    f.write(f"{mantid_version}\n")

with open(os.path.join(out_log_dir_2, log_file_name), "w") as f:

    f.write(f"{mantid_version}\n")

# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+

#                     ^^^^^ 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

# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+

print("[Info] Generating diagnostics plots for calibration...")

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)



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



print("\n======================")

print("Done with calibration.")

print("======================")



# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+

#                      |||||      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

# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+