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

/** Converts an EventWorkspace to an equivalent Workspace2D
 * @param inputMatrixW :: input event workspace
 * @return a MatrixWorkspace_sptr
 */
MatrixWorkspace_sptr
EventWorkspaceHelpers::convertEventTo2D(MatrixWorkspace_sptr inputMatrixW) {
  EventWorkspace_sptr inputW =
      boost::dynamic_pointer_cast<EventWorkspace>(inputMatrixW);
  if (!inputW)
    throw std::invalid_argument("EventWorkspaceHelpers::convertEventTo2D(): "
                                "Input workspace is not an EventWorkspace.");

  size_t numBins = inputW->blocksize();

  // Make a workspace 2D version of it
  MatrixWorkspace_sptr outputW;
  outputW = WorkspaceFactory::Instance().create(
      "Workspace2D", inputW->getNumberHistograms(), numBins + 1, numBins);
  WorkspaceFactory::Instance().initializeFromParent(inputW, outputW, false);

  // Now let's set all the X bins and values
  for (size_t i = 0; i < inputW->getNumberHistograms(); i++) {
    outputW->getSpectrum(i).copyInfoFrom(inputW->getSpectrum(i));
    outputW->setX(i, inputW->refX(i));

    MantidVec &Yout = outputW->dataY(i);
    const MantidVec &Yin = inputW->readY(i);
    for (size_t j = 0; j < numBins; j++)
      Yout[j] = Yin[j];

    MantidVec &Eout = outputW->dataE(i);
    const MantidVec &Ein = inputW->readE(i);
    for (size_t j = 0; j < numBins; j++)
      Eout[j] = Ein[j];
  }

  return outputW;
}

/** Execute the algorithm for a EventWorkspace input
 * @param ws :: EventWorkspace
 */
void SmoothNeighbours::execEvent(Mantid::DataObjects::EventWorkspace_sptr ws) {
  m_progress->resetNumSteps(inWS->getNumberHistograms(), 0.5, 1.0);

  // Get some stuff from the input workspace
  const size_t numberOfSpectra = outWI;
  const int YLength = static_cast<int>(inWS->blocksize());

  EventWorkspace_sptr outWS;
  // Make a brand new EventWorkspace
  outWS = boost::dynamic_pointer_cast<EventWorkspace>(
      API::WorkspaceFactory::Instance().create(
          "EventWorkspace", numberOfSpectra, YLength + 1, YLength));
  // Copy geometry over.
  API::WorkspaceFactory::Instance().initializeFromParent(*ws, *outWS, false);
  // Ensure thread-safety
  outWS->sortAll(TOF_SORT, nullptr);

  this->setProperty("OutputWorkspace",
                    boost::dynamic_pointer_cast<MatrixWorkspace>(outWS));

  // Go through all the output workspace
  PARALLEL_FOR_IF(Kernel::threadSafe(*ws, *outWS))
  for (int outWIi = 0; outWIi < int(numberOfSpectra); outWIi++) {
    PARALLEL_START_INTERUPT_REGION

    // Create the output event list (empty)
    EventList &outEL = outWS->getSpectrum(outWIi);

    // Which are the neighbours?
    std::vector<weightedNeighbour> &neighbours = m_neighbours[outWIi];
    std::vector<weightedNeighbour>::iterator it;
    for (it = neighbours.begin(); it != neighbours.end(); ++it) {
      size_t inWI = it->first;
      // if(sum)outEL.copyInfoFrom(*ws->getSpectrum(inWI));
      double weight = it->second;
      // Copy the event list
      EventList tmpEL = ws->getSpectrum(inWI);
      // Scale it
      tmpEL *= weight;
      // Add it
      outEL += tmpEL;
    }

    // Copy the single detector ID (of the center) and spectrum number from the
    // input workspace
    // if (!sum) outEL.copyInfoFrom(*ws->getSpectrum(outWIi));

    m_progress->report("Summing");
    PARALLEL_END_INTERUPT_REGION
  }
  PARALLEL_CHECK_INTERUPT_REGION

  // Give the 0-th X bins to all the output spectra.
  outWS->setAllX(inWS->binEdges(0));
  if (expandSumAllPixels)
    spreadPixels(outWS);
}

