Skip to content
Snippets Groups Projects
BinaryOperation.cpp 42 KiB
Newer Older
  int rhs_nhist = static_cast<int>(rhs->getNumberHistograms());
  int lhs_nhist = static_cast<int>(lhs->getNumberHistograms());

  // Initialize the table; filled with -1 meaning no match
  table->resize(lhs_nhist, -1);

  const detid2index_map rhs_det_to_wi = rhs->getDetectorIDToWorkspaceIndexMap();

  PARALLEL_FOR_NO_WSP_CHECK()
  for (int lhsWI = 0; lhsWI < lhs_nhist; lhsWI++) {
    bool done = false;

    // List of detectors on lhs side
    const std::set<detid_t> &lhsDets =
        lhs->getSpectrum(lhsWI)->getDetectorIDs();

    // ----------------- Matching Workspace Indices and Detector IDs
    // --------------------------------------
    // First off, try to match the workspace indices. Most times, this will be
    // ok right away.
    int64_t rhsWI = lhsWI;
    if (rhsWI < rhs_nhist) // don't go out of bounds
      // Get the detector IDs at that workspace index.
      const std::set<detid_t> &rhsDets =
          rhs->getSpectrum(rhsWI)->getDetectorIDs();

      // Checks that lhsDets is a subset of rhsDets
      if (std::includes(rhsDets.begin(), rhsDets.end(), lhsDets.begin(),
                        lhsDets.end())) {
        // We found the workspace index right away. No need to keep looking
        (*table)[lhsWI] = rhsWI;
        done = true;
    // ----------------- Scrambled Detector IDs with one Detector per Spectrum
    // --------------------------------------
    if (!done && (lhsDets.size() == 1)) {
      // Didn't find it. Try to use the RHS map.

      // First, we have to get the (single) detector ID of the LHS
      std::set<detid_t>::const_iterator lhsDets_it = lhsDets.begin();
      detid_t lhs_detector_ID = *lhsDets_it;

      // Now we use the RHS map to find it. This only works if both the lhs and
      // rhs have 1 detector per pixel
      detid2index_map::const_iterator map_it =
          rhs_det_to_wi.find(lhs_detector_ID);
      if (map_it != rhs_det_to_wi.end()) {
        rhsWI = map_it->second; // This is the workspace index in the RHS that
                                // matched lhs_detector_ID
      } else {
        // Did not find it!
        rhsWI = -1; // Marker to mean its not in the LHS.

        //            std::ostringstream mess;
        //            mess << "BinaryOperation: cannot find a RHS spectrum that
        //            contains the detectors in LHS workspace index " << lhsWI
        //            << "\n";
        //            throw std::runtime_error(mess.str());
      (*table)[lhsWI] = rhsWI;
      done = true; // Great, we did it.
    // ----------------- LHS detectors are subset of RHS, which are Grouped
    // --------------------------------------
    if (!done) {

      // Didn't find it? Now we need to iterate through the output workspace to
      //  match the detector ID.
      // NOTE: This can be SUPER SLOW!
      for (rhsWI = 0; rhsWI < static_cast<int64_t>(rhs_nhist); rhsWI++) {
        const std::set<detid_t> &rhsDets =
            rhs->getSpectrum(rhsWI)->getDetectorIDs();

        // Checks that lhsDets is a subset of rhsDets
        if (std::includes(rhsDets.begin(), rhsDets.end(), lhsDets.begin(),
                          lhsDets.end())) {
          // This one is right. Now we can stop looking.
          done = true;
          continue;
    // ------- Still nothing ! -----------
    if (!done) {
      (*table)[lhsWI] = -1;

      //          std::ostringstream mess;
      //          mess << "BinaryOperation: cannot find a RHS spectrum that
      //          contains the detectors in LHS workspace index " << lhsWI <<
      //          "\n";
      //          throw std::runtime_error(mess.str());
} // namespace Algorithms