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

std::string makeAxisName(const Kernel::V3D &Dir,
                         const std::vector<std::string> &QNames) {
  double eps(1.e-3);
  Kernel::V3D absDir(fabs(Dir.X()), fabs(Dir.Y()), fabs(Dir.Z()));
  std::string mainName;

  if ((absDir[0] >= absDir[1]) && (absDir[0] >= absDir[2])) {
    mainName = QNames[0];
  } else if (absDir[1] >= absDir[2]) {
    mainName = QNames[1];
  } else {
    mainName = QNames[2];
  }

  std::string name("["), separator = ",";
  for (size_t i = 0; i < 3; i++) {

    if (i == 2)
      separator = "]";
    if (absDir[i] < eps) {
      name += "0" + separator;
      continue;
    }
    if (Dir[i] < 0) {
      name += "-";
    }
    if (std::fabs(absDir[i] - 1) < eps) {
      name.append(mainName).append(separator);
      continue;
    }
    name.append(sprintfd(absDir[i], eps)).append(mainName).append(separator);
  }

  return name;
}

/**
 * Return the XML for a sphere.
 */
std::string sphereXML(double radius, const Kernel::V3D &centre, const std::string &id) {
    std::ostringstream xml;
    xml << "<sphere id=\"" << id << "\">"
        << "<centre x=\"" << centre.X() << "\"  y=\"" << centre.Y() << "\" z=\""
        << centre.Z() << "\" />"
        << "<radius val=\"" << radius << "\" />"
        << "</sphere>";
    return xml.str();
}

/**
 * Return a local point in a cylinder shape
 *
 * @param basis a basis vector
 * @param alongAxis symmetry axis vector of a cylinder
 * @param polarAngle a polar angle (in radians) of a point in a cylinder
 * @param radialLength radial position of point in a cylinder
 * @return a local point inside the cylinder
 */
Kernel::V3D localPointInCylinder(const Kernel::V3D &basis,
                                 const Kernel::V3D &alongAxis,
                                 double polarAngle, double radialLength) {
  // Use basis to get a second perpendicular vector to define basis2
  Kernel::V3D basis2;
  if (basis.X() == 0) {
    basis2.setX(1.);
  } else if (basis.Y() == 0) {
    basis2.setY(1.);
  } else if (basis.Z() == 0) {
    basis2.setZ(1.);
  } else {
    basis2.setX(-basis.Y());
    basis2.setY(basis.X());
    basis2.normalize();
  }
  const Kernel::V3D basis3{basis.cross_prod(basis2)};
  const Kernel::V3D localPoint{
      ((basis2 * std::cos(polarAngle) + basis3 * std::sin(polarAngle)) *
       radialLength) +
      alongAxis};
  return localPoint;
}

/**
 * Returns the symmetry axis for the given matrix
 *
 * According to ITA, 11.2 the axis of a symmetry operation can be determined by
 * solving the Eigenvalue problem \f$Wu = u\f$ for rotations or \f$Wu = -u\f$
 * for rotoinversions. This is implemented using the general real non-symmetric
 * eigen-problem solver provided by the GSL.
 *
 * @param matrix :: Matrix of a SymmetryOperation
 * @return Axis of symmetry element.
 */
