diff --git a/Code/Mantid/Framework/MDEvents/inc/MantidMDEvents/Integrate3DEvents.h b/Code/Mantid/Framework/MDEvents/inc/MantidMDEvents/Integrate3DEvents.h index 473b1b8103e1..d70da28336da 100644 --- a/Code/Mantid/Framework/MDEvents/inc/MantidMDEvents/Integrate3DEvents.h +++ b/Code/Mantid/Framework/MDEvents/inc/MantidMDEvents/Integrate3DEvents.h @@ -55,19 +55,19 @@ namespace MDEvents { */ typedef Mantid::Kernel::Matrix DblMatrix; -typedef boost::unordered_map> EventListMap; +typedef boost::unordered_map > > EventListMap; typedef boost::unordered_map PeakQMap; class DLLExport Integrate3DEvents { public: /// Construct object to store events around peaks and integrate peaks - Integrate3DEvents(std::vector const &peak_q_list, DblMatrix const &UBinv, + Integrate3DEvents(std::vector > const &peak_q_list, DblMatrix const &UBinv, double radius); ~Integrate3DEvents(); /// Add event Q's to lists of events near peaks - void addEvents(std::vector const &event_qs); + void addEvents(std::vector > const &event_qs); /// Find the net integrated intensity of a peak, using ellipsoidal volumes boost::shared_ptr ellipseIntegrateEvents(Mantid::Kernel::V3D const &peak_q, bool specify_size, @@ -78,12 +78,12 @@ class DLLExport Integrate3DEvents { private: /// Calculate the number of events in an ellipsoid centered at 0,0,0 - static int numInEllipsoid(std::vector const &events, + static double numInEllipsoid(std::vector > const &events, std::vector const &directions, std::vector const &sizes); /// Calculate the 3x3 covariance matrix of a list of Q-vectors at 0,0,0 - static void makeCovarianceMatrix(std::vector const &events, + static void makeCovarianceMatrix(std::vector > const &events, DblMatrix &matrix, double radius); /// Calculate the eigen vectors of a 3x3 real symmetric matrix @@ -91,7 +91,7 @@ class DLLExport Integrate3DEvents { std::vector &eigen_vectors); /// Calculate the standard deviation of 3D events in a specified direction - static double stdDev(std::vector const &events, Mantid::Kernel::V3D const &direction, + static double stdDev(std::vector > const &events, Mantid::Kernel::V3D const &direction, double radius); /// Form a map key as 10^12*h + 10^6*k + l from the integers h, k, l @@ -101,11 +101,11 @@ class DLLExport Integrate3DEvents { int64_t getHklKey(Mantid::Kernel::V3D const &q_vector); /// Add an event to the vector of events for the closest h,k,l - void addEvent(Mantid::Kernel::V3D event_Q); + void addEvent(std::pair event_Q); /// Find the net integrated intensity of a list of Q's using ellipsoids boost::shared_ptr ellipseIntegrateEvents( - std::vector const &ev_list, std::vector const &directions, + std::vector > const &ev_list, std::vector const &directions, std::vector const &sigmas, bool specify_size, double peak_radius, double back_inner_radius, double back_outer_radius, std::vector &axes_radii, double &inti, double &sigi); diff --git a/Code/Mantid/Framework/MDEvents/src/Integrate3DEvents.cpp b/Code/Mantid/Framework/MDEvents/src/Integrate3DEvents.cpp index acdee6cc6b26..be5badfbf2b9 100644 --- a/Code/Mantid/Framework/MDEvents/src/Integrate3DEvents.cpp +++ b/Code/Mantid/Framework/MDEvents/src/Integrate3DEvents.cpp @@ -33,16 +33,16 @@ using Mantid::Kernel::V3D; * an event to be stored in the list associated with * that peak. */ -Integrate3DEvents::Integrate3DEvents(std::vector const &peak_q_list, +Integrate3DEvents::Integrate3DEvents(std::vector> const &peak_q_list, DblMatrix const &UBinv, double radius) { this->UBinv = UBinv; this->radius = radius; int64_t hkl_key; - for (auto it = peak_q_list.begin(); it != peak_q_list.end(); ++it) { - hkl_key = getHklKey(*it); + for (size_t it = 0; it != peak_q_list.size(); ++it) { + hkl_key = getHklKey(peak_q_list[it].second); if (hkl_key != 0) // only save if hkl != (0,0,0) - peak_qs[hkl_key] = *it; + peak_qs[hkl_key] = peak_q_list[it].second; } } @@ -67,7 +67,7 @@ Integrate3DEvents::~Integrate3DEvents() {} * @param event_qs List of event Q vectors to add to lists of Q's associated * with peaks. */ -void Integrate3DEvents::addEvents(std::vector const &event_qs) { +void Integrate3DEvents::addEvents(std::vector > const &event_qs) { for (size_t i = 0; i < event_qs.size(); i++) { addEvent(event_qs[i]); } @@ -120,7 +120,12 @@ Mantid::Geometry::PeakShape_const_sptr Integrate3DEvents::ellipseIntegrateEvents return boost::make_shared(); } - std::vector &some_events = event_lists[hkl_key]; + std::vector > &some_events = event_lists[hkl_key]; + for (size_t it = 0; it != some_events.size(); ++it) { + hkl_key = getHklKey(some_events[it].second); + if (hkl_key != 0) // only save if hkl != (0,0,0) + peak_qs[hkl_key] = some_events[it].second; + } if (some_events.size() < 3) // if there are not enough events to { // find covariance matrix, return @@ -170,18 +175,18 @@ Mantid::Geometry::PeakShape_const_sptr Integrate3DEvents::ellipseIntegrateEvents * of the three axes of the ellisoid. * @return Then number of events that are in or on the specified ellipsoid. */ -int Integrate3DEvents::numInEllipsoid(std::vector const &events, +double Integrate3DEvents::numInEllipsoid(std::vector> const &events, std::vector const &directions, std::vector const &sizes) { - int count = 0; + double count = 0; for (size_t i = 0; i < events.size(); i++) { double sum = 0; for (size_t k = 0; k < 3; k++) { - double comp = events[i].scalar_prod(directions[k]) / sizes[k]; + double comp = events[i].second.scalar_prod(directions[k]) / sizes[k]; sum += comp * comp; } if (sum <= 1) - count++; + count += events[i].first; } return count; @@ -209,14 +214,16 @@ int Integrate3DEvents::numInEllipsoid(std::vector const &events, * peak center (0,0,0) will be used for * calculating the covariance matrix. */ -void Integrate3DEvents::makeCovarianceMatrix(std::vector const &events, + +void Integrate3DEvents::makeCovarianceMatrix(std::vector > const &events, DblMatrix &matrix, double radius) { for (int row = 0; row < 3; row++) { for (int col = 0; col < 3; col++) { double sum = 0; for (size_t i = 0; i < events.size(); i++) { - if (events[i].norm() <= radius) { - sum += events[i][row] * events[i][col]; + auto event = events[i].second; + if (event.norm() <= radius) { + sum += event[row] * event[col]; } } if (events.size() > 1) @@ -275,7 +282,7 @@ void Integrate3DEvents::getEigenVectors(DblMatrix const &cov_matrix, * @param radius Maximun size of event vectors that will be used * in calculating the standard deviation. */ -double Integrate3DEvents::stdDev(std::vector const &events, +double Integrate3DEvents::stdDev(std::vector > const &events, V3D const &direction, double radius) { double sum = 0; double sum_sq = 0; @@ -283,8 +290,9 @@ double Integrate3DEvents::stdDev(std::vector const &events, int count = 0; for (size_t i = 0; i < events.size(); i++) { - if (events[i].norm() <= radius) { - double dot_prod = events[i].scalar_prod(direction); + auto event = events[i].second; + if (event.norm() <= radius) { + double dot_prod = event.scalar_prod(direction); sum += dot_prod; sum_sq += dot_prod * dot_prod; count++; @@ -343,8 +351,8 @@ int64_t Integrate3DEvents::getHklKey(V3D const &q_vector) { * @param event_Q The Q-vector for the event that may be added to the * event_lists map, if it is close enough to some peak */ -void Integrate3DEvents::addEvent(V3D event_Q) { - int64_t hkl_key = getHklKey(event_Q); +void Integrate3DEvents::addEvent(std::pair event_Q) { + int64_t hkl_key = getHklKey(event_Q.second); if (hkl_key == 0) // don't keep events associated with 0,0,0 return; @@ -353,8 +361,8 @@ void Integrate3DEvents::addEvent(V3D event_Q) { if (peak_it != peak_qs.end()) { if (!peak_it->second.nullVector()) { - event_Q = event_Q - peak_it->second; - if (event_Q.norm() < radius) { + event_Q.second = event_Q.second - peak_it->second; + if (event_Q.second.norm() < radius) { event_lists[hkl_key].push_back(event_Q); } } @@ -397,7 +405,7 @@ void Integrate3DEvents::addEvent(V3D event_Q) { * */ PeakShapeEllipsoid_const_sptr Integrate3DEvents::ellipseIntegrateEvents( - std::vector const &ev_list, std::vector const &directions, + std::vector > const &ev_list, std::vector const &directions, std::vector const &sigmas, bool specify_size, double peak_radius, double back_inner_radius, double back_outer_radius, std::vector &axes_radii, double &inti, double &sigi) { diff --git a/Code/Mantid/Framework/MDEvents/src/IntegrateEllipsoids.cpp b/Code/Mantid/Framework/MDEvents/src/IntegrateEllipsoids.cpp index 211fc8d3e2f8..ddc312d7ea57 100644 --- a/Code/Mantid/Framework/MDEvents/src/IntegrateEllipsoids.cpp +++ b/Code/Mantid/Framework/MDEvents/src/IntegrateEllipsoids.cpp @@ -76,7 +76,7 @@ void qListFromEventWS(Integrate3DEvents &integrator, Progress &prog, double errorSq(1.); // ignorable garbage const std::vector &raw_events = events.getWeightedEventsNoTime(); - std::vector qList; + std::vector > qList; for (auto event = raw_events.begin(); event != raw_events.end(); ++event) { double val = unitConverter.convertUnits(event->tof()); qConverter->calcMatrixCoord(val, locCoord, signal, errorSq); @@ -84,7 +84,7 @@ void qListFromEventWS(Integrate3DEvents &integrator, Progress &prog, buffer[dim] = locCoord[dim]; } V3D qVec(buffer[0], buffer[1], buffer[2]); - qList.push_back(qVec); + qList.push_back(std::make_pair(event->m_weight, qVec)); } // end of loop over events in list integrator.addEvents(qList); @@ -125,7 +125,7 @@ void qListFromHistoWS(Integrate3DEvents &integrator, Progress &prog, double signal(1.); // ignorable garbage double errorSq(1.); // ignorable garbage - std::vector qList; + std::vector > qList; // TODO. we should be able to do this in an OMP loop. for (size_t j = 0; j < yVals.size(); ++j) { @@ -149,10 +149,8 @@ void qListFromHistoWS(Integrate3DEvents &integrator, Progress &prog, V3D qVec(buffer[0], buffer[1], buffer[2]); int yValCounts = int(yVal); // we deliberately truncate. // Account for counts in histograms by increasing the qList with the - // same q-poin - for (int k = 0; k < yValCounts; ++k) { - qList.push_back(qVec); // Not ideal to control the size dynamically? - } + // same q-point + qList.push_back(std::make_pair(yValCounts,qVec)); // Not ideal to control the size dynamically? } integrator.addEvents(qList); // We would have to put a lock around this. prog.report(); @@ -299,6 +297,7 @@ void IntegrateEllipsoids::exec() { size_t n_peaks = peak_ws->getNumberPeaks(); size_t indexed_count = 0; std::vector peak_q_list; + std::vector > qList; std::vector hkl_vectors; for (size_t i = 0; i < n_peaks; i++) // Note: we skip un-indexed peaks { @@ -307,6 +306,7 @@ void IntegrateEllipsoids::exec() { // just check for (0,0,0) { peak_q_list.push_back(V3D(peaks[i].getQLabFrame())); + qList.push_back(std::make_pair(1., V3D(peaks[i].getQLabFrame()))); V3D miller_ind((double)boost::math::iround(hkl[0]), (double)boost::math::iround(hkl[1]), (double)boost::math::iround(hkl[2])); @@ -345,7 +345,7 @@ void IntegrateEllipsoids::exec() { } // make the integrator - Integrate3DEvents integrator(peak_q_list, UBinv, radius); + Integrate3DEvents integrator(qList, UBinv, radius); // get the events and add // them to the inegrator @@ -379,10 +379,11 @@ void IntegrateEllipsoids::exec() { double inti; double sigi; std::vector principalaxis1,principalaxis2,principalaxis3; + V3D peak_q; for (size_t i = 0; i < n_peaks; i++) { V3D hkl(peaks[i].getH(), peaks[i].getK(), peaks[i].getL()); if (Geometry::IndexingUtils::ValidIndex(hkl, 1.0)) { - V3D peak_q(peaks[i].getQLabFrame()); + peak_q = peaks[i].getQLabFrame(); std::vector axes_radii; Mantid::Geometry::PeakShape_const_sptr shape = integrator.ellipseIntegrateEvents( @@ -451,10 +452,11 @@ void IntegrateEllipsoids::exec() { back_inner_radius = peak_radius; back_outer_radius = peak_radius * 1.25992105; // A factor of 2 ^ (1/3) will make the background // shell volume equal to the peak region volume. + V3D peak_q; for (size_t i = 0; i < n_peaks; i++) { V3D hkl(peaks[i].getH(), peaks[i].getK(), peaks[i].getL()); if (Geometry::IndexingUtils::ValidIndex(hkl, 1.0)) { - V3D peak_q(peaks[i].getQLabFrame()); + peak_q = peaks[i].getQLabFrame(); std::vector axes_radii; integrator.ellipseIntegrateEvents(peak_q, specify_size, peak_radius, back_inner_radius, back_outer_radius, diff --git a/Code/Mantid/Framework/MDEvents/test/Integrate3DEventsTest.h b/Code/Mantid/Framework/MDEvents/test/Integrate3DEventsTest.h index 56c313eced2b..687180da6092 100644 --- a/Code/Mantid/Framework/MDEvents/test/Integrate3DEventsTest.h +++ b/Code/Mantid/Framework/MDEvents/test/Integrate3DEventsTest.h @@ -35,14 +35,14 @@ class Integrate3DEventsTest : public CxxTest::TestSuite double sigi_some[] = { 27.4773, 26.533, 24.5561 }; // synthesize three peaks - std::vector peak_q_list; + std::vector > peak_q_list; V3D peak_1( 10, 0, 0 ); V3D peak_2( 0, 5, 0 ); V3D peak_3( 0, 0, 4 ); - peak_q_list.push_back( peak_1 ); - peak_q_list.push_back( peak_2 ); - peak_q_list.push_back( peak_3 ); + peak_q_list.push_back( std::make_pair( 1., peak_1 ) ); + peak_q_list.push_back( std::make_pair( 1., peak_2 ) ); + peak_q_list.push_back( std::make_pair( 1., peak_3 ) ); // synthesize a UB-inverse to map DblMatrix UBinv(3,3,false); // Q to h,k,l UBinv.setRow( 0, V3D( .1, 0, 0 ) ); @@ -55,31 +55,31 @@ class Integrate3DEventsTest : public CxxTest::TestSuite // around peak 1, 704 events around // peak 2, and 603 events around // peak 3. - std::vector event_Qs; + std::vector > event_Qs; for ( int i = -100; i <= 100; i++ ) { - event_Qs.push_back( V3D( peak_1 + V3D( (double)i/100.0, 0, 0 ) ) ); - event_Qs.push_back( V3D( peak_2 + V3D( (double)i/100.0, 0, 0 ) ) ); - event_Qs.push_back( V3D( peak_3 + V3D( (double)i/100.0, 0, 0 ) ) ); + event_Qs.push_back ( std::make_pair( 1., V3D( peak_1 + V3D( (double)i/100.0, 0, 0 ) ) ) ); + event_Qs.push_back ( std::make_pair( 1., V3D( peak_2 + V3D( (double)i/100.0, 0, 0 ) ) ) ); + event_Qs.push_back ( std::make_pair( 1., V3D( peak_3 + V3D( (double)i/100.0, 0, 0 ) ) ) ); - event_Qs.push_back( V3D( peak_1 + V3D( 0, (double)i/200.0, 0 ) ) ); - event_Qs.push_back( V3D( peak_2 + V3D( 0, (double)i/200.0, 0 ) ) ); - event_Qs.push_back( V3D( peak_3 + V3D( 0, (double)i/200.0, 0 ) ) ); + event_Qs.push_back ( std::make_pair( 1., V3D( peak_1 + V3D( 0, (double)i/200.0, 0 ) ) ) ); + event_Qs.push_back ( std::make_pair( 1., V3D( peak_2 + V3D( 0, (double)i/200.0, 0 ) ) ) ); + event_Qs.push_back ( std::make_pair( 1., V3D( peak_3 + V3D( 0, (double)i/200.0, 0 ) ) ) ); - event_Qs.push_back( V3D( peak_1 + V3D( 0, 0, (double)i/300.0 ) ) ); - event_Qs.push_back( V3D( peak_2 + V3D( 0, 0, (double)i/300.0 ) ) ); - event_Qs.push_back( V3D( peak_3 + V3D( 0, 0, (double)i/300.0 ) ) ); + event_Qs.push_back ( std::make_pair( 1., V3D( peak_1 + V3D( 0, 0, (double)i/300.0 ) ) ) ); + event_Qs.push_back ( std::make_pair( 1., V3D( peak_2 + V3D( 0, 0, (double)i/300.0 ) ) ) ); + event_Qs.push_back ( std::make_pair( 1., V3D( peak_3 + V3D( 0, 0, (double)i/300.0 ) ) ) ); } for ( int i = -50; i <= 50; i++ ) { - event_Qs.push_back( V3D( peak_1 + V3D( 0, (double)i/147.0, 0 ) ) ); - event_Qs.push_back( V3D( peak_2 + V3D( 0, (double)i/147.0, 0 ) ) ); + event_Qs.push_back ( std::make_pair( 1., V3D( peak_1 + V3D( 0, (double)i/147.0, 0 ) ) ) ); + event_Qs.push_back ( std::make_pair( 1., V3D( peak_2 + V3D( 0, (double)i/147.0, 0 ) ) ) ); } for ( int i = -25; i <= 25; i++ ) { - event_Qs.push_back( V3D( peak_1 + V3D( 0, 0, (double)i/61.0 ) ) ); + event_Qs.push_back ( std::make_pair(1., V3D( peak_1 + V3D( 0, 0, (double)i/61.0 ) ) ) ); } double radius = 1.3; @@ -98,7 +98,7 @@ class Integrate3DEventsTest : public CxxTest::TestSuite double sigi; for ( size_t i = 0; i < peak_q_list.size(); i++ ) { - auto shape = integrator.ellipseIntegrateEvents( peak_q_list[i], specify_size, + auto shape = integrator.ellipseIntegrateEvents( 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); @@ -114,7 +114,7 @@ class Integrate3DEventsTest : public CxxTest::TestSuite specify_size = false; for ( size_t i = 0; i < peak_q_list.size(); i++ ) { - integrator.ellipseIntegrateEvents( peak_q_list[i], specify_size, + integrator.ellipseIntegrateEvents( 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);