Skip to content

Commit

Permalink
Enabled to skip spectra without detectors. Refs #9976.
Browse files Browse the repository at this point in the history
  • Loading branch information
@wdzhou
wdzhou committed Jul 22, 2014
1 parent d2f2a67 commit 393d940
Showing 1 changed file with 131 additions and 127 deletions.
258 changes: 131 additions & 127 deletions Code/Mantid/Framework/Algorithms/src/FilterEvents.cpp
View file
Expand Up @@ -315,20 +315,37 @@ namespace Algorithms

for (size_t i = 0; i < numhist; ++i)
{
const EventList& elist = m_eventWS->getEventList(i);
const std::set<detid_t>& detids = elist.getDetectorIDs();
if (detids.size() == 0)
bool skip = false;

// Access detector of the spectrum
try
{
IDetector_const_sptr tempdet = m_eventWS->getDetector(i);
if (!tempdet)
skip = true;
}
catch (const Kernel::Exception::NotFoundError& err)
{
// No detector found
skip = true;
}

// Output
if (skip)
{
++ numskipspec;
m_vecSkip[i] = true;

++ numskipspec;
const EventList& elist = m_eventWS->getEventList(i);
numeventsskip += elist.getNumberEvents();
msgss << i;
if (numhist % 10 == 0)
if (numskipspec % 10 == 0)
msgss << "\n";
else
msgss << ",";
}
}

} // ENDFOR

if (numskipspec > 0)
{
Expand All @@ -341,7 +358,8 @@ namespace Algorithms
{
g_log.notice("There is no spectrum that does not have detectors.");
}
}

} // END-IF-ELSE

return;
}
Expand Down Expand Up @@ -609,10 +627,9 @@ namespace Algorithms

// Get
size_t numhist = m_eventWS->getNumberHistograms();
size_t numhistnodet = 0;
for (size_t i = 0; i < numhist; ++i)
{
try
if (!m_vecSkip[i])
{
IComponent_const_sptr tmpdet = boost::dynamic_pointer_cast<const IComponent>(m_eventWS->getDetector(i));
double l2 = instrument->getDistance(*tmpdet);
Expand All @@ -622,37 +639,13 @@ namespace Algorithms
m_detTofOffsets[i] = corrfactor;
corrws->dataY(i)[0] = corrfactor;
}
catch (const Kernel::Exception::NotFoundError& err)
{
const EventList& evlist = m_eventWS->getEventList(i);
const set<detid_t> detset = evlist.getDetectorIDs();

stringstream msgss;
for (set<detid_t>::iterator it = detset.begin(); it != detset.end(); ++it)
{
msgss << "DetID = " << *it << ", ";
IDetector_const_sptr det = m_eventWS->getDetectorByID(*it);
// if (!det)
// g_log.error() << "There is no detector defined for detector ID = " << *it << "\n";
}
g_log.error() << "Unable to find detector of spectrum " << i << "; EventList shows " << detset.size()
<< " detectors. They are " << msgss.str() << "\n";
++ numhistnodet;
}
}

if (numhistnodet > 0)
{
stringstream err;
err << "There are " << numhistnodet << " that do not have detectors. Unable to create to create correction.";
throw runtime_error(err.str());
}

return;
}

//----------------------------------------------------------------------------------------------
/**
/** Calculate TOF correction for direct geometry inelastic instrument
* Time = T_pulse + TOF*0 + L1/sqrt(E*2/m)
*/
void FilterEvents::setupDirectTOFCorrection(API::MatrixWorkspace_sptr corrws)
Expand Down Expand Up @@ -700,7 +693,7 @@ namespace Algorithms
}

//----------------------------------------------------------------------------------------------
/**
/** Calculate TOF correction for indirect geometry inelastic instrument
* Time = T_pulse + TOF - L2/sqrt(E_fix * 2 * meV / mass)
*/
void FilterEvents::setupIndirectTOFCorrection(API::MatrixWorkspace_sptr corrws)
Expand All @@ -721,58 +714,61 @@ namespace Algorithms

