C++ V3D::norm方法代码示例

    double Cone::distance(const Kernel::V3D& Pt) const
    /**
     Calculates the distance from the point to the Cone
     does not calculate the point on the cone that is closest
     @param Pt :: Point to calcuate from

     - normalise to a cone vertex at the origin
     - calculate the angle between the axis and the Point
     - Calculate the distance to P
     @return distance to Pt
     */
    {
      const Kernel::V3D Px = Pt - Centre;
      // test is the centre to point distance is zero
      if (Px.norm() < Tolerance)
        return Px.norm();
      double Pangle = Px.scalar_prod(Normal) / Px.norm();
      if (Pangle < 0.0)
        Pangle = acos(-Pangle);
      else
        Pangle = acos(Pangle);

      Pangle -= M_PI * alpha / 180.0;
      return Px.norm() * sin(Pangle);
    }

/**
 * Returns selected information for a "peak" at QLabFrame.
 *
 * @param qFrame      An arbitrary position in Q-space.  This does not have to
 *be the
 *                    position of a peak.
 * @param labCoords  Set true if the position is in the lab coordinate system,
 *false if
 *                    it is in the sample coordinate system.
 * @return a vector whose elements contain different information about the
 *"peak" at that position.
 *         each element is a pair of description of information and the string
 *form for the corresponding
 *         value.
 */
int PeaksWorkspace::peakInfoNumber(Kernel::V3D qFrame, bool labCoords) const {
  std::vector<std::pair<std::string, std::string>> Result;
  std::ostringstream oss;
  oss << std::setw(12) << std::fixed << std::setprecision(3) << (qFrame.norm());
  std::pair<std::string, std::string> QMag("|Q|", oss.str());
  Result.push_back(QMag);

  oss.str("");
  oss.clear();
  oss << std::setw(12) << std::fixed << std::setprecision(3)
      << (2.0 * M_PI / qFrame.norm());

  std::pair<std::string, std::string> dspc("d-spacing", oss.str());
  oss.str("");
  oss.clear();
  Result.push_back(dspc);

  int seqNum = -1;
  double minDist = 10000000;

  for (int i = 0; i < getNumberPeaks(); i++) {
    Peak pk = getPeak(i);
    V3D Q = pk.getQLabFrame();
    if (!labCoords)
      Q = pk.getQSampleFrame();
    double D = qFrame.distance(Q);
    if (D < minDist) {
      minDist = D;
      seqNum = i + 1;
    }
  }
  return seqNum;
}

double
Torus::distance(const Kernel::V3D& Pt) const
  /**
    Calculates the distance from the point to the Torus
    does not calculate the point on the Torusthat is closest
    @param Pt :: Point to calcuate from

    - normalise to a cone vertex at the origin
    - calculate the angle between the axis and the Point
    - Calculate the distance to P
    @return distance to Pt
  */
{
  const Kernel::V3D Px=Pt-Centre;
  // test is the centre to point distance is zero
  if(Px.norm()<Tolerance)
    return Px.norm();
  return Px.norm();
}

double Line::distance(const Kernel::V3D &A) const
/**
Distance of a point from the line
@param A :: test Point
@return absolute distance (not signed)
*/
{
    const double lambda = Direct.scalar_prod(A - Origin);
    Kernel::V3D L = getPoint(lambda);
    L -= A;
    return L.norm();
}

void Cone::setNorm(const Kernel::V3D &A)
/**
 Sets the Normal and the Base Equation
 @param A :: New Normal direction
 */
{
    if (A.norm() > Tolerance) {
        Normal = A;
        Normal.normalize();
        setBaseEqn();
    }
    return;
}

