ManagedRawFileWorkspace2D.cpp

#include "MantidDataHandling/ManagedRawFileWorkspace2D.h"
#include "MantidKernel/Exception.h"
#include "MantidAPI/RefAxis.h"
#include "MantidAPI/LocatedDataRef.h"
#include "MantidAPI/WorkspaceIterator.h"
#include "MantidAPI/WorkspaceIteratorCode.h"
#include "MantidAPI/WorkspaceProperty.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidKernel/UnitFactory.h"
#include "MantidKernel/ConfigService.h"
#include "MantidKernel/System.h"

#include "LoadRaw/isisraw2.h"
#include <cstdio>
#include <boost/timer.hpp>
#include <Poco/File.h>
#include <Poco/Path.h>
#include <Poco/Exception.h>
#include "MantidDataObjects/ManagedHistogram1D.h"

using Mantid::DataObjects::ManagedHistogram1D;

//DECLARE_WORKSPACE(ManagedRawFileWorkspace2D)

namespace Mantid
{
  namespace DataHandling
  {

    // Get a reference to the logger
    Kernel::Logger& ManagedRawFileWorkspace2D::g_log = Kernel::Logger::get("ManagedRawFileWorkspace2D");

    // Initialise the instance count
    int64_t ManagedRawFileWorkspace2D::g_uniqueID = 1;

    /// Constructor
    ManagedRawFileWorkspace2D::ManagedRawFileWorkspace2D(const std::string& fileName, int opt) :
      ManagedWorkspace2D(),
      isisRaw(new ISISRAW2),m_filenameRaw(fileName),m_fileRaw(NULL),m_readIndex(0),m_nmonitorSkipCounter(0)
    {
      this->setRawFile(opt);
    }

    ///Destructor
    ManagedRawFileWorkspace2D::~ManagedRawFileWorkspace2D()
    {
      if (m_fileRaw) fclose(m_fileRaw);
      removeTempFile();
    }

    /** Sets the RAW file for this workspace.
        @param opt :: Caching option.  0 - cache on local drive if raw file is very slow to read.
        1 - cache anyway, 2 - never cache.
    */
    void ManagedRawFileWorkspace2D::setRawFile(const int opt)
    {
      m_fileRaw = fopen(m_filenameRaw.c_str(),"rb");
      if (m_fileRaw == NULL)
      {
        g_log.error("Unable to open file " + m_filenameRaw);
        throw Kernel::Exception::FileError("Unable to open File:" , m_filenameRaw);
      }

      if ( needCache(opt) )
      {
        openTempFile();
      }

      isisRaw->ioRAW(m_fileRaw, true);

      m_numberOfBinBoundaries = isisRaw->t_ntc1 + 1;    
      m_numberOfPeriods = isisRaw->t_nper;
      initialize(isisRaw->t_nsp1,m_numberOfBinBoundaries,isisRaw->t_ntc1);
      int64_t noOfBlocks = m_noVectors / m_vectorsPerBlock;
      if ( noOfBlocks * m_vectorsPerBlock != m_noVectors ) ++noOfBlocks;
      m_changedBlock.resize(noOfBlocks,false);

      isisRaw->skipData(m_fileRaw,0);
      fgetpos(m_fileRaw, &m_data_pos); //< Save the data start position.

      getTimeChannels();
      getAxis(0)->unit() = Kernel::UnitFactory::Instance().create("TOF");
      setYUnit("Counts");
      setTitle(std::string(isisRaw->r_title, 80));
    }

