SaveNexusProcessed.cpp

// SaveNexusProcessed
// @author Ronald Fowler, based on SaveNexus
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidDataHandling/SaveNexusProcessed.h"
#include "MantidNexus/NexusFileIO.h"
#include "MantidAPI/WorkspaceOpOverloads.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidKernel/ConfigService.h"
#include "MantidAPI/FileProperty.h"

#include <Poco/File.h>
#include <Poco/Path.h>

#include <cmath>
#include <boost/shared_ptr.hpp>
#include "MantidNexusCPP/NeXusFile.hpp"

using namespace Mantid::API;

namespace Mantid
{
namespace DataHandling
{

  using namespace Kernel;
  using namespace API;
  using namespace DataObjects;
  using Geometry::Instrument_const_sptr;

  // Register the algorithm into the algorithm factory
  DECLARE_ALGORITHM(SaveNexusProcessed)
  
  /// Sets documentation strings for this algorithm
  void SaveNexusProcessed::initDocs()
  {
    this->setWikiSummary("The SaveNexusProcessed algorithm will write the given Mantid workspace to a Nexus file. SaveNexusProcessed may be invoked by [[SaveNexus]]. ");
    this->setOptionalMessage("The SaveNexusProcessed algorithm will write the given Mantid workspace to a Nexus file. SaveNexusProcessed may be invoked by SaveNexus.");
  }
  


  /// Empty default constructor
  SaveNexusProcessed::SaveNexusProcessed() : Algorithm()
  {
  }

  //-----------------------------------------------------------------------------------------------
  /** Initialisation method.
   *
   */
  void SaveNexusProcessed::init()
  {
    declareProperty(new WorkspaceProperty<Workspace>("InputWorkspace","",Direction::Input),
        "Name of the workspace to be saved");
    // Declare required input parameters for algorithm
    std::vector<std::string> exts;
    exts.push_back(".nxs");
    exts.push_back(".nx5");
    exts.push_back(".xml");

    declareProperty(new FileProperty("Filename", "", FileProperty::Save, exts),
        "The name of the Nexus file to write, as a full or relative\n"
        "path");

    // Declare optional parameters (title now optional, was mandatory)
    declareProperty("Title", "", new NullValidator<std::string>,
        "A title to describe the saved workspace");
    BoundedValidator<int> *mustBePositive = new BoundedValidator<int>();
    mustBePositive->setLower(0);

    declareProperty("WorkspaceIndexMin", 0, mustBePositive->clone(), 
        "Index number of first spectrum to write, only for single\n"
        "period data. Not yet implemented");
    declareProperty("WorkspaceIndexMax", Mantid::EMPTY_INT(), mustBePositive->clone(),
        "Index of last spectrum to write, only for single period\n"
        "data. Not yet implemented");
    declareProperty(new ArrayProperty<int>("WorkspaceIndexList"),
        "List of spectrum numbers to read, only for single period\n"
        "data. Not yet implemented");
    //declareProperty("EntryNumber", 0, mustBePositive);
    declareProperty("Append",false,"Determines whether .nxs file needs to be\n"
        "over written or appended");
  }


  /** Get the list of workspace indices to use
   *
   * @param spec :: returns the list of workspace indices
   */
  void SaveNexusProcessed::getSpectrumList(std::vector<int> & spec, MatrixWorkspace_const_sptr matrixWorkspace)
  {
    std::vector<int> spec_list = getProperty("WorkspaceIndexList");
    int spec_min = getProperty("WorkspaceIndexMin");
    int spec_max = getProperty("WorkspaceIndexMax");
    const bool list = !spec_list.empty();
    const bool interval = (spec_max != Mantid::EMPTY_INT());
    if ( spec_max == Mantid::EMPTY_INT() ) spec_max = 0;
    const int numberOfHist = static_cast<int>(matrixWorkspace->getNumberHistograms());

    if( interval )
    {
      if ( spec_max < spec_min || spec_max > numberOfHist-1 )
      {
        g_log.error("Invalid WorkspaceIndex min/max properties");
        throw std::invalid_argument("Inconsistent properties defined");
      }
      spec.reserve(1+spec_max-spec_min);
      for(int i=spec_min;i<=spec_max;i++)
        spec.push_back(i);
      if (list)
      {
        for(size_t i=0;i<spec_list.size();i++)
        {
          int s = spec_list[i];
          if ( s < 0 ) continue;
          if (s < spec_min || s > spec_max)
            spec.push_back(s);
        }
      }
    }
    else if (list)
    {
      spec_max=0;
      spec_min=numberOfHist-1;
      for(size_t i=0;i<spec_list.size();i++)
      {
        int s = spec_list[i];
        if ( s < 0 ) continue;
        spec.push_back(s);
        if (s > spec_max) spec_max = s;
        if (s < spec_min) spec_min = s;
      }
    }
    else
    {
      spec_min=0;
      spec_max=numberOfHist-1;
      spec.reserve(1+spec_max-spec_min);
      for(int i=spec_min;i<=spec_max;i++)
        spec.push_back(i);
    }
  }



