Core/src/DDRotation.cc

0001 #include <cstdio>
0002 #include <atomic>
0003 #include <cmath>
0004 #include <sstream>
0005 #include <string>
0006
0007 #include "DetectorDescription/Core/interface/DDRotationMatrix.h"
0008 #include "DetectorDescription/Core/interface/DDTranslation.h"
0009 #include "DetectorDescription/Core/interface/DDBase.h"
0010 #include "DetectorDescription/Core/interface/DDName.h"
0011 #include "DetectorDescription/Core/interface/DDTransform.h"
0012 #include "DataFormats/Math/interface/GeantUnits.h"
0013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0014 #include "FWCore/Utilities/interface/Exception.h"
0015 #include "Math/GenVector/AxisAngle.h"
0016 #include "Math/GenVector/Cartesian3D.h"
0017 #include "Math/GenVector/DisplacementVector3D.h"
0018
0019 using namespace geant_units::operators;
0020
0021 std::ostream& operator<<(std::ostream& os, const DDRotation& r) {
0022   DDBase<DDName, DDRotationMatrix*>::def_type defined(r.isDefined());
0023   if (defined.first) {
0024     os << *(defined.first) << " ";
0025     if (defined.second) {
0026       const DDRotationMatrix& rm = r.rotation();
0027       DDAxisAngle ra(rm);
0028       os << "t=" << convertRadToDeg(ra.Axis().Theta()) << "deg "
0029          << "p=" << convertRadToDeg(ra.Axis().Phi()) << "deg "
0030          << "a=" << convertRadToDeg(ra.Angle()) << "deg";
0031     } else {
0032       os << "* rotation not defined * ";
0033     }
0034   } else {
0035     os << "* rotation not declared * ";
0036   }
0037   return os;
0038 }
0039
0040 DDRotation::DDRotation() : DDBase<DDName, std::unique_ptr<DDRotationMatrix>>() {
0041   constexpr char const* baseName = "DdBlNa";
0042   // In this particular case, we do not really care about multiple threads
0043   // using the same counter, we simply need to have a unique id for the
0044   // blank matrix being created, so just making this static an atomic should do
0045   // the trick. In order to ensure repeatibility one should also include some
0046   // some run specific Id, I guess. Not sure it really matters.
0047   static std::atomic<int> countBlank;
0048   char buf[64];
0049   snprintf(buf, 64, "%s%i", baseName, countBlank++);
0050   create(DDName(buf, baseName), std::make_unique<DDRotationMatrix>());
0051 }
0052
0053 DDRotation::DDRotation(const DDName& name) : DDBase<DDName, std::unique_ptr<DDRotationMatrix>>() { create(name); }
0054
0055 DDRotation::DDRotation(const DDName& name, std::unique_ptr<DDRotationMatrix> rot)
0056     : DDBase<DDName, std::unique_ptr<DDRotationMatrix>>() {
0057   create(name, std::move(rot));
0058 }
0059
0060 DDRotation::DDRotation(std::unique_ptr<DDRotationMatrix> rot) : DDBase<DDName, std::unique_ptr<DDRotationMatrix>>() {
0061   static std::atomic<int> countNN;
0062   char buf[64];
0063   snprintf(buf, 64, "DdNoNa%i", countNN++);
0064   create(DDName(buf, "DdNoNa"), std::move(rot));
0065 }
0066
0067 DDRotation DDrot(const DDName& ddname, std::unique_ptr<DDRotationMatrix> rot) {
0068   // memory of rot goes sto DDRotationImpl!!
0069   return DDRotation(ddname, std::move(rot));
0070 }
0071
0072 std::unique_ptr<DDRotation> DDrotPtr(const DDName& ddname, std::unique_ptr<DDRotationMatrix> rot) {
0073   // memory of rot goes sto DDRotationImpl!!
0074   return std::make_unique<DDRotation>(ddname, std::move(rot));
0075 }
0076
0077 // makes sure that the DDRotationMatrix constructed is right-handed and orthogonal.
0078 DDRotation DDrot(
0079     const DDName& ddname, double thetaX, double phiX, double thetaY, double phiY, double thetaZ, double phiZ) {
0080   // define 3 unit std::vectors
0081   DD3Vector x(cos(phiX) * sin(thetaX), sin(phiX) * sin(thetaX), cos(thetaX));
0082   DD3Vector y(cos(phiY) * sin(thetaY), sin(phiY) * sin(thetaY), cos(thetaY));
0083   DD3Vector z(cos(phiZ) * sin(thetaZ), sin(phiZ) * sin(thetaZ), cos(thetaZ));
0084
0085   double tol = 1.0e-3;                 // Geant4 compatible
0086   double check = (x.Cross(y)).Dot(z);  // in case of a LEFT-handed orthogonal system this must be -1
0087   if (fabs(1. - check) > tol) {
0088     edm::LogError("DDRotation") << ddname << " is not a RIGHT-handed orthonormal matrix!" << std::endl;
0089     throw cms::Exception("DDException") << ddname.name() << " is not RIGHT-handed!";
0090   }
0091
0092   return DDRotation(ddname,
0093                     std::make_unique<DDRotationMatrix>(x.x(), y.x(), z.x(), x.y(), y.y(), z.y(), x.z(), y.z(), z.z()));
0094 }
0095
0096 DDRotation DDrotReflect(const DDName& ddname, std::unique_ptr<DDRotationMatrix> rot) {
0097   return DDRotation(ddname, std::move(rot));
0098 }
0099
0100 // makes sure that the DDRotationMatrix built is LEFT-handed coordinate system (i.e. reflected)
0101 DDRotation DDrotReflect(
0102     const DDName& ddname, double thetaX, double phiX, double thetaY, double phiY, double thetaZ, double phiZ) {
0103   // define 3 unit std::vectors forming the new left-handed axes
0104   DD3Vector x(cos(phiX) * sin(thetaX), sin(phiX) * sin(thetaX), cos(thetaX));
0105   DD3Vector y(cos(phiY) * sin(thetaY), sin(phiY) * sin(thetaY), cos(thetaY));
0106   DD3Vector z(cos(phiZ) * sin(thetaZ), sin(phiZ) * sin(thetaZ), cos(thetaZ));
0107
0108   double tol = 1.0e-3;                 // Geant4 compatible
0109   double check = (x.Cross(y)).Dot(z);  // in case of a LEFT-handed orthogonal system this must be -1
0110   if (fabs(1. + check) > tol) {
0111     edm::LogError("DDRotation") << ddname << " is not a LEFT-handed orthonormal matrix!" << std::endl;
0112     throw cms::Exception("DDException") << ddname.name() << " is not LEFT-handed!";
0113   }
0114
0115   return DDRotation(ddname,
0116                     std::make_unique<DDRotationMatrix>(x.x(), y.x(), z.x(), x.y(), y.y(), z.y(), x.z(), y.z(), z.z()));
0117 }
0118
0119 // does NOT check LEFT or Right handed coordinate system takes either.
0120 std::unique_ptr<DDRotationMatrix> DDcreateRotationMatrix(
0121     double thetaX, double phiX, double thetaY, double phiY, double thetaZ, double phiZ) {
0122   // define 3 unit std::vectors forming the new left-handed axes
0123   DD3Vector x(cos(phiX) * sin(thetaX), sin(phiX) * sin(thetaX), cos(thetaX));
0124   DD3Vector y(cos(phiY) * sin(thetaY), sin(phiY) * sin(thetaY), cos(thetaY));
0125   DD3Vector z(cos(phiZ) * sin(thetaZ), sin(phiZ) * sin(thetaZ), cos(thetaZ));
0126
0127   double tol = 1.0e-3;                 // Geant4 compatible
0128   double check = (x.Cross(y)).Dot(z);  // in case of a LEFT-handed orthogonal system this must be -1, RIGHT-handed: +1
0129   if ((1. - fabs(check)) > tol) {
0130     std::ostringstream o;
0131     o << "matrix is not an (left or right handed) orthonormal matrix! (in deg)" << std::endl
0132       << " thetaX=" << convertRadToDeg(thetaX) << " phiX=" << convertRadToDeg(phiX) << std::endl
0133       << " thetaY=" << convertRadToDeg(thetaY) << " phiY=" << convertRadToDeg(phiY) << std::endl
0134       << " thetaZ=" << convertRadToDeg(thetaZ) << " phiZ=" << convertRadToDeg(phiZ) << std::endl;
0135     edm::LogError("DDRotation") << o.str() << std::endl;
0136
0137     throw cms::Exception("DDException") << o.str();
0138   }
0139
0140   return std::make_unique<DDRotationMatrix>(x.x(), y.x(), z.x(), x.y(), y.y(), z.y(), x.z(), y.z(), z.z());
0141 }
0142
0143 DDRotation DDanonymousRot(std::unique_ptr<DDRotationMatrix> rot) { return DDRotation(std::move(rot)); }