Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Sign in
Toggle navigation
Menu
Open sidebar
mantidproject
mantid
Commits
264344dc
Commit
264344dc
authored
Jun 05, 2020
by
Anthony Lim
Browse files
updated fitting range for PSI background subtraction - clang format
parent
fe5bd167
Changes
2
Hide whitespace changes
Inline
Side-by-side
Framework/Muon/src/PSIBackgroundSubtraction.cpp
View file @
264344dc
...
...
@@ -8,10 +8,10 @@
#include
"MantidAPI/FunctionFactory.h"
#include
"MantidAPI/IFunction.h"
#include
"MantidAPI/MatrixWorkspace.h"
#include
"MantidAPI/Run.h"
#include
"MantidAPI/WorkspaceFactory.h"
#include
"MantidDataObjects/WorkspaceSingleValue.h"
#include
"MantidKernel/BoundedValidator.h"
#include
"MantidAPI/Run.h"
#include
<math.h>
...
...
@@ -31,11 +31,9 @@ const std::string LAST_GOOD = "Last good spectra ";
std
::
pair
<
double
,
double
>
getRange
(
const
MatrixWorkspace_sptr
&
inputWorkspace
,
const
size_t
&
index
)
{
auto
firstGoodIndex
=
std
::
stoi
(
inputWorkspace
->
getLog
(
FIRST_GOOD
+
std
::
to_string
(
index
))
->
value
());
inputWorkspace
->
getLog
(
FIRST_GOOD
+
std
::
to_string
(
index
))
->
value
());
auto
lastGoodIndex
=
std
::
stoi
(
inputWorkspace
->
getLog
(
LAST_GOOD
+
std
::
to_string
(
index
))
->
value
());
inputWorkspace
->
getLog
(
LAST_GOOD
+
std
::
to_string
(
index
))
->
value
());
auto
midGoodIndex
=
int
(
floor
(((
double
(
lastGoodIndex
)
-
double
(
firstGoodIndex
))
/
2.
)));
...
...
@@ -83,10 +81,11 @@ std::map<std::string, std::string> PSIBackgroundSubtraction::validateInputs() {
const
Mantid
::
API
::
Run
&
run
=
inputWS
->
run
();
for
(
size_t
index
=
0
;
index
<
inputWS
->
getNumberHistograms
();
index
++
)
{
int
firstGood
=
0
;
int
lastGood
=
0
;
try
{
firstGood
=
std
::
stod
(
run
.
getProperty
(
FIRST_GOOD
+
std
::
to_string
(
index
))
->
value
());
int
firstGood
=
0
;
int
lastGood
=
0
;
try
{
firstGood
=
std
::
stod
(
run
.
getProperty
(
FIRST_GOOD
+
std
::
to_string
(
index
))
->
value
());
}
catch
(
Kernel
::
Exception
::
NotFoundError
)
{
errors
[
"InputWorkspace"
]
=
"Input Workspace should should contain first food data. "
;
...
...
@@ -101,109 +100,106 @@ std::map<std::string, std::string> PSIBackgroundSubtraction::validateInputs() {
if
(
lastGood
<=
firstGood
)
{
errors
[
"InputWorkspace"
]
+=
"
\n
Input Workspace should have last good data > first good data. "
;
}
}
if
(
firstGood
<=
0
)
{
errors
[
"InputWorkspace"
]
+=
"
\n
Input Workspace should have first good data > 0. "
;
}
errors
[
"InputWorkspace"
]
+=
"
\n
Input Workspace should have first good data > 0. "
;
}
if
(
lastGood
>
inputWS
->
readX
(
index
).
size
())
{
errors
[
"InputWorkspace"
]
+=
"
\n
Input Workspace should have last good data < number of bins. "
;
}
}
return
errors
;
}
return
errors
;
}
void
PSIBackgroundSubtraction
::
exec
()
{
MatrixWorkspace_sptr
inputWS
=
getProperty
(
"InputWorkspace"
);
// Caclulate and subtract background from inputWS
calculateBackgroundUsingFit
(
inputWS
);
}
/**
* Calculate the background of a PSI workspace by performing a fit, comprising
of a FlatBackground and ExpDecayMuon, on the second half of the PSI data.
* @param inputWorkspace :: Input PSI workspace which is modified inplace
*/
void
PSIBackgroundSubtraction
::
calculateBackgroundUsingFit
(
MatrixWorkspace_sptr
&
inputWorkspace
)
{
IAlgorithm_sptr
fit
=
setupFitAlgorithm
(
inputWorkspace
->
getName
());
auto
numberOfHistograms
=
inputWorkspace
->
getNumberHistograms
();
std
::
vector
<
double
>
backgroundValues
(
numberOfHistograms
);
for
(
size_t
index
=
0
;
index
<
numberOfHistograms
;
++
index
)
{
std
::
pair
<
double
,
double
>
range
=
getRange
(
inputWorkspace
,
index
);
auto
[
background
,
fitQuality
]
=
calculateBackgroundFromFit
(
fit
,
range
,
static_cast
<
int
>
(
index
));
// If fit quality is poor, do not subtract the background and instead log
// a warning
if
(
fitQuality
>
FIT_TOLERANCE
)
{
g_log
.
warning
()
<<
"Fit quality obtained in PSIBackgroundSubtraction is poor. "
<<
"Skipping background calculation for WorkspaceIndex: "
<<
static_cast
<
int
>
(
index
)
<<
"
\n
"
;
}
else
{
backgroundValues
[
index
]
=
background
;
}
void
PSIBackgroundSubtraction
::
exec
()
{
MatrixWorkspace_sptr
inputWS
=
getProperty
(
"InputWorkspace"
);
// Caclulate and subtract background from inputWS
calculateBackgroundUsingFit
(
inputWS
);
}
/**
* Calculate the background of a PSI workspace by performing a fit, comprising
of a FlatBackground and ExpDecayMuon, on the second half of the PSI data.
* @param inputWorkspace :: Input PSI workspace which is modified inplace
*/
void
PSIBackgroundSubtraction
::
calculateBackgroundUsingFit
(
MatrixWorkspace_sptr
&
inputWorkspace
)
{
IAlgorithm_sptr
fit
=
setupFitAlgorithm
(
inputWorkspace
->
getName
());
auto
numberOfHistograms
=
inputWorkspace
->
getNumberHistograms
();
std
::
vector
<
double
>
backgroundValues
(
numberOfHistograms
);
for
(
size_t
index
=
0
;
index
<
numberOfHistograms
;
++
index
)
{
std
::
pair
<
double
,
double
>
range
=
getRange
(
inputWorkspace
,
index
);
auto
[
background
,
fitQuality
]
=
calculateBackgroundFromFit
(
fit
,
range
,
static_cast
<
int
>
(
index
));
// If fit quality is poor, do not subtract the background and instead log
// a warning
if
(
fitQuality
>
FIT_TOLERANCE
)
{
g_log
.
warning
()
<<
"Fit quality obtained in PSIBackgroundSubtraction is poor. "
<<
"Skipping background calculation for WorkspaceIndex: "
<<
static_cast
<
int
>
(
index
)
<<
"
\n
"
;
}
else
{
backgroundValues
[
index
]
=
background
;
}
// Create background workspace
IAlgorithm_sptr
wsAlg
=
createChildAlgorithm
(
"CreateWorkspace"
,
0.7
,
1.0
);
wsAlg
->
setProperty
<
std
::
vector
<
double
>>
(
"DataX"
,
std
::
vector
<
double
>
(
2
,
0.0
));
wsAlg
->
setProperty
<
std
::
vector
<
double
>>
(
"DataY"
,
std
::
move
(
backgroundValues
));
wsAlg
->
setProperty
<
int
>
(
"NSpec"
,
static_cast
<
int
>
(
numberOfHistograms
));
wsAlg
->
execute
();
MatrixWorkspace_sptr
backgroundWS
=
wsAlg
->
getProperty
(
"OutputWorkspace"
);
backgroundWS
->
setYUnit
(
"Counts"
);
// Subtract background from inputworkspace
auto
minusAlg
=
createChildAlgorithm
(
"Minus"
);
minusAlg
->
setProperty
(
"LHSWorkspace"
,
inputWorkspace
);
minusAlg
->
setProperty
(
"RHSWorkspace"
,
backgroundWS
);
minusAlg
->
setProperty
(
"OutputWorkspace"
,
inputWorkspace
);
minusAlg
->
execute
();
}
/**
* Setup the fit algorithm used to obtain the background from PSI data
of a FlatBackground and ExpDecayMuon, on the second half of the PSI data.
* @param wsName :: Input workspace name
* @return Initalised fitting algorithm
*/
IAlgorithm_sptr
PSIBackgroundSubtraction
::
setupFitAlgorithm
(
const
std
::
string
&
wsName
)
{
std
::
string
functionstring
=
"name=FlatBackground,A0=0;name=ExpDecayMuon"
;
IFunction_sptr
func
=
FunctionFactory
::
Instance
().
createInitialized
(
functionstring
);
IAlgorithm_sptr
fit
=
createChildAlgorithm
(
"Fit"
);
int
maxIterations
=
getProperty
(
"MaxIterations"
);
fit
->
initialize
();
fit
->
setProperty
(
"Function"
,
func
);
fit
->
setProperty
(
"MaxIterations"
,
maxIterations
);
fit
->
setPropertyValue
(
"Minimizer"
,
MINIMISER
);
fit
->
setProperty
(
"CreateOutput"
,
false
);
fit
->
setProperty
(
"InputWorkspace"
,
wsName
);
return
fit
;
}
/**
* Runs the fit algorithm used to obtain the background from PSI data.
* @param fit :: Reference to the child fit algorithm used to obtain the
* background
* @param maxTime :: The maximum value of time (x-axis), used for fit range
* @param workspaceIndex :: Index which the fit will be performed on
* @return A tuple of background and fit quality
*/
std
::
tuple
<
double
,
double
>
PSIBackgroundSubtraction
::
calculateBackgroundFromFit
(
IAlgorithm_sptr
&
fit
,
const
std
::
pair
<
double
,
double
>
&
range
,
const
int
&
workspaceIndex
)
{
fit
->
setProperty
(
"StartX"
,
range
.
first
);
fit
->
setProperty
(
"EndX"
,
range
.
second
);
fit
->
setProperty
(
"WorkspaceIndex"
,
workspaceIndex
);
fit
->
execute
();
IFunction_sptr
func
=
fit
->
getProperty
(
"Function"
);
double
flatbackground
=
func
->
getParameter
(
"f0.A0"
);
double
chi2
=
std
::
stod
(
fit
->
getPropertyValue
(
"OutputChi2overDof"
));
return
std
::
make_tuple
(
flatbackground
,
chi2
);
}
// Create background workspace
IAlgorithm_sptr
wsAlg
=
createChildAlgorithm
(
"CreateWorkspace"
,
0.7
,
1.0
);
wsAlg
->
setProperty
<
std
::
vector
<
double
>>
(
"DataX"
,
std
::
vector
<
double
>
(
2
,
0.0
));
wsAlg
->
setProperty
<
std
::
vector
<
double
>>
(
"DataY"
,
std
::
move
(
backgroundValues
));
wsAlg
->
setProperty
<
int
>
(
"NSpec"
,
static_cast
<
int
>
(
numberOfHistograms
));
wsAlg
->
execute
();
MatrixWorkspace_sptr
backgroundWS
=
wsAlg
->
getProperty
(
"OutputWorkspace"
);
backgroundWS
->
setYUnit
(
"Counts"
);
// Subtract background from inputworkspace
auto
minusAlg
=
createChildAlgorithm
(
"Minus"
);
minusAlg
->
setProperty
(
"LHSWorkspace"
,
inputWorkspace
);
minusAlg
->
setProperty
(
"RHSWorkspace"
,
backgroundWS
);
minusAlg
->
setProperty
(
"OutputWorkspace"
,
inputWorkspace
);
minusAlg
->
execute
();
}
/**
* Setup the fit algorithm used to obtain the background from PSI data
of a FlatBackground and ExpDecayMuon, on the second half of the PSI data.
* @param wsName :: Input workspace name
* @return Initalised fitting algorithm
*/
IAlgorithm_sptr
PSIBackgroundSubtraction
::
setupFitAlgorithm
(
const
std
::
string
&
wsName
)
{
std
::
string
functionstring
=
"name=FlatBackground,A0=0;name=ExpDecayMuon"
;
IFunction_sptr
func
=
FunctionFactory
::
Instance
().
createInitialized
(
functionstring
);
IAlgorithm_sptr
fit
=
createChildAlgorithm
(
"Fit"
);
int
maxIterations
=
getProperty
(
"MaxIterations"
);
fit
->
initialize
();
fit
->
setProperty
(
"Function"
,
func
);
fit
->
setProperty
(
"MaxIterations"
,
maxIterations
);
fit
->
setPropertyValue
(
"Minimizer"
,
MINIMISER
);
fit
->
setProperty
(
"CreateOutput"
,
false
);
fit
->
setProperty
(
"InputWorkspace"
,
wsName
);
return
fit
;
}
/**
* Runs the fit algorithm used to obtain the background from PSI data.
* @param fit :: Reference to the child fit algorithm used to obtain the
* background
* @param maxTime :: The maximum value of time (x-axis), used for fit range
* @param workspaceIndex :: Index which the fit will be performed on
* @return A tuple of background and fit quality
*/
std
::
tuple
<
double
,
double
>
PSIBackgroundSubtraction
::
calculateBackgroundFromFit
(
IAlgorithm_sptr
&
fit
,
const
std
::
pair
<
double
,
double
>
&
range
,
const
int
&
workspaceIndex
)
{
fit
->
setProperty
(
"StartX"
,
range
.
first
);
fit
->
setProperty
(
"EndX"
,
range
.
second
);
fit
->
setProperty
(
"WorkspaceIndex"
,
workspaceIndex
);
fit
->
execute
();
IFunction_sptr
func
=
fit
->
getProperty
(
"Function"
);
double
flatbackground
=
func
->
getParameter
(
"f0.A0"
);
double
chi2
=
std
::
stod
(
fit
->
getPropertyValue
(
"OutputChi2overDof"
));
return
std
::
make_tuple
(
flatbackground
,
chi2
);
}
}
// namespace Muon
}
// namespace M
uon
}
// namespace M
antid
Framework/Muon/test/PSIBackgroundSubtractionTest.h
View file @
264344dc
...
...
@@ -20,7 +20,8 @@ namespace {
constexpr
char
*
WORKSPACE_NAME
=
"DummyWS"
;
MatrixWorkspace_sptr
createCountsTestWorkspace
(
const
size_t
numberOfHistograms
,
const
size_t
numberOfBins
,
bool
addLogs
=
true
)
{
const
size_t
numberOfBins
,
bool
addLogs
=
true
)
{
MatrixWorkspace_sptr
ws
=
WorkspaceCreationHelper
::
create2DWorkspace
(
numberOfHistograms
,
numberOfBins
);
...
...
@@ -67,8 +68,10 @@ public:
}
private:
std
::
tuple
<
double
,
double
>
calculateBackgroundFromFit
(
IAlgorithm_sptr
&
,
const
std
::
pair
<
double
,
double
>
&
range
,
const
int
&
index
)
override
{
std
::
tuple
<
double
,
double
>
calculateBackgroundFromFit
(
IAlgorithm_sptr
&
,
const
std
::
pair
<
double
,
double
>
&
range
,
const
int
&
index
)
override
{
return
std
::
make_tuple
(
m_background
,
m_fitQuality
);
}
double
m_background
{
0.00
};
...
...
@@ -121,11 +124,10 @@ public:
auto
ws
=
createCountsTestWorkspace
(
numberOfHistograms
,
numberOfBins
,
false
);
for
(
int
index
=
0
;
index
<
numberOfHistograms
;
index
++
)
{
ws
->
mutableRun
().
addProperty
(
"First good spectra "
+
std
::
to_string
(
index
),
-
1
);
ws
->
mutableRun
().
addProperty
(
"First good spectra "
+
std
::
to_string
(
index
),
-
1
);
ws
->
mutableRun
().
addProperty
(
"Last good spectra "
+
std
::
to_string
(
index
),
numberOfBins
-
10
);
numberOfBins
-
10
);
}
alg
.
initialize
();
...
...
@@ -135,7 +137,7 @@ public:
clearADS
();
}
void
test_that_algorithm_does_not_execute_if_bad_last_good_data
()
{
void
test_that_algorithm_does_not_execute_if_bad_last_good_data
()
{
PSIBackgroundSubtraction
alg
;
int
numberOfHistograms
=
2
;
int
numberOfBins
=
100
;
...
...
@@ -145,7 +147,7 @@ public:
ws
->
mutableRun
().
addProperty
(
"First good spectra "
+
std
::
to_string
(
index
),
1
);
ws
->
mutableRun
().
addProperty
(
"Last good spectra "
+
std
::
to_string
(
index
),
numberOfBins
*
2
);
numberOfBins
*
2
);
}
alg
.
initialize
();
...
...
@@ -155,7 +157,7 @@ public:
clearADS
();
}
void
test_that_algorithm_does_not_execute_if_last_before_first_good_data
()
{
void
test_that_algorithm_does_not_execute_if_last_before_first_good_data
()
{
PSIBackgroundSubtraction
alg
;
int
numberOfHistograms
=
2
;
int
numberOfBins
=
100
;
...
...
Write
Preview
Supports
Markdown
0%
Try again
or
attach a new file
.
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment