Skip to content

Commit

Permalink
Refs #5021: method in EventList to create it from a histogram
Browse files Browse the repository at this point in the history 
for code reuse from ConvertToEventWorkspace
  • Loading branch information
Janik Zikovsky committed Mar 23, 2012
1 parent d456ed7 commit e3edcdf
Show file tree
Hide file tree
Showing 6 changed files with 155 additions and 69 deletions.
76 changes: 8 additions & 68 deletions Code/Mantid/Framework/Algorithms/src/ConvertToEventWorkspace.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,6 +64,9 @@ namespace Algorithms
declareProperty("GenerateMultipleEvents", false,
"Generate a number of evenly spread events in each bin. See the help for details.\n"
"Warning! This may use significantly more memory.");
declareProperty("GenerateZeros", false,
"Generate an event even for empty bins\n"
"Warning! This may use significantly more memory.");
declareProperty("MaxEventsPerBin", 10,
"If GenerateMultipleEvents is true, specifies a maximum number of events to generate in a single bin.\n"
"Use a value that matches your instrument's TOF resolution. Default 10.");
Expand All @@ -79,6 +82,7 @@ namespace Algorithms
Workspace2D_const_sptr inWS = getProperty("InputWorkspace");

bool GenerateMultipleEvents = getProperty("GenerateMultipleEvents");
bool GenerateZeros = getProperty("GenerateZeros");
int MaxEventsPerBin = getProperty("MaxEventsPerBin");

//Create the event workspace
Expand All @@ -88,85 +92,21 @@ namespace Algorithms
//Copy geometry, etc. over.
API::WorkspaceFactory::Instance().initializeFromParent(inWS, outWS, false);

// Cached values for later checks
double inf = std::numeric_limits<double>::infinity();
double ninf = -inf;

Progress prog(this, 0.0, 1.0, inWS->getNumberHistograms());
PARALLEL_FOR1(inWS)
for (int iwi=0; iwi<int(inWS->getNumberHistograms()); iwi++)
{
PARALLEL_START_INTERUPT_REGION
size_t wi = size_t(iwi);

// The input spectrum (a histogram)
const ISpectrum * inSpec = inWS->getSpectrum(wi);
const MantidVec & X = inSpec->readX();
const MantidVec & Y = inSpec->readY();
const MantidVec & E = inSpec->readE();
if (Y.size()+1 != X.size())
throw std::runtime_error("Expected a histogram (X vector should be 1 longer than the Y vector)");

// The output event list
EventList & el = outWS->getEventList(wi);
// Copy detector IDs and spectra
el.copyInfoFrom( *inSpec );
// We need weights but have no way to set the time. So use weighted, no time
el.switchTo(WEIGHTED_NOTIME);

for (size_t i=0; i<X.size()-1; i++)
{
double weight = Y[i];
if ((weight != 0.0) && (weight == weight) /*NAN check*/
&& (weight != inf) && (weight != ninf))
{
double error = E[i];
// Also check that the error is not a bad number
if ((error == error) /*NAN check*/
&& (error != inf) && (error != ninf))
{
if (GenerateMultipleEvents)
{
// --------- Multiple events per bin ----------
double errorSquared = error * error;
// Find how many events to fake
double val = weight / E[i];
val *= val;
// Convert to int with slight rounding up. This is to avoid rounding errors
int numEvents = int(val + 0.2);
if (numEvents < 1) numEvents = 1;
if (numEvents > MaxEventsPerBin) numEvents = MaxEventsPerBin;
// Scale the weight and error for each
weight /= numEvents;
errorSquared /= numEvents;

// Spread the TOF. e.g. 2 events = 0.25, 0.75.
double tofStep = (X[i+1] - X[i]) / (numEvents);
for (size_t j=0; j<size_t(numEvents); j++)
{
double tof = X[i] + tofStep * (0.5 + double(j));
// Create and add the event
el.addEventQuickly( WeightedEventNoTime(tof, weight, errorSquared) );
}
}
else
{
// --------- Single event per bin ----------
// TOF = midpoint of the bin
double tof = (X[i] + X[i+1]) / 2.0;
// Error squared is carried in the event
double errorSquared = E[i];
errorSquared *= errorSquared;
// Create and add the event
el.addEventQuickly( WeightedEventNoTime(tof, weight, errorSquared) );
}
} // error is nont NAN or infinite
} // weight is non-zero, not NAN, and non-infinite
} // (each bin)

// Set the X binning parameters
el.setX( inSpec->ptrX() );
// Manually set that this is sorted by TOF, since it is. This will make it "threadSafe" in other algos.
el.setSortOrder( TOF_SORT );

// This method fills in the events
el.createFromHistogram(inSpec, GenerateZeros, GenerateMultipleEvents, MaxEventsPerBin);

prog.report("Converting");
PARALLEL_END_INTERUPT_REGION
Expand Down
38 changes: 37 additions & 1 deletion Code/Mantid/Framework/Algorithms/test/ConvertToEventWorkspaceTest.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -161,7 +161,7 @@ class ConvertToEventWorkspaceTest : public CxxTest::TestSuite
AnalysisDataService::Instance().remove(outWSName);
AnalysisDataService::Instance().remove(inWSName);
}


