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

void ScaleX::execEvent()
{
  g_log.information("Processing event workspace");
  const MatrixWorkspace_const_sptr matrixInputWS = this->getProperty("InputWorkspace");
  const std::string op = getPropertyValue("Operation");
  EventWorkspace_const_sptr inputWS=boost::dynamic_pointer_cast<const EventWorkspace>(matrixInputWS);
  // generate the output workspace pointer
  API::MatrixWorkspace_sptr matrixOutputWS = this->getProperty("OutputWorkspace");
  EventWorkspace_sptr outputWS;
  if (matrixOutputWS == matrixInputWS)
    outputWS = boost::dynamic_pointer_cast<EventWorkspace>(matrixOutputWS);
  else
  {
    //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);
    //You need to copy over the data as well.
    outputWS->copyDataFrom( (*inputWS) );
    //Cast to the matrixOutputWS and save it
    matrixOutputWS = boost::dynamic_pointer_cast<MatrixWorkspace>(outputWS);
    this->setProperty("OutputWorkspace", matrixOutputWS);
  }
  int numHistograms = static_cast<int>(inputWS->getNumberHistograms());
  PARALLEL_FOR1(outputWS)
  for (int i=0; i < numHistograms; ++i)
  {
    PARALLEL_START_INTERUPT_REGION
    //Do the offsetting
    if ((i >= m_wi_min) && (i <= m_wi_max))
    {
      if(op=="Multiply")
      {
        outputWS->getEventList(i).scaleTof(getScaleFactor(inputWS, i));
        if( m_algFactor < 0 )
        {
          outputWS->getEventList(i).reverse();
        }
      }
      else if(op=="Add")
      {
        outputWS->getEventList(i).addTof(getScaleFactor(inputWS, i));
      }
    }
    m_progress->report("Scaling X");
    PARALLEL_END_INTERUPT_REGION
  }
  PARALLEL_CHECK_INTERUPT_REGION
  outputWS->clearMRU();
}

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

  }

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

/** Executes the algorithm
 */
void FilterByTime::exec()
{
  EventWorkspace_const_sptr inputWS = this->getProperty("InputWorkspace");

  // ---- Find the start/end times ----
  DateAndTime start, stop;

  double start_dbl, stop_dbl;
  start_dbl = getProperty("StartTime");
  stop_dbl = getProperty("StopTime");

  std::string start_str, stop_str;
  start_str = getPropertyValue("AbsoluteStartTime");
  stop_str = getPropertyValue("AbsoluteStopTime");

  if ( (start_str != "") && (stop_str != "") && (start_dbl <= 0.0) && (stop_dbl <= 0.0) )
  {
    // Use the absolute string
    start = DateAndTime( start_str );
    stop  = DateAndTime( stop_str );
  }
  else if ( (start_str == "") && (stop_str == "") && ((start_dbl > 0.0) || (stop_dbl > 0.0)) )
  {
    // Use the relative times in seconds.
    DateAndTime first = inputWS->getFirstPulseTime();
    DateAndTime last = inputWS->getLastPulseTime();
    start = first + start_dbl;
    if (stop_dbl > 0.0)
    {
      stop = first + stop_dbl;
    }
    else
    {
      this->getLogger().debug() << "No end filter time specified - assuming last pulse" << std::endl;
      stop = last + 10000.0;   // so we get all events - needs to be past last pulse
    }
  }
  else
  {
    //Either both or none were specified
    throw std::invalid_argument("You need to specify either the StartTime or StopTime parameters; or both the AbsoluteStartTime and AbsoluteStopTime parameters; but not other combinations.");
  }

  if (stop <= start)
    throw std::invalid_argument("The stop time should be larger than the start time.");

  // Make a brand new EventWorkspace
  EventWorkspace_sptr 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);
  // But we don't copy the data.

  setProperty("OutputWorkspace", outputWS);

  size_t numberOfSpectra = inputWS->getNumberHistograms();

  // Initialise the progress reporting object
  Progress prog(this,0.0,1.0,numberOfSpectra);

  // Loop over the histograms (detector spectra)
  PARALLEL_FOR_NO_WSP_CHECK()
  for (int64_t i = 0; i < int64_t(numberOfSpectra); ++i)
  {
    PARALLEL_START_INTERUPT_REGION

    //Get the output event list (should be empty)
    EventList& output_el = outputWS->getEventList(i);
    //and this is the input event list
    const EventList& input_el = inputWS->getEventList(i);

    //Perform the filtering
    input_el.filterByPulseTime(start, stop, output_el);

    prog.report();
    PARALLEL_END_INTERUPT_REGION
  }
  PARALLEL_CHECK_INTERUPT_REGION


  //Now filter out the run, using the DateAndTime type.
  outputWS->mutableRun().filterByTime(start, stop);

}