V3R SymmetryElementWithAxisGenerator::determineAxis(
    const Kernel::IntMatrix &matrix) const {
  gsl_matrix *eigenMatrix = getGSLMatrix(matrix);
  gsl_matrix *identityMatrix =
      getGSLIdentityMatrix(matrix.numRows(), matrix.numCols());

  gsl_eigen_genv_workspace *eigenWs = gsl_eigen_genv_alloc(matrix.numRows());

  gsl_vector_complex *alpha = gsl_vector_complex_alloc(3);
  gsl_vector *beta = gsl_vector_alloc(3);
  gsl_matrix_complex *eigenVectors = gsl_matrix_complex_alloc(3, 3);

  gsl_eigen_genv(eigenMatrix, identityMatrix, alpha, beta, eigenVectors,
                 eigenWs);
  gsl_eigen_genv_sort(alpha, beta, eigenVectors, GSL_EIGEN_SORT_ABS_DESC);

  double determinant = matrix.determinant();

  Kernel::V3D eigenVector;

  for (size_t i = 0; i < matrix.numCols(); ++i) {
    double eigenValue = GSL_REAL(gsl_complex_div_real(
        gsl_vector_complex_get(alpha, i), gsl_vector_get(beta, i)));

    if (fabs(eigenValue - determinant) < 1e-9) {
      for (size_t j = 0; j < matrix.numRows(); ++j) {
        double element = GSL_REAL(gsl_matrix_complex_get(eigenVectors, j, i));

        eigenVector[j] = element;
      }
    }
  }

  eigenVector *= determinant;

  double sumOfElements = eigenVector.X() + eigenVector.Y() + eigenVector.Z();

  if (sumOfElements < 0) {
    eigenVector *= -1.0;
  }

  gsl_matrix_free(eigenMatrix);
  gsl_matrix_free(identityMatrix);
  gsl_eigen_genv_free(eigenWs);
  gsl_vector_complex_free(alpha);
  gsl_vector_free(beta);
  gsl_matrix_complex_free(eigenVectors);

  double min = 1.0;
  for (size_t i = 0; i < 3; ++i) {
    double absoluteValue = fabs(eigenVector[i]);
    if (absoluteValue != 0.0 &&
        (eigenVector[i] < min && (absoluteValue - fabs(min)) < 1e-9)) {
      min = eigenVector[i];
    }
  }

  V3R axis;
  for (size_t i = 0; i < 3; ++i) {
    axis[i] = static_cast<int>(boost::math::round(eigenVector[i] / min));
  }

  return axis;
}

/** method to cacluate the detectors parameters and add them to the detectors
*averages
*@param spDet    -- shared pointer to the Mantid Detector
*@param Observer -- sample position or the centre of the polar system of
*coordinates to calculate detector's parameters.
*/
void AvrgDetector::addDetInfo(const Geometry::IDetector_const_sptr &spDet,
                              const Kernel::V3D &Observer) {
  m_nComponents++;
  Kernel::V3D detPos = spDet->getPos();
  Kernel::V3D toDet = (detPos - Observer);

  double dist2Det, Polar, Azimut, ringPolar, ringAzim;
  // identify the detector' position in the beam coordinate system:
  toDet.getSpherical(dist2Det, Polar, Azimut);
  if (m_nComponents <= 1) {
    m_FlightPathSum = dist2Det;
    m_PolarSum = Polar;
    m_AzimutSum = Azimut;

    m_AzimBase = Polar;
    m_PolarBase = Azimut;
    ringPolar = Polar;
    ringAzim = Azimut;
  } else {
    ringPolar = nearAngle(m_AzimBase, Polar);
    ringAzim = nearAngle(m_PolarBase, Azimut);
    m_FlightPathSum += dist2Det;
    m_PolarSum += ringPolar;
    m_AzimutSum += ringAzim;
  }

  // centre of the azimuthal ring (the ring  detectors form around the beam)
  Kernel::V3D ringCentre(0, 0, toDet.Z());

  // Get the bounding box
  Geometry::BoundingBox bbox;
  std::vector<Kernel::V3D> coord(3);

  Kernel::V3D er(0, 1, 0), e_th,
      ez(0, 0, 1); // ez along beamline, which is always oz;
  if (dist2Det)
    er = toDet / dist2Det; // direction to the detector
  Kernel::V3D e_tg =
      er.cross_prod(ez); // tangential to the ring and anticloakwise;
  e_tg.normalize();
  // make orthogonal -- projections are calculated in this coordinate system
  ez = e_tg.cross_prod(er);

  coord[0] = er;   // new X
  coord[1] = ez;   // new y
  coord[2] = e_tg; // new z
  bbox.setBoxAlignment(ringCentre, coord);

  spDet->getBoundingBox(bbox);

  // linear extensions of the bounding box orientied tangentially to the equal
  // scattering angle circle
  double azimMin = bbox.zMin();
  double azimMax = bbox.zMax();
  double polarMin = bbox.yMin(); // bounding box has been rotated according to
                                 // coord above, so z is along e_tg
  double polarMax = bbox.yMax();

  if (m_useSphericalSizes) {
    if (dist2Det == 0)
      dist2Det = 1;

    // convert to angular units
    double polarHalfSize =
        rad2deg * atan2(0.5 * (polarMax - polarMin), dist2Det);
    double azimHalfSize = rad2deg * atan2(0.5 * (azimMax - azimMin), dist2Det);

    polarMin = ringPolar - polarHalfSize;
    polarMax = ringPolar + polarHalfSize;
    azimMin = ringAzim - azimHalfSize;
    azimMax = ringAzim + azimHalfSize;
  }
  if (m_AzimMin > azimMin)
    m_AzimMin = azimMin;
  if (m_AzimMax < azimMax)
    m_AzimMax = azimMax;

  if (m_PolarMin > polarMin)
    m_PolarMin = polarMin;
  if (m_PolarMax < polarMax)
    m_PolarMax = polarMax;
}

