Refs #4693 add OptimizeExtinctionParameters algorithm

Code/Mantid/Framework/Crystal/CMakeLists.txt

@@ -15,6 +15,7 @@ set ( SRC_FILES
 	src/LoadIsawUB.cpp
 	src/MaskPeaksWorkspace.cpp
 	src/NormaliseVanadium.cpp
+	src/OptimizeExtinctionParameters.cpp
 	src/PeakIntegration.cpp
 	src/PeakIntensityVsRadius.cpp
 	src/PredictPeaks.cpp
@@ -40,12 +41,14 @@ set ( INC_FILES
 	inc/MantidCrystal/FindUBUsingMinMaxD.h
 	inc/MantidCrystal/IndexPeaks.h
     inc/MantidCrystal/IndexSXPeaks.h
+        inc/MantidCrystal/GSLFunctions.h
 	inc/MantidCrystal/IntegratePeakTimeSlices.h
 	inc/MantidCrystal/LoadHKL.h
 	inc/MantidCrystal/LoadIsawPeaks.h
 	inc/MantidCrystal/LoadIsawUB.h
 	inc/MantidCrystal/MaskPeaksWorkspace.h
 	inc/MantidCrystal/NormaliseVanadium.h
+	inc/MantidCrystal/OptimizeExtinctionParameters.h
 	inc/MantidCrystal/PeakIntegration.h
 	inc/MantidCrystal/PeakIntensityVsRadius.h
 	inc/MantidCrystal/PredictPeaks.h

Code/Mantid/Framework/Crystal/inc/MantidCrystal/OptimizeExtinctionParameters.h

@@ -0,0 +1,74 @@
+#ifndef MANTID_CRYSTAL_OptimizeExtinctionParameters_H_
+#define MANTID_CRYSTAL_OptimizeExtinctionParameters_H_
+
+#include "MantidAPI/Algorithm.h"
+#include "MantidKernel/System.h"
+#include "MantidDataObjects/OffsetsWorkspace.h"
+#include "MantidGeometry/Crystal/PointGroup.h"
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_multifit_nlin.h>
+#include <gsl/gsl_multimin.h>
+#include <gsl/gsl_statistics.h>
+
+namespace Mantid
+{
+namespace Algorithms
+{
+/**
+ Find the offsets for each detector
+
+ @author Vickie Lynch, SNS
+ @date 02/06/2012
+
+ Copyright &copy; 2009 ISIS Rutherford Appleton Laboratory & NScD Oak Ridge National Laboratory
+
+ This file is part of Mantid.
+
+ Mantid is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ Mantid is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+ File change history is stored at: <https://svn.mantidproject.org/mantid/trunk/Code/Mantid>
+ Code Documentation is available at: <http://doxygen.mantidproject.org>
+ */
+class DLLExport OptimizeExtinctionParameters: public API::Algorithm
+{
+public:
+  /// Default constructorMatrix
+  OptimizeExtinctionParameters();
+  /// Destructor
+  virtual ~OptimizeExtinctionParameters();
+  /// Algorithm's name for identification overriding a virtual method
+  virtual const std::string name() const { return "OptimizeExtinctionParameters"; }
+  /// Algorithm's version for identification overriding a virtual method
+  virtual int version() const { return 1; }
+  /// Algorithm's category for identification overriding a virtual method
+  virtual const std::string category() const { return "Crystal"; }
+  /// Call TOFExtinction as a sub-algorithm to get statistics of the peaks
+  double fitMosaic(double mosaic, double rcrystallite, std::string inname, std::string corrOption, std::string pointOption, std::string tofParams);
+
+private:
+   /// Point Groups possible
+  std::vector<Mantid::Geometry::PointGroup_sptr> m_pointGroups;
+
+  /// Sets documentation strings for this algorithm
+  virtual void initDocs();
+  // Overridden Algorithm methods
+  void init();
+  void exec();
+
+};
+
+} // namespace Algorithm
+} // namespace Mantid
+
+#endif /*MANTID_CRYSTAL_OptimizeExtinctionParameters_H_*/

Code/Mantid/Framework/Crystal/src/OptimizeExtinctionParameters.cpp

@@ -0,0 +1,238 @@
+/*WIKI* 
+
+
+*WIKI*/
+#include "MantidCrystal/OptimizeExtinctionParameters.h"
+#include "MantidCrystal/GSLFunctions.h"
+#include "MantidDataObjects/PeaksWorkspace.h"
+#include "MantidGeometry/Crystal/PointGroup.h"
+#include "MantidAPI/FileProperty.h"
+#include "MantidAPI/FunctionFactory.h"
+#include "MantidAPI/SpectraDetectorMap.h"
+#include "MantidAPI/WorkspaceValidators.h"
+#include "MantidAPI/IPeakFunction.h"
+#include "MantidAPI/IBackgroundFunction.h"
+#include "MantidAPI/CompositeFunctionMW.h"
+#include "MantidKernel/VectorHelper.h"
+#include <boost/math/special_functions/fpclassify.hpp>
+#include <fstream>
+#include <iomanip>
+#include <ostream>
+#include <sstream>
+
+using namespace Mantid::Geometry;
+
+namespace Mantid
+{
+  namespace Algorithms
+  {
+
+    // Register the class into the algorithm factory
+    DECLARE_ALGORITHM(OptimizeExtinctionParameters)
+
+    /// Sets documentation strings for this algorithm
+    void OptimizeExtinctionParameters::initDocs()
+    {
+      this->setWikiSummary("Finds optimal mosaic and r_crystallite parameters for extinction correction.");
+      this->setOptionalMessage("Finds optimal mosaic and r_crystallite parameters for extinction correction.");
+    }
+
+    using namespace Kernel;
+    using namespace API;
+    using std::size_t;
+    using namespace DataObjects;
+
+    /// Constructor
+    OptimizeExtinctionParameters::OptimizeExtinctionParameters()
+    {
+      m_pointGroups = getAllPointGroups();
+    }
+
+    /// Destructor
+    OptimizeExtinctionParameters::~OptimizeExtinctionParameters()
+    {}
+
+    static double gsl_costFunction(const gsl_vector *v, void *params)
+    {
+      std::string *p = (std::string *)params;
+      std::string inname = p[0];
+      std::string corrOption = p[1];
+      std::string pointOption = p[2];
+      std::string tofParams = p[3];
+      double mosaic = gsl_vector_get(v,0);
+      std::vector<double> tofParam = Kernel::VectorHelper::splitStringIntoVector<double>(tofParams);
+      double rcrystallite = tofParam[4];
+      if(v->size>1)rcrystallite = gsl_vector_get(v,1);
+      Mantid::Algorithms::OptimizeExtinctionParameters u;
+      return u.fitMosaic(mosaic, rcrystallite, inname, corrOption, pointOption, tofParams);
+    }
+
+    //-----------------------------------------------------------------------------------------
+    /** Initialisation method. Declares properties to be used in algorithm.
+     */
+    void OptimizeExtinctionParameters::init()
+    {
+
+    declareProperty(new WorkspaceProperty<PeaksWorkspace>("InputWorkspace","",Direction::InOut),
+        "An input PeaksWorkspace with an instrument.");
+    //declareProperty(new WorkspaceProperty<PeaksWorkspace>("OutputWorkspace","",Direction::Output));
+    std::vector<std::string> corrOptions;
+    corrOptions.push_back("Type I Zachariasen" );
+    corrOptions.push_back("Type I Gaussian" );
+    corrOptions.push_back("Type I Lorentzian" );
+    corrOptions.push_back("Type II Zachariasen" );
+    corrOptions.push_back("Type II Gaussian" );
+    corrOptions.push_back("Type II Lorentzian" );
+    corrOptions.push_back("Type I&II Zachariasen" );
+    corrOptions.push_back("Type I&II Gaussian" );
+    corrOptions.push_back("Type I&II Lorentzian" );
+    corrOptions.push_back( "None, Scaling Only" );
+    declareProperty("ExtinctionCorrectionType", corrOptions[0],new ListValidator(corrOptions),
+      "Select the type of extinction correction.");
+
+    BoundedValidator<double> *mustBePositive = new BoundedValidator<double> ();
+    mustBePositive->setLower(0.0);
+    declareProperty("LinearScatteringCoef", -1.0, mustBePositive,
+      "Linear scattering coefficient in 1/cm");
+    declareProperty("LinearAbsorptionCoef", -1.0, mustBePositive->clone(),
+      "Linear absorption coefficient at 1.8 Angstroms in 1/cm");
+    declareProperty("Radius", 0.10, "Radius of spherical crystal in cm");
+    declareProperty("Cell", 255.0, "Unit Cell Volume (Angstroms^3)");
+    declareProperty("Mosaic", 0.262, "Mosaic Spread (FWHM) (Degrees)");
+    declareProperty("RCrystallite", 6.0, "Becker-Coppens Crystallite Radius (micron)");
+    std::vector<std::string> propOptions;
+    for (size_t i=0; i<m_pointGroups.size(); ++i)
+      propOptions.push_back( m_pointGroups[i]->getName() );
+    declareProperty("PointGroup", propOptions[0],new ListValidator(propOptions),
+      "Which point group applies to this crystal?");
+
+      //Disable default gsl error handler (which is to call abort!)
+      gsl_set_error_handler_off();
+    }
+
+    //-----------------------------------------------------------------------------------------
+    /** Executes the algorithm
+     *
+     *  @throw Exception::FileError If the grouping file cannot be opened or read successfully
+     */
+    void OptimizeExtinctionParameters::exec()
+    {
+      std::string par[4];
+      std::string inname = getProperty("InputWorkspace");
+      par[0] = inname;
+      std::string type = getProperty("ExtinctionCorrectionType");
+      par[1] = type;
+      std::string group = getProperty("PointGroup");
+      par[2] = group;
+      double smu = getProperty("LinearScatteringCoef");
+      double amu = getProperty("LinearAbsorptionCoef");
+      double radius = getProperty("Radius");
+      double mosaic = getProperty("Mosaic");
+      double cell = getProperty("Cell");
+      double r_crystallite = getProperty("RCrystallite");
+      std::ostringstream strwi;
+      strwi<<smu<<","<<amu<<","<<radius<<","<<cell<<","<<r_crystallite;
+      par[3] = strwi.str();
+
+      const gsl_multimin_fminimizer_type *T =
+      gsl_multimin_fminimizer_nmsimplex;
+      gsl_multimin_fminimizer *s = NULL;
+      gsl_vector *ss, *x;
+      gsl_multimin_function minex_func;
+
+      // finally do the fitting
+
+      size_t nopt = 1;
+      if(type.compare("II")>0)nopt = 2;
+      size_t iter = 0;
+      int status = 0;
+      double size;
+
+      /* Starting point */
+      x = gsl_vector_alloc (nopt);
+      gsl_vector_set (x, 0, mosaic);
+      if(nopt>1)gsl_vector_set (x, 1, r_crystallite);
+
+      /* Set initial step sizes to 0.001 */
+      ss = gsl_vector_alloc (nopt);
+      gsl_vector_set_all (ss, 0.001);
+
+      /* Initialize method and iterate */
+      minex_func.n = nopt;
+      minex_func.f = &Mantid::Algorithms::gsl_costFunction;
+      minex_func.params = &par;
+
+      s = gsl_multimin_fminimizer_alloc (T, nopt);
+      gsl_multimin_fminimizer_set (s, &minex_func, x, ss);
+
+      do
+      {
+        iter++;
+        status = gsl_multimin_fminimizer_iterate(s);
+        if (status)
+          break;
+
+        size = gsl_multimin_fminimizer_size (s);
+        status = gsl_multimin_test_size (size, 1e-4);
+
+      }
+      while (status == GSL_CONTINUE && iter < 50);
+
+      // Output summary to log file
+      std::string reportOfDiffractionEventCalibrateDetectors = gsl_strerror(status);
+      //g_log.debug() << 
+      mosaic = gsl_vector_get (s->x, 0);
+      if(nopt>1)r_crystallite = gsl_vector_get (s->x, 1);
+      std::cout <<
+        " Method used = " << " Simplex" << 
+        " Iteration = " << iter << 
+        " Status = " << reportOfDiffractionEventCalibrateDetectors << 
+        " Minimize Sum = " << s->fval << 
+        " Mosaic   = " << mosaic << 
+        " R_crystallite  = " << r_crystallite << "  \n";
+      gsl_vector_free(x);
+      gsl_vector_free(ss);
+      gsl_multimin_fminimizer_free (s);
+
+    }
+
+
+    //-----------------------------------------------------------------------------------------
+   /** Calls Gaussian1D as a child algorithm to fit the offset peak in a spectrum
+    *
+    *  @param s :: The Workspace Index to fit
+    *  @return The calculated offset value
+    */
+
+    double OptimizeExtinctionParameters::fitMosaic(double mosaic, double rcrystallite, std::string inname, std::string corrOption, std::string pointOption, std::string tofParams)
+    {
+      PeaksWorkspace_sptr inputW = boost::dynamic_pointer_cast<PeaksWorkspace>
+           (AnalysisDataService::Instance().retrieve(inname));
+      std::vector<double> tofParam = Kernel::VectorHelper::splitStringIntoVector<double>(tofParams);
+
+      API::IAlgorithm_sptr tofextinction = createSubAlgorithm("TOFExtinction",0.0,0.2);
+      tofextinction->setProperty("InputWorkspace", inputW);
+      tofextinction->setProperty("OutputWorkspace", "tmp");
+      tofextinction->setProperty("ExtinctionCorrectionType", corrOption);
+      tofextinction->setProperty<double>("LinearScatteringCoef", tofParam[0]);
+      tofextinction->setProperty<double>("LinearAbsorptionCoef", tofParam[1]);
+      tofextinction->setProperty<double>("Radius", tofParam[2]);
+      tofextinction->setProperty<double>("Mosaic", mosaic);
+      tofextinction->setProperty<double>("Cell", tofParam[3]);
+      tofextinction->setProperty<double>("RCrystallite", rcrystallite);
+      tofextinction->executeAsSubAlg();
+      PeaksWorkspace_sptr peaksW = tofextinction->getProperty("OutputWorkspace");
+
+      API::IAlgorithm_sptr sorthkl = createSubAlgorithm("SortHKL",0.0,0.2);
+      sorthkl->setProperty("InputWorkspace", peaksW);
+      sorthkl->setProperty("OutputWorkspace", peaksW);
+      sorthkl->setProperty("PointGroup", pointOption);
+      sorthkl->executeAsSubAlg();
+      double Chisq = sorthkl->getProperty("OutputChi2");
+      std::cout << mosaic<<"  "<<rcrystallite<<"  "<<Chisq<<"\n";
+      return Chisq;
+    }
+
+
+  } // namespace Algorithm
+} // namespace Mantid

Code/Mantid/Framework/Crystal/src/SortHKL.cpp

@@ -191,7 +191,6 @@ namespace Crystal
     sig2.clear();
     setProperty<PeaksWorkspace_sptr>("OutputWorkspace", peaksW);
     setProperty("OutputChi2", Chisq);
-    std::cout << "Chisq = "<<Chisq<<"\n";
 
   }
   void SortHKL::Outliers(std::vector<double>& data, std::vector<double>& sig2)

Code/Mantid/Framework/Crystal/src/TOFExtinction.cpp

@@ -105,7 +105,7 @@ namespace Crystal
     declareProperty("Mosaic", 0.262, "Mosaic Spread (FWHM) (Degrees)");
     declareProperty("Cell", 255.0, "Unit Cell Volume (Angstroms^3)");
     declareProperty("RCrystallite", 6.0, "Becker-Coppens Crystallite Radius (micron)");
-    declareProperty("ScaleFactor", 1.2, "Multiply FSQ and sig(FSQ) by scaleFactor");
+    declareProperty("ScaleFactor", 1.0, "Multiply FSQ and sig(FSQ) by scaleFactor");
 
   }