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

//.........这里部分代码省略.........
      m_outputW = boost::dynamic_pointer_cast<MatrixWorkspace>(m_outputEW);
    } else {
      g_log.information() << "Not compressing event list\n";
      doSortEvents(m_outputW); // still sort to help some thing out
    }
  }
  m_progress->report();

  if (xmin > 0. || xmax > 0.) {
    double tempmin;
    double tempmax;
    m_outputW->getXMinMax(tempmin, tempmax);

    g_log.information() << "running CropWorkspace(TOFmin=" << xmin
                        << ", TOFmax=" << xmax << ") started at "
                        << Kernel::DateAndTime::getCurrentTime() << "\n";
    API::IAlgorithm_sptr cropAlg = createChildAlgorithm("CropWorkspace");
    cropAlg->setProperty("InputWorkspace", m_outputW);
    cropAlg->setProperty("OutputWorkspace", m_outputW);
    if ((xmin > 0.) && (xmin > tempmin))
      cropAlg->setProperty("Xmin", xmin);
    if ((xmax > 0.) && (xmax < tempmax))
      cropAlg->setProperty("Xmax", xmax);
    cropAlg->executeAsChildAlg();
    m_outputW = cropAlg->getProperty("OutputWorkspace");
    m_outputEW = boost::dynamic_pointer_cast<EventWorkspace>(m_outputW);
  }
  m_progress->report();

  // filter the input events if appropriate
  double removePromptPulseWidth = getProperty("RemovePromptPulseWidth");
  if (removePromptPulseWidth > 0.) {
    m_outputEW = boost::dynamic_pointer_cast<EventWorkspace>(m_outputW);
    if (m_outputEW->getNumberEvents() > 0) {
      g_log.information() << "running RemovePromptPulse(Width="
                          << removePromptPulseWidth << ") started at "
                          << Kernel::DateAndTime::getCurrentTime() << "\n";
      API::IAlgorithm_sptr filterPAlg =
          createChildAlgorithm("RemovePromptPulse");
      filterPAlg->setProperty("InputWorkspace", m_outputW);
      filterPAlg->setProperty("OutputWorkspace", m_outputW);
      filterPAlg->setProperty("Width", removePromptPulseWidth);
      filterPAlg->executeAsChildAlg();
      m_outputW = filterPAlg->getProperty("OutputWorkspace");
      m_outputEW = boost::dynamic_pointer_cast<EventWorkspace>(m_outputW);
    } else {
      g_log.information("skipping RemovePromptPulse on empty EventWorkspace");
    }
  }
  m_progress->report();

  if (maskBinTableWS) {
    g_log.information() << "running MaskBinsFromTable started at "
                        << Kernel::DateAndTime::getCurrentTime() << "\n";
    API::IAlgorithm_sptr alg = createChildAlgorithm("MaskBinsFromTable");
    alg->setProperty("InputWorkspace", m_outputW);
    alg->setProperty("OutputWorkspace", m_outputW);
    alg->setProperty("MaskingInformation", maskBinTableWS);
    alg->executeAsChildAlg();
    m_outputW = alg->getProperty("OutputWorkspace");
    m_outputEW = boost::dynamic_pointer_cast<EventWorkspace>(m_outputW);
  }
  m_progress->report();

  if (m_maskWS) {
    g_log.information() << "running MaskDetectors started at "

/// Executes the algorithm for events
void UnaryOperation::execEvent() {
  g_log.information("Processing event workspace");

  const MatrixWorkspace_const_sptr matrixInputWS =
      this->getProperty(inputPropName());
  EventWorkspace_const_sptr inputWS =
      boost::dynamic_pointer_cast<const EventWorkspace>(matrixInputWS);

  // generate the output workspace pointer
  API::MatrixWorkspace_sptr matrixOutputWS =
      this->getProperty(outputPropName());
  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);
  }

  // Now fetch any properties defined by concrete algorithm
  retrieveProperties();

  int64_t numHistograms = static_cast<int64_t>(inputWS->getNumberHistograms());
  API::Progress prog = API::Progress(this, 0.0, 1.0, numHistograms);
  PARALLEL_FOR1(outputWS)
  for (int64_t i = 0; i < numHistograms; ++i) {
    PARALLEL_START_INTERUPT_REGION
    // switch to weighted events if needed, and use the appropriate helper
    // function
    EventList *evlist = outputWS->getEventListPtr(i);
    switch (evlist->getEventType()) {
    case TOF:
      // Switch to weights if needed.
      evlist->switchTo(WEIGHTED);
    /* no break */
    // Fall through

    case WEIGHTED:
      unaryOperationEventHelper(evlist->getWeightedEvents());
      break;

    case WEIGHTED_NOTIME:
      unaryOperationEventHelper(evlist->getWeightedEventsNoTime());
      break;
    }

    prog.report();
    PARALLEL_END_INTERUPT_REGION
  }
  PARALLEL_CHECK_INTERUPT_REGION

  outputWS->clearMRU();
  if (inputWS->getNumberEvents() != outputWS->getNumberEvents()) {
    g_log.information() << "Number of events has changed!!!" << std::endl;
  }
}

