From feec9d47e04e8cd85c0b4851313eec5ba6804b16 Mon Sep 17 00:00:00 2001
From: Frank Gaede <frank.gaede@desy.de>
Date: Fri, 11 Dec 2015 19:15:35 +0000
Subject: [PATCH] - updated cone implmentation in Surfaces - fixed module id
for calo barrel surfaces - ...
---
DDDetectors/src/CaloFaceBarrel_surfaces.cpp | 4 +-
DDDetectors/src/CaloFaceEndcap_surfaces.cpp | 4 +-
DDRec/include/DDRec/Surface.h | 6 +++
DDRec/src/Surface.cpp | 53 +++++++++++++++++----
DDSegmentation/src/Segmentation.cpp | 1 +
DDSegmentation/src/TiledLayerGridXY.cpp | 4 +-
DDSurfaces/include/DDSurfaces/ISurface.h | 6 +++
7 files changed, 64 insertions(+), 14 deletions(-)
diff --git a/DDDetectors/src/CaloFaceBarrel_surfaces.cpp b/DDDetectors/src/CaloFaceBarrel_surfaces.cpp
index 59fb56de5..59daa32f2 100644
--- a/DDDetectors/src/CaloFaceBarrel_surfaces.cpp
+++ b/DDDetectors/src/CaloFaceBarrel_surfaces.cpp
@@ -97,7 +97,7 @@ namespace{
else if( name=="phi0" ) data.phi0 = value ;
else if( name=="symmetry") data.symmetry = value ;
else if( name=="systemID") data.systemID = value ;
- else if( name=="encoding") data.encoding = value ;
+ else if( name=="encoding") data.encoding = ptr ;
else {
std::cout << "DD4hep_CaloFaceBarrelSurfacePlugin: WARNING unknown parameter: " << name << std::endl ;
}
@@ -138,7 +138,7 @@ namespace{
for(unsigned i=0 ; i < symmetry ; ++i){
- bf["module"] = 1 ;
+ bf["module"] = i ;
double gam = phi0 + alpha/2. + i*alpha;
diff --git a/DDDetectors/src/CaloFaceEndcap_surfaces.cpp b/DDDetectors/src/CaloFaceEndcap_surfaces.cpp
index 13fa9c729..d48925899 100644
--- a/DDDetectors/src/CaloFaceEndcap_surfaces.cpp
+++ b/DDDetectors/src/CaloFaceEndcap_surfaces.cpp
@@ -102,7 +102,7 @@ namespace{
else if( name=="phi0" ) data.phi0 = value ;
else if( name=="symmetry") data.symmetry = value ;
else if( name=="systemID") data.systemID = value ;
- else if( name=="encoding") data.encoding = value ;
+ else if( name=="encoding") data.encoding = ptr ;
else {
std::cout << "DD4hep_CaloFaceEndcapSurfacePlugin: WARNING unknown parameter: " << name << std::endl ;
}
@@ -133,7 +133,7 @@ namespace{
DD4hep::DDSegmentation::BitField64 bf( "system:5,side:-2,layer:9,module:8,sensor:8" ) ;
- // DD4hep::DDSegmentation::BitField64 bf( data.encoding ) ;
+ //DD4hep::DDSegmentation::BitField64 bf( data.encoding ) ;
bf["system"] = data.systemID ;
bf["side"] = 1 ;
diff --git a/DDRec/include/DDRec/Surface.h b/DDRec/include/DDRec/Surface.h
index b4e484907..2b44c5f8f 100644
--- a/DDRec/include/DDRec/Surface.h
+++ b/DDRec/include/DDRec/Surface.h
@@ -413,6 +413,12 @@ namespace DD4hep {
*/
class VolConeImpl : public VolSurfaceBase {
+ //internal helper variables
+ double _ztip ; // z position of the tip in the volume coordinate system
+ double _zt0 ; // z distance of the front face from the tip
+ double _zt1 ; // z distance of the back face from the tip
+ double _tanTheta ; // tan of half the openeing angle
+
public:
/// default c'tor
diff --git a/DDRec/src/Surface.cpp b/DDRec/src/Surface.cpp
index 74375d750..6741386f5 100644
--- a/DDRec/src/Surface.cpp
+++ b/DDRec/src/Surface.cpp
@@ -336,13 +336,31 @@ namespace DD4hep {
throw std::runtime_error( sst.str() ) ;
}
- Vector3D n( 1. , v_val.phi() , ( v_val.theta() + M_PI/2. ) , Vector3D::spherical ) ;
+ double theta = v_val.theta() ;
+
+ Vector3D n( 1. , v_val.phi() , ( theta + M_PI/2. ) , Vector3D::spherical ) ;
Vector3D u_val = v_val.cross( n ) ;
setU( u_val ) ;
setOrigin( o_rphi ) ;
setNormal( n ) ;
+ // set helper variable for faster computations (describe cone with tip at origin)
+ _tanTheta = std::tan( theta ) ;
+ double tipoffset = o_val.rho() / _tanTheta ; // distance from tip to origin.z()
+ _ztip = o_val.z() - tipoffset ;
+
+ double dist_p = vol->GetShape()->DistFromInside( const_cast<double*> ( o_val.const_array() ) ,
+ const_cast<double*> ( v_val.const_array() ) ) ;
+ Vector3D vm = -1. * v_val ;
+ double dist_m = vol->GetShape()->DistFromInside( const_cast<double*> ( o_val.const_array() ) ,
+ const_cast<double*> ( vm.array() ) ) ;
+
+ double costh = std::cos( theta) ;
+ _zt0 = tipoffset - dist_m * costh ;
+ _zt1 = tipoffset + dist_p * costh ;
+
+
_type.setProperty( SurfaceType::Plane , false ) ;
_type.setProperty( SurfaceType::Cylinder, false ) ;
_type.setProperty( SurfaceType::Cone , true ) ;
@@ -378,7 +396,10 @@ namespace DD4hep {
while( phi < -M_PI ) phi += 2.*M_PI ;
while( phi > M_PI ) phi -= 2.*M_PI ;
- return Vector2D( origin().rho()*phi, ( point.z() - origin().z() ) / cos( _v.theta() ) ) ;
+
+ double r = ( point.z() - _ztip ) * _tanTheta ;
+
+ return Vector2D( r*phi, ( point.z() - origin().z() ) / cos( _v.theta() ) ) ;
}
@@ -386,26 +407,40 @@ namespace DD4hep {
double z = point.v() * cos( _v.theta() ) + origin().z() ;
- double phi = point.u() / origin().rho() + origin().phi() ;
-
+ double r = ( z - _ztip ) * _tanTheta ;
+
+ double phi = point.u() / r + origin().phi() ;
+
while( phi < -M_PI ) phi += 2.*M_PI ;
while( phi > M_PI ) phi -= 2.*M_PI ;
- return Vector3D( origin().rho() , phi, z , Vector3D::cylindrical ) ;
+ return Vector3D( r , phi, z , Vector3D::cylindrical ) ;
}
/** Distance to surface */
double VolConeImpl::distance(const Vector3D& point ) const {
+ // // if the point is in the other hemispere we return the distance to origin
+ // // -> this assumes that the cones do not cross the xy-plane ...
+ // // otherwise we get the distance to the mirrored part of the cone
+ // // needs more thought ..
+ // if( origin().z() * point.z() < 0. )
+ // return point.r() ;
+
//fixme: there are probably faster ways to compute this
// e.g by using the intercept theorem - tbd. ...
- const Vector2D& lp = globalToLocal( point ) ;
- const Vector3D& gp = localToGlobal( lp ) ;
+ // const Vector2D& lp = globalToLocal( point ) ;
+ // const Vector3D& gp = localToGlobal( lp ) ;
+
+ // Vector3D dz = point - gp ;
+
+ //return dz * normal( point ) ;
- Vector3D dz = point - gp ;
+ double zp = point.z() - _ztip ;
+ double r = point.rho() - zp * _tanTheta ;
+ return r * std::cos( _v.theta() ) ;
- return dz * normal( point ) ;
}
/// create outer bounding lines for the given symmetry of the polyhedron
diff --git a/DDSegmentation/src/Segmentation.cpp b/DDSegmentation/src/Segmentation.cpp
index dfabc427c..27be7fc59 100644
--- a/DDSegmentation/src/Segmentation.cpp
+++ b/DDSegmentation/src/Segmentation.cpp
@@ -27,6 +27,7 @@ using std::vector;
/// Default constructor used by derived classes passing the encoding string
Segmentation::Segmentation(const std::string& cellEncoding) :
_name("Segmentation"), _type("Segmentation"), _decoder(new BitField64(cellEncoding)), _ownsDecoder(true) {
+
}
/// Default constructor used by derived classes passing an existing decoder
diff --git a/DDSegmentation/src/TiledLayerGridXY.cpp b/DDSegmentation/src/TiledLayerGridXY.cpp
index 6070a084e..0c65656b9 100644
--- a/DDSegmentation/src/TiledLayerGridXY.cpp
+++ b/DDSegmentation/src/TiledLayerGridXY.cpp
@@ -21,7 +21,9 @@ TiledLayerGridXY::TiledLayerGridXY(const std::string& cellEncoding) :
// define type and description
_type = "TiledLayerGridXY";
_description = "Cartesian segmentation in the local XY-plane using optimal tiling depending on the layer dimensions";
-
+
+ std::cout << " ######### DD4hep::DDSegmentation::TiledLayerGridXY() " << std::endl ;
+
// register all necessary parameters
registerParameter("grid_size_x", "Cell size in X", _gridSizeX, 1., SegmentationParameter::LengthUnit);
registerParameter("grid_size_y", "Cell size in Y", _gridSizeY, 1., SegmentationParameter::LengthUnit);
diff --git a/DDSurfaces/include/DDSurfaces/ISurface.h b/DDSurfaces/include/DDSurfaces/ISurface.h
index 38d7bfdc5..49c693595 100644
--- a/DDSurfaces/include/DDSurfaces/ISurface.h
+++ b/DDSurfaces/include/DDSurfaces/ISurface.h
@@ -235,6 +235,9 @@ namespace DDSurfaces {
/** True if surface is parallel to Z with accuracy epsilon - result is cached in bit SurfaceType::ParallelToZ */
bool checkParallelToZ( const ISurface& surf , double epsilon=1.e-6 ) const {
+ if ( _bits[ SurfaceType::ParallelToZ ] ) // set in specific implementation
+ return true ;
+
double proj = std::fabs( surf.normal() * Vector3D(0.,0.,1.) ) ;
_bits.set( SurfaceType::ParallelToZ , ( proj < epsilon ) ) ;
@@ -249,6 +252,9 @@ namespace DDSurfaces {
/** True if surface is orthogonal to Z with accuracy epsilon - result is cached in bit SurfaceType::OrthogonalToZ */
bool checkOrthogonalToZ( const ISurface& surf , double epsilon=1.e-6 ) const {
+ if ( _bits[ SurfaceType::OrthogonalToZ ] ) // set in specific implementation
+ return true ;
+
double proj = std::fabs( surf.normal() * Vector3D(0.,0.,1.) ) ;
_bits.set( SurfaceType::OrthogonalToZ , ( std::fabs( proj - 1. ) < epsilon ) ) ;
--
GitLab