//---------------------------------------------------------------------- // Includes //---------------------------------------------------------------------- #include "MantidDataHandling/LoadMuonNexus2.h" #include "MantidDataHandling/LoadMuonNexus1.h" #include "MantidAPI/Axis.h" #include "MantidAPI/FileProperty.h" #include "MantidAPI/Progress.h" #include "MantidAPI/RegisterFileLoader.h" #include "MantidAPI/WorkspaceFactory.h" #include "MantidGeometry/Instrument/Detector.h" #include "MantidKernel/TimeSeriesProperty.h" #include "MantidKernel/UnitFactory.h" #include "MantidKernel/ConfigService.h" #include "MantidKernel/ArrayProperty.h" #include "MantidKernel/UnitLabelTypes.h" #include "MantidNexus/NexusClasses.h" #include "MantidDataObjects/Workspace2D.h" #include <nexus/NeXusFile.hpp> #include <nexus/NeXusException.hpp> #include <Poco/Path.h> #include <boost/lexical_cast.hpp> #include <boost/shared_ptr.hpp> #include <cmath> #include <numeric> namespace Mantid { namespace DataHandling { // Register the algorithm into the algorithm factory DECLARE_NEXUS_FILELOADER_ALGORITHM(LoadMuonNexus2) using namespace Kernel; using namespace API; using Geometry::Instrument; using namespace Mantid::NeXus; /// Empty default constructor LoadMuonNexus2::LoadMuonNexus2() : LoadMuonNexus() {} /** Executes the right version of the muon nexus loader: versions 1 or 2. * * @throw Exception::FileError If the Nexus file cannot be found/opened * @throw std::invalid_argument If the optional properties are set to invalid *values */ void LoadMuonNexus2::exec() { std::string filePath = getPropertyValue("Filename"); LoadMuonNexus1 load1; load1.initialize(); Kernel::NexusDescriptor descriptor(filePath); int confidence1 = load1.confidence(descriptor); int confidence2 = this->confidence(descriptor); // if none can load the file throw if (confidence1 < 80 && confidence2 < 80) { throw Kernel::Exception::FileError("Cannot open the file ", filePath); } if (confidence2 > confidence1) { // this loader doExec(); } else { // version 1 loader IAlgorithm_sptr childAlg = createChildAlgorithm("LoadMuonNexus", 0, 1, true, 1); auto version1Loader = boost::dynamic_pointer_cast<API::Algorithm>(childAlg); version1Loader->copyPropertiesFrom(*this); version1Loader->executeAsChildAlg(); this->copyPropertiesFrom(*version1Loader); API::Workspace_sptr outWS = version1Loader->getProperty("OutputWorkspace"); setProperty("OutputWorkspace", outWS); } } /** Read in a muon nexus file of the version 2. * * @throw Exception::FileError If the Nexus file cannot be found/opened * @throw std::invalid_argument If the optional properties are set to invalid *values */ void LoadMuonNexus2::doExec() { // Create the root Nexus class NXRoot root(getPropertyValue("Filename")); int64_t iEntry = getProperty("EntryNumber"); if (iEntry >= static_cast<int64_t>(root.groups().size())) { throw std::invalid_argument("EntryNumber is out of range"); } // Open the data entry m_entry_name = root.groups()[iEntry].nxname; NXEntry entry = root.openEntry(m_entry_name); // Read in the instrument name from the Nexus file m_instrument_name = entry.getString("instrument/name"); // Read the number of periods in this file if (entry.containsGroup("run")) { try { m_numberOfPeriods = entry.getInt("run/number_periods"); } catch (::NeXus::Exception &) { // assume 1 m_numberOfPeriods = 1; } } else { m_numberOfPeriods = 1; } // Need to extract the user-defined output workspace name Property *ws = getProperty("OutputWorkspace"); std::string localWSName = ws->value(); // If multiperiod, will need to hold the Instrument & Sample for copying boost::shared_ptr<Instrument> instrument; boost::shared_ptr<Sample> sample; std::string detectorName; // Only the first NXdata found for (auto &group : entry.groups()) { std::string className = group.nxclass; if (className == "NXdata") { detectorName = group.nxname; break; } } NXData dataGroup = entry.openNXData(detectorName); Mantid::NeXus::NXInt spectrum_index = dataGroup.openNXInt("spectrum_index"); spectrum_index.load(); m_numberOfSpectra = spectrum_index.dim0(); // Call private method to validate the optional parameters, if set checkOptionalProperties(); NXFloat raw_time = dataGroup.openNXFloat("raw_time"); raw_time.load(); int nBins = raw_time.dim0(); std::vector<double> timeBins; timeBins.assign(raw_time(), raw_time() + nBins); timeBins.push_back(raw_time[nBins - 1] + raw_time[1] - raw_time[0]); // Calculate the size of a workspace, given its number of periods & spectra to // read int total_specs; if (m_interval || m_list) { total_specs = static_cast<int>(m_spec_list.size()); if (m_interval) { total_specs += static_cast<int>((m_spec_max - m_spec_min + 1)); } else { m_spec_max = -1; // to stop entering the min - max loop } } else { total_specs = static_cast<int>(m_numberOfSpectra); // for nexus return all spectra m_spec_min = 1; m_spec_max = m_numberOfSpectra; // was +1? } // Create the 2D workspace for the output DataObjects::Workspace2D_sptr localWorkspace = boost::dynamic_pointer_cast<DataObjects::Workspace2D>( WorkspaceFactory::Instance().create("Workspace2D", total_specs, nBins + 1, nBins)); // Set the unit on the workspace to muon time, for now in the form of a Label // Unit boost::shared_ptr<Kernel::Units::Label> lblUnit = boost::dynamic_pointer_cast<Kernel::Units::Label>( UnitFactory::Instance().create("Label")); lblUnit->setLabel("Time", Units::Symbol::Microsecond); localWorkspace->getAxis(0)->unit() = lblUnit; // Set y axis unit localWorkspace->setYUnit("Counts"); // g_log.error()<<" number of perioids= "<<m_numberOfPeriods<<std::endl; WorkspaceGroup_sptr wsGrpSptr = WorkspaceGroup_sptr(new WorkspaceGroup); if (entry.containsDataSet("title")) { wsGrpSptr->setTitle(entry.getString("title")); } if (entry.containsDataSet("notes")) { wsGrpSptr->setComment(entry.getString("notes")); } if (m_numberOfPeriods > 1) { setProperty("OutputWorkspace", boost::dynamic_pointer_cast<Workspace>(wsGrpSptr)); } // period_index is currently unused // Mantid::NeXus::NXInt period_index = dataGroup.openNXInt("period_index"); // period_index.load(); Mantid::NeXus::NXInt counts = dataGroup.openIntData(); counts.load(); NXInstrument instr = entry.openNXInstrument("instrument"); if (instr.containsGroup("detector_fb")) { NXDetector detector = instr.openNXDetector("detector_fb"); if (detector.containsDataSet("time_zero")) { double dum = detector.getFloat("time_zero"); setProperty("TimeZero", dum); } if (detector.containsDataSet("first_good_time")) { double dum = detector.getFloat("first_good_time"); setProperty("FirstGoodData", dum); } } API::Progress progress(this, 0., 1., m_numberOfPeriods * total_specs); // Loop over the number of periods in the Nexus file, putting each period in a // separate workspace for (int period = 0; period < m_numberOfPeriods; ++period) { if (period == 0) { // Only run the Child Algorithms once loadRunDetails(localWorkspace); runLoadInstrument(localWorkspace); loadLogs(localWorkspace, entry, period); } else // We are working on a higher period of a multiperiod raw file { localWorkspace = boost::dynamic_pointer_cast<DataObjects::Workspace2D>( WorkspaceFactory::Instance().create(localWorkspace)); } std::string outws(""); if (m_numberOfPeriods > 1) { std::string outputWorkspace = "OutputWorkspace"; std::stringstream suffix; suffix << (period + 1); outws = outputWorkspace + "_" + suffix.str(); std::string WSName = localWSName + "_" + suffix.str(); declareProperty( new WorkspaceProperty<Workspace>(outws, WSName, Direction::Output)); if (wsGrpSptr) wsGrpSptr->addWorkspace(localWorkspace); } // create spectrum -> index correspondence std::map<int, int> index_spectrum; for (int i = 0; i < m_numberOfSpectra; ++i) { index_spectrum[spectrum_index[i]] = i; } int counter = 0; for (int spec = static_cast<int>(m_spec_min); spec <= static_cast<int>(m_spec_max); ++spec) { int i = index_spectrum[spec]; // if spec not found i is 0 loadData(counts, timeBins, counter, period, i, localWorkspace); localWorkspace->getSpectrum(counter)->setSpectrumNo(spectrum_index[i]); counter++; progress.report(); } // Read in the spectra in the optional list parameter, if set if (m_list) { for (auto spec : m_spec_list) { int k = index_spectrum[spec]; // if spec not found k is 0 loadData(counts, timeBins, counter, period, k, localWorkspace); localWorkspace->getSpectrum(counter)->setSpectrumNo(spectrum_index[k]); counter++; progress.report(); } } // Just a sanity check assert(counter == total_specs); bool autogroup = getProperty("AutoGroup"); if (autogroup) { g_log.warning( "Autogrouping is not implemented for muon NeXus version 2 files"); } // Assign the result to the output workspace property if (m_numberOfPeriods > 1) setProperty(outws, boost::static_pointer_cast<Workspace>(localWorkspace)); else { setProperty("OutputWorkspace", boost::dynamic_pointer_cast<Workspace>(localWorkspace)); } } // loop over periods } /** loadData * Load the counts data from an NXInt into a workspace */ void LoadMuonNexus2::loadData(const Mantid::NeXus::NXInt &counts, const std::vector<double> &timeBins, int wsIndex, int period, int spec, API::MatrixWorkspace_sptr localWorkspace) { MantidVec &X = localWorkspace->dataX(wsIndex); MantidVec &Y = localWorkspace->dataY(wsIndex); MantidVec &E = localWorkspace->dataE(wsIndex); X.assign(timeBins.begin(), timeBins.end()); int nBins = 0; int *data = nullptr; if (counts.rank() == 3) { nBins = counts.dim2(); data = &counts(period, spec, 0); } else if (counts.rank() == 2) { nBins = counts.dim1(); data = &counts(spec, 0); } else { throw std::runtime_error("Data have unsupported dimansionality"); } assert(nBins + 1 == static_cast<int>(timeBins.size())); Y.assign(data, data + nBins); typedef double (*uf)(double); uf dblSqrt = std::sqrt; std::transform(Y.begin(), Y.end(), E.begin(), dblSqrt); } /** Load logs from Nexus file. Logs are expected to be in * /run/sample group of the file. * @param ws :: The workspace to load the logs to. * @param entry :: The Nexus entry * @param period :: The period of this workspace */ void LoadMuonNexus2::loadLogs(API::MatrixWorkspace_sptr ws, NXEntry &entry, int period) { // Avoid compiler warning (void)period; std::string start_time = entry.getString("start_time"); std::string sampleName = entry.getString("sample/name"); NXMainClass runlogs = entry.openNXClass<NXMainClass>("sample"); ws->mutableSample().setName(sampleName); for (std::vector<NXClassInfo>::const_iterator it = runlogs.groups().begin(); it != runlogs.groups().end(); ++it) { NXLog nxLog = runlogs.openNXLog(it->nxname); Kernel::Property *logv = nxLog.createTimeSeries(start_time); if (!logv) continue; ws->mutableRun().addLogData(logv); } ws->setTitle(entry.getString("title")); if (entry.containsDataSet("notes")) { ws->setComment(entry.getString("notes")); } std::string run_num = boost::lexical_cast<std::string>(entry.getInt("run_number")); // The sample is left to delete the property ws->mutableRun().addLogData( new PropertyWithValue<std::string>("run_number", run_num)); ws->populateInstrumentParameters(); } /** Log the run details from the file * @param localWorkspace :: The workspace details to use */ void LoadMuonNexus2::loadRunDetails( DataObjects::Workspace2D_sptr localWorkspace) { API::Run &runDetails = localWorkspace->mutableRun(); runDetails.addProperty("run_title", localWorkspace->getTitle(), true); int numSpectra = static_cast<int>(localWorkspace->getNumberHistograms()); runDetails.addProperty("nspectra", numSpectra); m_filename = getPropertyValue("Filename"); NXRoot root(m_filename); NXEntry entry = root.openEntry(m_entry_name); std::string start_time = entry.getString("start_time"); runDetails.addProperty("run_start", start_time); std::string stop_time = entry.getString("end_time"); runDetails.addProperty("run_end", stop_time); if (entry.containsGroup("run")) { NXClass runRun = entry.openNXGroup("run"); if (runRun.containsDataSet("good_total_frames")) { int dum = runRun.getInt("good_total_frames"); runDetails.addProperty("goodfrm", dum); } if (runRun.containsDataSet("number_periods")) { int dum = runRun.getInt("number_periods"); runDetails.addProperty("nperiods", dum); } } { // Duration taken to be stop_time minus stat_time DateAndTime start(start_time); DateAndTime end(stop_time); double duration_in_secs = DateAndTime::secondsFromDuration(end - start); runDetails.addProperty("dur_secs", duration_in_secs); } } /** * Return the confidence with with this algorithm can load the file * @param descriptor A descriptor for the file * @returns An integer specifying the confidence level. 0 indicates it will not * be used */ int LoadMuonNexus2::confidence(Kernel::NexusDescriptor &descriptor) const { const auto &firstEntryNameType = descriptor.firstEntryNameType(); const std::string root = "/" + firstEntryNameType.first; if (!descriptor.pathExists(root + "/definition")) return 0; bool upperIDF(true); if (descriptor.pathExists(root + "/IDF_version")) upperIDF = true; else { if (descriptor.pathExists(root + "/idf_version")) upperIDF = false; else return 0; } try { std::string versionField = "idf_version"; if (upperIDF) versionField = "IDF_version"; auto &file = descriptor.data(); file.openPath(root + "/" + versionField); int32_t version = 0; file.getData(&version); if (version != 2) return 0; file.openPath(root + "/definition"); std::string def = file.getStrData(); if (def == "muonTD" || def == "pulsedTD") { // If all this succeeded then we'll assume this is an ISIS Muon NeXus file // version 2 return 81; } } catch (...) { } return 0; } } // namespace DataHandling } // namespace Mantid