/** Execute the algorithm.
 */
void IntegrateEllipsoids::exec() {
  // get the input workspace
  MatrixWorkspace_sptr wksp = getProperty("InputWorkspace");

  EventWorkspace_sptr eventWS =
      boost::dynamic_pointer_cast<EventWorkspace>(wksp);
  Workspace2D_sptr histoWS = boost::dynamic_pointer_cast<Workspace2D>(wksp);
  if (!eventWS && !histoWS) {
    throw std::runtime_error("IntegrateEllipsoids needs either a "
                             "EventWorkspace or Workspace2D as input.");
  }

  // error out if there are not events
  if (eventWS && eventWS->getNumberEvents() <= 0) {
    throw std::runtime_error(
        "IntegrateEllipsoids does not work for empty event lists");
  }

  PeaksWorkspace_sptr in_peak_ws = getProperty("PeaksWorkspace");
  if (!in_peak_ws) {
    throw std::runtime_error("Could not read the peaks workspace");
  }

  double radius = getProperty("RegionRadius");
  int numSigmas = getProperty("NumSigmas");
  double cutoffIsigI = getProperty("CutoffIsigI");
  bool specify_size = getProperty("SpecifySize");
  double peak_radius = getProperty("PeakSize");
  double back_inner_radius = getProperty("BackgroundInnerSize");
  double back_outer_radius = getProperty("BackgroundOuterSize");
  bool hkl_integ = getProperty("IntegrateInHKL");
  bool integrateEdge = getProperty("IntegrateIfOnEdge");
  if (!integrateEdge) {
    // This only fails in the unit tests which say that MaskBTP is not
    // registered
    try {
      runMaskDetectors(in_peak_ws, "Tube", "edges");
      runMaskDetectors(in_peak_ws, "Pixel", "edges");
    } catch (...) {
      g_log.error("Can't execute MaskBTP algorithm for this instrument to set "
                  "edge for IntegrateIfOnEdge option");
    }
    calculateE1(in_peak_ws->detectorInfo()); // fill E1Vec for use in detectorQ
  }

  Mantid::DataObjects::PeaksWorkspace_sptr peak_ws =
      getProperty("OutputWorkspace");
  if (peak_ws != in_peak_ws) {
    peak_ws = in_peak_ws->clone();
  }

  // get UBinv and the list of
  // peak Q's for the integrator
  std::vector<Peak> &peaks = peak_ws->getPeaks();
  size_t n_peaks = peak_ws->getNumberPeaks();
  size_t indexed_count = 0;
  std::vector<V3D> peak_q_list;
  std::vector<std::pair<double, V3D>> qList;
  std::vector<V3D> hkl_vectors;
  for (size_t i = 0; i < n_peaks; i++) // Note: we skip un-indexed peaks
  {
    V3D hkl(peaks[i].getH(), peaks[i].getK(), peaks[i].getL());
    if (Geometry::IndexingUtils::ValidIndex(hkl, 1.0)) // use tolerance == 1 to
                                                       // just check for (0,0,0)
    {
      peak_q_list.emplace_back(peaks[i].getQLabFrame());
      qList.emplace_back(1., V3D(peaks[i].getQLabFrame()));
      V3D miller_ind(static_cast<double>(boost::math::iround<double>(hkl[0])),
                     static_cast<double>(boost::math::iround<double>(hkl[1])),
                     static_cast<double>(boost::math::iround<double>(hkl[2])));
      hkl_vectors.push_back(miller_ind);
      indexed_count++;
    }
  }

  if (indexed_count < 3) {
    throw std::runtime_error(
        "At least three linearly independent indexed peaks are needed.");
  }
  // Get UB using indexed peaks and
  // lab-Q vectors
  Matrix<double> UB(3, 3, false);
  Geometry::IndexingUtils::Optimize_UB(UB, hkl_vectors, peak_q_list);
  Matrix<double> UBinv(UB);
  UBinv.Invert();
  UBinv *= (1.0 / (2.0 * M_PI));

  std::vector<double> PeakRadiusVector(n_peaks, peak_radius);
  std::vector<double> BackgroundInnerRadiusVector(n_peaks, back_inner_radius);
  std::vector<double> BackgroundOuterRadiusVector(n_peaks, back_outer_radius);
  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");
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  }

  const std::string emodeStr = getProperty("EMode");
  double efixedProp = getProperty("EFixed"),efixed;

  if( efixedProp == EMPTY_DBL() )
  {
    if (emodeStr == "Direct")
    {
      // Check if it has been store on the run object for this workspace
      if( this->inputWS->run().hasProperty("Ei"))
      {
        Kernel::Property* eiprop = this->inputWS->run().getProperty("Ei");
        efixedProp = boost::lexical_cast<double>(eiprop->value());
        g_log.debug() << "Using stored Ei value " << efixedProp << "\n";
      }
      else
      {
        throw std::invalid_argument("No Ei value has been set or stored within the run information.");
      }
    }
    else
    {
      // If not specified, will try to get Ef from the parameter file for indirect geometry, 
      // but it will be done for each spectrum separately, in case of different analyzer crystals
    }
  }

  // Get the parameter map
  const ParameterMap& pmap = outputWS->constInstrumentParameters();

  int64_t numHistograms = static_cast<int64_t>(inputWS->getNumberHistograms());
  API::Progress prog = API::Progress(this, 0.0, 1.0, numHistograms);
  PARALLEL_FOR1(outputWS)
  for (int64_t i=0; i < numHistograms; ++i)
  {
    PARALLEL_START_INTERUPT_REGION
    
    double Efi = 0;
    // Now get the detector object for this histogram to check if monitor
    // or to get Ef for indirect geometry
    if (emodeStr == "Indirect") 
    {
      if ( efixedProp != EMPTY_DBL()) Efi = efixedProp;
      else try 
      {
        IDetector_const_sptr det = inputWS->getDetector(i);
        if (!det->isMonitor())
        {
            try
            {
              Parameter_sptr par = pmap.getRecursive(det.get(),"Efixed");
              if (par)
              {
                Efi = par->value<double>();
                g_log.debug() << "Detector: " << det->getID() << " EFixed: " << Efi << "\n";
              }
            }
            catch (std::runtime_error&) { /* Throws if a DetectorGroup, use single provided value */ }
        }

      }
      catch(std::runtime_error&) { g_log.information() << "Workspace Index " << i << ": cannot find detector" << "\n"; }
    }

    if (emodeStr == "Indirect") efixed=Efi;
    else efixed=efixedProp;

    //Do the correction
		EventList *evlist=outputWS->getEventListPtr(i);
    switch (evlist->getEventType())
    {
      case TOF:
        //Switch to weights if needed.
        evlist->switchTo(WEIGHTED);
        /* no break */
        // Fall through
  
      case WEIGHTED:
        correctKiKfEventHelper(evlist->getWeightedEvents(), efixed,emodeStr);
        break;
  
      case WEIGHTED_NOTIME:
        correctKiKfEventHelper(evlist->getWeightedEventsNoTime(), efixed,emodeStr);
        break;
    }


    prog.report();
    PARALLEL_END_INTERUPT_REGION
  }
  PARALLEL_CHECK_INTERUPT_REGION

  outputWS->clearMRU();
  if (inputWS->getNumberEvents( ) != outputWS->getNumberEvents( ))
  {
    g_log.information() <<"Ef <= 0 or Ei <= 0 for "<<inputWS->getNumberEvents( )-outputWS->getNumberEvents( )<<" events, out of "<<inputWS->getNumberEvents( )<<std::endl;
    if ( efixedProp == EMPTY_DBL()) g_log.information()<<"Try to set fixed energy"<<std::endl ;
  }
}

  /** Execute the algorithm.
   */
  void IntegrateEllipsoids::exec()
  {
    // get the input workspace
    EventWorkspace_sptr wksp = getProperty("InputWorkspace");

    // this only works for unweighted events
    if (wksp->getEventType() != API::TOF)
    {
      throw std::runtime_error("IntegrateEllipsoids only works for raw events");
    }
    // error out if there are not events
    if (wksp->getNumberEvents() <= 0)
    {
      throw std::runtime_error("IntegrateEllipsoids does not work for empty event lists");
    }

    PeaksWorkspace_sptr in_peak_ws;
    in_peak_ws = boost::dynamic_pointer_cast<PeaksWorkspace>(
         AnalysisDataService::Instance().retrieve( getProperty("PeaksWorkspace")) );

    if (!in_peak_ws)
    {
      throw std::runtime_error("Could not read the peaks workspace");
    }

    Mantid::DataObjects::PeaksWorkspace_sptr peak_ws = getProperty("OutputWorkspace");
    if ( peak_ws != in_peak_ws )
    {
      peak_ws = in_peak_ws->clone();
    }
                                                   // get UBinv and the list of
                                                   // peak Q's for the integrator
    std::vector<Peak> & peaks = peak_ws->getPeaks();
    size_t n_peaks            = peak_ws->getNumberPeaks();
    size_t indexed_count      = 0;
    std::vector<V3D> peak_q_list;
    std::vector<V3D> hkl_vectors;
    for ( size_t i = 0; i < n_peaks; i++ )         // Note: we skip un-indexed peaks
    {
      V3D hkl( peaks[i].getH(), peaks[i].getK(), peaks[i].getL() );  
      if ( Geometry::IndexingUtils::ValidIndex( hkl, 1.0 ) )    // use tolerance == 1 to 
                                                                // just check for (0,0,0) 
      {
        peak_q_list.push_back( V3D( peaks[i].getQLabFrame() ) );
        V3D miller_ind( (double)boost::math::iround<double>(hkl[0]), 
                        (double)boost::math::iround<double>(hkl[1]),
                        (double)boost::math::iround<double>(hkl[2]) );
        hkl_vectors.push_back( V3D(miller_ind) );
        indexed_count++;
      }
    }

    if ( indexed_count < 3 )
    {
      throw std::runtime_error(
            "At least three linearly independent indexed peaks are needed.");
    }
                                             // Get UB using indexed peaks and
                                             // lab-Q vectors
    Matrix<double> UB(3,3,false);
    Geometry::IndexingUtils::Optimize_UB( UB, hkl_vectors, peak_q_list );
    Matrix<double> UBinv( UB );
    UBinv.Invert();
    UBinv *= (1.0/(2.0 * M_PI));

    double radius            = getProperty( "RegionRadius" );
    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);

