Skip to content

Commit

Permalink
refs #4201 A crude support for events
Browse files Browse the repository at this point in the history 
Using internal dynamic rebinning of the events itself (which may be incorrect)
  • Loading branch information
@abuts
abuts committed Jan 6, 2012
1 parent 087b2a8 commit 1312946
Show file tree
Hide file tree
Showing 6 changed files with 298 additions and 58 deletions.
10 changes: 8 additions & 2 deletions Code/Mantid/Framework/MDAlgorithms/inc/MantidMDAlgorithms/ConvertToMDEvents.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -92,8 +92,14 @@ namespace MDAlgorithms
EventWSType, //< Event worskapce
NInWSTypes
};
// way to treat the X-coorinate in the workspace:
enum XCoordType
{
Histohram, // typical for Matrix workspace -- deploys central average 0.5(X[i]+X[i+1]); other types of averaging are possible if needed
Axis // typical for events
};
/// predefenition of the class, which does all coordinate transformations, Linux compilers need this.
template<Q_state Q, AnalMode MODE, CnvrtUnits CONV>
template<Q_state Q, AnalMode MODE, CnvrtUnits CONV,XCoordType XTYPE>
struct COORD_TRANSFORMER;

class DLLExport ConvertToMDEvents : public MDEvents::BoxControllerSettingsAlgorithm
Expand Down Expand Up @@ -202,7 +208,7 @@ namespace MDAlgorithms

