Refs #4998 warn or omit when outer radius is off edge

Code/Mantid/Framework/MDEvents/inc/MantidMDEvents/IntegratePeaksMD.h

@@ -44,6 +44,12 @@ namespace MDEvents
     /// Input MDEventWorkspace
     Mantid::API::IMDEventWorkspace_sptr inWS;
 
+    /// Calculate if this Q is on a detector
+    bool detectorQ(Mantid::Kernel::V3D QLabFrame, double PeakRadius);
+
+    /// Instrument reference
+    Geometry::Instrument_const_sptr inst;
+
   };
 
 

Code/Mantid/Framework/MDEvents/src/IntegratePeaksMD.cpp

@@ -148,6 +148,9 @@ namespace MDEvents
     declareProperty("ReplaceIntensity", true, "Always replace intensity in PeaksWorkspacem (default).\n"
         "If false, then do not replace intensity if calculated value is 0 (used for SNSSingleCrystalReduction)");
 
+    declareProperty("IntegrateIfOnEdge", true, "Only warning if all of peak outer radius is not on detector (default).\n"
+        "If false, do not integrate if the outer radius is not on a detector.");
+
   }
 
   //----------------------------------------------------------------------------------------------
@@ -181,6 +184,7 @@ namespace MDEvents
     double BackgroundInnerRadius = getProperty("BackgroundInnerRadius");
     /// Replace intensity with 0
     bool replaceIntensity = getProperty("ReplaceIntensity");
+    bool integrateEdge = getProperty("IntegrateIfOnEdge");
     if (BackgroundInnerRadius < PeakRadius)
       BackgroundInnerRadius = PeakRadius;
 
@@ -200,6 +204,26 @@ namespace MDEvents
       else if (CoordinatesToUse == "HKL")
         pos = p.getHKL();
 
+      // Get the instrument and its detectors
+      inst = peakWS->getInstrument();
+      // Do not integrate if sphere is off edge of detector
+      if (BackgroundOuterRadius > PeakRadius)
+      {
+    	  if (!detectorQ(p.getQLabFrame(), BackgroundOuterRadius))
+          {
+             g_log.warning() << "Warning: sphere for integration is off edge of detector for peak " << i << std::endl;
+             if (!integrateEdge)continue;
+          }
+      }
+      else
+      {
+    	  if (!detectorQ(p.getQLabFrame(), PeakRadius))
+          {
+             g_log.warning() << "Warning: sphere for integration is off edge of detector for peak " << i << std::endl;
+             if (!integrateEdge)continue;
+          }
+      }
+
       // Build the sphere transformation
       bool dimensionsUsed[nd];
       coord_t center[nd];
@@ -280,6 +304,39 @@ namespace MDEvents
     setProperty("OutputWorkspace", peakWS);
 
   }
