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

int Line::intersect(std::list<Kernel::V3D> &PntOut, const Cylinder &Cyl) const
/**
For the line that intersects the cylinder generate
add the point to the VecOut, return number of points
added. It does not check the points for validity.

@param PntOut :: Vector of points found by the line/cylinder intersection
@param Cyl :: Cylinder to intersect line with
@return Number of points found by intersection
*/
{
    const Kernel::V3D Cent = Cyl.getCentre();
    const Kernel::V3D Ax = Origin - Cent;
    const Kernel::V3D N = Cyl.getNormal();
    const double R = Cyl.getRadius();
    const double vDn = N.scalar_prod(Direct);
    const double vDA = N.scalar_prod(Ax);
    // First solve the equation of intersection
    double C[3];
    C[0] = 1.0 - (vDn * vDn);
    C[1] = 2.0 * (Ax.scalar_prod(Direct) - vDA * vDn);
    C[2] = Ax.scalar_prod(Ax) - (R * R + vDA * vDA);
    std::pair<std::complex<double>, std::complex<double>> SQ;
    const int ix = solveQuadratic(C, SQ);
    // This takes the centre displacement into account:
    return lambdaPair(ix, SQ, PntOut);
}

/**
 Calculate if the point R is on
 the cone (Note: have to be careful here
 since angle calcuation calcuates an angle.
 We need a distance for tolerance!)
 @param R :: Point to check
 @return 1 if on surface and 0 if not not on surface
 */
int Cone::onSurface(const Kernel::V3D &R) const {

    const Kernel::V3D cR = R - Centre;
    double rptAngle = cR.scalar_prod(Normal);
    rptAngle *= rptAngle / cR.scalar_prod(cR);
    const double eqn(sqrt(rptAngle));

    return (fabs(eqn - cangle) > Tolerance) ? 0 : 1;
}

    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);
    }

double Plane::distance(const Kernel::V3D &A) const
/**
  Determine the distance of point A from the plane
  returns a value relative to the normal
  @param A :: point to get distance from
  @return singed distance from point
*/
{
  return A.scalar_prod(NormV) - Dist;
}

int Plane::setPlane(const Kernel::V3D &P, const Kernel::V3D &N)
/**
  Given a point and a normal direction set the plane
  @param P :: Point for plane to pass thought
  @param N :: Normal for the plane
  @retval 0 :: success
*/
{
  try {
    NormV = normalize(N);
  } catch (std::runtime_error &) {
    throw std::invalid_argument("Attempt to create Plane with zero normal");
  }
  Dist = P.scalar_prod(NormV);
  setBaseEqn();
  return 0;
}

/**
 * Returns whether the given reflection is allowed or not in this space group
 *
 * Space groups that contain translational symmetry cause certain reflections
 * to be absent due to the contributions of symmetry equivalent atoms to the
 * structure factor cancelling out. This method implements the procedure
 * described in ITA [1] to check whether a reflection is allowed or not
 * according to the symmetry operations in the space group. Please note that
 * certain arrangements of atoms can lead to additional conditions that can not
 * be determined using a space group's symmetry operations alone. For these
 * situations, Geometry::CrystalStructure can help.
 *
 * [1] International Tables for Crystallography (2006). Vol. A, ch. 12.3, p. 832
 *
 * @param hkl :: HKL to be checked.
 * @return :: true if the reflection is allowed, false otherwise.
 */
bool SpaceGroup::isAllowedReflection(const Kernel::V3D &hkl) const {
  for (const auto &operation : m_allOperations) {
    if (operation.hasTranslation()) {
      /* Floating point precision problem:
       *    (H . v) % 1.0 is not always exactly 0, so instead:
       *    | [(H . v) + delta] % 1.0 | > 1e-14 is checked
       * The transformation is only performed if necessary.
       */
      if ((fabs(fmod(fabs(hkl.scalar_prod(operation.reducedVector())) + 1e-15,
                     1.0)) > 1e-14) &&
          (operation.transformHKL(hkl) == hkl)) {
        return false;
      }
    }
  }

  return true;
}

double Cylinder::distance(const Kernel::V3D &A) const
/**
 Calculates the distance of point A from
 the surface of the  cylinder.
 @param A :: Point to calculate distance from
 @return :: +ve distance to the surface.

 \todo INCOMPLETE AS Does not deal with the case of
 non axis centre line  (has been updated?? )
 */
{
  // First find the normal going to the point
  const Kernel::V3D Amov = A - Centre;
  double lambda = Amov.scalar_prod(Normal);
  const Kernel::V3D Ccut = Normal * lambda;
  // The distance is from the centre line to the
  return fabs(Ccut.distance(Amov) - Radius);
}

/**
 * This function calculates the exponential contribution to the He3 tube
 * efficiency.
 * @param spectraIndex :: the current index to calculate
 * @param idet :: the current detector pointer
 * @throw out_of_range if twice tube thickness is greater than tube diameter
 * @return the exponential contribution for the given detector
 */
