Commit 2961863e authored by Zhang, Yuanpeng's avatar Zhang, Yuanpeng ⭐
Browse files

much faster calibration

parent 7f2edaad

cal_config.json

+1 −1
Original line number Diff line number Diff line 
{

    "Diamond": "/SNS/NOM/IPTS-28922/nexus/NOM_228263.nxs.h5",

    "Instrument": "NOM",

    "Date": "2025-09-25",

    "Date": "2025-09-27",

    "SampleEnv": "cryostat",

    "OutputDir": "/SNS/NOM/shared/autoreduce/calibration",

    "GenShadowMask": "shadow_mask_cryostat_228263.in",

utils/nom_cal.py

+193 −90
Original line number Diff line number Diff line
@@ -296,13 +296,10 @@ for wksp_index in range(num_det): 
# shot but many left-overs to process). Here we did not perform a systematic

# survey to see where the balancing point is. The chunk size of 50 is just what

# we found gives reasonable success rate and overall processing time.

chi_s_min = 5.

first_round_failed = list()

for i in range(int(38024 / 50) + 1):

ExtractSpectra(InputWorkspace="dia_calibrant",

               OutputWorkspace="dia_lr_tmp",

                   StartWorkspaceIndex=i * 50,

                   EndWorkspaceIndex=min((i + 1) * 50, 38024))

               StartWorkspaceIndex=0,

               EndWorkspaceIndex=38024)