    /// Constructs the time channel (X) vector(s)
    void ManagedRawFileWorkspace2D::getTimeChannels()
    {
      float* const timeChannels = new float[m_numberOfBinBoundaries];
      isisRaw->getTimeChannels(timeChannels, static_cast<int>(m_numberOfBinBoundaries));

      const int regimes = isisRaw->daep.n_tr_shift;
      if ( regimes >=2 )
      {
        g_log.debug() << "Raw file contains " << regimes << " time regimes\n"; 
        // If more than 1 regime, create a timeChannelsVec for each regime
        for (int i=0; i < regimes; ++i)
        {
          // Create a vector with the 'base' time channels
          boost::shared_ptr<MantidVec> channelsVec(
            new MantidVec(timeChannels,timeChannels + m_numberOfBinBoundaries));
          const double shift = isisRaw->daep.tr_shift[i];
          g_log.debug() << "Time regime " << i+1 << " shifted by " << shift << " microseconds\n";
          // Add on the shift for this vector
          std::transform(channelsVec->begin(), channelsVec->end(),
                         channelsVec->begin(), std::bind2nd(std::plus<double>(),shift));
          m_timeChannels.push_back(channelsVec);
        }
        // In this case, also need to populate the map of spectrum-regime correspondence
        const int64_t ndet = static_cast<int64_t>(isisRaw->i_det);
        std::map<int64_t,int64_t>::iterator hint = m_specTimeRegimes.begin();
        for (int64_t j=0; j < ndet; ++j)
        {
          // No checking for consistency here - that all detectors for given spectrum
          // are declared to use same time regime. Will just use first encountered
          hint = m_specTimeRegimes.insert(hint,std::make_pair(isisRaw->spec[j],isisRaw->timr[j]));
        }
      }
      else // Just need one in this case
      {
        boost::shared_ptr<MantidVec> channelsVec(
          new MantidVec(timeChannels,timeChannels + m_numberOfBinBoundaries));
        m_timeChannels.push_back(channelsVec);
      }
      // Done with the timeChannels C array so clean up
      delete[] timeChannels;
    }

    // Pointer to sqrt function (used in calculating errors below)
    typedef double (*uf)(double);
    static uf dblSqrt = std::sqrt;

