C++ Workspace2D_sptr::getAxis方法代码示例

/** Finalizes the calculation of the correlation spectrum
  *
  * This method offers a variable way of using the correlation spectrum
  *calculated previously.
  * The base version converts to Q and creates an appropriate output workspace.
  *
  * @param correctedCorrelatedIntensities :: Intensities of correlation
  *spectrum.
  * @param dValues :: d-spacings at which the spectrum was calculated.
  * @return A workspace containing the correlation spectrum.
  */
DataObjects::Workspace2D_sptr PoldiAutoCorrelationCore::finalizeCalculation(
    const std::vector<double> &correctedCorrelatedIntensities,
    const std::vector<double> &dValues) const {
  /* Finally, the d-Values are converted to q-Values for plotting etc. and
   * inserted into the output workspace. */
  size_t dCount = dValues.size();
  std::vector<double> qValues(dCount);

  PARALLEL_FOR_NO_WSP_CHECK()
  for (int i = 0; i < static_cast<int>(dCount); ++i) {
    qValues[dCount - i - 1] = Conversions::dToQ(dValues[i]);
  }

  m_logger.information() << "  Setting result..." << std::endl;
  DataObjects::Workspace2D_sptr outputWorkspace =
      boost::dynamic_pointer_cast<Mantid::DataObjects::Workspace2D>(
          WorkspaceFactory::Instance().create("Workspace2D", 1, dValues.size(),
                                              dValues.size()));

  outputWorkspace->getAxis(0)->setUnit("MomentumTransfer");

  outputWorkspace->dataY(0) = correctedCorrelatedIntensities;

  outputWorkspace->setX(0, qValues);

  return outputWorkspace;
}

/**
 * Create workspace to store the structure factor.
 * First spectrum is the real part, second spectrum is the imaginary part
 * X values are the modulus of the Q-vectors
 * @param h5file file identifier
 * @param gws pointer to WorkspaceGroup being filled
 * @param setName string name of dataset
 * @param qvmod vector of Q-vectors' moduli
 * @param sorting_indexes permutation of qvmod indexes to render it in increasing order of momemtum transfer
 */
void LoadSassena::loadFQ(const hid_t& h5file, API::WorkspaceGroup_sptr gws, const std::string setName, const MantidVec &qvmod, const std::vector<int> &sorting_indexes)
{
  const std::string gwsName = this->getPropertyValue("OutputWorkspace");
  int nq = static_cast<int>( qvmod.size() ); //number of q-vectors

  DataObjects::Workspace2D_sptr ws = boost::dynamic_pointer_cast<DataObjects::Workspace2D>(API::WorkspaceFactory::Instance().create("Workspace2D", 2, nq, nq));
  const std::string wsName = gwsName + std::string("_") + setName;
  ws->setTitle(wsName);

  double* buf = new double[nq*2];
  this->dataSetDouble(h5file,setName,buf);
  MantidVec& re = ws->dataY(0); // store the real part
  ws->dataX(0) = qvmod;  //X-axis values are the modulus of the q vector
  MantidVec& im = ws->dataY(1); // store the imaginary part
  ws->dataX(1) = qvmod;
  double *curr = buf;
  for(int iq=0; iq<nq; iq++){
    const int index=sorting_indexes[iq];
    re[index]=curr[0];
    im[index]=curr[1];
    curr+=2;
  }
  delete[] buf;

  // Set the Units
  ws->getAxis(0)->unit() = Kernel::UnitFactory::Instance().create("MomentumTransfer");

  this->registerWorkspace(gws,wsName,ws, "X-axis: Q-vector modulus; Y-axis: intermediate structure factor");
}

/**
 * load vectors onto a Workspace2D with 3 bins (the three components of the
 * vectors)
 * dataX for the origin of the vector (assumed (0,0,0) )
 * dataY for the tip of the vector
 * dataE is assumed (0,0,0), no errors
 * @param h5file file identifier
 * @param gws pointer to WorkspaceGroup being filled
 * @param sorting_indexes permutation of qvmod indexes to render it in
 * increasing order of momemtum transfer
 */