/// Workspace with infinity or NAN = don't create events there.
void test_with_nan_and_inf()
Expand Down Expand Up @@ -212,6 +212,42 @@ class ConvertToEventWorkspaceTest : public CxxTest::TestSuite
TS_ASSERT_EQUALS( outWS->getNumberEvents(), 1);
}



/// Create events for zero-weight bins
void test_GenerateZeros()
{
// Create the input
Workspace2D_sptr inWS = WorkspaceCreationHelper::Create2DWorkspace(1, 10);

// Clear the vector
inWS->dataY(0).assign(10, 0.0);

// Name of the output workspace.
std::string outWSName("ConvertToEventWorkspaceTest_OutputWS");

ConvertToEventWorkspace alg;
TS_ASSERT_THROWS_NOTHING( alg.initialize() )
TS_ASSERT( alg.isInitialized() )
TS_ASSERT_THROWS_NOTHING( alg.setProperty("InputWorkspace", inWS) );
TS_ASSERT_THROWS_NOTHING( alg.setProperty("GenerateMultipleEvents", true) );
TS_ASSERT_THROWS_NOTHING( alg.setProperty("GenerateZeros", true) );
TS_ASSERT_THROWS_NOTHING( alg.setPropertyValue("OutputWorkspace", outWSName) );
TS_ASSERT_THROWS_NOTHING( alg.execute(); );
TS_ASSERT( alg.isExecuted() );

// Retrieve the workspace from data service.
EventWorkspace_sptr outWS;
TS_ASSERT_THROWS_NOTHING( outWS = AnalysisDataService::Instance().retrieveWS<EventWorkspace>(outWSName) );
TS_ASSERT(outWS);
if (!outWS) return;

// Every bin is created, even though they were zeros
TS_ASSERT_EQUALS( outWS->getNumberEvents(), 10);
}



};


Expand Down
2 changes: 2 additions & 0 deletions Code/Mantid/Framework/DataObjects/inc/MantidDataObjects/EventList.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -132,6 +132,8 @@ class DLLExport EventList : public Mantid::API::IEventList

virtual ~EventList();

void createFromHistogram(const ISpectrum * spec, bool GenerateZeros, bool GenerateMultipleEvents, int MaxEventsPerBin);

EventList& operator=(const EventList&);

EventList& operator+=(const TofEvent& event);
Expand Down
91 changes: 91 additions & 0 deletions Code/Mantid/Framework/DataObjects/src/EventList.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -138,6 +138,97 @@ namespace DataObjects
//std::vector<TofEvent>().swap(events); //Trick to release the vector memory.
}