/// shalow class which is invoked from processQND procedure and describes the transformation from workspace coordinates to target coordinates
/// presumably will be completely inlined
template<Q_state Q, AnalMode MODE, CnvrtUnits CONV>
template<Q_state Q, AnalMode MODE, CnvrtUnits CONV,XCoordType XTYPE>
friend struct COORD_TRANSFORMER;
/// helper class to orginize metaloop on various algorithm options
template<Q_state Q,size_t N_ALGORITHMS >
Expand Down
34 changes: 17 additions & 17 deletions Code/Mantid/Framework/MDAlgorithms/inc/MantidMDAlgorithms/ConvertToMDEventsCoordTransf.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,7 @@ namespace MDAlgorithms
*
* Generic template defines interface to 3 functions which perform these three steps.
*/
template<Q_state Q,AnalMode MODE,CnvrtUnits CONV>
template<Q_state Q,AnalMode MODE,CnvrtUnits CONV,XCoordType Type>
struct COORD_TRANSFORMER
{
COORD_TRANSFORMER(ConvertToMDEvents *){}
Expand Down Expand Up @@ -103,8 +103,8 @@ struct COORD_TRANSFORMER
// ----> NoQ
// NoQ,ANY_Mode -- no units conversion. This templates just copies the data into MD events and not doing any momentum transformations
//
template<AnalMode MODE,CnvrtUnits CONV>
struct COORD_TRANSFORMER<NoQ,MODE,CONV>
template<AnalMode MODE,CnvrtUnits CONV,XCoordType Type>
struct COORD_TRANSFORMER<NoQ,MODE,CONV,Type>
{
inline bool calcGenericVariables(std::vector<coord_t> &Coord, size_t nd)
{
Expand Down Expand Up @@ -149,7 +149,7 @@ struct COORD_TRANSFORMER<NoQ,MODE,CONV>
// pointer to MD workspace convertor
ConvertToMDEvents *pHost;
// class which would convert units
UNITS_CONVERSION<CONV> CONV_UNITS_FROM;
UNITS_CONVERSION<CONV,Type> CONV_UNITS_FROM;
};
//
////----------------------------------------------------------------------------------------------------------------------
Expand All @@ -173,8 +173,8 @@ inline double k_trans<Indir>(double Ei, double E_tr){

// -----> modQ
// modQ,Inelastic
template<AnalMode MODE,CnvrtUnits CONV>
struct COORD_TRANSFORMER<modQ,MODE,CONV>
template<AnalMode MODE,CnvrtUnits CONV,XCoordType Type>
struct COORD_TRANSFORMER<modQ,MODE,CONV,Type>
{
inline bool calcGenericVariables(std::vector<coord_t> &Coord, size_t nd)
{
Expand Down Expand Up @@ -247,11 +247,11 @@ struct COORD_TRANSFORMER<modQ,MODE,CONV>
// Calling Mantid algorithm
ConvertToMDEvents *pHost;
// class that performs untis conversion;
UNITS_CONVERSION<CONV> CONV_UNITS_FROM;
UNITS_CONVERSION<CONV,Type> CONV_UNITS_FROM;
};
// modQ,Elastic
template<CnvrtUnits CONV>
struct COORD_TRANSFORMER<modQ,Elastic,CONV>
template<CnvrtUnits CONV,XCoordType Type>
struct COORD_TRANSFORMER<modQ,Elastic,CONV,Type>
{
inline bool calcGenericVariables(std::vector<coord_t> &Coord, size_t nd)
{
Expand Down Expand Up @@ -314,13 +314,13 @@ struct COORD_TRANSFORMER<modQ,Elastic,CONV>
// Calling Mantid algorithm
ConvertToMDEvents *pHost;
// class that performs untis conversion;
UNITS_CONVERSION<CONV> CONV_UNITS_FROM;
UNITS_CONVERSION<CONV,Type> CONV_UNITS_FROM;
};


// Direct/Indirect tramsformatiom, this template describes 3D Q analysis mode.
template<AnalMode MODE,CnvrtUnits CONV>
struct COORD_TRANSFORMER<Q3D,MODE,CONV>
template<AnalMode MODE,CnvrtUnits CONV,XCoordType Type>
struct COORD_TRANSFORMER<Q3D,MODE,CONV,Type>
{
inline bool calcGenericVariables(std::vector<coord_t> &Coord, size_t nd)
{
Expand Down Expand Up @@ -386,16 +386,16 @@ struct COORD_TRANSFORMER<Q3D,MODE,CONV>
// Calling Mantid algorithm
ConvertToMDEvents *pHost;
// class that performs untis conversion;
UNITS_CONVERSION<CONV> CONV_UNITS_FROM;
UNITS_CONVERSION<CONV,Type> CONV_UNITS_FROM;
};

// Elastic
template<CnvrtUnits CONV>
struct COORD_TRANSFORMER<Q3D,Elastic,CONV>
template<CnvrtUnits CONV,XCoordType Type>
struct COORD_TRANSFORMER<Q3D,Elastic,CONV,Type>
{
inline bool calcGenericVariables(std::vector<coord_t> &Coord, size_t nd)
{
// tree inital coordinates came from workspace and all are interconnnected. All additional are defined by properties:
// three inital coordinates came from workspace and all are interconnnected. All additional are defined by properties:
if(!pHost->fillAddProperties(Coord,nd,3))return false;
//
rotMat = pHost->getTransfMatrix();
Expand Down Expand Up @@ -445,7 +445,7 @@ struct COORD_TRANSFORMER<Q3D,Elastic,CONV>
// pointer to the algoritm, which calls all these transformations
ConvertToMDEvents *pHost;
// class that performs untis conversion;
UNITS_CONVERSION<CONV> CONV_UNITS_FROM;
UNITS_CONVERSION<CONV,Type> CONV_UNITS_FROM;
};


Expand Down
236 changes: 231 additions & 5 deletions Code/Mantid/Framework/MDAlgorithms/inc/MantidMDAlgorithms/ConvertToMDEventsMethods.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,8 +51,8 @@ void ConvertToMDEvents::processQNDHWS()
pProg = std::auto_ptr<API::Progress>(new API::Progress(this,0.0,1.0,numSpec));

// initiate the templated class which does the conversion of workspace data into MD WS coordinates;
COORD_TRANSFORMER<Q,MODE,CONV> trn(this);
// one of the dimensions has to be X-ws dimension -> need to add check for that;
COORD_TRANSFORMER<Q,MODE,CONV,Histohram> trn(this);


// copy experiment info into target workspace
API::ExperimentInfo_sptr ExperimentInfo(inWS2D->cloneExperimentInfo());
Expand All @@ -65,18 +65,22 @@ void ConvertToMDEvents::processQNDHWS()
const size_t specSize = this->inWS2D->blocksize();
size_t nValidSpectra = det_loc.det_id.size();

std::vector<coord_t> Coord(n_dims); // coordinates for single event
// if any property dimension is outside of the data range requested, the job is done;
if(!trn.calcGenericVariables(Coord,n_dims))return;

// take at least bufSize amout of data in one run for efficiency
size_t buf_size = ((specSize>SPLIT_LEVEL)?specSize:SPLIT_LEVEL);
// allocate temporary buffer for MD Events data
std::vector<coord_t> allCoord(0); // MD events coordinates buffer
allCoord.reserve(n_dims*buf_size);
std::vector<coord_t> Coord(n_dims); // coordinates for single event

std::vector<float> sig_err(2*buf_size); // array for signal and error.
std::vector<uint16_t> run_index(buf_size); // Buffer run index for each event
std::vector<uint32_t> det_ids(buf_size); // Buffer of det Id-s for each event


if(!trn.calcGenericVariables(Coord,n_dims))return; // if any property dimension is outside of the data range requested

//External loop over the spectra:
for (int64_t i = 0; i < int64_t(nValidSpectra); ++i)
{
Expand Down Expand Up @@ -132,8 +136,230 @@ void ConvertToMDEvents::processQNDHWS()
template<Q_state Q, AnalMode MODE, CnvrtUnits CONV>
void ConvertToMDEvents::processQNDEWS()
{
}
DataObjects::EventWorkspace_const_sptr pEventWS = boost::static_pointer_cast<const DataObjects::EventWorkspace>(inWS2D);

size_t lastNumBoxes = this->pWSWrapper->pWorkspace()->getBoxController()->getTotalNumMDBoxes();
// counder for the number of events
size_t n_added_events(0);
// amount of work
const size_t numSpec = inWS2D->getNumberHistograms();
size_t nValidSpectra = det_loc.det_id.size();
// progress reporter
pProg = std::auto_ptr<API::Progress>(new API::Progress(this,0.0,1.0,numSpec));

// initiate the templated class which does the conversion of workspace data into MD WS coordinates;
COORD_TRANSFORMER<Q,MODE,CONV,Histohram> trn(this);


// copy experiment info into target workspace
API::ExperimentInfo_sptr ExperimentInfo(inWS2D->cloneExperimentInfo());
// run index;
uint16_t runIndex = this->pWSWrapper->pWorkspace()->addExperimentInfo(ExperimentInfo);
// number of dimesnions
size_t n_dims = this->pWSWrapper->nDimensions();
// coordinates for single event
std::vector<coord_t> Coord(n_dims);
// if any property dimension is outside of the data range requested, the job is done;
if(!trn.calcGenericVariables(Coord,n_dims))return;

// take at least bufSize amout of data in one run for efficiency
size_t buf_size = SPLIT_LEVEL;
// allocate temporary buffer for MD Events data
std::vector<coord_t> allCoord(0); // MD events coordinates buffer
allCoord.reserve(n_dims*buf_size);

std::vector<float> sig_err(2*buf_size); // array for signal and error.
std::vector<uint16_t> run_index(buf_size); // Buffer run index for each event
std::vector<uint32_t> det_ids(buf_size); // Buffer of det Id-s for each event


for (size_t wi=0; wi < nValidSpectra; wi++)
{
size_t ic = det_loc.detIDMap[wi];
int32_t det_id = det_loc.det_id[wi];

const EventList & el = pEventWS->getEventList(ic);
size_t numEvents = el.getNumberEvents();


const MantidVec& X = el.dataX();
const MantidVec& Signal = el.dataY();
const MantidVec& Error = el.dataE();
if(!trn.calcYDepCoordinates(Coord,ic))continue; // skip y outsize of the range;

//=> START INTERNAL LOOP OVER THE "TIME"
for (size_t j = 0; j < numEvents-1; ++j)
{

if(!trn.calcMatrixCoord(X,ic,j,Coord))continue; // skip ND outside the range
// ADD RESULTING EVENTS TO THE WORKSPACE
float ErrSq = float(Error[j]*Error[j]);

// coppy all data into data buffer for future transformation into events;
sig_err[2*n_added_events+0]=float(Signal[j]);
sig_err[2*n_added_events+1]=ErrSq;
run_index[n_added_events] = runIndex;
det_ids[n_added_events] = det_id;
allCoord.insert(allCoord.end(),Coord.begin(),Coord.end());

n_added_events++;
} // end spectra loop
if(n_added_events>=buf_size){
pWSWrapper->addMDData(sig_err,run_index,det_ids,allCoord,n_added_events);

n_added_events=0;
pProg->report(wi);
}
} // end detectors loop;

if(n_added_events>0){
pWSWrapper->addMDData(sig_err,run_index,det_ids,allCoord,n_added_events);

n_added_events=0;
}

pWSWrapper->refreshCache();



//// Equivalent of: this->convertEventList(wi);
//EventList & el = in_ws->getEventList(wi);

//// We want to bind to the right templated function, so we have to know the type of TofEvent contained in the EventList.
//boost::function<void ()> func;
//switch (el.getEventType())
//{
//case TOF:
// func = boost::bind(&ConvertToDiffractionMDWorkspace::convertEventList<TofEvent>, &*this, static_cast<int>(wi));
// break;
//case WEIGHTED:
// func = boost::bind(&ConvertToDiffractionMDWorkspace::convertEventList<WeightedEvent>, &*this, static_cast<int>(wi));
// break;
//case WEIGHTED_NOTIME:
// func = boost::bind(&ConvertToDiffractionMDWorkspace::convertEventList<WeightedEventNoTime>, &*this, static_cast<int>(wi));
// break;
//default:
// throw std::runtime_error("EventList had an unexpected data type!");
//}

//// Give this task to the scheduler
//double cost = double(el.getNumberEvents());
//ts->push( new FunctionTask( func, cost) );

//// Keep a running total of how many events we've added
//eventsAdded += el.getNumberEvents();
//if (bc->shouldSplitBoxes(eventsAdded, lastNumBoxes))
//{
// // Do all the adding tasks
// tp.joinAll();
// // Now do all the splitting tasks
// ws->splitAllIfNeeded(ts);
// tp.joinAll();

// // Count the new # of boxes.
// lastNumBoxes = ws->getBoxController()->getTotalNumMDBoxes();

// eventsAdded = 0;
//}


// tp.joinAll();
// // Do a final splitting of everything
// ws->splitAllIfNeeded(ts);
// tp.joinAll();
//
// // Recount totals at the end.
//// cputim.reset();
// ws->refreshCache();
// //TODO: Centroid in parallel, maybe?
// ws->getBox()->refreshCentroid(NULL);



}


// // This little dance makes the getting vector of events more general (since you can't overload by return type).
// typename std::vector<T> * events_ptr;
// getEventsFrom(el, events_ptr);
// typename std::vector<T> & events = *events_ptr;
//
// // Iterators to start/end
// typename std::vector<T>::iterator it = events.begin();
// typename std::vector<T>::iterator it_end = events.end();
//
// for (; it != it_end; it++)
// {
// // Get the wavenumber in ang^-1 using the previously calculated constant.
// double wavenumber = wavenumber_in_angstrom_times_tof_in_microsec / it->tof();
//
//
// // Q vector = K_final - K_initial = wavenumber * (output_direction - input_direction)
// coord_t center[3] = {Q_dir_x * wavenumber, Q_dir_y * wavenumber, Q_dir_z * wavenumber};
//
//
// // Push the MDLeanEvent with the same weight
// out_events.push_back( MDE(float(it->weight()), float(it->errorSquared()), center) );
//
// }


//template <class T>
//void ConvertToDiffractionMDWorkspace::convertEventList(int workspaceIndex)
// {
// EventList & el = in_ws->getEventList(workspaceIndex);
// size_t numEvents = el.getNumberEvents();
//
// // Get the position of the detector there.
// std::set<detid_t>& detectors = el.getDetectorIDs();
// if (detectors.size() > 0)
// {
// // The 3D MDEvents that will be added into the MDEventWorkspce
// std::vector<MDE> out_events;
// out_events.reserve( el.getNumberEvents() );
//
//
// // This little dance makes the getting vector of events more general (since you can't overload by return type).
// typename std::vector<T> * events_ptr;
// getEventsFrom(el, events_ptr);
// typename std::vector<T> & events = *events_ptr;
//
// // Iterators to start/end
// typename std::vector<T>::iterator it = events.begin();
// typename std::vector<T>::iterator it_end = events.end();
//
// for (; it != it_end; it++)
// {
// // Get the wavenumber in ang^-1 using the previously calculated constant.
// double wavenumber = wavenumber_in_angstrom_times_tof_in_microsec / it->tof();
//
//
// // Q vector = K_final - K_initial = wavenumber * (output_direction - input_direction)
// coord_t center[3] = {Q_dir_x * wavenumber, Q_dir_y * wavenumber, Q_dir_z * wavenumber};
//
//
// // Push the MDLeanEvent with the same weight
// out_events.push_back( MDE(float(it->weight()), float(it->errorSquared()), center) );
//
// }
//
// // Clear out the EventList to save memory
// if (ClearInputWorkspace)
// {
// // Track how much memory you cleared
// size_t memoryCleared = el.getMemorySize();
// // Clear it now
// el.clear();
// // For Linux with tcmalloc, make sure memory goes back, if you've cleared 200 Megs
// MemoryManager::Instance().releaseFreeMemoryIfAccumulated(memoryCleared, (size_t)2e8);
// }
//
// // Add them to the MDEW
// ws->addEvents(out_events);
// }
//
// }




Expand Down

0 comments on commit 1312946

Please sign in to comment.