double He3TubeEfficiency::calculateExponential(
    std::size_t spectraIndex,
    boost::shared_ptr<const Geometry::IDetector> idet) {
  // Get the parameters for the current associated tube
  double pressure =
      this->getParameter("TubePressure", spectraIndex, "tube_pressure", idet);
  double tubethickness =
      this->getParameter("TubeThickness", spectraIndex, "tube_thickness", idet);
  double temperature = this->getParameter("TubeTemperature", spectraIndex,
                                          "tube_temperature", idet);

  double detRadius(0.0);
  Kernel::V3D detAxis;
  this->getDetectorGeometry(idet, detRadius, detAxis);
  double detDiameter = 2.0 * detRadius;
  double twiceTubeThickness = 2.0 * tubethickness;

  // now get the sin of the angle, it's the magnitude of the cross product of
  // unit vector along the detector tube axis and a unit vector directed from
  // the sample to the detector center
  Kernel::V3D vectorFromSample = idet->getPos() - this->samplePos;
  vectorFromSample.normalize();
  Kernel::Quat rot = idet->getRotation();
  // rotate the original cylinder object axis to get the detector axis in the
  // actual instrument
  rot.rotate(detAxis);
  detAxis.normalize();
  // Scalar product is quicker than cross product
  double cosTheta = detAxis.scalar_prod(vectorFromSample);
  double sinTheta = std::sqrt(1.0 - cosTheta * cosTheta);

  const double straight_path = detDiameter - twiceTubeThickness;
  if (std::fabs(straight_path - 0.0) < TOL) {
    throw std::out_of_range("Twice tube thickness cannot be greater than "
                            "or equal to the tube diameter");
  }

  const double pathlength = straight_path / sinTheta;
  return EXP_SCALAR_CONST * (pressure / temperature) * pathlength;
}

/**
 * Calculate volume.
 * @return The volume.
 */
double MeshObject::volume() const {
  // Select centre of bounding box as centre point.
  // For each triangle calculate the signed volume of
  // the tetrahedron formed by the triangle and the
  // centre point. Then add to total.

  BoundingBox bb = getBoundingBox();
  double cX = 0.5 * (bb.xMax() + bb.xMin());
  double cY = 0.5 * (bb.yMax() + bb.yMin());
  double cZ = 0.5 * (bb.zMax() + bb.zMin());
  Kernel::V3D centre(cX, cY, cZ);

  double volumeTimesSix(0.0);

  Kernel::V3D vertex1, vertex2, vertex3;
  for (size_t i = 0; getTriangle(i, vertex1, vertex2, vertex3); ++i) {
    Kernel::V3D a = vertex1 - centre;
    Kernel::V3D b = vertex2 - centre;
    Kernel::V3D c = vertex3 - centre;
    volumeTimesSix += a.scalar_prod(b.cross_prod(c));
  }

  return volumeTimesSix / 6.0;
}

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

int Plane::setSurface(const std::string &Pstr)
/**
   processes a standard MCNPX plane string:
   There are three types :
   - (A) px Distance
   - (B) p A B C D (equation Ax+By+Cz=D)
   - (C) p V3D V3D V3D
   @param Pstr :: String to make into a plane of type p{xyz} or p
   @return 0 on success, -ve of failure
*/
{
  // Two types of plane string p[x-z]  and p
  std::string Line = Pstr;
  std::string item;

  if (!Mantid::Kernel::Strings::section(Line, item) || tolower(item[0]) != 'p')
    return -1;
  // Only 3 need to be declared
  double surf[9] = {0.0, 0, 0, 0, 0};

  if (item.size() == 1) // PROCESS BASIC PLANE:
  {
    int cnt;
    for (cnt = 0; cnt < 9 && Mantid::Kernel::Strings::section(Line, surf[cnt]);
         cnt++)
      ;

    if (cnt != 4 && cnt != 9)
      return -3;
    if (cnt == 9) // V3D type
    {
      Kernel::V3D A = Kernel::V3D(surf[0], surf[1], surf[2]);
      Kernel::V3D B = Kernel::V3D(surf[3], surf[4], surf[5]);
      Kernel::V3D C = Kernel::V3D(surf[6], surf[7], surf[8]);
      B -= A;
      C -= A;
      NormV = B * C;
      NormV.normalize();
      Dist = A.scalar_prod(NormV);
    } else // Norm Equation:
    {
      NormV = Kernel::V3D(surf[0], surf[1], surf[2]);
      const double ll = NormV.normalize();
      if (ll < Tolerance) // avoid
        return -4;
      Dist = surf[3] / ll;
    }
  } else if (item.size() == 2) //  PROCESS px type PLANE
  {
    const int ptype = static_cast<int>(tolower(item[1]) - 'x');
    if (ptype < 0 || ptype > 2) // Not x,y,z
      return -5;
    surf[ptype] = 1.0;
    if (!Mantid::Kernel::Strings::convert(Line, Dist))
      return -6; // Too short or no number
    NormV = Kernel::V3D(surf[0], surf[1], surf[2]);
  } else
    return -3; // WRONG NAME

  setBaseEqn();
  return 0;
}