LoadQKK.cpp

/*WIKI* 


*WIKI*/
//---------------------------------------------------
// Includes
//---------------------------------------------------
#include "MantidDataHandling/LoadQKK.h"
#include "MantidAPI/FileProperty.h"
#include "MantidAPI/WorkspaceValidators.h"
#include "MantidKernel/UnitFactory.h"
#include "MantidAPI/LoadAlgorithmFactory.h"
#include "MantidGeometry/Instrument.h"
#include "MantidGeometry/Instrument/RectangularDetector.h"
#include "MantidGeometry/Objects/ShapeFactory.h"

#include "MantidNexus/NexusClasses.h"

#include <boost/math/special_functions/fpclassify.hpp>
#include <Poco/File.h>
#include <iostream>
#include <fstream>
#include <iomanip>

using namespace Mantid::DataHandling;
using namespace Mantid::API;
using namespace Mantid::Kernel;

namespace Mantid
{
  namespace DataHandling
  {

    // Register the algorithm into the AlgorithmFactory
    DECLARE_ALGORITHM(LoadQKK)

    //register the algorithm into loadalgorithm factory
    DECLARE_LOADALGORITHM(LoadQKK)

    /// Sets documentation strings for this algorithm
    void LoadQKK::initDocs()
    {
      this->setWikiSummary(
          "Loads a ANSTO QKK file. ");
      this->setOptionalMessage(
          "Loads a ANSTO QKK file. ");
    }

    /**
     * Initialise the algorithm. Declare properties which can be set before execution (input) or 
     * read from after the execution (output).
     */
    void LoadQKK::init()
    {
      // Specify file extensions which can be associated with a QKK file.
      std::vector<std::string> exts;
      exts.push_back(".nx.hdf");
      // Declare the Filename algorithm property. Mandatory. Sets the path to the file to load.
      declareProperty(new API::FileProperty("Filename", "", API::FileProperty::Load, exts),
          "The input filename of the stored data");
      // Declare the OutputWorkspace property. This sets the name of the workspace to be filled with the data
      // from the file. 
      declareProperty(new API::WorkspaceProperty<>("OutputWorkspace", "", Kernel::Direction::Output));
    }

