Parser/src/DDLRotationAndReflection.cc

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;  // this could become a throwWarning or something.
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;  // this could be a throwWarning or something.
0059   } else {
0060     std::string msg = "\nDDLRotationAndReflection::processElement tried to process wrong element.";
0061     throwError(msg);
0062   }
0063   // after a rotation or reflection rotation has been processed, clear it
0064   clear();
0065 }
0066
0067 // returns 1 if it is a left-handed CLHEP rotation matrix, 0 if not, but is okay, -1 if
0068 // it is not an orthonormal matrix.
0069 //
0070 // Upon encountering the end tag of a Rotation element, we've got to feed
0071 // the appropriate rotation in to the DDCore.  This is an attempt to do so.
0072 //
0073 // Basically, I cannibalized code from g3tog4 (see http link below) and then
0074 // provided the information from our DDL to the same calls.  Tim Cox showed me
0075 // how to build the rotation matrix (mathematically) and the g3tog4 code basically
0076 // did the rest.
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   // check for orthonormality and left-handedness
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     // colx
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     // coly
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     // colz
0155     z.SetX(sin(thetaZ) * cos(phiZ));
0156     z.SetY(sin(thetaZ) * sin(phiZ));
0157     z.SetZ(cos(thetaZ));
0158   }
0159   return z;
0160 }