/**
 * Find maximum angular half width of the bounding box from the observer, that
 * is
 * the greatest angle between the centre point and any corner point
 * @param observer :: Viewing point
 * @returns The value of the angular half-width
*/
double BoundingBox::angularWidth(const Kernel::V3D &observer) const {
  Kernel::V3D centre = centrePoint() - observer;
  std::vector<Kernel::V3D> pts;
  this->getFullBox(pts, observer);

  std::vector<Kernel::V3D>::const_iterator ip;
  double centre_norm_inv = 1.0 / centre.norm();
  double thetaMax(-1.0);
  for (ip = pts.begin(); ip != pts.end(); ++ip) {
    double theta = acos(ip->scalar_prod(centre) * centre_norm_inv / ip->norm());
    if (theta > thetaMax)
      thetaMax = theta;
  }
  return thetaMax;
}

/**
 * Establish if all vertices are coplanar
 * @param vertices : All vertices to check
 * @param normal : Surface normal
 * @return True only if all vertices are coplanar
 */
bool allCoplanar(const std::vector<Kernel::V3D> &vertices,
                 const Kernel::V3D &normal) {
  bool in_plane = true;
  auto v0 = vertices[0];
  const auto nx = normal[0];
  const auto ny = normal[1];
  const auto nz = normal[2];
  const auto k = nx * v0.X() + ny * v0.Y() + nz * v0.Z();
  const auto denom = normal.norm();
  const static double tolerance =
      1e-9; // Fixed Tolerance. Too expensive to calculate
            // based on machine uncertaintly for each
            // vertex.

  for (const auto &vertex : vertices) {
    auto d = (nx * vertex.X() + ny * vertex.Y() + nz * vertex.Z() - k) / denom;
    if (d > tolerance || d < -1 * tolerance) {
      in_plane = false;
      break;
    }
  }
  return in_plane;
}

//.........这里部分代码省略.........
                        << "\n";
    // Move in cm for small shifts
    g_log.information() << "Move (X)   = " << gsl_vector_get(s->x, 0) * 0.01
                        << "  \n";
    g_log.information() << "Move (Y)   = " << gsl_vector_get(s->x, 1) * 0.01
                        << "  \n";
    g_log.information() << "Move (Z)   = " << gsl_vector_get(s->x, 2) * 0.01
                        << "  \n";
    g_log.information() << "Rotate (X) = " << gsl_vector_get(s->x, 3) << "  \n";
    g_log.information() << "Rotate (Y) = " << gsl_vector_get(s->x, 4) << "  \n";
    g_log.information() << "Rotate (Z) = " << gsl_vector_get(s->x, 5) << "  \n";

    Kernel::V3D CalCenter =
        V3D(gsl_vector_get(s->x, 0) * 0.01, gsl_vector_get(s->x, 1) * 0.01,
            gsl_vector_get(s->x, 2) * 0.01);
    Kernel::V3D Center = detList[det]->getPos() + CalCenter;
    int pixmax = detList[det]->xpixels() - 1;
    int pixmid = (detList[det]->ypixels() - 1) / 2;
    BoundingBox box;
    detList[det]->getAtXY(pixmax, pixmid)->getBoundingBox(box);
    double baseX = box.xMax();
    double baseY = box.yMax();
    double baseZ = box.zMax();
    Kernel::V3D Base = V3D(baseX, baseY, baseZ) + CalCenter;
    pixmid = (detList[det]->xpixels() - 1) / 2;
    pixmax = detList[det]->ypixels() - 1;
    detList[det]->getAtXY(pixmid, pixmax)->getBoundingBox(box);
    double upX = box.xMax();
    double upY = box.yMax();
    double upZ = box.zMax();
    Kernel::V3D Up = V3D(upX, upY, upZ) + CalCenter;
    Base -= Center;
    Up -= Center;
    // Rotate around x
    baseX = Base[0];
    baseY = Base[1];
    baseZ = Base[2];
    double deg2rad = M_PI / 180.0;
    double angle = gsl_vector_get(s->x, 3) * deg2rad;
    Base = V3D(baseX, baseY * cos(angle) - baseZ * sin(angle),
               baseY * sin(angle) + baseZ * cos(angle));
    upX = Up[0];
    upY = Up[1];
    upZ = Up[2];
    Up = V3D(upX, upY * cos(angle) - upZ * sin(angle),
             upY * sin(angle) + upZ * cos(angle));
    // Rotate around y
    baseX = Base[0];
    baseY = Base[1];
    baseZ = Base[2];
    angle = gsl_vector_get(s->x, 4) * deg2rad;
    Base = V3D(baseZ * sin(angle) + baseX * cos(angle), baseY,
               baseZ * cos(angle) - baseX * sin(angle));
    upX = Up[0];
    upY = Up[1];
    upZ = Up[2];
    Up = V3D(upZ * cos(angle) - upX * sin(angle), upY,
             upZ * sin(angle) + upX * cos(angle));
    // Rotate around z
    baseX = Base[0];
    baseY = Base[1];
    baseZ = Base[2];
    angle = gsl_vector_get(s->x, 5) * deg2rad;
    Base = V3D(baseX * cos(angle) - baseY * sin(angle),
               baseX * sin(angle) + baseY * cos(angle), baseZ);
    upX = Up[0];
    upY = Up[1];
    upZ = Up[2];
    Up = V3D(upX * cos(angle) - upY * sin(angle),
             upX * sin(angle) + upY * cos(angle), upZ);
    Base.normalize();
    Up.normalize();
    Center *= 100.0;
    // << det+1  << "  "
    outfile << "5  " << detList[det]->getName().substr(4) << "  "
            << detList[det]->xpixels() << "  " << detList[det]->ypixels()
            << "  " << 100.0 * detList[det]->xsize() << "  "
            << 100.0 * detList[det]->ysize() << "  "
            << "0.2000"
            << "  " << Center.norm() << "  ";
    Center.write(outfile);
    outfile << "  ";
    Base.write(outfile);
    outfile << "  ";
    Up.write(outfile);
    outfile << "\n";

    // clean up dynamically allocated gsl stuff
    gsl_vector_free(x);
    gsl_vector_free(ss);
    gsl_multimin_fminimizer_free(s);

    // Remove the now-unneeded grouping workspace
    AnalysisDataService::Instance().remove(groupWSName);
    prog.report(detList[det]->getName());
  }

  // Closing
  outfile.close();
}

