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

/** Load in a single spectrum taken from a NeXus file
*  @param hist ::     The workspace index
*  @param i ::        The spectrum number
*  @param specNo ::   The spectrum number
*  @param nxload ::   A reference to the MuonNeXusReader object
*  @param lengthIn :: The number of elements in a spectrum
*  @param localWorkspace :: A pointer to the workspace in which the data will be
* stored
*/
void LoadMuonNexus1::loadData(size_t hist, specnum_t &i, specnum_t specNo,
                              MuonNexusReader &nxload, const int64_t lengthIn,
                              DataObjects::Workspace2D_sptr localWorkspace) {
  // Read in a spectrum
  // Put it into a vector, discarding the 1st entry, which is rubbish
  // But note that the last (overflow) bin is kept
  // For Nexus, not sure if above is the case, hence give all data for now
  MantidVec &Y = localWorkspace->dataY(hist);
  Y.assign(nxload.counts + i * lengthIn,
           nxload.counts + i * lengthIn + lengthIn);

  // Create and fill another vector for the errors, containing sqrt(count)
  MantidVec &E = localWorkspace->dataE(hist);
  typedef double (*uf)(double);
  uf dblSqrt = std::sqrt;
  std::transform(Y.begin(), Y.end(), E.begin(), dblSqrt);
  // Populate the workspace. Loop starts from 1, hence i-1

  // Create and fill another vector for the X axis
  auto timeChannels = new float[lengthIn + 1]();
  nxload.getTimeChannels(timeChannels, static_cast<const int>(lengthIn + 1));
  // Put the read in array into a vector (inside a shared pointer)
  auto timeChannelsVec = boost::make_shared<HistogramData::HistogramX>(
      timeChannels, timeChannels + lengthIn + 1);

  localWorkspace->setX(hist, timeChannelsVec);
  localWorkspace->getSpectrum(hist).setSpectrumNo(specNo);
  // Muon v1 files: always a one-to-one mapping between spectra and detectors
  localWorkspace->getSpectrum(hist).setDetectorID(static_cast<detid_t>(specNo));
  // Clean up
  delete[] timeChannels;
}

/**
* Perform a call to nxgetslab, via the NexusClasses wrapped methods for a given
* block-size
* @param data :: The NXDataSet object
* @param blocksize :: The block-size 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_loadBlockInfo.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_detBlockInfo.numberOfChannels;
    data_end += m_detBlockInfo.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);
    if (m_load_selected_spectra) {
      // local_workspace->getAxis(1)->setValue(hist,
      // static_cast<specid_t>(spec_num));
      auto spec = local_workspace->getSpectrum(hist);
      specid_t specID = m_specInd2specNum_map.at(hist);
      // set detectors corresponding to spectra Number
      spec->setDetectorIDs(m_spec2det_map.getDetectorIDsForSpectrumNo(specID));
      // set correct spectra Number
      spec->setSpectrumNo(specID);
    }

    ++hist;
    ++spec_num;
  }
}

/**
 * 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) {
  auto &X = ws->mutableX(specID);
  auto &Y = ws->mutableY(specID);
  auto &E = ws->mutableE(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->getSpectrum(specID).setSpectrumNo(specID);
}

/**
 * Populate spectra mapping to detector IDs
 *
 * TODO: Get the detector size information from the workspace directly
 *
 * @param localWorkspace: Workspace2D object
 * @param nxbins: number of bins in X
 * @param nybins: number of bins in Y
 */
void LoadSpice2D::runLoadMappingTable(
    DataObjects::Workspace2D_sptr localWorkspace, int nxbins, int nybins) {
  // Get the number of monitor channels
  boost::shared_ptr<const Geometry::Instrument> instrument =
      localWorkspace->getInstrument();
  std::vector<detid_t> monitors = instrument->getMonitors();
  const int nMonitors = static_cast<int>(monitors.size());

  // Number of monitors should be consistent with data file format
  if (nMonitors != LoadSpice2D::nMonitors) {
    std::stringstream error;
    error << "Geometry error for " << instrument->getName()
          << ": Spice data format defines " << LoadSpice2D::nMonitors
          << " monitors, " << nMonitors << " were/was found";
    throw std::runtime_error(error.str());
  }

  // Generate mapping of detector/channel IDs to spectrum ID

  // Detector/channel counter
  int icount = 0;

  // Monitor: IDs start at 1 and increment by 1
  for (int i = 0; i < nMonitors; i++) {
    localWorkspace->getSpectrum(icount)->setDetectorID(icount + 1);
    icount++;
  }

  // Detector pixels
  for (int ix = 0; ix < nxbins; ix++) {
    for (int iy = 0; iy < nybins; iy++) {
      localWorkspace->getSpectrum(icount)
          ->setDetectorID(1000000 + iy * 1000 + ix);
      icount++;
    }
  }
}

/**
* 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
* @param update_spectra2det_mapping :: reset spectra-detector map to the one
* calculated earlier. (Warning! -- this map has to be calculated correctly!)
*/
void
LoadISISNexus2::loadPeriodData(int64_t period, NXEntry &entry,
                               DataObjects::Workspace2D_sptr &local_workspace,
                               bool update_spectra2det_mapping) {
  int64_t hist_index = 0;
  int64_t period_index(period - 1);
  // int64_t first_monitor_spectrum = 0;

  for (auto block = m_spectraBlocks.begin(); block != m_spectraBlocks.end();
       ++block) {
    if (block->isMonitor) {
      NXData monitor = entry.openNXData(block->monName);
      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_monBlockInfo.numberOfChannels);
      MantidVec &E = local_workspace->dataE(hist_index);
      std::transform(Y.begin(), Y.end(), E.begin(), dblSqrt);

      if (update_spectra2det_mapping) {
        // local_workspace->getAxis(1)->setValue(hist_index,
        // static_cast<specid_t>(it->first));
        auto spec = local_workspace->getSpectrum(hist_index);
        specid_t specID = m_specInd2specNum_map.at(hist_index);
        spec->setDetectorIDs(
            m_spec2det_map.getDetectorIDsForSpectrumNo(specID));
        spec->setSpectrumNo(specID);
      }

      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 the list 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 block-size
      // 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;
  }
}