Update PaalmanPingsAbsorptionCorrection

[searscr]

Description of work

Do not merge before #39208

Updates PaalmanPingsAbsorptionCorrection to determine the GaugeVolume from the beam parameters provided from SetBeam. The following priority will be used to determine the integration volume for L1 distances:

1. Using DefineGaugeVolume or setting a GaugeVolume property for the workspace with a valid XML string defining shape.
2. Using SetBeam to find the intersection with the beam properties and the sample/container shape
3. Use the entire sample/container shape.

Summary of work

#39135 Created the functionality to get the intersection volume of the beam and the sample shape. This updates PaalmanPingsAbsorptionCorrection to use the new functionality to provide another method of determining the GaugeVolume.

Relates to EWM 10071

Story Summary:
modify the function and add in tests to demonstrate that it works. The tests are defined in the epic/parent story.

The gauge volume will be determined by the first one of the following that exists:

1. If the gauge volume is defined by DefineGaugeVolume
2. If the beam is defined, use the new functionality in MantidGeometry to calculate the gauge volume
3. Use the whole IObject as the gauge volume

Further detail of work

To test:

from mantid.simpleapi import *
import matplotlib.pyplot as plt
import numpy as np

ws = CreateWorkspace(DataX=[1.7981, 1.7983], DataY=[0.0]*4, NSpec=4, UnitX="Wavelength", YUnitLabel="Counts")
EditInstrumentGeometry(ws,
    PrimaryFlightPath=5.0,
    SpectrumIDs=[1, 2, 3, 4],
    L2=[2.0, 2.0, 2.0, 2.0],
    Polar=[10.0, 90.0, 170.0, 90.0],
    Azimuthal=[0.0, 0.0, 0.0, 45.0],
    DetectorIDs=[1, 2, 3, 4],
    InstrumentName="Instrument")
                  
SetSample(ws,
    Geometry={"Shape": "Cylinder", "Height": 5.68, "Radius": 0.295},
    Material={"ChemicalFormula": "La-(B11)5.94-(B10)0.06", "SampleNumberDensity": 0.1},
    ContainerGeometry={"Shape": 'HollowCylinder', 'Height': 5.68, 'InnerRadius': 0.295, 'OuterRadius': 0.315},
    ContainerMaterial={"ChemicalFormula": "V", "SampleNumberDensity": 0.0721})

# SetBeam to define the gauge volume
SetBeam(ws, Geometry={'Shape': 'Slit', 'Width': 3.0, 'Height': 3.0})
abs = PaalmanPingsAbsorptionCorrection(ws, ElementSize=0.1)
for a in ["ass", "assc", "acc", "acsc"]:
    print("{:4} {:.9f}(θ=10) {:.9f}(θ=90) {:.9f}(θ=170) {:.9f}(θ=90,φ=45)".format(a, *(mtd["abs_"+a].readY(i)[0] for i in range(4))))

---

Reviewer

Please comment on the points listed below (full description).
Your comments will be used as part of the gatekeeper process, so please comment clearly on what you have checked during your review. If changes are made to the PR during the review process then your final comment will be the most important for gatekeepers. In this comment you should make it clear why any earlier review is still valid, or confirm that all requested changes have been addressed.

Code Review

• Is the code of an acceptable quality?
• Does the code conform to the coding standards?
• Are the unit tests small and test the class in isolation?
• If there is GUI work does it follow the GUI standards?
• If there are changes in the release notes then do they describe the changes appropriately?
• Do the release notes conform to the release notes guide?

Functional Tests

• Do changes function as described? Add comments below that describe the tests performed?
• Do the changes handle unexpected situations, e.g. bad input?
• Has the relevant (user and developer) documentation been added/updated?

Does everything look good? Mark the review as Approve. A member of @mantidproject/gatekeepers will take care of it.

Gatekeeper

If you need to request changes to a PR then please add a comment and set the review status to "Request changes". This will stop the PR from showing up in the list for other gatekeepers.

[Kvieta1990]

I was using the following script (which I think is pretty much the same as the one put above in the PR) for testing,

ws = CreateWorkspace(DataX=[1.7981, 1.7983], DataY=[0.0]*4, NSpec=4, UnitX="Wavelength", YUnitLabel="Counts")
EditInstrumentGeometry(ws,
    PrimaryFlightPath=5.0,
    SpectrumIDs=[1, 2, 3, 4],
    L2=[2.0, 2.0, 2.0, 2.0],
    Polar=[10.0, 90.0, 170.0, 90.0],
    Azimuthal=[0.0, 0.0, 0.0, 45.0],
    DetectorIDs=[1, 2, 3, 4],
    InstrumentName="Instrument")
SetSample(ws,
    Geometry={"Shape": "Cylinder", "Height": 5.68, "Radius": 0.295},
    Material={"ChemicalFormula": "La-(B11)5.94-(B10)0.06", "SampleNumberDensity": 0.1},
    ContainerGeometry={"Shape": 'HollowCylinder', 'Height': 5.68, 'InnerRadius': 0.295, 'OuterRadius': 0.315},
    ContainerMaterial={"ChemicalFormula": "V", "SampleNumberDensity": 0.0721})
SetBeam(ws, Geometry={'Shape': 'Slit', 'Width': 100.0, 'Height': 200.0})
abs = PaalmanPingsAbsorptionCorrection(ws)
for a in ["ass", "assc", "acc", "acsc"]:
    print("{:4} {:.4f}(θ=10) {:.4f}(θ=90) {:.4f}(θ=170) {:.4f}(θ=90,φ=45)".format(a, *(mtd["abs_"+a].readY(i)[0] for i in range(4))))

