Skip to content

Commit

Permalink
Merge branch 'master' into edm4hep-mcp-momentum-double
Browse files Browse the repository at this point in the history
  • Loading branch information
@MarkusFrankATcernch
MarkusFrankATcernch authored Jan 17, 2024
2 parents 5c42c9f + c0120cd commit b8f2934
Show file tree
Hide file tree
Showing 10 changed files with 610 additions and 220 deletions.
4 changes: 2 additions & 2 deletions DDCAD/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,8 +37,8 @@ elseif(DD4HEP_LOAD_ASSIMP)

ExternalProject_Add(
assimp_project
SOURCE_DIR ${CMAKE_CURRENT_BINARY_DIR}/external/assimp
INSTALL_DIR ${CMAKE_INSTALL_PREFIX}
SOURCE_DIR ${CMAKE_CURRENT_BINARY_DIR}/external/assimp
INSTALL_DIR ${CMAKE_INSTALL_PREFIX}
GIT_REPOSITORY https://github.com/assimp/assimp.git
GIT_TAG v5.0.1
CMAKE_ARGS -DCMAKE_INSTALL_PREFIX:PATH=${CMAKE_INSTALL_PREFIX}
Expand Down
74 changes: 45 additions & 29 deletions DDCAD/include/DDCAD/Utilities.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,8 @@
#define DDCAD_UTILITIES_H

#include <vector>
#include <cmath>
#include <limits>

#include <TGeoTessellated.h>
#include <TGeoVector3.h>
Expand All @@ -24,52 +26,66 @@ namespace dd4hep {
/// Namespace for implementation details of the AIDA detector description toolkit
namespace cad {

inline std::stringstream streamFacet(TGeoFacet const& facet,
TGeoTessellated const& shape) {
#if ROOT_VERSION_CODE >= ROOT_VERSION(6,31,1)
inline std::string streamFacet(TGeoFacet const& facet,
TGeoTessellated const& shape) {
using ::operator<<;
std::stringstream str;
#if ROOT_VERSION_CODE >= ROOT_VERSION(6,31,1)
str << "{";
for (int i = 0; i < facet.GetNvert(); ++i) {
str << shape.GetVertex(facet[i]);
if (i != facet.GetNvert() - 1)
str << ", ";
str << shape.GetVertex(facet[i]);
if (i != facet.GetNvert() - 1)
str << ", ";
}
str << "}";
return str.str();
}
#else
inline std::string streamFacet(TGeoFacet const& facet,
TGeoTessellated const& /* shape */) {
using ::operator<<;
std::stringstream str;
str << facet;
#endif
return str;
return str.str();
}

inline std::stringstream streamVertices(ROOT::Geom::Vertex_t const& v1,
#endif

inline std::string streamVertices(ROOT::Geom::Vertex_t const& v1,
ROOT::Geom::Vertex_t const& v2,
ROOT::Geom::Vertex_t const& v3) {
using ::operator<<;
std::stringstream str;
str << "{" << v1 << ", " << v2 << ", " << v3 << "}";
return str;
return str.str();
}

// Determine if the facet is degenerated by calculating its determinant
inline bool facetIsDegenerated(std::vector<ROOT::Geom::Vertex_t> const& vertices){
const ROOT::Geom::Vertex_t& v1 = vertices[0];
const ROOT::Geom::Vertex_t& v2 = vertices[1];
const ROOT::Geom::Vertex_t& v3 = vertices[2];
constexpr double epsilon = 1.e-20;
// v1.x v2.x v3.x v1.x v2.x
//
// v1.y v2.y v3.y v1.y v2.y
//
// v1.z v2.z v3.z v1.z v2.z
double det = 0.0
+ v1.x() * v2.y() * v3.z()
+ v2.x() * v3.y() * v1.z()
+ v3.x() * v1.y() * v2.z()
- v1.z() * v2.y() * v3.x()
- v2.z() * v3.y() * v1.x()
- v3.z() * v1.y() * v2.x();
return det < epsilon;
using ROOT::Geom::Vertex_t;
// Compute normal using non-zero segments
constexpr double kTolerance = 1.e-20;
bool degenerated = true;
Vertex_t normal;
int fNvert = int(vertices.size());
for (int i = 0; i < fNvert - 1; ++i) {
Vertex_t e1 = vertices[i + 1] - vertices[i];
if (e1.Mag2() < kTolerance)
continue;
for (int j = i + 1; j < fNvert; ++j) {
Vertex_t e2 = vertices[(j + 1) % fNvert] - vertices[j];
if (e2.Mag2() < kTolerance)
continue;
normal = Vertex_t::Cross(e1, e2);
// e1 and e2 may be colinear
if (normal.Mag2() < kTolerance)
continue;
normal.Normalize();
degenerated = false;
break;
}
if (!degenerated)
break;
}
return degenerated;
}
}
}
Expand Down
55 changes: 43 additions & 12 deletions DDCAD/src/ASSIMPReader.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,9 +63,9 @@ ASSIMPReader::readShapes(const std::string& source, double unit_length) const
if ( dump_facets ) {
for( size_t i=0, n=shape->GetNfacets(); i < n; ++i ) {
const auto& facet = shape->GetFacet(i);
std::stringstream str = dd4hep::cad::streamFacet(facet, shape);
std::string str = dd4hep::cad::streamFacet(facet, shape);
printout(ALWAYS,"ASSIMPReader","++ Facet %4ld : %s",
i, str.str().c_str());
i, str.c_str());
}
}
shape->SetTitle(TESSELLATEDSOLID_TAG);
Expand All @@ -90,7 +90,6 @@ ASSIMPReader::readVolumes(const std::string& source, double unit_length) const
bool dump_facets = ((flags>>8)&0x1) == 1;
int aiflags = aiProcess_Triangulate|aiProcess_JoinIdenticalVertices|aiProcess_CalcTangentSpace;
auto scene = importer->ReadFile( source.c_str(), aiflags);
char text[1048];

if ( !scene ) {
except("ASSIMPReader","+++ FileNotFound: %s",source.c_str());
Expand All @@ -108,9 +107,45 @@ ASSIMPReader::readVolumes(const std::string& source, double unit_length) const
vertices.emplace_back(Vertex(v[i].x*unit, v[i].y*unit, v[i].z*unit));
}
TessellatedSolid shape(name,vertices);
if ( name.empty() ) {
name = _toString(result.size(), "tessellated_%ld");
}

/// NOTE: IMPORTANT!
/// ALWAYS add facets using the physical vertices!
/// TGeoTessellated takes care that the vertex map is unique and
/// assigns the proper indices to the facet.
for(unsigned int i=0; i < mesh->mNumFaces; i++) {
const unsigned int* idx = mesh->mFaces[i].mIndices;
shape->AddFacet(idx[0], idx[1], idx[2]);
bool degenerated = false;
if ( mesh->mFaces[i].mNumIndices == 3 )
degenerated = dd4hep::cad::facetIsDegenerated({vertices[idx[0]], vertices[idx[1]], vertices[idx[2]]});
else if ( mesh->mFaces[i].mNumIndices == 4 )
degenerated = dd4hep::cad::facetIsDegenerated({vertices[idx[0]], vertices[idx[1]], vertices[idx[2]], vertices[idx[3]]});

if ( degenerated ) {
printout(DEBUG, "ASSIMPReader", "+++ %s: Drop degenerated facet: %d %d %d",
name.c_str(), idx[0], idx[1], idx[2]);
}
#if ROOT_VERSION_CODE >= ROOT_VERSION(6,31,1)
else if ( mesh->mFaces[i].mNumIndices == 3 ) {
shape->AddFacet(vertices[idx[0]], vertices[idx[1]], vertices[idx[2]]);
}
else if ( mesh->mFaces[i].mNumIndices == 4 ) {
shape->AddFacet(vertices[idx[0]], vertices[idx[1]], vertices[idx[2]], vertices[idx[3]]);
}
#else
else if ( mesh->mFaces[i].mNumIndices == 3 ) {
shape->AddFacet(idx[0], idx[1], idx[2]);
}
else if ( mesh->mFaces[i].mNumIndices == 4 ) {
shape->AddFacet(idx[0], idx[1], idx[2], idx[3]);
}
#endif
else {
printout(INFO, "ASSIMPReader", "+++ %s: Fancy facet with %d indices.",
name.c_str(), mesh->mFaces[i].mNumIndices);
}
}
if ( shape->GetNfacets() > 2 ) {
std::string mat_name;
Expand All @@ -124,14 +159,9 @@ ASSIMPReader::readVolumes(const std::string& source, double unit_length) const
if ( !mat.isValid() ) {
printout(ERROR, "ASSIMPReader",
"+++ %s: No material named '%s' FOUND. Will use Air. [Missing material]",
text, mat_name.c_str());
name.c_str(), mat_name.c_str());
mat = detector.air();
}
if ( name.empty() ) {
::snprintf(text,sizeof(text),"tessellated_%ld", result.size());
text[sizeof(text)-1] = 0;
name = text;
}
Volume vol(name, Solid(shape.ptr()), mat);
if ( mesh->HasVertexColors(0) ) {
const aiColor4D* col = mesh->mColors[0];
Expand All @@ -146,6 +176,7 @@ ASSIMPReader::readVolumes(const std::string& source, double unit_length) const
}
}
if ( !vis.isValid() ) {
char text[1024];
::snprintf(text,sizeof(text),"vis_%s_%p", name.c_str(), (void*)vol.ptr());
text[sizeof(text)-1] = 0;
vis = VisAttr(text);
Expand All @@ -162,9 +193,9 @@ ASSIMPReader::readVolumes(const std::string& source, double unit_length) const
if ( dump_facets ) {
for( size_t i=0, n=shape->GetNfacets(); i < n; ++i ) {
const auto& facet = shape->GetFacet(i);
std::stringstream str = dd4hep::cad::streamFacet(facet, shape);
std::string str = dd4hep::cad::streamFacet(facet, shape);
printout(ALWAYS,"ASSIMPReader","++ Facet %4ld : %s",
i, str.str().c_str());
i, str.c_str());
}
}
result.emplace_back(std::unique_ptr<TGeoVolume>(vol.ptr()));
Expand Down
6 changes: 3 additions & 3 deletions DDCAD/src/ASSIMPWriter.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -159,7 +159,7 @@ namespace {
continue;
}
#if ROOT_VERSION_CODE >= ROOT_VERSION(6,31,1)
bool degenerated = dd4hep::cad::facetIsDegenerated(vertices);
bool degenerated = dd4hep::cad::facetIsDegenerated({vertices[vv0],vertices[vv1],vertices[vv2]});
#else
bool degenerated = true;
TGeoFacet f(&vertices, 3, vv0, vv1, vv2);
Expand Down Expand Up @@ -426,8 +426,8 @@ int ASSIMPWriter::write(const std::string& file_name,
ROOT::Geom::Vertex_t v1(vv[id[0]].x, vv[id[0]].y, vv[id[0]].z);
ROOT::Geom::Vertex_t v2(vv[id[1]].x, vv[id[1]].y, vv[id[1]].z);
ROOT::Geom::Vertex_t v3(vv[id[2]].x, vv[id[2]].y, vv[id[2]].z);
std::stringstream str = dd4hep::cad::streamVertices(v1, v2, v3);
printout(ALWAYS,"ASSIMPWriter","++ Facet %4ld : %s", j, str.str().c_str());
std::string str = dd4hep::cad::streamVertices(v1, v2, v3);
printout(ALWAYS,"ASSIMPWriter","++ Facet %4ld : %s", j, str.c_str());
}
}
else {
Expand Down
4 changes: 3 additions & 1 deletion DDCore/include/DD4hep/Handle.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -239,10 +239,12 @@ namespace dd4hep {

/// String conversions: boolean value to string \ingroup DD4HEP_CORE
std::string _toString(bool value);
/// String conversions: integer value to string \ingroup DD4HEP_CORE
/// String conversions: short integer value to string \ingroup DD4HEP_CORE
std::string _toString(short value, const char* fmt = "%d");
/// String conversions: integer value to string \ingroup DD4HEP_CORE
std::string _toString(int value, const char* fmt = "%d");
/// String conversions: unsigned long integer value to string \ingroup DD4HEP_CORE
std::string _toString(unsigned long value, const char* fmt = "%ld");
/// String conversions: float value to string \ingroup DD4HEP_CORE
std::string _toString(float value, const char* fmt = "%.17e");
/// String conversions: double value to string \ingroup DD4HEP_CORE
Expand Down
4 changes: 4 additions & 0 deletions DDCore/src/Handle.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -344,6 +344,10 @@ namespace dd4hep {
return __to_string(value, fmt);
}

string _toString(unsigned long value, const char* fmt) {
return __to_string(value, fmt);
}

string _toString(float value, const char* fmt) {
return __to_string(value, fmt);
}
Expand Down
73 changes: 41 additions & 32 deletions DDCore/src/plugins/DetectorCheck.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -97,6 +97,7 @@ namespace {
bool check_placements { false };
bool check_volmgr { false };
bool check_sensitive { false };
bool ignore_detector { false };

SensitiveDetector get_current_sensitive_detector();

Expand Down Expand Up @@ -449,7 +450,7 @@ void DetectorCheck::checkManagerSingleVolume(DetElement detector, PlacedVolume p
if ( pv.volume().isSensitive() ) {
PlacedVolume det_place = m_volMgr.lookupDetElementPlacement(vid);
++m_sens_counters.elements;
if ( pv.ptr() != det_place.ptr() ) {
if ( !ignore_detector && pv.ptr() != det_place.ptr() ) {
err << "VolumeMgrTest: Wrong placement "
<< " got " << det_place.name() << " (" << (void*)det_place.ptr() << ")"
<< " instead of " << pv.name() << " (" << (void*)pv.ptr() << ") "
Expand Down Expand Up @@ -551,33 +552,35 @@ void DetectorCheck::checkManagerSingleVolume(DetElement detector, PlacedVolume p
printout(ERROR, m_det.name(), "DETELEMENT_PERSISTENCY FAILED: World transformation have DIFFERET pointer!");
++m_place_counters.errors;
}

if ( pv.ptr() == det_elem.placement().ptr() ) {
// The computed transformation 'trafo' MUST be equal to:
// m_volMgr.worldTransformation(vid) AND det_elem.nominal().worldTransformation()
int res1 = detail::matrix::_matrixEqual(trafo, det_elem.nominal().worldTransformation());
int res2 = detail::matrix::_matrixEqual(trafo, m_volMgr.worldTransformation(m_mapping,vid));
if ( res1 != detail::matrix::MATRICES_EQUAL || res2 != detail::matrix::MATRICES_EQUAL ) {
printout(ERROR, m_det.name(), "DETELEMENT_PLACEMENT FAILED: World transformation DIFFER.");
++m_place_counters.errors;
}
else {
printout(INFO, m_det.name(), "DETELEMENT_PLACEMENT: PASSED. All matrices equal: %s",
volumeID(vid).c_str());
}
}
else {
// The computed transformation 'trafo' MUST be equal to:
// m_volMgr.worldTransformation(vid)
// The det_elem.nominal().worldTransformation() however is DIFFERENT!
int res2 = detail::matrix::_matrixEqual(trafo, m_volMgr.worldTransformation(m_mapping,vid));
if ( res2 != detail::matrix::MATRICES_EQUAL ) {
printout(ERROR, m_det.name(), "VOLUME_PLACEMENT FAILED: World transformation DIFFER.");
++m_place_counters.errors;

if ( !ignore_detector ) {
if ( pv.ptr() == det_elem.placement().ptr() ) {
// The computed transformation 'trafo' MUST be equal to:
// m_volMgr.worldTransformation(vid) AND det_elem.nominal().worldTransformation()
int res1 = detail::matrix::_matrixEqual(trafo, det_elem.nominal().worldTransformation());
int res2 = detail::matrix::_matrixEqual(trafo, m_volMgr.worldTransformation(m_mapping,vid));
if ( res1 != detail::matrix::MATRICES_EQUAL || res2 != detail::matrix::MATRICES_EQUAL ) {
printout(ERROR, m_det.name(), "DETELEMENT_PLACEMENT FAILED: World transformation DIFFER.");
++m_place_counters.errors;
}
else {
printout(INFO, m_det.name(), "DETELEMENT_PLACEMENT: PASSED. All matrices equal: %s",
volumeID(vid).c_str());
}
}
else {
printout(INFO, m_det.name(), "VOLUME_PLACEMENT: PASSED. All matrices equal: %s",
volumeID(vid).c_str());
// The computed transformation 'trafo' MUST be equal to:
// m_volMgr.worldTransformation(vid)
// The det_elem.nominal().worldTransformation() however is DIFFERENT!
int res2 = detail::matrix::_matrixEqual(trafo, m_volMgr.worldTransformation(m_mapping,vid));
if ( res2 != detail::matrix::MATRICES_EQUAL ) {
printout(ERROR, m_det.name(), "VOLUME_PLACEMENT FAILED: World transformation DIFFER.");
++m_place_counters.errors;
}
else {
printout(INFO, m_det.name(), "VOLUME_PLACEMENT: PASSED. All matrices equal: %s",
volumeID(vid).c_str());
}
}
}
}
Expand Down Expand Up @@ -643,6 +646,7 @@ void DetectorCheck::help(int argc,char** argv) {
" sensitive volume placements. \n\n"
" NOTE: Option requires proper PhysVolID setup \n"
" of the sensitive volume placements ! \n"
" -ignore_detector Ignore DetElement placement check for -volmgr \n"
<< std::endl;
std::cout << "Arguments: " << std::endl;
for(int iarg=0; iarg<argc;++iarg) {
Expand All @@ -659,6 +663,7 @@ long DetectorCheck::run(Detector& description,int argc,char** argv) {
bool structure = false;
bool sensitive = false;
bool placements = false;
bool ignore_de = false;
printout(ALWAYS, "DetectorCheck", "++ Processing plugin...");
for(int iarg=0; iarg<argc;++iarg) {
if ( argv[iarg] == 0 ) break;
Expand All @@ -674,6 +679,8 @@ long DetectorCheck::run(Detector& description,int argc,char** argv) {
geometry = true;
else if ( ::strncasecmp(argv[iarg], "-sensitive",4) == 0 )
sensitive = true;
else if ( ::strncasecmp(argv[iarg], "-ignore_detelement",4) == 0 )
ignore_de = true;
else if ( ::strncasecmp(argv[iarg], "-help",4) == 0 )
help(argc, argv);
else
Expand All @@ -685,12 +692,13 @@ long DetectorCheck::run(Detector& description,int argc,char** argv) {
if ( name == "all" || name == "All" || name == "ALL" ) {
for (const auto& det : description.detectors() ) {
printout(INFO, "DetectorCheck", "++ Processing subdetector: %s", det.second.name());
test.check_structure = structure;
test.check_placements = placements;
test.check_volmgr = volmgr;
test.check_geometry = geometry;
test.check_sensitive = sensitive;
test.execute(det.second, 9999);
test.check_structure = structure;
test.check_placements = placements;
test.check_volmgr = volmgr;
test.check_geometry = geometry;
test.check_sensitive = sensitive;
test.ignore_detector = ignore_de;
test.execute(det.second, 9999);
}
return 1;
}
Expand All @@ -702,6 +710,7 @@ long DetectorCheck::run(Detector& description,int argc,char** argv) {
test.check_volmgr = volmgr;
test.check_geometry = geometry;
test.check_sensitive = sensitive;
test.ignore_detector = ignore_de;
test.execute(det, 9999);
}
return 1;
Expand Down
Loading

0 comments on commit b8f2934

Please sign in to comment.