PDCalibration(InputWorkspace="dia_lr_tmp",

              TofBinning="300, -0.01, 16667",

              PreviousCalibrationTable="_cal",
@@ -316,8 +313,8 @@ for i in range(int(38024 / 50) + 1): 


CreateDetectorTable(InputWorkspace="dia_lr_tmp",

                    DetectorTableWorkspace="dia_lr_tmp_table")



    # Populate the offset and collect those failing first-shot fitting.

chi_s_min = 5.

first_round_failed = list()

for wksp_index in range(mtd["dia_lr_tmp"].getNumberHistograms() * 2):

    if wksp_index % 2 == 0:

        wksp_index_tmp = wksp_index // 2
@@ -331,13 +328,11 @@ for i in range(int(38024 / 50) + 1): 
            else:

                first_round_failed.append(det_id_tmp)



# Here we are performing the second-shot fitting for those failing pixels in the

# first round. If succeeded, we collect the offset, otherwise if it is stilling

# failing, we mask out the pixel.

for det_id in first_round_failed:

    ExtractSingleSpectrum(InputWorkspace="dia_calibrant",

ExtractSpectra(

    InputWorkspace="dia_calibrant",

    OutputWorkspace="dia_lr_snd_trial",

                          WorkspaceIndex=det_id)

    DetectorList=first_round_failed

)

PDCalibration(InputWorkspace="dia_lr_snd_trial",

              TofBinning="300, -0.005, 16667",

              PreviousCalibrationTable="_cal",
@@ -348,22 +343,74 @@ for det_id in first_round_failed: 
              CalibrationParameters="DIFC",

              OutputCalibrationTable="pdcal_table_snd_tmp",

              DiagnosticWorkspaces="pdcal_diag_snd_tmp")

    if float(mtd["pdcal_diag_snd_tmp_fitparam"].row(0)["chi2"]) < chi_s_min:

        tof_fitted_pos = float(mtd["pdcal_diag_snd_tmp_fitparam"].row(0)["centre"])



second_round_failed = []

for row_i, det_id in enumerate(first_round_failed):

    success = -1

    for j in range(2):

        if float(mtd["pdcal_diag_snd_tmp_fitparam"].row(row_i * 2 + j)["chi2"]) < chi_s_min:

            success = j

            if success == 0:

                break

    if success != -1:

        tof_fitted_pos = float(mtd["pdcal_diag_snd_tmp_fitparam"].row(row_i * 2 + success)["centre"])

        tof_fitted_pos = tof_fitted_pos / dia220_c_norm[success] * dia220_c_norm[0]

        difc_eng_tmp = difc_eng_init[det_id]

        dia220_c_fitted = tof_fitted_pos / difc_eng_tmp

        offset[det_id] = 1. + (dia220_c_norm[0] - dia220_c_fitted) / dia220_c_fitted

    else:

        second_round_failed.append(det_id)

        mtd["mask"].setY(det_id, [1.])



dia220_c_norm = (dia_a / np.sqrt(11.), dia_a / np.sqrt(8.), dia_a / np.sqrt(3.))



ExtractSpectra(

    InputWorkspace="dia_calibrant",

    OutputWorkspace="dia_lr_thrd_trial",

    DetectorList=second_round_failed

)

PDCalibration(InputWorkspace="dia_lr_thrd_trial",

              TofBinning="300, -0.001, 16667",

              PreviousCalibrationTable="_cal",

              PeakFunction="BackToBackExponential",

              BackgroundType="Linear",

              PeakPositions=dia220_c_norm,

              PeakWindow="1.0, 1.2, 1.2, 1.35, 1.9, 2.1",

              CalibrationParameters="DIFC",

              OutputCalibrationTable="pdcal_table_thrd_tmp",

              DiagnosticWorkspaces="pdcal_diag_thrd_tmp")



third_round_failed = []

for row_i, det_id in enumerate(second_round_failed):

    success = -1

    for j in range(3):

        if float(mtd["pdcal_diag_thrd_tmp_fitparam"].row(row_i * 3 + j)["chi2"]) < chi_s_min:

            success = j

            if success == 1:

                break

    if success != -1:

        tof_fitted_pos = float(mtd["pdcal_diag_thrd_tmp_fitparam"].row(row_i * 3 + success)["centre"])

        tof_fitted_pos = tof_fitted_pos / dia220_c_norm[success] * dia220_c_norm[1]

        difc_eng_tmp = difc_eng_init[det_id]

        dia220_c_fitted = tof_fitted_pos / difc_eng_tmp

        offset[det_id] = 1. + (dia220_c_norm[1] - dia220_c_fitted) / dia220_c_fitted

    else:

        third_round_failed.append(det_id)

        mtd["mask"].setY(det_id, [1.])



dia220_c_norm = (dia_a / np.sqrt(8.), dia_a / np.sqrt(3.))



# Same logic as above, but specifically for bank #3-#5, we will be throwing all

# involved pixels to `PDCalibration` in a single shot and collect those failing

# ones and perform a second trial.

ExtractSpectra(InputWorkspace="dia_calibrant",

ExtractSpectra(

    InputWorkspace="dia_calibrant",

    OutputWorkspace="dia_hr_tmp",

    StartWorkspaceIndex=38025,

               EndWorkspaceIndex=82943)

PDCalibration(InputWorkspace="dia_hr_tmp",

    EndWorkspaceIndex=82943

)

PDCalibration(

    InputWorkspace="dia_hr_tmp",

    TofBinning="300, -0.0008, 16667",

    PreviousCalibrationTable="_cal",

    PeakFunction="BackToBackExponential",
@@ -372,12 +419,14 @@ PDCalibration(InputWorkspace="dia_hr_tmp", 
    PeakWindow="1.1, 1.45, 1.9, 2.1",

    CalibrationParameters="DIFC",

    OutputCalibrationTable="pdcal_table_tmp",

              DiagnosticWorkspaces="pdcal_diag_tmp")

    DiagnosticWorkspaces="pdcal_diag_tmp"

)



CreateDetectorTable(InputWorkspace="dia_hr_tmp",

                    DetectorTableWorkspace="dia_hr_tmp_table")



chi_s_min = 5.

first_round_failed = list()

first_round_failed = []

for wksp_index in range(mtd["dia_hr_tmp"].getNumberHistograms() * 2):

    if wksp_index % 2 == 0:

        wksp_index_tmp = wksp_index // 2
@@ -391,12 +440,14 @@ for wksp_index in range(mtd["dia_hr_tmp"].getNumberHistograms() * 2): 
            else:

                first_round_failed.append(det_id_tmp)



for det_id in first_round_failed:

    ExtractSingleSpectrum(InputWorkspace="dia_calibrant",

ExtractSpectra(

    InputWorkspace="dia_calibrant",

    OutputWorkspace="dia_hr_snd_trial",

                          WorkspaceIndex=det_id)

    PDCalibration(InputWorkspace="dia_hr_snd_trial",

                  TofBinning="300, -0.0008, 16667",

    DetectorList=first_round_failed

)

PDCalibration(

    InputWorkspace="dia_hr_snd_trial",

    TofBinning="6000, -0.0005, 16667",

    PreviousCalibrationTable="_cal",

    PeakFunction="BackToBackExponential",

    BackgroundType="Linear",
@@ -404,18 +455,70 @@ for det_id in first_round_failed: 
    PeakWindow="1.1, 1.45, 1.9, 2.1",

    CalibrationParameters="DIFC",

    OutputCalibrationTable="pdcal_table_snd_tmp",

                  DiagnosticWorkspaces="pdcal_diag_snd_tmp")

    if float(mtd["pdcal_diag_snd_tmp_fitparam"].row(0)["chi2"]) < chi_s_min:

        tof_fitted_pos = float(mtd["pdcal_diag_snd_tmp_fitparam"].row(0)["centre"])

    DiagnosticWorkspaces="pdcal_diag_snd_tmp"

)



second_round_failed = []

for row_i, det_id in enumerate(first_round_failed):

    success = -1

    for j in range(2):

        if float(mtd["pdcal_diag_snd_tmp_fitparam"].row(row_i * 2 + j)["chi2"]) < chi_s_min:

            success = j

            if success == 0:

                break

    if success != -1:

        tof_fitted_pos = float(mtd["pdcal_diag_snd_tmp_fitparam"].row(row_i * 2 + success)["centre"])

        tof_fitted_pos = tof_fitted_pos / dia220_c_norm[success] * dia220_c_norm[0]

        difc_eng_tmp = difc_eng_init[det_id]

        dia220_c_fitted = tof_fitted_pos / difc_eng_tmp

        offset[det_id] = 1. + (dia220_c_norm[0] - dia220_c_fitted) / dia220_c_fitted

    else:

        second_round_failed.append(det_id)

        mtd["mask"].setY(det_id, [1.])



dia220_c_norm = (dia_a / np.sqrt(11.), dia_a / np.sqrt(8.), dia_a / np.sqrt(3.))



ExtractSpectra(

    InputWorkspace="dia_calibrant",

    OutputWorkspace="dia_hr_thrd_trial",

    DetectorList=second_round_failed

)

PDCalibration(

    InputWorkspace="dia_hr_thrd_trial",

    TofBinning="3000, -0.001, 16667",

    PreviousCalibrationTable="_cal",

    PeakFunction="BackToBackExponential",

    BackgroundType="Linear",

    PeakPositions=dia220_c_norm,

    PeakWindow="1.0, 1.2, 1.2, 1.35, 1.9, 2.1",

    CalibrationParameters="DIFC",

    OutputCalibrationTable="pdcal_table_thrd_tmp",

    DiagnosticWorkspaces="pdcal_diag_thrd_tmp"

)



third_round_failed = []

for row_i, det_id in enumerate(second_round_failed):

    success = -1

    for j in range(3):

        if float(mtd["pdcal_diag_thrd_tmp_fitparam"].row(row_i * 3 + j)["chi2"]) < chi_s_min:

            success = j

            if success == 1:

                break

    if success != -1:

        tof_fitted_pos = float(mtd["pdcal_diag_thrd_tmp_fitparam"].row(row_i * 3 + success)["centre"])

        tof_fitted_pos = tof_fitted_pos / dia220_c_norm[success] * dia220_c_norm[1]

        difc_eng_tmp = difc_eng_init[det_id]

        dia220_c_fitted = tof_fitted_pos / difc_eng_tmp

        offset[det_id] = 1. + (dia220_c_norm[1] - dia220_c_fitted) / dia220_c_fitted

    else:

        third_round_failed.append(det_id)

        mtd["mask"].setY(det_id, [1.])



# Make sure we mask out those pixels involved in the new generated mask list.

MaskDetectors(Workspace="dia_calibrant",

              MaskedWorkspace="mask")

MaskDetectors(

    Workspace="dia_calibrant",

    MaskedWorkspace="mask"

)



# Write out the offset value. If a pixel is masked, we will assign 0 as the

# offset.

utils/nom_cal_new.py

deleted100755 → 0
+0 −1381

File deleted.

Preview size limit exceeded, changes collapsed.

utils/nom_cal_test.py

deleted100755 → 0
+0 −1288

File deleted.

Preview size limit exceeded, changes collapsed.

utils/test_nom_calib.py

deleted100644 → 0
+0 −143
Original line number Diff line number Diff line 
# import mantid algorithms, numpy and matplotlib

from mantid.simpleapi import *

import matplotlib.pyplot as plt

import numpy as np



dia_a = 3.5671299351

dia220_c_norm = (dia_a / np.sqrt(8.), dia_a / np.sqrt(3.))

num_monitor = 2

chi_s_min = 5

info = mtd["dia_calibrant"].detectorInfo()

offset = {}



ExtractSpectra(InputWorkspace="dia_calibrant",

               OutputWorkspace="dia_hr_tmp",

               StartWorkspaceIndex=38025,

               EndWorkspaceIndex=82943)

PDCalibration(InputWorkspace="dia_hr_tmp",

              TofBinning="300, -0.0008, 16667",

              PreviousCalibrationTable="_cal",

              PeakFunction="BackToBackExponential",

              BackgroundType="Linear",

              PeakPositions=dia220_c_norm,

              PeakWindow="1.1, 1.45, 1.9, 2.1",

              CalibrationParameters="DIFC",

              OutputCalibrationTable="pdcal_table_tmp",

              DiagnosticWorkspaces="pdcal_diag_tmp")



CreateDetectorTable(InputWorkspace="dia_hr_tmp",

                    DetectorTableWorkspace="dia_hr_tmp_table")

chi_s_min = 5.

first_round_failed = list()

for wksp_index in range(mtd["dia_hr_tmp"].getNumberHistograms() * 2):

    if wksp_index % 2 == 0:

        wksp_index_tmp = wksp_index // 2

        det_id_tmp = int(mtd["dia_hr_tmp_table"].row(wksp_index_tmp)["Detector ID(s)"])

        if not info.isMasked(det_id_tmp + num_monitor):

            if float(mtd["pdcal_diag_tmp_fitparam"].row(wksp_index)["chi2"]) < chi_s_min:

                tof_fitted_pos = float(mtd["pdcal_diag_tmp_fitparam"].row(wksp_index)["centre"])

                difc_eng_tmp = 5000

                dia220_c_fitted = tof_fitted_pos / difc_eng_tmp

                offset[det_id_tmp] = 1. + (dia220_c_norm[0] - dia220_c_fitted) / dia220_c_fitted

            else:

                first_round_failed.append(det_id_tmp)



print("Done", len(first_round_failed))



ExtractSpectra(

    InputWorkspace="dia_calibrant",

    OutputWorkspace="dia_hr_snd_trial",

    DetectorList=first_round_failed

)

PDCalibration(InputWorkspace="dia_hr_snd_trial",

              TofBinning="6000, -0.0005, 16667",

              PreviousCalibrationTable="_cal",

              PeakFunction="BackToBackExponential",

              BackgroundType="Linear",

              PeakPositions=dia220_c_norm,

              PeakWindow="1.1, 1.45, 1.9, 2.1",

              CalibrationParameters="DIFC",

              OutputCalibrationTable="pdcal_table_snd_tmp",

              DiagnosticWorkspaces="pdcal_diag_snd_tmp")



second_round_failed = []

for row_i, det_id in enumerate(first_round_failed):

    success = 0

    for j in range(2):

        if float(mtd["pdcal_diag_snd_tmp_fitparam"].row(row_i * 2 + j)["chi2"]) < chi_s_min:

            success += 1

    if success > 0:

        pass

    else:

        second_round_failed.append(det_id)

        mtd["mask"].setY(det_id, [1.])



print("Done", len(second_round_failed))



dia220_c_norm = (dia_a / np.sqrt(11.), dia_a / np.sqrt(8.), dia_a / np.sqrt(3.))



ExtractSpectra(

    InputWorkspace="dia_calibrant",

    OutputWorkspace="dia_hr_thrd_trial",

    DetectorList=second_round_failed

)

PDCalibration(InputWorkspace="dia_hr_thrd_trial",

              TofBinning="3000, -0.001, 16667",

              PreviousCalibrationTable="_cal",

              PeakFunction="BackToBackExponential",

              BackgroundType="Linear",

              PeakPositions=dia220_c_norm,

              PeakWindow="1.0, 1.2, 1.2, 1.35, 1.9, 2.1",

              CalibrationParameters="DIFC",

              OutputCalibrationTable="pdcal_table_thrd_tmp",

              DiagnosticWorkspaces="pdcal_diag_thrd_tmp")



third_round_failed = []

for row_i, det_id in enumerate(second_round_failed):

    success = 0

    for j in range(2):

        if float(mtd["pdcal_diag_snd_tmp_fitparam"].row(row_i * 2 + j)["chi2"]) < chi_s_min:

            success += 1

    if success > 0:

        pass

    else:

        third_round_failed.append(det_id)

        if len(third_round_failed) <= 10:

            print(row_i)

        mtd["mask"].setY(det_id, [1.])



print("Done", len(third_round_failed))

# 

# dia220_c_norm = (dia_a / np.sqrt(11.), dia_a / np.sqrt(8.), dia_a / np.sqrt(3.))

# 

# ExtractSpectra(

#     InputWorkspace="dia_calibrant",

#     OutputWorkspace="dia_hr_fth_trial",

#     DetectorList=third_round_failed

# )

# PDCalibration(InputWorkspace="dia_hr_fth_trial",

#               TofBinning="3000, -0.001, 16667",

#               PreviousCalibrationTable="_cal",

#               PeakFunction="BackToBackExponential",

#               BackgroundType="Linear",

#               PeakPositions=dia220_c_norm,

#               PeakWindow="1.0, 1.2, 1.2, 1.35, 1.9, 2.1",

#               CalibrationParameters="DIFC",

#               OutputCalibrationTable="pdcal_table_fth_tmp",

#               DiagnosticWorkspaces="pdcal_diag_fth_tmp")

# 

# fourth_round_failed = []

# for row_i, det_id in enumerate(third_round_failed):

#     success = 0

#     for j in range(3):

#         if float(mtd["pdcal_diag_fth_tmp_fitparam"].row(row_i * 3 + j)["chi2"]) < chi_s_min:

#             success += 1

#     if success > 0:

#         pass

#     else:

#         fourth_round_failed.append(det_id)

#         if len(fourth_round_failed) <= 10:

#             print(row_i)

#         mtd["mask"].setY(det_id, [1.])

# 

# print("Done", len(fourth_round_failed))