Skip to content

Commit

Permalink
Refs #11239 Use weights for events
Browse files Browse the repository at this point in the history
  • Loading branch information
@VickieLynch
VickieLynch committed Mar 11, 2015
1 parent 40e5686 commit 13b0032
Show file tree
Hide file tree
Showing 4 changed files with 71 additions and 59 deletions.
18 changes: 9 additions & 9 deletions Code/Mantid/Framework/MDEvents/inc/MantidMDEvents/Integrate3DEvents.h
View file
Expand Up @@ -55,43 +55,43 @@ namespace MDEvents {
*/

typedef Mantid::Kernel::Matrix<double> DblMatrix;
typedef boost::unordered_map<int64_t, std::vector<Mantid::Kernel::V3D>> EventListMap;
typedef boost::unordered_map<int64_t, std::vector<std::pair<double, Mantid::Kernel::V3D> > > EventListMap;
typedef boost::unordered_map<int64_t, Mantid::Kernel::V3D> PeakQMap;

class DLLExport Integrate3DEvents {
public:
/// Construct object to store events around peaks and integrate peaks
Integrate3DEvents(std::vector<Mantid::Kernel::V3D> const &peak_q_list, DblMatrix const &UBinv,
Integrate3DEvents(std::vector<std::pair<double, Mantid::Kernel::V3D> > const &peak_q_list, DblMatrix const &UBinv,
double radius);

~Integrate3DEvents();

/// Add event Q's to lists of events near peaks
void addEvents(std::vector<Mantid::Kernel::V3D> const &event_qs);
void addEvents(std::vector<std::pair<double, Mantid::Kernel::V3D> > const &event_qs);

/// 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::pair<double, 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,
double &sigi);

private:
/// Calculate the number of events in an ellipsoid centered at 0,0,0
static int numInEllipsoid(std::vector<Mantid::Kernel::V3D> const &events,
static double numInEllipsoid(std::vector<std::pair<double, Mantid::Kernel::V3D> > const &events,
std::vector<Mantid::Kernel::V3D> const &directions,
std::vector<double> const &sizes);

/// Calculate the 3x3 covariance matrix of a list of Q-vectors at 0,0,0
static void makeCovarianceMatrix(std::vector<Mantid::Kernel::V3D> const &events,
static void makeCovarianceMatrix(std::vector<std::pair<double, Mantid::Kernel::V3D> > const &events,
DblMatrix &matrix, double radius);

/// Calculate the eigen vectors of a 3x3 real symmetric matrix
static void getEigenVectors(DblMatrix const &cov_matrix,
std::vector<Mantid::Kernel::V3D> &eigen_vectors);

/// Calculate the standard deviation of 3D events in a specified direction
static double stdDev(std::vector<Mantid::Kernel::V3D> const &events, Mantid::Kernel::V3D const &direction,
static double stdDev(std::vector<std::pair<double, Mantid::Kernel::V3D> > 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
Expand All @@ -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<double, Mantid::Kernel::V3D> event_Q);

/// Find the net integrated intensity of a list of Q's using ellipsoids
boost::shared_ptr<const Mantid::DataObjects::PeakShapeEllipsoid> ellipseIntegrateEvents(
std::vector<Mantid::Kernel::V3D> const &ev_list, std::vector<Mantid::Kernel::V3D> const &directions,
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);
Expand Down
54 changes: 31 additions & 23 deletions Code/Mantid/Framework/MDEvents/src/Integrate3DEvents.cpp
View file
Expand Up @@ -33,16 +33,16 @@ using Mantid::Kernel::V3D;
* an event to be stored in the list associated with
* that peak.
*/
Integrate3DEvents::Integrate3DEvents(std::vector<V3D> const &peak_q_list,
Integrate3DEvents::Integrate3DEvents(std::vector<std::pair<double, V3D>> 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;
}
}

Expand All @@ -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<V3D> const &event_qs) {
void Integrate3DEvents::addEvents(std::vector<std::pair<double, V3D> > const &event_qs) {
for (size_t i = 0; i < event_qs.size(); i++) {
addEvent(event_qs[i]);
}
Expand Down Expand Up @@ -109,18 +109,23 @@ void Integrate3DEvents::addEvents(std::vector<V3D> const &event_qs) {
*
*/
Mantid::Geometry::PeakShape_const_sptr Integrate3DEvents::ellipseIntegrateEvents(
V3D const &peak_q, bool specify_size, double peak_radius,
std::pair<double, 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
sigi = 0.0; // is wrong with the peak.

int64_t hkl_key = getHklKey(peak_q);
int64_t hkl_key = getHklKey(peak_q.second);
if (hkl_key == 0) {
return boost::make_shared<NoShape>();
}

std::vector<V3D> &some_events = event_lists[hkl_key];
std::vector<std::pair<double, V3D> > &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
Expand Down Expand Up @@ -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<V3D> const &events,
double Integrate3DEvents::numInEllipsoid(std::vector<std::pair<double, V3D>> const &events,
std::vector<V3D> const &directions,
std::vector<double> 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;
Expand Down Expand Up @@ -209,14 +214,16 @@ int Integrate3DEvents::numInEllipsoid(std::vector<V3D> const &events,
* peak center (0,0,0) will be used for
* calculating the covariance matrix.
*/
void Integrate3DEvents::makeCovarianceMatrix(std::vector<V3D> const &events,

void Integrate3DEvents::makeCovarianceMatrix(std::vector<std::pair<double, V3D> > 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)
Expand Down Expand Up @@ -275,16 +282,17 @@ 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<V3D> const &events,
double Integrate3DEvents::stdDev(std::vector<std::pair<double, V3D> > const &events,
V3D const &direction, double radius) {
double sum = 0;
double sum_sq = 0;
double stdev = 0;
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++;
Expand Down Expand Up @@ -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<double, V3D> event_Q) {
int64_t hkl_key = getHklKey(event_Q.second);

if (hkl_key == 0) // don't keep events associated with 0,0,0
return;
Expand All @@ -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);
}
}
Expand Down Expand Up @@ -397,7 +405,7 @@ void Integrate3DEvents::addEvent(V3D event_Q) {
*
*/
PeakShapeEllipsoid_const_sptr Integrate3DEvents::ellipseIntegrateEvents(
std::vector<V3D> const &ev_list, std::vector<V3D> const &directions,
std::vector<std::pair<double, V3D> > const &ev_list, std::vector<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) {
Expand Down
24 changes: 14 additions & 10 deletions Code/Mantid/Framework/MDEvents/src/IntegrateEllipsoids.cpp
View file
Expand Up @@ -76,15 +76,15 @@ void qListFromEventWS(Integrate3DEvents &integrator, Progress &prog,
double errorSq(1.); // ignorable garbage
const std::vector<WeightedEventNoTime> &raw_events =
events.getWeightedEventsNoTime();
std::vector<Mantid::Kernel::V3D> qList;
std::vector<std::pair<double, V3D> > qList;
for (auto event = raw_events.begin(); event != raw_events.end(); ++event) {
double val = unitConverter.convertUnits(event->tof());
qConverter->calcMatrixCoord(val, locCoord, signal, errorSq);
for (size_t dim = 0; dim < DIMS; ++dim) {
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);
Expand Down Expand Up @@ -125,7 +125,7 @@ void qListFromHistoWS(Integrate3DEvents &integrator, Progress &prog,
double signal(1.); // ignorable garbage
double errorSq(1.); // ignorable garbage

std::vector<V3D> qList;
std::vector<std::pair<double, V3D> > qList;

// TODO. we should be able to do this in an OMP loop.
for (size_t j = 0; j < yVals.size(); ++j) {
Expand All @@ -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();
Expand Down Expand Up @@ -299,6 +297,7 @@ void IntegrateEllipsoids::exec() {
size_t n_peaks = peak_ws->getNumberPeaks();
size_t indexed_count = 0;
std::vector<V3D> peak_q_list;
std::vector<std::pair<double, V3D> > qList;
std::vector<V3D> hkl_vectors;
for (size_t i = 0; i < n_peaks; i++) // Note: we skip un-indexed peaks
{
Expand All @@ -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<double>(hkl[0]),
(double)boost::math::iround<double>(hkl[1]),
(double)boost::math::iround<double>(hkl[2]));
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -379,10 +379,12 @@ void IntegrateEllipsoids::exec() {
double inti;
double sigi;
std::vector<double> principalaxis1,principalaxis2,principalaxis3;
std::pair <double, 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.first = 1.;
peak_q.second = peaks[i].getQLabFrame();
std::vector<double> axes_radii;
Mantid::Geometry::PeakShape_const_sptr shape =
integrator.ellipseIntegrateEvents(
Expand Down Expand Up @@ -451,10 +453,12 @@ 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.
std::pair <double, 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.first = 1.;
peak_q.second = peaks[i].getQLabFrame();
std::vector<double> axes_radii;
integrator.ellipseIntegrateEvents(peak_q, specify_size, peak_radius,
back_inner_radius, back_outer_radius,
Expand Down
34 changes: 17 additions & 17 deletions Code/Mantid/Framework/MDEvents/test/Integrate3DEventsTest.h
View file
Expand Up @@ -35,14 +35,14 @@ class Integrate3DEventsTest : public CxxTest::TestSuite
double sigi_some[] = { 27.4773, 26.533, 24.5561 };

// synthesize three peaks
std::vector<V3D> peak_q_list;
std::vector<std::pair<double, V3D> > 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 ) );
Expand All @@ -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<V3D> event_Qs;
std::vector<std::pair<double, V3D> > 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;
Expand Down

0 comments on commit 13b0032

Please sign in to comment.