Merge pull request #18891 from mantidproject/SNSPowderReduction_pchar…

…ge_fix

Fix bug in gd_prtn_chrg for chunked data

Framework/Algorithms/inc/MantidAlgorithms/NormaliseByCurrent.h

@@ -71,8 +71,8 @@ class DLLExport NormaliseByCurrent : public API::Algorithm {
   void init() override;
   void exec() override;
   // Extract the charge value from the logs.
-  double
-  extractCharge(boost::shared_ptr<Mantid::API::MatrixWorkspace> inputWS) const;
+  double extractCharge(boost::shared_ptr<Mantid::API::MatrixWorkspace> inputWS,
+                       const bool integratePCharge) const;
 };
 
 } // namespace Algorithm

Framework/Algorithms/src/NormaliseByCurrent.cpp

@@ -23,20 +23,26 @@ void NormaliseByCurrent::init() {
   declareProperty(make_unique<WorkspaceProperty<MatrixWorkspace>>(
                       "OutputWorkspace", "", Direction::Output),
                   "Name of the output workspace");
+  declareProperty(Kernel::make_unique<Kernel::PropertyWithValue<bool>>(
+                      "RecalculatePCharge", false, Kernel::Direction::Input),
+                  "Re-integrates the proton charge. This will modify the "
+                  "gd_prtn_chrg. Does nothing for multi-period data");
 }
 
 /**
  * Extract a value for the charge from the input workspace. Handles either
  * single period or multi-period data.
  *
  * @param inputWS :: The input workspace to extract the log details from.
- *
+ * @param integratePCharge :: recalculate the integrated proton charge if true
+
  * @throws Exception::NotFoundError, std::domain_error or
  * std::runtime_error if the charge value(s) are not set in the
  * workspace logs or if the values are invalid (0)
  */
 double NormaliseByCurrent::extractCharge(
-    boost::shared_ptr<Mantid::API::MatrixWorkspace> inputWS) const {
+    boost::shared_ptr<Mantid::API::MatrixWorkspace> inputWS,
+    const bool integratePCharge) const {
   // Get the good proton charge and check it's valid
   double charge(-1.0);
   const Run &run = inputWS->run();
@@ -82,6 +88,9 @@ double NormaliseByCurrent::extractCharge(
 
   } else {
     try {
+      if (integratePCharge) {
+        inputWS->run().integrateProtonCharge();
+      }
       charge = inputWS->run().getProtonCharge();
     } catch (Exception::NotFoundError &) {
       g_log.error() << "The proton charge is not set for the run attached to "
@@ -101,9 +110,10 @@ void NormaliseByCurrent::exec() {
   // Get the input workspace
   MatrixWorkspace_sptr inputWS = getProperty("InputWorkspace");
   MatrixWorkspace_sptr outputWS = getProperty("OutputWorkspace");
+  const bool integratePCharge = getProperty("RecalculatePCharge");
 
   // Get the good proton charge and check it's valid
-  double charge = extractCharge(inputWS);
+  double charge = extractCharge(inputWS, integratePCharge);
 
   g_log.information() << "Normalisation current: " << charge << " uamps\n";
 

Framework/Algorithms/test/NormaliseByCurrentTest.h

@@ -10,6 +10,7 @@
 #include "MantidAlgorithms/NormaliseByCurrent.h"
 #include "MantidDataObjects/EventWorkspace.h"
 #include "MantidKernel/ArrayProperty.h"
+#include "MantidKernel/TimeSeriesProperty.h"
 #include "MantidKernel/UnitFactory.h"
 
 using namespace Mantid;
@@ -21,7 +22,8 @@ using namespace Mantid::DataObjects;
 namespace {
 MatrixWorkspace_const_sptr doTest(MatrixWorkspace_sptr inWS,
                                   std::string wsNameOut, double expectedY,
-                                  double expectedE, bool performance = false) {
+                                  double expectedE, bool calcPcharge = false,
+                                  bool performance = false) {
   NormaliseByCurrent norm1;
   if (!norm1.isInitialized())
     norm1.initialize();
@@ -37,6 +39,8 @@ MatrixWorkspace_const_sptr doTest(MatrixWorkspace_sptr inWS,
   TS_ASSERT_THROWS_NOTHING(norm1.setProperty("InputWorkspace", inWS));
   TS_ASSERT_THROWS_NOTHING(
       norm1.setPropertyValue("OutputWorkspace", wsNameOut));
+  TS_ASSERT_THROWS_NOTHING(
+      norm1.setProperty("RecalculatePCharge", calcPcharge));
 
   TS_ASSERT_THROWS_NOTHING(norm1.execute());
   TS_ASSERT(norm1.isExecuted());
@@ -79,18 +83,18 @@ MatrixWorkspace_const_sptr doTest(MatrixWorkspace_sptr inWS,
 }
 
 MatrixWorkspace_const_sptr doTest(std::string wsNameIn, std::string wsNameOut,
-                                  double expectedY, double expectedE,
-                                  bool performance = false) {
+                                  const double pCharge, double expectedY,
+                                  double expectedE, bool performance = false) {
   MatrixWorkspace_sptr inWS =
       AnalysisDataService::Instance().retrieveWS<MatrixWorkspace>(wsNameIn);
 
   // Now set the charge
-  inWS->mutableRun().setProtonCharge(2.0);
+  inWS->mutableRun().setProtonCharge(pCharge);
   inWS->getAxis(0)->unit() =
       Mantid::Kernel::UnitFactory::Instance().create("TOF");
   inWS->setYUnit("Counts");
 
-  return doTest(inWS, wsNameOut, expectedY, expectedE, performance);
+  return doTest(inWS, wsNameOut, expectedY, expectedE, true, performance);
 }
 
 /// Helper method to add necessary log values to simulate multi-period data.
@@ -109,6 +113,25 @@ void addMultiPeriodLogsTo(MatrixWorkspace_sptr ws, int period,
   ws->mutableRun().addLogData(nperiodsProp);
   ws->mutableRun().addLogData(currentPeriodProp);
 }
+
+void addPChargeLogTo(MatrixWorkspace_sptr ws, const double pChargeAccum) {
+  auto pchargeLog =
+      Kernel::make_unique<Kernel::TimeSeriesProperty<double>>("proton_charge");
+
+  const Kernel::DateAndTime runstart(20000000000);
+  const int64_t pulsedt = 100 * 1000 * 1000;
+  const size_t numpulses = 100;
+  const double pCharge = pChargeAccum / static_cast<double>(numpulses);
+
+  for (int64_t pid = 0; pid < static_cast<int64_t>(numpulses); pid++) {
+    const int64_t pulsetime_i64 = pulsedt + runstart.totalNanoseconds();
+    const Kernel::DateAndTime pulsetime(pulsetime_i64);
+    pchargeLog->addValue(pulsetime, pCharge);
+  } // FOR each pulse
+
+  ws->mutableRun().addLogData(pchargeLog.release());
+  // ws->mutableRun().integrateProtonCharge(); // TODO
+}
 }
 class NormaliseByCurrentTest : public CxxTest::TestSuite {
 public:
@@ -140,7 +163,7 @@ class NormaliseByCurrentTest : public CxxTest::TestSuite {
   void test_exec() {
     AnalysisDataService::Instance().add(
         "normIn", WorkspaceCreationHelper::create2DWorkspaceBinned(10, 3, 1));
-    doTest("normIn", "normOut", 1.0, 0.5 * M_SQRT2);
+    doTest("normIn", "normOut", 2.0, 1.0, 0.5 * M_SQRT2);
     AnalysisDataService::Instance().remove("normIn");
     AnalysisDataService::Instance().remove("normOut");
   }
@@ -278,7 +301,7 @@ class NormaliseByCurrentTest : public CxxTest::TestSuite {
   void test_execInPlace() {
     AnalysisDataService::Instance().add(
         "normIn", WorkspaceCreationHelper::create2DWorkspaceBinned(10, 3, 1));
-    doTest("normIn", "normIn", 1.0, 0.5 * M_SQRT2);
+    doTest("normIn", "normIn", 2.0, 1.0, 0.5 * M_SQRT2);
     AnalysisDataService::Instance().remove("normIn");
   }
 
@@ -289,7 +312,25 @@ class NormaliseByCurrentTest : public CxxTest::TestSuite {
 
     EventWorkspace_const_sptr outputEvent;
     outputEvent = boost::dynamic_pointer_cast<const EventWorkspace>(
-        doTest("normInEvent", "normOutEvent", 1.0, 0.5 * M_SQRT2));
+        doTest("normInEvent", "normOutEvent", 2.0, 1.0, 0.5 * M_SQRT2));
+    // Output is an event workspace
+    TS_ASSERT(outputEvent);
+
+    AnalysisDataService::Instance().remove("normInEvent");
+    AnalysisDataService::Instance().remove("normOutEvent");
+  }
+
+  void test_execEventFunnyPCharge() {
+    auto wksp =
+        WorkspaceCreationHelper::createEventWorkspace(10, 3, 100, 0.0, 1.0, 2);
+    addPChargeLogTo(wksp, 2.); // pcharge intentionally doesn't match
+    AnalysisDataService::Instance().add("normInEvent", wksp);
+
+    // intentionally set the wrong `gd_prtn_chrg` to stress getting the right
+    // answer
+    EventWorkspace_const_sptr outputEvent;
+    outputEvent = boost::dynamic_pointer_cast<const EventWorkspace>(
+        doTest("normInEvent", "normOutEvent", 100.0, 1.0, 0.5 * M_SQRT2));
     // Output is an event workspace
     TS_ASSERT(outputEvent);
 
@@ -304,7 +345,7 @@ class NormaliseByCurrentTest : public CxxTest::TestSuite {
 
     EventWorkspace_const_sptr outputEvent;
     outputEvent = boost::dynamic_pointer_cast<const EventWorkspace>(
-        doTest("normInEvent", "normInEvent", 1.0, 0.5 * M_SQRT2));
+        doTest("normInEvent", "normInEvent", 2.0, 1.0, 0.5 * M_SQRT2));
     // Output is an event workspace
     TS_ASSERT(outputEvent);
 

Framework/PythonInterface/plugins/algorithms/SNSPowderReduction.py

@@ -707,7 +707,8 @@ def _loadAndSum(self, filename_list, outName, **filterWall):
             # temp_ws = self.get_workspace(sample_ws_name)
             if not (is_event_workspace(sample_ws_name) and get_workspace(sample_ws_name).getNumberEvents() == 0):
                 api.NormaliseByCurrent(InputWorkspace=sample_ws_name,
-                                       OutputWorkspace=sample_ws_name)
+                                       OutputWorkspace=sample_ws_name,
+                                       RecalculatePCharge=True)
                 get_workspace(sample_ws_name).getRun()['gsas_monitor'] = 1
             # END-IF
         # ENDI-IF
@@ -772,7 +773,8 @@ def _focusAndSum(self, filenames, filterWall, calib, preserveEvents=True):
 
         if self._normalisebycurrent is True:
             api.NormaliseByCurrent(InputWorkspace=sumRun,
-                                   OutputWorkspace=sumRun)
+                                   OutputWorkspace=sumRun,
+                                   RecalculatePCharge=True)
             get_workspace(sumRun).getRun()['gsas_monitor'] = 1
 
         return sumRun
@@ -945,7 +947,8 @@ def _focusChunks(self, filename, filterWall, calib,
             try:
                 if normalisebycurrent is True:
                     api.NormaliseByCurrent(InputWorkspace=output_wksp_list[split_index],
-                                           OutputWorkspace=output_wksp_list[split_index])
+                                           OutputWorkspace=output_wksp_list[split_index],
+                                           RecalculatePCharge=True)
                     get_workspace(output_wksp_list[split_index]).getRun()['gsas_monitor'] = 1
             except RuntimeError as e:
                 self.log().warning(str(e))