Re #7385. Merge branch 'master' of github.com:mantidproject/mantid

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

@@ -77,11 +77,11 @@ namespace Mantid
       double calculateQ(const double efixed,const int emode, const double deltaE,
                         const double twoTheta, const double azimuthal) const;
       /// Init variables cache base on the given workspace
-      void initCachedValues(API::MatrixWorkspace_const_sptr workspace);
+      void initCachedValues(const API::MatrixWorkspace_const_sptr & workspace);
       /// Init the theta index
-      void initThetaCache(API::MatrixWorkspace_const_sptr workspace);
+      void initAngularCachesNonPSD(const API::MatrixWorkspace_const_sptr & workspace);
       /// Get angles and calculate angular widths.
-      void getValuesAndWidths(API::MatrixWorkspace_const_sptr workspace);
+      void initAngularCachesPSD(const API::MatrixWorkspace_const_sptr &workspace);
 
       /// Create the output workspace
       DataObjects::RebinnedOutput_sptr setUpOutputWorkspace(API::MatrixWorkspace_const_sptr inputWorkspace,

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

@@ -23,6 +23,7 @@ two pixels at the same vertical position in adjacent tubes.
 #include "MantidAPI/SpectraAxis.h"
 #include "MantidAPI/WorkspaceFactory.h"
 #include "MantidGeometry/Instrument/DetectorGroup.h"
+#include "MantidGeometry/Instrument/ReferenceFrame.h"
 #include "MantidGeometry/Math/LaszloIntersection.h"
 #include "MantidGeometry/Math/Quadrilateral.h"
 #include "MantidGeometry/Math/Vertex2D.h"
@@ -41,6 +42,7 @@ namespace Algorithms
   DECLARE_ALGORITHM(SofQW3)
 
   using namespace Mantid::Kernel;
+  using namespace Mantid::Geometry;
   using namespace Mantid::API;
   using namespace Mantid::DataObjects;
   using Geometry::IDetector_const_sptr;
@@ -130,13 +132,13 @@ namespace Algorithms
     if (par.empty())
     {
       // Index theta cache
-      this->initThetaCache(inputWS);
+      this->initAngularCachesNonPSD(inputWS);
     }
     else
     {
       g_log.debug() << "Offset: " << par[0] << std::endl;
       this->m_detNeighbourOffset = static_cast<int>(par[0]);
-      this->getValuesAndWidths(inputWS);
+      this->initAngularCachesPSD(inputWS);
     }
 
     const MantidVec & X = inputWS->readX(0);
@@ -154,22 +156,11 @@ namespace Algorithms
       }
 
       double theta = this->m_theta[i];
-      double phi = 0.0;
-      double thetaWidth = 0.0;
-      double phiWidth = 0.0;
-      // Non-PSD mode
-      if (par.empty())
-      {
-        thetaWidth = this->m_thetaWidth;
-      }
-      // PSD mode
-      else
-      {
-        phi = this->m_phi[i];
-        thetaWidth = this->m_thetaWidths[i];
-        phiWidth = this->m_phiWidths[i];
-      }
+      double phi = this->m_phi[i];
+      double thetaWidth = this->m_thetaWidths[i];
+      double phiWidth = this->m_phiWidths[i];
 
+      // Compute polygon points
       double thetaHalfWidth = 0.5 * thetaWidth;
       double phiHalfWidth = 0.5 * phiWidth;
 
@@ -180,7 +171,7 @@ namespace Algorithms
       const double phiUpper = phi + phiHalfWidth;
 
       const double efixed = m_EmodeProperties.getEFixed(detector);
