Newer
Older
Campbell, Stuart
committed
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
Campbell, Stuart
committed
#include "MantidDataHandling/LoadEventNexus.h"
Campbell, Stuart
committed
#include "MantidGeometry/Instrument/CompAssembly.h"
#include "MantidKernel/ConfigService.h"
#include "MantidKernel/DateAndTime.h"
Janik Zikovsky
committed
#include "MantidKernel/ThreadPool.h"
#include "MantidKernel/FunctionTask.h"
Campbell, Stuart
committed
#include "MantidAPI/FileProperty.h"
#include "MantidKernel/UnitFactory.h"
#include "MantidKernel/Timer.h"
#include "MantidAPI/MemoryManager.h"
#include "MantidAPI/LoadAlgorithmFactory.h" // For the DECLARE_LOADALGORITHM macro
Gigg, Martyn Anthony
committed
#include "MantidAPI/SpectraAxis.h"
Campbell, Stuart
committed
#include <fstream>
#include <sstream>
#include <boost/algorithm/string/replace.hpp>
#include <Poco/File.h>
#include <Poco/Path.h>
Janik Zikovsky
committed
#include "MantidKernel/VisibleWhenProperty.h"
Gigg, Martyn Anthony
committed
#include "MantidKernel/EnabledWhenProperty.h"
Campbell, Stuart
committed
using std::endl;
using std::map;
using std::string;
using std::vector;
using namespace ::NeXus;
Campbell, Stuart
committed
using namespace Mantid::Geometry;
using namespace Mantid::DataObjects;
Janik Zikovsky
committed
using namespace Mantid::Kernel;
Campbell, Stuart
committed
namespace Mantid
{
Campbell, Stuart
committed
namespace DataHandling
Campbell, Stuart
committed
{
DECLARE_ALGORITHM(LoadEventNexus)
DECLARE_LOADALGORITHM(LoadEventNexus)
Janik Zikovsky
committed
Campbell, Stuart
committed
using namespace Kernel;
using namespace API;
using Geometry::Instrument;
Janik Zikovsky
committed
//===============================================================================================
//===============================================================================================
/** This task does the disk IO from loading the NXS file,
* and so will be on a disk IO mutex */
class ProcessBankData : public Task
{
public:
/**
*
* @param alg :: LoadEventNexus
* @param entry_name :: name of the bank
* @param pixelID_to_wi_map :: map pixel ID to Workspace Index
Janik Zikovsky
committed
* @param prog :: Progress reporter
* @param scheduler :: ThreadScheduler running this task
* @param event_id :: array with event IDs
* @param event_time_of_flight :: array with event TOFS
* @param numEvents :: how many events in the arrays
* @param startAt :: index of the first event from event_index
* @param event_index_ptr :: ptr to a vector of event index (length of # of pulses)
* @return
*/
ProcessBankData(LoadEventNexus * alg, std::string entry_name, detid2index_map * pixelID_to_wi_map,
Janik Zikovsky
committed
Progress * prog, ThreadScheduler * scheduler,
uint32_t * event_id, float * event_time_of_flight,
size_t numEvents, size_t startAt, std::vector<uint64_t> * event_index_ptr)
: Task(),
alg(alg), entry_name(entry_name), pixelID_to_wi_map(pixelID_to_wi_map), prog(prog), scheduler(scheduler),
Janik Zikovsky
committed
event_id(event_id), event_time_of_flight(event_time_of_flight), numEvents(numEvents), startAt(startAt),
event_index_ptr(event_index_ptr), event_index(*event_index_ptr)
{
// Cost is approximately proportional to the number of events to process.
Janik Zikovsky
committed
}
//----------------------------------------------------
// Run the data processing
void run()
{
//Local tof limits
double my_shortest_tof, my_longest_tof;
my_shortest_tof = static_cast<double>(std::numeric_limits<uint32_t>::max()) * 0.1;
my_longest_tof = 0.;
prog->report(entry_name + ": precount");
// ---- Pre-counting events per pixel ID ----
if (alg->precount)
{
Janik Zikovsky
committed
std::vector<size_t> counts;
// key = pixel ID, value = count
Gigg, Martyn Anthony
committed
counts.resize(alg->eventid_max+1);
Janik Zikovsky
committed
for (size_t i=0; i < numEvents; i++)
{
Janik Zikovsky
committed
detid_t thisId = detid_t(event_id[i]);
Gigg, Martyn Anthony
committed
if (thisId <= alg->eventid_max)
Janik Zikovsky
committed
counts[thisId]++;
Janik Zikovsky
committed
}
// Now we pre-allocate (reserve) the vectors of events in each pixel counted
Gigg, Martyn Anthony
committed
for (detid_t pixID = 0; pixID <= alg->eventid_max; pixID++)
Janik Zikovsky
committed
{
Janik Zikovsky
committed
if (counts[pixID] > 0)
{
//Find the the workspace index corresponding to that pixel ID
size_t wi = static_cast<size_t>((*pixelID_to_wi_map)[ pixID ]);
// Allocate it
alg->WS->getEventList(wi).reserve( counts[pixID] );
if (alg->getCancel()) break; // User cancellation
}
Janik Zikovsky
committed
}
}
// Check for cancelled algorithm
if (alg->getCancel())
{ delete [] event_id; delete [] event_time_of_flight; return; }
//Default pulse time (if none are found)
Mantid::Kernel::DateAndTime pulsetime;
// Index into the pulse array
int pulse_i = 0;
// And there are this many pulses
int numPulses = static_cast<int>(alg->pulseTimes.size());
if (numPulses > static_cast<int>(event_index.size()))
{
alg->getLogger().warning() << "Entry " << entry_name << "'s event_index vector is smaller than the proton_charge DAS log. This is inconsistent, so we cannot find pulse times for this entry.\n";
//This'll make the code skip looking for any pulse times.
pulse_i = numPulses + 1;
}
prog->report(entry_name + ": filling events");
Janik Zikovsky
committed
// The workspace
EventWorkspace_sptr WS = alg->WS;
// Will we need to compress?
bool compress = (alg->compressTolerance >= 0);
Janik Zikovsky
committed
// Which detector IDs were touched?
Gigg, Martyn Anthony
committed
std::vector<bool> usedDetIds(alg->eventid_max+1, false);
Janik Zikovsky
committed
Janik Zikovsky
committed
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
//Go through all events in the list
for (std::size_t i = 0; i < numEvents; i++)
{
//------ Find the pulse time for this event index ---------
if (pulse_i < numPulses-1)
{
bool breakOut = false;
//Go through event_index until you find where the index increases to encompass the current index. Your pulse = the one before.
while ( !((i+startAt >= event_index[pulse_i]) && (i+startAt < event_index[pulse_i+1])))
{
pulse_i++;
// Check once every new pulse if you need to cancel (checking on every event might slow things down more)
if (alg->getCancel()) breakOut = true;
if (pulse_i >= (numPulses-1))
break;
}
//Save the pulse time at this index for creating those events
pulsetime = alg->pulseTimes[pulse_i];
// Flag to break out of the event loop with using goto ;)
if (breakOut)
break;
}
//Create the tofevent
double tof = static_cast<double>( event_time_of_flight[i] );
if ((tof >= alg->filter_tof_min) && (tof <= alg->filter_tof_max))
{
//The event TOF passes the filter.
TofEvent event(tof, pulsetime);
Janik Zikovsky
committed
// We cached a pointer to the vector<tofEvent> -> so retrieve it and add the event
detid_t detId = event_id[i];
Gigg, Martyn Anthony
committed
if (detId <= alg->eventid_max)
Janik Zikovsky
committed
{
alg->eventVectors[detId]->push_back( event );
Janik Zikovsky
committed
Janik Zikovsky
committed
// //Find the the workspace index corresponding to that pixel ID
// size_t wi = static_cast<size_t>((*pixelID_to_wi_map)[event_id[i]]);
// // Add it to the list at that workspace index
// WS->getEventList(wi).addEventQuickly( event );
Janik Zikovsky
committed
Janik Zikovsky
committed
//Local tof limits
if (tof < my_shortest_tof) { my_shortest_tof = tof;}
if (tof > my_longest_tof) { my_longest_tof = tof;}
// Track all the touched wi
if (compress)
{
usedDetIds[detId] = true;
}
} // valid detector IDs
Janik Zikovsky
committed
Janik Zikovsky
committed
}
} //(for each event)
Janik Zikovsky
committed
//------------ Compress Events ------------------
if (compress)
{
Janik Zikovsky
committed
// Do it on all the detector IDs we touched
Gigg, Martyn Anthony
committed
for (detid_t pixID = 0; pixID <= alg->eventid_max; pixID++)
Janik Zikovsky
committed
{
Janik Zikovsky
committed
if (usedDetIds[pixID])
{
//Find the the workspace index corresponding to that pixel ID
size_t wi = static_cast<size_t>((*pixelID_to_wi_map)[ pixID ]);
EventList * el = WS->getEventListPtr(wi);
el->compressEvents(alg->compressTolerance, el);
}
Janik Zikovsky
committed
}
}
Janik Zikovsky
committed
//Join back up the tof limits to the global ones
PARALLEL_CRITICAL(tof_limits)
{
//This is not thread safe, so only one thread at a time runs this.
if (my_shortest_tof < alg->shortest_tof) { alg->shortest_tof = my_shortest_tof;}
if (my_longest_tof > alg->longest_tof ) { alg->longest_tof = my_longest_tof;}
}
// Free Memory
delete [] event_id;
delete [] event_time_of_flight;
delete event_index_ptr;
Janik Zikovsky
committed
// For Linux with tcmalloc, make sure memory goes back;
// but don't call if more than 15% of memory is still available, since that slows down the loading.
MemoryManager::Instance().releaseFreeMemoryIfAbove(0.85);
Janik Zikovsky
committed
}
private:
Janik Zikovsky
committed
LoadEventNexus * alg;
Janik Zikovsky
committed
std::string entry_name;
detid2index_map * pixelID_to_wi_map;
Janik Zikovsky
committed
Progress * prog;
Janik Zikovsky
committed
ThreadScheduler * scheduler;
Janik Zikovsky
committed
uint32_t * event_id;
Janik Zikovsky
committed
float * event_time_of_flight;
Janik Zikovsky
committed
size_t numEvents;
Janik Zikovsky
committed
size_t startAt;
/// ptr to a vector of event index vs time (length of # of pulses)
Janik Zikovsky
committed
std::vector<uint64_t> * event_index_ptr;
Janik Zikovsky
committed
std::vector<uint64_t> & event_index;
};
//===============================================================================================
//===============================================================================================
/** This task does the disk IO from loading the NXS file,
* and so will be on a disk IO mutex */
class LoadBankFromDiskTask : public Task
{
public:
//---------------------------------------------------------------------------------------------------
/** Constructor
*
* @param top_entry_name :: The pathname of the top level NXentry to use
Janik Zikovsky
committed
* @param entry_name :: The pathname of the bank to load
Gigg, Martyn Anthony
committed
* @param entry_type :: The classtype of the entry to load
* @param pixelID_to_wi_map :: a map where key = pixelID and value = the workpsace index to use.
Janik Zikovsky
committed
* @param prog :: an optional Progress object
* @param ioMutex :: a mutex shared for all Disk I-O tasks
* @param scheduler :: the ThreadScheduler that runs this task.
*/
Gigg, Martyn Anthony
committed
LoadBankFromDiskTask(LoadEventNexus * alg, const std::string& top_entry_name, const std::string& entry_name, const std::string & entry_type, detid2index_map * pixelID_to_wi_map,
Janik Zikovsky
committed
Progress * prog, Mutex * ioMutex, ThreadScheduler * scheduler)
: Task(),
Gigg, Martyn Anthony
committed
alg(alg), top_entry_name(top_entry_name), entry_name(entry_name), entry_type(entry_type),
pixelID_to_wi_map(pixelID_to_wi_map), prog(prog), scheduler(scheduler)
Janik Zikovsky
committed
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
{
setMutex(ioMutex);
}
//---------------------------------------------------------------------------------------------------
void run()
{
//The vectors we will be filling
std::vector<uint64_t> * event_index_ptr = new std::vector<uint64_t>();
std::vector<uint64_t> & event_index = *event_index_ptr;
// These give the limits in each file as to which events we actually load (when filtering by time).
std::vector<int> load_start(1); //TODO: Should this be size_t?
std::vector<int> load_size(1);
// Data arrays
uint32_t * event_id = NULL;
float * event_time_of_flight = NULL;
bool loadError = false ;
prog->report(entry_name + ": load from disk");
// Open the file
::NeXus::File file(alg->m_filename);
try
{
file.openGroup(top_entry_name, "NXentry");
Janik Zikovsky
committed
//Open the bankN_event group
Gigg, Martyn Anthony
committed
file.openGroup(entry_name, entry_type);
Janik Zikovsky
committed
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
// Get the event_index (a list of size of # of pulses giving the index in the event list for that pulse)
file.openData("event_index");
//Must be uint64
if (file.getInfo().type == ::NeXus::UINT64)
file.getData(event_index);
else
{
alg->getLogger().warning() << "Entry " << entry_name << "'s event_index field is not UINT64! It will be skipped.\n";
loadError = true;
}
file.closeData();
// Look for the sign that the bank is empty
if (event_index.size()==1)
{
if (event_index[0] == 0)
{
//One entry, only zero. This means NO events in this bank.
loadError = true;
alg->getLogger().debug() << "Bank " << entry_name << " is empty.\n";
}
}
if (event_index.size() != alg->pulseTimes.size())
{
alg->getLogger().debug() << "Bank " << entry_name << " has a mismatch between the number of event_index entries and the number of pulse times.\n";
Janik Zikovsky
committed
}
if (!loadError)
{
bool old_nexus_file_names = false;
// Get the list of pixel ID's
try
{
file.openData("event_id");
}
catch (::NeXus::Exception& )
{
//Older files (before Nov 5, 2010) used this field.
file.openData("event_pixel_id");
old_nexus_file_names = true;
}
// By default, use all available indices
Janik Zikovsky
committed
::NeXus::Info id_info = file.getInfo();
Janik Zikovsky
committed
//TODO: Handle the time filtering by changing the start/end offsets.
for (size_t i=0; i < alg->pulseTimes.size(); i++)
{
if (alg->pulseTimes[i] >= alg->filter_time_start)
{
start_event = event_index[i];
break; // stop looking
}
}
if (start_event > static_cast<size_t>(id_info.dims[0]))
{
// For bad file around SEQ_7872, Jul 15, 2011, Janik Zikovsky
alg->getLogger().information() << this->entry_name << "'s field 'event_index' seem to be invalid (> than the number of events in the bank). Filtering by time ignored.\n";
start_event = 0;
stop_event = static_cast<size_t>(id_info.dims[0]);
}
else
Janik Zikovsky
committed
{
for (size_t i=0; i < alg->pulseTimes.size(); i++)
Janik Zikovsky
committed
{
if (alg->pulseTimes[i] > alg->filter_time_stop)
{
stop_event = event_index[i];
break;
}
Janik Zikovsky
committed
}
}
// Make sure it is within range
Janik Zikovsky
committed
stop_event = id_info.dims[0];
alg->getLogger().debug() << entry_name << ": start_event " << start_event << " stop_event "<< stop_event << std::endl;
// These are the arguments to getSlab()
load_start[0] = static_cast<int>(start_event);
load_size[0] = static_cast<int>(stop_event - start_event);
Janik Zikovsky
committed
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
if ((load_size[0] > 0) && (load_start[0]>=0) )
{
// Now we allocate the required arrays
event_id = new uint32_t[load_size[0]];
event_time_of_flight = new float[load_size[0]];
// Check that the required space is there in the file.
if (id_info.dims[0] < load_size[0]+load_start[0])
{
alg->getLogger().warning() << "Entry " << entry_name << "'s event_id field is too small (" << id_info.dims[0]
<< ") to load the desired data size (" << load_size[0]+load_start[0] << ").\n";
loadError = true;
}
if (alg->getCancel()) loadError = true; //To allow cancelling the algorithm
if (!loadError)
{
//Must be uint32
if (id_info.type == ::NeXus::UINT32)
file.getSlab(event_id, load_start, load_size);
else
{
alg->getLogger().warning() << "Entry " << entry_name << "'s event_id field is not UINT32! It will be skipped.\n";
loadError = true;
}
file.closeData();
}
if (alg->getCancel()) loadError = true; //To allow cancelling the algorithm
if (!loadError)
{
// Get the list of event_time_of_flight's
if (!old_nexus_file_names)
file.openData("event_time_offset");
else
file.openData("event_time_of_flight");
// Check that the required space is there in the file.
::NeXus::Info tof_info = file.getInfo();
if (tof_info.dims[0] < load_size[0]+load_start[0])
{
alg->getLogger().warning() << "Entry " << entry_name << "'s event_time_offset field is too small to load the desired data.\n";
loadError = true;
}
//Check that the type is what it is supposed to be
if (tof_info.type == ::NeXus::FLOAT32)
file.getSlab(event_time_of_flight, load_start, load_size);
else
{
alg->getLogger().warning() << "Entry " << entry_name << "'s event_time_offset field is not FLOAT32! It will be skipped.\n";
loadError = true;
}
if (!loadError)
{
std::string units;
file.getAttr("units", units);
if (units != "microsecond")
{
alg->getLogger().warning() << "Entry " << entry_name << "'s event_time_offset field's units are not microsecond. It will be skipped.\n";
loadError = true;
}
file.closeData();
} //no error
} //no error
} // Size is at least 1
else
{
// Found a size that was 0 or less; stop processign
loadError=true;
}
} //no error
} // try block
Janik Zikovsky
committed
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
{
alg->getLogger().error() << "Error while loading bank " << entry_name << ":" << std::endl;
alg->getLogger().error() << e.what() << std::endl;
loadError = true;
}
catch (...)
{
alg->getLogger().error() << "Unspecified error while loading bank " << entry_name << std::endl;
loadError = true;
}
//Close up the file even if errors occured.
file.closeGroup();
file.close();
//Abort if anything failed
if (loadError)
{
prog->reportIncrement(2, entry_name + ": skipping");
delete [] event_id;
delete [] event_time_of_flight;
delete event_index_ptr;
return;
}
// No error? Launch a new task to process that data.
size_t numEvents = load_size[0];
size_t startAt = load_start[0];
ProcessBankData * newTask = new ProcessBankData(alg, entry_name,pixelID_to_wi_map,prog,scheduler,
Janik Zikovsky
committed
event_id,event_time_of_flight, numEvents, startAt, event_index_ptr);
scheduler->push(newTask);
}
private:
Janik Zikovsky
committed
LoadEventNexus * alg;
/// NXS name for top level NXentry
std::string top_entry_name;
Janik Zikovsky
committed
std::string entry_name;
Gigg, Martyn Anthony
committed
/// NXS type
std::string entry_type;
detid2index_map * pixelID_to_wi_map;
Janik Zikovsky
committed
Progress * prog;
Janik Zikovsky
committed
ThreadScheduler * scheduler;
};
//===============================================================================================
//===============================================================================================
Campbell, Stuart
committed
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
/// Empty default constructor
LoadEventNexus::LoadEventNexus() : IDataFileChecker()
{}
/**
* Do a quick file type check by looking at the first 100 bytes of the file
* @param filePath :: path of the file including name.
* @param nread :: no.of bytes read
* @param header :: The first 100 bytes of the file as a union
* @return true if the given file is of type which can be loaded by this algorithm
*/
bool LoadEventNexus::quickFileCheck(const std::string& filePath,size_t nread, const file_header& header)
{
std::string ext = this->extension(filePath);
// If the extension is nxs then give it a go
if( ext.compare("nxs") == 0 ) return true;
// If not then let's see if it is a HDF file by checking for the magic cookie
if ( nread >= sizeof(int32_t) && (ntohl(header.four_bytes) == g_hdf_cookie) ) return true;
return false;
}
/**
* Checks the file by opening it and reading few lines
* @param filePath :: name of the file inluding its path
* @return an integer value how much this algorithm can load the file
*/
int LoadEventNexus::fileCheck(const std::string& filePath)
{
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
int confidence(0);
typedef std::map<std::string,std::string> string_map_t;
try
{
string_map_t::const_iterator it;
::NeXus::File file = ::NeXus::File(filePath);
string_map_t entries = file.getEntries();
for(string_map_t::const_iterator it = entries.begin(); it != entries.end(); ++it)
{
if ( ((it->first == "entry") || (it->first == "raw_data_1")) && (it->second == "NXentry") )
{
file.openGroup(it->first, it->second);
string_map_t entries2 = file.getEntries();
for(string_map_t::const_iterator it2 = entries2.begin(); it2 != entries2.end(); ++it2)
{
if (it2->second == "NXevent_data")
{
confidence = 80;
}
}
file.closeGroup();
}
}
}
catch(::NeXus::Exception&)
{
}
return confidence;
Campbell, Stuart
committed
}
/// Sets documentation strings for this algorithm
void LoadEventNexus::initDocs()
{
this->setWikiSummary("Loads Event NeXus files (produced by the SNS) and stores it in an [[EventWorkspace]]. Optionally, you can filter out events falling outside a range of times-of-flight and/or a time interval. ");
this->setOptionalMessage("Loads Event NeXus files (produced by the SNS) and stores it in an EventWorkspace. Optionally, you can filter out events falling outside a range of times-of-flight and/or a time interval.");
}
Campbell, Stuart
committed
/// Initialisation method.
void LoadEventNexus::init()
{
std::vector<std::string> exts;
exts.push_back("_event.nxs");
exts.push_back(".nxs");
this->declareProperty(new FileProperty("Filename", "", FileProperty::Load, exts),
"The name of the Event NeXus file to read, including its full or relative path. \n"
"The file name is typically of the form INST_####_event.nxs (N.B. case sensitive if running on Linux)." );
Campbell, Stuart
committed
this->declareProperty(
new WorkspaceProperty<IEventWorkspace>("OutputWorkspace", "", Direction::Output),
"The name of the output EventWorkspace in which to load the EventNexus file." );
declareProperty(
new PropertyWithValue<double>("FilterByTofMin", EMPTY_DBL(), Direction::Input),
Campbell, Stuart
committed
"Optional: To exclude events that do not fall within a range of times-of-flight.\n"\
"This is the minimum accepted value in microseconds. Keep blank to load all events." );
Campbell, Stuart
committed
declareProperty(
new PropertyWithValue<double>("FilterByTofMax", EMPTY_DBL(), Direction::Input),
Campbell, Stuart
committed
"Optional: To exclude events that do not fall within a range of times-of-flight.\n"\
"This is the maximum accepted value in microseconds. Keep blank to load all events." );
Campbell, Stuart
committed
declareProperty(
new PropertyWithValue<double>("FilterByTimeStart", EMPTY_DBL(), Direction::Input),
Campbell, Stuart
committed
"Optional: To only include events after the provided start time, in seconds (relative to the start of the run).");
declareProperty(
new PropertyWithValue<double>("FilterByTimeStop", EMPTY_DBL(), Direction::Input),
Campbell, Stuart
committed
"Optional: To only include events before the provided stop time, in seconds (relative to the start of the run).");
Janik Zikovsky
committed
std::string grp1 = "Filter Events";
setPropertyGroup("FilterByTofMin", grp1);
setPropertyGroup("FilterByTofMax", grp1);
setPropertyGroup("FilterByTimeStart", grp1);
setPropertyGroup("FilterByTimeStop", grp1);
Janik Zikovsky
committed
Campbell, Stuart
committed
declareProperty(
new PropertyWithValue<string>("BankName", "", Direction::Input),
"Optional: To only include events from one bank. Any bank whose name does not match the given string will have no events.");
declareProperty(
new PropertyWithValue<bool>("SingleBankPixelsOnly", true, Direction::Input),
"Optional: Only applies if you specified a single bank to load with BankName.\n"
"Only pixels in the specified bank will be created if true; all of the instrument's pixels will be created otherwise.");
Janik Zikovsky
committed
setPropertySettings("SingleBankPixelsOnly", new VisibleWhenProperty(this, "BankName", IS_NOT_DEFAULT) );
Campbell, Stuart
committed
Janik Zikovsky
committed
std::string grp2 = "Loading a Single Bank";
setPropertyGroup("BankName", grp2);
setPropertyGroup("SingleBankPixelsOnly", grp2);
Campbell, Stuart
committed
declareProperty(
new PropertyWithValue<bool>("Precount", false, Direction::Input),
"Pre-count the number of events in each pixel before allocating memory (optional, default False). \n"
"This can significantly reduce memory use and memory fragmentation; it may also speed up loading.");
Janik Zikovsky
committed
declareProperty(
new PropertyWithValue<double>("CompressTolerance", -1.0, Direction::Input),
"Run CompressEvents while loading (optional, leave blank or negative to not do). \n"
"This specified the tolerance to use (in microseconds) when compressing.");
Janik Zikovsky
committed
std::string grp3 = "Reduce Memory Use";
setPropertyGroup("Precount", grp3);
setPropertyGroup("CompressTolerance", grp3);
declareProperty(
new PropertyWithValue<bool>("LoadMonitors", false, Direction::Input),
"Load the monitors from the file (optional, default False).");
Gigg, Martyn Anthony
committed
declareProperty(new PropertyWithValue<bool>("MonitorsAsEvents", false, Direction::Input),
"If present, load the monitors as events.\nWARNING: WILL SIGNIFICANTLY INCREASE MEMORY USAGE (optional, default False). \n");
Gigg, Martyn Anthony
committed
declareProperty(
new PropertyWithValue<double>("FilterMonByTofMin", EMPTY_DBL(), Direction::Input),
"Optional: To exclude events from monitors that do not fall within a range of times-of-flight.\n"\
"This is the minimum accepted value in microseconds." );
declareProperty(
new PropertyWithValue<double>("FilterMonByTofMax", EMPTY_DBL(), Direction::Input),
"Optional: To exclude events from monitors that do not fall within a range of times-of-flight.\n"\
"This is the maximum accepted value in microseconds." );
declareProperty(
new PropertyWithValue<double>("FilterMonByTimeStart", EMPTY_DBL(), Direction::Input),
"Optional: To only include events from monitors after the provided start time, in seconds (relative to the start of the run).");
declareProperty(
new PropertyWithValue<double>("FilterMonByTimeStop", EMPTY_DBL(), Direction::Input),
"Optional: To only include events from monitors before the provided stop time, in seconds (relative to the start of the run).");
Janik Zikovsky
committed
setPropertySettings("MonitorsAsEvents", new VisibleWhenProperty(this, "LoadMonitors", IS_EQUAL_TO, "1") );
Gigg, Martyn Anthony
committed
IPropertySettings *asEventsIsOn = new VisibleWhenProperty(this, "MonitorsAsEvents", IS_EQUAL_TO, "1");
setPropertySettings("FilterMonByTofMin", asEventsIsOn);
setPropertySettings("FilterMonByTofMax", asEventsIsOn->clone());
setPropertySettings("FilterMonByTimeStart", asEventsIsOn->clone());
setPropertySettings("FilterMonByTimeStop", asEventsIsOn->clone());
Janik Zikovsky
committed
std::string grp4 = "Monitors";
setPropertyGroup("LoadMonitors", grp4);
setPropertyGroup("MonitorsAsEvents", grp4);
Gigg, Martyn Anthony
committed
setPropertyGroup("FilterMonByTofMin", grp4);
setPropertyGroup("FilterMonByTofMax", grp4);
setPropertyGroup("FilterMonByTimeStart", grp4);
setPropertyGroup("FilterMonByTimeStop", grp4);
Janik Zikovsky
committed
Campbell, Stuart
committed
}
/// set the name of the top level NXentry m_top_entry_name
void LoadEventNexus::setTopEntryName()
{
typedef std::map<std::string,std::string> string_map_t;
try
{
string_map_t::const_iterator it;
::NeXus::File file = ::NeXus::File(m_filename);
string_map_t entries = file.getEntries();
for (it = entries.begin(); it != entries.end(); ++it)
{
if ( ((it->first == "entry") || (it->first == "raw_data_1")) && (it->second == "NXentry") )
{
m_top_entry_name = it->first;
break;
}
}
}
catch(const std::exception&)
{
g_log.error() << "Unable to determine name of top level NXentry - assuming \"entry\"." << std::endl;
m_top_entry_name = "entry";
}
Campbell, Stuart
committed
//------------------------------------------------------------------------------------------------
/** Executes the algorithm. Reading in the file and creating and populating
* the output workspace
*/
void LoadEventNexus::exec()
{
// Retrieve the filename from the properties
m_filename = getPropertyValue("Filename");
precount = getProperty("Precount");
Janik Zikovsky
committed
compressTolerance = getProperty("CompressTolerance");
Campbell, Stuart
committed
loadlogs = true;
Gigg, Martyn Anthony
committed
// //Get the limits to the filter
// filter_tof_min = getProperty("FilterByTofMin");
// filter_tof_max = getProperty("FilterByTofMax");
// if ( (filter_tof_min == EMPTY_DBL()) || (filter_tof_max == EMPTY_DBL()))
// {
// //Nothing specified. Include everything
// filter_tof_min = -1e20;
// filter_tof_max = +1e20;
// }
// else if ( (filter_tof_min != EMPTY_DBL()) || (filter_tof_max != EMPTY_DBL()))
// {
// //Both specified. Keep these values
// }
// else
// throw std::invalid_argument("You must specify both the min and max of time of flight to filter, or neither!");
Campbell, Stuart
committed
// Check to see if the monitors need to be loaded later
bool load_monitors = this->getProperty("LoadMonitors");
setTopEntryName();
Campbell, Stuart
committed
//Initialize progress reporting.
Janik Zikovsky
committed
int reports = 3;
Campbell, Stuart
committed
if (load_monitors)
reports++;
Progress prog(this,0.0,0.3, reports);
Gigg, Martyn Anthony
committed
// Load the detector events
WS = createEmptyEventWorkspace(); // Algorithm currently relies on an object-level workspace ptr
loadEvents(&prog, false); // Do not load monitor blocks
//Save output
this->setProperty<IEventWorkspace_sptr>("OutputWorkspace", WS);
// Load the monitors
if (load_monitors)
{
prog.report("Loading monitors");
const bool eventMonitors = getProperty("MonitorsAsEvents");
if( eventMonitors && this->hasEventMonitors() )
{
WS = createEmptyEventWorkspace(); // Algorithm currently relies on an object-level workspace ptr
Gigg, Martyn Anthony
committed
// // Set filter member variable attributes for
// filter_tof_min = getProperty("FilterMonByTofMin");
// filter_tof_max = getProperty("FilterMonByTofMax");
// if ( (filter_tof_min == EMPTY_DBL()) || (filter_tof_max == EMPTY_DBL()))
// {
// //Nothing specified. Include everything
// filter_tof_min = -1e20;
// filter_tof_max = +1e20;
// }
// else if ( (filter_tof_min != EMPTY_DBL()) || (filter_tof_max != EMPTY_DBL()))
// {
// //Both specified. Keep these values
// }
// else
// throw std::invalid_argument("You must specify both the min and max or neither of time of flight to filter the monitor");
// Perform the load
Gigg, Martyn Anthony
committed
loadEvents(&prog, true);
std::string mon_wsname = this->getProperty("OutputWorkspace");
mon_wsname.append("_monitors");
this->declareProperty(new WorkspaceProperty<IEventWorkspace>
("MonitorWorkspace", mon_wsname, Direction::Output), "Monitors from the Event NeXus file");
this->setProperty<IEventWorkspace_sptr>("MonitorWorkspace", WS);
}
else
{
this->runLoadMonitors();
}
}
// Clear any large vectors to free up memory.
this->pulseTimes.clear();
// Some memory feels like it sticks around (on Linux). Free it.
MemoryManager::Instance().releaseFreeMemory();
Campbell, Stuart
committed
Gigg, Martyn Anthony
committed
return;
}
Janik Zikovsky
committed
//-----------------------------------------------------------------------------
/** Generate a look-up table where the index = the pixel ID of an event
* and the value = a pointer to the EventList in the workspace
*/
void LoadEventNexus::makeMapToEventLists()
{
Gigg, Martyn Anthony
committed
eventid_max = 0; // seems like a safe lower bound
if( this->event_id_is_spec )
{
// Maximum possible spectrum number
detid2index_map::const_iterator it = pixelID_to_wi_map->end();
--it;
eventid_max = it->first;
}
else
{
// We want to pad out empty pixels.
detid2det_map detector_map;
WS->getInstrument()->getDetectors(detector_map);
Janik Zikovsky
committed
Gigg, Martyn Anthony
committed
// determine maximum pixel id
detid2det_map::iterator it;
for (it = detector_map.begin(); it != detector_map.end(); it++)
{
if (it->first > eventid_max) eventid_max = it->first;
}
}
Janik Zikovsky
committed
// Make an array where index = pixel ID
// Set the value to the 0th workspace index by default
Gigg, Martyn Anthony
committed
eventVectors.resize(eventid_max+1, &WS->getEventList(0).getEvents() );
for (detid_t j=0; j<eventid_max+1; j++)
Janik Zikovsky
committed
{
size_t wi = (*pixelID_to_wi_map)[j];
// Save a POINTER to the vector<tofEvent>
eventVectors[j] = &WS->getEventList(wi).getEvents();
}
}
Gigg, Martyn Anthony
committed
//-----------------------------------------------------------------------------
/**
* Load events from the file
* @param prop :: A pointer to the progress reporting object
* @param monitors :: If true the events from the monitors are loaded and not the main banks
*/
void LoadEventNexus::loadEvents(API::Progress * const prog, const bool monitors)
{
Gigg, Martyn Anthony
committed
// Get the time filters
setTimeFilters(monitors);
// The run_start will be loaded from the pulse times.
DateAndTime run_start(0,0);
Campbell, Stuart
committed
if (loadlogs)
{
Gigg, Martyn Anthony
committed
prog->doReport("Loading DAS logs");
Janik Zikovsky
committed
runLoadNexusLogs(m_filename, WS, pulseTimes, this);
run_start = WS->getFirstPulseTime();
Campbell, Stuart
committed
}
else
{
g_log.information() << "Skipping the loading of sample logs!" << endl;
}
Gigg, Martyn Anthony
committed
Campbell, Stuart
committed
//Load the instrument
Janik Zikovsky
committed
prog->report("Loading instrument");
instrument_loaded_correctly = runLoadInstrument(m_filename, WS, m_top_entry_name, this);
Campbell, Stuart
committed
if (!this->instrument_loaded_correctly)
throw std::runtime_error("Instrument was not initialized correctly! Loading cannot continue.");
// top level file information
::NeXus::File file(m_filename);
//Start with the base entry
file.openGroup(m_top_entry_name, "NXentry");
Campbell, Stuart
committed
//Now we want to go through all the bankN_event entries
Gigg, Martyn Anthony
committed
vector<string> bankNames;
Campbell, Stuart
committed
map<string, string> entries = file.getEntries();
map<string,string>::const_iterator it = entries.begin();
Gigg, Martyn Anthony
committed
std::string classType = monitors ? "NXmonitor" : "NXevent_data";
Campbell, Stuart
committed
for (; it != entries.end(); it++)
{
std::string entry_name(it->first);
std::string entry_class(it->second);
Gigg, Martyn Anthony
committed
if ( entry_class == classType )
Campbell, Stuart
committed
{
bankNames.push_back( entry_name );
}
}
//Close up the file
file.closeGroup();
file.close();
Gigg, Martyn Anthony
committed
// --------- Loading only one bank ? ----------------------------------
Campbell, Stuart
committed
std::string onebank = getProperty("BankName");
bool doOneBank = (onebank != "");
bool SingleBankPixelsOnly = getProperty("SingleBankPixelsOnly");
Gigg, Martyn Anthony
committed
if (doOneBank && !monitors)
Campbell, Stuart
committed
{
bool foundIt = false;
for (std::vector<string>::iterator it=bankNames.begin(); it!= bankNames.end(); it++)
Gigg, Martyn Anthony
committed
{
if (*it == ( onebank + "_events") )
Gigg, Martyn Anthony
committed
{
foundIt = true;
Gigg, Martyn Anthony
committed
break;
}
}
if (!foundIt)
{
throw std::invalid_argument("No entry named '" + onebank + "_events'" + " was found in the .NXS file.\n");
}
Campbell, Stuart
committed
bankNames.clear();
bankNames.push_back( onebank + "_events" );
Gigg, Martyn Anthony
committed
if( !SingleBankPixelsOnly ) onebank = ""; // Marker to load all pixels
}
else
{
onebank = "";
Campbell, Stuart
committed
}
// Delete the output workspace name if it existed
std::string outName = getPropertyValue("OutputWorkspace");
if (AnalysisDataService::Instance().doesExist(outName))
AnalysisDataService::Instance().remove( outName );
Campbell, Stuart
committed
Gigg, Martyn Anthony
committed
prog->report("Initializing all pixels");
Campbell, Stuart
committed
//----------------- Pad Empty Pixels -------------------------------
Gigg, Martyn Anthony
committed
// Create the required spectra mapping so that the workspace knows what to pad to
Gigg, Martyn Anthony
committed
createSpectraMapping(m_filename, WS, monitors, onebank);
Gigg, Martyn Anthony
committed
WS->padSpectra();
Janik Zikovsky
committed
//This map will be used to find the workspace index
Gigg, Martyn Anthony
committed
if( this->event_id_is_spec ) pixelID_to_wi_map = WS->getSpectrumToWorkspaceIndexMap();
else pixelID_to_wi_map = WS->getDetectorIDToWorkspaceIndexMap(true);
Janik Zikovsky
committed
// Cache a map for speed.
this->makeMapToEventLists();
// --------------------------- Time filtering ------------------------------------
Campbell, Stuart
committed
double filter_time_start_sec, filter_time_stop_sec;
filter_time_start_sec = getProperty("FilterByTimeStart");
filter_time_stop_sec = getProperty("FilterByTimeStop");
Campbell, Stuart
committed
//Default to ALL pulse times
bool is_time_filtered = false;