Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

[DD4hep] Increase precision of rotation matching in making DD4hep big XML file for DB #35299

Merged
merged 1 commit into from Sep 17, 2021
Merged

[DD4hep] Increase precision of rotation matching in making DD4hep big XML file for DB #35299

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
19 changes: 15 additions & 4 deletions DetectorDescription/DDCMS/src/DDNamespace.cc
View file
Open in desktop
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
Open in desktop
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