void GatherWorkspaces::execEvent() {

  // Every process in an MPI job must hit this next line or everything hangs!
  mpi::communicator included; // The communicator containing all processes
  // The root process needs to create a workspace of the appropriate size
  EventWorkspace_sptr outputWorkspace;
  if (included.rank() == 0) {
    g_log.debug() << "Total number of spectra is " << totalSpec << "\n";
    // Create the workspace for the output
    outputWorkspace = boost::dynamic_pointer_cast<EventWorkspace>(
        API::WorkspaceFactory::Instance().create("EventWorkspace", sumSpec,
                                                 numBins + hist, numBins));
    // Copy geometry over.
    API::WorkspaceFactory::Instance().initializeFromParent(
        eventW, outputWorkspace, true);
    setProperty("OutputWorkspace", outputWorkspace);
    ExperimentInfo_sptr inWS = inputWorkspace;
    outputWorkspace->copyExperimentInfoFrom(inWS.get());
  }

  for (size_t wi = 0; wi < totalSpec; wi++) {
    if (included.rank() == 0) {
      // How do we accumulate the data?
      std::string accum = this->getPropertyValue("AccumulationMethod");
      std::vector<Mantid::DataObjects::EventList> out_values;
      gather(included, eventW->getEventList(wi), out_values, 0);
      for (int i = 0; i < included.size(); i++) {
        size_t index = wi; // accum == "Add"
        if (accum == "Append")
          index = wi + i * totalSpec;
        outputWorkspace->dataX(index) = eventW->readX(wi);
        outputWorkspace->getOrAddEventList(index) += out_values[i];
        const ISpectrum *inSpec = eventW->getSpectrum(wi);
        ISpectrum *outSpec = outputWorkspace->getSpectrum(index);
        outSpec->clearDetectorIDs();
        outSpec->addDetectorIDs(inSpec->getDetectorIDs());
      }
    } else {
      gather(included, eventW->getEventList(wi), 0);
    }
  }
}