    // readData(int) should be changed to readNextSpectrum() returning the spectrum index
    // and skipData to skipNextSpectrum()
    void ManagedRawFileWorkspace2D::readDataBlock(DataObjects::ManagedDataBlock2D *newBlock,size_t startIndex)const
    {
      Poco::ScopedLock<Poco::FastMutex> mutex(m_mutex);
      if (!m_fileRaw)
      {
        g_log.error("Raw file was not open.");
        throw std::runtime_error("Raw file was not open.");
      }
      int64_t blockIndex = startIndex / m_vectorsPerBlock;

      // Modified data is stored in ManagedWorkspace2D flat file.
      if (m_changedBlock[blockIndex])
      {
        ManagedWorkspace2D::readDataBlock(newBlock,startIndex);
        return;
      }

      if(m_monitorList.size()>0)
      { 
        if ( static_cast<int>(startIndex) > m_readIndex)
        {
          while(static_cast<int>(startIndex) >  m_readIndex-m_nmonitorSkipCounter)
          {
            isisRaw->skipData(m_fileRaw,m_readIndex+1);// Adding 1 because we dropped the first spectrum.
            ++m_readIndex;
            if(isMonitor(m_readIndex))
              ++m_nmonitorSkipCounter;
          }
        }
        else
        {
          int nwords = 0;
          while(static_cast<int>(startIndex)+ m_nmonitorSkipCounter+1 <  m_readIndex)
          {
            if(isMonitor(m_readIndex))
              --m_nmonitorSkipCounter;
            --m_readIndex;
            nwords += 4*isisRaw->ddes[m_readIndex+1].nwords;
          }
          if (fseek(m_fileRaw,-nwords,SEEK_CUR) != 0)
          {
            fclose(m_fileRaw);
            removeTempFile();
            g_log.error("Error reading RAW file.");
            throw std::runtime_error("ManagedRawFileWorkspace2D: Error reading RAW file.");
          }
        }
        int64_t endIndex = startIndex+m_vectorsPerBlock < m_noVectors?startIndex+m_vectorsPerBlock:m_noVectors;
        if (endIndex >= static_cast<int64_t>(m_noVectors)) endIndex = static_cast<int64_t>(m_noVectors);
        int64_t index=startIndex;
        while(index<endIndex)
        {
          if(isMonitor(m_readIndex))
          {     isisRaw->skipData(m_fileRaw,m_readIndex+1);
            //g_log.error()<<"skipData called for monitor index"<<m_readIndex<<std::endl;
            ++m_nmonitorSkipCounter;
            ++m_readIndex;
          }
          else
          {
            isisRaw->readData(m_fileRaw,m_readIndex+1);
            //g_log.error()<<"readData called for spectrum index"<<m_readIndex<< " and wsIndex is "<<index<< std::endl;
            if( m_readIndex == static_cast<int64_t>(m_noVectors+m_monitorList.size()) )
              break;

            // The managed histogram we are modifying
            ManagedHistogram1D * spec = dynamic_cast<ManagedHistogram1D *>(newBlock->getSpectrum(index));

            MantidVec& y = spec->directDataY();
            y.assign(isisRaw->dat1 + 1, isisRaw->dat1 + m_numberOfBinBoundaries);  
            //g_log.error()<<"readData called for m_readIndex"<<m_readIndex<< " and wsIndex is "<<index<< "Y value at 0  column is "<<y[0]<<std::endl;
            MantidVec& e = spec->directDataE();
            e.resize(y.size(), 0);
            std::transform(y.begin(), y.end(), e.begin(), dblSqrt);
            if (m_timeChannels.size() == 1)
              spec->directSetX(m_timeChannels[0]);
            else
            {
              // std::map<int,int>::const_iterator regime = m_specTimeRegimes.find(index+1);
              std::map<int64_t,int64_t>::const_iterator regime = m_specTimeRegimes.find(m_readIndex+1);
              if ( regime == m_specTimeRegimes.end() ) 
              {
                g_log.error() << "Spectrum " << index << " not present in spec array:\n";
                g_log.error(" Assuming time regime of spectrum 1");
                regime = m_specTimeRegimes.begin();
              }
              spec->directSetX(m_timeChannels[(*regime).second-1]);
            }
            spec->setLoaded(true);
            ++index;
            ++m_readIndex;
          }
        }
 
      }
      else
      { 
        if ( static_cast<int>(startIndex) > m_readIndex)
        {
          while( static_cast<int>(startIndex) > m_readIndex)
          {
            isisRaw->skipData(m_fileRaw,m_readIndex+1);// Adding 1 because we dropped the first spectrum.
            ++m_readIndex;
          }
        }
        else
        {
          int nwords = 0;
          while( static_cast<int>(startIndex) < m_readIndex)
          {                      
            --m_readIndex;
            nwords += 4*isisRaw->ddes[m_readIndex+1].nwords;
          }
          if (fseek(m_fileRaw,-nwords,SEEK_CUR) != 0)
          {
            fclose(m_fileRaw);
            removeTempFile();
            g_log.error("Error reading RAW file.");
            throw std::runtime_error("ManagedRawFileWorkspace2D: Error reading RAW file.");
          }
        }
        int64_t endIndex = startIndex+m_vectorsPerBlock < m_noVectors?startIndex+m_vectorsPerBlock:m_noVectors;
        if (endIndex >=  static_cast<int64_t>(m_noVectors)) endIndex = m_noVectors;
        for(int64_t index = startIndex;index<endIndex;index++,m_readIndex++)
        {
          isisRaw->readData(m_fileRaw,m_readIndex+1);
          // g_log.error()<<"counter is "<<counter<<std::endl;

          // The managed histogram we are modifying
          ManagedHistogram1D * spec = dynamic_cast<ManagedHistogram1D *>(newBlock->getSpectrum(index));

          MantidVec& y = spec->directDataY();
          y.assign(isisRaw->dat1 + 1, isisRaw->dat1 + m_numberOfBinBoundaries);   
          MantidVec& e = spec->directDataE();
          e.resize(y.size(), 0);
          std::transform(y.begin(), y.end(), e.begin(), dblSqrt);
          if (m_timeChannels.size() == 1)
            spec->directSetX(m_timeChannels[0]);
          else
          {
            std::map<int64_t,int64_t>::const_iterator regime = m_specTimeRegimes.find(index+1);
            if ( regime == m_specTimeRegimes.end() ) 
            {
              g_log.error() << "Spectrum " << index << " not present in spec array:\n";
              g_log.error(" Assuming time regime of spectrum 1");
              regime = m_specTimeRegimes.begin();
            }
            spec->directSetX(m_timeChannels[(*regime).second-1]);
          }
          spec->setLoaded(true);
        }
      }

      newBlock->hasChanges(false);
      newBlock->setLoaded(true);
    }
    /** This method checks given spectrum is a monitor
     * @param readIndex :: a spectrum index
     * @return true if it's a monitor ,otherwise false
     */
    bool ManagedRawFileWorkspace2D::isMonitor(const specid_t readIndex)const
    {
      std::vector<specid_t>::const_iterator itr;
      for(itr=m_monitorList.begin();itr!=m_monitorList.end();++itr)
      {
        if((*itr)==readIndex)
          return true;
      }
      return false;
    }