  //-----------------------------------------------------------------------------------------------
  /** Executes the algorithm.
   *
   *  @throw runtime_error Thrown if algorithm cannot execute
   */
  void SaveNexusProcessed::exec()
  {
    //TODO: Remove?
    NXMEnableErrorReporting();

    Workspace_sptr inputWorkspace = getProperty("InputWorkspace");
    boost::shared_ptr<WorkspaceGroup> wsGrpSptr =
                  boost::dynamic_pointer_cast<WorkspaceGroup>(inputWorkspace);
    if(wsGrpSptr)
    {
      processGroups(wsGrpSptr, this->getProperties());
      return;
    }

    Progress prog_init(this, 0.0, 0.3, 5);

    // Retrieve the filename from the properties
    m_filename = getPropertyValue("Filename");
    //m_entryname = getPropertyValue("EntryName");
    m_title = getPropertyValue("Title");
    MatrixWorkspace_const_sptr matrixWorkspace = boost::dynamic_pointer_cast<const MatrixWorkspace>(inputWorkspace);
    ITableWorkspace_const_sptr tableWorkspace = boost::dynamic_pointer_cast<const ITableWorkspace>(inputWorkspace);
    // check if inputWorkspace is something we know how to save
    if (!matrixWorkspace && !tableWorkspace) 
      return;
    m_eventWorkspace = boost::dynamic_pointer_cast<const EventWorkspace>(matrixWorkspace);

    // If no title's been given, use the workspace title field
    if (m_title.empty()) 
      m_title = inputWorkspace->getTitle();

    // If we don't want to append then remove the file if it already exists
    bool append_to_file = getProperty("Append");
    if( !append_to_file )
    {
      Poco::File file(m_filename);
      if( file.exists() )
        file.remove();
    }


    const std::string workspaceID = inputWorkspace->id();
    if ((workspaceID.find("Workspace2D") == std::string::npos) &&
        !m_eventWorkspace && !tableWorkspace)
      throw Exception::NotImplementedError("SaveNexusProcessed passed invalid workspaces. Must be Workspace2D, EventWorkspace or ITableWorkspace.");


    Mantid::NeXus::NexusFileIO *nexusFile= new Mantid::NeXus::NexusFileIO();

    if( nexusFile->openNexusWrite( m_filename ) != 0 )
      throw Exception::FileError("Failed to open file", m_filename);

    // Equivalent C++ API handle
    ::NeXus::File * cppFile = new ::NeXus::File(nexusFile->fileID);

    prog_init.reportIncrement(1, "Opening file");
    if( nexusFile->writeNexusProcessedHeader( m_title ) != 0 )
      throw Exception::FileError("Failed to write to file", m_filename);

    prog_init.reportIncrement(1, "Writing header");

    // write instrument data, if present and writer enabled
    if (matrixWorkspace) 
    { 
      // Save the instrument names, ParameterMap, sample, run
      matrixWorkspace->saveExperimentInfoNexus(cppFile);
      prog_init.reportIncrement(1, "Writing sample and instrument");

      // check if all X() are in fact the same array
      const bool uniformSpectra = API::WorkspaceHelpers::commonBoundaries(matrixWorkspace);

      // Retrieve the workspace indices (from params)
      std::vector<int> spec;
      this->getSpectrumList(spec, matrixWorkspace);

      // Write out the data (2D or event)
      if (m_eventWorkspace)
      {
        this->execEvent(nexusFile,uniformSpectra,spec);
      }
      else
      {
        nexusFile->writeNexusProcessedData2D(matrixWorkspace,uniformSpectra,spec, "workspace", true);
      }

      // MW 27/10/10 - don't try and save the spectra-detector map if there isn't one
      if ( matrixWorkspace->getAxis(1)->isSpectra() )
      {
        cppFile->openGroup("instrument", "NXinstrument");
        matrixWorkspace->saveSpectraMapNexus(cppFile, spec, ::NeXus::LZW);
        cppFile->closeGroup();
      }

    }  // finish matrix workspace specifics 
    
    if (tableWorkspace)
    {
      nexusFile->writeNexusTableWorkspace(tableWorkspace,"table_workspace");
    }
 
    nexusFile->writeNexusProcessedProcess(inputWorkspace);
    nexusFile->closeNexusFile();

    delete nexusFile;

    return;
  }






  //-------------------------------------------------------------------------------------
  /** Append out each field of a vector of events to separate array.
   *
   * @param events :: vector of TofEvent or WeightedEvent, etc.
   * @param offset :: where the first event goes in the array
   * @param tofs, weights, errorSquareds, pulsetimes :: arrays to write to.
   *        Must be initialized and big enough,
   *        or NULL if they are not meant to be written to.
   */
  template<class T>
  void SaveNexusProcessed::appendEventListData( std::vector<T> events, size_t offset, double * tofs, float * weights, float * errorSquareds, int64_t * pulsetimes)
  {
    // Do nothing if there are no events.
    size_t num = events.size();
    if (num <= 0)
      return;

    typename std::vector<T>::const_iterator it;
    typename std::vector<T>::const_iterator it_end = events.end();
    size_t i = offset;

    // Fill the C-arrays with the fields from all the events, as requested.
    for (it = events.begin(); it != it_end; it++)
    {
      if (tofs) tofs[i] = it->tof();
      // FIXME: ticket #3051
      if (weights) weights[i] = static_cast<float>(it->weight());
      if (errorSquareds) errorSquareds[i] = static_cast<float>(it->errorSquared());
      if (pulsetimes) pulsetimes[i] = it->pulseTime().total_nanoseconds();
      i++;
    }
  }



