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

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

void CompressEvents::exec()
{
  // Get the input workspace
  EventWorkspace_sptr inputWS = getProperty("InputWorkspace");
  EventWorkspace_sptr outputWS = getProperty("OutputWorkspace");
  double tolerance = getProperty("Tolerance");

  // Some starting things
  bool inplace = (inputWS == outputWS);
  const int noSpectra = static_cast<int>(inputWS->getNumberHistograms());
  Progress prog(this,0.0,1.0, noSpectra*2);

  // Sort the input workspace in-place by TOF. This can be faster if there are few event lists.
  inputWS->sortAll(TOF_SORT, &prog);

  // Are we making a copy of the input workspace?
  if (!inplace)
  {
    //Make a brand new EventWorkspace
    outputWS = boost::dynamic_pointer_cast<EventWorkspace>(
        API::WorkspaceFactory::Instance().create("EventWorkspace", inputWS->getNumberHistograms(), 2, 1));
    //Copy geometry over.
    API::WorkspaceFactory::Instance().initializeFromParent(inputWS, outputWS, false);
    // We DONT copy the data though

    // Do we want to parallelize over event lists, or in each event list
    bool parallel_in_each = noSpectra < PARALLEL_GET_MAX_THREADS;
    //parallel_in_each = false;

    // Loop over the histograms (detector spectra)
    // Don't parallelize the loop if we are going to parallelize each event list.
    // cppcheck-suppress syntaxError
    PRAGMA_OMP( parallel for schedule(dynamic) if (!parallel_in_each) )
    for (int i = 0; i < noSpectra; ++i)
    {
      PARALLEL_START_INTERUPT_REGION
      //the loop variable i can't be signed because of OpenMp rules inforced in Linux. Using this signed type suppresses warnings below
      const size_t index = static_cast<size_t>(i);
      // The input event list
      EventList& input_el = inputWS->getEventList(index);
      // And on the output side
      EventList & output_el = outputWS->getOrAddEventList(index);
      // Copy other settings into output
      output_el.setX( input_el.ptrX() );

      // The EventList method does the work.
      input_el.compressEvents(tolerance, &output_el, parallel_in_each);

      prog.report("Compressing");
      PARALLEL_END_INTERUPT_REGION
    }
    PARALLEL_CHECK_INTERUPT_REGION

  }

    /** Executes the rebin algorithm
    *
    *  @throw runtime_error Thrown if the bin range does not intersect the range of the input workspace
    */
    void SortEvents::exec()
    {
      // Get the input workspace
      EventWorkspace_sptr eventW = getProperty("InputWorkspace");
      //And other properties
      std::string sortoption = getPropertyValue("SortBy");

      //------- EventWorkspace ---------------------------
      const size_t histnumber = eventW->getNumberHistograms();

      //Initialize progress reporting.
      Progress prog(this,0.0,1.0, histnumber);

      DataObjects::EventSortType sortType = DataObjects::TOF_SORT;
      if (sortoption == "Pulse Time")
        sortType = DataObjects::PULSETIME_SORT;
      else if (sortoption == "Pulse Time + TOF")
        sortType = DataObjects::PULSETIMETOF_SORT;

      //This runs the SortEvents algorithm in parallel
      eventW->sortAll(sortType, &prog);

      return;
    }

/** Carries out the bin-by-bin normalisation
 *  @param inputWorkspace The input workspace
 *  @param outputWorkspace The result workspace
 */
void NormaliseToMonitor::normaliseBinByBin(API::MatrixWorkspace_sptr inputWorkspace,
                                           API::MatrixWorkspace_sptr& outputWorkspace)
{ 
  EventWorkspace_sptr inputEvent = boost::dynamic_pointer_cast<EventWorkspace>(inputWorkspace);
  EventWorkspace_sptr outputEvent;

  // Only create output workspace if different to input one
  if (outputWorkspace != inputWorkspace )
  {
    if (inputEvent)
    {
      //Make a brand new EventWorkspace
      outputEvent = boost::dynamic_pointer_cast<EventWorkspace>(
          API::WorkspaceFactory::Instance().create("EventWorkspace", inputEvent->getNumberHistograms(), 2, 1));
      //Copy geometry and data
      API::WorkspaceFactory::Instance().initializeFromParent(inputEvent, outputEvent, false);
      outputEvent->copyDataFrom( (*inputEvent) );
      outputWorkspace = boost::dynamic_pointer_cast<MatrixWorkspace>(outputEvent);
    }
    else
      outputWorkspace = WorkspaceFactory::Instance().create(inputWorkspace);
  }

  // Get hold of the monitor spectrum
  const MantidVec& monX = m_monitor->readX(0);
  MantidVec& monY = m_monitor->dataY(0);
  MantidVec& monE = m_monitor->dataE(0);
  // Calculate the overall normalisation just the once if bins are all matching
  if (m_commonBins) this->normalisationFactor(m_monitor->readX(0),&monY,&monE);


  const size_t numHists = inputWorkspace->getNumberHistograms();
  MantidVec::size_type specLength = inputWorkspace->blocksize();
  Progress prog(this,0.0,1.0,numHists);
  // Loop over spectra
  PARALLEL_FOR3(inputWorkspace,outputWorkspace,m_monitor)
  for (int64_t i = 0; i < int64_t(numHists); ++i)
  {
    PARALLEL_START_INTERUPT_REGION
    prog.report();

    const MantidVec& X = inputWorkspace->readX(i);
    // If not rebinning, just point to our monitor spectra, otherwise create new vectors
    MantidVec* Y = ( m_commonBins ? &monY : new MantidVec(specLength) );
    MantidVec* E = ( m_commonBins ? &monE : new MantidVec(specLength) );

    if (!m_commonBins)
    {
      // ConvertUnits can give X vectors of all zeroes - skip these, they cause problems
      if (X.back() == 0.0 && X.front() == 0.0) continue;
      // Rebin the monitor spectrum to match the binning of the current data spectrum
      VectorHelper::rebinHistogram(monX,monY,monE,X,*Y,*E,false);
      // Recalculate the overall normalisation factor
      this->normalisationFactor(X,Y,E);
    }

    if (inputEvent)
    {
      // ----------------------------------- EventWorkspace ---------------------------------------
      EventList & outEL = outputEvent->getEventList(i);
      outEL.divide(X, *Y, *E);
    }
    else
    {
      // ----------------------------------- Workspace2D ---------------------------------------
      const MantidVec& inY = inputWorkspace->readY(i);
      const MantidVec& inE = inputWorkspace->readE(i);
      MantidVec& YOut = outputWorkspace->dataY(i);
      MantidVec& EOut = outputWorkspace->dataE(i);
      outputWorkspace->dataX(i) = inputWorkspace->readX(i);
      // The code below comes more or less straight out of Divide.cpp
      for (MantidVec::size_type k = 0; k < specLength; ++k)
      {
        // Get references to the input Y's
        const double& leftY = inY[k];
        const double& rightY = (*Y)[k];

        // Calculate result and store in local variable to avoid overwriting original data if
        // output workspace is same as one of the input ones
        const double newY = leftY/rightY;

        if (fabs(rightY)>1.0e-12 && fabs(newY)>1.0e-12)
        {
          const double lhsFactor = (inE[k]<1.0e-12|| fabs(leftY)<1.0e-12) ? 0.0 : pow((inE[k]/leftY),2);
          const double rhsFactor = (*E)[k]<1.0e-12 ? 0.0 : pow(((*E)[k]/rightY),2);
          EOut[k] = std::abs(newY) * sqrt(lhsFactor+rhsFactor);
        }

        // Now store the result
        YOut[k] = newY;
      } // end Workspace2D case
    } // end loop over current spectrum

    if (!m_commonBins) { delete Y; delete E; }
    PARALLEL_END_INTERUPT_REGION
  } // end loop over spectra
//.........这里部分代码省略.........

void SANSSolidAngleCorrection::execEvent() {
  MatrixWorkspace_sptr inputWS = getProperty("InputWorkspace");
  EventWorkspace_sptr inputEventWS =
      boost::dynamic_pointer_cast<EventWorkspace>(inputWS);

  const int numberOfSpectra =
      static_cast<int>(inputEventWS->getNumberHistograms());
  Progress progress(this, 0.0, 1.0, inputEventWS->getNumberHistograms());

  // generate the output workspace pointer
  MatrixWorkspace_sptr outputWS = this->getProperty("OutputWorkspace");
  EventWorkspace_sptr outputEventWS;
  if (outputWS == inputWS)
    outputEventWS = boost::dynamic_pointer_cast<EventWorkspace>(outputWS);
  else {
    // Make a brand new EventWorkspace
    outputEventWS = boost::dynamic_pointer_cast<EventWorkspace>(
        WorkspaceFactory::Instance().create(
            "EventWorkspace", inputEventWS->getNumberHistograms(), 2, 1));
    // Copy geometry over.
    WorkspaceFactory::Instance().initializeFromParent(inputEventWS,
                                                      outputEventWS, false);
    // You need to copy over the data as well.
    outputEventWS->copyDataFrom((*inputEventWS));

    // Cast to the matrixOutputWS and save it
    outputWS = boost::dynamic_pointer_cast<MatrixWorkspace>(outputEventWS);
    this->setProperty("OutputWorkspace", outputWS);
  }

  progress.report("Solid Angle Correction");

  PARALLEL_FOR2(inputEventWS, outputEventWS)
  for (int i = 0; i < numberOfSpectra; i++) {
    PARALLEL_START_INTERUPT_REGION
    IDetector_const_sptr det;
    try {
      det = inputEventWS->getDetector(i);
    } catch (Exception::NotFoundError &) {
      g_log.warning() << "Spectrum index " << i
                      << " has no detector assigned to it - discarding"
                      << std::endl;
      // Catch if no detector. Next line tests whether this happened - test
      // placed
      // outside here because Mac Intel compiler doesn't like 'continue' in a
      // catch
      // in an openmp block.
    }
    if (!det)
      continue;

    // Skip if we have a monitor or if the detector is masked.
    if (det->isMonitor() || det->isMasked())
      continue;

    // Compute solid angle correction factor
    const bool is_tube = getProperty("DetectorTubes");
    const double tanTheta = tan(inputEventWS->detectorTwoTheta(det));
    const double theta_term = sqrt(tanTheta * tanTheta + 1.0);
    double corr;
    if (is_tube) {
      const double tanAlpha = tan(getYTubeAngle(det, inputWS));
      const double alpha_term = sqrt(tanAlpha * tanAlpha + 1.0);
      corr = alpha_term * theta_term * theta_term;
    } else {
      corr = theta_term * theta_term * theta_term;
    }
    EventList &el = outputEventWS->getEventList(i);
    el *= corr;
    progress.report("Solid Angle Correction");
    PARALLEL_END_INTERUPT_REGION
  }
  PARALLEL_CHECK_INTERUPT_REGION

  setProperty("OutputMessage", "Solid angle correction applied");
}

//.........这里部分代码省略.........
    bool   specify_size      = getProperty( "SpecifySize" );
    double peak_radius       = getProperty( "PeakSize" );
    double back_inner_radius = getProperty( "BackgroundInnerSize" );
    double back_outer_radius = getProperty( "BackgroundOuterSize" );
    if ( specify_size )
    {
      if ( back_outer_radius > radius )
        throw std::runtime_error("BackgroundOuterSize must be less than or equal to the RegionRadius");

      if ( back_inner_radius >= back_outer_radius )
        throw std::runtime_error("BackgroundInnerSize must be less BackgroundOuterSize");

      if ( peak_radius > back_inner_radius )
        throw std::runtime_error("PeakSize must be less than or equal to the BackgroundInnerSize");
    }
                    
                                                  // make the integrator
    Integrate3DEvents integrator( peak_q_list, UBinv, radius );

                                                  // get the events and add
                                                  // them to the inegrator
    // set up a descripter of where we are going
    this->initTargetWSDescr(wksp);

    // units conersion helper
    UnitsConversionHelper unitConv;
    unitConv.initialize(m_targWSDescr, "Momentum");

    // initialize the MD coordinates conversion class
    MDTransf_sptr q_converter = MDTransfFactory::Instance().create(m_targWSDescr.AlgID);
    q_converter->initialize(m_targWSDescr);

    // set up the progress bar
    const size_t numSpectra = wksp->getNumberHistograms();
    Progress prog(this, 0.5, 1.0, numSpectra);

    // loop through the eventlists
    std::vector<double> buffer(DIMS);

    std::vector<V3D>  event_qs;
    for (std::size_t i = 0; i < numSpectra; ++i)
    {
      // get a reference to the event list
      const EventList& events = wksp->getEventList(i);

      // check to see if the event list is empty
      if (events.empty())
      {
        prog.report();
        continue; // nothing to do
      }

      // update which pixel is being converted
      std::vector<coord_t>locCoord(DIMS, 0.);
      unitConv.updateConversion(i);
      q_converter->calcYDepCoordinates(locCoord, i);

      // loop over the events
      double signal(1.);  // ignorable garbage
      double errorSq(1.); // ignorable garbage
      const std::vector<TofEvent>& raw_events = events.getEvents();
      event_qs.clear();
      for (auto event = raw_events.begin(); event != raw_events.end(); ++event)
      {
        double val = unitConv.convertUnits( event->tof() );
        q_converter->calcMatrixCoord( val, locCoord, signal, errorSq );