From 7dca07b8fa62b3e374612543d30cc53674968f43 Mon Sep 17 00:00:00 2001
From: Frank Gaede <frank.gaede@desy.de>
Date: Tue, 20 May 2014 08:43:08 +0000
Subject: [PATCH] - added new MaterialManager class for access to material
properties - added simple utility print_materials.cc - started to implement
averaging materials for surfaces - issue with ILDExDet geometry and
assemblies ...
---
DDRec/include/DDRec/MaterialManager.h | 77 ++++++++++
DDRec/src/MaterialManager.cpp | 152 +++++++++++++++++++
DDRec/src/Surface.cpp | 194 ++----------------------
UtilityApps/CMakeLists.txt | 7 +-
UtilityApps/src/print_materials.cc | 70 +++++++++
examples/ILDExSimu/src/test_surfaces.cc | 2 +
6 files changed, 322 insertions(+), 180 deletions(-)
create mode 100644 DDRec/include/DDRec/MaterialManager.h
create mode 100644 DDRec/src/MaterialManager.cpp
create mode 100644 UtilityApps/src/print_materials.cc
diff --git a/DDRec/include/DDRec/MaterialManager.h b/DDRec/include/DDRec/MaterialManager.h
new file mode 100644
index 000000000..5626dd2d1
--- /dev/null
+++ b/DDRec/include/DDRec/MaterialManager.h
@@ -0,0 +1,77 @@
+#ifndef DDRec_MaterialManager_H_
+#define DDRec_MaterialManager_H_
+
+#include "DD4hep/Objects.h"
+#include "DDSurfaces/Vector3D.h"
+
+#include <vector>
+
+
+class TGeoManager ;
+
+namespace DD4hep {
+ namespace DDRec {
+
+ typedef Geometry::Material Material ;
+
+ typedef std::vector< std::pair< Material, double > > MaterialVec ;
+
+ /** Material manager provides access to the material properties of the detector.
+ * Material can be accessed either for a given point or as a list of materials along a straight
+ * line between two points.
+ *
+ * @author F.Gaede, DESY
+ * @date May, 19 2014
+ * @version $Id:$
+ */
+ class MaterialManager {
+
+ public:
+
+ MaterialManager();
+
+ ~MaterialManager();
+
+ /** Get a vector with all the materials between the two points p0 and p1 with the corresponding thicknesses -
+ * element type is std::pair< Material, double >.
+ */
+ const MaterialVec& materials(const DDSurfaces::Vector3D& p0, const DDSurfaces::Vector3D& p1 ) ;
+
+ /** Get the material at the given position
+ */
+ const Material& material(const DDSurfaces::Vector3D& pos );
+
+ protected :
+
+ //cached materials
+ MaterialVec _mV ;
+ Material _m ;
+
+ // cached last points
+ DDSurfaces::Vector3D _p0 , _p1, _pos ;
+
+ TGeoManager* _tgeoMgr ;
+ };
+
+
+ /// dump Material operator
+ inline std::ostream& operator<<( std::ostream& os , const Material& m ) {
+ os << " " << m.name() << " Z: " << m.Z() << " A: " << m.A() << " densitiy: " << m.density() << " radiationLength: " << m.radLength()
+ << " interactionLength: " << m.intLength() ;
+ return os ;
+ }
+
+
+ /// dump MaterialVec operator
+ inline std::ostream& operator<<( std::ostream& os , const MaterialVec& m ) {
+
+ for( unsigned i=0,n=m.size() ; i<n ; ++i ) {
+ os << " material: " << m[i].first << " thickness : " << m[i].second << std::endl ;
+ }
+ return os ;
+ }
+
+ } /* namespace DDRec */
+} /* namespace DD4hep */
+
+#endif // DDRec_MaterialManager_H_
diff --git a/DDRec/src/MaterialManager.cpp b/DDRec/src/MaterialManager.cpp
new file mode 100644
index 000000000..ef115088d
--- /dev/null
+++ b/DDRec/src/MaterialManager.cpp
@@ -0,0 +1,152 @@
+#include "DDRec/MaterialManager.h"
+#include "DD4hep/Exceptions.h"
+#include "DD4hep/LCDD.h"
+
+#include "TGeoVolume.h"
+#include "TGeoManager.h"
+#include "TGeoNode.h"
+#include "TVirtualGeoTrack.h"
+
+#define MINSTEP 1.e-5
+
+namespace DD4hep {
+ namespace DDRec {
+
+
+ MaterialManager::MaterialManager() : _mV(0), _m( Material() ), _p0(),_p1(),_pos() {
+
+ _tgeoMgr = Geometry::LCDD::getInstance().world().volume()->GetGeoManager();
+ }
+
+ MaterialManager::~MaterialManager(){
+
+ }
+
+ const MaterialVec&MaterialManager:: materials(const DDSurfaces::Vector3D& p0, const DDSurfaces::Vector3D& p1 ) {
+
+ if( ( p0 != _p0 ) && ( p1 != _p1 ) ) {
+
+ //---------------------------------------
+ _mV.clear() ;
+
+ //
+ // algorithm copied from TGeoGearDistanceProperties.cc :
+ //
+
+ double startpoint[3], endpoint[3], direction[3];
+ double L=0;
+ for(unsigned int i=0; i<3; i++) {
+ startpoint[i] = p0[i];
+ endpoint[i] = p1[i];
+ direction[i] = endpoint[i] - startpoint[i];
+ L+=direction[i]*direction[i];
+ }
+ double totDist = sqrt( L ) ;
+
+ //normalize direction
+ for(unsigned int i=0; i<3; i++)
+ direction[i]=direction[i]/sqrt(L);
+
+ _tgeoMgr->AddTrack(0,12);
+ TGeoNode *node1 = _tgeoMgr->InitTrack(startpoint, direction);
+ //check if there is a node at startpoint
+ if(!node1)
+ throw std::runtime_error("No geometry node found at given location. Either there is no node placed here or position is outside of top volume.");
+
+ while ( !_tgeoMgr->IsOutside() ) {
+ TGeoNode *node2;
+ TVirtualGeoTrack *track;
+ node2 = _tgeoMgr->FindNextBoundaryAndStep(500, 1) ;
+
+ if(!node2 || _tgeoMgr->IsOutside())
+ break;
+
+ const double *position =(double*) _tgeoMgr->GetCurrentPoint();
+ const double *previouspos =(double*) _tgeoMgr->GetLastPoint();
+ double length=_tgeoMgr->GetStep();
+ track = _tgeoMgr->GetLastTrack();
+ //protection against infinitive loop in root which should not happen, but well it does...
+ //work around until solution within root can be found when the step gets very small e.g. 1e-10
+ //and the next boundary is never reached
+
+ if(length<MINSTEP) {
+
+ _tgeoMgr->SetCurrentPoint(position[0]+MINSTEP*direction[0], position[1]+MINSTEP*direction[1], position[2]+MINSTEP*direction[2]);
+ length=_tgeoMgr->GetStep();
+ node2 = _tgeoMgr->FindNextBoundaryAndStep(500, 1) ;
+
+ //fg: need to update positions as well !?
+ position = (const double*) _tgeoMgr->GetCurrentPoint();
+ previouspos = (const double*) _tgeoMgr->GetLastPoint();
+ }
+
+ // printf( " -- step length : %1.8e %1.8e %1.8e %1.8e %1.8e %1.8e %1.8e - %s \n" , length ,
+ // position[0], position[1], position[2], previouspos[0], previouspos[1], previouspos[2] , node1->GetMedium()->GetMaterial()->GetName() ) ;
+
+ DDSurfaces::Vector3D posV( position ) ;
+ double currDistance = ( posV - p0 ).r() ;
+
+ // //if the next boundary is further than end point
+ // if(fabs(position[0])>fabs(endpoint[0]) || fabs(position[1])>fabs(endpoint[1])
+ // || fabs(position[2])>fabs(endpoint[2]))
+
+ //if we travelled too far:
+ if( currDistance > totDist ) {
+
+ length=sqrt(pow(endpoint[0]-previouspos[0],2)
+ +pow(endpoint[1]-previouspos[1],2)
+ +pow(endpoint[2]-previouspos[2],2));
+
+ track->AddPoint(endpoint[0], endpoint[1], endpoint[2], 0.);
+
+
+ _mV.push_back( std::make_pair( Material( node1->GetMedium() ) , length ) ) ;
+
+ break;
+ }
+
+ track->AddPoint(position[0], position[1], position[2], 0.);
+
+ _mV.push_back( std::make_pair( Material( node1->GetMedium() ), length ) ) ;
+
+ node1=node2;
+ }
+
+ _tgeoMgr->ClearTracks();
+ _tgeoMgr->CleanGarbage();
+
+ //---------------------------------------
+
+ _p0 = p0 ;
+ _p1 = p1 ;
+ }
+
+ return _mV ; ;
+ }
+
+
+ const Geometry::Material& MaterialManager::material(const DDSurfaces::Vector3D& pos ){
+
+ if( pos != _pos ) {
+
+ TGeoNode *node=_tgeoMgr->FindNode( pos[0], pos[1], pos[2] ) ;
+
+ if( ! node ) {
+ std::stringstream err ;
+ err << " MaterialManager::material: No geometry node found at location: " << pos ;
+ throw std::runtime_error( err.str() );
+ }
+
+ // std::cout << " @@@ MaterialManager::material @ " << pos << " found volume : " << node->GetName() << std::endl ;
+
+ _m = Material( node->GetMedium() ) ;
+
+ _pos = pos ;
+ }
+
+ return _m ; ;
+ }
+
+
+ } /* namespace DDRec */
+} /* namespace DD4hep */
diff --git a/DDRec/src/Surface.cpp b/DDRec/src/Surface.cpp
index a386b8aa3..a11ddc344 100644
--- a/DDRec/src/Surface.cpp
+++ b/DDRec/src/Surface.cpp
@@ -1,6 +1,8 @@
#include "DDRec/Surface.h"
#include "DD4hep/objects/DetectorInterna.h"
+#include "DDRec/MaterialManager.h"
+
#include <math.h>
#include <memory>
#include <exception>
@@ -223,6 +225,18 @@ namespace DD4hep {
// neeed averaged material in case of several volumes...
// _volSurf.setInnerMaterial( _volSurf.volume().material() ) ;
+ MaterialManager matMgr ;
+ Vector3D p = _o - innerThickness() * _n ;
+
+ const MaterialVec& materials = matMgr.materials( _o , p ) ;
+
+ // std::cout << " ####### found materials between points : " << _o << " and " << p << " : " ;
+ // for( unsigned i=0,n=materials.size();i<n;++i){
+ // std::cout << materials[i].first.name() << "[" << materials[i].second << "], " ;
+ // }
+ // std::cout << std::endl ;
+ //std::cout << " #### material at origin : " << matMgr.material( _o ).name() << " material at endpoint : " << matMgr.material( p ).name() << std::endl ;
+
mat = _volSurf.volume().material() ;
// std::cout << " **** Surface::innerMaterial() - assigning volume material to surface : " << mat.name() << std::endl ;
@@ -243,7 +257,7 @@ namespace DD4hep {
mat = _volSurf.volume().material() ;
- //std::cout << " **** Surface::outerMaterial() - assigning volume material to surface : " << mat.name() << std::endl ;
+ // std::cout << " **** Surface::outerMaterial() - assigning volume material to surface : " << mat.name() << std::endl ;
}
return _volSurf.outerMaterial() ;
@@ -631,184 +645,6 @@ 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.const_array() , luRot.const_array() , 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 if( type().isCylinder() ) {
-
- // // if( shape->IsA() == TGeoTube::Class() ) {
- // if( shape->IsA() == TGeoConeSeg::Class() ) {
-
- // _vert.resize( nMax ) ;
-
- // TGeoTube* tube = ( TGeoTube* ) shape ;
-
- // double zHalf = tube->GetDZ() ;
-
- // Vector3D zv( 0. , 0. , zHalf ) ;
-
- // double r = lo.rho() ;
-
- // unsigned n = nMax / 4 ;
- // n = 8 ;
-
- // double dPhi = 2.* ROOT::Math::Pi() / double( n ) ;
-
- // for( unsigned i = 0 ; i < n ; ++i ) {
-
- // Vector3D rv0( r*sin( i *dPhi ) , r*cos( i *dPhi ) , 0. ) ;
- // Vector3D rv1( r*sin( (i+1)*dPhi ) , r*cos( (i+1)*dPhi ) , 0. ) ;
-
- // // 4 points on local cylinder
-
- // Vector3D pl0 = zv + rv0 ;
- // Vector3D pl1 = zv + rv1 ;
- // Vector3D pl2 = -zv + rv1 ;
- // Vector3D pl3 = -zv + rv0 ;
-
- // Vector3D pg0,pg1,pg2,pg3 ;
-
- // _wtM->LocalToMaster( pl0, pg0.array() ) ;
- // _wtM->LocalToMaster( pl1, pg1.array() ) ;
- // _wtM->LocalToMaster( pl2, pg2.array() ) ;
- // _wtM->LocalToMaster( pl3, pg3.array() ) ;
-
- // _vert[ i*5 ] = pg0 ;
- // _vert[ i*5 + 1 ] = pg1 ;
- // _vert[ i*5 + 2 ] = pg2 ;
- // _vert[ i*5 + 3 ] = pg3 ;
- // _vert[ i*5 + 3 ] = pg0 ;
-
- // }
- // }
- // }
- // return _vert ;
-
- // }
-
- //===================================================================================================================
-
- // double CylinderSurfaceSurface::radius() const {
-
- // return _volSurf.origin().rho() ;
- // }
-
- // //===================================================================================================================
} // namespace
} // namespace
diff --git a/UtilityApps/CMakeLists.txt b/UtilityApps/CMakeLists.txt
index 81a86d252..ab578f37c 100644
--- a/UtilityApps/CMakeLists.txt
+++ b/UtilityApps/CMakeLists.txt
@@ -18,6 +18,11 @@ target_link_libraries(geoConverter DD4hepCore)
add_executable(geoPluginRun src/plugin_runner.cpp)
target_link_libraries(geoPluginRun DD4hepCore)
#-----------------------------------------------------------------------------------
+add_executable( print_materials src/print_materials.cc)
+target_link_libraries(print_materials DD4hepCore DD4hepRec)
+#-----------------------------------------------------------------------------------
+
+
root_generate_dictionary( G__teve src/EvNavHandler.h LINKDEF src/LinkDef.h)
if(DD4HEP_USE_LCIO)
@@ -31,7 +36,7 @@ target_link_libraries( teveDisplay DD4hepCore ${ROOT_EVE_LIBRARIES} DD4hepRec ${
#--- install target-------------------------------------
-install(TARGETS geoDisplay geoConverter geoPluginRun teveDisplay
+install(TARGETS geoDisplay geoConverter geoPluginRun teveDisplay print_materials
RUNTIME DESTINATION bin
LIBRARY DESTINATION lib
)
diff --git a/UtilityApps/src/print_materials.cc b/UtilityApps/src/print_materials.cc
new file mode 100644
index 000000000..8a8a0c74d
--- /dev/null
+++ b/UtilityApps/src/print_materials.cc
@@ -0,0 +1,70 @@
+// $Id:$
+//====================================================================
+// AIDA Detector description implementation for LCD
+//--------------------------------------------------------------------
+//
+// Simple program to print all the materials in a detector on
+// a straight line between two given points
+//
+// Author : F.Gaede, DESY
+//
+//====================================================================
+#include "DD4hep/LCDD.h"
+#include "DD4hep/TGeoUnits.h"
+
+#include "DDRec/MaterialManager.h"
+
+using namespace std ;
+using namespace DD4hep ;
+using namespace DD4hep::Geometry;
+using namespace DD4hep::DDRec;
+using namespace DDSurfaces ;
+using namespace tgeo ;
+
+//=============================================================================
+
+int main(int argc, char** argv ){
+
+ if( argc != 8 ) {
+ std::cout << " usage: print_materials compact.xml x0 y0 z0 x1 y1 z1 " << std::endl
+ << " -> prints the materials on a straight line between the two given points ( unit is cm) " ;
+ exit(1) ;
+ }
+
+ std::string inFile = argv[1] ;
+
+ std::stringstream sstr ;
+ sstr << argv[2] << " " << argv[3] << " " << argv[4] << " " << argv[5] << " " << argv[6] << " " << argv[7] ;
+
+ double x0,y0,z0,x1,y1,z1 ;
+ sstr >> x0 ;
+ sstr >> y0 ;
+ sstr >> z0 ;
+ sstr >> x1 ;
+ sstr >> y1 ;
+ sstr >> z1 ;
+
+ LCDD& lcdd = LCDD::getInstance();
+
+ lcdd.fromCompact( inFile );
+
+ Vector3D p0( x0, y0, z0 ) ;
+ Vector3D p1( x1, y1, z1 ) ;
+
+ MaterialManager matMgr ;
+
+ const MaterialVec& materials = matMgr.materials( p0 , p1 ) ;
+
+ std::cout << std::endl << " ####### materials between the two points : " << p0 << "*cm and " << p1 << "*cm : " << std::endl ;
+
+ for( unsigned i=0,n=materials.size();i<n;++i){
+
+ std::cout << " " << materials[i].first << " \t thickness: " << materials[i].second / tgeo::mm << " mm " << std::endl ;
+ }
+
+ std::cout << "############################################################################### " << std::endl << std::endl ;
+
+ return 0;
+}
+
+//=============================================================================
diff --git a/examples/ILDExSimu/src/test_surfaces.cc b/examples/ILDExSimu/src/test_surfaces.cc
index e5b79cce6..c98e6da17 100644
--- a/examples/ILDExSimu/src/test_surfaces.cc
+++ b/examples/ILDExSimu/src/test_surfaces.cc
@@ -118,6 +118,8 @@ int main(int argc, char** argv ){
if( surf != 0 ){
+ // std::cout << " found surface " << *surf << std::endl ;
+
Vector3D point( sHit->getPosition()[0]* tgeo::mm , sHit->getPosition()[1]* tgeo::mm , sHit->getPosition()[2]* tgeo::mm ) ;
double dist = surf->distance( point ) ;
--
GitLab