Merge pull request #13312 from mantidproject/11534_calculate_q_resolution_in_Q1D

11534 calculate q resolution in Q1D

Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/Qhelper.h

@@ -37,6 +37,11 @@ namespace Algorithms {
 */
 class Qhelper {
 public:
+  void examineInput(API::MatrixWorkspace_const_sptr dataWS,
+                    API::MatrixWorkspace_const_sptr binAdj,
+                    API::MatrixWorkspace_const_sptr detectAdj,
+                    API::MatrixWorkspace_const_sptr qResolution);
+
   void examineInput(API::MatrixWorkspace_const_sptr dataWS,
                     API::MatrixWorkspace_const_sptr binAdj,
                     API::MatrixWorkspace_const_sptr detectAdj);

Code/Mantid/Framework/Algorithms/src/Q1D2.cpp

@@ -86,6 +86,11 @@ void Q1D2::init() {
   declareProperty(
       "ExtraLength", 0.0, mustBePositive,
       "Additional length for gravity correction.");
+
+  declareProperty(
+      new WorkspaceProperty<>("QResolution", "", Direction::Input,
+                              PropertyMode::Optional, dataVal),
+      "Workspace to calculate the Q resolution.\n");
 }
 /**
   @ throw invalid_argument if the workspaces are not mututially compatible
@@ -95,12 +100,14 @@ void Q1D2::exec() {
   MatrixWorkspace_const_sptr waveAdj = getProperty("WavelengthAdj");
   MatrixWorkspace_const_sptr pixelAdj = getProperty("PixelAdj");
   MatrixWorkspace_const_sptr wavePixelAdj = getProperty("WavePixelAdj");
+  MatrixWorkspace_const_sptr qResolution = getProperty("QResolution");
+
   const bool doGravity = getProperty("AccountForGravity");
   m_doSolidAngle = getProperty("SolidAngleWeighting");
 
   // throws if we don't have common binning or another incompatibility
   Qhelper helper;
-  helper.examineInput(m_dataWS, waveAdj, pixelAdj);
+  helper.examineInput(m_dataWS, waveAdj, pixelAdj, qResolution);
   // FIXME: how to examine the wavePixelAdj?
   g_log.debug() << "All input workspaces were found to be valid\n";
   // normalization as a function of wavelength (i.e. centers of x-value bins)
@@ -122,10 +129,15 @@ void Q1D2::exec() {
   MantidVec normSum(YOut.size(), 0.0);
   // the error on the normalisation
   MantidVec normError2(YOut.size(), 0.0);
+  // the averaged Q resolution
+  MantidVec &qResolutionOut = outputWS->dataDx(0);
 
   const int numSpec = static_cast<int>(m_dataWS->getNumberHistograms());
   Progress progress(this, 0.05, 1.0, numSpec + 1);
 
+  // Flag to decide if Q Resolution is to be used
+  auto useQResolution = qResolution ? true : false; 
+
   PARALLEL_FOR3(m_dataWS, outputWS, pixelAdj)
   for (int i = 0; i < numSpec; ++i) {
     PARALLEL_START_INTERUPT_REGION
@@ -160,6 +172,7 @@ void Q1D2::exec() {
       continue;
     }
 
+
     const size_t numWavbins = m_dataWS->readY(i).size() - wavStart;
     // make just one call to new to reduce CPU overhead on each thread, access
     // to these
@@ -184,6 +197,12 @@ void Q1D2::exec() {
     MantidVec::const_iterator YIn = m_dataWS->readY(i).begin() + wavStart;
     MantidVec::const_iterator EIn = m_dataWS->readE(i).begin() + wavStart;
 
+    // Pointers to the QResolution data. Note that the xdata was initially the same, hence
+    // the same indexing applies to the y values of m_dataWS and qResolution
+    // If we want to use it set it to the correct value, else to YIN, although that does not matter, as
+    // we won't use it
+    MantidVec::const_iterator QResIn = useQResolution ? (qResolution->readY(i).begin() + wavStart) : YIn;
+
     // when finding the output Q bin remember that the input Q bins (from the
     // convert to wavelength) start high and reduce
     MantidVec::const_iterator loc = QOut.end();
@@ -205,13 +224,20 @@ void Q1D2::exec() {
           // at the end
           EOutTo2[bin] += (*EIn) * (*EIn);
           normError2[bin] += *normETo2s;
+          if (useQResolution) {
+            qResolutionOut[bin] += (*YIn) * (*QResIn);
+          }
         }
       }
+
+      // Increment the QResolution iterator
+      if (useQResolution) {
+        ++QResIn;
+      }
     }
 
     PARALLEL_CRITICAL(q1d_spectra_map) {
       progress.report("Computing I(Q)");
-
       // Add up the detector IDs in the output spectrum at workspace index 0
       const ISpectrum *inSpec = m_dataWS->getSpectrum(i);
       ISpectrum *outSpec = outputWS->getSpectrum(0);
@@ -222,12 +248,30 @@ void Q1D2::exec() {
   }
   PARALLEL_CHECK_INTERUPT_REGION
 
+  if (useQResolution) {
+    // The number of Q (x)_ values is N, while the number of DeltaQ values is N-1,
+    // Richard Heenan suggested to duplicate the last entry of DeltaQ
+    Mantid::MantidVec::const_iterator countsIterator = YOut.begin();
+    Mantid::MantidVec::iterator qResolutionIterator = qResolutionOut.begin();
+    for (;countsIterator!= YOut.end(); ++countsIterator, ++qResolutionIterator) {
+      // Divide by the counts of the Qbin, if the counts are 0, the the qresolution will also be 0
+      if ((*countsIterator) > 0.0) {
+        *qResolutionIterator = (*qResolutionIterator)/(*countsIterator);
+      }
+    }
+    // Now duplicate write the second to last element into the last element of the deltaQ vector
+    if (qResolutionOut.size() > 1) {
+      qResolutionOut.rbegin()[0] = qResolutionOut.rbegin()[1];
+    }
+  }
+
   bool doOutputParts = getProperty("OutputParts");
   if (doOutputParts) {
     MatrixWorkspace_sptr ws_sumOfCounts =
         WorkspaceFactory::Instance().create(outputWS);
     ws_sumOfCounts->dataX(0) = outputWS->dataX(0);
     ws_sumOfCounts->dataY(0) = outputWS->dataY(0);
+    ws_sumOfCounts->dataDx(0) = outputWS->dataDx(0);
     for (size_t i = 0; i < outputWS->dataE(0).size(); i++) {
       ws_sumOfCounts->dataE(0)[i] = sqrt(outputWS->dataE(0)[i]);
     }
@@ -235,6 +279,7 @@ void Q1D2::exec() {
     MatrixWorkspace_sptr ws_sumOfNormFactors =
         WorkspaceFactory::Instance().create(outputWS);
     ws_sumOfNormFactors->dataX(0) = outputWS->dataX(0);
+    ws_sumOfNormFactors->dataDx(0) = outputWS->dataDx(0);
     for (size_t i = 0; i < ws_sumOfNormFactors->dataY(0).size(); i++) {
       ws_sumOfNormFactors->dataY(0)[i] = normSum[i];
       ws_sumOfNormFactors->dataE(0)[i] = sqrt(normError2[i]);

Code/Mantid/Framework/Algorithms/src/Qhelper.cpp

@@ -17,6 +17,46 @@ using namespace Kernel;
 using namespace API;
 using namespace Geometry;
 
+/** Checks if workspaces input to Q1D or Qxy are reasonable
+  @param dataWS data workspace
+  @param binAdj (WavelengthAdj) workpace that will be checked to see if it has
+  one spectrum and the same number of bins as dataWS
+  @param detectAdj (PixelAdj) passing NULL for this wont raise an error, if set
+  it will be checked this workspace has as many histograms as dataWS each with
+  one bin
+  @param qResolution: the QResolution workspace
+  @throw invalid_argument if the workspaces are not mututially compatible
+*/
+void Qhelper::examineInput(API::MatrixWorkspace_const_sptr dataWS,
+                           API::MatrixWorkspace_const_sptr binAdj,
+                           API::MatrixWorkspace_const_sptr detectAdj,
+                           API::MatrixWorkspace_const_sptr qResolution) {
+
+  // Check the compatibility of dataWS, binAdj and detectAdj
+  examineInput(dataWS, binAdj, detectAdj);
+
+  // Check the compatibility of the QResolution workspace
+  if (qResolution) {
+    // We require the same number of histograms
+    if (qResolution->getNumberHistograms() != dataWS->getNumberHistograms()) {
+      throw std::invalid_argument("The QResolution should have one spectrum"
+                                   "per spectrum of the input workspace");
+    }
+
+    // We require the same binning for the input workspace and the q resolution
+    // workspace
+    MantidVec::const_iterator reqX = dataWS->readX(0).begin();
+    MantidVec::const_iterator qResX = qResolution->readX(0).begin();
+    for (; reqX != dataWS->readX(0).end(); ++reqX, ++qResX) {
+      if (*reqX != *qResX) {
+        throw std::invalid_argument("The QResolution needs to have the same binning as"
+                                    "as the input workspace.");
+      }
+    }
+  }
+}
+
+
 /** Checks if workspaces input to Q1D or Qxy are reasonable
   @param dataWS data workspace
   @param binAdj (WavelengthAdj) workpace that will be checked to see if it has
@@ -103,6 +143,7 @@ void Qhelper::examineInput(API::MatrixWorkspace_const_sptr dataWS,
   }
 }
 
+
 /** Finds the first index number of the first wavelength bin that should
 * included based on the
 *  the calculation: W = Wcut (Rcut-R)/Rcut

Code/Mantid/Framework/Algorithms/test/Q1D2Test.h

@@ -13,6 +13,8 @@
 
 #include "MantidTestHelpers/WorkspaceCreationHelper.h"
 
+#include <algorithm>
+
 using namespace Mantid::API;
 using namespace Mantid::Kernel;
 using namespace Mantid::DataHandling;
@@ -24,6 +26,8 @@ static double flat_cell061Es [] = {8.140295E-03, 8.260089E-03, 8.198814E-03, 8.3
 /// defined below this creates some input data
 void createInputWorkspaces(int start, int end, Mantid::API::MatrixWorkspace_sptr & input, Mantid::API::MatrixWorkspace_sptr & wave, Mantid::API::MatrixWorkspace_sptr & pixels);
 void createInputWorkSpacesForMasking ( Mantid::API::MatrixWorkspace_sptr & input, Mantid::API::MatrixWorkspace_sptr & wave, Mantid::API::MatrixWorkspace_sptr & pixels);
+void createQResolutionWorkspace(Mantid::API::MatrixWorkspace_sptr & qResolution, Mantid::API::MatrixWorkspace_sptr & alteredInput, Mantid::API::MatrixWorkspace_sptr & input,  double value1, double value2);
+
 
 class Q1D2Test : public CxxTest::TestSuite
 {
@@ -341,7 +345,7 @@ public:
 
     Mantid::API::AnalysisDataService::Instance().remove(m_noGrav);
   }
-    
+
   void testInvalidInput()
   {
     Mantid::Algorithms::Q1D2 Q1D;
@@ -362,8 +366,94 @@ public:
       Q1D.execute();
 
     TS_ASSERT( ! Q1D.isExecuted() )
+
+    // Restore the old binning
+    xData[15] -= 0.001;
   }
-  
+
+  void testRunsWithQResolution() {
+    // Arrange
+    Mantid::API::MatrixWorkspace_sptr qResolution, alteredInput;
+    const double value1 = 1;
+    const double value2 = 2;
+    createQResolutionWorkspace(qResolution, alteredInput, m_inputWS, value1, value2);
+
+    // Act
+    Mantid::Algorithms::Q1D2 Q1D;
+    TS_ASSERT_THROWS_NOTHING( Q1D.initialize() );
+    TS_ASSERT( Q1D.isInitialized() )
+
+    const std::string outputWS("Q1D2Test_result");
+    Q1D.setProperty("DetBankWorkspace", alteredInput); 
+    Q1D.setProperty("WavelengthAdj", m_wavNorm);
+    Q1D.setProperty("PixelAdj", m_pixel);
+    Q1D.setPropertyValue("OutputWorkspace", outputWS);
+    Q1D.setPropertyValue("OutputBinning", "0.1,-0.02,0.5");
+    Q1D.setProperty("QResolution", qResolution);
+    Q1D.setProperty("OutputParts", true);
+
+    // Assert
+    TS_ASSERT_THROWS_NOTHING(Q1D.execute());
+
+    TS_ASSERT( Q1D.isExecuted() )
+
+    Mantid::API::MatrixWorkspace_sptr result;
+    TS_ASSERT_THROWS_NOTHING( result = boost::dynamic_pointer_cast<Mantid::API::MatrixWorkspace>
+                              (Mantid::API::AnalysisDataService::Instance().retrieve(outputWS)) )
+
+    Mantid::API::MatrixWorkspace_sptr sumOfCounts;
+    TS_ASSERT_THROWS_NOTHING( sumOfCounts = boost::dynamic_pointer_cast<Mantid::API::MatrixWorkspace>
+                              (Mantid::API::AnalysisDataService::Instance().retrieve(outputWS+"_sumOfCounts")) )
+
+    // That output will be SUM_i(Yin_i*QRES_in_i)/(SUM_i(Y_in_i)) for each q value
+    // In our test workspace we set QRes_in_1 to 1, this means that all DX values should be SUM_i(Y_in_i)/(SUM_i(Y_in_i))
+    // which is either 1 or 0. It can be 0 if no data falls into this bin. We make sure that there is at least one bin
+    // with a count of 1
+    auto& dataDX = result->dataDx(0);
+    unsigned int counter = 0;
+    for (auto it = dataDX.begin(); it != dataDX.end(); ++it) {
+      if (*it==1.0) {
+        counter++;
+      }
+      // Make sure that the value is either 1 or 0
+      TSM_ASSERT("The DX value should be either 1 or 0", (*it==1.0) || (*it==0.0));
+    }
+    TSM_ASSERT("There should be at least one bin which is not 0", counter > 0);
+
+    // Clean Up
+    Mantid::API::AnalysisDataService::Instance().remove(outputWS);
+    Mantid::API::AnalysisDataService::Instance().remove(outputWS+"_sumOfCounts");
+    Mantid::API::AnalysisDataService::Instance().remove(outputWS+"_sumOfNormFactors");
+  }
+
+    void testThatDxIsNotPopulatedWhenQResolutionIsNotChoosen() {
+    // Arrange
+    Mantid::Algorithms::Q1D2 Q1D;
+    TS_ASSERT_THROWS_NOTHING( Q1D.initialize() );
+    TS_ASSERT( Q1D.isInitialized() )
+
+    const std::string outputWS("Q1D2Test_result");
+    Q1D.setProperty("DetBankWorkspace", m_inputWS); 
+    Q1D.setProperty("WavelengthAdj", m_wavNorm);
+    Q1D.setProperty("PixelAdj", m_pixel);
+    Q1D.setPropertyValue("OutputWorkspace", outputWS);
+    Q1D.setPropertyValue("OutputBinning", "0.1,-0.02,0.5");
+    TS_ASSERT_THROWS_NOTHING(Q1D.execute());
+
+    TS_ASSERT( Q1D.isExecuted() )
+
+    Mantid::API::MatrixWorkspace_sptr result;
+    TS_ASSERT_THROWS_NOTHING( result = boost::dynamic_pointer_cast<Mantid::API::MatrixWorkspace>
+                              (Mantid::API::AnalysisDataService::Instance().retrieve(outputWS)) )
+
+    // Make sure that the Q resolution is not calculated
+    auto& dataDX = result->dataDx(0);
+    for (auto it = dataDX.begin(); it != dataDX.end(); ++it) {
+      TSM_ASSERT("All Dx values should be 0, as we didn't use the QResolution ooption", (*it==0.0));
+    }
+    Mantid::API::AnalysisDataService::Instance().remove(outputWS);
+  }
+
   ///stop the constructor from being run every time algorithms test suite is initialised
   static Q1D2Test *createSuite() { return new Q1D2Test(); }
   static void destroySuite(Q1D2Test *suite) { delete suite; }
@@ -483,4 +573,25 @@ void createInputWorkSpacesForMasking ( Mantid::API::MatrixWorkspace_sptr & input
 {
    createInputWorkspaces ( 9001, 9030, input, wave, pixels );
 }
+
+void createQResolutionWorkspace(Mantid::API::MatrixWorkspace_sptr & qResolution, Mantid::API::MatrixWorkspace_sptr & alteredInput, Mantid::API::MatrixWorkspace_sptr & input, double value1, double value2) {
+  //The q resolution workspace is almost the same to the input workspace, except for the y value, we set all Y values to 1
+  qResolution = Mantid::API::MatrixWorkspace_sptr(input->clone().release());
+  alteredInput = Mantid::API::MatrixWorkspace_sptr(input->clone().release());
+
+  // Populate Y with Value1
+  for (size_t i = 0; i < qResolution->getNumberHistograms(); ++i) {
+    auto& data = qResolution->dataY(i);
+    std::fill(data.begin(), data.end(), value1);
+  }
+
+  // Populate Y with Value2
+  for (size_t i = 0; i < alteredInput->getNumberHistograms(); ++i) {
+    auto& data = alteredInput->dataY(i);
+    std::fill(data.begin(), data.end(), value2);
+  }
+
+
+}
+
 #endif /*Q1D2Test_H_*/