Skip to content

Commit

Permalink
Merge pull request #35299 from cvuosalo/rotReorder
Browse files Browse the repository at this point in the history 
[DD4hep] Increase precision of rotation matching in making DD4hep big XML file for DB
  • Loading branch information
@cmsbuild
cmsbuild committed Sep 17, 2021
2 parents 7cee8f4 + 2fd71aa commit e662c63
Show file tree
Hide file tree
Showing 2 changed files with 61 additions and 38 deletions.
19 changes: 15 additions & 4 deletions DetectorDescription/DDCMS/src/DDNamespace.cc
View file
Expand Up @@ -6,6 +6,8 @@
#include "XML/XML.h"

#include <TClass.h>
#include <iomanip>
#include <sstream>
#include <unordered_map>
#include <vector>

Expand All @@ -14,7 +16,7 @@ using namespace cms;

double cms::rotation_utils::roundBinary(double value) {
value = cms_rounding::roundIfNear0(value);
static constexpr double roundingVal = 1 << 14;
static constexpr double roundingVal = 1 << 24;
value = (round(value * roundingVal) / roundingVal);
// Set -0 to 0
return (cms_rounding::roundIfNear0(value));
Expand All @@ -24,7 +26,10 @@ std::string cms::rotation_utils::rotHash(const Double_t* rot) {
std::string hashVal;
for (int row = 0; row <= 2; ++row) {
for (int col = 0; col <= 2; ++col) {
hashVal += std::to_string(roundBinary(rot[(3 * row) + col]));
std::ostringstream numStream;
numStream << std::fixed << std::setprecision(7);
numStream << roundBinary(rot[(3 * row) + col]);
hashVal += numStream.str();
}
}
return (hashVal);
Expand All @@ -36,7 +41,10 @@ std::string cms::rotation_utils::rotHash(const dd4hep::Rotation3D& rot) {
matrix.assign(9, 0.);
rot.GetComponents(matrix.begin());
for (double val : matrix) {
hashVal += std::to_string(roundBinary(val));
std::ostringstream numStream;
numStream << std::fixed << std::setprecision(7);
numStream << roundBinary(val);
hashVal += numStream.str();
}
return (hashVal);
}
Expand Down Expand Up @@ -152,7 +160,10 @@ void DDNamespace::addRotation(const string& name, const dd4hep::Rotation3D& rot)
m_context->rotations[n] = rot;
if (m_context->makePayload) {
string hashVal = cms::rotation_utils::rotHash(rot);
m_context->rotRevMap[hashVal] = n;
if (m_context->rotRevMap.find(hashVal) == m_context->rotRevMap.end()) {
// Only set a rotation that is not already in the map
m_context->rotRevMap[hashVal] = n;
}
}
}

Expand Down
80 changes: 46 additions & 34 deletions DetectorDescription/OfflineDBLoader/src/DDCoreToDDXMLOutput.cc
View file
Expand Up @@ -39,6 +39,9 @@ static inline constexpr NumType convertGPerCcToMgPerCc(NumType gPerCc) // g/cm^
return (gPerCc * 1000.);
}

static constexpr double tol0 = 1.e-11; // Tiny values to be considered equal to 0
static constexpr double reflectTol = 1.0e-3; // Tolerance for recognizing reflections; Geant4-compatible

