diff --git a/DDRec/include/DDRec/Material.h b/DDRec/include/DDRec/Material.h
index 01215ad0d43b4112b7d2e887acbb1fc7390cbacd..0e9d713963789e4e140fba6f64a8fd7b5b39dfb2 100644
--- a/DDRec/include/DDRec/Material.h
+++ b/DDRec/include/DDRec/Material.h
@@ -10,7 +10,8 @@ namespace DD4hep {
namespace DDRec {
- /** Simple data class that implements the DDSurfaces::IMaterial interface.
+ /** Simple data class that implements the DDSurfaces::IMaterial interface
+ * and is used in the Surface implementation.
*
* @author F.Gaede, DESY
* @date May, 20 2014
@@ -28,18 +29,27 @@ namespace DD4hep {
public:
- /** Instantiate from Geometry::Material handle */
- MaterialData( Geometry::Material m ) {
- // handle data can only be accessed if the handle is valid
+ /** Instantiate from Geometry::Material - default initialization if handle is not valid */
+ MaterialData( Geometry::Material m ) :
+
+ _name("unknown"),
+ _Z( -1. ),
+ _A( 0. ),
+ _rho( 0. ),
+ _x0( 0. ),
+ _lambda( 0.) {
+
if( m.isValid() ) {
+
_name= m.name() ;
- _Z= m.Z() ;
- _A= m.A() ;
- _rho= m.density() ;
- _x0= m.radLength() ;
- _lambda= m.intLength() ;
+ _Z = m.Z() ;
+ _A = m.A() ;
+ _rho = m.density() ;
+ _x0 = m.radLength() ;
+ _lambda = m.intLength() ;
+
}
- }
+ }
/** Default c'tor .*/
MaterialData() : _name("unknown"),
@@ -66,8 +76,31 @@ namespace DD4hep {
_lambda( m.interactionLength() ) {}
+ /// assignment from Geometry::Material
+ MaterialData& operator=(const Geometry::Material& m){
+
+ if( m.isValid() ) {
+
+ _name = m.name() ;
+ _Z = m.Z() ;
+ _A = m.A() ;
+ _rho = m.density() ;
+ _x0 = m.radLength() ;
+ _lambda = m.intLength() ;
+
+ } else {
+
+ _name= "unknown";
+ _Z = -1. ;
+ _A = 0. ;
+ _rho = 0. ;
+ _x0 = 0. ;
+ _lambda = 0. ;
+ }
+ }
+
/// true if initialized
- bool isValid() const { return true ; } //( _Z > 0. ) ; }
+ bool isValid() const { return ( _Z > 0. ) ; }
/** D'tor.*/
virtual ~MaterialData() {}
diff --git a/DDRec/include/DDRec/MaterialManager.h b/DDRec/include/DDRec/MaterialManager.h
index a115997dff6c4bcd3bad8a88bdd6780ff301d8ba..4d033da8213a45d64dad8684c3fcafa79886ce69 100644
--- a/DDRec/include/DDRec/MaterialManager.h
+++ b/DDRec/include/DDRec/MaterialManager.h
@@ -4,6 +4,7 @@
#include "DD4hep/Objects.h"
#include "DDSurfaces/Vector3D.h"
#include "DDRec/Material.h"
+#include "DD4hep/TGeoUnits.h"
#include <vector>
@@ -34,16 +35,22 @@ namespace DD4hep {
~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 >.
+ * element type is std::pair< Material, double >. Materials with a thickness smaller than epsilon (default 1e-4=1mu)
+ * are ignored. Avoid calling this method in inner loops as the computation is not cheap. Ideally the result should be cached,
+ * for example as an averaged material @see createAveragedMaterial().
*/
- const MaterialVec& materials(const DDSurfaces::Vector3D& p0, const DDSurfaces::Vector3D& p1 ) ;
+ const MaterialVec& materialsBetween(const DDSurfaces::Vector3D& p0, const DDSurfaces::Vector3D& p1 , double epsilon=1e-4 ) ;
- /** Get the material at the given position
+ /** Get the material at the given position.
*/
- const Material& material(const DDSurfaces::Vector3D& pos );
+ const Material& materialAt(const DDSurfaces::Vector3D& pos );
- void createAveragedMaterial( const MaterialVec& materials, MaterialData& mData ) ;
+ /** Create a material with averaged properties from all materials in the list.
+ * A and Z are averaged by relative number of atoms(molecules), rho is averaged by relative volume
+ * and the inverse radiation and interaction lengths are averaged by relative weight.
+ */
+ MaterialData createAveragedMaterial( const MaterialVec& materials ) ;
protected :
@@ -60,7 +67,8 @@ namespace DD4hep {
/// 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()
+ os << " " << m.name() << " Z: " << m.Z() << " A: " << m.A() << " densitiy: " << m.density()
+ << " radiationLength: " << m.radLength()
<< " interactionLength: " << m.intLength() ;
return os ;
}
@@ -70,7 +78,7 @@ namespace DD4hep {
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 ;
+ os << " material: " << m[i].first << " thickness: " << m[i].second << std::endl ;
}
return os ;
}
diff --git a/DDRec/src/MaterialManager.cpp b/DDRec/src/MaterialManager.cpp
index ef115088d15ec35af952fb053f75653ffa75ba25..90a9deedf50d004eab784ba8c4c98258e61b734e 100644
--- a/DDRec/src/MaterialManager.cpp
+++ b/DDRec/src/MaterialManager.cpp
@@ -22,7 +22,7 @@ namespace DD4hep {
}
- const MaterialVec&MaterialManager:: materials(const DDSurfaces::Vector3D& p0, const DDSurfaces::Vector3D& p1 ) {
+ const MaterialVec&MaterialManager:: materialsBetween(const DDSurfaces::Vector3D& p0, const DDSurfaces::Vector3D& p1 , double epsilon) {
if( ( p0 != _p0 ) && ( p1 != _p1 ) ) {
@@ -30,7 +30,7 @@ namespace DD4hep {
_mV.clear() ;
//
- // algorithm copied from TGeoGearDistanceProperties.cc :
+ // algorithm copied from TGeoGearDistanceProperties.cc (A.Munnich):
//
double startpoint[3], endpoint[3], direction[3];
@@ -100,14 +100,16 @@ namespace DD4hep {
track->AddPoint(endpoint[0], endpoint[1], endpoint[2], 0.);
- _mV.push_back( std::make_pair( Material( node1->GetMedium() ) , length ) ) ;
+ if( length > epsilon )
+ _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 ) ) ;
+ if( length > epsilon )
+ _mV.push_back( std::make_pair( Material( node1->GetMedium() ), length ) ) ;
node1=node2;
}
@@ -125,7 +127,7 @@ namespace DD4hep {
}
- const Geometry::Material& MaterialManager::material(const DDSurfaces::Vector3D& pos ){
+ const Geometry::Material& MaterialManager::materialAt(const DDSurfaces::Vector3D& pos ){
if( pos != _pos ) {
@@ -147,6 +149,60 @@ namespace DD4hep {
return _m ; ;
}
+ MaterialData MaterialManager::createAveragedMaterial( const MaterialVec& materials ) {
+
+ std::stringstream sstr ;
+
+ double sum_l = 0 ;
+ double sum_rho_l = 0 ;
+ double sum_rho_l_over_A = 0 ;
+ double sum_rho_l_Z_over_A = 0 ;
+ //double sum_rho_l_over_x = 0 ;
+ double sum_l_over_x = 0 ;
+ //double sum_rho_l_over_lambda = 0 ;
+ double sum_l_over_lambda = 0 ;
+
+ for(unsigned i=0,n=materials.size(); i<n ; ++i){
+
+ Material mat = materials[i].first ;
+ double l = materials[i].second ;
+
+ if( i != 0 ) sstr << "_" ;
+ sstr << mat.name() << "_" << l ;
+
+ double rho = mat.density() ;
+ double A = mat.A() ;
+ double Z = mat.Z() ;
+ double x = mat.radLength() ;
+ double lambda = mat.intLength() ;
+
+ sum_l += l ;
+ sum_rho_l += rho * l ;
+ sum_rho_l_over_A += rho * l / A ;
+ sum_rho_l_Z_over_A += rho * l * Z / A ;
+ sum_l_over_x += l / x ;
+ sum_l_over_lambda += l / lambda ;
+ // sum_rho_l_over_x += rho * l / x ;
+ // sum_rho_l_over_lambda += rho * l / lambda ;
+ }
+
+ double rho = sum_rho_l / sum_l ;
+
+ double A = sum_rho_l / sum_rho_l_over_A ;
+ double Z = sum_rho_l_Z_over_A / sum_rho_l_over_A ;
+
+ // radiation and interaction lengths already given in cm - average by length
+
+ // double x = sum_rho_l / sum_rho_l_over_x ;
+ double x = sum_l / sum_l_over_x ;
+
+ // double lambda = sum_rho_l / sum_rho_l_over_lambda ;
+ double lambda = sum_l / sum_l_over_lambda ;
+
+
+ return MaterialData( sstr.str() , Z, A, rho, x, lambda ) ;
+
+ }
} /* namespace DDRec */
} /* namespace DD4hep */
diff --git a/DDRec/src/Surface.cpp b/DDRec/src/Surface.cpp
index facff75633e6b4825a309718b2f2b90eb49d85e4..5b02db720dd16418e633026e95f9725eef8a3344 100644
--- a/DDRec/src/Surface.cpp
+++ b/DDRec/src/Surface.cpp
@@ -37,8 +37,8 @@ namespace DD4hep {
_o( Vector3D() ) ,
_th_i( 0. ),
_th_o( 0. ),
- _innerMat( Geometry::Material() ),
- _outerMat( Geometry::Material() ) {
+ _innerMat( MaterialData() ),
+ _outerMat( MaterialData() ) {
}
@@ -50,8 +50,8 @@ namespace DD4hep {
_o( o ),
_th_i( thickness_inner ),
_th_o( thickness_outer ),
- _innerMat( Geometry::Material() ),
- _outerMat( Geometry::Material() ) {
+ _innerMat( MaterialData() ),
+ _outerMat( MaterialData() ) {
}
@@ -220,30 +220,29 @@ namespace DD4hep {
SurfaceMaterial& mat = _volSurf->_innerMat ;
+ // std::cout << " **** Surface::innerMaterial() " << mat << std::endl ;
+
if( ! mat.isValid() ) {
- // fixme: for now just set the material of the volume holding the surface
- // 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 ) ;
+ const MaterialVec& materials = matMgr.materialsBetween( _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 << " ####### 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 ;
+ // const MaterialData& matAvg = matMgr.createAveragedMaterial( materials ) ;
+ // mat = matAvg ;
+ // std::cout << " **** Surface::innerMaterial() - assigning averaged material to surface : " << mat << std::endl ;
- std::cout << " #### material at origin : " << matMgr.material( _o ).name() << " material at endpoint : " << matMgr.material( p ).name() << std::endl ;
-
- mat = _volSurf.volume().material() ;
+ mat = ( materials.size() > 1 ? matMgr.createAveragedMaterial( materials ) : materials[0].first ) ;
- std::cout << " **** Surface::innerMaterial() - assigning volume material to surface : " << mat.name() << std::endl ;
}
- return _volSurf.innerMaterial() ;
+ return mat ;
}
@@ -253,16 +252,16 @@ namespace DD4hep {
if( ! mat.isValid() ) {
- // fixme: for now just set the material of the volume holding the surface
- // neeed averaged material in case of several volumes...
- // _volSurf.setOuterMaterial( _volSurf.volume().material() ) ;
+ MaterialManager matMgr ;
+
+ Vector3D p = _o + outerThickness() * _n ;
+
+ const MaterialVec& materials = matMgr.materialsBetween( _o , p ) ;
- mat = _volSurf.volume().material() ;
+ mat = ( materials.size() > 1 ? matMgr.createAveragedMaterial( materials ) : materials[0].first ) ;
- // std::cout << " **** Surface::outerMaterial() - assigning volume material to surface : " << mat.name() << std::endl ;
}
-
- return _volSurf.outerMaterial() ;
+ return mat ;
}
double Surface::distance(const Vector3D& point ) const {
diff --git a/UtilityApps/src/print_materials.cc b/UtilityApps/src/print_materials.cc
index 8a8a0c74de7a8b5525ec3caac888992c151b2b48..9d10af238ee4a0fe3c71e28108496764b68a57ca 100644
--- a/UtilityApps/src/print_materials.cc
+++ b/UtilityApps/src/print_materials.cc
@@ -53,17 +53,48 @@ int main(int argc, char** argv ){
MaterialManager matMgr ;
- const MaterialVec& materials = matMgr.materials( p0 , p1 ) ;
+ const MaterialVec& materials = matMgr.materialsBetween( p0 , p1 ) ;
std::cout << std::endl << " ####### materials between the two points : " << p0 << "*cm and " << p1 << "*cm : " << std::endl ;
+ double sum_x0 = 0 ;
+ double sum_lambda = 0 ;
+ double path_length = 0 ;
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 ;
+ Material mat = materials[i].first ;
+ double length = materials[i].second ;
+
+ double nx0 = length / mat.radLength() ;
+ sum_x0 += nx0 ;
+
+ double nLambda = length / mat.intLength() ;
+ sum_lambda += nLambda ;
+
+ path_length += length ;
+
+ std::cout << " " << mat << " thickness: " << length << " path_length:" << path_length << " integrated_X0: " << sum_x0 << " integrated_lambda: " << sum_lambda << std::endl ;
}
std::cout << "############################################################################### " << std::endl << std::endl ;
+
+ const MaterialData& avMat = matMgr.createAveragedMaterial( materials ) ;
+
+ std::cout << " averaged Material : " << " Z: " << avMat.Z() << " A: " << avMat.A() << " densitiy: " << avMat.density()
+ << " radiationLength: " << avMat.radiationLength()
+ << " interactionLength: " << avMat.interactionLength() << std::endl << std::endl ;
+
+
+ std::cout << " Total length : " << path_length / tgeo::mm << " mm " << std::endl ;
+
+ std::cout << " Integrated radiation lengths : " << path_length / avMat.radiationLength() << " X0 " << std::endl ;
+
+ std::cout << " Integrated interaction lengths : " << path_length / avMat.interactionLength() << " lambda " << std::endl ;
+
+ std::cout << "############################################################################### " << std::endl << std::endl ;
+
+
return 0;
}
diff --git a/doc/release.notes b/doc/release.notes
index ccda81c525348ab596a39eb65008bf748df58260..3f4948d9ea22a55a5523ea556142d229724be6f9 100644
--- a/doc/release.notes
+++ b/doc/release.notes
@@ -1,6 +1,32 @@
DD4hep ---- Release Notes
=================================
+2014/05/21 Frank Gaede
+-----------------------
+ - add MaterialManager class providing
+ - access to materials at any point or on straight
+ line between two points
+ - creation of material with averaged properties (A,Z,rho,x0,Lambda)
+
+ - added utility print_materials to print material properties along
+ a straight line between two points including integrated radiation and
+ interaction lengths (useful for debugging geometries and materials)
+
+ - use avaeraged material for Surfaces where the thickness extends
+ beyond the volume boundaries
+ - introduced new simple data class MaterialData for this
+
+
+ *** known issues
+ ** materials don't work for detectors with assemblies in assemblies
+ as the TGeo navigation dose not seem to work:
+ Error in <TGeoVoxelFinder::SortAll>: Wrong bounding box for volume SIT_assembly
+ -> ROOT bug or feature ?
+
+ ** using <composite/> in compound materials results in incorrect material properties
+ see ILDExDet/compact/materials Polysterene as example
+
+
2014/05/06 Frank Gaede
-----------------------
- DDSurfaces/DDRec:
diff --git a/examples/ILDExDet/compact/materials.xml b/examples/ILDExDet/compact/materials.xml
index 907706d5ae8c8282fba14315acc43eeb8d6d9710..59499dc2f266a06d7eafb2d6c58b79ed0bca5747 100644
--- a/examples/ILDExDet/compact/materials.xml
+++ b/examples/ILDExDet/compact/materials.xml
@@ -48,6 +48,19 @@
</material>
<material name="Polystyrene">
+ <D value="1.032" unit="g/cm3"/>
+ <fraction n="0.077418" ref="H"/>
+ <fraction n="0.922582" ref="C"/>
+ </material>
+
+ <!-- FIXME: definition with composites results
+ in invalid material properties:
+Polystyrene Z: 135 A: 249.37 densitiy: 1.032 radiationLength: 1.19602 interactionLength: 1.1644
+ - the above results in:
+Polystyrene Z: 5.61291 A: 11.1589 densitiy: 1.032 radiationLength: 42.0422 interactionLength: 70.7405
+
+ -->
+ <material name="Polystyrene_BAD">
<D value="1.032" unit="g/cm3"/>
<composite n="19" ref="C"/>
<composite n="21" ref="H"/>