Refs #14181 error of volume

Framework/Crystal/src/OptimizeLatticeForCellType.cpp

@@ -92,7 +92,7 @@ void OptimizeLatticeForCellType::exec() {
   runWS.push_back(ws);
 
   if (perRun) {
-    std::vector<std::pair<std::string, bool>> criteria;
+    std::vector<std::pair<std::string, bool> > criteria;
     // Sort by run number
     criteria.push_back(std::pair<std::string, bool>("runnumber", true));
     ws->sort(criteria);
@@ -130,7 +130,8 @@ void OptimizeLatticeForCellType::exec() {
     IAlgorithm_sptr fit_alg;
     try {
       fit_alg = createChildAlgorithm("Fit", -1, -1, false);
-    } catch (Exception::NotFoundError &) {
+    }
+    catch (Exception::NotFoundError &) {
       g_log.error("Can't locate Fit algorithm");
       throw;
     }
@@ -146,15 +147,13 @@ void OptimizeLatticeForCellType::exec() {
 
     double chisq = fit_alg->getProperty("OutputChi2overDoF");
     Geometry::UnitCell refinedCell = latticeFunction->getUnitCell();
-    /*std::vector<double> sigabc;
-    for (size_t i = 0; i < latticeFunction->nParams(); i++)
-      sigabc.push_back(latticeFunction->getError(i));*/
 
     IAlgorithm_sptr ub_alg;
     try {
       ub_alg =
           createChildAlgorithm("FindUBUsingLatticeParameters", -1, -1, false);
-    } catch (Exception::NotFoundError &) {
+    }
+    catch (Exception::NotFoundError &) {
       g_log.error("Can't locate FindUBUsingLatticeParameters algorithm");
       throw;
     }
@@ -178,8 +177,44 @@ void OptimizeLatticeForCellType::exec() {
                        refinedCell.errorc(), refinedCell.erroralpha(),
                        refinedCell.errorbeta(), refinedCell.errorgamma());
 
+    // From latcon.py by Art Schultz
+    double alpha1 = refinedCell.alpha() - 0.5 * refinedCell.erroralpha();
+    double Va1 =
+        UnitCell(refinedCell.a(), refinedCell.b(), refinedCell.c(), alpha1,
+                 refinedCell.beta(), refinedCell.gamma()).volume();
+    double alpha2 = refinedCell.alpha() + 0.5 * refinedCell.erroralpha();
+    double Va2 =
+        UnitCell(refinedCell.a(), refinedCell.b(), refinedCell.c(), alpha2,
+                 refinedCell.beta(), refinedCell.gamma()).volume();
+    double delta_V_alpha = Va2 - Va1;
+
+    double beta1 = refinedCell.beta() - 0.5 * refinedCell.errorbeta();
+    Va1 = UnitCell(refinedCell.a(), refinedCell.b(), refinedCell.c(),
+                   refinedCell.alpha(), beta1, refinedCell.gamma()).volume();
+    double beta2 = refinedCell.beta() + 0.5 * refinedCell.errorbeta();
+    Va2 = UnitCell(refinedCell.a(), refinedCell.b(), refinedCell.c(),
+                   refinedCell.alpha(), beta2, refinedCell.gamma()).volume();
+    double delta_V_beta = Va2 - Va1;
+
+    double gamma1 = refinedCell.gamma() - 0.5 * refinedCell.errorgamma();
+    Va1 = UnitCell(refinedCell.a(), refinedCell.b(), refinedCell.c(),
+                   refinedCell.alpha(), refinedCell.beta(), gamma1).volume();
+    double gamma2 = refinedCell.gamma() + 0.5 * refinedCell.errorgamma();
+    Va2 = UnitCell(refinedCell.a(), refinedCell.b(), refinedCell.c(),
+                   refinedCell.alpha(), refinedCell.beta(), gamma2).volume();
+    double delta_V_gamma = Va2 - Va1;
+
+    double sigV = refinedCell.volume() *
+                  sqrt(std::pow(refinedCell.errora() / refinedCell.a(), 2) +
+                       std::pow(refinedCell.errorb() / refinedCell.b(), 2) +
+                       std::pow(refinedCell.errorc() / refinedCell.c(), 2) +
+                       std::pow(delta_V_alpha / refinedCell.volume(), 2) +
+                       std::pow(delta_V_beta / refinedCell.volume(), 2) +
+                       std::pow(delta_V_gamma / refinedCell.volume(), 2));
+
     // Show the modified lattice parameters
-    g_log.notice() << runWS[i_run]->getName() << "  " << o_lattice << "\n";
+    g_log.notice() << runWS[i_run]->getName() << "  " << o_lattice << "  "
+                   << sigV << "\n";
 
     runWS[i_run]->mutableSample().setOrientedLattice(&o_lattice);
 
@@ -205,9 +240,9 @@ void OptimizeLatticeForCellType::exec() {
                                                runWS[i_run]->getName() +
                                                ".integrate");
       savePks_alg->executeAsChildAlg();
-      g_log.notice() << "See output file: "
-                     << outputdir + "ls" + runWS[i_run]->getName() +
-                            ".integrate"
+      g_log.notice() << "See output file: " << outputdir + "ls" +
+                                                   runWS[i_run]->getName() +
+                                                   ".integrate"
                      << "\n";
       // Save UB
       Mantid::API::IAlgorithm_sptr saveUB_alg =

Framework/Geometry/inc/MantidGeometry/Crystal/UnitCell.h

@@ -121,6 +121,7 @@ public:
   void setErroralpha(double _alphaerr, const int angleunit = angDegrees);
   void setErrorbeta(double _betaerr, const int angleunit = angDegrees);
   void setErrorgamma(double _gammaerr, const int angleunit = angDegrees);
+  void setErrorvolume(double _volerr);
   // Get errors in latice parameters
   double errora() const;
   double errorb() const;
@@ -128,6 +129,7 @@ public:
   double erroralpha(const int angleunit = angDegrees) const;
   double errorbeta(const int angleunit = angDegrees) const;
   double errorgamma(const int angleunit = angDegrees) const;
+  double errorvolume() const;
   // Access private variables
   const Kernel::DblMatrix &getG() const;
   const Kernel::DblMatrix &getGstar() const;

Framework/Geometry/src/Crystal/UnitCell.cpp

@@ -21,7 +21,7 @@ UnitCell::UnitCell() : da(6), ra(6), errorda(6), G(3, 3), Gstar(3, 3), B(3, 3) {
   da[0] = da[1] = da[2] = 1.;
   da[3] = da[4] = da[5] = deg2rad * 90.0;
   errorda[0] = errorda[1] = errorda[2] = errorda[3] = errorda[4] = errorda[5] =
-      0.0;
+      errorda[6] = 0.0;
   recalculate();
 }
 
@@ -50,7 +50,7 @@ UnitCell::UnitCell(double _a, double _b, double _c)
   // Angles are 90 degrees in radians ->Pi/2
   da[3] = da[4] = da[5] = 0.5 * M_PI;
   errorda[0] = errorda[1] = errorda[2] = errorda[3] = errorda[4] = errorda[5] =
-      0.0;
+      errorda[6] = 0.0;
   recalculate();
 }
 
@@ -75,7 +75,7 @@ UnitCell::UnitCell(double _a, double _b, double _c, double _alpha, double _beta,
     da[5] = _gamma;
   }
   errorda[0] = errorda[1] = errorda[2] = errorda[3] = errorda[4] = errorda[5] =
-      0.0;
+      errorda[6] = 0.0;
   recalculate();
 }
 
@@ -266,6 +266,11 @@ double UnitCell::errorgamma(const int angleunit) const {
   }
 }
 
+/** Get lattice parameter error
+@return errorc :: errorlattice parameter \f$ volume \f$ (in \f$ \mbox{\AA} \f$ )
+*/
+double UnitCell::errorvolume() const { return errorda[6]; }
+
 /** Set lattice parameters
   @param _a, _b, _c, _alpha, _beta, _gamma :: lattice parameters\n
 @param angleunit :: units for angle, of type #AngleUnits . Default is degrees.
@@ -405,6 +410,11 @@ void UnitCell::setErrorgamma(double _gammaerr, const int angleunit) {
     errorda[5] = _gammaerr;
 }
 
+/** Set lattice parameter error
+@param _cerr :: lattice parameter \f$ volume \f$ error (in \f$ \mbox{\AA} \f$
+)*/
+void UnitCell::setErrorvolume(double _volerr) { errorda[6] = _volerr; }
+
 /// Return d-spacing (\f$ \mbox{ \AA } \f$) for a given h,k,l coordinate
 double UnitCell::d(double h, double k, double l) const {
   return 1.0 / dstar(V3D(h, k, l));
@@ -575,7 +585,8 @@ std::ostream &operator<<(std::ostream &out, const UnitCell &unitCell) {
       << std::setw(12) << unitCell.c() << std::fixed << std::setprecision(6)
       << std::setw(12) << unitCell.alpha() << std::fixed << std::setprecision(6)
       << std::setw(12) << unitCell.beta() << std::fixed << std::setprecision(6)
-      << std::setw(12) << unitCell.gamma();
+      << std::setw(12) << unitCell.gamma() << std::fixed << std::setprecision(6)
+      << std::setw(12) << unitCell.volume();
 
   // write out the uncertainty if there is a positive one somewhere
   if ((unitCell.errora() > 0) || (unitCell.errorb() > 0) ||
@@ -589,7 +600,8 @@ std::ostream &operator<<(std::ostream &out, const UnitCell &unitCell) {
         << std::setw(12) << unitCell.erroralpha() << std::fixed
         << std::setprecision(6) << std::setw(12) << unitCell.errorbeta()
         << std::fixed << std::setprecision(6) << std::setw(12)
-        << unitCell.errorgamma();
+        << unitCell.errorgamma() << std::fixed << std::setprecision(6)
+        << std::setw(12) << unitCell.errorvolume();
 
   return out;
 }
@@ -607,12 +619,12 @@ std::string unitCellToStr(const UnitCell &unitCell) {
 
 UnitCell strToUnitCell(const std::string &unitCellString) {
   boost::char_separator<char> separator(" ");
-  boost::tokenizer<boost::char_separator<char>> cellTokens(unitCellString,
-                                                           separator);
+  boost::tokenizer<boost::char_separator<char> > cellTokens(unitCellString,
+                                                            separator);
 
   std::vector<double> components;
 
-  for (boost::tokenizer<boost::char_separator<char>>::iterator token =
+  for (boost::tokenizer<boost::char_separator<char> >::iterator token =
            cellTokens.begin();
        token != cellTokens.end(); ++token) {
     components.push_back(boost::lexical_cast<double>(*token));