//.........这里部分代码省略.........
  // TODO: Check why this parallel stuff breaks
  // Loop over the histograms (detector spectra)
  // PARALLEL_FOR1(outputWS)
  for (int64_t i = 0; i < numberOfSpectra_i; ++i) {

    // Lets find what row this spectrum ID 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 specid = det->getID();

      // int spectraNumber = static_cast<int>(spectraColumn->toDouble(i));
      // wsid = outputWS->getIndexFromSpectrumNumber(spectraNumber);
      g_log.debug() << "###### Spectra #" << specid
                    << " ==> Workspace ID:" << wsid << std::endl;

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

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

        g_log.debug() << "specId from detector table = "
                      << spectraColumn[detectorRow] << std::endl;

        // l1 = l1Column->toDouble(detectorRow);
        // l2 = l2Column->toDouble(detectorRow);
        // twoTheta = deg2rad * twoThetaColumn->toDouble(detectorRow);
        // efixed = efixedColumn->toDouble(detectorRow);
        // emode = static_cast<int>(emodeColumn->toDouble(detectorRow));

        g_log.debug() << "###### Spectra #" << specid
                      << " ==> Det Table Row:" << detectorRow << std::endl;

        g_log.debug() << "\tL1=" << l1 << ",L2=" << l2 << ",TT=" << twoTheta
                      << ",EF=" << efixed << ",EM=" << emode << std::endl;

        // Make local copies of the units. This allows running the loop in
        // parallel
        Unit *localFromUnit = fromUnit->clone();
        Unit *localOutputUnit = outputUnit->clone();
        /// @todo Don't yet consider hold-off (delta)
        const double delta = 0.0;
        // Convert the input unit to time-of-flight
        localFromUnit->toTOF(outputWS->dataX(wsid), emptyVec, l1, l2, twoTheta,
                             emode, efixed, delta);
        // Convert from time-of-flight to the desired unit
        localOutputUnit->fromTOF(outputWS->dataX(wsid), emptyVec, l1, l2,
                                 twoTheta, emode, efixed, delta);
        // EventWorkspace part, modifying the EventLists.
        if (m_inputEvents) {
          eventWS->getEventList(wsid)
              .convertUnitsViaTof(localFromUnit, localOutputUnit);
        }
        // Clear unit memory
        delete localFromUnit;
        delete localOutputUnit;

      } else {
        // Not found
        g_log.debug() << "Spectrum " << specid << " not found!" << std::endl;
        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." << std::endl;
  }
  if (m_inputEvents)
    eventWS->clearMRU();
}