-
+      const specid_t specNo = inputWS->getSpectrum(i)->getSpectrumNo();
       for(size_t j = 0; j < nEnergyBins; ++j)
       {
         m_progress->report("Computing polygon intersections");
@@ -193,6 +184,13 @@ namespace Algorithms
         const V2D lr(dE_jp1, this->calculateQ(efixed,emode, dE_jp1, thetaLower, phiLower));
         const V2D ur(dE_jp1, this->calculateQ(efixed,emode, dE_jp1, thetaUpper, phiUpper));
         const V2D ul(dE_j, this->calculateQ(efixed,emode, dE_j, thetaUpper, phiUpper));
+        if(g_log.is(Logger::Priority::PRIO_DEBUG))
+        {
+          g_log.debug() << "Spectrum=" << specNo << ", theta=" << theta << ",thetaWidth=" << thetaWidth
+                              << ", phi=" << phi << ", phiWidth=" << phiWidth
+                              << ". QE polygon: ll=" << ll << ", lr=" << lr << ", ur=" << ur << ", ul=" << ul << "\n";
+        }
+
         Quadrilateral inputQ = Quadrilateral(ll, lr, ur, ul);
 
         this->rebinToFractionalOutput(inputQ, inputWS, i, j, outputWS, m_Qout);
@@ -245,16 +243,21 @@ namespace Algorithms
    * values are required very frequently so the total time is more than
    * offset by this precaching step
    */
-  void SofQW3::initThetaCache(API::MatrixWorkspace_const_sptr workspace)
+  void SofQW3::initAngularCachesNonPSD(const API::MatrixWorkspace_const_sptr & workspace)
   {
     const size_t nhist = workspace->getNumberHistograms();
     this->m_theta = std::vector<double>(nhist);
-    size_t ndets(0);
-    double minTheta(DBL_MAX), maxTheta(-DBL_MAX);
+    this->m_thetaWidths = std::vector<double>(nhist);
+    // Force phi widths to zero
+    this->m_phi = std::vector<double>(nhist, 0.0);
+    this->m_phiWidths = std::vector<double>(nhist, 0.0);
 
-    for(int64_t i = 0 ; i < (int64_t)nhist; ++i) //signed for OpenMP
-    {
+    auto inst = workspace->getInstrument();
+    const auto samplePos = inst->getSample()->getPos();
+    const PointingAlong upDir = inst->getReferenceFrame()->pointingUp();
 
+    for(size_t i = 0 ; i < nhist; ++i) //signed for OpenMP
+    {
       m_progress->report("Calculating detector angles");
       IDetector_const_sptr det;
       try
@@ -280,25 +283,42 @@ namespace Algorithms
       if( !det || det->isMonitor() )
       {
         this->m_theta[i] = -1.0; // Indicates a detector to skip
+        this->m_thetaWidths[i] = -1.0;
+        continue;
       }
-      else
+      const double theta = workspace->detectorTwoTheta(det);
+      this->m_theta[i] = theta;
+
+      /**
+       * Determine width from shape geometry. A group is assumed to contain
+       * detectors with the same shape & r, theta value, i.e. a ring mapped-group
+       * The shape is retrieved and rotated to match the rotation of the detector.
+       * The angular width is computed using the l2 distance from the sample
+       */
+      if(auto group = boost::dynamic_pointer_cast<const DetectorGroup>(det))
       {
-        ++ndets;
-        const double theta = workspace->detectorTwoTheta(det);
-        this->m_theta[i] = theta;
-        if( theta < minTheta )
-        {
-          minTheta = theta;
-        }
-        else if( theta > maxTheta )
-        {
-          maxTheta = theta;
-        }
+        // assume they all have same shape and same r,theta
+        auto dets = group->getDetectors();
+        det = dets[0];
+      }
+      const auto pos = det->getPos();
+      double l2(0.0),t(0.0),p(0.0);
+      pos.getSpherical(l2,t,p);
+      // Get the shape
+      auto shape = det->shape(); // Defined in its own reference frame with centre at 0,0,0
+      auto rot = det->getRotation();
+      BoundingBox bbox = shape->getBoundingBox();
+      auto maxPoint(bbox.maxPoint());
+      rot.rotate(maxPoint);
+      double boxWidth = maxPoint[upDir];
+
+      m_thetaWidths[i] = std::fabs(2.0*std::atan(boxWidth/l2));
+      if(g_log.is(Logger::Priority::PRIO_DEBUG))
+      {
+        g_log.debug() << "Detector at spectrum =" << workspace->getSpectrum(i)->getSpectrumNo() << ", width=" << m_thetaWidths[i]*180.0/M_PI << " degrees\n";
       }
     }
 
-    this->m_thetaWidth = (maxTheta - minTheta)/static_cast<double>(ndets);
-    g_log.information() << "Calculated detector width in theta=" << (m_thetaWidth*180.0/M_PI) << " degrees.\n";
   }
 
   /**
@@ -307,7 +327,7 @@ namespace Algorithms
    * it calculates the two-theta and azimuthal angle widths.
    * @param workspace : the workspace containing the needed detector information
    */
-  void SofQW3::getValuesAndWidths(API::MatrixWorkspace_const_sptr workspace)
+  void SofQW3::initAngularCachesPSD(const API::MatrixWorkspace_const_sptr &workspace)
   {
     // Trigger a build of the nearst neighbors outside the OpenMP loop
     const int numNeighbours = 4;
@@ -352,20 +372,20 @@ namespace Algorithms
             spec == deltaPlusT || spec == deltaMinusT)
         {
           DetConstPtr detector_n = workspace->getDetector(spec - 1);
-          double theta_n = workspace->detectorTwoTheta(detector_n);
+          double theta_n = workspace->detectorTwoTheta(detector_n)/2.0;
           double phi_n = detector_n->getPhi();
 
           double dTheta = std::fabs(theta - theta_n);
           double dPhi = std::fabs(phi - phi_n);
           if (dTheta > thetaWidth)
           {
             thetaWidth = dTheta;
-            //g_log.debug() << "Current ThetaWidth: " << thetaWidth << std::endl;
+            g_log.information() << "Current ThetaWidth: " << thetaWidth*180/M_PI << std::endl;
           }
           if (dPhi > phiWidth)
           {
             phiWidth = dPhi;
-            //g_log.debug() << "Current PhiWidth: " << phiWidth << std::endl;
+            g_log.information() << "Current PhiWidth: " << phiWidth*180/M_PI << std::endl;
           }
         }
       }

Code/Mantid/Framework/Geometry/inc/MantidGeometry/Math/PolygonEdge.h

@@ -11,7 +11,6 @@ namespace Mantid
 {
 namespace Geometry
 {
-
   /** PolygonEdge 
     Defines a directed edge between two points on a polygon
 
@@ -73,8 +72,8 @@ namespace Geometry
 
   /// Enumeration for point type w.r.t an edge
   enum PointClassification {
-    Left,         /**< Point is to left of edge */
-    Right,        /**< Point is to right of edge */
+    OnLeft,         /**< Point is to left of edge */
+    OnRight,        /**< Point is to right of edge */
     Beyond,       /**< Point is right of edge destination */
     Behind,       /**< Point is left of edge origin */
     Between,      /**< Point is between edge origin and destination */

Code/Mantid/Framework/Geometry/src/Math/ConvexPolygon.cpp

@@ -124,7 +124,7 @@ namespace Mantid
       do
       {
         PolygonEdge edge(v->point(), v->next()->point());
-        if( classify(point, edge) == Left ) return false;
+        if( classify(point, edge) == OnLeft ) return false;
         v = v->next();
       }
       while( v != m_head );

Code/Mantid/Framework/Geometry/src/Math/LaszloIntersection.cpp

@@ -133,11 +133,11 @@ namespace Mantid
               }
             }
           }
