Merge pull request #15029 from mantidproject/14857_integrate_edge

Estimate missing intensity of peaks off edge of detector

Framework/MDAlgorithms/inc/MantidMDAlgorithms/IntegratePeaksMD2.h

@@ -48,7 +48,7 @@ private:
 
   /// Calculate if this Q is on a detector
   void calculateE1(Geometry::Instrument_const_sptr inst);
-  bool detectorQ(Mantid::Kernel::V3D QLabFrame, double PeakRadius);
+  double detectorQ(Mantid::Kernel::V3D QLabFrame, double PeakRadius);
   void runMaskDetectors(Mantid::DataObjects::PeaksWorkspace_sptr peakWS,
                         std::string property, std::string values);
 

Framework/MDAlgorithms/src/IntegratePeaksMD2.cpp

@@ -150,6 +150,12 @@ void IntegratePeaksMD2::init() {
                   "PeakRadius plus this value multiplied"
                   "by the magnitude of Q at the peak center so each peak has a "
                   "different integration radius.  Q includes the 2*pi factor.");
+
+  declareProperty(
+      "CorrectIfOnEdge", false,
+      "Only warning if all of peak outer radius is not on detector (default).\n"
+      "If false, correct for volume off edge for both background and "
+      "intensity.");
 }
 
 //----------------------------------------------------------------------------------------------
