From d3019abf890cee7736f809f1c88f8585b0bf1952 Mon Sep 17 00:00:00 2001
From: Markus Frank <Markus.Frank@cern.ch>
Date: Tue, 16 Jan 2024 13:52:34 +0100
Subject: [PATCH] Fix facetIsDegenerated utility function
---
DDCAD/CMakeLists.txt | 4 +-
DDCAD/include/DDCAD/Utilities.h | 15 +-
examples/AlignDet/src/AlephTPC_geo.cpp | 282 ++++++++++++-------------
3 files changed, 151 insertions(+), 150 deletions(-)
diff --git a/DDCAD/CMakeLists.txt b/DDCAD/CMakeLists.txt
index a85572aff..8e9df077d 100644
--- a/DDCAD/CMakeLists.txt
+++ b/DDCAD/CMakeLists.txt
@@ -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}
diff --git a/DDCAD/include/DDCAD/Utilities.h b/DDCAD/include/DDCAD/Utilities.h
index 6e12f0401..effbbeb5d 100644
--- a/DDCAD/include/DDCAD/Utilities.h
+++ b/DDCAD/include/DDCAD/Utilities.h
@@ -14,6 +14,7 @@
#define DDCAD_UTILITIES_H
#include <vector>
+#include <cmath>
#include <TGeoTessellated.h>
#include <TGeoVector3.h>
@@ -63,13 +64,13 @@ namespace dd4hep {
//
// 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;
+ + 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 std::fabs(det) < epsilon;
}
}
}
diff --git a/examples/AlignDet/src/AlephTPC_geo.cpp b/examples/AlignDet/src/AlephTPC_geo.cpp
index 0c0fe041d..b4a089d03 100644
--- a/examples/AlignDet/src/AlephTPC_geo.cpp
+++ b/examples/AlignDet/src/AlephTPC_geo.cpp
@@ -164,133 +164,133 @@ static Ref_t create_element(Detector& description, xml_h e, SensitiveDetector se
if ( layer_vis.empty() ) layer_vis = sector_vis;
if ( sector_type == 'K' ) {
- double angle = M_PI/12.0;
- double angle1 = std::tan(angle);
- double v[8][2];
- v[0][0] = rmin;
- v[0][1] = 0;
- v[1][0] = rmin;
- v[1][1] = rmin*angle1;
- v[2][0] = rmax;
- v[2][1] = rmax*angle1;
- v[3][0] = rmax;
- v[3][1] = 0;
- ::memcpy(&v[4][0],&v[0][0],8*sizeof(double));
- EightPointSolid upper(layer_thickness/2,&v[0][0]);
-
- v[0][0] = rmin;
- v[0][1] = 0;
- v[1][0] = rmax;
- v[1][1] = 0;
- v[2][0] = rmax;
- v[2][1] = -rmax*angle1;
- v[3][0] = rmin;
- v[3][1] = -rmin*angle1;
- ::memcpy(&v[4][0],&v[0][0],8*sizeof(double));
- EightPointSolid lower(layer_thickness/2,&v[0][0]);
-
- v[0][0] = rmin;
- v[0][1] = gap_half;
- v[1][0] = rmin;
- v[1][1] = rmin*angle1;
- v[2][0] = rmax;
- v[2][1] = rmax*angle1;
- v[3][0] = rmax;
- v[3][1] = gap_half;
- ::memcpy(&v[4][0],&v[0][0],8*sizeof(double));
- EightPointSolid top(layer_thickness/2,&v[0][0]);
-
- v[0][0] = rmin;
- v[0][1] = -gap_half;
- v[1][0] = rmax;
- v[1][1] = -gap_half;
- v[2][0] = rmax;
- v[2][1] = -rmax*angle1;
- v[3][0] = rmin;
- v[3][1] = -rmin*angle1;
- ::memcpy(&v[4][0],&v[0][0],8*sizeof(double));
- EightPointSolid bottom(layer_thickness/2,&v[0][0]);
-
- UnionSolid u(UnionSolid(upper,lower),top,Rotation3D(RotationZ(2*(angle))));
- tm = UnionSolid(u,bottom,Rotation3D(RotationZ(-2*(angle))));
+ double angle = M_PI/12.0;
+ double angle1 = std::tan(angle);
+ double v[8][2];
+ v[0][0] = rmin;
+ v[0][1] = 0;
+ v[1][0] = rmin;
+ v[1][1] = rmin*angle1;
+ v[2][0] = rmax;
+ v[2][1] = rmax*angle1;
+ v[3][0] = rmax;
+ v[3][1] = 0;
+ ::memcpy(&v[4][0],&v[0][0],8*sizeof(double));
+ EightPointSolid upper(layer_thickness/2,&v[0][0]);
+
+ v[0][0] = rmin;
+ v[0][1] = 0;
+ v[1][0] = rmax;
+ v[1][1] = 0;
+ v[2][0] = rmax;
+ v[2][1] = -rmax*angle1;
+ v[3][0] = rmin;
+ v[3][1] = -rmin*angle1;
+ ::memcpy(&v[4][0],&v[0][0],8*sizeof(double));
+ EightPointSolid lower(layer_thickness/2,&v[0][0]);
+
+ v[0][0] = rmin;
+ v[0][1] = gap_half;
+ v[1][0] = rmin;
+ v[1][1] = rmin*angle1;
+ v[2][0] = rmax;
+ v[2][1] = rmax*angle1;
+ v[3][0] = rmax;
+ v[3][1] = gap_half;
+ ::memcpy(&v[4][0],&v[0][0],8*sizeof(double));
+ EightPointSolid top(layer_thickness/2,&v[0][0]);
+
+ v[0][0] = rmin;
+ v[0][1] = -gap_half;
+ v[1][0] = rmax;
+ v[1][1] = -gap_half;
+ v[2][0] = rmax;
+ v[2][1] = -rmax*angle1;
+ v[3][0] = rmin;
+ v[3][1] = -rmin*angle1;
+ ::memcpy(&v[4][0],&v[0][0],8*sizeof(double));
+ EightPointSolid bottom(layer_thickness/2,&v[0][0]);
+
+ UnionSolid u(UnionSolid(upper,lower),top,Rotation3D(RotationZ(2*(angle))));
+ tm = UnionSolid(u,bottom,Rotation3D(RotationZ(-2*(angle))));
}
else {
- double overlap = sector_type=='W' ? 20 : -20;
- double angle = M_PI/12.0;
- double angle1 = std::tan(angle);
- double v[8][2], dr = 0;
-
- if ( sector_type == 'W' ) {
- rmax += overlap*std::tan(angle/2);
- dr = overlap*std::tan(angle/2);
- }
-
- v[0][0] = rmin;
- v[0][1] = -rmin*angle1+gap_half;
- v[1][0] = rmin;
- v[1][1] = rmin*angle1-gap_half;
- v[2][0] = rmax;
- v[2][1] = rmax*angle1-gap_half;
- v[3][0] = rmax;
- v[3][1] = -rmax*angle1+gap_half;
- ::memcpy(&v[4][0],&v[0][0],8*sizeof(double));
- printCoordinates(sector_type,"t2:",v);
- EightPointSolid sectorSolid(layer_thickness/2,&v[0][0]);
-
- if ( sector_type=='W' ) {
- v[0][0] = (rmax+rmin)/2-dr;
- v[0][1] = (rmax+rmin)/2*angle1-gap_half;
- v[1][0] = rmax+0.0001;
- v[1][1] = rmax*angle1-gap_half;
- v[2][0] = rmax+0.0001;
- v[2][1] = (rmax-overlap)*angle1-gap_half;
- v[3][0] = (rmax+rmin)/2-dr;
- v[3][1] = (rmax+rmin-overlap)/2*angle1-gap_half;
- }
- else {
- v[0][0] = (rmax+rmin)/2-dr;
- v[0][1] = (rmax+rmin)/2*angle1-gap_half;
- v[3][0] = rmax+0.0001;
- v[3][1] = rmax*angle1-gap_half;
- v[2][0] = rmax+0.0001;
- v[2][1] = (rmax-overlap)*angle1-gap_half;
- v[1][0] = (rmax+rmin)/2-dr;
- v[1][1] = (rmax+rmin-overlap)/2*angle1-gap_half;
- }
- printCoordinates(sector_type,"upper_boolean:",v);
- ::memcpy(&v[4][0],&v[0][0],8*sizeof(double));
- EightPointSolid upper_boolean((sector_type == 'W' ? 1.1 : 1.0) * layer_thickness/2,&v[0][0]);
-
- if ( sector_type=='W' ) {
- v[0][0] = (rmax+rmin)/2-dr;
- v[0][1] = -((rmax+rmin)/2*angle1-gap_half);
- v[1][0] = (rmax+rmin)/2-dr;
- v[1][1] = -((rmax+rmin-overlap)/2*angle1-gap_half);
- v[2][0] = rmax+0.0001;
- v[2][1] = -((rmax-overlap)*angle1-gap_half);
- v[3][0] = rmax+0.0001;
- v[3][1] = -(rmax*angle1-gap_half);
- }
- else {
- v[0][0] = (rmax+rmin)/2-dr;
- v[0][1] = -((rmax+rmin)/2*angle1-gap_half);
- v[3][0] = (rmax+rmin)/2-dr;
- v[3][1] = -((rmax+rmin-overlap)/2*angle1-gap_half);
- v[2][0] = rmax+0.0001;
- v[2][1] = -((rmax-overlap)*angle1-gap_half);
- v[1][0] = rmax+0.0001;
- v[1][1] = -(rmax*angle1-gap_half);
- }
- printCoordinates(sector_type,"lower_boolean:",v);
- ::memcpy(&v[4][0],&v[0][0],8*sizeof(double));
-
- // For W sectors make the subtraction solid slightly thicker to ensure everything is cut off
- // and no left-overs from numerical precision are left.
- EightPointSolid lower_boolean((sector_type == 'W' ? 1.1 : 1.0) * layer_thickness/2,&v[0][0]);
- if ( sector_type == 'W' )
- tm = SubtractionSolid(SubtractionSolid(sectorSolid,upper_boolean),lower_boolean);
- else
- tm = UnionSolid(UnionSolid(sectorSolid,upper_boolean),lower_boolean);
+ double overlap = sector_type=='W' ? 20 : -20;
+ double angle = M_PI/12.0;
+ double angle1 = std::tan(angle);
+ double v[8][2], dr = 0;
+
+ if ( sector_type == 'W' ) {
+ rmax += overlap*std::tan(angle/2);
+ dr = overlap*std::tan(angle/2);
+ }
+
+ v[0][0] = rmin;
+ v[0][1] = -rmin*angle1+gap_half;
+ v[1][0] = rmin;
+ v[1][1] = rmin*angle1-gap_half;
+ v[2][0] = rmax;
+ v[2][1] = rmax*angle1-gap_half;
+ v[3][0] = rmax;
+ v[3][1] = -rmax*angle1+gap_half;
+ ::memcpy(&v[4][0],&v[0][0],8*sizeof(double));
+ printCoordinates(sector_type,"t2:",v);
+ EightPointSolid sectorSolid(layer_thickness/2,&v[0][0]);
+
+ if ( sector_type=='W' ) {
+ v[0][0] = (rmax+rmin)/2-dr;
+ v[0][1] = (rmax+rmin)/2*angle1-gap_half;
+ v[1][0] = rmax+0.0001;
+ v[1][1] = rmax*angle1-gap_half;
+ v[2][0] = rmax+0.0001;
+ v[2][1] = (rmax-overlap)*angle1-gap_half;
+ v[3][0] = (rmax+rmin)/2-dr;
+ v[3][1] = (rmax+rmin-overlap)/2*angle1-gap_half;
+ }
+ else {
+ v[0][0] = (rmax+rmin)/2-dr;
+ v[0][1] = (rmax+rmin)/2*angle1-gap_half;
+ v[3][0] = rmax+0.0001;
+ v[3][1] = rmax*angle1-gap_half;
+ v[2][0] = rmax+0.0001;
+ v[2][1] = (rmax-overlap)*angle1-gap_half;
+ v[1][0] = (rmax+rmin)/2-dr;
+ v[1][1] = (rmax+rmin-overlap)/2*angle1-gap_half;
+ }
+ printCoordinates(sector_type,"upper_boolean:",v);
+ ::memcpy(&v[4][0],&v[0][0],8*sizeof(double));
+ EightPointSolid upper_boolean((sector_type == 'W' ? 1.1 : 1.0) * layer_thickness/2,&v[0][0]);
+
+ if ( sector_type=='W' ) {
+ v[0][0] = (rmax+rmin)/2-dr;
+ v[0][1] = -((rmax+rmin)/2*angle1-gap_half);
+ v[1][0] = (rmax+rmin)/2-dr;
+ v[1][1] = -((rmax+rmin-overlap)/2*angle1-gap_half);
+ v[2][0] = rmax+0.0001;
+ v[2][1] = -((rmax-overlap)*angle1-gap_half);
+ v[3][0] = rmax+0.0001;
+ v[3][1] = -(rmax*angle1-gap_half);
+ }
+ else {
+ v[0][0] = (rmax+rmin)/2-dr;
+ v[0][1] = -((rmax+rmin)/2*angle1-gap_half);
+ v[3][0] = (rmax+rmin)/2-dr;
+ v[3][1] = -((rmax+rmin-overlap)/2*angle1-gap_half);
+ v[2][0] = rmax+0.0001;
+ v[2][1] = -((rmax-overlap)*angle1-gap_half);
+ v[1][0] = rmax+0.0001;
+ v[1][1] = -(rmax*angle1-gap_half);
+ }
+ printCoordinates(sector_type,"lower_boolean:",v);
+ ::memcpy(&v[4][0],&v[0][0],8*sizeof(double));
+
+ // For W sectors make the subtraction solid slightly thicker to ensure everything is cut off
+ // and no left-overs from numerical precision are left.
+ EightPointSolid lower_boolean((sector_type == 'W' ? 1.1 : 1.0) * layer_thickness/2,&v[0][0]);
+ if ( sector_type == 'W' )
+ tm = SubtractionSolid(SubtractionSolid(sectorSolid,upper_boolean),lower_boolean);
+ else
+ tm = UnionSolid(UnionSolid(sectorSolid,upper_boolean),lower_boolean);
}
Volume secVol(name+"_sector_"+sector_type+_toString(i_layer,"_layer%d"),tm,description.material(layer_mat));
@@ -304,22 +304,22 @@ static Ref_t create_element(Detector& description, xml_h e, SensitiveDetector se
int j = i + (sector_type=='W' ? 1 : 0);
double phi = dphi*j+shift + (sector_type=='K' ? 0 : M_PI/12.0);
if ( sector_type == 'K' || (i%2)==0 ) {
- Transform3D trA(RotationZYX(phi,0,0),Position(0,0,0.00001));
- Transform3D trB(RotationZYX(phi,0,0),Position(0,0,0.00001));
-
- pv = endCapAVol.placeVolume(sector,trA);
- pv.addPhysVolID("type", sector_type=='K' ? 1 : sector_type=='M' ? 2 : 3);
- pv.addPhysVolID("sector",j);
- DetElement sectorA(detA,detA.name()+_toString(sector_count,"_sector%02d"),1);
- sectorA.setPlacement(pv);
-
- pv = endCapBVol.placeVolume(sector,trB);
- pv.addPhysVolID("type", sector_type=='K' ? 1 : sector_type=='M' ? 2 : 3);
- pv.addPhysVolID("sector",j);
- DetElement sectorB(detB,detB.name()+_toString(sector_count,"_sector%02d"),1);
- sectorB.setPlacement(pv);
- cout << "Placed " << sector_type << " sector at phi=" << phi << endl;
- ++sector_count;
+ Transform3D trA(RotationZYX(phi,0,0),Position(0,0,0.00001));
+ Transform3D trB(RotationZYX(phi,0,0),Position(0,0,0.00001));
+
+ pv = endCapAVol.placeVolume(sector,trA);
+ pv.addPhysVolID("type", sector_type=='K' ? 1 : sector_type=='M' ? 2 : 3);
+ pv.addPhysVolID("sector",j);
+ DetElement sectorA(detA,detA.name()+_toString(sector_count,"_sector%02d"),1);
+ sectorA.setPlacement(pv);
+
+ pv = endCapBVol.placeVolume(sector,trB);
+ pv.addPhysVolID("type", sector_type=='K' ? 1 : sector_type=='M' ? 2 : 3);
+ pv.addPhysVolID("sector",j);
+ DetElement sectorB(detB,detB.name()+_toString(sector_count,"_sector%02d"),1);
+ sectorB.setPlacement(pv);
+ cout << "Placed " << sector_type << " sector at phi=" << phi << endl;
+ ++sector_count;
}
}
}
--
GitLab