    void ManagedRawFileWorkspace2D::writeDataBlock(DataObjects::ManagedDataBlock2D *toWrite) const
    {
      Poco::ScopedLock<Poco::FastMutex> mutex(m_mutex);
      // ManagedWorkspace2D resets the hasChanges flag but we need to make sure we keep track of it here as well
      const bool blockHasChanged = toWrite->hasChanges();
      ManagedWorkspace2D::writeDataBlock(toWrite);
      int blockIndex = static_cast<int>(toWrite->minIndex() / m_vectorsPerBlock);
      m_changedBlock[blockIndex] = blockHasChanged;
    }

    /**
       Decides if the raw file must be copied to a cache file on the local drive to improve reading time.
    */
    bool ManagedRawFileWorkspace2D::needCache(const int opt)
    {
      if (opt == 1) return true;
      if (opt == 2) return false;

      return Kernel::ConfigService::Instance().isNetworkDrive(m_filenameRaw);
    }

    /**   Opens a temporary file
     */
    void ManagedRawFileWorkspace2D::openTempFile()
    {
      // Look for the (optional) path from the configuration file
      std::string path = Kernel::ConfigService::Instance().getString("ManagedWorkspace.FilePath");
      if( path.empty() || !Poco::File(path).exists() || !Poco::File(path).canWrite() )
      {
        path = Mantid::Kernel::ConfigService::Instance().getUserPropertiesDir();
        g_log.debug() << "Temporary file written to " << path << std::endl;
      }
      if ( ( *(path.rbegin()) != '/' ) && ( *(path.rbegin()) != '\\' ) )
      {
        path.push_back('/');
      }

      std::stringstream filename;
      filename << "WS2D_" << Poco::Path(m_filenameRaw).getBaseName() <<'_'<< ManagedRawFileWorkspace2D::g_uniqueID <<".raw";
      // Increment the instance count
      ++ManagedRawFileWorkspace2D::g_uniqueID;
      m_tempfile = path + filename.str();
      Poco::File(m_filenameRaw).copyTo(m_tempfile);

      FILE *fileRaw = fopen(m_tempfile.c_str(),"rb");
      if (fileRaw)
      {
        fclose(m_fileRaw);
        m_fileRaw = fileRaw;
      }

    }

    /**   Removes the temporary file
     */
    void ManagedRawFileWorkspace2D::removeTempFile() const
    {
      if (!m_tempfile.empty()) Poco::File(m_tempfile).remove();
    }
        

  } // namespace DataHandling
} //NamespaceMantid