  //-----------------------------------------------------------------------------------------------
  /** Execute the saving of event data.
   * This will make one long event list for all events contained.
   * */
  void SaveNexusProcessed::execEvent(Mantid::NeXus::NexusFileIO * nexusFile,const bool uniformSpectra,const std::vector<int> spec)
  {
    prog = new Progress(this, 0.3, 1.0, m_eventWorkspace->getNumberEvents()*2);

    // Start by writing out the axes and crap
    nexusFile->writeNexusProcessedData2D(m_eventWorkspace, uniformSpectra, spec, "event_workspace", false);

    // Make a super long list of tofs, weights, etc.
    std::vector<int64_t> indices;
    indices.reserve( m_eventWorkspace->getNumberHistograms()+1 );
    // First we need to index the events in each spectrum
    size_t index = 0;
    for (int wi =0; wi < static_cast<int>(m_eventWorkspace->getNumberHistograms()); wi++)
    {
      indices.push_back(index);
      // Track the total # of events
      index += m_eventWorkspace->getEventList(wi).getNumberEvents();
    }
    indices.push_back(index);

    // Initialize all the arrays
    int64_t num = index;
    double * tofs = NULL;
    float * weights = NULL;
    float * errorSquareds = NULL;
    int64_t * pulsetimes = NULL;

    // overall event type.
    EventType type = m_eventWorkspace->getEventType();
    bool writeTOF = true;
    bool writePulsetime = false;
    bool writeWeight = false;
    bool writeError = false;

    switch (type)
    {
    case TOF:
      writePulsetime = true;
      break;
    case WEIGHTED:
      writePulsetime = true;
      writeWeight = true;
      writeError = true;
      break;
    case WEIGHTED_NOTIME:
      writeWeight = true;
      writeError = true;
      break;
    }

    // --- Initialize the combined event arrays ----
    if (writeTOF)
      tofs = new double[num];
    if (writeWeight)
      weights = new float[num];
    if (writeError)
      errorSquareds = new float[num];
    if (writePulsetime)
      pulsetimes = new int64_t[num];

    // --- Fill in the combined event arrays ----
    PARALLEL_FOR_NO_WSP_CHECK()
    for (int wi=0; wi < static_cast<int>(m_eventWorkspace->getNumberHistograms()); wi++)
    {
      PARALLEL_START_INTERUPT_REGION
      const DataObjects::EventList & el = m_eventWorkspace->getEventList(wi);

      // This is where it will land in the output array.
      // It is okay to write in parallel since none should step on each other.
      size_t offset = indices[wi];
      std::string eventType;

      switch (el.getEventType())
      {
      case TOF:
        eventType = "TOF";
        appendEventListData( el.getEvents(), offset, tofs, weights, errorSquareds, pulsetimes);
        break;
      case WEIGHTED:
        eventType = "WEIGHTED";
        appendEventListData( el.getWeightedEvents(), offset, tofs, weights, errorSquareds, pulsetimes);
        break;
      case WEIGHTED_NOTIME:
        eventType = "WEIGHTED_NOTIME";
        appendEventListData( el.getWeightedEventsNoTime(), offset, tofs, weights, errorSquareds, pulsetimes);
        break;
      }
      prog->reportIncrement(el.getNumberEvents(), "Copying EventList");

      PARALLEL_END_INTERUPT_REGION
    }
    PARALLEL_CHECK_INTERUPT_REGION


    // Write out to the NXS file.
    nexusFile->writeNexusProcessedDataEventCombined(m_eventWorkspace, indices, tofs, weights, errorSquareds, pulsetimes,
        false /*DONT compress - much faster*/);

    // Free mem.
    delete [] tofs;
    delete [] weights;
    delete [] errorSquareds;
    delete [] pulsetimes;
  }

  //-----------------------------------------------------------------------------------------------
  /** virtual method to set the non workspace properties for this algorithm
   *  @param alg :: pointer to the algorithm
   *  @param propertyName :: name of the property
   *  @param propertyValue :: value  of the property
   *  @param perioidNum :: period number
   */
  void SaveNexusProcessed::setOtherProperties(IAlgorithm* alg,const std::string& propertyName,const std::string& propertyValue,int perioidNum)
  {
    if(!propertyName.compare("Append"))
    {	if(perioidNum!=1)
    { alg->setPropertyValue(propertyName,"1");
    }
    else alg->setPropertyValue(propertyName,propertyValue);
    }
    else
      Algorithm::setOtherProperties(alg,propertyName,propertyValue,perioidNum);
  }

} // namespace DataHandling
} // namespace Mantid