@@ -257,6 +263,7 @@ void IntegratePeaksMD2::integrate(typename MDEventWorkspace<MDE, nd>::sptr ws) {
   /// Replace intensity with 0
   bool replaceIntensity = getProperty("ReplaceIntensity");
   bool integrateEdge = getProperty("IntegrateIfOnEdge");
+  bool correctEdge = getProperty("CorrectIfOnEdge");
 
   std::string profileFunction = getProperty("ProfileFunction");
   std::string integrationOption = getProperty("IntegrationOption");
@@ -270,6 +277,11 @@ void IntegratePeaksMD2::integrate(typename MDEventWorkspace<MDE, nd>::sptr ws) {
     outFile = save_path + outFile;
     out.open(outFile.c_str(), std::ofstream::out);
   }
+  // volume of Background sphere with inner volume subtracted
+  double volumeBkg = 4.0 / 3.0 * M_PI * (std::pow(BackgroundOuterRadius, 3) -
+                                         std::pow(BackgroundOuterRadius, 3));
+  // volume of PeakRadius sphere
+  double volumeRadius = 4.0 / 3.0 * M_PI * std::pow(PeakRadius, 3);
   //
   // If the following OMP pragma is included, this algorithm seg faults
   // sporadically when processing multiple TOPAZ runs in a script, on
@@ -302,10 +314,12 @@ void IntegratePeaksMD2::integrate(typename MDEventWorkspace<MDE, nd>::sptr ws) {
 
     // Do not integrate if sphere is off edge of detector
 
-    if (!detectorQ(p.getQLabFrame(),
-                   std::max(BackgroundOuterRadius, PeakRadius))) {
+    double edge = detectorQ(p.getQLabFrame(),
+                            std::max(BackgroundOuterRadius, PeakRadius));
+    if (edge < std::max(BackgroundOuterRadius, PeakRadius)) {
       g_log.warning() << "Warning: sphere/cylinder for integration is off edge "
-                         "of detector for peak " << i << std::endl;
+                         "of detector for peak " << i
+                      << "; radius of edge =  " << edge << std::endl;
       if (!integrateEdge) {
         if (replaceIntensity) {
           p.setIntensity(0.0);
@@ -620,10 +634,32 @@ void IntegratePeaksMD2::integrate(typename MDEventWorkspace<MDE, nd>::sptr ws) {
         2.0 * std::max(PeakRadiusVector[i], BackgroundOuterRadiusVector[i]));
     // Save it back in the peak object.
     if (signal != 0. || replaceIntensity) {
-      p.setIntensity(signal - ratio * background_total - bgSignal);
-      p.setSigmaIntensity(sqrt(errorSquared +
-                               ratio * ratio * std::fabs(background_total) +
-                               bgErrorSquared));
+      double edgeMultiplier = 1.0;
+      double peakMultiplier = 1.0;
+      if (correctEdge) {
+        if (edge < BackgroundOuterRadius) {
+          double e1 = BackgroundOuterRadius - edge;
+          // volume of cap of sphere with h = edge
+          double f1 =
+              M_PI * std::pow(e1, 2) / 3 * (3 * BackgroundOuterRadius - e1);
+          edgeMultiplier = volumeBkg / (volumeBkg - f1);
+        }
+        if (edge < PeakRadius) {
+          double sigma = PeakRadius / 3.0;
+          // assume gaussian peak
+          double e1 =
+              std::exp(-std::pow(edge, 2) / (2 * sigma * sigma)) * PeakRadius;
+          // volume of cap of sphere with h = edge
+          double f1 = M_PI * std::pow(e1, 2) / 3 * (3 * PeakRadius - e1);
+          peakMultiplier = volumeRadius / (volumeRadius - f1);
+        }
+      }
+      p.setIntensity(peakMultiplier * signal -
+                     edgeMultiplier * (ratio * background_total + bgSignal));
+      p.setSigmaIntensity(
+          sqrt(peakMultiplier * errorSquared +
+               edgeMultiplier * (ratio * ratio * std::fabs(background_total) +
+                                 bgErrorSquared)));
     }
 
     g_log.information() << "Peak " << i << " at " << pos << ": signal "
@@ -703,27 +739,45 @@ void IntegratePeaksMD2::calculateE1(Geometry::Instrument_const_sptr inst) {
  * @param QLabFrame: The Peak center.
  * @param r: Peak radius.
  */
-bool IntegratePeaksMD2::detectorQ(Mantid::Kernel::V3D QLabFrame, double r) {
-
+double IntegratePeaksMD2::detectorQ(Mantid::Kernel::V3D QLabFrame, double r) {
+  double edge = r;
   for (auto E1 = E1Vec.begin(); E1 != E1Vec.end(); ++E1) {
     V3D distv = QLabFrame -
                 *E1 * (QLabFrame.scalar_prod(
                           *E1)); // distance to the trajectory as a vector
-    if (distv.norm() < r) {
-      return false;
+    if (distv.norm() < std::min(r, edge)) {
+      edge = distv.norm();
     }
   }
-  return true;
+  return edge;
 }
 void IntegratePeaksMD2::runMaskDetectors(
     Mantid::DataObjects::PeaksWorkspace_sptr peakWS, std::string property,
     std::string values) {
-  IAlgorithm_sptr alg = createChildAlgorithm("MaskBTP");
-  alg->setProperty<Workspace_sptr>("Workspace", peakWS);
-  alg->setProperty(property, values);
-  if (!alg->execute())
-    throw std::runtime_error(
-        "MaskDetectors Child Algorithm has not executed successfully");
+  // For CORELLI do not count as edge if next to another detector bank
+  if (property == "Tube" && peakWS->getInstrument()->getName() == "CORELLI") {
+    IAlgorithm_sptr alg = createChildAlgorithm("MaskBTP");
+    alg->setProperty<Workspace_sptr>("Workspace", peakWS);
+    alg->setProperty("Bank", "1,7,12,17,22,27,30,59,63,69,74,79,84,89");
+    alg->setProperty(property, "1");
+    if (!alg->execute())
+      throw std::runtime_error(
+          "MaskDetectors Child Algorithm has not executed successfully");
+    IAlgorithm_sptr alg2 = createChildAlgorithm("MaskBTP");
+    alg2->setProperty<Workspace_sptr>("Workspace", peakWS);
+    alg2->setProperty("Bank", "6,11,16,21,26,29,58,62,68,73,78,83,88,91");
+    alg2->setProperty(property, "16");
+    if (!alg2->execute())
+      throw std::runtime_error(
+          "MaskDetectors Child Algorithm has not executed successfully");
+  } else {
+    IAlgorithm_sptr alg = createChildAlgorithm("MaskBTP");
+    alg->setProperty<Workspace_sptr>("Workspace", peakWS);
+    alg->setProperty(property, values);
+    if (!alg->execute())
+      throw std::runtime_error(
+          "MaskDetectors Child Algorithm has not executed successfully");
+  }
 }
 
 void IntegratePeaksMD2::checkOverlap(

docs/source/algorithms/IntegratePeaksMD-v2.rst

@@ -97,7 +97,7 @@ IntegrateIfOnEdge option
 ###################################
 
 Edges for each bank or pack of tubes of the instrument are defined by masking the edges in the PeaksWorkspace instrument. 
-e.g. For CORELLI, tubes 1 and 16, and pixels 0 and 255.
+e.g. For TOPAZ pixels 0 and 255 in both directions for the Rectangular Detector.
 Q in the lab frame for every peak is calculated, call it C
 For every point on the edge, the trajectory in reciprocal space is a straight line, going through:
 
@@ -124,6 +124,53 @@ for the integration, one of the detector trajectories on the edge is too close t
 This method is also applied to all masked pixels.  If there are masked pixels trajectories inside an integration volume, the peak must be rejected.
 
    
+CorrectIfOnEdge option
+###################################
+
+This is an extension of what was calculated for the IntegrateIfOnEdge option.  It will only be calculated if this option  
+is true and the minimum dv is less than PeakRadius or BackgroundOuterRadius.  
+
+For the background if
+
+:math:`\left|dv\right|_{min}<BackgroundOuterRadius` 
+
+:math:`h = BackgroundOuterRadius - \left|dv\right|_{min}`
+
+From the minimum of dv the volume of the cap of the sphere is found:
+
+:math:`V_{cap} = \pi h^2 / 3 (3 * BackgroundOuterRadius - h)`
+
+The volume of the total sphere is calculated and for the background the volume of the inner radius must be subtracted:
+
+:math:`V_{shell} = 4/3 \pi (BackgroundOuterRadius^3 - BackgroundInnerRadius^3)`
+
+The integrated intensity is multiplied by the ratio of the volume of the sphere divided by the volume where data was collected
+
+:math:`I_{bkgMultiplier} = V_{shell} / (V_{shell} - V_{cap})`
+
+
+For the peak assume that the shape is Gaussian.  If
+
+:math:`\left|dv\right|_{min}<PeakRadius`
+
+:math:`\sigma = PeakRadius / 3`
+
+:math:`h = PeakRadius * exp(-\left|dv\right|_{min}^2 / (2 \sigma^2)`
+
+From the minimum of dv the volume of the cap of the sphere is found:
+
+:math:`V_{cap} = \pi h^2 / 3 (3 * PeakRadius - h)`
+
+and the volume of the sphere is calculated
+
+:math:`V_{sphere} = 4/3 \pi PeakRadius^3`
+
+The integrated intensity is multiplied by the ratio of the volume of the sphere divided by the volume where data was collected
+
+:math:`I_{peakMultiplier} = V_{sphere} / (V_{sphere} - V_{cap})`
+
+
+   
 Usage
 ------