diff --git a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksProfileFitting.py b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksProfileFitting.py
index a3418ba4f677720fbf4eaac5417a8e341ddf0af5..111b06c4dc39c8492bbc30c1b89a32867e1733d6 100644
--- a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksProfileFitting.py
+++ b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksProfileFitting.py
@@ -43,12 +43,9 @@ class IntegratePeaksProfileFitting(PythonAlgorithm):
self.declareProperty("RunNumber", defaultValue=0,
doc="Run Number to integrate")
- self.declareProperty("DQPixel", defaultValue=0.003, validator=FloatBoundedValidator(lower=0., exclusive=True),
- doc="The side length of each voxel in the non-MD histogram used for fitting (1/Angstrom)")
-
self.declareProperty(FileProperty(name="UBFile",defaultValue="",action=FileAction.OptionalLoad,
extensions=[".mat"]),
- doc="File containing the UB Matrix in ISAW format.")
+ doc="File containing the UB Matrix in ISAW format. Leave blank to use loaded UB Matrix.")
self.declareProperty(FileProperty(name="ModeratorCoefficientsFile",
defaultValue="",action=FileAction.OptionalLoad,
extensions=[".dat"]),
@@ -59,22 +56,13 @@ class IntegratePeaksProfileFitting(PythonAlgorithm):
self.declareProperty("IntensityCutoff", defaultValue=0., doc="Minimum number of counts to force a profile")
edgeDocString = 'Pixels within EdgeCutoff from a detector edge will be have a profile forced. Currently for 256x256 cameras only.'
self.declareProperty("EdgeCutoff", defaultValue=0., doc=edgeDocString)
- self.declareProperty("FracHKL", defaultValue=0.5, validator=FloatBoundedValidator(lower=0., exclusive=True),
- doc="Fraction of HKL to consider for profile fitting.")
self.declareProperty("FracStop", defaultValue=0.05, validator=FloatBoundedValidator(lower=0., exclusive=True),
doc="Fraction of max counts to include in peak selection.")
self.declareProperty("MinpplFrac", defaultValue=0.9, doc="Min fraction of predicted background level to check")
self.declareProperty("MaxpplFrac", defaultValue=1.1, doc="Max fraction of predicted background level to check")
- mindtBinWidthDocString = "Smallest spacing (in microseconds) between data points for TOF profile fitting."
- self.declareProperty("MindtBinWidth", defaultValue=15, doc=mindtBinWidthDocString)
-
- self.declareProperty("NTheta", defaultValue=50, doc="Number of bins for bivarite Gaussian along the scattering angle.")
- self.declareProperty("NPhi", defaultValue=50, doc="Number of bins for bivariate Gaussian along the azimuthal angle.")
self.declareProperty("DQMax", defaultValue=0.15, doc="Largest total side length (in Angstrom) to consider for profile fitting.")
- self.declareProperty("DtSpread", defaultValue=0.03, validator=FloatBoundedValidator(lower=0., exclusive=True),
- doc="The fraction of the peak TOF to consider for TOF profile fitting.")
self.declareProperty("PeakNumber", defaultValue=-1, doc="Which Peak to fit. Leave negative for all.")
def PyExec(self):
@@ -83,10 +71,8 @@ class IntegratePeaksProfileFitting(PythonAlgorithm):
from mantid.simpleapi import LoadIsawUB
import pickle
from scipy.ndimage.filters import convolve
-
MDdata = self.getProperty('InputWorkspace').value
peaks_ws = self.getProperty('PeaksWorkspace').value
- fracHKL = self.getProperty('FracHKL').value
fracStop = self.getProperty('FracStop').value
dQMax = self.getProperty('DQMax').value
UBFile = self.getProperty('UBFile').value
@@ -96,11 +82,17 @@ class IntegratePeaksProfileFitting(PythonAlgorithm):
edgeCutoff = self.getProperty('EdgeCutoff').value
peakNumberToFit = self.getProperty('PeakNumber').value
- LoadIsawUB(InputWorkspace=peaks_ws, FileName=UBFile)
+ if UBFile == '' and peaks_ws.sample().hasOrientedLattice():
+ logger.information("Using UB file already available in PeaksWorkspace")
+ else:
+ try:
+ LoadIsawUB(InputWorkspace=peaks_ws, FileName=UBFile)
+ except:
+ logger.error("peaks_ws does not have a UB matrix loaded. Must provide a file")
+
UBMatrix = peaks_ws.sample().getOrientedLattice().getUB()
dQ = np.abs(ICCFT.getDQFracHKL(UBMatrix, frac=0.5))
dQ[dQ>dQMax] = dQMax
- dQPixel = self.getProperty('DQPixel').value
q_frame='lab'
mtd['MDdata'] = MDdata
@@ -113,21 +105,40 @@ class IntegratePeaksProfileFitting(PythonAlgorithm):
else:
strongPeakParams = None #This will not force any profiles
- nTheta = self.getProperty('NTheta').value
- nPhi = self.getProperty('NPhi').value
zBG = 1.96
- mindtBinWidth = self.getProperty('MindtBinWidth').value
pplmin_frac = self.getProperty('MinpplFrac').value
pplmax_frac = self.getProperty('MaxpplFrac').value
sampleRun = self.getProperty('RunNumber').value
+
+ # There are a few instrument specific parameters that we define here. In some cases,
+ # it may improve fitting to set tweak these parameters, but for simplicity we define these here
+ # The default values are good for MaNDi - new instruments can be added by adding a different elif
+ # statement.
+ # If you change these values or add an instrument, documentation should also be changed.
+ try:
+ nTheta = peaks_ws.getInstrument().getIntParameter("numBinsTheta")[0]
+ nPhi = peaks_ws.getInstrument().getIntParameter("numBinsPhi")[0]
+ mindtBinWidth = peaks_ws.getInstrument().getNumberParameter("mindtBinWidth")[0]
+ maxdtBinWidth = peaks_ws.getInstrument().getNumberParameter("maxdtBinWidth")[0]
+ fracHKL = peaks_ws.getInstrument().getNumberParameter("fracHKL")[0]
+ dQPixel = peaks_ws.getInstrument().getNumberParameter("dQPixel")[0]
+ peakMaskSize = peaks_ws.getInstrument().getIntParameter("peakMaskSize")[0]
+
+ except:
+ raise
+ logger.error("Cannot find all parameters in instrument parameters file.")
+ sys.exit(1)
+
neigh_length_m=3
qMask = ICCFT.getHKLMask(UBMatrix, frac=fracHKL, dQPixel=dQPixel,dQ=dQ)
+ iccFitDict = ICCFT.parseConstraints(peaks_ws) #Contains constraints and guesses for ICC Fitting
+
numgood = 0
numerrors = 0
# Create the parameters workspace
- keys = ['peakNumber','Alpha', 'Beta', 'R', 'T0', 'bgBVG', 'chiSq3d', 'dQ', 'KConv', 'MuPH',
+ keys = ['peakNumber','Alpha', 'Beta', 'R', 'T0', 'bgBVG', 'chiSq3d', 'chiSq', 'dQ', 'KConv', 'MuPH',
'MuTH', 'newQ', 'Scale', 'scale3d', 'SigP', 'SigX', 'SigY', 'Intens3d', 'SigInt3d']
datatypes = ['float']*len(keys)
datatypes[np.where(np.array(keys)=='newQ')[0][0]] = 'V3D'
@@ -153,13 +164,16 @@ class IntegratePeaksProfileFitting(PythonAlgorithm):
box = ICCFT.getBoxFracHKL(peak, peaks_ws, MDdata, UBMatrix, peakNumber,
dQ, fracHKL=0.5, dQPixel=dQPixel, q_frame=q_frame)
# Will force weak peaks to be fit using a neighboring peak profile
- Y3D, goodIDX, pp_lambda, params = BVGFT.get3DPeak(peak, box, padeCoefficients,qMask,
+ Y3D, goodIDX, pp_lambda, params = BVGFT.get3DPeak(peak, peaks_ws, box, padeCoefficients,qMask,
nTheta=nTheta, nPhi=nPhi, plotResults=False,
zBG=zBG,fracBoxToHistogram=1.0,bgPolyOrder=1,
strongPeakParams=strongPeakParams,
q_frame=q_frame, mindtBinWidth=mindtBinWidth,
+ maxdtBinWidth=maxdtBinWidth,
pplmin_frac=pplmin_frac, pplmax_frac=pplmax_frac,
- forceCutoff=forceCutoff, edgeCutoff=edgeCutoff)
+ forceCutoff=forceCutoff, edgeCutoff=edgeCutoff,
+ peakMaskSize=peakMaskSize,
+ iccFitDict=iccFitDict)
# First we get the peak intensity
peakIDX = Y3D/Y3D.max() > fracStop
diff --git a/docs/source/algorithms/IntegratePeaksProfileFitting-v1.rst b/docs/source/algorithms/IntegratePeaksProfileFitting-v1.rst
index 1efc25edf46d931d73494296237e68d0744c20c3..6641913f21ccc541071976363ecc71625cf2e5c6 100644
--- a/docs/source/algorithms/IntegratePeaksProfileFitting-v1.rst
+++ b/docs/source/algorithms/IntegratePeaksProfileFitting-v1.rst
@@ -34,6 +34,43 @@ The algorithms takes two input workspaces:
- The OutputParamsWorkspace is a TableWorkspace containing the fit parameters.
Peaks which could not be fit are omitted.
+Instrument-Defined Parameters
+-----------------------------
+In addition to the input parameters defined above, there are several other parameters
+to be aware of which are pre-defined for each instrument. The instrument is determined
+from the instrument that is loaded into PeaksWorkspace. If the instrument parameters file
+does not contain paramters, the algorithm defaults to MaNDi parameters. Default
+values are below:
+
++--------------+----------------------------+----------+----------+---------+
+| Parameter | Description | MaNDi | TOPAZ | CORELLI |
++==============+============================+==========+==========+=========+
+| DQPixel | The side length for each | | | |
+| | voxel used for fitting. | 0.003 | 0.01 | 0.007 |
+| | Units: 1/Angstrom | | | |
++--------------+----------------------------+----------+----------+---------+
+| FracHKL | The distance between peaks | | | |
+| | (in fraction of hkl) that | 0.4 | 0.4 | 0.4 |
+| | is used for fitting. | | | |
++--------------+----------------------------+----------+----------+---------+
+| MinDtBinWidth| The smallest time bin used | | | |
+| | for fitting the TOF profile| 15 | 2 | 2 |
+| | Units: microseconds | | | |
++--------------+----------------------------+----------+----------+---------+
+| MaxDtBinWidth| The largest time bin used | | | |
+| | for fitting the TOF profile| 50 | 15 | 60 |
+| | Units: microseconds | | | |
++--------------+----------------------------+----------+----------+---------+
+| NTheta | The number of bins along | | | |
+| | the scattering direction | 50 | 50 | 50 |
+| | used for BVG fitting. | | | |
++--------------+----------------------------+----------+----------+---------+
+| NPhi | The number of bins along | | | |
+| | the azimuthal direction | 50 | 50 | 50 |
+| | used for BVG fitting. | | | |
++--------------+----------------------------+----------+----------+---------+
+
+
Calculations
------------
This algorithm will fit a set of peaks in a PeaksWorkspace. The intensity profile
@@ -114,10 +151,9 @@ Usage
LoadIsawPeaks(Filename='/SNS/MANDI/shared/ProfileFitting/demo_5921.integrate', OutputWorkspace='peaks_ws')
IntegratePeaksProfileFitting(OutputPeaksWorkspace='peaks_ws_out', OutputParamsWorkspace='params_ws',
- InputWorkspace='MANDI_5921_md', PeaksWorkspace='peaks_ws', RunNumber=5921, DtSpread=0.015,
- UBFile='/SNS/MANDI/shared/ProfileFitting/demo_5921.mat',
+ InputWorkspace='MANDI_5921_md', PeaksWorkspace='peaks_ws', RunNumber=5921,
+ UBFile='/SNS/MANDI/shared/ProfileFitting/demo_5921.mat', MinpplFrac=0.9, MaxpplFrac=1.1,
ModeratorCoefficientsFile='/SNS/MANDI/shared/ProfileFitting/franz_coefficients_2017.dat',
- MinpplFrac=0.9, MaxpplFrac=1.1, MindtBinWidth=15,
StrongPeakParamsFile='/SNS/MANDI/shared/ProfileFitting/strongPeakParams_beta_lac_mut_mbvg.pkl',
peakNumber=30)
diff --git a/docs/source/release/v3.14.0/diffraction.rst b/docs/source/release/v3.14.0/diffraction.rst
index a1cc728f84ecdcb97b7557595cc0ca681f434918..a1d75415248c3d8d43beffa6e050af3a35a991e2 100644
--- a/docs/source/release/v3.14.0/diffraction.rst
+++ b/docs/source/release/v3.14.0/diffraction.rst
@@ -10,3 +10,12 @@ Diffraction Changes
improvements, followed by bug fixes.
:ref:`Release 3.14.0 <v3.14.0>`
+
+
+Single Crystal Diffraction
+--------------------------
+
+Improvements
+############
+
+- :ref:`IntegratePeaksProfileFitting <algm-IntegratePeaksProfileFitting>` now supports MaNDi, TOPAZ, and CORELLI. Other instruments can easily be added as well.
diff --git a/instrument/CORELLI_Parameters.xml b/instrument/CORELLI_Parameters.xml
index e9cf3f940d4ea06da9e75f25eb91b28d9f51c7c2..1fde9c8970ce8946ff84970c19103b0c7095d366 100644
--- a/instrument/CORELLI_Parameters.xml
+++ b/instrument/CORELLI_Parameters.xml
@@ -14,6 +14,88 @@
<value val="(101.9 * incidentEnergy^(-0.41) * exp(-incidentEnergy/282.0))" />
</parameter>
+
+ <!-- Need to fill in gaps for peak profile fitting -->
+ <parameter name="fitConvolvedPeak" type="bool">
+ <value val="true"/>
+ </parameter>
+
+ <!-- Multiplier for profile fitting for BVG polar angle -->
+ <parameter name="sigX0Scale">
+ <value val="2." />
+ </parameter>
+
+ <!-- Multiplier for profile fitting for BVG azimuthal angle -->
+ <parameter name="sigY0Scale">
+ <value val="2." />
+ </parameter>
+
+ <!-- Number of rows between detector gaps for profile fitting -->
+ <parameter name="numDetRows" type="int">
+ <value val="255" />
+ </parameter>
+
+ <!-- Number of cols between detector gaps for profile fitting -->
+ <parameter name="numDetCols" type="int">
+ <value val="16" />
+ </parameter>
+
+ <!-- Number of polar bins for BVG histogramming when profile fitting -->
+ <parameter name="numBinsTheta" type="int">
+ <value val="50" />
+ </parameter>
+
+ <!-- Number of azimuthal bins for BVG histogramming when profile fitting -->
+ <parameter name="numBinsPhi" type="int">
+ <value val="50" />
+ </parameter>
+
+ <!-- Fraction along (h,k,l) to use for profile fitting. 0.5 is the next peak. -->
+ <parameter name="fracHKL">
+ <value val="0.4" />
+ </parameter>
+
+ <!-- Side length of each voxel for fitting in units of angstrom^-1 -->
+ <parameter name="dQPixel">
+ <value val="0.007" />
+ </parameter>
+
+ <!-- Minimum spacing for profile fitting the TOF profile. Units of microseconds -->
+ <parameter name="mindtBinWidth">
+ <value val="2" />
+ </parameter>
+
+ <!-- Maximum spacing for profile fitting the TOF profile. Units of microseconds -->
+ <parameter name="maxdtBinWidth">
+ <value val="60" />
+ </parameter>
+
+ <!-- Size of peak mask for background calculation in units of dQPixel -->
+ <parameter name="peakMaskSize" type="int">
+ <value val="10" />
+ </parameter>
+
+ <!-- Constraints for ICC fitting. Valid names are iccA, iccB, iccR, iccT0, iccScale0
+ iccHatWidth and iccKConv. Inputs are strings with values separated by
+ spaces which are prased by the IntegratePeaksProfileFitting algorithm.
+ If two values are given they are treated as the lower and upper bounds. If
+ three are given they are the lower bound, upper bound, and initial guess.-->
+ <parameter name="iccA" type="string">
+ <value val="0.25 0.75 0.5" />
+ </parameter>
+
+ <parameter name="iccB" type="string">
+ <value val="0.001 0.3 0.005" />
+ </parameter>
+
+ <parameter name="iccR" type="string">
+ <value val="0.05 0.15 0.1" />
+ </parameter>
+
+ <parameter name="iccKConv" type="string">
+ <value val="10.0 800.0 100.0" />
+ </parameter>
+
</component-link>
</parameter-file>
diff --git a/instrument/MANDI_Parameters.xml b/instrument/MANDI_Parameters.xml
index a617566959d86d03cf807d8d316bf02b43a8db00..9ac7e699e5c294a906f17ebcabf9ddf4a6a87296 100644
--- a/instrument/MANDI_Parameters.xml
+++ b/instrument/MANDI_Parameters.xml
@@ -8,6 +8,66 @@
<value val="1"/>
</parameter>
+<!-- Need to fill in gaps for peak profile fitting -->
+<parameter name="fitConvolvedPeak" type="bool">
+ <value val="true"/>
+</parameter>
+
+<!-- Multiplier for profile fitting for BVG polar angle -->
+<parameter name="sigX0Scale">
+ <value val="2." />
+</parameter>
+
+<!-- Multiplier for profile fitting for BVG azimuthal angle -->
+<parameter name="sigY0Scale">
+ <value val="2." />
+</parameter>
+
+<!-- Number of rows between detector gaps for profile fitting -->
+<parameter name="numDetRows" type="int">
+ <value val="255" />
+</parameter>
+
+<!-- Number of cols between detector gaps for profile fitting -->
+<parameter name="numDetCols" type="int">
+ <value val="255" />
+</parameter>
+
+<!-- Number of polar bins for BVG histogramming when profile fitting -->
+<parameter name="numBinsTheta" type="int">
+ <value val="50" />
+</parameter>
+
+<!-- Number of azimuthal bins for BVG histogramming when profile fitting -->
+<parameter name="numBinsPhi" type="int">
+ <value val="50" />
+</parameter>
+
+<!-- Fraction along (h,k,l) to use for profile fitting. 0.5 is the next peak. -->
+<parameter name="fracHKL">
+ <value val="0.4" />
+</parameter>
+
+<!-- Side length of each voxel for fitting in units of angstrom^-1 -->
+<parameter name="dQPixel">
+ <value val="0.003" />
+</parameter>
+
+<!-- Minimum spacing for profile fitting the TOF profile. Units of microseconds -->
+<parameter name="mindtBinWidth">
+ <value val="15" />
+</parameter>
+
+<!-- Maximum spacing for profile fitting the TOF profile. Units of microseconds -->
+<parameter name="maxdtBinWidth">
+ <value val="50" />
+</parameter>
+
+<!-- Size of peak mask for background calculation in units of dQPixel -->
+<parameter name="peakMaskSize" type="int">
+ <value val="5" />
+</parameter>
+
</component-link>
</parameter-file>
diff --git a/instrument/MANDI_Parameters_2015_08_01.xml b/instrument/MANDI_Parameters_2015_08_01.xml
index 52ba77f6fdbe47c82b3dc7d2b3fc290aec796878..2de3d3e50778bb6b52a523910e4af3c8a03a82a0 100644
--- a/instrument/MANDI_Parameters_2015_08_01.xml
+++ b/instrument/MANDI_Parameters_2015_08_01.xml
@@ -8,4 +8,65 @@
</parameter>
</component-link>
+<!-- Need to fill in gaps for peak profile fitting -->
+<parameter name="fitConvolvedPeak" type="bool">
+ <value val="true"/>
+</parameter>
+
+<!-- Multiplier for profile fitting for BVG polar angle -->
+<parameter name="sigX0Scale">
+ <value val="2." />
+</parameter>
+
+<!-- Multiplier for profile fitting for BVG azimuthal angle -->
+<parameter name="sigY0Scale">
+ <value val="2." />
+</parameter>
+
+<!-- Number of rows between detector gaps for profile fitting -->
+<parameter name="numDetRows" type="int">
+ <value val="255" />
+</parameter>
+
+<!-- Number of cols between detector gaps for profile fitting -->
+<parameter name="numDetCols" type="int">
+ <value val="255" />
+</parameter>
+
+<!-- Number of polar bins for BVG histogramming when profile fitting -->
+<parameter name="numBinsTheta" type="int">
+ <value val="50" />
+</parameter>
+
+<!-- Number of azimuthal bins for BVG histogramming when profile fitting -->
+<parameter name="numBinsPhi" type="int">
+ <value val="50" />
+</parameter>
+
+<!-- Fraction along (h,k,l) to use for profile fitting. 0.5 is the next peak. -->
+<parameter name="fracHKL">
+ <value val="0.4" />
+</parameter>
+
+<!-- Side length of each voxel for fitting in units of angstrom^-1 -->
+<parameter name="dQPixel">
+ <value val="0.003" />
+</parameter>
+
+<!-- Minimum spacing for profile fitting the TOF profile. Units of microseconds -->
+<parameter name="mindtBinWidth">
+ <value val="15" />
+</parameter>
+
+<!-- Maximum spacing for profile fitting the TOF profile. Units of microseconds -->
+<parameter name="maxdtBinWidth">
+ <value val="50" />
+</parameter>
+
+<!-- Size of peak mask for background calculation in units of dQPixel -->
+<parameter name="peakMaskSize" type="int">
+ <value val="5" />
+</parameter>
+
+
</parameter-file>
diff --git a/instrument/MANDI_Parameters_2016_02_01.xml b/instrument/MANDI_Parameters_2016_02_01.xml
index 5172ef8dff0345ed75af05507ca19f83af787b98..5c875a8014aefca8597163bd2113d3bd668a3d31 100644
--- a/instrument/MANDI_Parameters_2016_02_01.xml
+++ b/instrument/MANDI_Parameters_2016_02_01.xml
@@ -6,6 +6,66 @@
<parameter name="T0">
<value val="-1.185000"/>
</parameter>
-</component-link>
+<!-- Need to fill in gaps for peak profile fitting -->
+<parameter name="fitConvolvedPeak" type="bool">
+ <value val="true"/>
+</parameter>
+
+<!-- Multiplier for profile fitting for BVG polar angle -->
+<parameter name="sigX0Scale">
+ <value val="2." />
+</parameter>
+
+<!-- Multiplier for profile fitting for BVG azimuthal angle -->
+<parameter name="sigY0Scale">
+ <value val="2." />
+</parameter>
+
+<!-- Number of rows between detector gaps for profile fitting -->
+<parameter name="numDetRows" type="int">
+ <value val="255" />
+</parameter>
+
+<!-- Number of cols between detector gaps for profile fitting -->
+<parameter name="numDetCols" type="int">
+ <value val="255" />
+</parameter>
+
+<!-- Number of polar bins for BVG histogramming when profile fitting -->
+<parameter name="numBinsTheta" type="int">
+ <value val="50" />
+</parameter>
+
+<!-- Number of azimuthal bins for BVG histogramming when profile fitting -->
+<parameter name="numBinsPhi" type="int">
+ <value val="50" />
+</parameter>
+
+<!-- Fraction along (h,k,l) to use for profile fitting. 0.5 is the next peak. -->
+<parameter name="fracHKL">
+ <value val="0.4" />
+</parameter>
+
+<!-- Side length of each voxel for fitting in units of angstrom^-1 -->
+<parameter name="dQPixel">
+ <value val="0.003" />
+</parameter>
+
+<!-- Minimum spacing for profile fitting the TOF profile. Units of microseconds -->
+<parameter name="mindtBinWidth">
+ <value val="15" />
+</parameter>
+
+<!-- Maximum spacing for profile fitting the TOF profile. Units of microseconds -->
+<parameter name="maxdtBinWidth">
+ <value val="50" />
+</parameter>
+
+<!-- Size of peak mask for background calculation in units of dQPixel -->
+<parameter name="peakMaskSize" type="int">
+ <value val="5" />
+</parameter>
+
+</component-link>
</parameter-file>
diff --git a/instrument/TOPAZ_Parameters.xml b/instrument/TOPAZ_Parameters.xml
index d0f57f5ff91415d011e034685985e5a84ee23d93..d2956eb6df187922bd173b30c6e5b95099a2aa24 100644
--- a/instrument/TOPAZ_Parameters.xml
+++ b/instrument/TOPAZ_Parameters.xml
@@ -83,6 +83,79 @@ detScale={13:1.046504,14:1.259293,16:1.02449,17:1.18898,18:0.88014,19:0.98665,\
<parameter name="detScale49">
<value val="0.79420"/>
</parameter>
+
+<!-- Need to fill in gaps for peak profile fitting -->
+<parameter name="fitConvolvedPeak" type="bool">
+ <value val="false"/>
+</parameter>
+
+<!-- Multiplier for profile fitting for BVG polar angle -->
+<parameter name="sigX0Scale">
+ <value val="3." />
+</parameter>
+
+<!-- Multiplier for profile fitting for BVG azimuthal angle -->
+<parameter name="sigY0Scale">
+ <value val="3." />
+</parameter>
+
+<!-- Number of rows between detector gaps for profile fitting -->
+<parameter name="numDetRows" type="int">
+ <value val="255" />
+</parameter>
+
+<!-- Number of cols between detector gaps for profile fitting -->
+<parameter name="numDetCols" type="int">
+ <value val="255" />
+</parameter>
+
+<!-- Number of polar bins for BVG histogramming when profile fitting -->
+<parameter name="numBinsTheta" type="int">
+ <value val="50" />
+</parameter>
+
+<!-- Number of azimuthal bins for BVG histogramming when profile fitting -->
+<parameter name="numBinsPhi" type="int">
+ <value val="50" />
+</parameter>
+
+<!-- Fraction along (h,k,l) to use for profile fitting. 0.5 is the next peak. -->
+<parameter name="fracHKL">
+ <value val="0.4" />
+</parameter>
+
+<!-- Side length of each voxel for fitting in units of angstrom^-1 -->
+<parameter name="dQPixel">
+ <value val="0.01" />
+</parameter>
+
+<!-- Minimum spacing for profile fitting the TOF profile. Units of microseconds -->
+<parameter name="mindtBinWidth">
+ <value val="2" />
+</parameter>
+
+<!-- Maximum spacing for profile fitting the TOF profile. Units of microseconds -->
+<parameter name="maxdtBinWidth">
+ <value val="15" />
+</parameter>
+
+<!-- Size of peak mask for background calculation in units of dQPixel -->
+<parameter name="peakMaskSize" type="int">
+ <value val="15" />
+</parameter>
+
+<!-- Constraints for ICC fitting. Valid names are iccA, iccB, iccR, iccT0, iccScale0
+ iccHatWidth and iccKConv. Inputs are strings with values separated by
+ spaces which are prased by the IntegratePeaksProfileFitting algorithm.
+ If two values are given they are treated as the lower and upper bounds. If
+ three are given they are the lower bound, upper bound, and initial guess.-->
+<parameter name="iccB" type="string">
+ <value val="0.001 0.3 0.005" />
+</parameter>
+
+<parameter name="iccKConv" type="string">
+ <value val="10.0 800.0 100.0" />
+</parameter>
</component-link>
</parameter-file>
diff --git a/scripts/SCD_Reduction/BVGFitTools.py b/scripts/SCD_Reduction/BVGFitTools.py
index 3a2aa938e19f7025e0dfb85720f9b37e16ba0d17..7fbf60c1295ecec73dbf5e3b70f4da661a252a01 100644
--- a/scripts/SCD_Reduction/BVGFitTools.py
+++ b/scripts/SCD_Reduction/BVGFitTools.py
@@ -11,11 +11,11 @@ import BivariateGaussian as BivariateGaussian
plt.ion()
-def get3DPeak(peak, box, padeCoefficients, qMask, nTheta=150, nPhi=150, fracBoxToHistogram=1.0,
+def get3DPeak(peak, peaks_ws, box, padeCoefficients, qMask, nTheta=150, nPhi=150, fracBoxToHistogram=1.0,
plotResults=False, zBG=1.96, bgPolyOrder=1, fICCParams=None, oldICCFit=None,
strongPeakParams=None, forceCutoff=250, edgeCutoff=15,
neigh_length_m=3, q_frame='sample', dtSpread=0.03, pplmin_frac=0.8, pplmax_frac=1.5, mindtBinWidth=1,
- figureNumber=2):
+ maxdtBinWidth=50, figureNumber=2, peakMaskSize=5, iccFitDict=None):
n_events = box.getNumEventsArray()
if q_frame == 'lab':
@@ -30,13 +30,14 @@ def get3DPeak(peak, box, padeCoefficients, qMask, nTheta=150, nPhi=150, fracBoxT
goodIDX, pp_lambda = ICCFT.getBGRemovedIndices(
n_events, peak=peak, box=box, qMask=qMask, calc_pp_lambda=True, padeCoefficients=padeCoefficients,
neigh_length_m=neigh_length_m, pp_lambda=None, pplmin_frac=pplmin_frac,
- pplmax_frac=pplmax_frac, mindtBinWidth=mindtBinWidth)
-
+ pplmax_frac=pplmax_frac, mindtBinWidth=mindtBinWidth, maxdtBinWidth=maxdtBinWidth,
+ peakMaskSize=peakMaskSize, iccFitDict=iccFitDict)
YTOF, fICC, x_lims = fitTOFCoordinate(
box, peak, padeCoefficients, dtSpread=dtSpread, qMask=qMask, bgPolyOrder=bgPolyOrder, zBG=zBG,
plotResults=plotResults, pp_lambda=pp_lambda, neigh_length_m=neigh_length_m, pplmin_frac=pplmin_frac,
- pplmax_frac=pplmax_frac, mindtBinWidth=mindtBinWidth)
-
+ pplmax_frac=pplmax_frac, mindtBinWidth=mindtBinWidth, maxdtBinWidth=maxdtBinWidth,
+ peakMaskSize=peakMaskSize, iccFitDict=iccFitDict)
+ chiSqTOF = mtd['fit_Parameters'].column(1)[-1]
else: # we already did I-C profile, so we'll just read the parameters
pp_lambda = fICCParams[-1]
fICC = ICC.IkedaCarpenterConvoluted()
@@ -49,7 +50,9 @@ def get3DPeak(peak, box, padeCoefficients, qMask, nTheta=150, nPhi=150, fracBoxT
fICC['HatWidth'] = fICCParams[10]
fICC['KConv'] = fICCParams[11]
goodIDX, _ = ICCFT.getBGRemovedIndices(
- n_events, pp_lambda=pp_lambda, qMask=qMask)
+ n_events, pp_lambda=pp_lambda, qMask=qMask, peakMaskSize=peakMaskSize,
+ iccFitDict=iccFitDict)
+ chiSqTOF = fICCParams[4] #Last entry
# Get the 3D TOF component, YTOF
if oldICCFit is not None:
@@ -71,8 +74,34 @@ def get3DPeak(peak, box, padeCoefficients, qMask, nTheta=150, nPhi=150, fracBoxT
goodIDX *= qMask
X = boxToTOFThetaPhi(box, peak)
dEdge = edgeCutoff
+
+ # This section defines detector size to determine if a peak is too
+ # close to the edge. Order is [NROWS, NCOLS].
+ try:
+ numDetRows = peaks_ws.getInstrument().getIntParameter("numDetRows")[0]
+ numDetCols = peaks_ws.getInstrument().getIntParameter("numDetCols")[0]
+ nPixels = [numDetRows, numDetCols]
+ except:
+ UserWarning('Detector size not found in instrument parameters file. Assuming a 255*255 detector!')
+ nPixels = [255,255]
+
useForceParams = peak.getIntensity() < forceCutoff or peak.getRow() <= dEdge or peak.getRow(
- ) >= 255 - dEdge or peak.getCol() <= dEdge or peak.getCol() >= 255 - dEdge
+ ) >= nPixels[0] - dEdge or peak.getCol() <= dEdge or peak.getCol() >= nPixels[1] - dEdge
+
+ #Here we retrieve some instrument specific parameters
+ try:
+ doPeakConvolution = peaks_ws.getInstrument().getBoolParameter("fitConvolvedPeak")[0]
+ except:
+ doPeakConvolution = False
+ try:
+ sigX0Scale = peaks_ws.getInstrument().getNumberParameter("sigX0Scale")[0]
+ except:
+ sigX0Scale = 1.0
+ try:
+ sigY0Scale = peaks_ws.getInstrument().getNumberParameter("sigY0Scale")[0]
+ except:
+ sigY0Scale = 1.0
+
if strongPeakParams is not None and useForceParams: # We will force parameters on this fit
ph = np.arctan2(q0[1], q0[0])
th = np.arctan2(q0[2], np.hypot(q0[0], q0[1]))
@@ -85,10 +114,12 @@ def get3DPeak(peak, box, padeCoefficients, qMask, nTheta=150, nPhi=150, fracBoxT
phthPeak[0],
phthPeak[1]))
params, h, t, p = doBVGFit(box, nTheta=nTheta, nPhi=nPhi, fracBoxToHistogram=fracBoxToHistogram,
- goodIDX=goodIDX, forceParams=strongPeakParams[nnIDX])
+ goodIDX=goodIDX, forceParams=strongPeakParams[nnIDX],
+ doPeakConvolution=doPeakConvolution, sigX0Scale=sigX0Scale, sigY0Scale=sigY0Scale)
else: # Just do the fit - no nearest neighbor assumptions
params, h, t, p = doBVGFit(
- box, nTheta=nTheta, nPhi=nPhi, fracBoxToHistogram=fracBoxToHistogram, goodIDX=goodIDX)
+ box, nTheta=nTheta, nPhi=nPhi, fracBoxToHistogram=fracBoxToHistogram, goodIDX=goodIDX,
+ doPeakConvolution=doPeakConvolution, sigX0Scale=sigX0Scale, sigY0Scale=sigY0Scale)
if plotResults:
compareBVGFitData(
@@ -113,7 +144,10 @@ def get3DPeak(peak, box, padeCoefficients, qMask, nTheta=150, nPhi=150, fracBoxT
YBVG = bvg(1.0, mu, sigma, XTHETA, XPHI, 0)
# Do scaling to the data
- Y, redChiSq, scaleFactor = fitScaling(n_events, box, YTOF, YBVG)
+ if doPeakConvolution: #This means peaks will have gaps, so we only use good data to scale
+ Y, redChiSq, scaleFactor = fitScaling(n_events, box, YTOF, YBVG, goodIDX=goodIDX)
+ else:
+ Y, redChiSq, scaleFactor = fitScaling(n_events, box, YTOF, YBVG)
YBVG2 = bvg(1.0, mu, sigma, XTHETA, XPHI, 0)
YTOF2 = getYTOF(fICC, XTOF, x_lims)
Y2 = YTOF2 * YBVG2
@@ -139,6 +173,7 @@ def get3DPeak(peak, box, padeCoefficients, qMask, nTheta=150, nPhi=150, fracBoxT
retParams['bgBVG'] = bgBVG
retParams['scale3d'] = scaleFactor
retParams['chiSq3d'] = redChiSq
+ retParams['chiSq'] = chiSqTOF
retParams['dQ'] = np.linalg.norm(newCenter - q0)
retParams['newQ'] = newCenter
@@ -176,8 +211,7 @@ def fitScaling(n_events, box, YTOF, YBVG, goodIDX=None, neigh_length_m=3):
goodIDX[max(fitMaxIDX[0] - dP, 0):min(fitMaxIDX[0] + dP, goodIDX.shape[0]),
max(fitMaxIDX[1] - dP, 0):min(fitMaxIDX[1] + dP, goodIDX.shape[1]),
max(fitMaxIDX[2] - dP, 0):min(fitMaxIDX[2] + dP, goodIDX.shape[2])] = True
- goodIDX = np.logical_and(goodIDX, conv_n_events > 0)
-
+ goodIDX = np.logical_and(goodIDX, conv_n_events > 0)
# A1 = slope, A0 = offset
scaleLinear = Polynomial(n=1)
scaleLinear.constrain("A1>0")
@@ -217,7 +251,8 @@ def getXTOF(box, peak):
def fitTOFCoordinate(box, peak, padeCoefficients, dtSpread=0.03, minFracPixels=0.01,
neigh_length_m=3, zBG=1.96, bgPolyOrder=1, qMask=None, plotResults=False,
- fracStop=0.01, pp_lambda=None, pplmin_frac=0.8, pplmax_frac=1.5, mindtBinWidth=1):
+ fracStop=0.01, pp_lambda=None, pplmin_frac=0.8, pplmax_frac=1.5, mindtBinWidth=1,
+ maxdtBinWidth=50, peakMaskSize=5, iccFitDict=None):
# Get info from the peak
tof = peak.getTOF() # in us
@@ -235,10 +270,13 @@ def fitTOFCoordinate(box, peak, padeCoefficients, dtSpread=0.03, minFracPixels=0
dtSpread=dtSpread, minFracPixels=minFracPixels,
neigh_length_m=neigh_length_m, zBG=zBG, pp_lambda=pp_lambda,
pplmin_frac=pplmin_frac, pplmax_frac=pplmax_frac,
- mindtBinWidth=mindtBinWidth)
+ mindtBinWidth=mindtBinWidth, maxdtBinWidth=maxdtBinWidth,
+ peakMaskSize=peakMaskSize,
+ iccFitDict=iccFitDict)
fitResults, fICC = ICCFT.doICCFit(tofWS, energy, flightPath,
- padeCoefficients, fitOrder=bgPolyOrder, constraintScheme=1)
+ padeCoefficients, fitOrder=bgPolyOrder, constraintScheme=1,
+ iccFitDict=iccFitDict)
for i, param in enumerate(['A', 'B', 'R', 'T0', 'Scale', 'HatWidth', 'KConv']):
fICC[param] = mtd['fit_Parameters'].row(i)['Value']
@@ -391,7 +429,8 @@ def compareBVGFitData(box, params, nTheta=200, nPhi=200, figNumber=2, fracBoxToH
def doBVGFit(box, nTheta=200, nPhi=200, zBG=1.96, fracBoxToHistogram=1.0, goodIDX=None,
- forceParams=None, forceTolerance=0.1, dth=10, dph=10):
+ forceParams=None, forceTolerance=0.1, dth=10, dph=10,
+ doPeakConvolution=False, sigX0Scale=1., sigY0Scale=1.):
"""
doBVGFit takes a binned MDbox and returns the fit of the peak shape along the non-TOF direction. This is done in one of two ways:
1) Standard least squares fit of the 2D histogram.
@@ -408,6 +447,8 @@ def doBVGFit(box, nTheta=200, nPhi=200, zBG=1.96, fracBoxToHistogram=1.0, goodID
forceParams: set of parameters to force. These are the same format as a row in strongPeaksParams
forceTolerance: the factor we allow sigX, sigY, sigP to change when forcing peaks. Not used if forceParams is None.
dth, dph: The peak center may move by (dth, dph) from predicted position (in units of histogram pixels).
+ doPeakConvolution: boolean stating whether we should fit a convolved (smoothed) peak. This is useful for filling in
+ gaps for 3He detector tube packs.
"""
h, thBins, phBins = getAngularHistogram(
@@ -454,6 +495,17 @@ def doBVGFit(box, nTheta=200, nPhi=200, zBG=1.96, fracBoxToHistogram=1.0, goodID
boundsDict['SigP'] = [-1., 1.]
boundsDict['Bg'] = [0, np.inf]
+ # Here we can make instrument-specific changes to our initial guesses and boundaries
+ sigX0 = sigX0*sigX0Scale
+ sigY0 = sigY0*sigY0Scale
+
+ if doPeakConvolution:
+ neigh_length_m = 5
+ convBox = 1.0*np.ones([neigh_length_m, neigh_length_m]) / neigh_length_m**2
+ conv_h = convolve(h, convBox)
+ H[:,:,0] = conv_h
+ H[:,:,1] = conv_h
+
# Set our initial guess
m = BivariateGaussian.BivariateGaussian()
m.init()
@@ -468,13 +520,15 @@ def doBVGFit(box, nTheta=200, nPhi=200, zBG=1.96, fracBoxToHistogram=1.0, goodID
m.setAttributeValue('nX', h.shape[0])
m.setAttributeValue('nY', h.shape[1])
m.setConstraints(boundsDict)
+
# Do the fit
+ #bvgWS = CreateWorkspace(OutputWorkspace='bvgWS', DataX=pos.ravel(
+ #), DataY=H.ravel(), DataE=np.sqrt(H.ravel()))
bvgWS = CreateWorkspace(OutputWorkspace='bvgWS', DataX=pos.ravel(
), DataY=H.ravel(), DataE=np.sqrt(H.ravel()))
fitResults = Fit(Function=m, InputWorkspace='bvgWS', Output='bvgfit',
Minimizer='Levenberg-MarquardtMD')
-
elif forceParams is not None:
p0 = np.zeros(7)
p0[0] = np.max(h)
@@ -508,6 +562,16 @@ def doBVGFit(box, nTheta=200, nPhi=200, zBG=1.96, fracBoxToHistogram=1.0, goodID
boundsDict['SigX'] = [bounds[0][3], bounds[1][3]]
boundsDict['SigY'] = [bounds[0][4], bounds[1][4]]
boundsDict['SigP'] = [bounds[0][5], bounds[1][5]]
+
+ # Here we can make instrument-specific changes to our initial guesses and boundaries
+
+ if doPeakConvolution:
+ neigh_length_m = 5
+ convBox = 1.0*np.ones([neigh_length_m, neigh_length_m]) / neigh_length_m**2
+ conv_h = convolve(h, convBox)
+ H[:,:,0] = conv_h
+ H[:,:,1] = conv_h
+
# Set our initial guess
m = BivariateGaussian.BivariateGaussian()
m.init()
diff --git a/scripts/SCD_Reduction/ICCFitTools.py b/scripts/SCD_Reduction/ICCFitTools.py
index e450536bfa6ddfebe0b90009dfc23122327869b1..a7186882ac598707262b638480986fbeec87ad34 100644
--- a/scripts/SCD_Reduction/ICCFitTools.py
+++ b/scripts/SCD_Reduction/ICCFitTools.py
@@ -13,6 +13,21 @@ from scipy.ndimage.filters import convolve
plt.ion()
+def parseConstraints(peaks_ws):
+ """
+ returns a dictionary containing parameters for ICC fitting. Parameters
+ are derived from instrument parameters files (see MANDI_Parameters.xml
+ for an example).
+ """
+ possibleKeys = ['iccA', 'iccB', 'iccR', 'iccT0', 'iccScale0', 'iccHatWidth', 'iccKConv']
+ d = {}
+ for paramName in possibleKeys:
+ if peaks_ws.getInstrument().hasParameter(paramName):
+ vals = np.array(peaks_ws.getInstrument().getStringParameter(paramName)[0].split(),dtype=float)
+ d[paramName] = vals
+ return d
+
+
def scatFun(x, A, bg):
"""
scatFun: returns A/x+bg. Used for background estimation.
@@ -92,7 +107,8 @@ def getQXQYQZ(box):
def getQuickTOFWS(box, peak, padeCoefficients, goodIDX=None, dtSpread=0.03, qMask=None,
- pp_lambda=None, minppl_frac=0.8, maxppl_frac=1.5, mindtBinWidth=1, constraintScheme=1):
+ pp_lambda=None, minppl_frac=0.8, maxppl_frac=1.5, mindtBinWidth=1, maxdtBinWidth=50,
+ constraintScheme=1, peakMaskSize=5, iccFitDict=None):
"""
getQuickTOFWS - generates a quick-and-dirty TOFWS. Useful for determining the background.
Input:
@@ -106,7 +122,9 @@ def getQuickTOFWS(box, peak, padeCoefficients, goodIDX=None, dtSpread=0.03, qMas
pp_lambda - nominal background level. Will be calculated if set to None.
minppl_frac, maxppl_frac: fraction of the predicted pp_lambda to try if calculating pp_lambda
mindtBinWidth - the minimum binwidth (in us) that we will allow when histogramming.
+ maxdtBinWidth - the maximum binwidth (in us) that we will allow when histogramming.
constraintScheme - which constraint scheme we use. Typically set to 1
+ iccFitDict - a dictionary containing ICC fit constraints and possibly initial guesses
Output:
chiSq - reduced chiSquared from fitting the TOF profile
h - list of [Y, X], with Y and X being numpy arrays of the Y and X of the tof profile
@@ -128,9 +146,11 @@ def getQuickTOFWS(box, peak, padeCoefficients, goodIDX=None, dtSpread=0.03, qMas
tofWS, ppl = getTOFWS(box, flightPath, scatteringHalfAngle, tof, peak, qMask, dtSpread=dtSpread,
minFracPixels=0.01, neigh_length_m=3, zBG=1.96, pp_lambda=pp_lambda,
calc_pp_lambda=calc_pp_lambda, pplmin_frac=minppl_frac, pplmax_frac=minppl_frac,
- mindtBinWidth=mindtBinWidth)
+ mindtBinWidth=mindtBinWidth, maxdtBinWidth=maxdtBinWidth,
+ peakMaskSize=peakMaskSize, iccFitDict=iccFitDict)
fitResults, fICC = doICCFit(
- tofWS, energy, flightPath, padeCoefficients, fitOrder=1, constraintScheme=constraintScheme)
+ tofWS, energy, flightPath, padeCoefficients, fitOrder=1, constraintScheme=constraintScheme,
+ iccFitDict=iccFitDict)
h = [tofWS.readY(0), tofWS.readX(0)]
chiSq = fitResults.OutputChi2overDoF
@@ -193,7 +213,8 @@ def getPoissionGoodIDX(n_events, zBG=1.96, neigh_length_m=3):
def getOptimizedGoodIDX(n_events, padeCoefficients, zBG=1.96, neigh_length_m=3, qMask=None,
peak=None, box=None, pp_lambda=None, peakNumber=-1, minppl_frac=0.8,
- maxppl_frac=1.5, mindtBinWidth=1, constraintScheme=1):
+ maxppl_frac=1.5, mindtBinWidth=1, maxdtBinWidth=50,
+ constraintScheme=1, peakMaskSize=5, iccFitDict=None):
"""
getOptimizedGoodIDX - returns a numpy arrays which is true if the voxel contains events at
the zBG z level (1.96=95%CI). Rather than using Poission statistics, this function
@@ -211,9 +232,11 @@ def getOptimizedGoodIDX(n_events, padeCoefficients, zBG=1.96, neigh_length_m=3,
pp_lambda - Currently unused. Leave as None. TODO: remove this.
peakNumber - currently unused. TODO: Remove this.
minppl_frac, maxppl_frac; range around predicted pp_lambda to check.
- mindtBinWidth - the small dt (in us) allowed for constructing the TOF profile.
+ mindtBinWidth - the smallest dt (in us) allowed for constructing the TOF profile.
+ mindtBinWidth - the largest dt (in us) allowed for constructing the TOF profile.
constraintScheme - sets the constraints for TOF profile fitting. Leave as 1 if you're
not sure how to modify this.
+ iccFitDict - a dictionary containing ICC fit constraints and possibly initial guesses
Output:
goodIDX: a numpy arrays the same size as n_events that is True of False for if it contains
@@ -239,7 +262,8 @@ def getOptimizedGoodIDX(n_events, padeCoefficients, zBG=1.96, neigh_length_m=3,
cX = nX//2
cY = nY//2
cZ = nZ//2
- dP = 5
+ dP = peakMaskSize
+
peakMask = qMask.copy()
peakMask[cX-dP:cX+dP, cY-dP:cY+dP, cZ-dP:cZ+dP] = 0
neigh_length_m=3
@@ -274,9 +298,10 @@ def getOptimizedGoodIDX(n_events, padeCoefficients, zBG=1.96, neigh_length_m=3,
try:
chiSq, h, intens, sigma = getQuickTOFWS(box, peak, padeCoefficients, goodIDX=goodIDX, qMask=qMask, pp_lambda=pp_lambda,
minppl_frac=minppl_frac, maxppl_frac=maxppl_frac, mindtBinWidth=mindtBinWidth,
- constraintScheme=constraintScheme)
+ maxdtBinWidth=maxdtBinWidth, constraintScheme=constraintScheme,
+ peakMaskSize=peakMaskSize, iccFitDict=iccFitDict)
except:
- # raise
+ #raise
break
chiSqList[i] = chiSq
ISIGList[i] = intens/sigma
@@ -288,6 +313,7 @@ def getOptimizedGoodIDX(n_events, padeCoefficients, zBG=1.96, neigh_length_m=3,
except RuntimeError:
# This is caused by there being fewer datapoints remaining than parameters. For now, we just hope
# we found a satisfactory answer.
+ raise
break
except KeyboardInterrupt:
sys.exit()
@@ -296,7 +322,9 @@ def getOptimizedGoodIDX(n_events, padeCoefficients, zBG=1.96, neigh_length_m=3,
goodIDX, _ = getBGRemovedIndices(n_events, pp_lambda=pp_lambda)
chiSq, h, intens, sigma = getQuickTOFWS(box, peak, padeCoefficients, goodIDX=goodIDX, qMask=qMask,
pp_lambda=pp_lambda, minppl_frac=minppl_frac, maxppl_frac=maxppl_frac,
- mindtBinWidth=mindtBinWidth)
+ mindtBinWidth=mindtBinWidth, maxdtBinWidth=maxdtBinWidth,
+ peakMaskSize=peakMaskSize,
+ iccFitDict=iccFitDict)
if qMask is not None:
return goodIDX*qMask, pp_lambda
return goodIDX, pp_lambda
@@ -304,7 +332,8 @@ def getOptimizedGoodIDX(n_events, padeCoefficients, zBG=1.96, neigh_length_m=3,
def getBGRemovedIndices(n_events, zBG=1.96, calc_pp_lambda=False, neigh_length_m=3, qMask=None,
peak=None, box=None, pp_lambda=None, peakNumber=-1, padeCoefficients=None,
- pplmin_frac=0.8, pplmax_frac=1.5, mindtBinWidth=1, constraintScheme=1):
+ pplmin_frac=0.8, pplmax_frac=1.5, mindtBinWidth=1, maxdtBinWidth=50,
+ constraintScheme=1, peakMaskSize=5, iccFitDict=None):
"""
getBGRemovedIndices - A wrapper for getOptimizedGoodIDX
Input:
@@ -319,9 +348,11 @@ def getBGRemovedIndices(n_events, zBG=1.96, calc_pp_lambda=False, neigh_length_m
pp_lambda - Currently unused. Leave as None. TODO: remove this.
peakNumber - currently unused. TODO: Remove this.
minppl_frac, maxppl_frac; range around predicted pp_lambda to check.
- mindtBinWidth - the small dt (in us) allowed for constructing the TOF profile.
+ mindtBinWidth - the smallest dt (in us) allowed for constructing the TOF profile.
+ maxdtBinWidth - the largest dt (in us) allowed for constructing the TOF profile.
constraintScheme - sets the constraints for TOF profile fitting. Leave as 1 if you're
not sure how to modify this.
+ iccFitDict - a dictionary containing ICC fit constraints and possibly initial guesses
Output:
goodIDX: a numpy arrays the same size as n_events that is True of False for if it contains
@@ -337,8 +368,6 @@ def getBGRemovedIndices(n_events, zBG=1.96, calc_pp_lambda=False, neigh_length_m
sys.exit(
'Error in ICCFT:getBGRemovedIndices: calc_pp_lambda is True, but no moderator coefficients are provided.')
- # TODO: this result should be multiplied by qMask if qMask is not None - but I need to check that that change won't affect
- # other workflows.
if pp_lambda is not None:
# Set up some things to only consider good pixels
hasEventsIDX = n_events > 0
@@ -361,7 +390,9 @@ def getBGRemovedIndices(n_events, zBG=1.96, calc_pp_lambda=False, neigh_length_m
return getOptimizedGoodIDX(n_events, padeCoefficients, zBG=1.96, neigh_length_m=neigh_length_m,
minppl_frac=pplmin_frac, maxppl_frac=pplmax_frac, qMask=qMask, peak=peak,
box=box, pp_lambda=pp_lambda, peakNumber=peakNumber,
- mindtBinWidth=mindtBinWidth, constraintScheme=constraintScheme)
+ mindtBinWidth=mindtBinWidth, maxdtBinWidth=maxdtBinWidth,
+ constraintScheme=constraintScheme,
+ peakMaskSize=peakMaskSize, iccFitDict=iccFitDict)
except KeyboardInterrupt:
sys.exit()
except:
@@ -521,8 +552,8 @@ def get_pp_lambda(n_events, hasEventsIDX):
def getTOFWS(box, flightPath, scatteringHalfAngle, tofPeak, peak, qMask, zBG=-1.0, dtSpread=0.02,
minFracPixels=0.005, workspaceNumber=None, neigh_length_m=0, pp_lambda=None, calc_pp_lambda=False,
- padeCoefficients=None, pplmin_frac=0.8, pplmax_frac=1.5,
- mindtBinWidth=1, constraintScheme=1):
+ padeCoefficients=None, pplmin_frac=0.8, pplmax_frac=1.5, peakMaskSize=5,
+ mindtBinWidth=1, maxdtBinWidth=50, constraintScheme=1, iccFitDict=None):
"""
Builds a TOF profile from the data in box which is nominally centered around a peak.
Input:
@@ -545,8 +576,10 @@ def getTOFWS(box, flightPath, scatteringHalfAngle, tofPeak, peak, qMask, zBG=-1.
want to, you can feed the value in as pp_lambda (calculated elsewhere).
minppl_frac, maxppl_frac; range around predicted pp_lambda to check.
mindtBinWidth - the small dt (in us) allowed for constructing the TOF profile.
+ maxdtBinWidth - the largest dt (in us) allowed for constructing the TOF profile.
constraintScheme - sets the constraints for TOF profile fitting. Leave as 1 if you're
not sure how to modify this.
+ iccFitDict - a dictionary containing ICC fit constraints and possibly initial guesses
Output:
tofWS - a mantid containing the TOF profile. X-axis is TOF (units: us) and
@@ -563,8 +596,9 @@ def getTOFWS(box, flightPath, scatteringHalfAngle, tofPeak, peak, qMask, zBG=-1.
goodIDX, pp_lambda = getBGRemovedIndices(n_events, box=box, qMask=qMask, peak=peak, pp_lambda=pp_lambda,
calc_pp_lambda=calc_pp_lambda, padeCoefficients=padeCoefficients,
pplmin_frac=pplmin_frac, pplmax_frac=pplmax_frac,
- mindtBinWidth=mindtBinWidth, constraintScheme=constraintScheme)
- # TODO bad naming, but a lot of the naming in this function assumes it
+ mindtBinWidth=mindtBinWidth, maxdtBinWidth=maxdtBinWidth,
+ constraintScheme=constraintScheme,
+ peakMaskSize=peakMaskSize, iccFitDict=iccFitDict)
hasEventsIDX = np.logical_and(goodIDX, qMask)
boxMeanIDX = np.where(hasEventsIDX)
else: # don't do background removal - just consider one pixel at a time
@@ -608,7 +642,7 @@ def getTOFWS(box, flightPath, scatteringHalfAngle, tofPeak, peak, qMask, zBG=-1.
[qx[qx.shape[0]//2 + 1], qy[qy.shape[0]//2+1], qz[qz.shape[0]//2+1]])
dtBinWidth = np.abs(tD-tC)
dtBinWidth = max(mindtBinWidth, dtBinWidth)
- dtBinWidth = min(50, dtBinWidth)
+ dtBinWidth = min(maxdtBinWidth, dtBinWidth)
tBins = np.arange(tMin, tMax, dtBinWidth)
weightList = n_events[hasEventsIDX] # - pp_lambda
h = np.histogram(tList, tBins, weights=weightList)
@@ -792,7 +826,8 @@ def getBoxFracHKL(peak, peaks_ws, MDdata, UBMatrix, peakNumber, dQ, dQPixel=0.00
return Box
-def doICCFit(tofWS, energy, flightPath, padeCoefficients, constraintScheme=None, outputWSName='fit', fitOrder=1):
+def doICCFit(tofWS, energy, flightPath, padeCoefficients, constraintScheme=None, outputWSName='fit', fitOrder=1,
+ iccFitDict=None):
"""
doICCFit - Carries out the actual least squares fit for the TOF workspace.
Intput:
@@ -811,6 +846,7 @@ def doICCFit(tofWS, energy, flightPath, padeCoefficients, constraintScheme=None,
outputWSName - the base name for output workspaces. Leave as 'fit' unless you are
doing multiple fits.
fitOrder - the background polynomial order
+ iccFitDict - a dictionary containing ICC fit constraints and possibly initial guesses
Returns:
fitResults - the output from Mantid's Fit() routine
fICC - an IkedaCarpenterConvoluted function with parameters set to the fit values.
@@ -833,12 +869,30 @@ def doICCFit(tofWS, energy, flightPath, padeCoefficients, constraintScheme=None,
fICC.setParameter(4, x0[4])
#fICC.setPenalizedConstraints(A0=[0.01, 1.0], B0=[0.005, 1.5], R0=[0.01, 1.0], T00=[0,1.0e10], KConv0=[10,500],penalty=1.0e20)
if constraintScheme == 1:
+ # Set these bounds as defaults - they can be changed for each instrument
+ # They can be changed by setting parameters in the INSTRUMENT_Parameters.xml file.
+ A0 = [0.5*x0[0], 1.5*x0[0]]
+ B0 = [0.5*x0[1], 1.5*x0[1]]
+ R0 = [0.5*x0[2], 1.5*x0[2]]
+ T00 = [0.,1.e10]
+ HatWidth0 = [0., 5.]
+ Scale0 = [0., np.inf]
+ KConv0 = [100, 140]
+
+ # Now we see what instrument specific parameters we have
+ if iccFitDict is not None:
+ possibleKeys = ['iccA', 'iccB', 'iccR', 'iccT0', 'iccScale0', 'iccHatWidth', 'iccKConv']
+ for keyIDX, (key, bounds) in enumerate(zip(possibleKeys, [A0, B0, R0, T00, Scale0, HatWidth0, KConv0])):
+ if key in iccFitDict:
+ bounds[0] = iccFitDict[key][0]
+ bounds[1] = iccFitDict[key][1]
+ if len(iccFitDict[key] == 3):
+ x0[keyIDX] = iccFitDict[key][2]
+ fICC.setParameter(keyIDX, x0[keyIDX])
try:
- fICC.setPenalizedConstraints(A0=[0.5*x0[0], 1.5*x0[0]], B0=[0.5*x0[1], 1.5*x0[1]], R0=[
- 0.5*x0[2], 1.5*x0[2]], T00=[0., 1.e10], KConv0=[100, 140], penalty=1.0e10)
+ fICC.setPenalizedConstraints(A0=A0, B0=B0, R0=R0, T00=T00, KConv0=KConv0, penalty=1.0e10)
except:
- fICC.setPenalizedConstraints(A0=[0.5*x0[0], 1.5*x0[0]], B0=[0.5*x0[1], 1.5*x0[1]], R0=[
- 0.5*x0[2], 1.5*x0[2]], T00=[0., 1.e10], KConv0=[100, 140], penalty=None)
+ fICC.setPenalizedConstraints(A0=A0, B0=B0, R0=R0, T00=T00, KConv0=KConv0, penalty=None)
if constraintScheme == 2:
try:
fICC.setPenalizedConstraints(A0=[0.0001, 1.0], B0=[0.005, 1.5], R0=[0.00, 1.], Scale0=[
@@ -862,7 +916,8 @@ def integrateSample(run, MDdata, peaks_ws, paramList, UBMatrix, dQ, qMask, padeC
figsFormat=None, dtSpread=0.02, fracHKL=0.5, minFracPixels=0.0000, fracStop=0.01,
dQPixel=0.005, p=None, neigh_length_m=0, zBG=-1.0, bgPolyOrder=1,
doIterativeBackgroundFitting=False, q_frame='sample',
- progressFile=None, minpplfrac=0.8, maxpplfrac=1.5, mindtBinWidth=1, keepFitDict=False, constraintScheme=1):
+ progressFile=None, minpplfrac=0.8, maxpplfrac=1.5, mindtBinWidth=1, maxdtBinWidth=50,
+ keepFitDict=False, constraintScheme=1, peakMaskSize=5, iccFitDict=None):
"""
integrateSample contains the loop that integrates over all of the peaks in a run and saves the results. Importantly, it also handles
errors (mostly by passing and recording special values for failed fits.)
@@ -894,9 +949,11 @@ def integrateSample(run, MDdata, peaks_ws, paramList, UBMatrix, dQ, qMask, padeC
minpplfrac, maxpplfrac - the range of pp_lambdas to check around the predicted pp_lambda as a fraction
of pp_lambda
mindtBinWidth - the smallest dt bin width (in us) allowed for TOF profile construction
+ mindtBinWidth - the largest dt bin width (in us) allowed for TOF profile construction
keepFitDict= bool; if True then each fit will be saved in a dictionary and returned. For large peak sets,
this can take a lot of memory.
constraintScheme - which constraint scheme we will use - leave as 1 if you're not sure what this does.
+ iccFitDict - a dictionary containing ICC fit constraints and possibly initial guesses
Returns:
peaks_ws - the peaks_ws with updated I, sig(I)
paramList - a list of fit parameters for each peak. Parameters are in the order:
@@ -934,13 +991,16 @@ def integrateSample(run, MDdata, peaks_ws, paramList, UBMatrix, dQ, qMask, padeC
goodIDX, pp_lambda = getBGRemovedIndices(n_events, peak=peak, box=Box, qMask=qMask,
calc_pp_lambda=True, padeCoefficients=padeCoefficients,
mindtBinWidth=mindtBinWidth,
+ maxdtBinWidth=maxdtBinWidth,
pplmin_frac=minpplfrac, pplmax_frac=maxpplfrac,
- constraintScheme=constraintScheme)
+ constraintScheme=constraintScheme,
+ peakMaskSize=peakMaskSize, iccFitDict=iccFitDict)
# --IN PRINCIPLE!!! WE CALCULATE THIS BEFORE GETTING HERE
tofWS = mtd['tofWS']
fitResults, fICC = doICCFit(
- tofWS, energy, flightPath, padeCoefficients, fitOrder=bgPolyOrder, constraintScheme=constraintScheme)
+ tofWS, energy, flightPath, padeCoefficients, fitOrder=bgPolyOrder, constraintScheme=constraintScheme,
+ iccFitDict=iccFitDict)
chiSq = fitResults.OutputChi2overDoF
r = mtd['fit_Workspace']