for (size_t i = 0; i < numhist; ++i)
{
double shift;

IDetector_const_sptr det = m_eventWS->getDetector(i);
if (!det->isMonitor())
if (!m_vecSkip[i])
{
// Get E_fix
double efix = 0.;
try
double shift;
IDetector_const_sptr det = m_eventWS->getDetector(i);
if (!det->isMonitor())
{
Parameter_sptr par = pmap.getRecursive(det.get(),"Efixed");
if (par)
// Get E_fix
double efix = 0.;
try
{
efix = par->value<double>();
g_log.debug() << "Detector: " << det->getID() << " of spectrum " << i << " EFixed: " << efix << "\n";
Parameter_sptr par = pmap.getRecursive(det.get(),"Efixed");
if (par)
{
efix = par->value<double>();
g_log.debug() << "Detector: " << det->getID() << " of spectrum " << i << " EFixed: " << efix << "\n";
}
else
{
g_log.warning() << "Detector: " << det->getID() << " of spectrum " << i << " does not have EFixed set up."
<< "\n";
}
}
else
catch (std::runtime_error&)
{
g_log.warning() << "Detector: " << det->getID() << " of spectrum " << i << " does not have EFixed set up."
<< "\n";
// Throws if a DetectorGroup, use single provided value
stringstream errmsg;
errmsg << "Inelastic instrument detector " << det->getID() << " of spectrum " << i
<< " does not have EFixed ";
throw runtime_error(errmsg.str());
}

// Get L2
double l2 = det->getPos().distance(samplepos);

// Calculate shift
shift = -1. * l2 / sqrt(efix * twomev_d_mass);

g_log.notice() << "Detector " << i << ": " << "L2 = " << l2 << ", EFix = " << efix << ", Shift = " << shift << "\n";
}
catch (std::runtime_error&)
else
{
// Throws if a DetectorGroup, use single provided value
stringstream errmsg;
errmsg << "Inelastic instrument detector " << det->getID() << " of spectrum " << i
<< " does not have EFixed ";
throw runtime_error(errmsg.str());
}
// Monitor:
g_log.warning() << "Spectrum " << i << " contains detector " << det->getID() << " is a monitor. " << "\n";

// Get L2
double l2 = det->getPos().distance(samplepos);
shift = 0.;
}

// Calculate shift
shift = -1. * l2 / sqrt(efix * twomev_d_mass);
// Set up the shifts
m_detTofOffsets[i] = 1.0;
m_detTofShifts[i] = shift;

g_log.notice() << "Detector " << i << ": " << "L2 = " << l2 << ", EFix = " << efix << ", Shift = " << shift << "\n";
}
else
{
// Monitor:
g_log.warning() << "Spectrum " << i << " contains detector " << det->getID() << " is a monitor. " << "\n";
corrws->dataY(i)[0] = 1.0;
corrws->dataY(i)[1] = shift;

shift = 0.;
}

// Set up the shifts
m_detTofOffsets[i] = 1.0;
m_detTofShifts[i] = shift;

corrws->dataY(i)[0] = 1.0;
corrws->dataY(i)[1] = shift;
} // ENDOF (all spectra)