namespace cms::rotation_utils {
/* For debugging
static double determinant(const dd4hep::Rotation3D &rot) {
Expand All @@ -51,7 +54,8 @@ namespace cms::rotation_utils {
}
*/

static const std::string identityHash("1.0000000.0000000.0000000.0000001.0000000.0000000.0000000.0000001.000000");
static const std::string identityHash(
"1.00000000.00000000.00000000.00000001.00000000.00000000.00000000.00000001.0000000");

static void addRotWithNewName(cms::DDNamespace& ns, std::string& name, const dd4hep::Rotation3D& rot) {
const dd4hep::Rotation3D& rot2 = rot;
Expand All @@ -61,15 +65,15 @@ namespace cms::rotation_utils {

static void addRotWithNewName(cms::DDNamespace& ns, std::string& name, const Double_t* rot) {
using namespace cms_rounding;
dd4hep::Rotation3D rot2(roundIfNear0(rot[0]),
roundIfNear0(rot[1]),
roundIfNear0(rot[2]),
roundIfNear0(rot[3]),
roundIfNear0(rot[4]),
roundIfNear0(rot[5]),
roundIfNear0(rot[6]),
roundIfNear0(rot[7]),
roundIfNear0(rot[8]));
dd4hep::Rotation3D rot2(roundIfNear0(rot[0], tol0),
roundIfNear0(rot[1], tol0),
roundIfNear0(rot[2], tol0),
roundIfNear0(rot[3], tol0),
roundIfNear0(rot[4], tol0),
roundIfNear0(rot[5], tol0),
roundIfNear0(rot[6], tol0),
roundIfNear0(rot[7], tol0),
roundIfNear0(rot[8], tol0));
addRotWithNewName(ns, name, rot2);
}

Expand Down Expand Up @@ -115,8 +119,12 @@ void DDCoreToDDXMLOutput::solid(const dd4hep::Solid& solid, const cms::DDParsing
xos << " y=\"" << trans[1] << "*mm\"";
xos << " z=\"" << trans[2] << "*mm\"";
xos << "/>" << std::endl;
std::string rotNameStr = cms::rotation_utils::rotName(rs.rightMatrix()->GetRotationMatrix(), context);
xos << "<rRotation name=\"" << rotNameStr << "\"/>" << std::endl;
auto rot = rs.rightMatrix()->GetRotationMatrix();
// The identity rotation can be omitted.
if (cms::rotation_utils::rotHash(rot) != cms::rotation_utils::identityHash) {
std::string rotNameStr = cms::rotation_utils::rotName(rot, context);
xos << "<rRotation name=\"" << rotNameStr << "\"/>" << std::endl;
}
if (shape == cms::DDSolidShape::ddunion) {
xos << "</UnionSolid>" << std::endl;
} else if (shape == cms::DDSolidShape::ddsubtraction) {
Expand Down Expand Up @@ -654,19 +662,21 @@ void DDCoreToDDXMLOutput::element(const TGeoMaterial* material, std::ostream& xo
}

void DDCoreToDDXMLOutput::rotation(const DDRotation& rotation, std::ostream& xos, const std::string& rotn) {
double tol = 1.0e-3; // Geant4 compatible
DD3Vector x, y, z;
rotation.rotation().GetComponents(x, y, z);
double a, b, c;
x.GetCoordinates(a, b, c);
x.SetCoordinates(cms_rounding::roundIfNear0(a), cms_rounding::roundIfNear0(b), cms_rounding::roundIfNear0(c));
x.SetCoordinates(
cms_rounding::roundIfNear0(a, tol0), cms_rounding::roundIfNear0(b, tol0), cms_rounding::roundIfNear0(c, tol0));
y.GetCoordinates(a, b, c);
y.SetCoordinates(cms_rounding::roundIfNear0(a), cms_rounding::roundIfNear0(b), cms_rounding::roundIfNear0(c));
y.SetCoordinates(
cms_rounding::roundIfNear0(a, tol0), cms_rounding::roundIfNear0(b, tol0), cms_rounding::roundIfNear0(c, tol0));
z.GetCoordinates(a, b, c);
z.SetCoordinates(cms_rounding::roundIfNear0(a), cms_rounding::roundIfNear0(b), cms_rounding::roundIfNear0(c));
z.SetCoordinates(
cms_rounding::roundIfNear0(a, tol0), cms_rounding::roundIfNear0(b, tol0), cms_rounding::roundIfNear0(c, tol0));
double check = (x.Cross(y)).Dot(z); // in case of a LEFT-handed orthogonal system
// this must be -1
bool reflection((1. - check) > tol);
bool reflection((1. - check) > reflectTol);
std::string rotName = rotation.toString();
if (rotName == ":") {
if (!rotn.empty()) {
Expand All @@ -685,44 +695,46 @@ void DDCoreToDDXMLOutput::rotation(const DDRotation& rotation, std::ostream& xos
}
using namespace cms_rounding;
xos << "name=\"" << rotName << "\""
<< " phiX=\"" << roundIfNear0(convertRadToDeg(x.phi()), 4.e-4) << "*deg\""
<< " thetaX=\"" << roundIfNear0(convertRadToDeg(x.theta()), 4.e-4) << "*deg\""
<< " phiY=\"" << roundIfNear0(convertRadToDeg(y.phi()), 4.e-4) << "*deg\""
<< " thetaY=\"" << roundIfNear0(convertRadToDeg(y.theta()), 4.e-4) << "*deg\""
<< " phiZ=\"" << roundIfNear0(convertRadToDeg(z.phi()), 4.e-4) << "*deg\""
<< " thetaZ=\"" << roundIfNear0(convertRadToDeg(z.theta()), 4.e-4) << "*deg\"/>" << std::endl;
<< " phiX=\"" << roundIfNear0(convertRadToDeg(x.phi()), tol0) << "*deg\""
<< " thetaX=\"" << roundIfNear0(convertRadToDeg(x.theta()), tol0) << "*deg\""
<< " phiY=\"" << roundIfNear0(convertRadToDeg(y.phi()), tol0) << "*deg\""
<< " thetaY=\"" << roundIfNear0(convertRadToDeg(y.theta()), tol0) << "*deg\""
<< " phiZ=\"" << roundIfNear0(convertRadToDeg(z.phi()), tol0) << "*deg\""
<< " thetaZ=\"" << roundIfNear0(convertRadToDeg(z.theta()), tol0) << "*deg\"/>" << std::endl;
}

void DDCoreToDDXMLOutput::rotation(const dd4hep::Rotation3D& rotation,
std::ostream& xos,
const cms::DDParsingContext& context,
const std::string& rotn) {
double tol = 1.0e-3; // Geant4 compatible
ROOT::Math::XYZVector x, y, z;
rotation.GetComponents(x, y, z);
double a, b, c;
x.GetCoordinates(a, b, c);
x.SetCoordinates(cms_rounding::roundIfNear0(a), cms_rounding::roundIfNear0(b), cms_rounding::roundIfNear0(c));
x.SetCoordinates(
cms_rounding::roundIfNear0(a, tol0), cms_rounding::roundIfNear0(b, tol0), cms_rounding::roundIfNear0(c, tol0));
y.GetCoordinates(a, b, c);
y.SetCoordinates(cms_rounding::roundIfNear0(a), cms_rounding::roundIfNear0(b), cms_rounding::roundIfNear0(c));
y.SetCoordinates(
cms_rounding::roundIfNear0(a, tol0), cms_rounding::roundIfNear0(b, tol0), cms_rounding::roundIfNear0(c, tol0));
z.GetCoordinates(a, b, c);
z.SetCoordinates(cms_rounding::roundIfNear0(a), cms_rounding::roundIfNear0(b), cms_rounding::roundIfNear0(c));
z.SetCoordinates(
cms_rounding::roundIfNear0(a, tol0), cms_rounding::roundIfNear0(b, tol0), cms_rounding::roundIfNear0(c, tol0));
double check = (x.Cross(y)).Dot(z); // in case of a LEFT-handed orthogonal system
// this must be -1
bool reflection((1. - check) > tol);
bool reflection((1. - check) > reflectTol);
if (!reflection) {
xos << "<Rotation ";
} else {
xos << "<ReflectionRotation ";
}
using namespace cms_rounding;
xos << "name=\"" << rotn << "\""
<< " phiX=\"" << roundIfNear0(convertRadToDeg(x.phi()), 4.e-4) << "*deg\""
<< " thetaX=\"" << roundIfNear0(convertRadToDeg(x.theta()), 4.e-4) << "*deg\""
<< " phiY=\"" << roundIfNear0(convertRadToDeg(y.phi()), 4.e-4) << "*deg\""
<< " thetaY=\"" << roundIfNear0(convertRadToDeg(y.theta()), 4.e-4) << "*deg\""
<< " phiZ=\"" << roundIfNear0(convertRadToDeg(z.phi()), 4.e-4) << "*deg\""
<< " thetaZ=\"" << roundIfNear0(convertRadToDeg(z.theta()), 4.e-4) << "*deg\"/>" << std::endl;
<< " phiX=\"" << roundIfNear0(convertRadToDeg(x.phi()), tol0) << "*deg\""
<< " thetaX=\"" << roundIfNear0(convertRadToDeg(x.theta()), tol0) << "*deg\""
<< " phiY=\"" << roundIfNear0(convertRadToDeg(y.phi()), tol0) << "*deg\""
<< " thetaY=\"" << roundIfNear0(convertRadToDeg(y.theta()), tol0) << "*deg\""
<< " phiZ=\"" << roundIfNear0(convertRadToDeg(z.phi()), tol0) << "*deg\""
<< " thetaZ=\"" << roundIfNear0(convertRadToDeg(z.theta()), tol0) << "*deg\"/>" << std::endl;
}

void DDCoreToDDXMLOutput::logicalPart(const DDLogicalPart& lp, std::ostream& xos) {
Expand Down

0 comments on commit e662c63

Please sign in to comment.