diff --git a/Code/Mantid/Framework/MDAlgorithms/inc/MantidMDAlgorithms/Integrate3DEvents.h b/Code/Mantid/Framework/MDAlgorithms/inc/MantidMDAlgorithms/Integrate3DEvents.h index 0a2623c8aeb1c3145bbf073514cc1d5bcdbfc55f..977aa62e935118570ea96689f16775ded461ccdc 100644 --- a/Code/Mantid/Framework/MDAlgorithms/inc/MantidMDAlgorithms/Integrate3DEvents.h +++ b/Code/Mantid/Framework/MDAlgorithms/inc/MantidMDAlgorithms/Integrate3DEvents.h @@ -70,7 +70,7 @@ public: void addEvents(std::vector<std::pair<double, Mantid::Kernel::V3D> > const &event_qs, bool hkl_integ); /// Find the net integrated intensity of a peak, using ellipsoidal volumes - boost::shared_ptr<const Mantid::Geometry::PeakShape> ellipseIntegrateEvents(Mantid::Kernel::V3D const &peak_q, bool specify_size, + boost::shared_ptr<const Mantid::Geometry::PeakShape> ellipseIntegrateEvents(std::vector<Kernel::V3D> E1Vec, Mantid::Kernel::V3D const &peak_q, bool specify_size, double peak_radius, double back_inner_radius, double back_outer_radius, std::vector<double> &axes_radii, double &inti, @@ -82,6 +82,12 @@ private: std::vector<Mantid::Kernel::V3D> const &directions, std::vector<double> const &sizes); + /// Calculate the number of events in an ellipsoid centered at 0,0,0 + static double numInEllipsoidBkg(std::vector<std::pair<double, Mantid::Kernel::V3D> > const &events, + std::vector<Mantid::Kernel::V3D> const &directions, + std::vector<double> const &sizes, + std::vector<double> const &sizesIn); + /// Calculate the 3x3 covariance matrix of a list of Q-vectors at 0,0,0 static void makeCovarianceMatrix(std::vector<std::pair<double, Mantid::Kernel::V3D> > const &events, Kernel::DblMatrix &matrix, double radius); @@ -106,11 +112,11 @@ private: /// Find the net integrated intensity of a list of Q's using ellipsoids boost::shared_ptr<const Mantid::DataObjects::PeakShapeEllipsoid> ellipseIntegrateEvents( - std::vector<std::pair<double, Mantid::Kernel::V3D> > const &ev_list, std::vector<Mantid::Kernel::V3D> const &directions, + std::vector<Kernel::V3D> E1Vec, Kernel::V3D const &peak_q, std::vector<std::pair<double, Mantid::Kernel::V3D> > const &ev_list, std::vector<Mantid::Kernel::V3D> const &directions, std::vector<double> const &sigmas, bool specify_size, double peak_radius, double back_inner_radius, double back_outer_radius, std::vector<double> &axes_radii, double &inti, double &sigi); - + double detectorQ(std::vector<Kernel::V3D> E1Vec, const Mantid::Kernel::V3D QLabFrame, std::vector<double>& r); // Private data members PeakQMap peak_qs; // hashtable with peak Q-vectors EventListMap event_lists; // hashtable with lists of events for each peak diff --git a/Code/Mantid/Framework/MDAlgorithms/inc/MantidMDAlgorithms/IntegrateEllipsoids.h b/Code/Mantid/Framework/MDAlgorithms/inc/MantidMDAlgorithms/IntegrateEllipsoids.h index e391ac2a93906cd0509a535d0b95df17723bb736..6dedaa0affec0ffc880276cb1a1596e4f623ae3a 100644 --- a/Code/Mantid/Framework/MDAlgorithms/inc/MantidMDAlgorithms/IntegrateEllipsoids.h +++ b/Code/Mantid/Framework/MDAlgorithms/inc/MantidMDAlgorithms/IntegrateEllipsoids.h @@ -9,6 +9,7 @@ #include "MantidMDAlgorithms/MDTransfInterface.h" #include "MantidDataObjects/Workspace2D.h" #include "MantidDataObjects/EventWorkspace.h" +#include "MantidDataObjects/PeaksWorkspace.h" namespace Mantid { namespace MDAlgorithms { @@ -37,6 +38,15 @@ private: DataObjects::Workspace2D_sptr &wksp, Kernel::DblMatrix const &UBinv, bool hkl_integ); + /// Calculate if this Q is on a detector + void calculateE1(Geometry::Instrument_const_sptr inst) ; + + void runMaskDetectors(Mantid::DataObjects::PeaksWorkspace_sptr peakWS, + std::string property, std::string values); + + /// save for all detector pixels + std::vector<Kernel::V3D> E1Vec; + MDWSDescription m_targWSDescr; void initTargetWSDescr(API::MatrixWorkspace_sptr &wksp); diff --git a/Code/Mantid/Framework/MDAlgorithms/src/Integrate3DEvents.cpp b/Code/Mantid/Framework/MDAlgorithms/src/Integrate3DEvents.cpp index 9b316b490e337cd364c951ffd09c16601f4b3342..eaeb223c84367dc734e61c6a8445b4fb4f92d4cf 100644 --- a/Code/Mantid/Framework/MDAlgorithms/src/Integrate3DEvents.cpp +++ b/Code/Mantid/Framework/MDAlgorithms/src/Integrate3DEvents.cpp @@ -84,6 +84,7 @@ void Integrate3DEvents::addEvents(std::vector<std::pair<double, V3D> > const &ev * half the major axis length of the ellipsoids, and the other axes of the * ellipsoids will be set proportionally, based on the standard deviations. * + * @param E1Vec Vector of values for calculating edge of detectors * @param peak_q The Q-vector for the peak center. * @param specify_size If true the integration will be done using the * ellipsoids with major axes determined by the @@ -111,7 +112,7 @@ void Integrate3DEvents::addEvents(std::vector<std::pair<double, V3D> > const &ev * */ Mantid::Geometry::PeakShape_const_sptr Integrate3DEvents::ellipseIntegrateEvents( - V3D const &peak_q, bool specify_size, double peak_radius, + std::vector<Kernel::V3D> E1Vec, V3D const &peak_q, bool specify_size, double peak_radius, double back_inner_radius, double back_outer_radius, std::vector<double> &axes_radii, double &inti, double &sigi) { inti = 0.0; // default values, in case something @@ -154,7 +155,7 @@ Mantid::Geometry::PeakShape_const_sptr Integrate3DEvents::ellipseIntegrateEvents return boost::make_shared<NoShape>(); // ellipsoids will be zero. } - return ellipseIntegrateEvents(some_events, eigen_vectors, sigmas, specify_size, + return ellipseIntegrateEvents(E1Vec, peak_q, some_events, eigen_vectors, sigmas, specify_size, peak_radius, back_inner_radius, back_outer_radius, axes_radii, inti, sigi); } @@ -188,7 +189,47 @@ double Integrate3DEvents::numInEllipsoid(std::vector<std::pair<double, V3D>> con return count; } +/** + * Calculate the number of events in an ellipsoid centered at 0,0,0 with + * the three specified axes and the three specified sizes in the direction + * of those axes. NOTE: The three axes must be mutually orthogonal unit + * vectors. + * + * @param events List of 3D events centered at 0,0,0 + * @param directions List of 3 orthonormal directions for the axes of + * the ellipsoid. + * @param sizes List of three values a,b,c giving half the length + * of the three axes of the ellisoid. + * @param sizesIn List of three values a,b,c giving half the length + * of the three inner axes of the ellisoid. + * @return Then number of events that are in or on the specified ellipsoid. + */ +double Integrate3DEvents::numInEllipsoidBkg(std::vector<std::pair<double, V3D>> const &events, + std::vector<V3D> const &directions, + std::vector<double> const &sizes, + std::vector<double> const &sizesIn) { + double count = 0; + std::vector<double> eventVec; + for (size_t i = 0; i < events.size(); i++) { + double sum = 0; + double sumIn = 0; + for (size_t k = 0; k < 3; k++) { + double comp = events[i].second.scalar_prod(directions[k]) / sizes[k]; + sum += comp * comp; + comp = events[i].second.scalar_prod(directions[k]) / sizesIn[k]; + sumIn += comp * comp; + } + if (sum <= 1 && sumIn >= 1) + eventVec.push_back(events[i].first); + } + std::sort(eventVec.begin(), eventVec.end()); + // Remove top 1% of background + for (size_t k = 0; k < static_cast<size_t>(0.99 * static_cast<double>(eventVec.size())); k++) { + count += eventVec[k]; + } + return count; +} /** * Given a list of events, associated with a particular peak * and already SHIFTED to be centered at (0,0,0), calculate the 3x3 @@ -387,6 +428,8 @@ void Integrate3DEvents::addEvent(std::pair<double, V3D> event_Q, bool hkl_integ) * axes for the events and the standard deviations in the the directions * of the principal axes. * + * @param E1Vec Vector of values for calculating edge of detectors + * @param peak_q The Q-vector for the peak center. * @param ev_list List of events centered at (0,0,0) for a * particular peak. * @param directions The three principal axes of the list of events @@ -418,7 +461,7 @@ void Integrate3DEvents::addEvent(std::pair<double, V3D> event_Q, bool hkl_integ) * */ PeakShapeEllipsoid_const_sptr Integrate3DEvents::ellipseIntegrateEvents( - std::vector<std::pair<double, Mantid::Kernel::V3D> > const &ev_list, + std::vector<Kernel::V3D> E1Vec, V3D const &peak_q, std::vector<std::pair<double, Mantid::Kernel::V3D> > const &ev_list, std::vector<Mantid::Kernel::V3D> const &directions, std::vector<double> const &sigmas, bool specify_size, double peak_radius, double back_inner_radius, double back_outer_radius, @@ -460,25 +503,35 @@ PeakShapeEllipsoid_const_sptr Integrate3DEvents::ellipseIntegrateEvents( } axes_radii.clear(); + std::vector<double> abcBackgroundOuterRadii; + std::vector<double> abcBackgroundInnerRadii; + std::vector<double> abcRadii; for (int i = 0; i < 3; i++) { - axes_radii.push_back(r3 * sigmas[i]); + abcBackgroundOuterRadii.push_back(r3 * sigmas[i]); + abcBackgroundInnerRadii.push_back(r2 * sigmas[i]); + abcRadii.push_back(r1 * sigmas[i]); + axes_radii.push_back(r1 * sigmas[i]); } - auto abcBackgroundOuterRadii = axes_radii; - double back2 = numInEllipsoid(ev_list, directions, axes_radii); - for (int i = 0; i < 3; i++) { - axes_radii[i] = r2 * sigmas[i]; + if(E1Vec.size() > 0){ + double h3 = 1.0 - detectorQ(E1Vec, peak_q, abcBackgroundOuterRadii); + // scaled from area of circle minus segment when r normalized to 1 + double m3= std::sqrt(1.0 - (std::acos(1.0-h3)-(1.0-h3)*std::sqrt(2.0*h3-h3*h3))/M_PI); + double h1 = 1.0 - detectorQ(E1Vec, peak_q, axes_radii); + // Do not use peak if edge of detector is inside integration radius + if (h1 > 0.0) return boost::make_shared<const PeakShapeEllipsoid>(directions, abcRadii, abcBackgroundInnerRadii, abcBackgroundOuterRadii, Mantid::Kernel::QLab, "IntegrateEllipsoids"); + r3 *= m3; + if (r2 != r1) { + double h2 = 1.0 - detectorQ(E1Vec, peak_q, abcBackgroundInnerRadii); + // scaled from area of circle minus segment when r normalized to 1 + double m2= std::sqrt(1.0 - (std::acos(1.0-h2)-(1.0-h2)*std::sqrt(2.0*h2-h2*h2))/M_PI); + r2 *= m2; + } } - auto abcBackgroundInnerRadii = axes_radii; - double back1 = numInEllipsoid(ev_list, directions, axes_radii); - for (int i = 0; i < 3; i++) { - axes_radii[i] = r1 * sigmas[i]; - } - auto abcRadii = axes_radii; - double peak_w_back = numInEllipsoid(ev_list, directions, axes_radii); + double backgrd = numInEllipsoidBkg(ev_list, directions, abcBackgroundOuterRadii, abcBackgroundInnerRadii); - double backgrd = back2 - back1; + double peak_w_back = numInEllipsoid(ev_list, directions, axes_radii); double ratio = pow(r1, 3) / (pow(r3, 3) - pow(r2, 3)); @@ -488,7 +541,28 @@ PeakShapeEllipsoid_const_sptr Integrate3DEvents::ellipseIntegrateEvents( // Make the shape and return it. return boost::make_shared<const PeakShapeEllipsoid>(directions, abcRadii, abcBackgroundInnerRadii, abcBackgroundOuterRadii, Mantid::Kernel::QLab, "IntegrateEllipsoids"); } - +/** Calculate if this Q is on a detector + * The distance from C to OE is given by dv=C-E*(C.scalar_prod(E)) + * If dv.norm<integration_radius, one of the detector trajectories on the edge + *is too close to the peak + * This method is applied to all masked pixels. If there are masked pixels + *trajectories inside an integration volume, the peak must be rejected. + * + * @param E1Vec Vector of values for calculating edge of detectors + * @param QLabFrame: The Peak center. + * @param r: Peak radius. + */ +double Integrate3DEvents::detectorQ(std::vector<Kernel::V3D> E1Vec, const Mantid::Kernel::V3D QLabFrame, std::vector<double>& r) { + double quot = 1.0; + for (auto E1 = E1Vec.begin(); E1 != E1Vec.end(); ++E1) { + V3D distv = QLabFrame - *E1 * (QLabFrame.scalar_prod(*E1)); // distance to the trajectory as a vector + double quot0 = distv.norm() / *( std::min_element(r.begin(), r.end())); + if (quot0 < quot) { + quot = quot0; + } + } +return quot; +} } // namespace MDAlgorithms } // namespace Mantid diff --git a/Code/Mantid/Framework/MDAlgorithms/src/IntegrateEllipsoids.cpp b/Code/Mantid/Framework/MDAlgorithms/src/IntegrateEllipsoids.cpp index 2b93b9dd1bbf12f5c9b703021e89f42867cf922a..1ae713b3f761cb23f1a8967b5b512118ec630e3c 100644 --- a/Code/Mantid/Framework/MDAlgorithms/src/IntegrateEllipsoids.cpp +++ b/Code/Mantid/Framework/MDAlgorithms/src/IntegrateEllipsoids.cpp @@ -67,6 +67,7 @@ void IntegrateEllipsoids::qListFromEventWS(Integrate3DEvents &integrator, Progre EventList &events = wksp->getEventList(i); events.switchTo(WEIGHTED_NOTIME); + events.compressEvents(1e-5, &events); // check to see if the event list is empty if (events.empty()) { @@ -275,6 +276,11 @@ void IntegrateEllipsoids::init() { declareProperty( "IntegrateInHKL", false, "If true, integrate in HKL space not Q space."); + + declareProperty( + "IntegrateIfOnEdge", true, + "Set to false to not integrate if peak radius is off edge of detector." + "Background will be scaled if background radius is off edge."); } //--------------------------------------------------------------------- @@ -303,11 +309,6 @@ void IntegrateEllipsoids::exec() { throw std::runtime_error("Could not read the peaks workspace"); } - Mantid::DataObjects::PeaksWorkspace_sptr peak_ws = - getProperty("OutputWorkspace"); - if (peak_ws != in_peak_ws) { - peak_ws.reset(in_peak_ws->clone().release()); - } double radius = getProperty("RegionRadius"); int numSigmas = getProperty("NumSigmas"); double cutoffIsigI = getProperty("CutoffIsigI"); @@ -316,6 +317,28 @@ void IntegrateEllipsoids::exec() { double back_inner_radius = getProperty("BackgroundInnerSize"); double back_outer_radius = getProperty("BackgroundOuterSize"); bool hkl_integ = getProperty("IntegrateInHKL"); + bool integrateEdge = getProperty("IntegrateIfOnEdge"); + if (!integrateEdge) + { + // This only fails in the unit tests which say that MaskBTP is not registered + try { + runMaskDetectors(in_peak_ws, "Tube", "edges"); + runMaskDetectors(in_peak_ws, "Pixel", "edges"); + } catch (...) { + g_log.error("Can't execute MaskBTP algorithm for this instrument to set " + "edge for IntegrateIfOnEdge option"); + } + // Get the instrument and its detectors + Geometry::Instrument_const_sptr inst = in_peak_ws->getInstrument(); + calculateE1(inst); //fill E1Vec for use in detectorQ + } + + Mantid::DataObjects::PeaksWorkspace_sptr peak_ws = + getProperty("OutputWorkspace"); + if (peak_ws != in_peak_ws) { + peak_ws.reset(in_peak_ws->clone().release()); + } + // get UBinv and the list of // peak Q's for the integrator std::vector<Peak> &peaks = peak_ws->getPeaks(); @@ -399,7 +422,7 @@ void IntegrateEllipsoids::exec() { peak_q = peaks[i].getQLabFrame(); std::vector<double> axes_radii; Mantid::Geometry::PeakShape_const_sptr shape = - integrator.ellipseIntegrateEvents( + integrator.ellipseIntegrateEvents(E1Vec, peak_q, specify_size, peak_radius, back_inner_radius, back_outer_radius, axes_radii, inti, sigi); peaks[i].setIntensity(inti); @@ -474,7 +497,7 @@ void IntegrateEllipsoids::exec() { if (Geometry::IndexingUtils::ValidIndex(hkl, 1.0)) { peak_q = peaks[i].getQLabFrame(); std::vector<double> axes_radii; - integrator.ellipseIntegrateEvents( + integrator.ellipseIntegrateEvents(E1Vec, peak_q, specify_size, peak_radius, back_inner_radius, back_outer_radius, axes_radii, inti, sigi); peaks[i].setIntensity(inti); @@ -554,5 +577,45 @@ void IntegrateEllipsoids::initTargetWSDescr(MatrixWorkspace_sptr &wksp) { m_targWSDescr.m_PreprDetTable = table; } +/* + * Define edges for each instrument by masking. For CORELLI, tubes 1 and 16, and + *pixels 0 and 255. + * Get Q in the lab frame for every peak, call it C + * For every point on the edge, the trajectory in reciprocal space is a straight + *line, going through O=V3D(0,0,0). + * Calculate a point at a fixed momentum, say k=1. Q in the lab frame + *E=V3D(-k*sin(tt)*cos(ph),-k*sin(tt)*sin(ph),k-k*cos(ph)). + * Normalize E to 1: E=E*(1./E.norm()) + * + * @param inst: instrument + */ +void IntegrateEllipsoids::calculateE1(Geometry::Instrument_const_sptr inst) { + std::vector<detid_t> detectorIDs = inst->getDetectorIDs(); + + for (auto detID = detectorIDs.begin(); detID != detectorIDs.end(); ++detID) { + Mantid::Geometry::IDetector_const_sptr det = inst->getDetector(*detID); + if (det->isMonitor()) + continue; // skip monitor + if (!det->isMasked()) + continue; // edge is masked so don't check if not masked + double tt1 = det->getTwoTheta(V3D(0, 0, 0), V3D(0, 0, 1)); // two theta + double ph1 = det->getPhi(); // phi + V3D E1 = V3D(-std::sin(tt1) * std::cos(ph1), -std::sin(tt1) * std::sin(ph1), + 1. - std::cos(tt1)); // end of trajectory + E1 = E1 * (1. / E1.norm()); // normalize + E1Vec.push_back(E1); + } + } + +void IntegrateEllipsoids::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"); +} } // namespace MDAlgorithms } // namespace Mantid diff --git a/Code/Mantid/Framework/MDAlgorithms/test/Integrate3DEventsTest.h b/Code/Mantid/Framework/MDAlgorithms/test/Integrate3DEventsTest.h index 7bdbc2900e398d299ff41e765e1cfc8252e4a103..2fd9f4c2c8a52227eb28fb6e1355331663e37927 100644 --- a/Code/Mantid/Framework/MDAlgorithms/test/Integrate3DEventsTest.h +++ b/Code/Mantid/Framework/MDAlgorithms/test/Integrate3DEventsTest.h @@ -27,8 +27,8 @@ public: double inti_all[] = { 755, 704, 603 }; double sigi_all[] = { 27.4773, 26.533, 24.5561 }; - double inti_some[] = { 691, 648, 603 }; - double sigi_some[] = { 27.4773, 26.533, 24.5561 }; + double inti_some[] = { 692, 649, 603 }; + double sigi_some[] = { 27.4590, 26.5141, 24.5561 }; // synthesize three peaks std::vector<std::pair<double, V3D> > peak_q_list; @@ -90,11 +90,12 @@ public: double back_inner_radius = 1.2; double back_outer_radius = 1.3; std::vector<double> new_sigma; + std::vector<Kernel::V3D> E1Vec; double inti; double sigi; for ( size_t i = 0; i < peak_q_list.size(); i++ ) { - auto shape = integrator.ellipseIntegrateEvents( peak_q_list[i].second, specify_size, + auto shape = integrator.ellipseIntegrateEvents( E1Vec, peak_q_list[i].second, specify_size, peak_radius, back_inner_radius, back_outer_radius, new_sigma, inti, sigi ); TS_ASSERT_DELTA( inti, inti_all[i], 0.1); @@ -110,7 +111,7 @@ public: specify_size = false; for ( size_t i = 0; i < peak_q_list.size(); i++ ) { - integrator.ellipseIntegrateEvents( peak_q_list[i].second, specify_size, + integrator.ellipseIntegrateEvents( E1Vec, peak_q_list[i].second, specify_size, peak_radius, back_inner_radius, back_outer_radius, new_sigma, inti, sigi ); TS_ASSERT_DELTA( inti, inti_some[i], 0.1); diff --git a/Code/Mantid/Framework/MDAlgorithms/test/IntegrateEllipsoidsTest.h b/Code/Mantid/Framework/MDAlgorithms/test/IntegrateEllipsoidsTest.h index 108ca61408f6b54514822dfa03e4f1354d38c9d0..5f62d987f9ea4662cdfc181886ad72326730d258 100644 --- a/Code/Mantid/Framework/MDAlgorithms/test/IntegrateEllipsoidsTest.h +++ b/Code/Mantid/Framework/MDAlgorithms/test/IntegrateEllipsoidsTest.h @@ -253,17 +253,17 @@ public: m_peaksWS->getNumberPeaks()); TSM_ASSERT_DELTA("Wrong intensity for peak 0", - integratedPeaksWS->getPeak(0).getIntensity(), -2, 0.01); + integratedPeaksWS->getPeak(0).getIntensity(), -1, 0.01); TSM_ASSERT_DELTA("Wrong intensity for peak 1", - integratedPeaksWS->getPeak(1).getIntensity(), 2, 0.01); + integratedPeaksWS->getPeak(1).getIntensity(), 3, 0.01); TSM_ASSERT_DELTA("Wrong intensity for peak 2", - integratedPeaksWS->getPeak(2).getIntensity(), -2, 0.01); + integratedPeaksWS->getPeak(2).getIntensity(), -1, 0.01); TSM_ASSERT_DELTA("Wrong intensity for peak 3", integratedPeaksWS->getPeak(3).getIntensity(), 15, 0.01); TSM_ASSERT_DELTA("Wrong intensity for peak 4", - integratedPeaksWS->getPeak(4).getIntensity(), 11, 0.01); + integratedPeaksWS->getPeak(4).getIntensity(), 12, 0.01); TSM_ASSERT_DELTA("Wrong intensity for peak 5", - integratedPeaksWS->getPeak(5).getIntensity(), 10, 0.01); + integratedPeaksWS->getPeak(5).getIntensity(), 11, 0.01); } @@ -283,17 +283,17 @@ public: integratedPeaksWS->getNumberPeaks(), m_peaksWS->getNumberPeaks()); TSM_ASSERT_DELTA("Wrong intensity for peak 0", - integratedPeaksWS->getPeak(0).getIntensity(), 1.0, 0.01); + integratedPeaksWS->getPeak(0).getIntensity(), 2.0, 0.01); TSM_ASSERT_DELTA("Wrong intensity for peak 1", - integratedPeaksWS->getPeak(1).getIntensity(), 0, 0.01); + integratedPeaksWS->getPeak(1).getIntensity(), 1.0, 0.01); TSM_ASSERT_DELTA("Wrong intensity for peak 2", - integratedPeaksWS->getPeak(2).getIntensity(), 1.0, 0.01); + integratedPeaksWS->getPeak(2).getIntensity(), 2.0, 0.01); TSM_ASSERT_DELTA("Wrong intensity for peak 3", - integratedPeaksWS->getPeak(3).getIntensity(), 10, 0.01); + integratedPeaksWS->getPeak(3).getIntensity(), 11, 0.01); TSM_ASSERT_DELTA("Wrong intensity for peak 4", - integratedPeaksWS->getPeak(4).getIntensity(), 13, 0.01); + integratedPeaksWS->getPeak(4).getIntensity(), 14, 0.01); TSM_ASSERT_DELTA("Wrong intensity for peak 5", - integratedPeaksWS->getPeak(5).getIntensity(), 12, 0.01); + integratedPeaksWS->getPeak(5).getIntensity(), 13, 0.01); } }; diff --git a/Code/Mantid/Testing/SystemTests/tests/analysis/EllipsoidIntegr.py b/Code/Mantid/Testing/SystemTests/tests/analysis/EllipsoidIntegr.py index 21b008adba94b701e1f2ccfbfcf0ee3d9a30c44e..38271684cb7604b016c06a97528b4b1227505ee6 100644 --- a/Code/Mantid/Testing/SystemTests/tests/analysis/EllipsoidIntegr.py +++ b/Code/Mantid/Testing/SystemTests/tests/analysis/EllipsoidIntegr.py @@ -18,11 +18,11 @@ class EllipsoidIntegr( stresstesting.MantidStressTest): def runTest(self): # expected results with size determined # automatically from projected event sigmas - inti_auto = [ 88, 99, 23, 33, 8, 8, 4 ] + inti_auto = [ 89, 101, 24, 34, 9, 9, 5 ] sigi_auto = [ 13.784, 18.1384, 13.1529, 9.94987, 5.83095, 10.2956, 10.2956] # expected results with fixed size # ellipsoids - inti_fixed = [ 87.541, 95.3934, 21.3607, 33.4262, 7.36066, 9.68852, 3.54098 ] + inti_fixed = [ 88.5902, 97.4918, 22.4098, 34.4754, 8.40984, 10.7377, 4.59016 ] sigi_fixed = [ 13.9656, 18.4523, 13.4335, 10.1106, 5.94223, 10.5231, 10.5375 ] # first, load peaks into a peaks workspace diff --git a/Code/Mantid/docs/source/algorithms/IntegrateEllipsoids-v1.rst b/Code/Mantid/docs/source/algorithms/IntegrateEllipsoids-v1.rst index 6c9b414ca4608dc14e753e00b664a578bc905f07..3befdb87212c9a5c39a02be24d6d0eca022ccd6a 100644 --- a/Code/Mantid/docs/source/algorithms/IntegrateEllipsoids-v1.rst +++ b/Code/Mantid/docs/source/algorithms/IntegrateEllipsoids-v1.rst @@ -165,6 +165,43 @@ surface. This algorithm uses principle component analysis to determine the principle axis for each peak. For the event list (QLab) associated with each peak, the algorithm determines a covariance matrix, and uses that to establish eigenvectors corresponding to the principle axis (all orthogonal). The sizes of each principle axis are used define the region of which events will be counted/integrated from those already associated with each peak. +IntegrateIfOnEdge=False 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. +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: + +:math:`\vec{O}=(0,0,0)` + +Calculate a point at a fixed momentum, say k=1. +Q in the lab frame: + +:math:`\vec{E}=(-k*sin(\theta)*cos(\phi),-k*sin(\theta)*sin(\phi),k-k*cos(\phi))` + +Normalize E to 1: + +:math:`\vec{E}=\vec{E}*(1./\left|\vec{E}\right|)` + +The distance from C to OE is given by: + +:math:`dv=\vec{C}-\vec{E}*(\vec{C} \cdot \vec{E})` + +If: + +:math:`\left|dv\right|<PeakRadius` + +for the integration, one of the detector trajectories on the edge is too close to the peak +This method is also applied to all masked pixels. If there are masked pixels trajectories inside an integration volume, the peak must be rejected. +If there are masked pixel trajectories inside the background volume, the background events are scaled by estimating the volume of the ellipsoid +on the detector. + +Sigma from the background +################################### +The sigma from the background could be too small because the background contains events from other peaks. In an effort to reduce this, all the +background events are sorted and the top 1% are removed. + Usage ------