    /**
     * Execute the algorithm.
     */
    void LoadQKK::exec()
    {
      using namespace Mantid::API;
      // Get the name of the file.
      std::string filename = getPropertyValue("Filename");

      // Open the root.
      NeXus::NXRoot root(filename);
      // Open the first NXentry found in the file.
      NeXus::NXEntry entry = root.openFirstEntry();
      // Open NXdata group with name "data"
      NeXus::NXData data = entry.openNXData("data");
      // Read in wavelength value
      double wavelength = static_cast<double>(data.getFloat("wavelength"));
      // open the data set with the counts. It is identified by the signal=1 attribute
      NeXus::NXInt hmm = data.openIntData();
      // Read the data into memory
      hmm.load();

      // Get the wavelength spread
      double wavelength_spread = static_cast<double>(entry.getFloat("instrument/velocity_selector/wavelength_spread"));
      double wavelength0 = wavelength - wavelength_spread / 2;
      double wavelength1 = wavelength + wavelength_spread / 2;

      // hmm is a 3d array with axes: sample_x, y_pixel_offset, x_pixel_offset
      size_t ny = hmm.dim1(); // second dimension
      size_t nx = hmm.dim2(); // third dimension
      size_t nHist = ny * nx; // number of spectra in the dataset
      if (nHist == 0)
      {
        throw std::runtime_error("Error in data dimensions: " + boost::lexical_cast<std::string>(ny) + " X " + boost::lexical_cast<std::string>(nx));
      }

      // Set the workspace structure. The workspace will contain nHist spectra each having a single wavelength bin.
      const size_t xWidth = 2; // number of wavelength bin boundaries
      const size_t yWidth = 1; // number of bins

      // Create a workspace with nHist spectra and a single y bin.
      MatrixWorkspace_sptr outputWorkspace = boost::dynamic_pointer_cast<MatrixWorkspace>(
          WorkspaceFactory::Instance().create("Workspace2D", nHist, xWidth, yWidth));
      // Set the units of the x axis as Wavelength
      outputWorkspace->getAxis(0)->unit() = UnitFactory::Instance().create("Wavelength");
      // Set the units of the data as Counts
      outputWorkspace->setYUnitLabel("Counts");

      //  Put the data into outputWorkspace
      size_t count = 0;
      for (size_t i = 0; i < ny; ++i)
      for (size_t j = 0; j < nx; ++j)
      {
        // Move data across
        double c = hmm(0,int(i),int(j));
        outputWorkspace->dataX(count)[0] = wavelength0;
        outputWorkspace->dataX(count)[1] = wavelength1;
        outputWorkspace->dataY(count)[0] = c;
        outputWorkspace->dataE(count)[0] = sqrt(c);
        ++count;
      }

      // Build instrument geometry

      // Create a new instrument and set its name
      std::string instrumentname = "QUOKKA";
      Geometry::Instrument_sptr instrument(new Geometry::Instrument(instrumentname));
      outputWorkspace->setInstrument(instrument);

      // Add dummy source and samplepos to instrument

      // Create an instrument component wich will represent the sample position.
      Geometry::ObjComponent *samplepos = new Geometry::ObjComponent("Sample",instrument.get());
      instrument->add(samplepos);
      instrument->markAsSamplePos(samplepos);
      // Put the sample in the centre of the coordinate system
      samplepos->setPos(0.0,0.0,0.0);

      // Create a component to represent the source
      Geometry::ObjComponent *source = new Geometry::ObjComponent("Source",instrument.get());
      instrument->add(source);
      instrument->markAsSource(source);

      // Read in the L1 value and place the source at (0,0,-L1)
      double l1 = static_cast<double>(entry.getFloat("instrument/parameters/L1"));
      source->setPos(0.0,0.0,-1.0*l1);

      // Create a component for the detector. 

      // We assumed that these are the dimensions of the detector, and height is in y direction and width is in x direction
      double height = static_cast<double>(entry.getFloat("instrument/detector/active_height"));
      double width = static_cast<double>(entry.getFloat("instrument/detector/active_width"));
      // Convert them to metres
      height /= 1000;
      width /= 1000;

      // We assumed that individual pixels have the same size and shape of a cuboid with dimensions:
      double pixel_height = height / static_cast<double>(ny);
      double pixel_width = width / static_cast<double>(nx);
      // Create size strings for shape creation
      std::string pixel_height_str = boost::lexical_cast<std::string>(pixel_height / 2);
      std::string pixel_width_str = boost::lexical_cast<std::string>(pixel_width / 2);
      // Set the depth of a pixel to a very small number
      std::string pixel_depth_str = "0.00001";

      // Create a RectangularDetector which represents a rectangular array of pixels
      Geometry::RectangularDetector* bank = new Geometry::RectangularDetector("bank",instrument.get());
      // Define shape of a pixel as an XML string. See http://www.mantidproject.org/HowToDefineGeometricShape for details
      // on shapes in Mantid.
      std::string detXML =   "<cuboid id=\"pixel\">"
        "<left-front-bottom-point   x= \""+pixel_width_str+"\" y=\"-"+pixel_height_str+"\" z=\"0\"  />"
        "<left-front-top-point      x= \""+pixel_width_str+"\" y=\"-"+pixel_height_str+"\" z=\""+pixel_depth_str+"\"  />"
        "<left-back-bottom-point    x=\"-"+pixel_width_str+"\" y=\"-"+pixel_height_str+"\" z=\"0\"  />"
        "<right-front-bottom-point  x= \""+pixel_width_str+"\" y= \""+pixel_height_str+"\" z=\"0\"  />"
        "</cuboid>";
      // Create a shape object which will be shared by all pixels.
      Geometry::Object_sptr shape = Geometry::ShapeFactory().createShape(detXML);
      // Initialise the detector specifying the sizes.
      bank->initialize(shape,int(nx),0,pixel_width,int(ny),0,pixel_height,1,true,int(nx));
      for (int i = 0; i < static_cast<int>(ny); ++i)
      for (int j = 0; j < static_cast<int>(nx); ++j)
      {
        instrument->markAsDetector(bank->getAtXY(j,i).get());
      }
      // Position the detector so the z axis goes through its centre
      bank->setPos(-width / 2, -height / 2, 0);

      // Set the workspace title
      outputWorkspace->setTitle(entry.getString("experiment/title"));
      // Attach the created workspace to the OutputWorkspace property. The workspace will also be saved in AnalysisDataService
      // and can be retrieved by its name.
      setProperty("OutputWorkspace", outputWorkspace);

    }

    /**This method does a quick file type check by checking 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 LoadQKK::quickFileCheck(const std::string& filePath, size_t nread, const file_header& header)
    {
      std::string extn=extension(filePath);
      bool bnexs(false);
      (!extn.compare("nxs")||!extn.compare("nx5")||!extn.compare("nx.hdf"))?bnexs=true:bnexs=false;
      /*
      * HDF files have magic cookie in the first 4 bytes
      */
      if ( ((nread >= sizeof(unsigned)) && (ntohl(header.four_bytes) == g_hdf_cookie)) || bnexs )
      {
        //hdf
        return true;
      }
      else if ( (nread >= sizeof(g_hdf5_signature)) && 
                (!memcmp(header.full_hdr, g_hdf5_signature, sizeof(g_hdf5_signature))) )
      { 
        //hdf5
        return true;
      }
      return false;
    }

    /** Checks the file by opening it and reading few lines
     *  @param filePath :: name of the file including its path
     *  @return an integer value how much this algorithm can load the file
     */
    int LoadQKK::fileCheck(const std::string& filePath)
    {
      NeXus::NXRoot root(filePath);

      try
      {
        // Check if there exists a data set with name hmm_xy
        NeXus::NXEntry entry = root.openFirstEntry();
        if (entry.containsGroup("data"))
        {
          NeXus::NXData data = entry.openNXData("data");
          if (data.getDataSetInfo("hmm_xy").stat != NX_ERROR)
          {
            return 80; // Give it an 80% chance that the file is readable
          }
        }
      }
      catch(...)
      {
        return 0;
      }
      return 0;
    }

  }//namespace
}//namespace