Skip to content
Snippets Groups Projects
Commit 47e1dfba authored by Doucet, Mathieu's avatar Doucet, Mathieu
Browse files

Re #9159 first cut for the USANS geometry

parent 0541871e
No related branches found
No related tags found
No related merge requests found
"""*WIKI*
Simulate a USANS workspace
*WIKI*"""
from mantid.simpleapi import *
from mantid.api import *
from mantid.kernel import *
import math
import numpy
class USANSSimulation(PythonAlgorithm):
def category(self):
return "SANS"
def name(self):
return "USANSSimulation"
def PyInit(self):
self.declareProperty("TwoTheta", 0.01, "Scattering angle in degrees")
self.declareProperty(FloatArrayProperty("WavelengthPeaks", values=[0.9, 1.2, 1.8, 3.6],
direction=Direction.Input), "Wavelength peaks out of the monochromator")
# Model parameters
self.declareProperty("SphereRadius", 200.0, "Radius for the sphere model (Angstrom)")
self.declareProperty("Background", 0.0, "Background")
self.declareProperty("SigmaPeak", 0.01, "Width of the wavelength peaks")
self.declareProperty(MatrixWorkspaceProperty("OutputWorkspace", "", Direction.Output), "Output workspace")
self.declareProperty(MatrixWorkspaceProperty("MonitorWorkspace", "", Direction.Output), "Output monitor workspace")
def PyExec(self):
workspace = self.getPropertyValue("OutputWorkspace")
out_ws = CreateSimulationWorkspace(Instrument="USANS",
BinParams="0,50,32000",
UnitX="TOF",
OutputWorkspace=workspace)
data_x = out_ws.dataX(0)
mon_ws_name = self.getPropertyValue("MonitorWorkspace")
mon_ws = CreateWorkspace(dataX=data_x, dataY=numpy.zeros(len(data_x)-1),
UnitX="TOF", OutputWorkspace=mon_ws_name)
mon_y = mon_ws.dataY(0)
mon_e = mon_ws.dataE(0)
# Number of pixels for the main detector
n_pixels = out_ws.getNumberHistograms()/2
# Clean up the workspace
for j in range(n_pixels):
data_y = out_ws.dataY(j)
for i in range(len(data_y)):
data_y[i] = 0.0
# Fill monitor workspace with fake beam profile
for i in range(len(data_x)-1):
wl_i = 0.0039560/30.0*(data_x[i]+data_x[i+1])/2.0
mon_y[i] = 1000.0*math.exp(-wl_i)
mon_e[i] = math.sqrt(mon_y[i])
# Add analyzer theta value and monochromator angle theta_b in logs
two_theta = self.getProperty("TwoTheta").value
theta_b = 70.0
theta = theta_b + two_theta
out_ws.getRun().addProperty("AnalyzerTheta", theta, 'degree', True)
out_ws.getRun().addProperty("TwoTheta", two_theta, 'degree', True)
out_ws.getRun().addProperty("MonochromatorTheta", theta_b, 'degree', True)
# List of wavelength peaks, and width of the peaks
wl_peaks = self.getProperty("WavelengthPeaks").value
sigma = self.getProperty("SigmaPeak").value
for wl in wl_peaks:
q = 6.28*math.sin(two_theta)/wl
Logger.get("USANS").notice( "wl = %g; Q = %g" % (wl, q))
for i in range(len(data_x)-1):
wl_i = 0.0039560/30.0*(data_x[i]+data_x[i+1])/2.0
# Scale the I(q) by a Gaussian to simulate the wavelength peaks selected by the monochromator
flux = 1.0e6/(sigma*math.sqrt(2.0*math.pi))*math.exp(-(wl_i-wl)*(wl_i-wl)/(2.0*sigma*sigma))
# Multiply by beam profile
flux *= mon_y[i]
# Account for transmission
flux *= math.exp(-wl_i/2.0)
# Transmission detector
for j in range(n_pixels, 2*n_pixels):
det_pos = out_ws.getInstrument().getDetector(j).getPos()
r = math.sqrt(det_pos.Y()*det_pos.Y()+det_pos.X()*det_pos.X())
sigma = 0.01
scale = math.exp(-r*r/(2.0*sigma*sigma))
data_y = out_ws.dataY(j)
data_y[i] += int(scale*flux)
data_e = out_ws.dataE(j)
data_e[i] = math.sqrt(data_e[i]*data_e[i]+scale*scale*flux*flux)
# Compute I(q) and store the results
q_i = q*wl/wl_i
i_q = self._sphere_model(q_i, scale=flux)
for j in range(n_pixels):
det_pos = out_ws.getInstrument().getDetector(j).getPos()
r = math.sqrt(det_pos.Y()*det_pos.Y()+det_pos.X()*det_pos.X())
sigma = 0.01
scale = math.exp(-r*r/(2.0*sigma*sigma))
data_y = out_ws.dataY(j)
data_y[i] += int(i_q*scale)
data_e = out_ws.dataE(j)
data_e[i] = math.sqrt(data_e[i]*data_e[i]+i_q*i_q*scale*scale)
self.setProperty("OutputWorkspace", out_ws)
self.setProperty("MonitorWorkspace", mon_ws)
def _sphere_model(self, q, scale):
"""
Return I(q) for a sphere model
@param q: q-value
@param scale: normalization factor to give I(q)
"""
radius = self.getProperty("SphereRadius").value
bck = self.getProperty("Background").value
qr = q*radius
bes = 3.0*(math.sin(qr)-qr*math.cos(qr))/(qr*qr*qr) if not qr == 0.0 else 1.0
vol = 4.0*math.pi/3.0*radius*radius*radius
f2 = vol*bes*bes*1.0e-6
return(scale*f2+bck)
#############################################################################################
AlgorithmFactory.subscribe(USANSSimulation())
...@@ -360,6 +360,10 @@ ...@@ -360,6 +360,10 @@
<technique>Small Angle Scattering</technique> <technique>Small Angle Scattering</technique>
</instrument> </instrument>
<instrument name="USANS" shortname="USANS" beamline="1A">
<technique>Small Angle Scattering</technique>
</instrument>
<instrument name="VULCAN" beamline="7"> <instrument name="VULCAN" beamline="7">
<technique>Neutron Diffraction</technique> <technique>Neutron Diffraction</technique>
</instrument> </instrument>
......
<?xml version="1.0" encoding="utf-8"?>
<instrument xmlns="http://www.mantidproject.org/IDF/1.0" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.mantidproject.org/IDF/1.0 Schema/IDFSchema.xsd"
name="USANS" valid-from="1900-01-31 23:59:59" valid-to="2100-01-31 23:59:59" last-modified="2014-03-10 08:53:09">
<defaults>
<length unit="meter" />
<angle unit="degree" />
<reference-frame>
<along-beam axis="z" />
<pointing-up axis="y" />
<handedness val="right" />
</reference-frame>
</defaults>
<component type="moderator">
<location z="-28.011" />
</component>
<type name="moderator" is="Source"></type>
<!-- Sample position -->
<component type="sample-position">
<location y="0.0" x="0.0" z="0.0" />
</component>
<type name="sample-position" is="SamplePos" />
<!-- Moderator - monitor distance: 27606 mm -->
<component type="monitors" idlist="monitors">
<location/>
</component>
<type name="monitors">
<component type="monitor">
<location z="-0.405" name="monitor1"/>
</component>
</type>
<!-- Detector Panels -->
<component type="main_detector" idlist="main_detector-idlist">
<!-- location x="0.0" y="0.0" z="0.0" rot="0.0" axis-x="0.0" axis-y="0.0" axis-z="1.0"/ -->
<location/>
</component>
<type name="main_detector">
<component type="bank1">
<location/>
</component>
<component type="bank2">
<location/>
</component>
</type>
<type name="bank1">
<component type="fourpack">
<location y="0.0" x="0.0" z="2.03175" />
</component>
</type>
<type name="bank2">
<component type="fourpack">
<location y="0.0" x="-0.0096266" z="2.04175" />
</component>
</type>
<component type="transmission_detector" idlist="transmission_detector-idlist">
<location/>
</component>
<type name="transmission_detector">
<component type="fourpack">
<location y="0.0" x="0.0" z="0.710316" />
</component>
<component type="fourpack">
<location y="0.0" x="-0.0096266" z="0.720316" />
</component>
</type>
<!-- 4-PACK -->
<type name="fourpack">
<properties/>
<component type="tube">
<location x="-0.0288798" name="tube1"/>
<location x="-0.0096266" name="tube2"/>
<location x="0.0096266" name="tube3"/>
<location x="0.0288798" name="tube4"/>
</component>
</type>
<!--TUBE-->
<type name="tube" outline="yes">
<properties/>
<component type="pixel">
<location y="-0.10998200" name="pixel0"/>
<location y="-0.10183519" name="pixel1"/>
<location y="-0.09368837" name="pixel2"/>
<location y="-0.08554156" name="pixel3"/>
<location y="-0.07739474" name="pixel4"/>
<location y="-0.06924793" name="pixel5"/>
<location y="-0.06110111" name="pixel6"/>
<location y="-0.05295430" name="pixel7"/>
<location y="-0.04480748" name="pixel8"/>
<location y="-0.03666067" name="pixel9"/>
<location y="-0.02851385" name="pixel10"/>
<location y="-0.02036704" name="pixel11"/>
<location y="-0.01222022" name="pixel12"/>
<location y="-0.00407341" name="pixel13"/>
<location y="0.00407341" name="pixel14"/>
<location y="0.01222022" name="pixel15"/>
<location y="0.02036704" name="pixel16"/>
<location y="0.02851385" name="pixel17"/>
<location y="0.03666067" name="pixel18"/>
<location y="0.04480748" name="pixel19"/>
<location y="0.05295430" name="pixel20"/>
<location y="0.06110111" name="pixel21"/>
<location y="0.06924793" name="pixel22"/>
<location y="0.07739474" name="pixel23"/>
<location y="0.08554156" name="pixel24"/>
<location y="0.09368837" name="pixel25"/>
<location y="0.10183519" name="pixel26"/>
<location y="0.10998200" name="pixel27"/>
</component>
</type>
<!--PIXEL FOR TUBE-->
<type is="detector" name="pixel">
<cylinder id="cyl-approx">
<centre-of-bottom-base p="0.0" r="0.0" t="0.0"/>
<axis y="1.0" x="0.0" z="0.0"/>
<radius val="0.00635"/>
<height val="0.0093741875"/>
</cylinder>
<algebra val="cyl-approx"/>
</type>
<!--MONITOR SHAPE-->
<type is="monitor" name="monitor">
<cuboid id="shape">
<left-front-bottom-point y="-0.08255" x="-0.0254" z="-0.01905"/>
<left-front-top-point y="0.08255" x="-0.0254" z="-0.01905"/>
<left-back-bottom-point y="-0.08255" x="-0.0254" z="0.01905"/>
<right-front-bottom-point y="-0.08255" x="0.0254" z="-0.01905"/>
</cuboid>
<algebra val="shape"/>
</type>
<!--DETECTOR IDs-->
<idlist idname="main_detector-idlist">
<id start="0" end="223"/>
</idlist>
<idlist idname="transmission_detector-idlist">
<id start="224" end="447"/>
</idlist>
<!--MONITOR IDs-->
<idlist idname="monitors">
<id val="-1"/>
</idlist>
<!--DETECTOR PARAMETERS-->
<component-link name="USANS">
<parameter name="wavelength_peaks" type="string">
<value val="0.9, 1.2, 1.8, 3.6" />
</parameter>
</component-link>
</instrument>
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment