Running simulations¶
Simulations in RAMP are executed using Python. The architecture of RAMP is contained in the mcramp Python package, and is responsible for handling all of the OpenCL operations required to parallelize instrument simulations.
A detailed example¶
The basic workflow for RAMP is to write an instrument definition file for your instrument, run the simulation in a Python script, and then analyze the resulting data using standard scientific computing libraries for Python (SciPy, etc.).
The following notebook demonstrates this basic workflow for a simplified powder diffractometer, consisting of a monochromated source, and a sample with a single Bragg peak, described in the instrument definition file powder.json.
The snippet below contains a small two line initialization of the OpenCL internals required to identify the device on which the simulations should run. For simulations being run on personal workstations, this will be sufficient and can be safely copied and pasted into any RAMP script.
[1]:
import mcramp as mcr
import pyopencl as cl
import numpy as np
ctx = cl.create_some_context()
queue = cl.CommandQueue(ctx)
Loading instruments¶
The instrument powder.json contains the variables Ei and Mono_angle which we can use to vary the wavelength of incident neutrons. First we shall calculate some values for these variables.
[2]:
# Physical constants
h=6.62607015e-34
mn=1.674929e-27
meV=1.602176634e-22
overAA = 1e10
Ei = 9.09 # meV - 3AA neutrons
Li = h / np.sqrt(2*mn*Ei*meV) * overAA
Mono_d_spacing = 3.355 # AA
Mono_angle = np.arcsin(Li / (2 * Mono_d_spacing))
The instrument definition file also contains the variable detector_binning
which we can use to vary the histogram bins in the detector output.
[3]:
detector_binning = [-40.0, 1.0, 140.0]
Instruments in RAMP are represented by the Instrument
class, which contains the methods responsible for instantiating, simulating, and obtaining the results of simulations. Instances of the Instrument class are instantiated as follows:
[4]:
inst = mcr.Instrument(
'powder.json',
ctx,
queue,
Ei=Ei,
Mono_angle=Mono_angle,
detector_binning=detector_binning
)
C:\Users\wyz83832\AppData\Local\Continuum\anaconda3\lib\site-packages\pyopencl\__init__.py:235: CompilerWarning: Non-empty compiler output encountered. Set the environment variable PYOPENCL_COMPILER_OUTPUT=1 to see more.
"to see more.", CompilerWarning)
The constructor for the Instrument
class takes three arguments by default: the filename of the instrument definition file, the OpenCL context, and the OpenCL queue (the internals we set up earlier). The additional keyword arguments Ei
and Mono_angle
correspond to variables in the instrument definition file.
Executing simulations¶
Simulations are executed by calling the execute(N)
method of the Instrument
class, where the argument N
is the desired number of neutron histories to be simulated.
[5]:
N = int(1e7)
inst.execute(N)
Analysing simulation results¶
The easiest way to get a quick visual confirmation that a simulation has returned sensible results is by plotting the output of the instrument monitors using the plot()
method of the Instrument
class. This will do some default plotting using matplotlib
of the histograms of all the detectors in the instrument.
[6]:
inst.plot()

We can see that there is a peak at approximately 60 degrees. Let’s adjust the detector binning and re-run the simulation to obtain a better view of the peak shape.
[7]:
detector_binning = [50.0, 0.02, 70.0]
inst = mcr.Instrument(
'powder.json',
ctx,
queue,
Ei=Ei,
Mono_angle=Mono_angle,
detector_binning=detector_binning
)
N = int(1e7)
inst.execute(N)
inst.plot()

Pretty ugly, but not surprising for such a crude instrument. To perform a more detailed analysis of the detector output, we can obtain the histogram using the data()
method of the Instrument
class. This method returns a dictionary whose keys are the names of components in the instrument definition and whose values are data produced by the corresponding component. Typically this will be non-zero for detectors which provide their axes and histogrammed intensities.
[8]:
data = inst.data()
detector_output = data['Detector']
theta_axis = detector_output[0]
intensity = detector_output[1]
With this data in hand, we can now perform some analysis. To find the peak centre, let’s naievely fit our lineshape to a Gaussian distribution and find the position of the mean.
[9]:
from scipy.optimize import curve_fit
func = lambda x, c, mu, sigma: c*np.exp(-(x - mu)**2 / (2*sigma**2))
popt, pcov = curve_fit(func, xdata=theta_axis, ydata=intensity, p0=[1e10, 60.0, 2.5])
print(popt)
[8.84435273e+09 6.03711688e+01 1.02435847e+00]
Giving a peak center of 60.4 degrees (as compared to the nominal answer of 60 degrees exactly for our chosen sample and instrument parameters). Plotting the resulting fit yields…
[10]:
import matplotlib.pyplot as plt
plt.plot(theta_axis, intensity, 'b-')
plt.plot(theta_axis, func(theta_axis, *popt), 'r-')
plt.ylabel("Intensity")
plt.xlabel("Theta [deg]")
plt.show()

To save the results of a simulation for analysis at a later date, the Instrument
class provides the save()
method. Any component with a filename
argument will save it’s data()
output as a numpy file filename.npy
when inst.save()
is called.
Instrument class reference¶
-
class
mcramp.
Instrument
(fn, ctx, queue, **kwargs)[source]¶ Core class which instantiates instrument definiton files, and provides public methods for executing simulations, and plotting, saving, and analyzing results
- Parameters
- fnstr
Instrument definition file name
- ctxpyopencl.Context
OpenCL context within which the simulation is to be executed
- queuepyopencl.CommandQueue
OpenCL command queue for enqueueing simulation kernels
- **kwargs
Values for variables contained in the instrument definition file
Methods
data
(self)Returns a dictionary whose keys are the names of components in the instrument definition file, and whose values are the data arrays returned by the respective component
execute
(self, N[, debug])Executes the instrument simulation
plot
(self)Calls the plotting function of each component and displays the output
save
(self)Saves the data array returned by each component’s data() function to a numpy file “filename.npy”, where filename is the value of the filename attribute in the instrument definition file