Here, you can see that I was using a huge beam size. The thing is, when the beam size is larger than the sample size in every dimension, the calculation result should be the same as for the case where we don't have the beam size defined at all. However, the result created above is not the same as with the script in which I commented out the SetBeam line. Another observation is, as long as the beam size covers the whole sample, indeed the result no longer changes no matter how large we make the beam size be.

I am suspecting this may be something to do with the integration volume setting for sample and container. Maybe we were missing the intersecting volume setting for the container?

[searscr]

@Kvieta1990 I made an update that should fix the issue you are seeing with the beam illuminating the whole sample not giving the same result.

[Kvieta1990]

@Kvieta1990 I made an update that should fix the issue you are seeing with the beam illuminating the whole sample not giving the same result.

Thank you! I will grab the change and do the test.

[Kvieta1990]

I check the latest push and now having a super big beam size does create consistent outputs compared to the no beam defined case. Thank you!

I have another thing. When setting the beam size to be,

SetBeam(ws, Geometry={'Shape': 'Slit', 'Width': 1.0, 'Height': 2.0})

I got the output like this,

ass  0.0387(θ=10) 0.0522(θ=90) 0.0653(θ=170) 0.0389(θ=90,φ=45)
assc 0.0374(θ=10) 0.0504(θ=90) 0.0632(θ=170) 0.0373(θ=90,φ=45)
acc  0.2781(θ=10) 0.2785(θ=90) 0.2790(θ=170) 0.2751(θ=90,φ=45)
acsc 0.0974(θ=10) 0.1362(θ=90) 0.1886(θ=170) 0.1274(θ=90,φ=45)

Same script as mentioned above was used here.

Compared to the no beam defined situation (i.e., assuming the whole sample exposed to the beam), the difference seems to be a bit too large to me. I don't have a good expectation based on a solid physics ground, since the evaluation of the coefficients involves a normalization over the integration volume and as we reduce the integration volume size, the value on the top gets smaller but the normalization factor at the bottom also gets smaller (see the Eqn. (8) here). But the observed difference being so big makes me wonder whether we do the normalization correctly here. If I look at the codes here, I found the division is done over the sample volume (from the variable used there m_sampleVolume and the comment, I think it is referring to the total volume of the sample). But from my understanding of the absorption coefficient evaluation, the normalization should be done over the integration volume but NOT the total sample volume. Before we put in the beam size definition and gauge volume implementation, this is not a problem since integration volume and the sample volume are the same. But now, they are no longer necessarily the same so we have to check and make sure we do the normalization correctly.

We can discuss about this further. Unfortunately I have a full morning event tomorrow so I cannot go to the status meeting. But I will be free in the afternoon. So if you want to discuss, I can try to tune in later tomorrow. Thanks!

[ekapadi]

I finished with the "non test code" part of the code review. There are only a few minor things that I found. I'll go over the test code, and actually build it and run the tests early tomorrow.

[Kvieta1990]

@searscr I am wondering your thoughts about my very last comment about the volume normalization? Please feel free to let me know if you want to discuss. Thank you!

[searscr]

@Kvieta1990 Here we are setting the "sample volume" to the calculated volume from the Raster. From what I can tell with the equation you provided, it seems to be correct.

[Kvieta1990]

@Kvieta1990 Here we are setting the "sample volume" to the calculated volume from the Raster. From what I can tell with the equation you provided, it seems to be correct.

I am not sure actually. When I see m_sampleVolume = raster.totalvolume;, I would assume it means the total volume of the whole sample? But what we actually should have is the intersection volume with the beam.

For example, we can do a little test to print out the m_sampleVolume for the two situations -- one with beam defined and the other without beam defined. The m_sampleVolume value should be different, if the defined beam size is smaller than the sample.

[searscr]

@Kvieta1990 You were right, I had the sample and the integration volume backwards in the function to calculate the Raster volume so it was using the full sample size. I have updated the algorithm and now running the same test you ran before with

SetBeam(ws, Geometry={'Shape': 'Slit', 'Width': 1.0, 'Height': 2.0})

The result is

ass  0.146831240(θ=10) 0.197244313(θ=90) 0.250620450(θ=170) 0.156742892(θ=90,φ=45)
assc 0.140854936(θ=10) 0.189831094(θ=90) 0.241162035(θ=170) 0.149544197(θ=90,φ=45)
acc  0.944069217(θ=10) 0.944432246(θ=90) 0.944598198(θ=170) 0.933454398(θ=90,φ=45)
acsc 0.322538050(θ=10) 0.419365979(θ=90) 0.576807534(θ=170) 0.392303011(θ=90,φ=45)

compared to the full sample

ass  0.146621897(θ=10) 0.197750467(θ=90) 0.251731382(θ=170) 0.162254627(θ=90,φ=45)
assc 0.140687118(θ=10) 0.190336654(θ=90) 0.242260149(θ=170) 0.155058082(θ=90,φ=45)
acc  0.942924283(θ=10) 0.942705426(θ=90) 0.943423084(θ=170) 0.932408403(θ=90,φ=45)
acsc 0.325109461(θ=10) 0.421832405(θ=90) 0.577852023(θ=170) 0.401417991(θ=90,φ=45)

Does that look more correct?

[Kvieta1990]

That looks more reasonable! Thank you so much @searscr!

I will check later today to see whether I have more comments and I will let you know.

[ekapadi]

I built the code and ran through all of the tests as described (, including the details in the discussion with the CIS about what's supposed to happen when the beam profile is larger than the sample). I am happy with the quality of the changes, and as far as I can tell, everything is working as expected.

[Kvieta1990]

@searscr Sorry for the delayed response. I did check the PR, running the tests we did before and it seems to be all fine. Thank you so much!

[github-actions]

👋 Hi, @searscr,

Conflicts have been detected against the base branch. Please rebase your branch against the base branch.

[ekapadi]

Re-approved after rebase.