File indexing completed on 2024-04-06 12:05:33
0001 #include "DetectorDescription/Parser/src/DDLRotationAndReflection.h"
0002 #include "DetectorDescription/Core/interface/DDRotationMatrix.h"
0003 #include "DetectorDescription/Core/interface/DDName.h"
0004 #include "DetectorDescription/Core/interface/DDTransform.h"
0005 #include "DetectorDescription/Core/interface/ClhepEvaluator.h"
0006 #include "DataFormats/Math/interface/GeantUnits.h"
0007 #include "DetectorDescription/Parser/interface/DDLElementRegistry.h"
0008 #include "DetectorDescription/Parser/src/DDXMLElement.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 #include "Math/GenVector/Cartesian3D.h"
0011 #include "Math/GenVector/DisplacementVector3D.h"
0012
0013 #include <cmath>
0014 #include <iostream>
0015 #include <map>
0016 #include <utility>
0017
0018 class DDCompactView;
0019
0020 using namespace geant_units::operators;
0021
0022 DDLRotationAndReflection::DDLRotationAndReflection(DDLElementRegistry* myreg) : DDXMLElement(myreg) {}
0023
0024 void DDLRotationAndReflection::processElement(const std::string& name, const std::string& nmspace, DDCompactView& cpv) {
0025 DD3Vector x = makeX(nmspace);
0026 DD3Vector y = makeY(nmspace);
0027 DD3Vector z = makeZ(nmspace);
0028
0029 DDXMLAttribute atts = getAttributeSet();
0030
0031 if ((name == "Rotation") && isLeftHanded(x, y, z, nmspace) == 0) {
0032 DDRotation ddrot = DDrot(getDDName(nmspace), std::make_unique<DDRotationMatrix>(x, y, z));
0033 } else if ((name == "Rotation") && isLeftHanded(x, y, z, nmspace) == 1) {
0034 std::string msg("\nDDLRotationAndReflection attempted to make a");
0035 msg += " left-handed rotation with a Rotation element. If";
0036 msg += " you meant to make a reflection, use ReflectionRotation";
0037 msg += " elements, otherwise, please check your matrix. Other";
0038 msg += " errors may follow. Rotation matrix not created.";
0039 edm::LogError("DetectorDescription_Parser_Rotation_and_Reflection")
0040 << msg << std::endl;
0041 } else if (name == "ReflectionRotation" && isLeftHanded(x, y, z, nmspace) == 1) {
0042 ClhepEvaluator& ev = myRegistry_->evaluator();
0043 DDRotation ddrot = DDrotReflect(getDDName(nmspace),
0044 ev.eval(nmspace, atts.find("thetaX")->second),
0045 ev.eval(nmspace, atts.find("phiX")->second),
0046 ev.eval(nmspace, atts.find("thetaY")->second),
0047 ev.eval(nmspace, atts.find("phiY")->second),
0048 ev.eval(nmspace, atts.find("thetaZ")->second),
0049 ev.eval(nmspace, atts.find("phiZ")->second));
0050 } else if (name == "ReflectionRotation" && isLeftHanded(x, y, z, nmspace) == 0) {
0051 std::string msg("WARNING: Attempted to make a right-handed");
0052 msg += " rotation using a ReflectionRotation element. ";
0053 msg += " If you meant to make a Rotation, use Rotation";
0054 msg += " elements, otherwise, please check your matrix.";
0055 msg += " Other errors may follow. ReflectionRotation";
0056 msg += " matrix not created.";
0057 edm::LogError("DetectorDescription_Parser_Rotation_and_Reflection")
0058 << msg << std::endl;
0059 } else {
0060 std::string msg = "\nDDLRotationAndReflection::processElement tried to process wrong element.";
0061 throwError(msg);
0062 }
0063
0064 clear();
0065 }
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077
0078
0079 int DDLRotationAndReflection::isLeftHanded(const DD3Vector& x,
0080 const DD3Vector& y,
0081 const DD3Vector& z,
0082 const std::string& nmspace) {
0083 int ret = 0;
0084
0085
0086
0087 double check = (x.Cross(y)).Dot(z);
0088 double tol = 1.0e-3;
0089 ClhepEvaluator& ev = myRegistry_->evaluator();
0090 DDXMLAttribute atts = getAttributeSet();
0091
0092 if (1.0 - std::abs(check) > tol) {
0093 std::cout << "DDLRotationAndReflection Coordinate axes forming rotation matrix " << getDDName(nmspace)
0094 << " are not orthonormal.(tolerance=" << tol << " check=" << std::abs(check) << ")" << std::endl
0095 << " thetaX=" << (atts.find("thetaX")->second) << ' '
0096 << convertRadToDeg(ev.eval(nmspace, atts.find("thetaX")->second)) << std::endl
0097 << " phiX=" << (atts.find("phiX")->second) << ' '
0098 << convertRadToDeg(ev.eval(nmspace, atts.find("phiX")->second)) << std::endl
0099 << " thetaY=" << (atts.find("thetaY")->second) << ' '
0100 << convertRadToDeg(ev.eval(nmspace, atts.find("thetaY")->second)) << std::endl
0101 << " phiY=" << (atts.find("phiY")->second) << ' '
0102 << convertRadToDeg(ev.eval(nmspace, atts.find("phiY")->second)) << std::endl
0103 << " thetaZ=" << (atts.find("thetaZ")->second) << ' '
0104 << convertRadToDeg(ev.eval(nmspace, atts.find("thetaZ")->second)) << std::endl
0105 << " phiZ=" << (atts.find("phiZ")->second) << ' '
0106 << convertRadToDeg(ev.eval(nmspace, atts.find("phiZ")->second)) << std::endl
0107 << " WAS NOT CREATED!" << std::endl;
0108 ret = -1;
0109 } else if (1.0 + check <= tol) {
0110 ret = 1;
0111 }
0112 return ret;
0113 }
0114
0115 DD3Vector DDLRotationAndReflection::makeX(const std::string& nmspace) {
0116 DD3Vector x;
0117 DDXMLAttribute atts = getAttributeSet();
0118 if (atts.find("thetaX") != atts.end()) {
0119 ClhepEvaluator& ev = myRegistry_->evaluator();
0120 double thetaX = ev.eval(nmspace, atts.find("thetaX")->second);
0121 double phiX = ev.eval(nmspace, atts.find("phiX")->second);
0122
0123 x.SetX(sin(thetaX) * cos(phiX));
0124 x.SetY(sin(thetaX) * sin(phiX));
0125 x.SetZ(cos(thetaX));
0126 }
0127 return x;
0128 }
0129
0130 DD3Vector DDLRotationAndReflection::makeY(const std::string& nmspace) {
0131 DD3Vector y;
0132 DDXMLAttribute atts = getAttributeSet();
0133 if (atts.find("thetaY") != atts.end()) {
0134 ClhepEvaluator& ev = myRegistry_->evaluator();
0135 double thetaY = ev.eval(nmspace, atts.find("thetaY")->second);
0136 double phiY = ev.eval(nmspace, atts.find("phiY")->second);
0137
0138
0139 y.SetX(sin(thetaY) * cos(phiY));
0140 y.SetY(sin(thetaY) * sin(phiY));
0141 y.SetZ(cos(thetaY));
0142 }
0143 return y;
0144 }
0145
0146 DD3Vector DDLRotationAndReflection::makeZ(const std::string& nmspace) {
0147 DD3Vector z;
0148 DDXMLAttribute atts = getAttributeSet();
0149 if (atts.find("thetaZ") != atts.end()) {
0150 ClhepEvaluator& ev = myRegistry_->evaluator();
0151 double thetaZ = ev.eval(nmspace, atts.find("thetaZ")->second);
0152 double phiZ = ev.eval(nmspace, atts.find("phiZ")->second);
0153
0154
0155 z.SetX(sin(thetaZ) * cos(phiZ));
0156 z.SetY(sin(thetaZ) * sin(phiZ));
0157 z.SetZ(cos(thetaZ));
0158 }
0159 return z;
0160 }