Skip to content

Commit

Permalink
Refs #4202. Fix normalisation of SofQW2
Browse files Browse the repository at this point in the history 
Also fix SofQW for new BASIS definition where some detectors have
no efixed.
  • Loading branch information
@martyngigg
martyngigg committed Jan 27, 2012
1 parent 5e3bf09 commit 2b1a6cb
Show file tree
Hide file tree
Showing 3 changed files with 76 additions and 24 deletions.
2 changes: 2 additions & 0 deletions Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/SofQW2.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -102,6 +102,8 @@ namespace Mantid
std::vector<double> m_thetaPts;
/// Theta width
double m_thetaWidth;
///
API::MatrixWorkspace_sptr m_numIntersectionsWS;
};


Expand Down
12 changes: 10 additions & 2 deletions Code/Mantid/Framework/Algorithms/src/SofQW.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -129,16 +129,20 @@ void SofQW::exec()
if (emode==2)
{
try {
Parameter_sptr par = pmap.get(spectrumDet.get(),"Efixed");
Parameter_sptr par = pmap.get(spectrumDet.get(),"EFixed");
if (par)
{
efixed = par->value<double>();
}
else
{
continue;
}
} catch (std::runtime_error&) { /* Throws if a DetectorGroup, use single provided value */ }
}

// For inelastic scattering the simple relationship q=4*pi*sinTheta/lambda does not hold. In order to
// be completely general wemust calculate the momentum transfer by calculating the incident and final
// be completely general we must calculate the momentum transfer by calculating the incident and final
// wave vectors and then use |q| = sqrt[(ki - kf)*(ki - kf)]
DetectorGroup_const_sptr detGroup = boost::dynamic_pointer_cast<const DetectorGroup>(spectrumDet);
std::vector<IDetector_const_sptr> detectors;
Expand Down Expand Up @@ -238,7 +242,11 @@ API::MatrixWorkspace_sptr SofQW::setUpOutputWorkspace(API::MatrixWorkspace_const

// Set the axis units
verticalAxis->unit() = UnitFactory::Instance().create("MomentumTransfer");
verticalAxis->title() = "|Q|";

// Set the X axis title (for conversion to MD)
outputWorkspace->getAxis(0)->title() = "Energy transfer";

setProperty("OutputWorkspace",outputWorkspace);
return outputWorkspace;
}
Expand Down
86 changes: 64 additions & 22 deletions Code/Mantid/Framework/Algorithms/src/SofQW2.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,8 +17,6 @@ Converts a 2D workspace from units of spectrum number/energy transfer to the in
#include "MantidGeometry/Math/Vertex2D.h"
#include "MantidGeometry/Instrument/DetectorGroup.h"

#include "MantidKernel/Timer.h"

namespace Mantid
{
namespace Algorithms
Expand All @@ -27,7 +25,7 @@ namespace Mantid
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(SofQW2)

using namespace Mantid::Kernel;
using namespace Mantid::Kernel;
using namespace Mantid::API;
using Geometry::IDetector_const_sptr;
using Geometry::DetectorGroup;
Expand Down Expand Up @@ -64,14 +62,20 @@ namespace Mantid
}

MatrixWorkspace_sptr outputWS = this->setUpOutputWorkspace(inputWS, m_Qout);
// We also need to keep track of how many intersections went into a pixel so that we can
// normalize correctly
m_numIntersectionsWS = WorkspaceFactory::Instance().create(outputWS);

const size_t nOutputHist(outputWS->getNumberHistograms());
const size_t nenergyBins = inputWS->blocksize();

//Progress reports & cancellation
const size_t nreports(static_cast<size_t>(inputWS->getNumberHistograms()*inputWS->blocksize()));
m_progress = boost::shared_ptr<API::Progress>(new API::Progress(this, 0.0, 1.0, nreports));

// Compute input caches
initCachedValues(inputWS);
const size_t nTheta = m_thetaPts.size();
const size_t nenergyBins = inputWS->blocksize();
const MantidVec & X = inputWS->readX(0);

// Select the calculate Q method based on the mode
Expand All @@ -93,6 +97,11 @@ namespace Mantid
PARALLEL_START_INTERUPT_REGION

const double theta = m_thetaPts[i];
if( theta < 0.0 ) // One to skip
{
continue;
}

const double efixed = getEFixed(inputWS->getDetector(i));
const double thetaLower = theta - 0.5*m_thetaWidth;
const double thetaUpper = theta + 0.5*m_thetaWidth;
Expand All @@ -118,25 +127,31 @@ namespace Mantid
}
PARALLEL_CHECK_INTERUPT_REGION