/** Build meaningful dimension names for different conversion modes
 * @param TargWSDescription the class-container to keep the dimension names and
 dimension unints
 * @param FrameID -- the ID describing the target transformation frame (lab,
 sample, hkl)
 * @param ScaleID -- the scale ID which define how the dimensions are scaled

*/
void MDWSTransform::setQ3DDimensionsNames(
    MDWSDescription &TargWSDescription, CnvrtToMD::TargetFrame FrameID,
    CnvrtToMD::CoordScaling ScaleID) const {

  std::vector<Kernel::V3D> dimDirections;
  // set default dimension names:
  std::vector<std::string> dimNames = TargWSDescription.getDimNames();

  // define B-matrix and Lattice parameters to one in case if no OrientedLattice
  // is there
  Kernel::DblMatrix Bm(3, 3, true);
  std::vector<double> LatPar(3, 1);
  if (TargWSDescription.hasLattice()) { // redefine B-matrix and Lattice
                                        // parameters from real oriented lattice
                                        // if there is one
    auto spLatt = TargWSDescription.getLattice();
    Bm = spLatt->getB();
    for (int i = 0; i < 3; i++)
      LatPar[i] = spLatt->a(i);
  }
  if (FrameID == CnvrtToMD::AutoSelect)
    FrameID = findTargetFrame(TargWSDescription);

  switch (FrameID) {
  case (CnvrtToMD::LabFrame): {
    dimNames[0] = "Q_lab_x";
    dimNames[1] = "Q_lab_y";
    dimNames[2] = "Q_lab_z";
    TargWSDescription.setCoordinateSystem(Mantid::Kernel::QLab);
    TargWSDescription.setFrame(Geometry::QLab::QLabName);
    break;
  }
  case (CnvrtToMD::SampleFrame): {
    dimNames[0] = "Q_sample_x";
    dimNames[1] = "Q_sample_y";
    dimNames[2] = "Q_sample_z";
    TargWSDescription.setCoordinateSystem(Mantid::Kernel::QSample);
    TargWSDescription.setFrame(Geometry::QSample::QSampleName);
    break;
  }
  case (CnvrtToMD::HKLFrame): {
    dimNames[0] = "H";
    dimNames[1] = "K";
    dimNames[2] = "L";

    Kernel::MDUnit_uptr mdUnit(new Kernel::InverseAngstromsUnit);
    TargWSDescription.setCoordinateSystem(Mantid::Kernel::HKL);
    TargWSDescription.setFrame(Geometry::HKL::HKLName);
    break;
  }
  default:
    throw(std::invalid_argument(" Unknown or undefined Target Frame ID"));
  }

  dimDirections.resize(3);
  dimDirections[0] = m_UProj;
  dimDirections[1] = m_VProj;
  dimDirections[2] = m_WProj;
  if (ScaleID == OrthogonalHKLScale) {
    std::vector<Kernel::V3D> uv(2);
    uv[0] = m_UProj;
    uv[1] = m_VProj;
    dimDirections = Kernel::V3D::makeVectorsOrthogonal(uv);
  }
  // axis names:
  if ((FrameID == CnvrtToMD::LabFrame) || (FrameID == CnvrtToMD::SampleFrame))
    for (int i = 0; i < 3; i++)
      TargWSDescription.setDimName(i, dimNames[i]);
  else
    for (int i = 0; i < 3; i++)
      TargWSDescription.setDimName(
          i, MDAlgorithms::makeAxisName(dimDirections[i], dimNames));

  if (ScaleID == NoScaling) {
    for (int i = 0; i < 3; i++)
      TargWSDescription.setDimUnit(i, "A^-1");
  }
  if (ScaleID == SingleScale) {
    double dMax(-1.e+32);
    for (int i = 0; i < 3; i++)
      dMax = (dMax > LatPar[i]) ? (dMax) : (LatPar[i]);
    for (int i = 0; i < 3; i++)
      TargWSDescription.setDimUnit(
          i, "in " + MDAlgorithms::sprintfd(2 * M_PI / dMax, 1.e-3) + " A^-1");
  }
  if ((ScaleID == OrthogonalHKLScale) || (ScaleID == HKLScale)) {
    // get the length along each of the axes
    std::vector<double> len;
    Kernel::V3D x;
    x = Bm * dimDirections[0];
    len.push_back(2 * M_PI * x.norm());
    x = Bm * dimDirections[1];
//.........这里部分代码省略.........

/**
 * Make a map of the conversion factors between tof and D-spacing
 * for all pixel IDs in a workspace.
 * map vulcan should contain the module/module and stack/stack offset
 *
 * @param vulcan :: map between detector ID and vulcan correction factor.
 * @param offsetsWS :: OffsetsWorkspace to be filled.
 */
void LoadDspacemap::CalculateOffsetsFromVulcanFactors(
    std::map<detid_t, double> &vulcan,
    Mantid::DataObjects::OffsetsWorkspace_sptr offsetsWS) {
  // Get a pointer to the instrument contained in the workspace
  // At this point, instrument VULCAN has been created?
  Instrument_const_sptr instrument = offsetsWS->getInstrument();

  g_log.notice() << "Name of instrument = " << instrument->getName()
                 << std::endl;
  g_log.notice() << "Input map (dict):  size = " << vulcan.size() << std::endl;

  // To get all the detector ID's
  detid2det_map allDetectors;
  instrument->getDetectors(allDetectors);

  detid2det_map::const_iterator it;
  int numfinds = 0;
  g_log.notice() << "Input number of detectors = " << allDetectors.size()
                 << std::endl;

  // Get detector information
  double l1, beamline_norm;
  Kernel::V3D beamline, samplePos;
  instrument->getInstrumentParameters(l1, beamline, beamline_norm, samplePos);

  /*** A survey of parent detector
  std::map<detid_t, bool> parents;
  for (it = allDetectors.begin(); it != allDetectors.end(); it++){
    int32_t detid = it->first;

    // def boost::shared_ptr<const Mantid::Geometry::IDetector>
  IDetector_const_sptr;

    std::string parentname =
  it->second->getParent()->getComponentID()->getName();
    g_log.notice() << "Name = " << parentname << std::endl;
    // parents.insert(parentid, true);
  }
  ***/

  /*** Here some special configuration for VULCAN is hard-coded here!
   *   Including (1) Super-Parent Information
   ***/
  Kernel::V3D referencePos;
  detid_t anydetinrefmodule = 21 * 1250 + 5;

  std::map<detid_t, Geometry::IDetector_const_sptr>::iterator det_iter =
      allDetectors.find(anydetinrefmodule);

  if (det_iter == allDetectors.end()) {
    throw std::invalid_argument("Any Detector ID is Instrument's detector");
  }
  referencePos = det_iter->second->getParent()->getPos();
  double refl2 = referencePos.norm();
  double halfcosTwoThetaRef =
      referencePos.scalar_prod(beamline) / (refl2 * beamline_norm);
  double sinThetaRef = sqrt(0.5 - halfcosTwoThetaRef);
  double difcRef = sinThetaRef * (l1 + refl2) / CONSTANT;

  // Loop over all detectors in instrument to find the offset
  for (it = allDetectors.begin(); it != allDetectors.end(); ++it) {
    int detectorID = it->first;
    Geometry::IDetector_const_sptr det = it->second;
    double offset = 0.0;

    // Find the vulcan factor;
    double vulcan_factor = 0.0;
    std::map<detid_t, double>::const_iterator vulcan_iter =
        vulcan.find(detectorID);
    if (vulcan_iter != vulcan.end()) {
      vulcan_factor = vulcan_iter->second;
      numfinds++;
    }

    // g_log.notice() << "Selected Detector with ID = " << detectorID << "  ID2
    // = " << id2 << std::endl; proved to be same

    double intermoduleoffset = 0;
    double interstackoffset = 0;

    detid_t intermoduleid = detid_t(detectorID / 1250) * 1250 + 1250 - 2;
    vulcan_iter = vulcan.find(intermoduleid);
    if (vulcan_iter == vulcan.end()) {
      g_log.error() << "Cannot find inter-module offset ID = " << intermoduleid
                    << std::endl;
    } else {
      intermoduleoffset = vulcan_iter->second;
    }

    detid_t interstackid = detid_t(detectorID / 1250) * 1250 + 1250 - 1;
    vulcan_iter = vulcan.find(interstackid);
    if (vulcan_iter == vulcan.end()) {
//.........这里部分代码省略.........

/** Read the scaling information from a file (e.g. merlin_detector.sca) or from
 * the RAW file (.raw)
 *  @param scalingFile :: Name of scaling file .sca
 *  @param truepos :: V3D vector of actual positions as read from the file
 *  @return False if unable to open file, True otherwise
 */
bool SetScalingPSD::processScalingFile(const std::string &scalingFile,
                                       std::vector<Kernel::V3D> &truepos) {
  // Read the scaling information from a text file (.sca extension) or from a
  // raw file (.raw)
  // This is really corrected positions as (r,theta,phi) for each detector
  // Compare these with the instrument values to determine the change in
  // position and the scaling
  // which may be necessary for each pixel if in a tube.
  // movePos is used to updated positions
  std::map<int, Kernel::V3D> posMap;
  std::map<int, double> scaleMap;
  std::map<int, double>::iterator its;

  Instrument_const_sptr instrument = m_workspace->getInstrument();
  if (scalingFile.find(".sca") != std::string::npos ||
      scalingFile.find(".SCA") != std::string::npos) {
    // read a .sca text format file
    // format consists of a short header followed by one line per detector

    std::ifstream sFile(scalingFile.c_str());
    if (!sFile) {
      g_log.error() << "Unable to open scaling file " << scalingFile
                    << std::endl;
      return false;
    }
    std::string str;
    getline(sFile,
            str); // skip header line should be <filename> generated by <prog>
    int detectorCount;
    getline(sFile, str); // get detector count line
    std::istringstream istr(str);
    istr >> detectorCount;
    if (detectorCount < 1) {
      g_log.error("Bad detector count in scaling file");
      throw std::runtime_error("Bad detector count in scaling file");
    }
    truepos.reserve(detectorCount);
    getline(sFile, str); // skip title line
    int detIdLast = -10;
    Kernel::V3D truPosLast, detPosLast;

    Progress prog(this, 0.0, 0.5, detectorCount);
    // Now loop through lines, one for each detector/monitor. The latter are
    // ignored.

    while (getline(sFile, str)) {
      if (str.empty() || str[0] == '#')
        continue;
      std::istringstream istr(str);

      // read 6 values from the line to get the 3 (l2,theta,phi) of interest
      int detIndex, code;
      double l2, theta, phi, offset;
      istr >> detIndex >> offset >> l2 >> code >> theta >> phi;

      // sanity check on angles - l2 should be +ve but sample file has a few -ve
      // values
      // on monitors
      if (theta > 181.0 || theta < -1 || phi < -181 || phi > 181) {
        g_log.error("Position angle data out of range in .sca file");
        throw std::runtime_error(
            "Position angle data out of range in .sca file");
      }
      Kernel::V3D truPos;
      // use abs as correction file has -ve l2 for first few detectors
      truPos.spherical(fabs(l2), theta, phi);
      truepos.push_back(truPos);
      //
      Geometry::IDetector_const_sptr det;
      try {
        det = instrument->getDetector(detIndex);
      } catch (Kernel::Exception::NotFoundError &) {
        continue;
      }
      Kernel::V3D detPos = det->getPos();
      Kernel::V3D shift = truPos - detPos;

      // scaling applied to dets that are not monitors and have sequential IDs
      if (detIdLast == detIndex - 1 && !det->isMonitor()) {
        Kernel::V3D diffI = detPos - detPosLast;
        Kernel::V3D diffT = truPos - truPosLast;
        double scale = diffT.norm() / diffI.norm();
        Kernel::V3D scaleDir = diffT / diffT.norm();
        // Wish to store the scaling in a map, if we already have a scaling
        // for this detector (i.e. from the other side) we average the two
        // values. End of tube detectors only have one scaling estimate.
        scaleMap[detIndex] = scale;
        its = scaleMap.find(detIndex - 1);
        if (its == scaleMap.end())
          scaleMap[detIndex - 1] = scale;
        else
          its->second = 0.5 * (its->second + scale);
        // std::cout << detIndex << scale << scaleDir << std::endl;
      }
//.........这里部分代码省略.........

/** Convert a SPICE 2D Det MatrixWorkspace to MDEvents and append to an
 * MDEventWorkspace
 * It is optional to use a virtual instrument or copy from input data workspace
 * @brief ConvertCWSDExpToMomentum::convertSpiceMatrixToMomentumMDEvents
 * @param dataws :: data matrix workspace
 * @param usevirtual :: boolean flag to use virtual instrument
 * @param startdetid :: starting detid for detectors from this workspace mapping
 * to virtual instrument in MDEventWorkspace
 * @param runnumber :: run number for all MDEvents created from this matrix
 * workspace
 */
void ConvertCWSDExpToMomentum::convertSpiceMatrixToMomentumMDEvents(
    MatrixWorkspace_sptr dataws, bool usevirtual, const detid_t &startdetid,
    const int runnumber) {
  // Create transformation matrix from which the transformation is
  Kernel::DblMatrix rotationMatrix;
  setupTransferMatrix(dataws, rotationMatrix);

  g_log.information() << "Before insert new event, output workspace has "
                      << m_outputWS->getNEvents() << "Events.\n";

  // Creates a new instance of the MDEventInserte to output workspace
  MDEventWorkspace<MDEvent<3>, 3>::sptr mdws_mdevt_3 =
      boost::dynamic_pointer_cast<MDEventWorkspace<MDEvent<3>, 3>>(m_outputWS);
  MDEventInserter<MDEventWorkspace<MDEvent<3>, 3>::sptr> inserter(mdws_mdevt_3);

  // Calcualte k_i: it is assumed that all k_i are same for one Pt.
  // number, i.e., one 2D XML file
  Kernel::V3D sourcePos = dataws->getInstrument()->getSource()->getPos();
  Kernel::V3D samplePos = dataws->getInstrument()->getSample()->getPos();
  if (dataws->readX(0).size() != 2)
    throw std::runtime_error(
        "Input matrix workspace has wrong dimension in X-axis.");
  double momentum = 0.5 * (dataws->readX(0)[0] + dataws->readX(0)[1]);
  Kernel::V3D ki = (samplePos - sourcePos) * (momentum / sourcePos.norm());

  g_log.debug() << "Source at " << sourcePos.toString()
                << ", Norm = " << sourcePos.norm()
                << ", momentum = " << momentum << "\n"
                << "k_i = " << ki.toString() << "\n";

  // Go though each spectrum to conver to MDEvent
  size_t numspec = dataws->getNumberHistograms();
  double maxsignal = 0;
  size_t nummdevents = 0;
  for (size_t iws = 0; iws < numspec; ++iws) {
    // Get detector positions and signal
    double signal = dataws->readY(iws)[0];
    // Skip event with 0 signal
    if (signal < 0.001)
      continue;
    double error = dataws->readE(iws)[0];
    Kernel::V3D detpos = dataws->getDetector(iws)->getPos();
    std::vector<Mantid::coord_t> q_sample(3);

    // Calculate Q-sample and new detector ID in virtual instrument.
    Kernel::V3D qlab = convertToQSample(samplePos, ki, detpos, momentum,
                                        q_sample, rotationMatrix);
    detid_t native_detid = dataws->getDetector(iws)->getID();
    detid_t detid = native_detid + startdetid;

    // Insert
    inserter.insertMDEvent(
        static_cast<float>(signal), static_cast<float>(error * error),
        static_cast<uint16_t>(runnumber), detid, q_sample.data());
    updateQRange(q_sample);

    g_log.debug() << "Q-lab = " << qlab.toString() << "\n";
    g_log.debug() << "Insert DetID " << detid << ", signal = " << signal
                  << ", with q_sample = " << q_sample[0] << ", " << q_sample[1]
                  << ", " << q_sample[2] << "\n";

    // Update some statistical inforamtion
    if (signal > maxsignal)
      maxsignal = signal;
    ++nummdevents;
  }

  g_log.information() << "Imported Matrixworkspace: Max. Signal = " << maxsignal
                      << ", Add " << nummdevents << " MDEvents "
                      << "\n";

  // Add experiment info including instrument, goniometer and run number
  ExperimentInfo_sptr expinfo = boost::make_shared<ExperimentInfo>();
  if (usevirtual)
    expinfo->setInstrument(m_virtualInstrument);
  else {
    Geometry::Instrument_const_sptr tmp_inst = dataws->getInstrument();
    expinfo->setInstrument(tmp_inst);
  }
  expinfo->mutableRun().setGoniometer(dataws->run().getGoniometer(), false);
  expinfo->mutableRun().addProperty("run_number", runnumber);
  // Add all the other propertys from original data workspace
  const std::vector<Kernel::Property *> vec_property =
      dataws->run().getProperties();
  for (auto property : vec_property) {
    expinfo->mutableRun().addProperty(property->clone());
  }

  m_outputWS->addExperimentInfo(expinfo);
//.........这里部分代码省略.........