void ModeratorTzero::execEvent()
{
  g_log.information("Processing event workspace");

  const MatrixWorkspace_const_sptr matrixInputWS = getProperty("InputWorkspace");
  EventWorkspace_const_sptr inputWS= boost::dynamic_pointer_cast<const EventWorkspace>(matrixInputWS);

  // generate the output workspace pointer
  const size_t numHists = static_cast<size_t>(inputWS->getNumberHistograms());
  Mantid::API::MatrixWorkspace_sptr matrixOutputWS = getProperty("OutputWorkspace");
  EventWorkspace_sptr outputWS;
  if (matrixOutputWS == matrixInputWS)
  {
    outputWS = boost::dynamic_pointer_cast<EventWorkspace>(matrixOutputWS);
  }
  else
  {
    //Make a brand new EventWorkspace
    outputWS = boost::dynamic_pointer_cast<EventWorkspace>(WorkspaceFactory::Instance().create("EventWorkspace", numHists, 2, 1));
    //Copy geometry over.
    WorkspaceFactory::Instance().initializeFromParent(inputWS, outputWS, false);
    //You need to copy over the data as well.
    outputWS->copyDataFrom( (*inputWS) );
    //Cast to the matrixOutputWS and save it
    matrixOutputWS = boost::dynamic_pointer_cast<MatrixWorkspace>(outputWS);
    setProperty("OutputWorkspace", matrixOutputWS);
  }

  //Get a pointer to the sample
  IComponent_const_sptr sample = outputWS->getInstrument()->getSample();

  // Loop over the spectra
  Progress prog(this,0.0,1.0,numHists); //report progress of algorithm
  PARALLEL_FOR1(outputWS)
  for (int i = 0; i < static_cast<int>(numHists); ++i)
  {
    PARALLEL_START_INTERUPT_REGION
    size_t wsIndex = static_cast<size_t>(i);
    EventList &evlist=outputWS->getEventList(wsIndex);
    if( evlist.getNumberEvents() > 0 ) //don't bother with empty lists
    {
      double L1=CalculateL1(matrixOutputWS, wsIndex); // distance from source to sample or monitor
      double t2=CalculateT2(matrixOutputWS, wsIndex); // time from sample to detector
      if(t2>=0) //t2 < 0 when no detector info is available
      {
        double tof, E1;
        mu::Parser parser;
        parser.DefineVar("incidentEnergy", &E1); // associate variable E1 to this parser
        parser.SetExpr(m_formula);
        E1=m_convfactor*(L1/m_t1min)*(L1/m_t1min);
        double min_t0_next=parser.Eval(); // fast neutrons are shifted by min_t0_next, irrespective of tof

        // fix the histogram bins
        MantidVec &x=evlist.dataX();
        for (MantidVec::iterator iter=x.begin(); iter!=x.end(); ++iter)
        {
          tof=*iter;
          if(tof<m_t1min+t2)
            tof-=min_t0_next;
          else
            tof-=CalculateT0(tof, L1, t2, E1, parser);
          *iter=tof;
        }

        MantidVec tofs=evlist.getTofs();
        for(unsigned int itof=0; itof<tofs.size(); itof++)
        {
          tof=tofs[itof]+0.002*(rand()%100 -50); // add a [-0.1,0.1] microsecond noise to avoid artifacts resulting from original tof data
          if(tof<m_t1min+t2)
            tof-=min_t0_next;
          else
            tof-=CalculateT0(tof, L1, t2, E1, parser);
          tofs[itof]=tof;
        }
        evlist.setTofs(tofs);
        evlist.setSortOrder(Mantid::DataObjects::EventSortType::UNSORTED);
      }
    }
    prog.report();
    PARALLEL_END_INTERUPT_REGION
  }
  PARALLEL_CHECK_INTERUPT_REGION
  outputWS->clearMRU(); // Clears the Most Recent Used lists */
} // end of void ModeratorTzero::execEvent()

/** Executes the algorithm
 *@param localworkspace :: the input workspace
 *@param indices :: set of indices to sum up
 */