return;
Expand Down Expand Up @@ -958,33 +954,37 @@ namespace Algorithms
{
PARALLEL_START_INTERUPT_REGION

// Get the output event lists (should be empty) to be a map
std::map<int, DataObjects::EventList* > outputs;
PARALLEL_CRITICAL(build_elist)
// Filter the non-skipped
if (!m_vecSkip[iws])
{
for (wsiter = m_outputWS.begin(); wsiter != m_outputWS.end(); ++ wsiter)
// Get the output event lists (should be empty) to be a map
std::map<int, DataObjects::EventList* > outputs;
PARALLEL_CRITICAL(build_elist)
{
int index = wsiter->first;
DataObjects::EventList* output_el = wsiter->second->getEventListPtr(iws);
outputs.insert(std::make_pair(index, output_el));
for (wsiter = m_outputWS.begin(); wsiter != m_outputWS.end(); ++ wsiter)
{
int index = wsiter->first;
DataObjects::EventList* output_el = wsiter->second->getEventListPtr(iws);
outputs.insert(std::make_pair(index, output_el));
}
}
}

// Get a holder on input workspace's event list of this spectrum
const DataObjects::EventList& input_el = m_eventWS->getEventList(iws);
// Get a holder on input workspace's event list of this spectrum
const DataObjects::EventList& input_el = m_eventWS->getEventList(iws);

// Perform the filtering (using the splitting function and just one output)
if (mFilterByPulseTime)
{
input_el.splitByPulseTime(m_splitters, outputs);
}
else if (m_tofCorrType != NoneCorrect)
{
input_el.splitByFullTime(m_splitters, outputs, true, m_detTofOffsets[iws], m_detTofShifts[iws]);
}
else
{
input_el.splitByFullTime(m_splitters, outputs, false, 1.0, 0.0);
// Perform the filtering (using the splitting function and just one output)
if (mFilterByPulseTime)
{
input_el.splitByPulseTime(m_splitters, outputs);
}
else if (m_tofCorrType != NoneCorrect)
{
input_el.splitByFullTime(m_splitters, outputs, true, m_detTofOffsets[iws], m_detTofShifts[iws]);
}
else
{
input_el.splitByFullTime(m_splitters, outputs, false, 1.0, 0.0);
}
}

PARALLEL_END_INTERUPT_REGION
Expand Down Expand Up @@ -1062,44 +1062,48 @@ namespace Algorithms
{
PARALLEL_START_INTERUPT_REGION

// Get the output event lists (should be empty) to be a map
map<int, DataObjects::EventList* > outputs;
PARALLEL_CRITICAL(build_elist)
// Filter the non-skipped spectrum
if (!m_vecSkip[iws])
{
for (wsiter = m_outputWS.begin(); wsiter != m_outputWS.end(); ++ wsiter)
// Get the output event lists (should be empty) to be a map
map<int, DataObjects::EventList* > outputs;
PARALLEL_CRITICAL(build_elist)
{
int index = wsiter->first;
DataObjects::EventList* output_el = wsiter->second->getEventListPtr(iws);
outputs.insert(std::make_pair(index, output_el));
for (wsiter = m_outputWS.begin(); wsiter != m_outputWS.end(); ++ wsiter)
{
int index = wsiter->first;
DataObjects::EventList* output_el = wsiter->second->getEventListPtr(iws);
outputs.insert(std::make_pair(index, output_el));
}
}
}

// Get a holder on input workspace's event list of this spectrum
const DataObjects::EventList& input_el = m_eventWS->getEventList(iws);
// Get a holder on input workspace's event list of this spectrum
const DataObjects::EventList& input_el = m_eventWS->getEventList(iws);

bool printdetail = false;
if (m_useDBSpectrum)
printdetail = (iws == static_cast<int64_t>(m_dbWSIndex));
bool printdetail = false;
if (m_useDBSpectrum)
printdetail = (iws == static_cast<int64_t>(m_dbWSIndex));

// Perform the filtering (using the splitting function and just one output)
std::string logmessage("");
if (mFilterByPulseTime)
{
throw runtime_error("It is not a good practice to split fast event by pulse time. ");
}
else if (m_tofCorrType != NoneCorrect)
{
logmessage = input_el.splitByFullTimeMatrixSplitter(m_vecSplitterTime, m_vecSplitterGroup, outputs,
true, m_detTofOffsets[iws], m_detTofShifts[iws]);
}
else
{
logmessage = input_el.splitByFullTimeMatrixSplitter(m_vecSplitterTime, m_vecSplitterGroup, outputs,
false, 1.0, 0.0);
}
// Perform the filtering (using the splitting function and just one output)
std::string logmessage("");
if (mFilterByPulseTime)
{
throw runtime_error("It is not a good practice to split fast event by pulse time. ");
}
else if (m_tofCorrType != NoneCorrect)
{
logmessage = input_el.splitByFullTimeMatrixSplitter(m_vecSplitterTime, m_vecSplitterGroup, outputs,
true, m_detTofOffsets[iws], m_detTofShifts[iws]);
}
else
{
logmessage = input_el.splitByFullTimeMatrixSplitter(m_vecSplitterTime, m_vecSplitterGroup, outputs,
false, 1.0, 0.0);
}

if (printdetail)
g_log.notice(logmessage);
if (printdetail)
g_log.notice(logmessage);
}

PARALLEL_END_INTERUPT_REGION
} // END FOR i = 0
Expand Down

0 comments on commit 393d940

Please sign in to comment.