DDLRotationAndReflection.cc

Go to the documentation of this file.

#include "DetectorDescription/Parser/src/DDLRotationAndReflection.h"
#include "DetectorDescription/Core/interface/DDRotationMatrix.h"
#include "DetectorDescription/Core/interface/DDName.h"
#include "DetectorDescription/Core/interface/DDTransform.h"
#include "DetectorDescription/Core/interface/ClhepEvaluator.h"
#include "DataFormats/Math/interface/GeantUnits.h"
#include "DetectorDescription/Parser/interface/DDLElementRegistry.h"
#include "DetectorDescription/Parser/src/DDXMLElement.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "Math/GenVector/Cartesian3D.h"
#include "Math/GenVector/DisplacementVector3D.h"

#include <cmath>
#include <iostream>
#include <map>
#include <utility>

class DDCompactView;

using namespace geant_units::operators;

DDLRotationAndReflection::DDLRotationAndReflection(DDLElementRegistry* myreg) : DDXMLElement(myreg) {}

void DDLRotationAndReflection::processElement(const std::string& name, const std::string& nmspace, DDCompactView& cpv) {
  DD3Vector x = makeX(nmspace);
  DD3Vector y = makeY(nmspace);
  DD3Vector z = makeZ(nmspace);

  DDXMLAttribute atts = getAttributeSet();

  if ((name == "Rotation") && isLeftHanded(x, y, z, nmspace) == 0) {
    DDRotation ddrot = DDrot(getDDName(nmspace), std::make_unique<DDRotationMatrix>(x, y, z));
  } else if ((name == "Rotation") && isLeftHanded(x, y, z, nmspace) == 1) {
    std::string msg("\nDDLRotationAndReflection attempted to make a");
    msg += " left-handed rotation with a Rotation element. If";
    msg += " you meant to make a reflection, use ReflectionRotation";
    msg += " elements, otherwise, please check your matrix.  Other";
    msg += " errors may follow.  Rotation  matrix not created.";
    edm::LogError("DetectorDescription_Parser_Rotation_and_Reflection")
        << msg << std::endl;  // this could become a throwWarning or something.
  } else if (name == "ReflectionRotation" && isLeftHanded(x, y, z, nmspace) == 1) {
    ClhepEvaluator& ev = myRegistry_->evaluator();
    DDRotation ddrot = DDrotReflect(getDDName(nmspace),
                                    ev.eval(nmspace, atts.find("thetaX")->second),
                                    ev.eval(nmspace, atts.find("phiX")->second),
                                    ev.eval(nmspace, atts.find("thetaY")->second),
                                    ev.eval(nmspace, atts.find("phiY")->second),
                                    ev.eval(nmspace, atts.find("thetaZ")->second),
                                    ev.eval(nmspace, atts.find("phiZ")->second));
  } else if (name == "ReflectionRotation" && isLeftHanded(x, y, z, nmspace) == 0) {
    std::string msg("WARNING:  Attempted to make a right-handed");
    msg += " rotation using a ReflectionRotation element. ";
    msg += " If you meant to make a Rotation, use Rotation";
    msg += " elements, otherwise, please check your matrix.";
    msg += "  Other errors may follow.  ReflectionRotation";
    msg += " matrix not created.";
    edm::LogError("DetectorDescription_Parser_Rotation_and_Reflection")
        << msg << std::endl;  // this could be a throwWarning or something.
  } else {
    std::string msg = "\nDDLRotationAndReflection::processElement tried to process wrong element.";
    throwError(msg);
  }
  // after a rotation or reflection rotation has been processed, clear it
  clear();
}

// returns 1 if it is a left-handed CLHEP rotation matrix, 0 if not, but is okay, -1 if
// it is not an orthonormal matrix.
//
// Upon encountering the end tag of a Rotation element, we've got to feed
// the appropriate rotation in to the DDCore.  This is an attempt to do so.
//
// Basically, I cannibalized code from g3tog4 (see http link below) and then
// provided the information from our DDL to the same calls.  Tim Cox showed me
// how to build the rotation matrix (mathematically) and the g3tog4 code basically
// did the rest.
//