+  /** Calculate if this Q is on a detector
+   *
+   * @param QLabFrame of radius of integration
+   */
+  bool IntegratePeaksMD::detectorQ(Mantid::Kernel::V3D QLabFrame, double r)
+  {
+    bool in = true;
+    int nAngles = 8;
+    double dAngles = static_cast<coord_t>(nAngles);
+    // check 64 points in theta and phi at outer radius
+    for (int i=0; i < nAngles; ++i)
+    {
+      double theta = 6.28318531/dAngles * i;
+      for (int j=0; j < nAngles; ++j)
+      {
+    	 double phi = 6.28318531/dAngles * j;
+         V3D edge = V3D(QLabFrame.X()+r*std::cos(theta)*std::sin(phi),
+           QLabFrame.Y()+r*std::sin(theta)*std::sin(phi), QLabFrame.Z()+r*std::cos(phi));
+         // Create the peak using the Q in the lab frame with all its info:
+         try
+         {
+			 Peak p(inst, edge);
+			 in = (in && p.findDetector());
+			 if (!in) return in;
+         }
+         catch (...)
+         {
+        	 return false;
+         }
+      }
+    }
+    return in;
+  }
 
   //----------------------------------------------------------------------------------------------
   /** Execute the algorithm.

Code/Mantid/Framework/MDEvents/test/IntegratePeaksMDTest.h

@@ -48,7 +48,7 @@ class IntegratePeaksMDTest : public CxxTest::TestSuite
   /** Run the IntegratePeaksMD with the given peak radius integration param */
   static void doRun(double PeakRadius, double BackgroundRadius,
       std::string OutputWorkspace = "IntegratePeaksMDTest_peaks",
-      double BackgroundStartRadius = 0.0)
+      double BackgroundStartRadius = 0.0, bool edge = true)
   {
     IntegratePeaksMD alg;
     TS_ASSERT_THROWS_NOTHING( alg.initialize() )
@@ -58,6 +58,7 @@ class IntegratePeaksMDTest : public CxxTest::TestSuite
     TS_ASSERT_THROWS_NOTHING( alg.setProperty("PeakRadius", PeakRadius ) );
     TS_ASSERT_THROWS_NOTHING( alg.setProperty("BackgroundOuterRadius", BackgroundRadius ) );
     TS_ASSERT_THROWS_NOTHING( alg.setProperty("BackgroundInnerRadius", BackgroundStartRadius ) );
+    TS_ASSERT_THROWS_NOTHING( alg.setProperty("IntegrateIfOnEdge", edge ) );
     TS_ASSERT_THROWS_NOTHING( alg.setPropertyValue("PeaksWorkspace", "IntegratePeaksMDTest_peaks" ) );
     TS_ASSERT_THROWS_NOTHING( alg.setPropertyValue("OutputWorkspace", OutputWorkspace) );
     TS_ASSERT_THROWS_NOTHING( alg.execute() );
@@ -108,13 +109,30 @@ class IntegratePeaksMDTest : public CxxTest::TestSuite
     TS_ASSERT_DELTA( mdews->getBox()->getSignal(), 3000.0, 1e-2);
 
     // Make a fake instrument - doesn't matter, we won't use it really
-    Instrument_sptr inst = ComponentCreationHelper::createTestInstrumentCylindrical(5);
+    Instrument_sptr inst = ComponentCreationHelper::createTestInstrumentRectangular(1, 100, 0.05);
+
+    // --- Make a fake PeaksWorkspace ---
+    PeaksWorkspace_sptr peakWS0(new PeaksWorkspace());
+    peakWS0->setInstrument(inst);
+    peakWS0->addPeak( Peak(inst, 15050, 1.0 ) );
+
+    TS_ASSERT_EQUALS( peakWS0->getPeak(0).getIntensity(), 0.0);
+    AnalysisDataService::Instance().add("IntegratePeaksMDTest_peaks",peakWS0);
+
+    // ------------- Integrate with 0.1 radius but IntegrateIfOnEdge false------------------------
+    doRun(0.1,0.0,"IntegratePeaksMDTest_peaks",0.0,false);
+
+    TS_ASSERT_DELTA( peakWS0->getPeak(0).getIntensity(), 2.0, 1e-2);
+
+    // Error is also calculated
+    TS_ASSERT_DELTA( peakWS0->getPeak(0).getSigmaIntensity(), sqrt(2.0), 1e-2);
+    AnalysisDataService::Instance().remove("IntegratePeaksMDTest_peaks");
 
     // --- Make a fake PeaksWorkspace ---
     PeaksWorkspace_sptr peakWS(new PeaksWorkspace());
-    peakWS->addPeak( Peak(inst, 1, 1.0, V3D(0., 0., 0.) ) );
-    peakWS->addPeak( Peak(inst, 1, 1.0, V3D(2., 3., 4.) ) );
-    peakWS->addPeak( Peak(inst, 1, 1.0, V3D(6., 6., 6.) ) );
+    peakWS->addPeak( Peak(inst, 15050, 1.0, V3D(0., 0., 0.) ) );
+    peakWS->addPeak( Peak(inst, 15050, 1.0, V3D(2., 3., 4.) ) );
+    peakWS->addPeak( Peak(inst, 15050, 1.0, V3D(6., 6., 6.) ) );
 
     TS_ASSERT_EQUALS( peakWS->getPeak(0).getIntensity(), 0.0);
     AnalysisDataService::Instance().add("IntegratePeaksMDTest_peaks",peakWS);
@@ -178,9 +196,10 @@ class IntegratePeaksMDTest : public CxxTest::TestSuite
     // These had no bg, so they are the same
     TS_ASSERT_DELTA( peakWS->getPeak(1).getIntensity(), 1000.0, 1e-2);
     TS_ASSERT_DELTA( peakWS->getPeak(2).getIntensity(), 125.0, 10.0);
-
+    
     AnalysisDataService::Instance().remove("IntegratePeaksMDTest_MDEWS");
     AnalysisDataService::Instance().remove("IntegratePeaksMDTest_peaks");
+
   }