Re #11868 Trying different options for the fit. This one seems to work

mantidproject · Sep 29, 2015 · 1b3a9de · 1b3a9de
1 parent 5e35b04
commit 1b3a9de
Show file tree

Hide file tree

Showing 2 changed files with 103 additions and 2 deletions.
diff --git a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/CalMuonDetectorPhases.h b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/CalMuonDetectorPhases.h
@@ -54,6 +54,9 @@ class DLLExport CalMuonDetectorPhases : public API::Algorithm {
   API::MatrixWorkspace_sptr
   prepareWorkspace(const API::MatrixWorkspace_sptr &ws, double startTime,
                    double endTime);
+  API::ITableWorkspace_sptr fitWorkspace(const API::MatrixWorkspace_sptr &ws,
+                                         double freq);
+  std::string createFittingFunction(int nspec, double freq);
 };
 } // namespace Algorithms
 } // namespace Mantid

diff --git a/Code/Mantid/Framework/Algorithms/src/CalMuonDetectorPhases.cpp b/Code/Mantid/Framework/Algorithms/src/CalMuonDetectorPhases.cpp
@@ -1,6 +1,9 @@
 #include "MantidAlgorithms/CalMuonDetectorPhases.h"
 
+#include "MantidAPI/FunctionFactory.h"
 #include "MantidAPI/ITableWorkspace.h"
+#include "MantidAPI/IFunction.h"
+#include "MantidAPI/MultiDomainFunction.h"
 
 namespace Mantid {
 namespace Algorithms {
@@ -61,11 +64,106 @@ void CalMuonDetectorPhases::exec() {
   // removes exponential decay
   API::MatrixWorkspace_sptr tempWS = prepareWorkspace(inputWS, startTime, endTime);
 
-  API::ITableWorkspace_sptr tab =
-      API::WorkspaceFactory::Instance().createTable("TableWorkspace");
+  int nspec = static_cast<int>(inputWS->getNumberHistograms());
+
+  // Create the fitting function f(x) = A * sin ( w * x + p )
+  std::string funcStr = createFittingFunction(nspec, freq);
+  // Create the function from string
+  API::IFunction_sptr func =
+      API::FunctionFactory::Instance().createInitialized(funcStr);
+  // Create the multi domain function
+  boost::shared_ptr<API::MultiDomainFunction> multi =
+      boost::dynamic_pointer_cast<API::MultiDomainFunction>(func);
+  // Set the domain indices
+  for (int i = 0; i < nspec; ++i) {
+    multi->setDomainIndex(i, i);
+  }
+
+  API::IAlgorithm_sptr fit = createChildAlgorithm("Fit");
+  fit->initialize();
+  fit->setProperty("Function",
+                   boost::dynamic_pointer_cast<API::IFunction>(multi));
+  fit->setProperty("InputWorkspace", inputWS);
+  fit->setProperty("WorkspaceIndex", 0);
+  for (int s = 1; s < nspec; s++) {
+    std::string suffix = boost::lexical_cast<std::string>(s);
+    fit->setProperty("InputWorkspace_" + suffix, inputWS);
+    fit->setProperty("WorkspaceIndex_" + suffix, s);
+  }
+  fit->setProperty("CreateOutput", true);
+  fit->execute();
+  API::ITableWorkspace_sptr tab = fit->getProperty("OutputParameters");
+
+  // Set the result table
   setProperty("DetectorTable", tab);
 }
 
+/** TODO Description
+*/
+API::ITableWorkspace_sptr CalMuonDetectorPhases::fitWorkspace(const API::MatrixWorkspace_sptr &ws, double freq) {
+
+  int nspec = static_cast<int>(ws->getNumberHistograms());
+
+  // Create the fitting function f(x) = A * sin ( w * x + p )
+  std::string funcStr = createFittingFunction(nspec, freq);
+  // Create the function from string
+  API::IFunction_sptr func =
+      API::FunctionFactory::Instance().createInitialized(funcStr);
+  // Create the multi domain function
+  boost::shared_ptr<API::MultiDomainFunction> multi =
+      boost::dynamic_pointer_cast<API::MultiDomainFunction>(func);
+  // Set the domain indices
+  for (int i = 0; i < nspec; ++i) {
+    multi->setDomainIndex(i, i);
+  }
+
+  API::IAlgorithm_sptr fit = createChildAlgorithm("Fit");
+  fit->initialize();
+  fit->setProperty("Function",
+                   boost::dynamic_pointer_cast<API::IFunction>(multi));
+  fit->setProperty("InputWorkspace", ws);
+  fit->setProperty("WorkspaceIndex", 0);
+  for (int s = 1; s < nspec; s++) {
+    std::string suffix = boost::lexical_cast<std::string>(s);
+    fit->setProperty("InputWorkspace_" + suffix, ws);
+    fit->setProperty("WorkspaceIndex_" + suffix, s);
+  }
+  fit->setProperty("CreateOutput", true);
+  fit->execute();
+  API::ITableWorkspace_sptr tab = fit->getProperty("OutputParameters");
+
+  return tab;
+}
+
+/** TODO Description
+*/
+std::string CalMuonDetectorPhases::createFittingFunction(int nspec, double freq) {
+
+  // The fitting function is:
+  // f(x) = A * sin ( w * x + p )
+  // where w is shared across workspaces
+  std::ostringstream ss;
+  ss << "composite=MultiDomainFunction,NumDeriv=true;";
+  for (int s = 0; s < nspec; s++) {
+    ss << "name=UserFunction,";
+    ss << "Formula=A*sin(w*x+p),";
+    ss << "A=1.0,";
+    ss << "w=" << freq << ",";
+    ss << "p=1.0;";
+  }
+  ss << "ties=(";
+  for (int s = 1; s < nspec - 1; s++) {
+    ss << "f";
+    ss << boost::lexical_cast<std::string>(s);
+    ss << ".w=f0.w,";
+  }
+  ss << "f";
+  ss << boost::lexical_cast<std::string>(nspec - 1);
+  ss << ".w=f0.w)";
+
+  return ss.str();
+}
+
 /** TODO Description
 */
 API::MatrixWorkspace_sptr