// --------------------------------------------------------------------------
/** Create an EventList from a histogram. This converts bins to weighted
* events.
* Any existing events are cleared.
*
* @param spec :: ISpectrum ptr to histogram data.
* @param GenerateZeros :: if true, generate event(s) for empty bins
* @param GenerateMultipleEvents :: if true, create several evenly-spaced fake events inside the bin
* @param MaxEventsPerBin :: max number of events to generate in one bin, if GenerateMultipleEvents
*/
void EventList::createFromHistogram(const ISpectrum * inSpec, bool GenerateZeros, bool GenerateMultipleEvents, int MaxEventsPerBin)
{
// Fresh start
this->clear(true);

// Cached values for later checks
double inf = std::numeric_limits<double>::infinity();
double ninf = -inf;

// Get the input histogram
const MantidVec & X = inSpec->readX();
const MantidVec & Y = inSpec->readY();
const MantidVec & E = inSpec->readE();
if (Y.size()+1 != X.size())
throw std::runtime_error("Expected a histogram (X vector should be 1 longer than the Y vector)");

// Copy detector IDs and spectra
this->copyInfoFrom( *inSpec );
// We need weights but have no way to set the time. So use weighted, no time
this->switchTo(WEIGHTED_NOTIME);

for (size_t i=0; i<X.size()-1; i++)
{
double weight = Y[i];
if ((weight != 0.0 || GenerateZeros)
&& (weight == weight) /*NAN check*/
&& (weight != inf) && (weight != ninf))
{
double error = E[i];
// Also check that the error is not a bad number
if ((error == error) /*NAN check*/
&& (error != inf) && (error != ninf))
{
if (GenerateMultipleEvents)
{
// --------- Multiple events per bin ----------
double errorSquared = error * error;
// Find how many events to fake
double val = weight / E[i];
val *= val;
// Convert to int with slight rounding up. This is to avoid rounding errors
int numEvents = int(val + 0.2);
if (numEvents < 1) numEvents = 1;
if (numEvents > MaxEventsPerBin) numEvents = MaxEventsPerBin;
// Scale the weight and error for each
weight /= numEvents;
errorSquared /= numEvents;

// Spread the TOF. e.g. 2 events = 0.25, 0.75.
double tofStep = (X[i+1] - X[i]) / (numEvents);
for (size_t j=0; j<size_t(numEvents); j++)
{
double tof = X[i] + tofStep * (0.5 + double(j));
// Create and add the event
this->addEventQuickly( WeightedEventNoTime(tof, weight, errorSquared) );
}
}
else
{
// --------- Single event per bin ----------
// TOF = midpoint of the bin
double tof = (X[i] + X[i+1]) / 2.0;
// Error squared is carried in the event
double errorSquared = E[i];
errorSquared *= errorSquared;
// Create and add the event
this->addEventQuickly( WeightedEventNoTime(tof, weight, errorSquared) );
}
} // error is nont NAN or infinite
} // weight is non-zero, not NAN, and non-infinite
} // (each bin)

// Set the X binning parameters
this->setX( inSpec->ptrX() );

// Manually set that this is sorted by TOF, since it is. This will make it "threadSafe" in other algos.
this->setSortOrder( TOF_SORT );
}


// --------------------------------------------------------------------------
// --- Operators -------------------------------------------------------------------

Expand Down
6 changes: 6 additions & 0 deletions Code/Mantid/Framework/MDEvents/inc/MantidMDEvents/ConvertToDiffractionMDWorkspace.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,8 +48,14 @@ namespace MDEvents
template <class T>
void convertEventList(int workspaceIndex);

void convertSpectrum(int workspaceIndex);

/// The input event workspace
DataObjects::EventWorkspace_sptr in_ws;

/// The input event workspace
DataObjects::EventWorkspace_sptr m_inMWS;

/// The output MDEventWorkspace<3>
MDEvents::MDEventWorkspace3Lean::sptr ws;
/// Do we clear events on the input during loading?
Expand Down
11 changes: 11 additions & 0 deletions Code/Mantid/Framework/MDEvents/src/ConvertToDiffractionMDWorkspace.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -128,6 +128,17 @@ namespace MDEvents
/// Our MDLeanEvent dimension
typedef MDLeanEvent<3> MDE;

//----------------------------------------------------------------------------------------------
/** Convert one spectrum (histogram) into an event list
* with one event per bin, and then pass through to
* another method to convert to 3D q-space and add it to the MDEventWorkspace
*
* @param workspaceIndex :: index into the workspace
*/
void ConvertToDiffractionMDWorkspace::convertSpectrum(int workspaceIndex)
{
// this->workspaceIndex
}

//----------------------------------------------------------------------------------------------
/** Convert an event list to 3D q-space and add it to the MDEventWorkspace
Expand Down

0 comments on commit e3edcdf

Please sign in to comment.