/** Create output workspace
 * @brief ConvertCWSDExpToMomentum::createExperimentMDWorkspace
 * @return
 */
API::IMDEventWorkspace_sptr ConvertCWSDMDtoHKL::createHKLMDWorkspace(
    const std::vector<Kernel::V3D> &vec_hkl,
    const std::vector<signal_t> &vec_signal,
    const std::vector<detid_t> &vec_detid) {
  // Check
  if (vec_hkl.size() != vec_signal.size() ||
      vec_signal.size() != vec_detid.size())
    throw std::invalid_argument("Input vectors for HKL, signal and detector "
                                "IDs are of different size!");

  // Create workspace in Q_sample with dimenion as 3
  size_t nDimension = 3;
  IMDEventWorkspace_sptr mdws =
      MDEventFactory::CreateMDWorkspace(nDimension, "MDEvent");

  // Extract Dimensions and add to the output workspace.
  std::vector<std::string> vec_ID(3);
  vec_ID[0] = "H";
  vec_ID[1] = "K";
  vec_ID[2] = "L";

  std::vector<std::string> dimensionNames(3);
  dimensionNames[0] = "H";
  dimensionNames[1] = "K";
  dimensionNames[2] = "L";

  Mantid::Kernel::SpecialCoordinateSystem coordinateSystem =
      Mantid::Kernel::HKL;

  // Add dimensions
  std::vector<double> m_extentMins(3);
  std::vector<double> m_extentMaxs(3);
  std::vector<size_t> m_numBins(3, 100);
  getRange(vec_hkl, m_extentMins, m_extentMaxs);

  // Get MDFrame of HKL type with RLU
  auto unitFactory = makeMDUnitFactoryChain();
  auto unit = unitFactory->create(Units::Symbol::RLU.ascii());
  Mantid::Geometry::HKL frame(unit);

  for (size_t i = 0; i < nDimension; ++i) {
    std::string id = vec_ID[i];
    std::string name = dimensionNames[i];
    // std::string units = "A^-1";
    mdws->addDimension(
        Geometry::MDHistoDimension_sptr(new Geometry::MDHistoDimension(
            id, name, frame, static_cast<coord_t>(m_extentMins[i]),
            static_cast<coord_t>(m_extentMaxs[i]), m_numBins[i])));
  }

  // Set coordinate system
  mdws->setCoordinateSystem(coordinateSystem);

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

  // Go though each spectrum to conver to MDEvent
  for (size_t iq = 0; iq < vec_hkl.size(); ++iq) {
    Kernel::V3D hkl = vec_hkl[iq];
    std::vector<Mantid::coord_t> millerindex(3);
    millerindex[0] = static_cast<float>(hkl.X());
    millerindex[1] = static_cast<float>(hkl.Y());
    millerindex[2] = static_cast<float>(hkl.Z());

    signal_t signal = vec_signal[iq];
    signal_t error = std::sqrt(signal);
    uint16_t runnumber = 1;
    detid_t detid = vec_detid[iq];

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

  return mdws;
}

double MeshObject2D::distanceToPlane(const Kernel::V3D &point) const {
  return ((point.X() * m_planeParameters.a) +
          (point.Y() * m_planeParameters.b) +
          (point.Z() * m_planeParameters.c) + m_planeParameters.k);
}