//.........这里部分代码省略.........
  partWorkspaces.resize(numThreads);
  buffers.resize(numThreads);
  eventVectors = new EventVector_pt *[numThreads];

  // cppcheck-suppress syntaxError
  PRAGMA_OMP( parallel for if (parallelProcessing) )
  for (int i = 0; i < int(numThreads); i++) {
    // This is the partial workspace we are about to create (if in parallel)
    EventWorkspace_sptr partWS;
    if (parallelProcessing) {
      prog->report("Creating Partial Workspace");
      // Create a partial workspace
      partWS = EventWorkspace_sptr(new EventWorkspace());
      // Make sure to initialize.
      partWS->initialize(1, 1, 1);
      // Copy all the spectra numbers and stuff (no actual events to copy
      // though).
      partWS->copyDataFrom(*workspace);
      // Push it in the array
      partWorkspaces[i] = partWS;
    } else
      partWS = workspace;

    // Allocate the buffers
    buffers[i] = new DasEvent[loadBlockSize];

    // For each partial workspace, make an array where index = detector ID and
    // value = pointer to the events vector
    eventVectors[i] = new EventVector_pt[detid_max + 1];
    EventVector_pt *theseEventVectors = eventVectors[i];
    for (detid_t j = 0; j < detid_max + 1; j++) {
      size_t wi = pixel_to_wkspindex[j];
      // Save a POINTER to the vector<tofEvent>
      theseEventVectors[j] = &partWS->getSpectrum(wi).getEvents();
    }
  }

  g_log.debug() << tim << " to create " << partWorkspaces.size()
                << " workspaces for parallel loading.\n";

  prog->resetNumSteps(numBlocks, 0.1, 0.8);

  // ---------------------------------- LOAD THE DATA --------------------------
  PRAGMA_OMP( parallel for schedule(dynamic, 1) if (parallelProcessing) )
  for (int blockNum = 0; blockNum < int(numBlocks); blockNum++) {
    PARALLEL_START_INTERUPT_REGION

    // Find the workspace for this particular thread
    EventWorkspace_sptr ws;
    size_t threadNum = 0;
    if (parallelProcessing) {
      threadNum = PARALLEL_THREAD_NUMBER;
      ws = partWorkspaces[threadNum];
    } else
      ws = workspace;

    // Get the buffer (for this thread)
    DasEvent *event_buffer = buffers[threadNum];

    // Get the speeding-up array of vector<tofEvent> where index = detid.
    EventVector_pt *theseEventVectors = eventVectors[threadNum];

    // Where to start in the file?
    size_t fileOffset = first_event + (loadBlockSize * blockNum);
    // May need to reduce size of last (or only) block
    size_t current_event_buffer_size =

/** Execute the algorithm.
 */
void ResampleX::exec() {
  // generically having access to the input workspace is a good idea
  MatrixWorkspace_sptr inputWS = getProperty("InputWorkspace");
  MatrixWorkspace_sptr outputWS = getProperty("OutputWorkspace");
  bool inPlace = (inputWS == outputWS); // Rebinning in-place
  m_isDistribution = inputWS->isDistribution();
  m_isHistogram = inputWS->isHistogramData();
  int numSpectra = static_cast<int>(inputWS->getNumberHistograms());

  // the easy parameters
  m_useLogBinning = getProperty("LogBinning");
  m_numBins = getProperty("NumberBins");
  m_preserveEvents = getProperty("PreserveEvents");

  // determine the xmin/xmax for the workspace
  vector<double> xmins = getProperty("XMin");
  vector<double> xmaxs = getProperty("XMax");
  string error = determineXMinMax(inputWS, xmins, xmaxs);
  if (!error.empty())
    throw std::runtime_error(error);

  bool common_limits = true;
  {
    double xmin_common = xmins[0];
    double xmax_common = xmaxs[0];
    for (size_t i = 1; i < xmins.size(); ++i) {
      if (xmins[i] != xmin_common) {
        common_limits = false;
        break;
      }
      if (xmaxs[i] != xmax_common) {
        common_limits = false;
        break;
      }
    }
  }
  if (common_limits) {
    g_log.debug() << "Common limits between all spectra\n";
  } else {
    g_log.debug() << "Does not have common limits between all spectra\n";
  }

  // start doing actual work
  EventWorkspace_const_sptr inputEventWS =
      boost::dynamic_pointer_cast<const EventWorkspace>(inputWS);
  if (inputEventWS != NULL) {
    if (m_preserveEvents) {
      EventWorkspace_sptr outputEventWS =
          boost::dynamic_pointer_cast<EventWorkspace>(outputWS);
      if (inPlace) {
        g_log.debug() << "Rebinning event workspace in place\n";
      } else {
        g_log.debug() << "Rebinning event workspace out of place\n";

        // copy the event workspace to a new EventWorkspace
        outputEventWS = boost::dynamic_pointer_cast<EventWorkspace>(
            API::WorkspaceFactory::Instance().create(
                "EventWorkspace", inputWS->getNumberHistograms(), 2, 1));
        // copy geometry over.
        API::WorkspaceFactory::Instance().initializeFromParent(
            inputEventWS, outputEventWS, false);
        // copy over the data as well.
        outputEventWS->copyDataFrom((*inputEventWS));
      }

      if (common_limits) {
        // get the delta from the first since they are all the same
        MantidVecPtr xValues;
        double delta =
            this->determineBinning(xValues.access(), xmins[0], xmaxs[0]);
        g_log.debug() << "delta = " << delta << "\n";
        outputEventWS->setAllX(xValues);
      } else {
        // initialize progress reporting.
        Progress prog(this, 0.0, 1.0, numSpectra);

        // do the rebinning
        PARALLEL_FOR2(inputEventWS, outputWS)
        for (int wkspIndex = 0; wkspIndex < numSpectra; ++wkspIndex) {
          PARALLEL_START_INTERUPT_REGION
          MantidVec xValues;
          double delta = this->determineBinning(xValues, xmins[wkspIndex],
                                                xmaxs[wkspIndex]);
          g_log.debug() << "delta[wkspindex=" << wkspIndex << "] = " << delta
                        << " xmin=" << xmins[wkspIndex]
                        << " xmax=" << xmaxs[wkspIndex] << "\n";
          outputEventWS->getSpectrum(wkspIndex)->setX(xValues);
          prog.report(name()); // Report progress
          PARALLEL_END_INTERUPT_REGION
        }
        PARALLEL_CHECK_INTERUPT_REGION
      }

      this->setProperty(
          "OutputWorkspace",
          boost::dynamic_pointer_cast<MatrixWorkspace>(outputEventWS));
    }    // end if (m_preserveEvents)
    else // event workspace -> matrix workspace
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  // Loop over the histograms (detector spectra)
  // PARALLEL_FOR_IF(Kernel::threadSafe(*outputWS))
  for (int64_t i = 0; i < numberOfSpectra_i; ++i) {

    // Lets find what row this spectrum Number appears in our detector table.

    // PARALLEL_START_INTERUPT_REGION

    std::size_t wsid = i;

    try {

      double deg2rad = M_PI / 180.;

      auto det = outputWS->getDetector(i);
      int specNo = det->getID();

      // int spectraNumber = static_cast<int>(spectraColumn->toDouble(i));
      // wsid = outputWS->getIndexFromSpectrumNumber(spectraNumber);
      g_log.debug() << "###### Spectra #" << specNo
                    << " ==> Workspace ID:" << wsid << '\n';

      // Now we need to find the row that contains this spectrum
      std::vector<int>::const_iterator specIter;

      specIter = std::find(spectraColumn.begin(), spectraColumn.end(), specNo);
      if (specIter != spectraColumn.end()) {
        const size_t detectorRow =
            std::distance(spectraColumn.begin(), specIter);
        const double l1 = l1Column[detectorRow];
        const double l2 = l2Column[detectorRow];
        const double twoTheta = twoThetaColumn[detectorRow] * deg2rad;
        const double efixed = efixedColumn[detectorRow];
        const int emode = emodeColumn[detectorRow];

        if (g_log.is(Logger::Priority::PRIO_DEBUG)) {
          g_log.debug() << "specNo from detector table = "
                        << spectraColumn[detectorRow] << '\n';

          g_log.debug() << "###### Spectra #" << specNo
                        << " ==> Det Table Row:" << detectorRow << '\n';

          g_log.debug() << "\tL1=" << l1 << ",L2=" << l2 << ",TT=" << twoTheta
                        << ",EF=" << efixed << ",EM=" << emode << '\n';
        }

        // Make local copies of the units. This allows running the loop in
        // parallel
        auto localFromUnit = std::unique_ptr<Unit>(fromUnit->clone());
        auto localOutputUnit = std::unique_ptr<Unit>(outputUnit->clone());
        /// @todo Don't yet consider hold-off (delta)
        const double delta = 0.0;
        std::vector<double> values(outputWS->x(wsid).begin(),
                                   outputWS->x(wsid).end());

        // Convert the input unit to time-of-flight
        localFromUnit->toTOF(values, emptyVec, l1, l2, twoTheta, emode, efixed,
                             delta);
        // Convert from time-of-flight to the desired unit
        localOutputUnit->fromTOF(values, emptyVec, l1, l2, twoTheta, emode,
                                 efixed, delta);

        outputWS->mutableX(wsid) = std::move(values);

        // EventWorkspace part, modifying the EventLists.
        if (m_inputEvents) {
          eventWS->getSpectrum(wsid)
              .convertUnitsViaTof(localFromUnit.get(), localOutputUnit.get());
        }

      } else {
        // Not found
        failedDetectorCount++;
        outputWS->maskWorkspaceIndex(wsid);
      }

    } catch (Exception::NotFoundError &) {
      // Get to here if exception thrown when calculating distance to detector
      failedDetectorCount++;
      // Since you usually (always?) get to here when there's no attached
      // detectors, this call is
      // the same as just zeroing out the data (calling clearData on the
      // spectrum)
      outputWS->maskWorkspaceIndex(i);
    }

    prog.report("Convert to " + m_outputUnit->unitID());
    // PARALLEL_END_INTERUPT_REGION
  } // loop over spectra
  // PARALLEL_CHECK_INTERUPT_REGION

  if (failedDetectorCount != 0) {
    g_log.information() << "Something went wrong for " << failedDetectorCount
                        << " spectra. Masking spectrum.\n";
  }
  if (m_inputEvents)
    eventWS->clearMRU();

  return outputWS;
}