void SumSpectra::execEvent(EventWorkspace_const_sptr localworkspace,
                           std::set<int> &indices) {
  // Make a brand new EventWorkspace
  EventWorkspace_sptr outputWorkspace =
      boost::dynamic_pointer_cast<EventWorkspace>(
          API::WorkspaceFactory::Instance().create("EventWorkspace", 1, 2, 1));
  // Copy geometry over.
  API::WorkspaceFactory::Instance().initializeFromParent(localworkspace,
                                                         outputWorkspace, true);

  Progress progress(this, 0, 1, indices.size());

  // Get the pointer to the output event list
  EventList &outEL = outputWorkspace->getEventList(0);
  outEL.setSpectrumNo(m_outSpecId);
  outEL.clearDetectorIDs();

  // Loop over spectra
  std::set<int>::iterator it;
  size_t numSpectra(0);
  size_t numMasked(0);
  size_t numZeros(0);
  // for (int i = m_minSpec; i <= m_maxSpec; ++i)
  for (it = indices.begin(); it != indices.end(); ++it) {
    int i = *it;
    // Don't go outside the range.
    if ((i >= m_numberOfSpectra) || (i < 0)) {
      g_log.error() << "Invalid index " << i
                    << " was specified. Sum was aborted.\n";
      break;
    }

    try {
      // Get the detector object for this spectrum
      Geometry::IDetector_const_sptr det = localworkspace->getDetector(i);
      // Skip monitors, if the property is set to do so
      if (!m_keepMonitors && det->isMonitor())
        continue;
      // Skip masked detectors
      if (det->isMasked()) {
        numMasked++;
        continue;
      }
    } catch (...) {
      // if the detector not found just carry on
    }
    numSpectra++;

    // Add the event lists with the operator
    const EventList &tOutEL = localworkspace->getEventList(i);
    if (tOutEL.empty()) {
      ++numZeros;
    }
    outEL += tOutEL;

    progress.report();
  }

  // Set all X bins on the output
  cow_ptr<MantidVec> XValues;
  XValues.access() = localworkspace->readX(0);
  outputWorkspace->setAllX(XValues);

  outputWorkspace->mutableRun().addProperty("NumAllSpectra", int(numSpectra),
                                            "", true);
  outputWorkspace->mutableRun().addProperty("NumMaskSpectra", int(numMasked),
                                            "", true);
  outputWorkspace->mutableRun().addProperty("NumZeroSpectra", int(numZeros), "",
                                            true);

  // Assign it to the output workspace property
  setProperty("OutputWorkspace",
              boost::dynamic_pointer_cast<MatrixWorkspace>(outputWorkspace));
}

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

  g_log.debug() << tim << " to create " << partWorkspaces.size()
                << " workspaces for parallel loading." << std::endl;

  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 align detectors algorithm for an event workspace.
 */
void AlignDetectors::execEvent() {
  // g_log.information("Processing event workspace");

  // the calibration information is already read in at this point

  // convert the input workspace into the event workspace we already know it is
  const MatrixWorkspace_const_sptr matrixInputWS =
      this->getProperty("InputWorkspace");
  EventWorkspace_const_sptr inputWS =
      boost::dynamic_pointer_cast<const EventWorkspace>(matrixInputWS);

  // generate the output workspace pointer
  API::MatrixWorkspace_sptr matrixOutputWS =
      this->getProperty("OutputWorkspace");
  EventWorkspace_sptr outputWS;
  if (matrixOutputWS == matrixInputWS)
    outputWS = boost::dynamic_pointer_cast<EventWorkspace>(matrixOutputWS);
  else {
    // 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);
    // You need to copy over the data as well.
    outputWS->copyDataFrom((*inputWS));

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

  // Set the final unit that our output workspace will have
  setXAxisUnits(outputWS);

  ConversionFactors converter = ConversionFactors(m_calibrationWS);

  Progress progress(this, 0.0, 1.0, m_numberOfSpectra);

  PARALLEL_FOR_NO_WSP_CHECK()
  for (int64_t i = 0; i < m_numberOfSpectra; ++i) {
    PARALLEL_START_INTERUPT_REGION

    auto toDspacing = converter.getConversionFunc(
        inputWS->getSpectrum(size_t(i))->getDetectorIDs());
    outputWS->getEventList(i).convertTof(toDspacing);

    progress.report();
    PARALLEL_END_INTERUPT_REGION
  }
  PARALLEL_CHECK_INTERUPT_REGION

  if (outputWS->getTofMin() < 0.) {
    std::stringstream msg;
    msg << "Something wrong with the calibration. Negative minimum d-spacing "
           "created. d_min = " << outputWS->getTofMin() << " d_max "
        << outputWS->getTofMax();
    g_log.warning(msg.str());
  }
  outputWS->clearMRU();
}

//.........这里部分代码省略.........
    {
      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 );
        for ( size_t dim = 0; dim < DIMS; ++dim )
        {
          buffer[dim] = locCoord[dim];
        }
        V3D q_vec( buffer[0], buffer[1], buffer[2] );
        event_qs.push_back( q_vec );
      } // end of loop over events in list

      integrator.addEvents( event_qs );

      prog.report();
    } // end of loop over spectra

    double inti;
    double sigi;
    for ( size_t i = 0; i < n_peaks; i++ )
    {
      V3D hkl( peaks[i].getH(), peaks[i].getK(), peaks[i].getL() );
      if ( Geometry::IndexingUtils::ValidIndex( hkl, 1.0 ) ) 
      {
        V3D peak_q( peaks[i].getQLabFrame() );
        integrator.ellipseIntegrateEvents( peak_q, 
          specify_size, peak_radius, back_inner_radius, back_outer_radius,
          inti, sigi );
        peaks[i].setIntensity( inti );
        peaks[i].setSigmaIntensity( sigi );
      }
      else
      {
        peaks[i].setIntensity( 0.0 );
        peaks[i].setSigmaIntensity( 0.0 );
      }
  
    }

    peak_ws->mutableRun().addProperty("PeaksIntegrated", 1, true);

    setProperty("OutputWorkspace", peak_ws);
  }

