Commit 9d4fb919 authored by Shang, Yingrui's avatar Shang, Yingrui Committed by Zhou, Wenduo
Browse files

change the calculation of sigma q for equation 11.26-11.28 in master doc

parent 84be8bdd

drtsans/iq.py

+10 −10
Original line number Diff line number Diff line
@@ -877,18 +877,18 @@ def _do_1d_weighted_binning( 


    General description of algorithm:



    Equation 11.22

    Equation 11.26

    I(Q') = sum_{Q, lambda}^{K} (I(Q, lambda) / sigma(Q, lambda)^2) /

            sum_{Q, lambda}^{K} (1 / sigma(Q, lambda)^2)



    Equation 11.23

    Equation 11.27

    sigmaI(Q') = sqrt(sum_{Q, lambda}^{K} (sigma(Q, lambda / sigma(Q, lambda)^2)^2) /

                 sum_{Q, lambda}^{K} (1 / sigma(Q, lambda)^2)

               = sqrt(sum_{Q, lambda}^{K} (1 / sigma(Q, lambda)^2)) /

                 sum_{Q, lambda}^{K}(1/sigma(Q, lambda)^2)

               = 1 / sqrt(sum_{Q, lambda}^{K} (1 / sigma(Q, lambda)^2))



    Equation 11.24

    Equation 11.28

    sigmaQ(Q') = sum_{Q, lambda}^{K}(sigmaQ(Q, lambda)/sigma^2(Q, lambda)^2) /

                 sum_{Q, lambda}^{K}(1/sigma(Q, lambda)^2)
@@ -926,17 +926,17 @@ def _do_1d_weighted_binning( 
            mod_q_array, bins=bin_edges, weights=invert_sigma2_array

        )



        # Calculate Equation 11.22: I(Q)

        # Calculate Equation 11.26: I(Q)

        #  I(Q') = sum_{Q, lambda}^{K} (I(Q, lambda) / sigma(Q, lambda)^2) /

        #              sum_{Q, lambda}^{K} (1 / sigma(Q, lambda)^2)

        # denominator in Equation 11.22: sum_{Q, lambda}^{K} (I(Q, lambda) / sigma(Q, lambda)^2)

        i_raw_array, _ = np.histogram(

            mod_q_array, bins=bin_edges, weights=intensity_array * invert_sigma2_array

        )

        # denominator divided by nominator (11.22)

        # numerator divided by denominator (11.26)

        binned_intensity_array = i_raw_array / w_array



        # Calculate equation 11.23: sigmaI(Q)

        # Calculate equation 11.27: sigmaI(Q)

        # sigmaI(Q') = sqrt(sum_{Q, lambda}^{K} (sigma(Q, lambda / sigma(Q, lambda)^2)^2) /

        #              sum_{Q, lambda}^{K} (1 / sigma(Q, lambda)^2)

        #            = sqrt(sum_{Q, lambda}^{K} (1 / sigma(Q, lambda)^2)) /
@@ -945,18 +945,18 @@ def _do_1d_weighted_binning( 
        # Thus histogrammed sigmaI can be obtained from histogrammed invert_sigma2_array directly

        binned_error_array = 1 / np.sqrt(w_array)



        # Calculate equation 11.24:  sigmaQ (i.e., Q resolution)

        # Calculate equation 11.28:  sigmaQ (i.e., Q resolution)

        # sigmaQ(Q') = sum_{Q, lambda}^{K}(sigmaQ(Q, lambda)/sigma^2(Q, lambda)^2) /

        #              sum_{Q, lambda}^{K}(1/sigma(Q, lambda)^2)

        # denominator in Equation 11.24: sum_{Q, lambda}^{K}(sigmaQ(Q, lambda)/sigma^2(Q, lambda)^2)

        # denominator in Equation 11.28: sum_{Q, lambda}^{K}(1/sigma(Q, lambda)^2)

        if delta_q_array is None:

            binned_dq_array = None

        else:

            binned_dq, _ = np.histogram(

                mod_q_array, bins=bin_edges, weights=delta_q_array * invert_sigma2_array

            )

            # denominator divided by nominator (11.24)

            binned_dq_array = binned_dq / i_raw_array

            # numerator divided by denominator (11.28)

            binned_dq_array = binned_dq / w_array



        return binned_intensity_array, binned_error_array, binned_dq_array

tests/integration/new/drtsans/tof/eqsans/test_inelastic_incoherence_correction.py

+5 −1
Original line number Diff line number Diff line
@@ -346,11 +346,15 @@ def test_incoherence_correction_elastic_normalization(reference_dir, generatecle 
    # Check output result

    iq1d_base_name = "EQSANS_125707__Iq.dat"

    test_iq1d_file = os.path.join(test_dir, iq1d_base_name)

    # FIXME: The gold data are not stored inside the repository so when

    # gold data are changed a version prefix is added with date and developer

    # information. The old data will be kept as it is.

    version = "20220321_rys_"

    assert os.path.exists(

        test_iq1d_file

    ), f"Expected test result {test_iq1d_file} does not exist"

    gold_iq1d_file = os.path.join(

        reference_dir.new.eqsans, "test_incoherence_correction", iq1d_base_name

        reference_dir.new.eqsans, "test_incoherence_correction", version + iq1d_base_name

    )

    assert os.path.exists(

        gold_iq1d_file

tests/integration/new/drtsans/tof/eqsans/test_reduction_api.py

+7 −1
Original line number Diff line number Diff line
@@ -220,9 +220,15 @@ def test_weighted_binning_setup(run_config, basename, generatecleanfile, referen 
    # Verify binned I(Q)

    gold_file_dict = dict()

    gold_dir = os.path.join(reference_dir.new.eqsans, "gold_data")



    # FIXME: The gold data are not stored inside the repository so when

    # gold data are changed a version prefix is added with date and developer

    # information. The old data will be kept as it is.

    version = "20220321_rys_"



    for frame_index in range(2):

        iq1d_h5_name = os.path.join(

            gold_dir, f"gold_88980_weighted_1d_{frame_index}.h5"

            gold_dir, f"{version}gold_88980_weighted_1d_{frame_index}.h5"

        )

        gold_file_dict[1, frame_index, 0] = iq1d_h5_name

        iq2d_h5_name = os.path.join(

tests/unit/new/drtsans/test_i_of_q_1d_binning.py

+30 −26
Original line number Diff line number Diff line
@@ -286,26 +286,26 @@ def test_1d_weighted_binning(): 
    # Generate test data

    test_qiew_array = np.array(

        [

            [0.013, np.nan, np.nan, 1.0],

            [0.021, np.nan, np.nan, 1.0],

            [0.032, np.nan, np.nan, 1.0],

            [0.041, 10.000, 1.1000, 1.0],

            [0.048, 10.000, 1.5000, 1.0],

            [0.060, 10.000, 1.1300, 1.0],

            [0.070, 10.000, 2.3000, 1.0],

            [0.080, 10.000, 0.7000, 1.0],

            [0.090, 10.000, 0.8000, 1.0],

            [0.099, 10.000, 1.2000, 1.0],

            [0.011, np.nan, np.nan, 2.0],

            [0.021, 12.000, 2.3000, 2.0],

            [0.031, 12.000, 1.5000, 2.0],

            [0.039, 12.000, 1.7000, 2.0],

            [0.049, 12.000, 0.5000, 2.0],

            [0.062, 12.000, 2.4000, 2.0],

            [0.073, 12.000, 0.6000, 2.0],

            [0.079, 12.000, 0.9000, 2.0],

            [0.091, 12.000, 1.3000, 2.0],

            [0.099, 12.000, 1.4000, 2.0],

            [0.013, np.nan, np.nan, 1.0, 0.0005],

            [0.021, np.nan, np.nan, 1.0, 0.0007],

            [0.032, np.nan, np.nan, 1.0, 0.0010],

            [0.041, 10.000, 1.1000, 1.0, 0.0002],

            [0.048, 10.000, 1.5000, 1.0, 0.0001],

            [0.060, 10.000, 1.1300, 1.0, 0.0004],

            [0.070, 10.000, 2.3000, 1.0, 0.0003],

            [0.080, 10.000, 0.7000, 1.0, 0.0006],

            [0.090, 10.000, 0.8000, 1.0, 0.0009],

            [0.099, 10.000, 1.2000, 1.0, 0.00035],

            [0.011, np.nan, np.nan, 2.0, 0.00015],

            [0.021, 12.000, 2.3000, 2.0, 0.0002],

            [0.031, 12.000, 1.5000, 2.0, 0.0006],

            [0.039, 12.000, 1.7000, 2.0, 0.0004],

            [0.049, 12.000, 0.5000, 2.0, 0.0005],

            [0.062, 12.000, 2.4000, 2.0, 0.0006],

            [0.073, 12.000, 0.6000, 2.0, 0.0007],

            [0.079, 12.000, 0.9000, 2.0, 0.0006],

            [0.091, 12.000, 1.3000, 2.0, 0.0009],

            [0.099, 12.000, 1.4000, 2.0, 0.00015],

        ]

    )
@@ -313,6 +313,7 @@ def test_1d_weighted_binning(): 
        intensity=test_qiew_array[:, 1],

        error=test_qiew_array[:, 2],

        mod_q=test_qiew_array[:, 0],

        delta_mod_q=test_qiew_array[:, 4],

        wavelength=test_qiew_array[:, 3],

    )
@@ -326,18 +327,18 @@ def test_1d_weighted_binning(): 
    # Expected binned I(Q, wl)

    expected_binned_qiew = np.array(

        [

            [0.025, 10, 0.887045495441276, 1.0],

            [0.075, 10, 0.435609182296287, 1.0],

            [0.025, 12, 0.448133161649441, 2.0],

            [0.075, 12, 0.434869885395938, 2.0],

            [0.025, 10, 0.887045495441276, 1.0, 0.00016502890173410406],

            [0.075, 10, 0.435609182296287, 1.0, 0.0006155217563843883],

            [0.025, 12, 0.448133161649441, 2.0, 0.0004905877304625013],

            [0.075, 12, 0.434869885395938, 2.0, 0.0006426826759447178],

        ]

    )



    # Expected binned I(Q)

    expected_binned_qie = np.array(

        [

            [0.025, 11.5933404636534, 0.399987461455165],

            [0.075, 11.0016985965419, 0.307760492837406],

            [0.025, 11.5933404636534, 0.399987461455165, 0.00042439192929037884],

            [0.075, 11.0016985965419, 0.307760492837406, 0.0006291252838865729],

        ]

    )
@@ -355,6 +356,7 @@ def test_1d_weighted_binning(): 
    np.testing.assert_allclose(binned_iq_all_wl.mod_q, expected_binned_qie[:, 0])

    np.testing.assert_allclose(binned_iq_all_wl.intensity, expected_binned_qie[:, 1])

    np.testing.assert_allclose(binned_iq_all_wl.error, expected_binned_qie[:, 2])

    np.testing.assert_allclose(binned_iq_all_wl.delta_mod_q, expected_binned_qie[:, 3])



    # Do weighted binning on each wavelength

    binned_iq_per_wl = bin_intensity_into_q1d(
@@ -369,6 +371,7 @@ def test_1d_weighted_binning(): 
    np.testing.assert_allclose(binned_iq_per_wl.intensity, expected_binned_qiew[:, 1])

    np.testing.assert_allclose(binned_iq_per_wl.error, expected_binned_qiew[:, 2])

    np.testing.assert_allclose(binned_iq_per_wl.wavelength, expected_binned_qiew[:, 3])

    np.testing.assert_allclose(binned_iq_per_wl.delta_mod_q, expected_binned_qiew[:, 4])



    # Do weighted binning on Q-binned I(Q, wl)

    binned_iq_all_wl = bin_intensity_into_q1d(
@@ -378,6 +381,7 @@ def test_1d_weighted_binning(): 
    np.testing.assert_allclose(binned_iq_all_wl.mod_q, expected_binned_qie[:, 0])

    np.testing.assert_allclose(binned_iq_all_wl.intensity, expected_binned_qie[:, 1])

    np.testing.assert_allclose(binned_iq_all_wl.error, expected_binned_qie[:, 2])

    np.testing.assert_allclose(binned_iq_all_wl.delta_mod_q, expected_binned_qie[:, 3])





if __name__ == "__main__":