int DDLRotationAndReflection::isLeftHanded(const DD3Vector& x,
                                           const DD3Vector& y,
                                           const DD3Vector& z,
                                           const std::string& nmspace) {
  int ret = 0;

  // check for orthonormality and left-handedness

  double check = (x.Cross(y)).Dot(z);
  double tol = 1.0e-3;
  ClhepEvaluator& ev = myRegistry_->evaluator();
  DDXMLAttribute atts = getAttributeSet();

  if (1.0 - std::abs(check) > tol) {
    std::cout << "DDLRotationAndReflection Coordinate axes forming rotation matrix " << getDDName(nmspace)
              << " are not orthonormal.(tolerance=" << tol << " check=" << std::abs(check) << ")" << std::endl
              << " thetaX=" << (atts.find("thetaX")->second) << ' '
              << convertRadToDeg(ev.eval(nmspace, atts.find("thetaX")->second)) << std::endl
              << " phiX=" << (atts.find("phiX")->second) << ' '
              << convertRadToDeg(ev.eval(nmspace, atts.find("phiX")->second)) << std::endl
              << " thetaY=" << (atts.find("thetaY")->second) << ' '
              << convertRadToDeg(ev.eval(nmspace, atts.find("thetaY")->second)) << std::endl
              << " phiY=" << (atts.find("phiY")->second) << ' '
              << convertRadToDeg(ev.eval(nmspace, atts.find("phiY")->second)) << std::endl
              << " thetaZ=" << (atts.find("thetaZ")->second) << ' '
              << convertRadToDeg(ev.eval(nmspace, atts.find("thetaZ")->second)) << std::endl
              << " phiZ=" << (atts.find("phiZ")->second) << ' '
              << convertRadToDeg(ev.eval(nmspace, atts.find("phiZ")->second)) << std::endl
              << "  WAS NOT CREATED!" << std::endl;
    ret = -1;
  } else if (1.0 + check <= tol) {
    ret = 1;
  }
  return ret;
}

DD3Vector DDLRotationAndReflection::makeX(const std::string& nmspace) {
  DD3Vector x;
  DDXMLAttribute atts = getAttributeSet();
  if (atts.find("thetaX") != atts.end()) {
    ClhepEvaluator& ev = myRegistry_->evaluator();
    double thetaX = ev.eval(nmspace, atts.find("thetaX")->second);
    double phiX = ev.eval(nmspace, atts.find("phiX")->second);
    // colx
    x.SetX(sin(thetaX) * cos(phiX));
    x.SetY(sin(thetaX) * sin(phiX));
    x.SetZ(cos(thetaX));
  }
  return x;
}

DD3Vector DDLRotationAndReflection::makeY(const std::string& nmspace) {
  DD3Vector y;
  DDXMLAttribute atts = getAttributeSet();
  if (atts.find("thetaY") != atts.end()) {
    ClhepEvaluator& ev = myRegistry_->evaluator();
    double thetaY = ev.eval(nmspace, atts.find("thetaY")->second);
    double phiY = ev.eval(nmspace, atts.find("phiY")->second);

    // coly
    y.SetX(sin(thetaY) * cos(phiY));
    y.SetY(sin(thetaY) * sin(phiY));
    y.SetZ(cos(thetaY));
  }
  return y;
}

DD3Vector DDLRotationAndReflection::makeZ(const std::string& nmspace) {
  DD3Vector z;
  DDXMLAttribute atts = getAttributeSet();
  if (atts.find("thetaZ") != atts.end()) {
    ClhepEvaluator& ev = myRegistry_->evaluator();
    double thetaZ = ev.eval(nmspace, atts.find("thetaZ")->second);
    double phiZ = ev.eval(nmspace, atts.find("phiZ")->second);

    // colz
    z.SetX(sin(thetaZ) * cos(phiZ));
    z.SetY(sin(thetaZ) * sin(phiZ));
    z.SetZ(cos(thetaZ));
  }
  return z;
}