const MantidVec LoadSassena::loadQvectors(const hid_t &h5file,
                                          API::WorkspaceGroup_sptr gws,
                                          std::vector<int> &sorting_indexes) {

  const std::string gwsName = this->getPropertyValue("OutputWorkspace");
  const std::string setName("qvectors");

  hsize_t dims[3];
  if (dataSetInfo(h5file, setName, dims) < 0) {
    throw Kernel::Exception::FileError(
        "Unable to read " + setName + " dataset info:", m_filename);
  }
  int nq = static_cast<int>(dims[0]); // number of q-vectors
  double *buf = new double[nq * 3];
  this->dataSetDouble(h5file, "qvectors", buf);

  MantidVec qvmod; // store the modulus of the vector
  double *curr = buf;
  for (int iq = 0; iq < nq; iq++) {
    qvmod.push_back(
        sqrt(curr[0] * curr[0] + curr[1] * curr[1] + curr[2] * curr[2]));
    curr += 3;
  }

  if (getProperty("SortByQVectors")) {
    std::vector<mypair> qvmodpair;
    for (int iq = 0; iq < nq; iq++)
      qvmodpair.push_back(mypair(qvmod[iq], iq));
    std::sort(qvmodpair.begin(), qvmodpair.end(), compare);
    for (int iq = 0; iq < nq; iq++)
      sorting_indexes.push_back(qvmodpair[iq].second);
    std::sort(qvmod.begin(), qvmod.end());
  } else
    for (int iq = 0; iq < nq; iq++)
      sorting_indexes.push_back(iq);

  DataObjects::Workspace2D_sptr ws =
      boost::dynamic_pointer_cast<DataObjects::Workspace2D>(
          API::WorkspaceFactory::Instance().create("Workspace2D", nq, 3, 3));
  std::string wsName = gwsName + std::string("_") + setName;
  ws->setTitle(wsName);

  for (int iq = 0; iq < nq; iq++) {
    MantidVec &Y = ws->dataY(iq);
    const int index = sorting_indexes[iq];
    curr = buf + 3 * index;
    Y.assign(curr, curr + 3);
  }
  delete[] buf;

  ws->getAxis(0)->unit() = Kernel::UnitFactory::Instance().create(
      "MomentumTransfer"); // Set the Units

  this->registerWorkspace(
      gws, wsName, ws, "X-axis: origin of Q-vectors; Y-axis: tip of Q-vectors");
  return qvmod;
}

 /**
  * Convenience function to store a detector value into a given spectrum.
  * Note that this type of data doesn't use TOD, so that we use a single dummy
  * bin in X. Each detector is defined as a spectrum of length 1.
  * @param ws: workspace
  * @param specID: ID of the spectrum to store the value in
  * @param value: value to store [count]
  * @param error: error on the value [count]
  * @param wavelength: wavelength value [Angstrom]
  * @param dwavelength: error on the wavelength [Angstrom]
  */
 void store_value(DataObjects::Workspace2D_sptr ws, int specID,
     double value, double error, double wavelength, double dwavelength)
 {
   MantidVec& X = ws->dataX(specID);
   MantidVec& Y = ws->dataY(specID);
   MantidVec& E = ws->dataE(specID);
   // The following is mostly to make Mantid happy by defining a histogram with
   // a single bin around the neutron wavelength
   X[0] = wavelength-dwavelength/2.0;
   X[1] = wavelength+dwavelength/2.0;
   Y[0] = value;
   E[0] = error;
   ws->getAxis(1)->setValue(specID, specID);
 }

    /** Read in a single spectrum from the raw file
     *  @param tcbs ::     The vector containing the time bin boundaries
     *  @param hist ::     The workspace index
     *  @param i ::        The spectrum number
     *  @param iraw ::     A reference to the ISISRAW object
     *  @param lengthIn :: The number of elements in a spectrum
     *  @param spectrum :: Pointer to the array into which the spectrum will be read
     *  @param localWorkspace :: A pointer to the workspace in which the data will be stored
     */
    void LoadRaw::loadData(const MantidVecPtr::ptr_type& tcbs, int32_t hist, specid_t& i, ISISRAW& iraw, const int& lengthIn, int* spectrum, DataObjects::Workspace2D_sptr localWorkspace)
    {
      // Read in a spectrum
      memcpy(spectrum, iraw.dat1 + i * lengthIn, lengthIn * sizeof(int));
      // Put it into a vector, discarding the 1st entry, which is rubbish
      // But note that the last (overflow) bin is kept
      MantidVec& Y = localWorkspace->dataY(hist);
      Y.assign(spectrum + 1, spectrum + lengthIn);
      // Create and fill another vector for the errors, containing sqrt(count)
      MantidVec& E = localWorkspace->dataE(hist);
      std::transform(Y.begin(), Y.end(), E.begin(), dblSqrt);
      // Populate the workspace. Loop starts from 1, hence i-1
      localWorkspace->setX(hist, tcbs);

      localWorkspace->getAxis(1)->setValue(hist, i);
      // NOTE: Raw numbers go straight into the workspace
      //     - no account taken of bin widths/units etc.
    }

 /**
 * Perform a call to nxgetslab, via the NexusClasses wrapped methods for a given blocksize
 * @param data :: The NXDataSet object
 * @param blocksize :: The blocksize to use
 * @param period :: The period number
 * @param start :: The index within the file to start reading from (zero based)
 * @param hist :: The workspace index to start reading into
 * @param spec_num :: The spectrum number that matches the hist variable
 * @param local_workspace :: The workspace to fill the data with
 */
 void LoadISISNexus2::loadBlock(NXDataSetTyped<int> & data, int64_t blocksize, int64_t period, int64_t start,
     int64_t &hist, int64_t& spec_num,
   DataObjects::Workspace2D_sptr local_workspace)
 {
   data.load(static_cast<int>(blocksize), static_cast<int>(period), static_cast<int>(start)); // TODO this is just wrong
   int *data_start = data();
   int *data_end = data_start + m_numberOfChannels;
   int64_t final(hist + blocksize);
   while( hist < final )
   {
     m_progress->report("Loading data");
     MantidVec& Y = local_workspace->dataY(hist);
     Y.assign(data_start, data_end);
     data_start += m_numberOfChannels; data_end += m_numberOfChannels;
     MantidVec& E = local_workspace->dataE(hist);
     std::transform(Y.begin(), Y.end(), E.begin(), dblSqrt);
     // Populate the workspace. Loop starts from 1, hence i-1
     local_workspace->setX(hist, m_tof_data);
     local_workspace->getAxis(1)->spectraNo(hist)= static_cast<specid_t>(spec_num);
     ++hist;
     ++spec_num;
   }
 }

    /**
    * Load a given period into the workspace
    * @param period :: The period number to load (starting from 1) 
    * @param entry :: The opened root entry node for accessing the monitor and data nodes
    * @param local_workspace :: The workspace to place the data in
    */
    void LoadISISNexus2::loadPeriodData(int64_t period, NXEntry & entry, DataObjects::Workspace2D_sptr local_workspace)
    {
      int64_t hist_index = 0;
      int64_t period_index(period - 1);
      int64_t first_monitor_spectrum = 0;

      if( !m_monitors.empty() )
      {
        first_monitor_spectrum = m_monitors.begin()->first;
        hist_index = first_monitor_spectrum - 1;
        for(std::map<int64_t,std::string>::const_iterator it = m_monitors.begin();
          it != m_monitors.end(); ++it)
        {
          NXData monitor = entry.openNXData(it->second);
          NXInt mondata = monitor.openIntData();
          m_progress->report("Loading monitor");
          mondata.load(1,static_cast<int>(period-1)); // TODO this is just wrong
          MantidVec& Y = local_workspace->dataY(hist_index);
          Y.assign(mondata(),mondata() + m_numberOfChannels);
          MantidVec& E = local_workspace->dataE(hist_index);
          std::transform(Y.begin(), Y.end(), E.begin(), dblSqrt);
          local_workspace->getAxis(1)->spectraNo(hist_index) = static_cast<specid_t>(it->first);

          NXFloat timeBins = monitor.openNXFloat("time_of_flight");
          timeBins.load();
          local_workspace->dataX(hist_index).assign(timeBins(),timeBins() + timeBins.dim0());
          hist_index++;
        }

        if (first_monitor_spectrum > 1)
        {
          hist_index = 0;
        }
      }
      
      if( m_have_detector )
      {
        NXData nxdata = entry.openNXData("detector_1");
        NXDataSetTyped<int> data = nxdata.openIntData();
        data.open();
        //Start with thelist members that are lower than the required spectrum
        const int * const spec_begin = m_spec.get();
        std::vector<int64_t>::iterator min_end = m_spec_list.end();
        if( !m_spec_list.empty() )
        {
          // If we have a list, by now it is ordered so first pull in the range below the starting block range
          // Note the reverse iteration as we want the last one
          if( m_range_supplied )
          {
            min_end = std::find_if(m_spec_list.begin(), m_spec_list.end(), std::bind2nd(std::greater<int>(), m_spec_min));
          }

          for( std::vector<int64_t>::iterator itr = m_spec_list.begin(); itr < min_end; ++itr )
          {
            // Load each
            int64_t spectra_no = (*itr);
            // For this to work correctly, we assume that the spectrum list increases monotonically
            int64_t filestart = std::lower_bound(spec_begin,m_spec_end,spectra_no) - spec_begin;
            m_progress->report("Loading data");
            loadBlock(data, static_cast<int64_t>(1), period_index, filestart, hist_index, spectra_no, local_workspace);
          }
        }    

        if( m_range_supplied )
        {
          // When reading in blocks we need to be careful that the range is exactly divisible by the blocksize
          // and if not have an extra read of the left overs
          const int64_t blocksize = 8;
          const int64_t rangesize = (m_spec_max - m_spec_min + 1) - m_monitors.size();
          const int64_t fullblocks = rangesize / blocksize;
          int64_t read_stop = 0;
          int64_t spectra_no = m_spec_min;
          if (first_monitor_spectrum == 1)
          {// this if crudely checks whether the monitors are at the begining or end of the spectra
            spectra_no += static_cast<int>(m_monitors.size());
          }
          // For this to work correctly, we assume that the spectrum list increases monotonically
          int64_t filestart = std::lower_bound(spec_begin,m_spec_end,spectra_no) - spec_begin;
          if( fullblocks > 0 )
          {
            read_stop = (fullblocks * blocksize);// + m_monitors.size(); //RNT: I think monitors are excluded from the data
            //for( ; hist_index < read_stop; )
            for(int64_t i = 0; i < fullblocks; ++i)
            {
              loadBlock(data, blocksize, period_index, filestart, hist_index, spectra_no, local_workspace);
              filestart += blocksize;
            }
          }
          int64_t finalblock = rangesize - (fullblocks * blocksize);
          if( finalblock > 0 )
          {
            loadBlock(data, finalblock, period_index, filestart, hist_index, spectra_no,  local_workspace);
          }
        }
//.........这里部分代码省略.........

/**
 * Create one workspace to hold the real part and another to hold the imaginary
* part.
 * We symmetrize the structure factor to negative times
 * Y-values are structure factor for each Q-value
 * X-values are time bins
 * @param h5file file identifier
 * @param gws pointer to WorkspaceGroup being filled
 * @param setName string name of dataset
 * @param qvmod vector of Q-vectors' moduli
* @param sorting_indexes permutation of qvmod indexes to render it in increasing
* order of momemtum transfer
*/
void LoadSassena::loadFQT(const hid_t &h5file, API::WorkspaceGroup_sptr gws,
                          const std::string setName, const MantidVec &qvmod,
                          const std::vector<int> &sorting_indexes) {
  const std::string gwsName = this->getPropertyValue("OutputWorkspace");
  int nq = static_cast<int>(qvmod.size()); // number of q-vectors
  const double dt =
      getProperty("TimeUnit"); // time unit increment, in picoseconds;
  hsize_t dims[3];
  if (dataSetInfo(h5file, setName, dims) < 0) {
    throw Kernel::Exception::FileError(
        "Unable to read " + setName + " dataset info:", m_filename);
  }
  int nnt = static_cast<int>(dims[1]); // number of non-negative time points
  int nt = 2 * nnt - 1;                // number of time points
  int origin = nnt - 1;
  double *buf = new double[nq * nnt * 2];
  this->dataSetDouble(h5file, setName, buf);

  DataObjects::Workspace2D_sptr wsRe =
      boost::dynamic_pointer_cast<DataObjects::Workspace2D>(
          API::WorkspaceFactory::Instance().create("Workspace2D", nq, nt, nt));
  const std::string wsReName =
      gwsName + std::string("_") + setName + std::string(".Re");
  wsRe->setTitle(wsReName);

  DataObjects::Workspace2D_sptr wsIm =
      boost::dynamic_pointer_cast<DataObjects::Workspace2D>(
          API::WorkspaceFactory::Instance().create("Workspace2D", nq, nt, nt));
  const std::string wsImName =
      gwsName + std::string("_") + setName + std::string(".Im");
  wsIm->setTitle(wsImName);

  for (int iq = 0; iq < nq; iq++) {
    MantidVec &reX = wsRe->dataX(iq);
    MantidVec &imX = wsIm->dataX(iq);
    MantidVec &reY = wsRe->dataY(iq);
    MantidVec &imY = wsIm->dataY(iq);
    const int index = sorting_indexes[iq];
    double *curr = buf + index * nnt * 2;
    for (int it = 0; it < nnt; it++) {
      reX[origin + it] = it * dt; // time point for the real part
      reY[origin + it] =
          *curr; // real part of the intermediate structure factor
      reX[origin - it] = -it * dt; // symmetric negative time
      reY[origin - it] = *curr;    // symmetric value for the negative time
      curr++;
      imX[origin + it] = it * dt;
      imY[origin + it] = *curr;
      imX[origin - it] = -it * dt;
      imY[origin - it] = -(*curr); // antisymmetric value for negative times
      curr++;
    }
  }
  delete[] buf;

  // Set the Time unit for the X-axis
  wsRe->getAxis(0)->unit() = Kernel::UnitFactory::Instance().create("Label");
  auto unitPtr = boost::dynamic_pointer_cast<Kernel::Units::Label>(
      wsRe->getAxis(0)->unit());
  unitPtr->setLabel("Time", "picoseconds");

  wsIm->getAxis(0)->unit() = Kernel::UnitFactory::Instance().create("Label");
  unitPtr = boost::dynamic_pointer_cast<Kernel::Units::Label>(
      wsIm->getAxis(0)->unit());
  unitPtr->setLabel("Time", "picoseconds");

  // Create a numeric axis to replace the default vertical one
  API::Axis *const verticalAxisRe = new API::NumericAxis(nq);
  API::Axis *const verticalAxisIm = new API::NumericAxis(nq);

  wsRe->replaceAxis(1, verticalAxisRe);
  wsIm->replaceAxis(1, verticalAxisIm);

  // Now set the axis values
  for (int i = 0; i < nq; ++i) {
    verticalAxisRe->setValue(i, qvmod[i]);
    verticalAxisIm->setValue(i, qvmod[i]);
  }

  // Set the axis units
  verticalAxisRe->unit() =
      Kernel::UnitFactory::Instance().create("MomentumTransfer");
  verticalAxisRe->title() = "|Q|";
  verticalAxisIm->unit() =
      Kernel::UnitFactory::Instance().create("MomentumTransfer");
  verticalAxisIm->title() = "|Q|";

//.........这里部分代码省略.........

    /** Executes the algorithm. Reading in the file and creating and populating
     *  the output workspace
     *
     *  @throw Exception::FileError If the RAW file cannot be found/opened
     *  @throw std::invalid_argument If the optional properties are set to invalid values
     */
    void LoadRaw::exec()
    {
      // Retrieve the filename from the properties
      m_filename = getPropertyValue("Filename");

      LoadRawHelper *helper = new LoadRawHelper;
      FILE* file = helper->openRawFile(m_filename);
      ISISRAW iraw;
      iraw.ioRAW(file, true);

      std::string title(iraw.r_title, 80);
      g_log.information("**** Run title: "+title+ "***");
      
      // Read in the number of spectra in the RAW file
      m_numberOfSpectra = iraw.t_nsp1;
      // Read the number of periods in this file
      m_numberOfPeriods = iraw.t_nper;
      // Need to extract the user-defined output workspace name
      Property *ws = getProperty("OutputWorkspace");
      std::string localWSName = ws->value();

      // Call private method to validate the optional parameters, if set
      checkOptionalProperties();

      // Read the number of time channels (i.e. bins) from the RAW file
      const int channelsPerSpectrum = iraw.t_ntc1;
      // Read in the time bin boundaries
      const int lengthIn = channelsPerSpectrum + 1;
      float* timeChannels = new float[lengthIn];
      iraw.getTimeChannels(timeChannels, lengthIn);
      // Put the read in array into a vector (inside a shared pointer)
      boost::shared_ptr<MantidVec> timeChannelsVec
                          (new MantidVec(timeChannels, timeChannels + lengthIn));
      // Create an array to hold the read-in data
      int* spectrum = new int[lengthIn];

      // Calculate the size of a workspace, given its number of periods & spectra to read
      specid_t total_specs;
      if( m_interval || m_list)
      {
        total_specs = static_cast<specid_t>(m_spec_list.size());
        if (m_interval)
        {
          total_specs += (m_spec_max-m_spec_min+1);
          m_spec_max += 1;
        }
      }
      else
      {
        total_specs = m_numberOfSpectra;
        // In this case want all the spectra, but zeroth spectrum is garbage so go from 1 to NSP1
        m_spec_min = 1;
        m_spec_max = m_numberOfSpectra + 1;
      }

      double histTotal = static_cast<double>(total_specs * m_numberOfPeriods);
      int32_t histCurrent = -1;

      // Create the 2D workspace for the output
      DataObjects::Workspace2D_sptr localWorkspace = boost::dynamic_pointer_cast<DataObjects::Workspace2D>
               (WorkspaceFactory::Instance().create("Workspace2D",total_specs,lengthIn,lengthIn-1));
      localWorkspace->getAxis(0)->unit() = UnitFactory::Instance().create("TOF");
      localWorkspace->setTitle(title);
      // Run parameters
      helper->loadRunParameters(localWorkspace, &iraw);
      delete helper;
      helper = NULL;


      // Loop over the number of periods in the raw file, putting each period in a separate workspace
      for (int period = 0; period < m_numberOfPeriods; ++period) {

        if ( period > 0 ) localWorkspace =  boost::dynamic_pointer_cast<DataObjects::Workspace2D>
                                              (WorkspaceFactory::Instance().create(localWorkspace));

        specid_t counter = 0;
        for (specid_t i = m_spec_min; i < m_spec_max; ++i)
        {
          // Shift the histogram to read if we're not in the first period
          int32_t histToRead = i + period*total_specs;
          loadData(timeChannelsVec,counter,histToRead,iraw,lengthIn,spectrum,localWorkspace );
          counter++;
          if (++histCurrent % 100 == 0) progress(static_cast<double>(histCurrent)/histTotal);
          interruption_point();
        }
        // Read in the spectra in the optional list parameter, if set
        if (m_list)
        {
          for(size_t i=0; i < m_spec_list.size(); ++i)
          {
            loadData(timeChannelsVec,counter,m_spec_list[i],iraw,lengthIn,spectrum, localWorkspace );
            counter++;
            if (++histCurrent % 100 == 0) progress(static_cast<double>(histCurrent)/histTotal);
            interruption_point();
//.........这里部分代码省略.........

//.........这里部分代码省略.........
      int numberXPixels = 0;
      int numberYPixels = 0;
      std::string data_type = element->getAttribute("type");
      boost::regex b_re_sig("INT\\d+\\[(\\d+),(\\d+)\\]");
      if (boost::regex_match(data_type, b_re_sig))
      {
        boost::match_results<std::string::const_iterator> match;
        boost::regex_search(data_type, match, b_re_sig);
        // match[0] is the full string
        Kernel::Strings::convert(match[1], numberXPixels);
        Kernel::Strings::convert(match[2], numberYPixels);
      }
      if (numberXPixels==0 || numberYPixels==0)
        g_log.notice() << "Could not read in the number of pixels!" << std::endl;

      // We no longer read from the meta data because that data is wrong
      //from_element<int>(numberXPixels, sasEntryElem, "Number_of_X_Pixels", fileName);
      //from_element<int>(numberYPixels, sasEntryElem, "Number_of_Y_Pixels", fileName);

      // Store sample-detector distance
      declareProperty("SampleDetectorDistance", distance, Kernel::Direction::Output);

      // Create the output workspace

      // Number of bins: we use a single dummy TOF bin
      int nBins = 1;
      // Number of detectors: should be pulled from the geometry description. Use detector pixels for now.
      // The number of spectram also includes the monitor and the timer.
      int numSpectra = numberXPixels*numberYPixels + LoadSpice2D::nMonitors;

      DataObjects::Workspace2D_sptr ws = boost::dynamic_pointer_cast<DataObjects::Workspace2D>(
          API::WorkspaceFactory::Instance().create("Workspace2D", numSpectra, nBins+1, nBins));
      ws->setTitle(wsTitle);
      ws->getAxis(0)->unit() = Kernel::UnitFactory::Instance().create("Wavelength");
      ws->setYUnit("");
      API::Workspace_sptr workspace = boost::static_pointer_cast<API::Workspace>(ws);
      setProperty("OutputWorkspace", workspace);

      // Parse out each pixel. Pixels can be separated by white space, a tab, or an end-of-line character
      Poco::StringTokenizer pixels(data_str, " \n\t", Poco::StringTokenizer::TOK_TRIM | Poco::StringTokenizer::TOK_IGNORE_EMPTY);
      Poco::StringTokenizer::Iterator pixel = pixels.begin();

      // Check that we don't keep within the size of the workspace
      size_t pixelcount = pixels.count();
      if( pixelcount != static_cast<size_t>(numberXPixels*numberYPixels) )
      {
        throw Kernel::Exception::FileError("Inconsistent data set: "
            "There were more data pixels found than declared in the Spice XML meta-data.", fileName);
      }
      if( numSpectra == 0 )
      {
        throw Kernel::Exception::FileError("Empty data set: the data file has no pixel data.", fileName);
      }

      // Go through all detectors/channels
      int ipixel = 0;

      // Store monitor count
      store_value(ws, ipixel++, monitorCounts, monitorCounts>0 ? sqrt(monitorCounts) : 0.0,
          wavelength, dwavelength);

      // Store counting time
      store_value(ws, ipixel++, countingTime, 0.0, wavelength, dwavelength);

      // Store detector pixels
      while (pixel != pixels.end())

/// Execute the algorithm
void SassenaFFT::exec()
{
  const std::string gwsName = this->getPropertyValue("InputWorkspace");
  API::WorkspaceGroup_sptr gws = this->getProperty("InputWorkspace");

  const std::string ftqReName = gwsName + "_fqt.Re";
  const std::string ftqImName = gwsName + "_fqt.Im";

  // Make sure the intermediate structure factor is there
  if(!gws->contains(ftqReName) )
  {
    const std::string errMessg = "workspace "+gwsName+" does not contain an intermediate structure factor";
    this->g_log.error(errMessg);
    throw Kernel::Exception::NotFoundError("group workspace does not contain",ftqReName);
  }

  // Retrieve the real and imaginary parts of the intermediate scattering function
  DataObjects::Workspace2D_sptr fqtRe = boost::dynamic_pointer_cast<DataObjects::Workspace2D>( gws->getItem( ftqReName ) );
  DataObjects::Workspace2D_sptr fqtIm = boost::dynamic_pointer_cast<DataObjects::Workspace2D>( gws->getItem( ftqImName ) );

  // Calculate the FFT for all spectra, retaining only the real part since F(q,-t) = F*(q,t)
  int part=3; // extract the real part of the transform, assuming I(Q,t) is real
  const std::string sqwName = gwsName + "_sqw";
  API::IAlgorithm_sptr fft = this->createChildAlgorithm("ExtractFFTSpectrum");
  fft->setProperty<DataObjects::Workspace2D_sptr>("InputWorkspace", fqtRe);
  if( !this->getProperty("FFTonlyRealPart") )
  {
    part=0; // extract the real part of the transform, assuming I(Q,t) is complex
    fft->setProperty<DataObjects::Workspace2D_sptr>("InputImagWorkspace", fqtIm);
  }
  fft->setPropertyValue("OutputWorkspace", sqwName );
  fft->setProperty<int>("FFTPart",part); // extract the real part
  fft->executeAsChildAlg();
  API::MatrixWorkspace_sptr sqw0 = fft->getProperty("OutputWorkspace");
  DataObjects::Workspace2D_sptr sqw = boost::dynamic_pointer_cast<DataObjects::Workspace2D>( sqw0 );
  API::AnalysisDataService::Instance().addOrReplace( sqwName, sqw );

  // Transform the X-axis to appropriate dimensions
  // We assume the units of the intermediate scattering function are in picoseconds
  // The resulting frequency unit is in mili-eV, thus use m_ps2meV
  API::IAlgorithm_sptr scaleX = this->createChildAlgorithm("ScaleX");
  scaleX->setProperty<DataObjects::Workspace2D_sptr>("InputWorkspace",sqw);
  scaleX->setProperty<double>("Factor", m_ps2meV);
  scaleX->setProperty<DataObjects::Workspace2D_sptr>("OutputWorkspace", sqw);
  scaleX->executeAsChildAlg();

  //Do we apply the detailed balance condition exp(E/(2*kT)) ?
  if( this->getProperty("DetailedBalance") )
  {
    double T = this->getProperty("Temp");
    // The ExponentialCorrection algorithm assumes the form C0*exp(-C1*x). Note the explicit minus in the exponent
    API::IAlgorithm_sptr ec = this->createChildAlgorithm("ExponentialCorrection");
    ec->setProperty<DataObjects::Workspace2D_sptr>("InputWorkspace", sqw);
    ec->setProperty<DataObjects::Workspace2D_sptr>("OutputWorkspace", sqw);
    ec->setProperty<double>("C0",1.0);
    ec->setProperty<double>("C1",-1.0/(2.0*T*m_T2ueV)); // Temperature in units of ueV
    ec->setPropertyValue("Operation","Multiply");
    ec->executeAsChildAlg();
  }

  // Set the Energy unit for the X-axis
  sqw->getAxis(0)->unit() = Kernel::UnitFactory::Instance().create("DeltaE");

  // Add to group workspace, except if we are replacing the workspace. In this case, the group workspace
  // is already notified of the changes by the analysis data service.
  if(!gws->contains(sqwName))
  {
    gws->add( sqwName );
  }
  else
  {
    this->g_log.information("Workspace "+sqwName+" replaced with new contents");
  }

}

    /**
    * Load a given period into the workspace
    * @param period :: The period number to load (starting from 1) 
    * @param entry :: The opened root entry node for accessing the monitor and data nodes
    * @param local_workspace :: The workspace to place the data in
    */
    void LoadISISNexus2::loadPeriodData(int64_t period, NXEntry & entry, DataObjects::Workspace2D_sptr local_workspace)
    {
      int64_t hist_index = 0;
      int64_t period_index(period - 1);
      //int64_t first_monitor_spectrum = 0;

      for(auto block = m_spectraBlocks.cbegin(); block != m_spectraBlocks.cend(); ++block)
      {
        if ( block->isMonitor )
        {
          auto it = m_monitors.find( block->first );
          assert( it != m_monitors.end() );
          NXData monitor = entry.openNXData(it->second);
          NXInt mondata = monitor.openIntData();
          m_progress->report("Loading monitor");
          mondata.load(1,static_cast<int>(period-1)); // TODO this is just wrong
          std::cerr << "monitor " << monitor.name() << std::endl;
          MantidVec& Y = local_workspace->dataY(hist_index);
          Y.assign(mondata(),mondata() + m_numberOfChannels);
          MantidVec& E = local_workspace->dataE(hist_index);
          std::transform(Y.begin(), Y.end(), E.begin(), dblSqrt);
          local_workspace->getAxis(1)->spectraNo(hist_index) = static_cast<specid_t>(it->first);

          NXFloat timeBins = monitor.openNXFloat("time_of_flight");
          timeBins.load();
          local_workspace->dataX(hist_index).assign(timeBins(),timeBins() + timeBins.dim0());
          hist_index++;
        }
        else if( m_have_detector )
        {
          NXData nxdata = entry.openNXData("detector_1");
          NXDataSetTyped<int> data = nxdata.openIntData();
          data.open();
          //Start with thelist members that are lower than the required spectrum
          const int * const spec_begin = m_spec.get();
          // When reading in blocks we need to be careful that the range is exactly divisible by the blocksize
          // and if not have an extra read of the left overs
          const int64_t blocksize = 8;
          const int64_t rangesize = block->last - block->first + 1;
          const int64_t fullblocks = rangesize / blocksize;
          int64_t spectra_no = block->first;

          // For this to work correctly, we assume that the spectrum list increases monotonically
          int64_t filestart = std::lower_bound(spec_begin,m_spec_end,spectra_no) - spec_begin;
          if( fullblocks > 0 )
          {
            for(int64_t i = 0; i < fullblocks; ++i)
            {
              loadBlock(data, blocksize, period_index, filestart, hist_index, spectra_no, local_workspace);
              filestart += blocksize;
            }
          }
          int64_t finalblock = rangesize - (fullblocks * blocksize);
          if( finalblock > 0 )
          {
            loadBlock(data, finalblock, period_index, filestart, hist_index, spectra_no,  local_workspace);
          }
        }
      }

      try
      {
        const std::string title = entry.getString("title");
        local_workspace->setTitle(title);
        // write the title into the log file (run object)
        local_workspace->mutableRun().addProperty("run_title", title, true);
      }
      catch (std::runtime_error &)
      {
        g_log.debug() << "No title was found in the input file, " << getPropertyValue("Filename") << std::endl;
      }
    }

    /** Executes the algorithm. Reading in the file and creating and populating
     *  the output workspace
     *
     *  @throw Exception::FileError If the RAW file cannot be found/opened
     *  @throw std::invalid_argument If the optional properties are set to invalid values
     */
    void LoadRaw2::exec()
    {
      // Retrieve the filename from the properties
      m_filename = getPropertyValue("Filename");
      LoadRawHelper *helper = new LoadRawHelper;
      FILE* file = helper->openRawFile(m_filename);
      isisRaw->ioRAW(file, true);
      
      std::string title(isisRaw->r_title, 80);
      g_log.information("**** Run title: "+ title + "***");

      // Read in the number of spectra in the RAW file
      m_numberOfSpectra = isisRaw->t_nsp1;
      // Read the number of periods in this file
      m_numberOfPeriods = isisRaw->t_nper;
      // Read the number of time channels (i.e. bins) from the RAW file
      const int channelsPerSpectrum = isisRaw->t_ntc1;
      // Read in the time bin boundaries
      const int lengthIn = channelsPerSpectrum + 1;

      // Call private method to validate the optional parameters, if set
      checkOptionalProperties();

      // Calculate the size of a workspace, given its number of periods & spectra to read
      specid_t total_specs;
      if( m_interval || m_list)
      {
        if (m_interval)
        {
          total_specs = (m_spec_max-m_spec_min+1);
          m_spec_max += 1;
        }
        else
            total_specs = 0;

        if (m_list)
        {
            if (m_interval)
            {
                for(std::vector<specid_t>::iterator it=m_spec_list.begin();it!=m_spec_list.end();)
                    if (*it >= m_spec_min && *it <m_spec_max)
                    {
                        it = m_spec_list.erase(it);
                    }
                    else
                        ++it;

            }
            if (m_spec_list.size() == 0) m_list = false;
            total_specs += static_cast<specid_t>(m_spec_list.size());
        }
      }
      else
      {
        total_specs = m_numberOfSpectra;
        // In this case want all the spectra, but zeroth spectrum is garbage so go from 1 to NSP1
        m_spec_min = 1;
        m_spec_max = m_numberOfSpectra + 1;
      }

      // If there is not enough memory use ManagedRawFileWorkspace2D.
      if ( ConfigService::Instance().getString("ManagedRawFileWorkspace.DoNotUse") != "1" &&
           m_numberOfPeriods == 1 && MemoryManager::Instance().goForManagedWorkspace(total_specs,lengthIn,channelsPerSpectrum) &&
          total_specs == m_numberOfSpectra)
      {
        const std::string cache_option = getPropertyValue("Cache");
        size_t option = find(m_cache_options.begin(),m_cache_options.end(),cache_option) - m_cache_options.begin();
        DataObjects::Workspace2D_sptr localWorkspace =
          DataObjects::Workspace2D_sptr(
          new ManagedRawFileWorkspace2D(m_filename, static_cast<int>(option)));
        progress(0.,"Reading raw file...");
        helper->loadRunParameters(localWorkspace, isisRaw.get());
        runLoadInstrument(localWorkspace );
        runLoadMappingTable(localWorkspace );
        runLoadLog(localWorkspace );
        const int period_number = 1;
        Property* log=createPeriodLog(period_number);
        if(log)
        {
          localWorkspace->mutableRun().addLogData(log);
          localWorkspace->mutableRun().addLogData(createCurrentPeriodLog(period_number));
        }
              localWorkspace->mutableRun().setProtonCharge(isisRaw->rpb.r_gd_prtn_chrg);
        for (int i = 0; i < m_numberOfSpectra; ++i)
          localWorkspace->getAxis(1)->setValue(i, i+1);
        localWorkspace->populateInstrumentParameters();
        setProperty("OutputWorkspace",localWorkspace);
        return;
      }

      float* timeChannels = new float[lengthIn];
      isisRaw->getTimeChannels(timeChannels, lengthIn);
      // Put the read in array into a vector (inside a shared pointer)
      boost::shared_ptr<MantidVec> timeChannelsVec
//.........这里部分代码省略.........