C++ EventWorkspace_const_sptr类代码示例

/** Execute the algorithm.
 *
 *  @throw runtime_error Thrown if algorithm cannot execute
 */
void SaveToSNSHistogramNexus::exec() {
  // NXMSetError(NULL, nexus_print_error);
  NXMEnableErrorReporting();

  // Retrieve the filename from the properties
  m_inputFilename = getPropertyValue("InputFileName");
  m_outputFilename = getPropertyValue("OutputFileName");
  m_compress = getProperty("Compress");

  inputWorkspace = getProperty("InputWorkspace");

  // We'll need to get workspace indices
  map = inputWorkspace->getDetectorIDToWorkspaceIndexMap();

  // Start the progress bar. 3 reports per histogram.
  prog = new Progress(this, 0, 1.0, inputWorkspace->getNumberHistograms() * 3);

  EventWorkspace_const_sptr eventWorkspace =
      boost::dynamic_pointer_cast<const EventWorkspace>(inputWorkspace);
  if (eventWorkspace) {
    eventWorkspace->sortAll(TOF_SORT, prog);
  }

  int ret;
  ret = this->copy_file(m_inputFilename.c_str(), NXACC_READ,
                        m_outputFilename.c_str(), NXACC_CREATE5);

  if (ret == NX_ERROR)
    throw std::runtime_error("Nexus error while copying the file.");
}