const size_t nOutputHist(outputWS->getNumberHistograms());
// The errors need square-rooting
/**
* Normalise the output by the total accumulated fraction
* and square root the errors
*/

PARALLEL_FOR1(outputWS)
for(int64_t i = 0; i < static_cast<int64_t>(nOutputHist); ++i)
{
PARALLEL_START_INTERUPT_REGION

MantidVec& errors = outputWS->dataE(i);
std::transform(errors.begin(), errors.end(), errors.begin(), (double (*)(double))std::sqrt);
MantidVec & outputY = outputWS->dataY(i);
MantidVec & outputE = outputWS->dataE(i);
const MantidVec & noIntersects = m_numIntersectionsWS->readY(i);
for(size_t j = 0; j < nenergyBins; ++j)
{
m_progress->report("Calculating errors and normalising");
const double numIntersects = noIntersects[j];
if(numIntersects > 0.0) outputY[j] /= numIntersects;
outputE[j] = std::sqrt(outputE[j]);
}

PARALLEL_END_INTERUPT_REGION
}
PARALLEL_CHECK_INTERUPT_REGION

// Divide by the bin width if the input is a distribution
if( inputWS->isDistribution() )
{
this->makeDistribution(outputWS, m_Qout);
}
}

/**
Expand Down Expand Up @@ -167,11 +182,11 @@ namespace Mantid
ConvexPolygon overlap = intersectionByLaszlo(outputQ, inputQ);
const double weight = overlap.area()/inputQ.area();
outputWS->dataY(qi)[ei] += inputWS->readY(i)[j] * weight;
m_numIntersectionsWS->dataY(qi)[ei] += weight;
outputWS->dataE(qi)[ei] += std::pow(inputWS->readE(i)[j] * weight, 2);
}
catch(Geometry::NoIntersectionException &)
{}

}
}
}
Expand Down Expand Up @@ -295,8 +310,6 @@ namespace Mantid
*/
void SofQW2::initCachedValues(API::MatrixWorkspace_const_sptr workspace)
{
m_progress->report("Initializing caches");

// Retrieve the emode & efixed properties
const std::string emode = getProperty("EMode");
// Convert back to an integer representation
Expand Down Expand Up @@ -355,19 +368,30 @@ namespace Mantid
*/
void SofQW2::initThetaCache(API::MatrixWorkspace_const_sptr workspace)
{
Mantid::Kernel::Timer clock;
const size_t nhist = workspace->getNumberHistograms();
m_thetaPts = std::vector<double>(nhist);
size_t ndets(0);
double minTheta(DBL_MAX), maxTheta(-DBL_MAX);

PARALLEL_FOR1(workspace)
for(int64_t i = 0 ; i < (int64_t)nhist; ++i) //signed for OpenMP
{
PARALLEL_START_INTERUPT_REGION

m_progress->report("Calculating detector angles");
IDetector_const_sptr det;
try
{
det = workspace->getDetector(i);
// Check to see if there is an EFixed, if not skip it
try
{
getEFixed(det);
}
catch(std::runtime_error&)
{
det.reset();
}
}
catch(Kernel::Exception::NotFoundError&)
{
Expand All @@ -376,16 +400,34 @@ namespace Mantid
// in an openmp block.
}
// If no detector found, skip onto the next spectrum
if( !det || det->isMonitor() ) continue;
m_thetaPts[i] = workspace->detectorTwoTheta(det);
if( !det || det->isMonitor() )
{
m_thetaPts[i] = -1.0; // Indicates a detector to skip
}
else
{
PARALLEL_ATOMIC
++ndets;
const double theta = workspace->detectorTwoTheta(det);
m_thetaPts[i] = theta;
if( theta < minTheta )
{
PARALLEL_CRITICAL(minTheta)
minTheta = theta;
}
else if( theta > maxTheta )
{
PARALLEL_CRITICAL(maxTheta)
maxTheta = theta;
}
}

PARALLEL_END_INTERUPT_REGION
}
PARALLEL_CHECK_INTERUPT_REGION

m_thetaWidth = std::abs(m_thetaPts.back() - m_thetaPts.front())/static_cast<double>(m_thetaPts.size());

g_log.debug() << "Initializing theta Cache took " << clock.elapsed() << " seconds\n";
m_thetaWidth = (maxTheta - minTheta)/static_cast<double>(ndets);
g_log.information() << "Calculated detector width in theta=" << (m_thetaWidth*180.0/M_PI) << " degrees.\n";
}


Expand Down

0 comments on commit 2b1a6cb

Please sign in to comment.