C++ MatrixWorkspace_sptr::maskWorkspaceIndex方法代码示例

    /**
     * Mask the outlier values to get a better median value.
     * @param medianvec The median value calculated from the current counts.
     * @param countsWS The counts workspace. Any outliers will be masked here.
     * @param indexmap Index map.
     * @returns The number failed.
     */
    int MedianDetectorTest::maskOutliers(const std::vector<double> medianvec, API::MatrixWorkspace_sptr countsWS,std::vector<std::vector<size_t> > indexmap)
    {

      // Fractions of the median
      const double out_lo = getProperty("LowOutlier");
      const double out_hi = getProperty("HighOutlier");

      int numFailed(0);

      bool checkForMask = false;
      Geometry::Instrument_const_sptr instrument = countsWS->getInstrument();
      if (instrument != NULL)
      {
        checkForMask = ((instrument->getSource() != NULL) && (instrument->getSample() != NULL));
      }

      for (size_t i=0; i<indexmap.size();  ++i)
      {
        std::vector<size_t> hists=indexmap.at(i);
        double median=medianvec.at(i);

        PARALLEL_FOR1(countsWS)
        for(int j = 0; j < static_cast<int>(hists.size()); ++j)
        {
          const double value = countsWS->readY(hists.at(j))[0];
          if ((value == 0.) && checkForMask)
          {
            const std::set<detid_t>& detids = countsWS->getSpectrum(hists.at(j))->getDetectorIDs();
            if (instrument->isDetectorMasked(detids))
            {
              numFailed -= 1; // it was already masked
            }
          }
          if( (value < out_lo*median) && (value > 0.0) )
          {
            countsWS->maskWorkspaceIndex(hists.at(j));
            PARALLEL_ATOMIC
            ++numFailed;
          }
          else if( value > out_hi*median )
          {
            countsWS->maskWorkspaceIndex(hists.at(j));
            PARALLEL_ATOMIC
            ++numFailed;
          }
        }
        PARALLEL_CHECK_INTERUPT_REGION
      }

      return numFailed;
    }

 /**
  * Checks if the spectra at the given index of either input workspace is masked. If so then the output spectra has zeroed data
  * and is also masked. 
  * @param lhs :: A pointer to the left-hand operand
  * @param rhs :: A pointer to the right-hand operand
  * @param index :: The workspace index to check
  * @param out :: A pointer to the output workspace
  * @returns True if further processing is not required on the spectra, false if the binary operation should be performed.
  */
 bool BinaryOperation::propagateSpectraMask(const API::MatrixWorkspace_const_sptr lhs, const API::MatrixWorkspace_const_sptr rhs, 
     const int64_t index, API::MatrixWorkspace_sptr out)
 {
   bool continueOp(true);
   IDetector_const_sptr det_lhs, det_rhs;
   try
   {
     det_lhs = lhs->getDetector(index);
     det_rhs = rhs->getDetector(index);
   }
   catch(std::runtime_error &)
   {
   }
   catch(std::domain_error &)
   {
     // try statement will throw a domain_error when the axis is not a spectra axis.
     return continueOp;
   }
   if( (det_lhs && det_lhs->isMasked()) || ( det_rhs && det_rhs->isMasked()) )
   {
     continueOp = false;
     out->maskWorkspaceIndex(index);
   }
   return continueOp;
 }

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