/** 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) {
  auto outputWorkspace = create<EventWorkspace>(*localworkspace, 1);

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

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

  const auto &spectrumInfo = localworkspace->spectrumInfo();
  // Loop over spectra
  size_t numSpectra(0);
  size_t numMasked(0);
  size_t numZeros(0);
  for (const auto i : indices) {
    // 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;
    }

    if (spectrumInfo.hasDetectors(i)) {
      // Skip monitors, if the property is set to do so
      if (!m_keepMonitors && spectrumInfo.isMonitor(i))
        continue;
      // Skip masked detectors
      if (spectrumInfo.isMasked(i)) {
        numMasked++;
        continue;
      }
    }
    numSpectra++;

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

    progress.report();
  }

  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", std::move(outputWorkspace));
}

void FindCenterOfMassPosition2::exec()
{
  MatrixWorkspace_sptr inputWSWvl = getProperty("InputWorkspace");
  MatrixWorkspace_sptr inputWS;

  // Option to exclude beam area
  bool direct_beam = getProperty("DirectBeam");

  //TODO: Need an input for the X bin to use, assume 0 for now
  int specID = 0;
  // Initial center location
  double center_x = getProperty("CenterX");
  double center_y = getProperty("CenterY");
  const double tolerance = getProperty("Tolerance");
  // Iteration cutoff
  int max_iteration = 200;
  // Radius of the beam area, in pixels
  double beam_radius = getProperty("BeamRadius");

  // Get the number of monitors. We assume that all monitors are stored in the first spectra
  const int numSpec = static_cast<int>(inputWSWvl->getNumberHistograms());

  // Set up the progress reporting object
  Progress progress(this,0.0,1.0,max_iteration);

  EventWorkspace_const_sptr inputEventWS = boost::dynamic_pointer_cast<const EventWorkspace>(inputWSWvl);
  if(inputEventWS)
  {
    std::vector<double> y_values(numSpec);
    std::vector<double> e_values(numSpec);

    PARALLEL_FOR_NO_WSP_CHECK()
    for (int i = 0; i < numSpec; i++)
    {
      double sum_i(0), err_i(0);
      progress.report("Integrating events");
      const EventList& el = inputEventWS->getEventList(i);
      el.integrate(0,0,true,sum_i,err_i);
      y_values[i] = sum_i;
      e_values[i] = err_i;
    }

    IAlgorithm_sptr algo = createChildAlgorithm("CreateWorkspace", 0.7, 1.0);
    algo->setProperty< std::vector<double> >("DataX", std::vector<double>(2,0.0) );
    algo->setProperty< std::vector<double> >("DataY", y_values );
    algo->setProperty< std::vector<double> >("DataE", e_values );
    algo->setProperty<int>("NSpec", numSpec );
    algo->execute();

    inputWS = algo->getProperty("OutputWorkspace");
    WorkspaceFactory::Instance().initializeFromParent(inputWSWvl, inputWS, false);
  }
  else
  {

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

      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);
      }

      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 >= wi_min) && (i <= wi_max))
        {
          outputWS->getEventList(i).scaleTof(factor);
          if( factor < 0 )
          {
            outputWS->getEventList(i).reverse();
          }
        }
        m_progress->report("Scaling X");
        PARALLEL_END_INTERUPT_REGION
      }
      PARALLEL_CHECK_INTERUPT_REGION

      outputWS->clearMRU();
    }

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

double DiffractionEventCalibrateDetectors::intensity(
    double x, double y, double z, double rotx, double roty, double rotz,
    std::string detname, std::string inname, std::string outname,
    std::string peakOpt, std::string rb_param, std::string groupWSName) {

  EventWorkspace_sptr inputW = boost::dynamic_pointer_cast<EventWorkspace>(
      AnalysisDataService::Instance().retrieve(inname));

  bool debug = true;
  CPUTimer tim;

  movedetector(x, y, z, rotx, roty, rotz, detname, inputW);
  if (debug)
    std::cout << tim << " to movedetector()\n";

  IAlgorithm_sptr alg3 = createChildAlgorithm("ConvertUnits");
  alg3->setProperty<EventWorkspace_sptr>("InputWorkspace", inputW);
  alg3->setPropertyValue("OutputWorkspace", outname);
  alg3->setPropertyValue("Target", "dSpacing");
  alg3->executeAsChildAlg();
  MatrixWorkspace_sptr outputW = alg3->getProperty("OutputWorkspace");

  if (debug)
    std::cout << tim << " to ConvertUnits\n";

  IAlgorithm_sptr alg4 = createChildAlgorithm("DiffractionFocussing");
  alg4->setProperty<MatrixWorkspace_sptr>("InputWorkspace", outputW);
  alg4->setProperty<MatrixWorkspace_sptr>("OutputWorkspace", outputW);
  alg4->setPropertyValue("GroupingFileName", "");
  alg4->setPropertyValue("GroupingWorkspace", groupWSName);
  alg4->executeAsChildAlg();
  outputW = alg4->getProperty("OutputWorkspace");

  // Remove file
  if (debug)
    std::cout << tim << " to DiffractionFocussing\n";

  IAlgorithm_sptr alg5 = createChildAlgorithm("Rebin");
  alg5->setProperty<MatrixWorkspace_sptr>("InputWorkspace", outputW);
  alg5->setProperty<MatrixWorkspace_sptr>("OutputWorkspace", outputW);
  alg5->setPropertyValue("Params", rb_param);
  alg5->executeAsChildAlg();
  outputW = alg5->getProperty("OutputWorkspace");

  if (debug)
    std::cout << tim << " to Rebin\n";

  // Find point of peak centre
  const MantidVec &yValues = outputW->readY(0);
  auto it = std::max_element(yValues.begin(), yValues.end());
  double peakHeight = *it;
  if (peakHeight == 0)
    return -0.000;
  double peakLoc = outputW->readX(0)[it - yValues.begin()];

  IAlgorithm_sptr fit_alg;
  try {
    // set the ChildAlgorithm no to log as this will be run once per spectra
    fit_alg = createChildAlgorithm("Fit", -1, -1, false);
  } catch (Exception::NotFoundError &) {
    g_log.error("Can't locate Fit algorithm");
    throw;
  }
  std::ostringstream fun_str;
  fun_str << "name=Gaussian,Height=" << peakHeight
          << ",Sigma=0.01,PeakCentre=" << peakLoc;
  fit_alg->setProperty("Function", fun_str.str());
  fit_alg->setProperty("InputWorkspace", outputW);
  fit_alg->setProperty("WorkspaceIndex", 0);
  fit_alg->setProperty("StartX", outputW->readX(0)[0]);
  fit_alg->setProperty("EndX", outputW->readX(0)[outputW->blocksize()]);
  fit_alg->setProperty("MaxIterations", 200);
  fit_alg->setProperty("Output", "fit");
  fit_alg->executeAsChildAlg();

  if (debug)
    std::cout << tim << " to Fit\n";

  std::vector<double> params; // = fit_alg->getProperty("Parameters");
  Mantid::API::IFunction_sptr fun_res = fit_alg->getProperty("Function");
  for (size_t i = 0; i < fun_res->nParams(); ++i) {
    params.push_back(fun_res->getParameter(i));
  }
  peakHeight = params[0];
  peakLoc = params[1];

  movedetector(-x, -y, -z, -rotx, -roty, -rotz, detname, inputW);

  if (debug)
    std::cout << tim << " to movedetector()\n";

  // Optimize C/peakheight + |peakLoc-peakOpt|  where C is scaled by number of
  // events
  EventWorkspace_const_sptr inputE =
      boost::dynamic_pointer_cast<const EventWorkspace>(inputW);
  return (static_cast<int>(inputE->getNumberEvents()) / 1.e6) / peakHeight +
         std::fabs(peakLoc - boost::lexical_cast<double>(peakOpt));
}

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

/** Executes the rebin algorithm
*
*  @throw runtime_error Thrown if the bin range does not intersect the range of
*the input workspace
*/
void Rebin::exec() {
  // Get the input workspace
  MatrixWorkspace_sptr inputWS = getProperty("InputWorkspace");
  MatrixWorkspace_sptr outputWS = getProperty("OutputWorkspace");

  // Are we preserving event workspace-iness?
  bool PreserveEvents = getProperty("PreserveEvents");

  // Rebinning in-place
  bool inPlace = (inputWS == outputWS);

  std::vector<double> rbParams =
      rebinParamsFromInput(getProperty("Params"), *inputWS, g_log);

  const bool dist = inputWS->isDistribution();
  const bool isHist = inputWS->isHistogramData();

  // workspace independent determination of length
  const int histnumber = static_cast<int>(inputWS->getNumberHistograms());

  //-------------------------------------------------------

  bool fullBinsOnly = getProperty("FullBinsOnly");

  MantidVecPtr XValues_new;
  // create new output X axis
  const int ntcnew = VectorHelper::createAxisFromRebinParams(
      rbParams, XValues_new.access(), true, fullBinsOnly);

  //---------------------------------------------------------------------------------
  // Now, determine if the input workspace is actually an EventWorkspace
  EventWorkspace_const_sptr eventInputWS =
      boost::dynamic_pointer_cast<const EventWorkspace>(inputWS);

  if (eventInputWS != NULL) {
    //------- EventWorkspace as input -------------------------------------
    EventWorkspace_sptr eventOutputWS =
        boost::dynamic_pointer_cast<EventWorkspace>(outputWS);

    if (inPlace && PreserveEvents) {
      // -------------Rebin in-place, preserving events
      // ----------------------------------------------
      // This only sets the X axis. Actual rebinning will be done upon data
      // access.
      eventOutputWS->setAllX(XValues_new);
      this->setProperty(
          "OutputWorkspace",
          boost::dynamic_pointer_cast<MatrixWorkspace>(eventOutputWS));
    } else if (!inPlace && PreserveEvents) {
      // -------- NOT in-place, but you want to keep events for some reason.
      // ----------------------
      // Must copy the event workspace to a new EventWorkspace (and bin that).

      // Make a brand new EventWorkspace
      eventOutputWS = boost::dynamic_pointer_cast<EventWorkspace>(
          API::WorkspaceFactory::Instance().create(
              "EventWorkspace", inputWS->getNumberHistograms(), 2, 1));
      // Copy geometry over.
      API::WorkspaceFactory::Instance().initializeFromParent(
          inputWS, eventOutputWS, false);
      // You need to copy over the data as well.
      eventOutputWS->copyDataFrom((*eventInputWS));

      // This only sets the X axis. Actual rebinning will be done upon data
      // access.
      eventOutputWS->setAllX(XValues_new);

      // Cast to the matrixOutputWS and save it
      this->setProperty(
          "OutputWorkspace",
          boost::dynamic_pointer_cast<MatrixWorkspace>(eventOutputWS));
    } else {
      //--------- Different output, OR you're inplace but not preserving Events
      //--- create a Workspace2D -------
      g_log.information() << "Creating a Workspace2D from the EventWorkspace "
                          << eventInputWS->getName() << ".\n";

      // Create a Workspace2D
      // This creates a new Workspace2D through a torturous route using the
      // WorkspaceFactory.
      // The Workspace2D is created with an EMPTY CONSTRUCTOR
      outputWS = WorkspaceFactory::Instance().create("Workspace2D", histnumber,
                                                     ntcnew, ntcnew - 1);
      WorkspaceFactory::Instance().initializeFromParent(inputWS, outputWS,
                                                        true);

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

      // Go through all the histograms and set the data
      PARALLEL_FOR3(inputWS, eventInputWS, outputWS)
      for (int i = 0; i < histnumber; ++i) {
        PARALLEL_START_INTERUPT_REGION

        // Set the X axis for each output histogram
//.........这里部分代码省略.........

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

double EQSANSTofStructure::getTofOffset(EventWorkspace_const_sptr inputWS,
                                        bool frame_skipping) {
  //# Storage for chopper information read from the logs
  double chopper_set_phase[4] = {0, 0, 0, 0};
  double chopper_speed[4] = {0, 0, 0, 0};
  double chopper_actual_phase[4] = {0, 0, 0, 0};
  double chopper_wl_1[4] = {0, 0, 0, 0};
  double chopper_wl_2[4] = {0, 0, 0, 0};
  double frame_wl_1 = 0;
  double frame_srcpulse_wl_1 = 0;
  double frame_wl_2 = 0;
  double chopper_srcpulse_wl_1[4] = {0, 0, 0, 0};
  double chopper_frameskip_wl_1[4] = {0, 0, 0, 0};
  double chopper_frameskip_wl_2[4] = {0, 0, 0, 0};
  double chopper_frameskip_srcpulse_wl_1[4] = {0, 0, 0, 0};

  // Calculate the frame width
  auto frequencyLog = dynamic_cast<TimeSeriesProperty<double> *>(
      inputWS->run().getLogData("frequency"));
  if (!frequencyLog) {
    throw std::runtime_error("Frequency log not found.");
  }
  double frequency = frequencyLog->getStatistics().mean;
  double tof_frame_width = 1.0e6 / frequency;

  double tmp_frame_width = tof_frame_width;
  if (frame_skipping)
    tmp_frame_width *= 2.0;

  // Choice of parameter set
  int m_set = 0;
  if (frame_skipping)
    m_set = 1;

  bool first = true;
  bool first_skip = true;
  double frameskip_wl_1 = 0;
  double frameskip_srcpulse_wl_1 = 0;
  double frameskip_wl_2 = 0;

  for (int i = 0; i < 4; i++) {
    // Read chopper information
    std::ostringstream phase_str;
    phase_str << "Phase" << i + 1;
    auto log = dynamic_cast<TimeSeriesProperty<double> *>(
        inputWS->run().getLogData(phase_str.str()));
    if (!log) {
      throw std::runtime_error("Phase log not found.");
    }
    chopper_set_phase[i] = log->getStatistics().mean;
    std::ostringstream speed_str;
    speed_str << "Speed" << i + 1;
    log = dynamic_cast<TimeSeriesProperty<double> *>(
        inputWS->run().getLogData(speed_str.str()));
    if (!log) {
      throw std::runtime_error("Speed log not found.");
    }
    chopper_speed[i] = log->getStatistics().mean;

    // Only process choppers with non-zero speed
    if (chopper_speed[i] <= 0)
      continue;

    chopper_actual_phase[i] =
        chopper_set_phase[i] - CHOPPER_PHASE_OFFSET[m_set][i];

    while (chopper_actual_phase[i] < 0)
      chopper_actual_phase[i] += tmp_frame_width;

    double x1 =
        (chopper_actual_phase[i] -
         (tmp_frame_width * 0.5 * CHOPPER_ANGLE[i] / 360.)); // opening edge
    double x2 =
        (chopper_actual_phase[i] +
         (tmp_frame_width * 0.5 * CHOPPER_ANGLE[i] / 360.)); // closing edge
    if (!frame_skipping)                                     // not skipping
    {
      while (x1 < 0) {
        x1 += tmp_frame_width;
        x2 += tmp_frame_width;
      }
    }

    if (x1 > 0) {
      chopper_wl_1[i] = 3.9560346 * x1 / CHOPPER_LOCATION[i];
      chopper_srcpulse_wl_1[i] =
          3.9560346 * (x1 - chopper_wl_1[i] * PULSEWIDTH) / CHOPPER_LOCATION[i];
    } else
      chopper_wl_1[i] = chopper_srcpulse_wl_1[i] = 0.;

    if (x2 > 0)
      chopper_wl_2[i] = 3.9560346 * x2 / CHOPPER_LOCATION[i];
    else
      chopper_wl_2[i] = 0.;

    if (first) {
      frame_wl_1 = chopper_wl_1[i];
      frame_srcpulse_wl_1 = chopper_srcpulse_wl_1[i];
      frame_wl_2 = chopper_wl_2[i];
      first = false;
//.........这里部分代码省略.........

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

/** 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();
  const 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 != nullptr) {
    if (m_preserveEvents) {
      if (inPlace) {
        g_log.debug() << "Rebinning event workspace in place\n";
      } else {
        g_log.debug() << "Rebinning event workspace out of place\n";
        outputWS = inputWS->clone();
      }
      auto outputEventWS =
          boost::dynamic_pointer_cast<EventWorkspace>(outputWS);

      if (common_limits) {
        // get the delta from the first since they are all the same
        BinEdges xValues(0);
        const double delta = this->determineBinning(xValues.mutableRawData(),
                                                    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_FOR_IF(Kernel::threadSafe(*inputEventWS, *outputWS))
        for (int wkspIndex = 0; wkspIndex < numSpectra; ++wkspIndex) {
          PARALLEL_START_INTERUPT_REGION
          BinEdges xValues(0);
          const double delta = this->determineBinning(
              xValues.mutableRawData(), xmins[wkspIndex], xmaxs[wkspIndex]);
          g_log.debug() << "delta[wkspindex=" << wkspIndex << "] = " << delta
                        << " xmin=" << xmins[wkspIndex]
                        << " xmax=" << xmaxs[wkspIndex] << "\n";
          outputEventWS->setHistogram(wkspIndex, xValues);
          prog.report(name()); // Report progress
          PARALLEL_END_INTERUPT_REGION
        }
        PARALLEL_CHECK_INTERUPT_REGION
      }
    }    // end if (m_preserveEvents)
    else // event workspace -> matrix workspace
    {
      //--------- Different output, OR you're inplace but not preserving Events
      g_log.information() << "Creating a Workspace2D from the EventWorkspace "
                          << inputEventWS->getName() << ".\n";
      outputWS = create<DataObjects::Workspace2D>(
          *inputWS, numSpectra, HistogramData::BinEdges(m_numBins + 1));

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

      // Go through all the histograms and set the data
      PARALLEL_FOR_IF(Kernel::threadSafe(*inputEventWS, *outputWS))
      for (int wkspIndex = 0; wkspIndex < numSpectra; ++wkspIndex) {
//.........这里部分代码省略.........

    /** Executes the algorithm
     *
     *  @throw Exception::FileError If the grouping file cannot be opened or read successfully
     */
    void GetDetOffsetsMultiPeaks::exec()
    {
      const double BAD_OFFSET(1000.); // mark things that didn't work with this

      MatrixWorkspace_sptr inputW=getProperty("InputWorkspace");
      double maxOffset=getProperty("MaxOffset");
      int nspec=static_cast<int>(inputW->getNumberHistograms());
      // Create the output OffsetsWorkspace
      OffsetsWorkspace_sptr outputW(new OffsetsWorkspace(inputW->getInstrument()));
      // determine min/max d-spacing of the workspace
      double wkspDmin, wkspDmax;
      inputW->getXMinMax(wkspDmin, wkspDmax);
      // Create the output MaskWorkspace
      MatrixWorkspace_sptr maskWS(new MaskWorkspace(inputW->getInstrument()));
      //To get the workspace index from the detector ID
      detid2index_map * pixel_to_wi = maskWS->getDetectorIDToWorkspaceIndexMap(true);
      // the peak positions and where to fit
      std::vector<double> peakPositions = getProperty("DReference");
      std::sort(peakPositions.begin(), peakPositions.end());
      std::vector<double> fitWindows = generateWindows(wkspDmin, wkspDmax, peakPositions, this->getProperty("FitWindowMaxWidth"));
      g_log.information() << "windows : ";
      if (fitWindows.empty())
      {
        g_log.information() << "empty\n";
      }
      else
      {
        for (std::vector<double>::const_iterator it = fitWindows.begin(); it != fitWindows.end(); ++it)
          g_log.information() << *it << " ";
        g_log.information() << "\n";
      }

      // some shortcuts for event workspaces
      EventWorkspace_const_sptr eventW = boost::dynamic_pointer_cast<const EventWorkspace>( inputW );
      bool isEvent = false;
      if (eventW)
        isEvent = true;

      // cache the peak and background function names
      m_peakType = this->getPropertyValue("PeakFunction");
      m_backType = this->getPropertyValue("BackgroundType");
      // the maximum allowable chisq value for an individual peak fit
      m_maxChiSq = this->getProperty("MaxChiSq");

      // Fit all the spectra with a gaussian
      Progress prog(this, 0, 1.0, nspec);
      // cppcheck-suppress syntaxError
      PRAGMA_OMP(parallel for schedule(dynamic, 1) )
      for (int wi=0;wi<nspec;++wi)
      {
        PARALLEL_START_INTERUPT_REGION
        double offset = 0.0;
        double fitSum = 0.0;
        // checks for dead detectors
        if ((isEvent) && (eventW->getEventList(wi).empty()))
        {
          // dead detector will be masked
          offset = BAD_OFFSET;
        }
        else {
          const MantidVec& Y = inputW->readY(wi);
          const int YLength = static_cast<int>(Y.size());
          double sumY = 0.0;
          for (int i = 0; i < YLength; i++) sumY += Y[i];
          if (sumY < 1.e-30)
          {
            // Dead detector will be masked
            offset=BAD_OFFSET;
          }
        }
        if (offset < 10.)
        {
          // Fit the peak
          std::vector<double> peakPosToFit, peakPosFitted, chisq;
          size_t nparams;
          double minD, maxD;
          fitSpectra(wi, inputW, peakPositions, fitWindows, nparams, minD, maxD, peakPosToFit, peakPosFitted, chisq);
          if (nparams > 0)
          {
          //double * params = new double[2*nparams+1];
          double params[153];
          if(nparams > 50) nparams = 50;
          params[0] = static_cast<double>(nparams);
          params[1] = minD;
          params[2] = maxD;
          for (size_t i = 0; i < nparams; i++)
          {
            params[i+3] = peakPosToFit[i];
          }
          for (size_t i = 0; i < nparams; i++)
          {
            params[i+3+nparams] = peakPosFitted[i];
          }
          for (size_t i = 0; i < nparams; i++)
          {
            params[i+3+2*nparams] = chisq[i];
//.........这里部分代码省略.........

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

}

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