15#include "MantidHistogramData/Histogram.h"
44 return "Inelastic\\Utility;Diagnostics;Events\\EventFiltering";
49 return "Calculates instrument count rate as the function of the "
50 "experiment time and adds CountRate log to the source workspace.";
60 "Name of the event workspace to calculate counting rate for.");
62 "Minimal value of X-range for the rate calculations. If left "
63 "to default, Workspace X-axis minimal value is used.");
65 "Maximal value of X-range for the rate calculations. If left "
66 "to default, Workspace X-axis maximal value is used.");
68 "RangeUnits",
"Energy",
71 "The units from Mantid Unit factory for calculating the "
72 "counting rate and XMin-XMax ranges are in. If the "
73 "X-axis of the input workspace is not expressed"
74 "in these units, unit conversion will be performed, so the "
75 "workspace should contain all necessary information for this "
76 "conversion. E.g. if *RangeUnits* is *EnergyTransfer*, Ei "
77 "log containing incident energy value should be attached to the "
78 "input workspace. See ConvertUnits algorithm for the details.");
79 std::vector<std::string> propOptions{
"Elastic",
"Direct",
"Indirect"};
80 declareProperty(
"EMode",
"Elastic", std::make_shared<Kernel::StringListValidator>(propOptions),
81 "The energy mode for 'RangeUnits' conversion mode (default: elastic)");
84 std::string used_logs_mode(
"Used normalization logs");
86 "Usually you want to normalize counting rate to some "
87 "rate related to the source beam intensity. Change this to "
88 "'false' if appropriate time series log is broken || not attached to "
89 "the input workspace.");
91 "If the normalization log contains "
92 "cumulative counting, derivative "
93 "of this log is necessary to get "
94 "correct normalization values.");
96 "The name of the log, used in the counting rate normalization. ");
103 "The name of the processed time series log with instrument "
104 "count rate to be added"
105 " to the source workspace");
107 "If true, the count rate log will have the normalization log "
108 "accuracy; If false, the 'NumTimeSteps' in the visualization "
109 "workspace below will be used for the target log granularity too.");
112 std::string spur_vis_mode(
"Spurion visualization");
115 "Optional name to build 2D matrix workspace for spurion visualization. "
116 "If name is provided, a 2D workspace with this name will be created "
117 "containing data to visualize counting rate as function of time in the "
120 auto mustBeReasonable = std::make_shared<Kernel::BoundedValidator<int>>();
121 mustBeReasonable->setLower(3);
124 "Number of time steps (time accuracy) the visualization workspace has. "
125 "Also number of steps in 'CountRateLogName' log if "
126 "'UseNormLogGranularity' is set to false. Should be bigger than 3");
129 "Number of steps (accuracy) of the visualization workspace has along "
144 throw std::runtime_error(
"Event workspace " + sourceWS->getName() +
145 " contains events without necessary frame "
146 "information. Can not process counting rate");
162 const std::string logname =
getProperty(
"CountRateLogName");
163 auto newlog = std::make_unique<Kernel::TimeSeriesProperty<double>>(logname);
166 sourceWS->mutableRun().addProperty(std::move(newlog),
true);
188 std::vector<double> countNormalization;
192 std::unique_ptr<std::mutex[]> pVisWS_locks;
194 pVisWS_locks.reset(
new std::mutex[
m_visWs->getNumberHistograms()]);
197 auto nHist =
static_cast<int64_t
>(InputWorkspace->getNumberHistograms());
201 auto dTRangeMin =
static_cast<double>(
m_TRangeMin.totalNanoseconds());
202 auto dTRangeMax =
static_cast<double>(
m_TRangeMax.totalNanoseconds());
203 std::vector<MantidVec> Buff;
211 Buff.resize(nThreads);
216 for (int64_t i = 0; i < nHist; ++i) {
236 for (
int i = 0; i < nThreads; i++) {
237 countRate[j] += Buff[i][j];
240 if (!countNormalization.empty()) {
241 countRate[j] /= countNormalization[j];
244 if (!countNormalization.empty() && this->buildVisWS()) {
246 for (int64_t j = 0; j < int64_t(
m_visNorm.size()); j++) {
254 double dt = (dTRangeMax - dTRangeMin) /
static_cast<double>(
m_numLogSteps);
257 times[i] = Types::Core::DateAndTime(t0 +
static_cast<int64_t
>((0.5 +
double(i)) * dt));
273 for (
const Types::Event::TofEvent &ev : events) {
274 double pulsetime =
static_cast<double>(ev.pulseTime().totalNanoseconds());
275 double tof = ev.tof();
276 if (pulsetime < m_visT0 || pulsetime >=
m_visTmax)
281 auto n_spec =
static_cast<size_t>((pulsetime -
m_visT0) /
m_visDT);
283 (spectraLocks + n_spec)->lock();
284 auto &
Y =
m_visWs->mutableY(n_spec);
286 (spectraLocks + n_spec)->unlock();
305 std::string NormLogName =
getProperty(
"NormalizationLogName");
306 std::string TargetLog =
getProperty(
"CountRateLogName");
307 if (NormLogName == TargetLog) {
308 throw std::invalid_argument(
"Target log name: " + TargetLog +
" and normalization log name: " + NormLogName +
309 " can not be the same");
314 bool useLogAccuracy =
getProperty(
"UseNormLogGranularity");
316 bool logPresent = InputWorkspace->run().hasProperty(NormLogName);
320 "' values requested but the log is not attached to the "
321 "workspace. Normalization disabled");
324 g_log.
warning() <<
"Normalization by log: '" << NormLogName
325 <<
"' -- log derivative requested but the source log is "
327 "the workspace. Log derivative will not be used.\n";
330 if (useLogAccuracy) {
331 g_log.
warning() <<
"Using accuracy of the log: '" << NormLogName
332 <<
"' is requested but the log is not attached to the "
333 "workspace. Will use accuracy defined by "
334 "'NumTimeSteps' property value.\n";
335 useLogAccuracy =
false;
352 if (!useLogAccuracy) {
353 g_log.
warning() <<
"Change of the counting log accuracy while "
354 "normalizing by log values is not implemented. Will "
355 "use log accuracy.\n";
356 useLogAccuracy =
true;
361 Types::Core::DateAndTime runTMin, runTMax;
362 InputWorkspace->getPulseTimeMinMax(runTMin, runTMax);
364 if (useLogAccuracy) {
365 Types::Core::DateAndTime tLogMin, tLogMax;
368 auto pSource = InputWorkspace->run().getTimeSeriesProperty<
double>(NormLogName);
369 tLogMin = pSource->firstTime();
370 tLogMax = pSource->lastTime();
376 if (tLogMin < runTMin || tLogMax > runTMax) {
377 if (tLogMin > runTMax || tLogMax < runTMin) {
380 " time lies outside of the the whole experiment time. "
381 "Log normalization impossible.");
382 useLogAccuracy =
false;
392 if (tLogMin > runTMin || tLogMax < runTMax) {
394 " time does not cover the whole experiment time. "
395 "Log normalization impossible.");
396 useLogAccuracy =
false;
401 if (useLogAccuracy) {
406 " smaller then 2. Can not normalize using this log.");
408 useLogAccuracy =
false;
414 double t_epsilon = double(runTMax.totalNanoseconds()) * (1 + std::numeric_limits<double>::epsilon());
415 auto eps_increment =
static_cast<int64_t
>(t_epsilon - double(runTMax.totalNanoseconds()));
418 if (useLogAccuracy) {
423 auto iTMax = runTMax.totalNanoseconds();
430 " is not consistent with number of log steps. "
431 "Can not use this log normalization");
432 useLogAccuracy =
false;
435 if (iTMax1 <= iTMax) {
443 if (!useLogAccuracy) {
463 std::string RangeUnits =
getProperty(
"RangeUnits");
464 auto axis = InputWorkspace->getAxis(0);
465 const auto unit = axis->unit();
468 if (unit->unitID() != RangeUnits) {
470 std::string wsName = InputWorkspace->getName();
471 if (wsName.empty()) {
472 wsName =
"_CountRate_UnitsConverted";
474 wsName =
"_" + wsName +
"_converted";
477 conv->setProperty(
"InputWorkspace", InputWorkspace);
478 conv->setPropertyValue(
"OutputWorkspace", wsName);
480 conv->setProperty(
"Emode", Emode);
481 conv->setProperty(
"Target", RangeUnits);
483 conv->setRethrows(
true);
485 wst = conv->getProperty(
"OutputWorkspace");
488 wst = InputWorkspace;
490 m_workingWS = std::dynamic_pointer_cast<DataObjects::EventWorkspace>(wst);
492 throw std::runtime_error(
"SetWSDataRanges:Can not retrieve EventWorkspace "
493 "after converting units");
505 double realMin, realMax;
510 m_XRangeMax = realMax * (1. + std::numeric_limits<double>::epsilon());
519 m_XRangeMax = realMax * (1. + std::numeric_limits<double>::epsilon());
522 g_log.
debug() <<
"Workspace constrain min range changed from: " <<
m_XRangeMin <<
" To: " << realMin << std::endl;
526 g_log.
debug() <<
"Workspace constrain max range changed from: " <<
m_XRangeMax <<
" To: " << realMax << std::endl;
527 m_XRangeMax = realMax * (1. + std::numeric_limits<double>::epsilon());
530 if (m_XRangeMax < realMin || m_XRangeMin > realMax) {
537 throw std::invalid_argument(
" Minimal spurion search data limit is bigger "
538 "than the maximal limit. ( Min: " +
548 std::string visWSName =
getProperty(
"VisualizationWs");
549 if (visWSName.empty()) {
560 "workspace exceeds the number of points in the "
561 "normalization log. This mode is not supported so "
562 "number of time steps decreased to be equal to "
563 "the number of normalization log points\n";
568 std::string RangeUnits =
getProperty(
"RangeUnits");
570 m_visWs = create<Workspace2D>(numTBins, BinEdges(numXBins + 1));
575 if (std::isinf(Xmax)) {
577 for (int64_t i = Xbin.size() - 1; i >= 0; --i) {
578 if (!std::isinf(Xbin[i])) {
583 if (std::isinf(Xmax)) {
584 g_log.
warning() <<
"All X-range for visualization workspace is infinity. "
585 "Can not build visualization workspace in the units "
595 std::vector<double> xx(numXBins);
596 for (
int i = 0; i < numXBins; ++i) {
597 xx[i] =
m_XRangeMin + (0.5 +
static_cast<double>(i)) * dX;
599 auto ax0 = std::make_unique<API::NumericAxis>(xx);
600 ax0->setUnit(RangeUnits);
601 m_visWs->replaceAxis(0, std::move(ax0));
605 static_cast<double>(numTBins)) *
608 for (
int i = 0; i < numTBins; i++) {
609 xx[i] = (0.5 +
static_cast<double>(i)) * dt;
611 auto ax1 = std::make_unique<API::NumericAxis>(xx);
613 labelY->setLabel(
"sec");
614 ax1->unit() = labelY;
615 m_visWs->replaceAxis(1, std::move(ax1));
653 g_log.
warning() <<
"CalculateCountRate::buildVisWSNormalization: No source "
654 "normalization log is found. Will not normalize "
655 "visualization workspace\n";
661 throw std::runtime_error(
"Can not retrieve Y-axis from visualization workspace");
663 normalization.assign(ax->length(), 0.);
#define DECLARE_ALGORITHM(classname)
#define PARALLEL_THREAD_NUMBER
#define PARALLEL_START_INTERRUPT_REGION
Begins a block to skip processing is the algorithm has been interupted Note the end of the block if n...
#define PARALLEL_END_INTERRUPT_REGION
Ends a block to skip processing is the algorithm has been interupted Note the start of the block if n...
#define PARALLEL_NUMBER_OF_THREADS
#define PARALLEL_CHECK_INTERRUPT_REGION
Adds a check after a Parallel region to see if it was interupted.
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
virtual std::shared_ptr< Algorithm > createChildAlgorithm(const std::string &name, const double startProgress=-1., const double endProgress=-1., const bool enableLogging=true, const int &version=-1)
Create a Child Algorithm.
Class to represent a numeric axis of a workspace.
Helper class for reporting progress from algorithms.
A property class for workspaces.
In normal circumstances an instrument in event mode counts neutrons with constant steady rate which d...
void setSourceWSandXRanges(DataObjects::EventWorkspace_sptr &InputWorkspace)
void calcRateLog(DataObjects::EventWorkspace_sptr &InputWorkspace, Kernel::TimeSeriesProperty< double > *const targLog)
Process input workspace to calculate instrument counting rate as function of experiment time.
bool m_doVis
should algo generate visualization VS
void setOutLogParameters(const DataObjects::EventWorkspace_sptr &InputWorkspace)
void checkAndInitVisWorkspace()
Check if visualization workspace is necessary and initiate it if requested.
bool m_normalizeResult
pointer to the log used to normalize results or NULL if no such log present on input workspace.
bool buildVisWS() const
helper function to test if visualization workspace is requested
DataObjects::EventWorkspace_sptr m_workingWS
temporary workspace used to keep intermediate results
double m_XRangeMin
spurion search ranges (TOF or other units requested)
bool normalizeCountRate() const
Helper function: true if count rate should be normalized and false otherwise.
std::vector< double > m_visNorm
void disableNormalization(const std::string &NormLogError)
Disable normalization using normalization log.
Types::Core::DateAndTime m_TRangeMin
experiment time ranges:
DataObjects::Workspace2D_sptr m_visWs
shared pointer to the optional visualization workspace
Types::Core::DateAndTime m_TRangeMax
Kernel::TimeSeriesProperty< double > const * m_pNormalizationLog
const std::string summary() const override
Algorithm's summary for use in the GUI and help.
const std::string category() const override
Algorithm's category for identification.
bool m_rangeExplicit
specify if rate is calculated in selected frame interval (range defined) or all frame should be used
const std::string name() const override
Algorithms name for identification.
void init() override
Declare the algorithm's properties.
void histogramEvents(const DataObjects::EventList &el, std::mutex *spectraLocks)
histogram event list into visualization workspace
std::unique_ptr< Kernel::TimeSeriesProperty< double > > m_tmpLogHolder
void buildVisWSNormalization(std::vector< double > &normalization)
method to prepare normalization vector for the visualisation workspace using data from normalization ...
int m_numLogSteps
default number of points in the target log
int version() const override
Algorithm's version for identification.
bool useLogDerivative() const
Helper function to test if log derivative is used.
void exec() override
Execute the algorithm.
std::vector< Types::Event::TofEvent > & getEvents()
Return the list of TofEvents contained.
void generateCountsHistogramPulseTime(const double &xMin, const double &xMax, MantidVec &Y, const double TofMin=std::numeric_limits< double >::lowest(), const double TofMax=std::numeric_limits< double >::max()) const
With respect to PulseTime fill a histogram given equal histogram bins.
bool empty() const
Much like stl containers, returns true if there is nothing in the event list.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void setPropertyGroup(const std::string &name, const std::string &group)
Set the group for a given property.
void debug(const std::string &msg)
Logs at debug level.
void warning(const std::string &msg)
Logs at warning level.
void information(const std::string &msg)
Logs at information level.
Validator to check that a property is not left empty.
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
The concrete, templated class for properties.
const std::string & name() const
Get the property's name.
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
A specialised Property class for holding a series of time-value pairs.
void replaceValues(const std::vector< Types::Core::DateAndTime > ×, const std::vector< TYPE > &values)
Replaces the time series with new values time series values.
TimeSeriesProperty< TYPE > * clone() const override
"Virtual" copy constructor
std::vector< TYPE > valuesAsVector() const
Return the time series's values (unfiltered) as a vector<TYPE>
Types::Core::DateAndTime firstTime() const
Returns the first time regardless of filter.
void histogramData(const Types::Core::DateAndTime &tMin, const Types::Core::DateAndTime &tMax, std::vector< double > &counts) const
generate constant time-step histogram from the property values
Types::Core::DateAndTime lastTime() const
Returns the last time.
std::unique_ptr< TimeSeriesProperty< double > > getDerivative() const
Return time series property, containing time derivative of current property.
int realSize() const override
Returns the real size of the time series property map:
EventType
What kind of event list is being stored.
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
std::shared_ptr< EventWorkspace > EventWorkspace_sptr
shared pointer to the EventWorkspace class
std::unique_ptr< T > create(const P &parent, const IndexArg &indexArg, const HistArg &histArg)
This is the create() method that all the other create() methods call.
std::vector< double > MantidVec
typedef for the data storage used in Mantid matrix workspaces
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
std::string to_string(const wide_integer< Bits, Signed > &n)
@ InOut
Both an input & output workspace.
@ Input
An input workspace.
@ Output
An output workspace.