-          if (pclass == Right)
+          if (pclass == OnRight)
           {
             inflag = PIsInside;
           }
-          else if (qclass == Right)
+          else if (qclass == OnRight)
           {
             inflag = QIsInside;
           }
@@ -164,7 +164,7 @@ namespace Mantid
 #endif
         if (pAIMSq && qAIMSp)
         {
-          if( (inflag == QIsInside) || ((inflag == Unknown) && (pclass == Left)) )
+          if( (inflag == QIsInside) || ((inflag == Unknown) && (pclass == OnLeft)) )
           {
 #ifdef VERBOSE
             std::cout << "Move edge on P" << std::endl;
@@ -195,7 +195,7 @@ namespace Mantid
         }
         else 
         {
-          if ((inflag == QIsInside) || ((inflag == Unknown) && (pclass == Left)))
+          if ((inflag == QIsInside) || ((inflag == Unknown) && (pclass == OnLeft)))
           {
 #ifdef VERBOSE
             std::cout << "Move edge on P" << std::endl;

Code/Mantid/Framework/Geometry/src/Math/PolygonEdge.cpp

@@ -49,11 +49,11 @@ namespace Mantid
       double sa = a.X() * b.Y() - b.X() * a.Y();
       if(sa > 0.0)
       {
-        return Left;
+        return OnLeft;
       }
       if(sa < 0.0)
       {
-        return Right;
+        return OnRight;
       }
       if((a.X() * b.X() < 0.0) || (a.Y() * b.Y() < 0.0))
       {
@@ -90,7 +90,7 @@ namespace Mantid
       if( Kernel::equals(denom, 0.0) )
       {
         PointClassification edgeClass = classify(focusEdge.start(), refEdge);
-        if( edgeClass == Left || edgeClass == Right )
+        if( edgeClass == OnLeft || edgeClass == OnRight )
         {
           return PolygonEdge::Parallel;
         }
@@ -178,11 +178,11 @@ namespace Mantid
         double cb = vb.X() * va.Y();
         if( Kernel::gtEquals(ca, cb) )
         {
-          return (aclass != Right);
+          return (aclass != OnRight);
         }
         else
         {
-          return (aclass != Left);
+          return (aclass != OnLeft);
         }
       }
       else

Code/Mantid/Framework/Geometry/test/PolygonEdgeTest.h

@@ -69,8 +69,8 @@ class PolygonEdgeTest : public CxxTest::TestSuite
   {
     const PolygonEdge edge(V2D(0.1,0.1), V2D(2.0,2.0));
 
-    TS_ASSERT_EQUALS(classify(V2D(0.05,0.1), edge), Mantid::Geometry::Left);
-    TS_ASSERT_EQUALS(classify(V2D(0.3,0.1), edge), Mantid::Geometry::Right);
+    TS_ASSERT_EQUALS(classify(V2D(0.05,0.1), edge), Mantid::Geometry::OnLeft);
+    TS_ASSERT_EQUALS(classify(V2D(0.3,0.1), edge), Mantid::Geometry::OnRight);
     TS_ASSERT_EQUALS(classify(V2D(-0.05,-0.05), edge), Mantid::Geometry::Behind);
     TS_ASSERT_EQUALS(classify(V2D(2.5,2.5), edge), Mantid::Geometry::Beyond);
     TS_ASSERT_EQUALS(classify(V2D(1.4,1.4), edge), Mantid::Geometry::Between);

Code/Mantid/Framework/PythonInterface/mantid/api/src/Exports/ITableWorkspace.cpp

@@ -30,7 +30,7 @@ namespace
   /// Boost macro for "looping" over builtin types
   #define BUILTIN_TYPES \
     BOOST_PP_TUPLE_TO_LIST( \
-        5, (double, std::string, int, int64_t, float) \
+      6, (double, std::string, int, int64_t, float, size_t)    \
     )
   #define USER_TYPES \
     BOOST_PP_TUPLE_TO_LIST( \