/**
 * Execute the align detectors algorithm for an event workspace.
 */
void AlignDetectors::execEvent()
{
    //g_log.information("Processing event workspace");

    // the calibration information is already read in at this point

    // convert the input workspace into the event workspace we already know it is
    const MatrixWorkspace_const_sptr matrixInputWS = this->getProperty("InputWorkspace");
    EventWorkspace_const_sptr inputWS
        = boost::dynamic_pointer_cast<const EventWorkspace>(matrixInputWS);

    // generate the output workspace pointer
    API::MatrixWorkspace_sptr matrixOutputWS = this->getProperty("OutputWorkspace");
    EventWorkspace_sptr outputWS;
    if (matrixOutputWS == matrixInputWS)
        outputWS = boost::dynamic_pointer_cast<EventWorkspace>(matrixOutputWS);
    else
    {
        //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);
        //outputWS->mutableSpectraMap().clear();
        //You need to copy over the data as well.
        outputWS->copyDataFrom( (*inputWS) );

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

    // Set the final unit that our output workspace will have
    outputWS->getAxis(0)->unit() = UnitFactory::Instance().create("dSpacing");

    const int64_t numberOfSpectra = static_cast<int64_t>(inputWS->getNumberHistograms());

    // Initialise the progress reporting object
    Progress progress(this,0.0,1.0,numberOfSpectra);

    PARALLEL_FOR_NO_WSP_CHECK()
    for (int64_t i = 0; i < int64_t(numberOfSpectra); ++i)
    {
        PARALLEL_START_INTERUPT_REGION
        // Compute the conversion factor
        double factor = calcConversionFromMap(this->tofToDmap, inputWS->getSpectrum(size_t(i))->getDetectorIDs());

        //Perform the multiplication on all events
        outputWS->getEventList(i).convertTof(factor);

        progress.report();
        PARALLEL_END_INTERUPT_REGION
    }
    PARALLEL_CHECK_INTERUPT_REGION

    if (outputWS->getTofMin() < 0.)
    {
        std::stringstream msg;
        msg << "Something wrong with the calibration. Negative minimum d-spacing created. d_min = "
            <<  outputWS->getTofMin() << " d_max " << outputWS->getTofMax();
        throw std::runtime_error(msg.str());
    }
    outputWS->clearMRU();
}

void ModeratorTzero::execEvent(const std::string &emode) {
  g_log.information("Processing event workspace");

  const MatrixWorkspace_const_sptr matrixInputWS =
      getProperty("InputWorkspace");
  EventWorkspace_const_sptr inputWS =
      boost::dynamic_pointer_cast<const EventWorkspace>(matrixInputWS);

  // generate the output workspace pointer
  const size_t numHists = static_cast<size_t>(inputWS->getNumberHistograms());
  Mantid::API::MatrixWorkspace_sptr matrixOutputWS =
      getProperty("OutputWorkspace");
  EventWorkspace_sptr outputWS;
  if (matrixOutputWS == matrixInputWS) {
    outputWS = boost::dynamic_pointer_cast<EventWorkspace>(matrixOutputWS);
  } else {
    // Make a brand new EventWorkspace
    outputWS = boost::dynamic_pointer_cast<EventWorkspace>(
        WorkspaceFactory::Instance().create("EventWorkspace", numHists, 2, 1));
    // Copy geometry over.
    WorkspaceFactory::Instance().initializeFromParent(inputWS, outputWS, false);
    // You need to copy over the data as well.
    outputWS->copyDataFrom((*inputWS));
    // Cast to the matrixOutputWS and save it
    matrixOutputWS = boost::dynamic_pointer_cast<MatrixWorkspace>(outputWS);
    setProperty("OutputWorkspace", matrixOutputWS);
  }

  // Get pointers to sample and source
  IComponent_const_sptr source = m_instrument->getSource();
  IComponent_const_sptr sample = m_instrument->getSample();
  double Lss = source->getDistance(*sample); // distance from source to sample

  // calculate tof shift once for all neutrons if emode==Direct
  double t0_direct(-1);
  if (emode == "Direct") {
    Kernel::Property *eiprop = inputWS->run().getProperty("Ei");
    double Ei = boost::lexical_cast<double>(eiprop->value());
    mu::Parser parser;
    parser.DefineVar("incidentEnergy", &Ei); // associate E1 to this parser
    parser.SetExpr(m_formula);
    t0_direct = parser.Eval();
  }

  // Loop over the spectra
  Progress prog(this, 0.0, 1.0, numHists); // report progress of algorithm
  PARALLEL_FOR1(outputWS)
  for (int i = 0; i < static_cast<int>(numHists); ++i) {
    PARALLEL_START_INTERUPT_REGION
    size_t wsIndex = static_cast<size_t>(i);
    EventList &evlist = outputWS->getEventList(wsIndex);
    if (evlist.getNumberEvents() > 0) // don't bother with empty lists
    {
      IDetector_const_sptr det;
      double L1(Lss); // distance from source to sample
      double L2(-1);  // distance from sample to detector

      try {
        det = inputWS->getDetector(i);
        if (det->isMonitor()) {
          // redefine the sample as the monitor
          L1 = source->getDistance(*det);
          L2 = 0;
        } else {
          L2 = sample->getDistance(*det);
        }
      } catch (Exception::NotFoundError &) {
        g_log.error() << "Unable to calculate distances to/from detector" << i
                      << std::endl;
      }

      if (L2 >= 0) {
        // One parser for each parallel processor needed (except Edirect mode)
        double E1;
        mu::Parser parser;
        parser.DefineVar("incidentEnergy", &E1); // associate E1 to this parser
        parser.SetExpr(m_formula);

        // fast neutrons are shifted by min_t0_next, irrespective of tof
        double v1_max = L1 / m_t1min;
        E1 = m_convfactor * v1_max * v1_max;
        double min_t0_next = parser.Eval();

        if (emode == "Indirect") {
          double t2(-1.0); // time from sample to detector. (-1) signals error
          if (det->isMonitor()) {
            t2 = 0.0;
          } else {
            static const double convFact =
                1.0e-6 * sqrt(2 * PhysicalConstants::meV /
                              PhysicalConstants::NeutronMass);
            std::vector<double> wsProp = det->getNumberParameter("Efixed");
            if (!wsProp.empty()) {
              double E2 = wsProp.at(0);        //[E2]=meV
              double v2 = convFact * sqrt(E2); //[v2]=meter/microsec
              t2 = L2 / v2;
            } else {
              // t2 is kept to -1 if no Efixed is found
              g_log.debug() << "Efixed not found for detector " << i
                            << std::endl;
//.........这里部分代码省略.........