From 8391c6f222ae91a8100b377b89a17951a0209f21 Mon Sep 17 00:00:00 2001
From: Frank Gaede <frank.gaede@desy.de>
Date: Thu, 17 Apr 2014 09:51:11 +0000
Subject: [PATCH] - added functionality to draw surfaces in teve_display -
new method Surface::getVertices() (only planes so far) - added some
convenient methods to Vector3D
---
DDRec/include/DDRec/Surface.h | 12 +-
DDRec/src/Surface.cpp | 129 ++++++++++++++++
DDSurfaces/include/DDSurfaces/Vector3D.h | 25 +++-
UtilityApps/CMakeLists.txt | 8 +-
UtilityApps/src/teve_display.cpp | 180 ++++++++++++++++++++++-
examples/ILDExSimu/src/test_surfaces.cc | 27 +++-
6 files changed, 365 insertions(+), 16 deletions(-)
diff --git a/DDRec/include/DDRec/Surface.h b/DDRec/include/DDRec/Surface.h
index 9f7ba349c..605947e0b 100644
--- a/DDRec/include/DDRec/Surface.h
+++ b/DDRec/include/DDRec/Surface.h
@@ -297,7 +297,6 @@ namespace DD4hep {
Vector3D _n ;
Vector3D _o ;
-
Surface() :_det( Geometry::DetElement() ), _volSurf( VolSurface() ), _wtM( 0 ) , _id( 0) { }
public:
@@ -314,6 +313,11 @@ namespace DD4hep {
*/
virtual const SurfaceType& type() const { return _type ; }
+ Geometry::Volume volume() const { return _volSurf.volume() ; }
+
+ VolSurface volSurface() const { return _volSurf ; }
+
+
//==== geometry ====
/** First direction of measurement U */
@@ -346,6 +350,12 @@ namespace DD4hep {
/// Checks if the given point lies within the surface
virtual bool insideBounds(const Vector3D& point, double epsilon=1.e-4) const ;
+ //---------------------------------------------------
+ /** Get vertices constraining the surface for drawing ( might not be exact boundaries) -
+ * at most nMax points are returned.
+ */
+ std::vector< Vector3D > getVertices( unsigned nMax=360 ) ;
+
protected:
void initialize() ;
diff --git a/DDRec/src/Surface.cpp b/DDRec/src/Surface.cpp
index 02dcf69cf..9f2081ec2 100644
--- a/DDRec/src/Surface.cpp
+++ b/DDRec/src/Surface.cpp
@@ -6,6 +6,8 @@
#include <exception>
#include "TGeoMatrix.h"
+#include "TGeoShape.h"
+#include "TRotation.h"
namespace DD4hep {
namespace DDRec {
@@ -374,6 +376,133 @@ namespace DD4hep {
}
//===================================================================================================================
+ std::vector< Vector3D > Surface::getVertices(unsigned nMax) {
+
+
+ const static double epsilon = 1e-6 ;
+
+ std::vector< Vector3D > _vert ;
+
+
+ // get local and global surface vectors
+ const DDSurfaces::Vector3D& lu = _volSurf.u() ;
+ const DDSurfaces::Vector3D& lv = _volSurf.v() ;
+ const DDSurfaces::Vector3D& ln = _volSurf.normal() ;
+ const DDSurfaces::Vector3D& lo = _volSurf.origin() ;
+
+ Volume vol = volume() ;
+ const TGeoShape* shape = vol->GetShape() ;
+
+
+ if( type().isPlane() ) {
+
+ if( shape->IsA() == TGeoBBox::Class() ) {
+
+ TGeoBBox* box = ( TGeoBBox* ) shape ;
+
+ DDSurfaces::Vector3D boxDim( box->GetDX() , box->GetDY() , box->GetDZ() ) ;
+
+
+ bool isYZ = std::fabs( ln.x() - 1.0 ) < epsilon ; // normal parallel to x
+ bool isXZ = std::fabs( ln.y() - 1.0 ) < epsilon ; // normal parallel to y
+ bool isXY = std::fabs( ln.z() - 1.0 ) < epsilon ; // normal parallel to z
+
+
+ if( isYZ || isXZ || isXY ) { // plane is parallel to one of the box' sides -> need 4 vertices from box dimensions
+
+ // if isYZ :
+ unsigned uidx = 1 ;
+ unsigned vidx = 2 ;
+
+ DDSurfaces::Vector3D ubl( 0., 1., 0. ) ;
+ DDSurfaces::Vector3D vbl( 0., 0., 1. ) ;
+
+ if( isXZ ) {
+
+ ubl.fill( 1., 0., 0. ) ;
+ vbl.fill( 0., 0., 1. ) ;
+ uidx = 0 ;
+ vidx = 2 ;
+
+ } else if( isXY ) {
+
+ ubl.fill( 1., 0., 0. ) ;
+ vbl.fill( 0., 1., 0. ) ;
+ uidx = 0 ;
+ vidx = 1 ;
+ }
+
+ DDSurfaces::Vector3D ub ;
+ DDSurfaces::Vector3D vb ;
+ _wtM->LocalToMasterVect( ubl , ub.array() ) ;
+ _wtM->LocalToMasterVect( vbl , vb.array() ) ;
+
+ _vert.resize(4) ;
+
+ _vert[0] = _o + boxDim[ uidx ] * ub + boxDim[ vidx ] * vb ;
+ _vert[1] = _o - boxDim[ uidx ] * ub + boxDim[ vidx ] * vb ;
+ _vert[2] = _o - boxDim[ uidx ] * ub - boxDim[ vidx ] * vb ;
+ _vert[3] = _o + boxDim[ uidx ] * ub - boxDim[ vidx ] * vb ;
+
+ return _vert ;
+ }
+ }
+
+ // ===== default for arbitrary planes in arbitrary shapes =================
+
+ // We create nMax vertices by rotating the local u vector around the normal
+ // and checking the distance to the volume boundary in that direction.
+ // This is brute force and not very smart, as many points are created on straight
+ // lines and the edges are still rounded.
+ // The alterative would be to compute the true intersections a plane and the most
+ // common shapes - at least for boxes that should be not too hard. To be done...
+
+ _vert.resize( nMax ) ;
+
+ double dAlpha = 2.* ROOT::Math::Pi() / double( nMax ) ;
+
+ TVector3 normal( ln.x() , ln.y() , ln.z() ) ;
+
+ for(unsigned i=0 ; i< nMax ; ++i ){
+
+ double alpha = double(i) * dAlpha ;
+
+ TVector3 vec( lu.x() , lu.y() , lu.z() ) ;
+
+ TRotation rot ;
+ rot.Rotate( alpha , normal );
+
+ TVector3 vecR = rot * vec ;
+
+ DDSurfaces::Vector3D luRot ;
+ luRot.fill( vecR ) ;
+
+ double dist = shape->DistFromInside( lo, luRot , 3, 0.1 ) ;
+
+ // local point at volume boundary
+ DDSurfaces::Vector3D lp = lo + dist * luRot ;
+
+ DDSurfaces::Vector3D gp ;
+
+ _wtM->LocalToMaster( lp , gp.array() ) ;
+
+ // std::cout << " **** normal:" << ln << " lu:" << lu << " alpha:" << alpha << " luRot:" << luRot << " lp :" << lp << " gp:" << gp << std::endl;
+
+ _vert[i] = gp ;
+ }
+
+
+ } else { // cylinder
+
+ // TODO
+ }
+
+ return _vert ;
+
+ }
+
+ //===================================================================================================================
+
} // namespace
} // namespace
diff --git a/DDSurfaces/include/DDSurfaces/Vector3D.h b/DDSurfaces/include/DDSurfaces/Vector3D.h
index 530d73634..dae2cc3f4 100644
--- a/DDSurfaces/include/DDSurfaces/Vector3D.h
+++ b/DDSurfaces/include/DDSurfaces/Vector3D.h
@@ -55,18 +55,28 @@ namespace DDSurfaces {
// }
+ /// fill vector from arbitrary class that defines operator[]
+ template <class T>
+ inline const Vector3D& fill( const T& v ) {
+
+ _x = v[0] ; _y = v[1] ; _z = v[2] ;
+ return *this ;
+ }
+
+ /// fill vector from double array
inline const Vector3D& fill( const double* v) {
_x = v[0] ; _y = v[1] ; _z = v[2] ;
return *this ;
}
+ /// fill from double values
inline const Vector3D& fill( double x, double y, double z) {
-
_x = x ; _y = y ; _z = z ;
return *this ;
}
+
/** Cartesian x coordinate */
inline double x() const { return _x ; }
@@ -86,7 +96,7 @@ namespace DDSurfaces {
inline double& z() { return _z ; }
- /** Accessing x,y,z with bracket operator */
+ /** Accessing x,y,z with bracket operator */
inline double operator[](int i) const {
switch(i) {
case 0: return _x ; break ;
@@ -176,6 +186,11 @@ namespace DDSurfaces {
return &_x ;
}
+ /// direct access to data as double* - allows modification
+ inline double* array() {
+ return &_x ;
+ }
+
/** Implicit templated conversion to anything that has a c'tor T(x,y,z)
* and accessor functions x(),y(),z(). For safety the result is checked which
@@ -249,6 +264,12 @@ namespace DDSurfaces {
return Vector3D( s * v.x() , s * v.y() , s * v.z() ) ;
}
+ /** Negative vector */
+ inline Vector3D operator-( const Vector3D& v) {
+
+ return Vector3D( -v.x(), - v.y(), - v.z() ) ;
+ }
+
/// operator for scalar product
inline double operator*( const Vector3D& v0, const Vector3D& v1 ){
return v0.dot( v1 ) ;
diff --git a/UtilityApps/CMakeLists.txt b/UtilityApps/CMakeLists.txt
index 5517871f9..300988a62 100644
--- a/UtilityApps/CMakeLists.txt
+++ b/UtilityApps/CMakeLists.txt
@@ -2,11 +2,13 @@ cmake_minimum_required(VERSION 2.8.3 FATAL_ERROR)
include_directories( ${CMAKE_SOURCE_DIR}/DDCore/include
${ROOT_INCLUDE_DIR}
- ${CMAKE_SOURCE_DIR}/DDSegmentation/include)
+ ${CMAKE_SOURCE_DIR}/DDSegmentation/include
+ ${CMAKE_SOURCE_DIR}/DDSurfaces/include
+ ${CMAKE_SOURCE_DIR}/DDRec/include )
#-----------------------------------------------------------------------------------
add_executable(geoDisplay src/display.cpp)
-target_link_libraries(geoDisplay DD4hepCore)
+target_link_libraries(geoDisplay DD4hepCore )
#-----------------------------------------------------------------------------------
add_executable(geoConverter src/converter.cpp)
target_link_libraries(geoConverter DD4hepCore)
@@ -15,7 +17,7 @@ add_executable(geoPluginRun src/plugin_runner.cpp)
target_link_libraries(geoPluginRun DD4hepCore)
#-----------------------------------------------------------------------------------
add_executable(teveDisplay src/teve_display.cpp)
-target_link_libraries(teveDisplay DD4hepCore ${ROOT_EVE_LIBRARIES})
+target_link_libraries(teveDisplay DD4hepCore ${ROOT_EVE_LIBRARIES} DD4hepRec )
#--- install target-------------------------------------
diff --git a/UtilityApps/src/teve_display.cpp b/UtilityApps/src/teve_display.cpp
index 5bec5df6d..3a2378cbe 100644
--- a/UtilityApps/src/teve_display.cpp
+++ b/UtilityApps/src/teve_display.cpp
@@ -10,6 +10,9 @@
//====================================================================
#include "DD4hep/Factories.h"
#include "DD4hep/LCDD.h"
+#include "DDRec/SurfaceManager.h"
+
+
#include "run_plugin.h"
#include "TRint.h"
#include "TEveGeoNode.h"
@@ -25,26 +28,194 @@
#define private public
#include "TEveManager.h"
+#include "TEveStraightLineSet.h"
+#include "TRandom.h"
+#include "TGeoShape.h"
+
+using namespace DD4hep ;
+using namespace DDRec ;
+using namespace Geometry ;
+using namespace DDSurfaces ;
+
+//=====================================================================================
+TEveStraightLineSet* drawSurfaceVectors() {
+
+ TEveStraightLineSet* ls = new TEveStraightLineSet("SurfaceVectors");
+
+ LCDD& lcdd = LCDD::getInstance();
+
+ DetElement world = lcdd.world() ;
+
+ // create a list of all surfaces in the detector:
+ SurfaceManager surfMan( world ) ;
+
+ const SurfaceList& sL = surfMan.surfaceList() ;
+
+ for( SurfaceList::const_iterator it = sL.begin() ; it != sL.end() ; ++it ){
+
+ ISurface* surf = *it ;
+ const DDSurfaces::Vector3D& u = surf->u() ;
+ const DDSurfaces::Vector3D& v = surf->v() ;
+ const DDSurfaces::Vector3D& n = surf->normal() ;
+ const DDSurfaces::Vector3D& o = surf->origin() ;
+
+ DDSurfaces::Vector3D ou = o + u ;
+ DDSurfaces::Vector3D ov = o + v ;
+ DDSurfaces::Vector3D on = o + n ;
+
+ ls->AddLine( o.x(), o.y(), o.z(), ou.x() , ou.y() , ou.z() ) ;
+ ls->AddLine( o.x(), o.y(), o.z(), ov.x() , ov.y() , ov.z() ) ;
+ ls->AddLine( o.x(), o.y(), o.z(), on.x() , on.y() , on.z() ) ;
+
+ ls->AddMarker( o.x(), o.y(), o.z() );
+ }
+ ls->SetLineColor( kGreen ) ;
+ ls->SetMarkerColor( kBlue ) ;
+ ls->SetMarkerSize(1);
+ ls->SetMarkerStyle(4);
+
+ gEve->AddElement(ls);
+
+ return ls;
+}
+//=====================================================================================
+TEveStraightLineSet* drawSurfaces() {
+
+ TEveStraightLineSet* ls = new TEveStraightLineSet("Surfaces");
+
+ LCDD& lcdd = LCDD::getInstance();
+
+ DetElement world = lcdd.world() ;
+
+ // create a list of all surfaces in the detector:
+ SurfaceManager surfMan( world ) ;
+
+ const SurfaceList& sL = surfMan.surfaceList() ;
+
+ for( SurfaceList::const_iterator it = sL.begin() ; it != sL.end() ; ++it ){
+
+ Surface* surf = *it ;
+
+ const std::vector< Vector3D > vert = surf->getVertices() ;
+
+ if( vert.empty() )
+ std::cout << " **** drawSurfaces() : empty vertices for surface " << *surf << std::endl ;
+
+ unsigned nV = vert.size() ;
+
+ for( unsigned i=0 ; i<nV ; ++i){
+
+ unsigned j = ( i < nV-1 ? i+1 : 0 ) ;
+
+ ls->AddLine( vert[i].x(), vert[i].y(), vert[i].z(), vert[j].x() , vert[j].y() , vert[j].z() ) ;
+ }
+
+
+ // Volume vol = surf->volume() ;
+ // const TGeoShape* shape = vol->GetShape() ;
+
+ // // get extend of surface volume along u and v
+ // VolSurface volSurf = surf->volSurface() ;
+
+ // const DDSurfaces::Vector3D& lu = volSurf.u() ;
+ // const DDSurfaces::Vector3D& lv = volSurf.v() ;
+ // const DDSurfaces::Vector3D& ln = volSurf.normal() ;
+ // const DDSurfaces::Vector3D& lo = volSurf.origin() ;
+
+
+ // // TGeoShape::DistFromInside(Double_t *point[3], Double_t *dir[3],
+ // // Int_t iact, Double_t step, Double_t *safe)
+
+ // double safe = 0. ;
+ // double lup = shape->DistFromInside( lo, lu, 2, 0.1 , &safe ) ;
+ // double lun = shape->DistFromInside( lo, -lu, 2, 0.1 , &safe ) ;
+ // double lvp = shape->DistFromInside( lo, lv, 2, 0.1 , &safe ) ;
+ // double lvn = shape->DistFromInside( lo, -lv, 2, 0.1 , &safe ) ;
+
+
+ // // std::cout << "*** lu : " << lu << " lv : " << lv << " ln : " << ln << " lo : " << lo << std::endl ;
+ // // std::cout << "*** lup : " << lup << " lun : " << lun << " lvp : " << lvp << " lvn : " << lvn << std::endl ;
+
+ // const DDSurfaces::Vector3D& u = surf->u() ;
+ // const DDSurfaces::Vector3D& v = surf->v() ;
+ // const DDSurfaces::Vector3D& n = surf->normal() ;
+ // const DDSurfaces::Vector3D& o = surf->origin() ;
+
+
+ // if( surf->type().isPlane() ) {
+
+ // DDSurfaces::Vector3D c0 = o + lup * u + lvp * v ;
+ // DDSurfaces::Vector3D c1 = o + lup * u - lvn * v ;
+ // DDSurfaces::Vector3D c2 = o - lun * u - lvn * v ;
+ // DDSurfaces::Vector3D c3 = o - lun * u + lvn * v ;
+
+
+ // ls->AddLine( c0.x(), c0.y(), c0.z(), c1.x() , c1.y() , c1.z() ) ;
+ // ls->AddLine( c1.x(), c1.y(), c1.z(), c2.x() , c2.y() , c2.z() ) ;
+ // ls->AddLine( c2.x(), c2.y(), c2.z(), c3.x() , c3.y() , c3.z() ) ;
+ // ls->AddLine( c3.x(), c3.y(), c3.z(), c0.x() , c0.y() , c0.z() ) ;
+
+ // // ls->AddMarker( o.x(), o.y(), o.z() );
+ // }
+
+ ls->SetLineColor( kRed ) ;
+ ls->SetMarkerColor( kRed ) ;
+ ls->SetMarkerSize(.1);
+ ls->SetMarkerStyle(4);
+
+ } // don't know what to do ...
+
+ gEve->AddElement(ls);
+
+ return ls;
+}
+//=====================================================================================
+
+
static long teve_display(LCDD& lcdd, int /* argc */, char** /* argv */) {
TGeoManager* mgr = &lcdd.manager();
TEveManager::Create();
gEve->fGeometries->Add(new TObjString("DefaultGeometry"),mgr);
TEveGeoTopNode* tn = new TEveGeoTopNode(mgr, mgr->GetTopNode());
tn->SetVisLevel(4);
- gEve->AddGlobalElement(tn);
+
+ // // ---- try to set transparency - does not seem to work ...
+ // TGeoNode* node1 = gGeoManager->GetTopNode();
+ // int k_nodes = node1->GetNdaughters();
+ // for(int i=0; i<k_nodes; i++) {
+ // TGeoNode * CurrentNode = node1->GetDaughter(i);
+ // CurrentNode->GetMedium()->GetMaterial()->SetTransparency(80);
+ // // TEveGeoNode CurrentNode( node1->GetDaughter(i) ) ;
+ // // CurrentNode.SetMainTransparency( 80 ) ;
+ // }
+ // this code seems to have no effect either ...
+ tn->CSCApplyMainTransparencyToAllChildren() ;
+ tn->SetMainTransparency( 80 ) ;
+
+
+ drawSurfaceVectors() ;
+ drawSurfaces() ;
+
+ gEve->AddGlobalElement(tn);
gEve->FullRedraw3D(kTRUE);
+
// EClipType not exported to CINT (see TGLUtil.h):
// 0 - no clip, 1 - clip plane, 2 - clip box
TGLViewer *v = gEve->GetDefaultGLViewer();
- v->GetClipSet()->SetClipType(TGLClip::kClipPlane);
- v->ColorSet().Background().SetColor(kMagenta+4);
+
+ // v->GetClipSet()->SetClipType(TGLClip::kClipPlane);
+ // v->ColorSet().Background().SetColor(kMagenta+4);
+ v->ColorSet().Background().SetColor(kWhite);
+
v->SetGuideState(TGLUtil::kAxesEdge, kTRUE, kFALSE, 0);
v->RefreshPadEditor(v);
- v->CurrentCamera().RotateRad(-1.2, 0.5);
+ // v->CurrentCamera().RotateRad(-1.2, 0.5);
v->DoDraw();
return 1;
+
+
}
DECLARE_APPLY(DD4hepTEveDisplay,teve_display)
@@ -52,3 +223,4 @@ DECLARE_APPLY(DD4hepTEveDisplay,teve_display)
int main(int argc,char** argv) {
return main_default("DD4hepTEveDisplay",argc,argv);
}
+
diff --git a/examples/ILDExSimu/src/test_surfaces.cc b/examples/ILDExSimu/src/test_surfaces.cc
index ff94f0b9c..e5b79cce6 100644
--- a/examples/ILDExSimu/src/test_surfaces.cc
+++ b/examples/ILDExSimu/src/test_surfaces.cc
@@ -84,7 +84,7 @@ int main(int argc, char** argv ){
std::vector< std::string > colNames ;
colNames.push_back( "VXDCollection" ) ;
- colNames.push_back( "SITCollection" ) ;
+ // colNames.push_back( "SITCollection" ) ;
while( ( evt = rdr->readNextEvent() ) != 0 ){
@@ -132,11 +132,15 @@ int main(int argc, char** argv ){
test( isInside , true , sst.str() ) ;
- // std::cout << " found surface with same id " << *surf << std::endl
- // << " point : " << point
- // << " is inside : " << isInside
- // << " distance from surface : " << dist/tgeo::mm << std::endl
- // << std::endl ;
+ if( ! isInside ) {
+
+ std::cout << " found surface " << *surf << std::endl
+ << " id : " << idDecoder
+ << " point : " << point
+ << " is inside : " << isInside
+ << " distance from surface : " << dist/tgeo::mm << std::endl
+ << std::endl ;
+ }
// ====== test that slightly moved hit points are inside their surface ================================
@@ -146,6 +150,17 @@ int main(int argc, char** argv ){
isInside = surf->insideBounds( point2 ) ;
test( isInside , true , sst.str() ) ;
+ if( ! isInside ) {
+
+ std::cout << " found surface " << *surf << std::endl
+ << " id : " << idDecoder
+ << " point : " << point
+ << " is inside : " << isInside
+ << " distance from surface : " << dist/tgeo::mm << std::endl
+ << std::endl ;
+
+ }
+
// ====== test that moved hit points are outside their surface ================================
Vector3D point3 = point + 1e-3 * surf->normal() ;
--
GitLab