//.........这里部分代码省略.........

/** Execute the algorithm.
 */
void IntegrateEllipsoids::exec() {
  // get the input workspace
  MatrixWorkspace_sptr wksp = getProperty("InputWorkspace");

  EventWorkspace_sptr eventWS =
      boost::dynamic_pointer_cast<EventWorkspace>(wksp);
  Workspace2D_sptr histoWS = boost::dynamic_pointer_cast<Workspace2D>(wksp);
  if (!eventWS && !histoWS) {
    throw std::runtime_error("IntegrateEllipsoids needs either a "
                             "EventWorkspace or Workspace2D as input.");
  }

  // error out if there are not events
  if (eventWS && eventWS->getNumberEvents() <= 0) {
    throw std::runtime_error(
        "IntegrateEllipsoids does not work for empty event lists");
  }

  PeaksWorkspace_sptr in_peak_ws = getProperty("PeaksWorkspace");
  if (!in_peak_ws) {
    throw std::runtime_error("Could not read the peaks workspace");
  }

  double radius_m = getProperty("RegionRadius");
  double radius_s = getProperty("SatelliteRegionRadius");
  int numSigmas = getProperty("NumSigmas");
  double cutoffIsigI = getProperty("CutoffIsigI");
  bool specify_size = getProperty("SpecifySize");
  double peak_radius = getProperty("PeakSize");
  double sate_peak_radius = getProperty("SatellitePeakSize");
  double back_inner_radius = getProperty("BackgroundInnerSize");
  double sate_back_inner_radius = getProperty("SatelliteBackgroundInnerSize");
  double back_outer_radius = getProperty("BackgroundOuterSize");
  double sate_back_outer_radius = getProperty("SatelliteBackgroundOuterSize");
  bool hkl_integ = getProperty("IntegrateInHKL");
  bool integrateEdge = getProperty("IntegrateIfOnEdge");
  bool adaptiveQBackground = getProperty("AdaptiveQBackground");
  double adaptiveQMultiplier = getProperty("AdaptiveQMultiplier");
  double adaptiveQBackgroundMultiplier = 0.0;
  bool useOnePercentBackgroundCorrection =
      getProperty("UseOnePercentBackgroundCorrection");
  if (adaptiveQBackground)
    adaptiveQBackgroundMultiplier = adaptiveQMultiplier;
  if (!integrateEdge) {
    // This only fails in the unit tests which say that MaskBTP is not
    // registered
    try {
      runMaskDetectors(in_peak_ws, "Tube", "edges");
      runMaskDetectors(in_peak_ws, "Pixel", "edges");
    } catch (...) {
      g_log.error("Can't execute MaskBTP algorithm for this instrument to set "
                  "edge for IntegrateIfOnEdge option");
    }
    calculateE1(in_peak_ws->detectorInfo()); // fill E1Vec for use in detectorQ
  }

  Mantid::DataObjects::PeaksWorkspace_sptr peak_ws =
      getProperty("OutputWorkspace");
  if (peak_ws != in_peak_ws)
    peak_ws = in_peak_ws->clone();

  // get UBinv and the list of
  // peak Q's for the integrator
  std::vector<Peak> &peaks = peak_ws->getPeaks();
  size_t n_peaks = peak_ws->getNumberPeaks();
  size_t indexed_count = 0;
  std::vector<V3D> peak_q_list;
  std::vector<std::pair<double, V3D>> qList;
  std::vector<V3D> hkl_vectors;
  std::vector<V3D> mnp_vectors;
  int ModDim = 0;
  for (size_t i = 0; i < n_peaks; i++) // Note: we skip un-indexed peaks
  {
    V3D hkl(peaks[i].getIntHKL());
    V3D mnp(peaks[i].getIntMNP());

    if (mnp[0] != 0 && ModDim == 0)
      ModDim = 1;
    if (mnp[1] != 0 && ModDim == 1)
      ModDim = 2;
    if (mnp[2] != 0 && ModDim == 2)
      ModDim = 3;

    // use tolerance == 1 to just check for (0,0,0,0,0,0)
    if (Geometry::IndexingUtils::ValidIndex(hkl, 1.0)) {
      peak_q_list.emplace_back(peaks[i].getQLabFrame());
      qList.emplace_back(1., V3D(peaks[i].getQLabFrame()));
      hkl_vectors.push_back(hkl);
      mnp_vectors.push_back(mnp);
      indexed_count++;
    }
  }

  if (indexed_count < 3)
    throw std::runtime_error(
        "At least three linearly independent indexed peaks are needed.");

  // Get UB using indexed peaks and
//.........这里部分代码省略.........