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

/**
* 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;
  }
}

/** Load in a single spectrum taken from a NeXus file
*  @param hist ::     The workspace index
*  @param i ::        The spectrum index
*  @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, specid_t &i, specid_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  
  float *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)
  boost::shared_ptr<MantidVec> timeChannelsVec(
    new MantidVec(timeChannels, timeChannels + lengthIn+1));

  localWorkspace->setX(hist, timeChannelsVec);
  localWorkspace->getSpectrum(hist)->setSpectrumNo(specNo);

  // Clean up
  delete[] timeChannels;

}

/**
 * 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->getSpectrum(specID)->setSpectrumNo(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;
   }
 }

/** Select background points
  */
void ProcessBackground::selectFromGivenXValues() {
  // Get special input properties
  std::vector<double> bkgdpoints = getProperty("BackgroundPoints");
  string mode = getProperty("BackgroundPointSelectMode");

  // Construct background workspace for fit
  std::vector<double> realx, realy, reale;
  const MantidVec &vecX = m_dataWS->readX(m_wsIndex);
  const MantidVec &vecY = m_dataWS->readY(m_wsIndex);
  const MantidVec &vecE = m_dataWS->readE(m_wsIndex);
  for (size_t i = 0; i < bkgdpoints.size(); ++i) {
    // Data range validation
    double bkgdpoint = bkgdpoints[i];
    if (bkgdpoint < vecX.front()) {
      g_log.warning() << "Input background point " << bkgdpoint
                      << " is out of lower boundary.  "
                      << "Use X[0] = " << vecX.front() << " instead."
                      << "\n";
      bkgdpoint = vecX.front();
    } else if (bkgdpoint > vecX.back()) {
      g_log.warning() << "Input background point " << bkgdpoint
                      << " is out of upper boundary.  Use X[-1] = "
                      << vecX.back() << " instead."
                      << "\n";
      bkgdpoint = vecX.back();
    }

    // Find the index in
    std::vector<double>::const_iterator it;
    it = std::lower_bound(vecX.begin(), vecX.end(), bkgdpoint);
    size_t index = size_t(it - vecX.begin());

    g_log.debug() << "DBx502 Background Points " << i << " Index = " << index
                  << " For TOF = " << bkgdpoints[i] << " in [" << vecX[0]
                  << ", " << vecX.back() << "] "
                  << "\n";

    // Add to list
    realx.push_back(vecX[index]);
    realy.push_back(vecY[index]);
    reale.push_back(vecE[index]);

  } // ENDFOR (i)

  DataObjects::Workspace2D_sptr bkgdWS =
      boost::dynamic_pointer_cast<DataObjects::Workspace2D>(
          API::WorkspaceFactory::Instance().create("Workspace2D", 1,
                                                   realx.size(), realy.size()));
  for (size_t i = 0; i < realx.size(); ++i) {
    bkgdWS->dataX(0)[i] = realx[i];
    bkgdWS->dataY(0)[i] = realy[i];
    bkgdWS->dataE(0)[i] = reale[i];
  }

  // Select background points according to mode
  if (mode.compare("All Background Points") == 0) {
    // Select (possibly) all background points
    m_outputWS = autoBackgroundSelection(bkgdWS);
  } else if (mode.compare("Input Background Points Only") == 0) {
    // Use the input background points only
    m_outputWS = bkgdWS;
  } else {
    stringstream errss;
    errss << "Background select mode " << mode
          << " is not supported by ProcessBackground.";
    g_log.error(errss.str());
    throw runtime_error(errss.str());
  }

  return;
}

    /**
    * 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);
          }
        }
//.........这里部分代码省略.........

/**
* 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;
  }
}

//.........这里部分代码省略.........
      // 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->setTitle(title);
      localWorkspace->getAxis(0)->unit() = UnitFactory::Instance().create("TOF");
      // Run parameters
      helper->loadRunParameters(localWorkspace, isisRaw.get());
      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));
        }

        isisRaw->skipData(file,period*(m_numberOfSpectra+1));
        int counter = 0;
        for (int i = 1; i <= m_numberOfSpectra; ++i)
        {
            int histToRead = i + period*(m_numberOfSpectra+1);
            if ((i >= m_spec_min && i < m_spec_max) ||
                (m_list && find(m_spec_list.begin(),m_spec_list.end(),i) != m_spec_list.end()))
            {
                isisRaw->readData(file,histToRead);
                // Copy the data into the workspace vector, discarding the 1st entry, which is rubbish
                // But note that the last (overflow) bin is kept
                MantidVec& Y = localWorkspace->dataY(counter);
                Y.assign(isisRaw->dat1 + 1, isisRaw->dat1 + lengthIn);
                // Fill the vector for the errors, containing sqrt(count)
                MantidVec& E = localWorkspace->dataE(counter);
                std::transform(Y.begin(), Y.end(), E.begin(), dblSqrt);
                // Set the X vector pointer and spectrum number
                localWorkspace->setX(counter, timeChannelsVec);
                localWorkspace->getAxis(1)->setValue(counter, i);
                // NOTE: Raw numbers go straight into the workspace
                //     - no account taken of bin widths/units etc.
                ++counter;
                pr.report();
            }
            else
            {
                isisRaw->skipData(file,histToRead);
            }
        }

        // Just a sanity check
        assert(counter == total_specs);

        std::string outputWorkspace = "OutputWorkspace";
        if (period == 0)
        {
          // Only run the Child Algorithms once
          runLoadInstrument(localWorkspace );
          runLoadMappingTable(localWorkspace );
          runLoadLog(localWorkspace );
          const int period_number = period + 1;
          Property* log=createPeriodLog(period_number);
          if(log)
          {
            localWorkspace->mutableRun().addLogData(log);
            localWorkspace->mutableRun().addLogData(createCurrentPeriodLog(period_number));
          }