Newer
Older
#include "DD4hep/objects/DetectorInterna.h"
#include "DDRec/MaterialManager.h"
#include <exception>
#include "TGeoMatrix.h"
#include "TGeoShape.h"
#include "TRotation.h"
//TGeoTrd1 is apparently not included by defautl
#include "TGeoTrd1.h"
namespace DD4hep {
namespace DDRec {
using namespace Geometry ;
//======================================================================================================
void VolSurfaceBase::setU(const Vector3D& u_val) { _u = u_val ; }
void VolSurfaceBase::setV(const Vector3D& v_val) { _v = v_val ; }
void VolSurfaceBase::setNormal(const Vector3D& n) { _n = n ; }
long64 VolSurfaceBase::id() const { return 0 ; }
const SurfaceType& VolSurfaceBase::type() const { return _type ; }
Vector3D VolSurfaceBase::u( const Vector3D& point ) const { point.x() ; return _u ; }
Vector3D VolSurfaceBase::v(const Vector3D& point ) const { point.x() ; return _v ; }
Vector3D VolSurfaceBase::normal(const Vector3D& point ) const { point.x() ; return _n ; }
const Vector3D& VolSurfaceBase::origin() const { return _o ;}
Vector2D VolSurfaceBase::globalToLocal( const Vector3D& point) const {
Vector3D p = point - origin() ;
// create new orthogonal unit vectors
// FIXME: these vectors should be cached really ...
double uv = u() * v() ;
Vector3D uprime = ( u() - uv * v() ).unit() ;
Vector3D vprime = ( v() - uv * u() ).unit() ;
double uup = u() * uprime ;
double vvp = v() * vprime ;
return Vector2D( p*uprime / uup , p*vprime / vvp ) ;
Vector3D VolSurfaceBase::localToGlobal( const Vector2D& point) const {
Vector3D g = origin() + point[0] * u() + point[1] * v() ;
return g ;
}
const IMaterial& VolSurfaceBase::innerMaterial() const{ return _innerMat ; }
const IMaterial& VolSurfaceBase::outerMaterial() const { return _outerMat ; }
double VolSurfaceBase::innerThickness() const { return _th_i ; }
double VolSurfaceBase::outerThickness() const { return _th_o ; }
double VolSurfaceBase::length_along_u() const {
const DDSurfaces::Vector3D& u_val = this->u( o ) ;
DDSurfaces::Vector3D um = -1. * u_val ;
double dist_p = 0. ;
double dist_m = 0. ;
// std::cout << " VolSurfaceBase::length_along_u() : o = " << o << " u = " << this->u( o )
// << " -u = " << um << std::endl ;
if( volume()->GetShape()->Contains( o.const_array() ) ){
dist_p = volume()->GetShape()->DistFromInside( const_cast<double*> ( o.const_array() ) ,
const_cast<double*> ( u_val.const_array() ) ) ;
dist_m = volume()->GetShape()->DistFromInside( const_cast<double*> ( o.const_array() ) ,
const_cast<double*> ( um.array() ) ) ;
// std::cout << " VolSurfaceBase::length_along_u() : shape contains(o) = " << volume()->GetShape()->Contains( o.const_array() )
// << " dist_p " << dist_p
// << " dist_m " << dist_m
// << std::endl ;
} else{
dist_p = volume()->GetShape()->DistFromOutside( const_cast<double*> ( o.const_array() ) ,
const_cast<double*> ( u_val.const_array() ) ) ;
dist_m = volume()->GetShape()->DistFromOutside( const_cast<double*> ( o.const_array() ) ,
const_cast<double*> ( um.array() ) ) ;
dist_p *= 1.0001 ;
dist_m *= 1.0001 ;
// std::cout << " VolSurfaceBase::length_along_u() : shape contains(o) = " << volume()->GetShape()->Contains( o.const_array() )
// << " dist_p " << dist_p
// << " dist_m " << dist_m
// << std::endl ;
DDSurfaces::Vector3D o_1 = this->origin() + dist_p * u_val ;
DDSurfaces::Vector3D o_2 = this->origin() + dist_m * um ;
dist_p += volume()->GetShape()->DistFromInside( const_cast<double*> ( o_1.const_array() ) ,
const_cast<double*> ( u_val.const_array() ) ) ;
dist_m += volume()->GetShape()->DistFromInside( const_cast<double*> ( o_2.const_array() ) ,
const_cast<double*> ( um.array() ) ) ;
// std::cout << " VolSurfaceBase::length_along_u() : shape contains(o) = " << volume()->GetShape()->Contains( o.const_array() )
// << " dist_p " << dist_p
// << " dist_m " << dist_m
// << std::endl ;
}
double VolSurfaceBase::length_along_v() const {
const DDSurfaces::Vector3D& v_val = this->v( o ) ;
DDSurfaces::Vector3D vm = -1. * v_val ;
double dist_p = 0. ;
double dist_m = 0. ;
// std::cout << " VolSurfaceBase::length_along_u() : o = " << o << " u = " << this->u( o )
// << " -u = " << vm << std::endl ;
if( volume()->GetShape()->Contains( o.const_array() ) ){
dist_p = volume()->GetShape()->DistFromInside( const_cast<double*> ( o.const_array() ) ,
const_cast<double*> ( v_val.const_array() ) ) ;
dist_m = volume()->GetShape()->DistFromInside( const_cast<double*> ( o.const_array() ) ,
const_cast<double*> ( vm.array() ) ) ;
// std::cout << " VolSurfaceBase::length_along_u() : shape contains(o) = " << volume()->GetShape()->Contains( o.const_array() )
// << " dist_p " << dist_p
// << " dist_m " << dist_m
// << std::endl ;
} else{
dist_p = volume()->GetShape()->DistFromOutside( const_cast<double*> ( o.const_array() ) ,
const_cast<double*> ( v_val.const_array() ) ) ;
dist_m = volume()->GetShape()->DistFromOutside( const_cast<double*> ( o.const_array() ) ,
const_cast<double*> ( vm.array() ) ) ;
dist_p *= 1.0001 ;
dist_m *= 1.0001 ;
// std::cout << " VolSurfaceBase::length_along_u() : shape contains(o) = " << volume()->GetShape()->Contains( o.const_array() )
// << " dist_p " << dist_p
// << " dist_m " << dist_m
// << std::endl ;
DDSurfaces::Vector3D o_1 = this->origin() + dist_p * v_val ;
DDSurfaces::Vector3D o_2 = this->origin() + dist_m * vm ;
dist_p += volume()->GetShape()->DistFromInside( const_cast<double*> ( o_1.const_array() ) ,
const_cast<double*> ( v_val.const_array() ) ) ;
dist_m += volume()->GetShape()->DistFromInside( const_cast<double*> ( o_2.const_array() ) ,
const_cast<double*> ( vm.array() ) ) ;
// std::cout << " VolSurfaceBase::length_along_u() : shape contains(o) = " << volume()->GetShape()->Contains( o.const_array() )
// << " dist_p " << dist_p
// << " dist_m " << dist_m
// << std::endl ;
}
double VolSurfaceBase::distance(const Vector3D& point ) const { point.x() ; return 1.e99 ; }
/// Checks if the given point lies within the surface
bool VolSurfaceBase::insideBounds(const Vector3D& point, double epsilon) const {
double dist = std::abs ( distance( point ) ) ;
bool inShape = volume()->GetShape()->Contains( point.const_array() ) ;
std::cout << " ** Surface::insideBound( " << point << " ) - distance = " << dist
<< " origin = " << origin() << " normal = " << normal()
<< " p * n = " << point * normal()
<< " isInShape : " << inShape << std::endl ;
//fixme: older versions of ROOT (~<5.34.10 ) take a non const pointer as argument - therefore use a const cast here for the time being ...
return ( std::abs ( distance( point ) ) < epsilon ) && volume()->GetShape()->Contains( const_cast<double*> (point.const_array() ) ) ;
std::vector< std::pair<Vector3D, Vector3D> > VolSurfaceBase::getLines(unsigned nMax) {
// dummy implementation returning empty set
nMax ;
std::vector< std::pair<Vector3D, Vector3D> > lines ;
return lines ;
}
//===================================================================
// simple wrapper methods forwarding the call to the implementation object
long64 VolSurface::id() const { return _surf->id() ; }
const SurfaceType& VolSurface::type() const { return _surf->type() ; }
Vector3D VolSurface::u( const Vector3D& point ) const { return _surf->u(point) ; }
Vector3D VolSurface::v(const Vector3D& point ) const { return _surf->v(point) ; }
Vector3D VolSurface::normal(const Vector3D& point ) const { return _surf->normal(point) ; }
const Vector3D& VolSurface::origin() const { return _surf->origin() ;}
Vector2D VolSurface::globalToLocal( const Vector3D& point) const { return _surf->globalToLocal( point ) ; }
Vector3D VolSurface::localToGlobal( const Vector2D& point) const { return _surf->localToGlobal( point) ; }
const IMaterial& VolSurface::innerMaterial() const{ return _surf->innerMaterial() ; }
const IMaterial& VolSurface::outerMaterial() const { return _surf->outerMaterial() ; }
double VolSurface::innerThickness() const { return _surf->innerThickness() ; }
double VolSurface::outerThickness() const { return _surf->outerThickness() ; }
double VolSurface::length_along_u() const { return _surf->length_along_u() ; }
double VolSurface::length_along_v() const { return _surf->length_along_v() ; }
double VolSurface::distance(const Vector3D& point ) const { return _surf->distance( point ) ; }
bool VolSurface::insideBounds(const Vector3D& point, double epsilon) const {
return _surf->insideBounds( point, epsilon ) ;
}
std::vector< std::pair<Vector3D, Vector3D> > VolSurface::getLines(unsigned nMax) {
return _surf->getLines(nMax) ;
}
//===================================================================
/** Distance to planar surface */
double VolPlaneImpl::distance(const Vector3D& point ) const {
return ( point - origin() ) * normal() ;
}
//======================================================================================================
VolCylinderImpl::VolCylinderImpl( Geometry::Volume vol, SurfaceType typ,
double thickness_inner ,double thickness_outer, Vector3D o ) :
VolSurfaceBase(typ, thickness_inner, thickness_outer, Vector3D() , Vector3D() , Vector3D() , o , vol, 0) {
Vector3D v_val( 0., 0., 1. ) ;
Vector3D o_rphi( o.x() , o.y() , 0. ) ;
_type.setProperty( SurfaceType::Plane , false ) ;
_type.setProperty( SurfaceType::Cylinder , true ) ;
_type.checkParallelToZ( *this ) ;
_type.checkOrthogonalToZ( *this ) ;
Vector3D VolCylinderImpl::u(const Vector3D& point ) const {
Vector3D n( 1. , point.phi() , 0. , Vector3D::cylindrical ) ;
Vector3D VolCylinderImpl::normal(const Vector3D& point ) const {
// normal is just given by phi of the point
return Vector3D( 1. , point.phi() , 0. , Vector3D::cylindrical ) ;
Vector2D VolCylinderImpl::globalToLocal( const Vector3D& point) const {
// cylinder is parallel to v here so u is Z and v is r *Phi
double phi = point.phi() - origin().phi() ;
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() ) ;
Vector3D VolCylinderImpl::localToGlobal( const Vector2D& point) const {
double z = point.v() + origin().z() ;
double phi = point.u() / origin().rho() + 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 ) ;
}
double VolCylinderImpl::distance(const Vector3D& point ) const {
//================================================================================================================
// // delete all surfaces attached to this volume
// for( VolSurfaceList::iterator i=begin(), n=end() ; i !=n ; ++i ) {
// i->clear() ;
// }
}
//=======================================================
SurfaceList::~SurfaceList(){
if( _isOwner ) {
// delete all surfaces attached to this volume
for( SurfaceList::iterator i=begin(), n=end() ; i !=n ; ++i ) {
delete (*i) ;
}
}
}
//================================================================================================================
VolSurfaceList* volSurfaceList( DetElement& det ) {
VolSurfaceList* list = 0 ;
try {
} catch(const std::runtime_error& e){
list = det.addExtension<VolSurfaceList >( new VolSurfaceList ) ;
//======================================================================================================================
bool findVolume( PlacedVolume pv, Volume theVol, std::list< PlacedVolume >& volList ) {
volList.push_back( pv ) ;
// unsigned count = volList.size() ;
// for(unsigned i=0 ; i < count ; ++i) {
// std::cout << " **" ;
// }
// std::cout << " searching for volume: " << theVol.name() << " " << std::hex << theVol.ptr() << " <-> pv.volume : " << pv.name() << " " << pv.volume().ptr()
// << " pv.volume().ptr() == theVol.ptr() " << (pv.volume().ptr() == theVol.ptr() )
// << std::endl ;
if( pv.volume().ptr() == theVol.ptr() ) {
} else {
// std::cout << " *** findVolume: Invalid placement: - node pointer Null for volume: " << pv.name() << std::endl ;
throw std::runtime_error("*** findVolume: Invalid placement: - node pointer Null ! " + std::string( pv.name() ) );
}
// Volume vol = pv.volume();
// std::cout << " ndau = " << node->GetNdaughters() << std::endl ;
for (Int_t idau = 0, ndau = node->GetNdaughters(); idau < ndau; ++idau) {
TGeoNode* daughter = node->GetDaughter(idau);
PlacedVolume placement( daughter );
if ( !placement.data() ) {
throw std::runtime_error("*** findVolume: Invalid not instrumented placement:"+std::string(daughter->GetName())
+" [Internal error -- bad detector constructor]");
}
PlacedVolume pv_dau = Ref_t(daughter); // why use a Ref_t here ???
// std::cout << " ----- found in daughter volume !!! " << std::hex << pv_dau.volume().ptr() << std::endl ;
return false ;
//--------------------------------
}
}
//======================================================================================================================
Surface::Surface( Geometry::DetElement det, VolSurface volSurf ) : _det( det) , _volSurf( volSurf ),
initialize() ;
}
long64 Surface::id() const { return _id ; }
const SurfaceType& Surface::type() const { return _type ; }
Vector3D Surface::u( const Vector3D& point ) const { point.x() ; return _u ; }
Vector3D Surface::v(const Vector3D& point ) const { point.x() ; return _v ; }
Vector3D Surface::normal(const Vector3D& point ) const { point.x() ; return _n ; }
const Vector3D& Surface::origin() const { return _o ;}
double Surface::innerThickness() const { return _volSurf.innerThickness() ; }
double Surface::outerThickness() const { return _volSurf.outerThickness() ; }
double Surface::length_along_u() const { return _volSurf.length_along_u() ; }
double Surface::length_along_v() const { return _volSurf.length_along_v() ; }
const IMaterial& Surface::innerMaterial() const {
const IMaterial& mat = _volSurf.innerMaterial() ;
if( ! ( mat.Z() > 0 ) ) {
Vector3D p = _o - innerThickness() * _n ;
const MaterialVec& materials = matMgr.materialsBetween( _o , p ) ;
_volSurf.setInnerMaterial( materials.size() > 1 ?
matMgr.createAveragedMaterial( materials ) :
materials[0].first ) ;
return mat ;
const IMaterial& Surface::outerMaterial() const {
const IMaterial& mat = _volSurf.outerMaterial() ;
if( ! ( mat.Z() > 0 ) ) {
Vector3D p = _o - outerThickness() * _n ;
const MaterialVec& materials = matMgr.materialsBetween( _o , p ) ;
_volSurf.setOuterMaterial( materials.size() > 1 ?
matMgr.createAveragedMaterial( materials ) :
materials[0].first ) ;
return mat ;
Vector2D Surface::globalToLocal( const Vector3D& point) const {
// create new orthogonal unit vectors
// FIXME: these vectors should be cached really ...
double uv = u() * v() ;
Vector3D uprime = ( u() - uv * v() ).unit() ;
Vector3D vprime = ( v() - uv * u() ).unit() ;
double uup = u() * uprime ;
double vvp = v() * vprime ;
return Vector2D( p*uprime / uup , p*vprime / vvp ) ;
Vector3D Surface::localToGlobal( const Vector2D& point) const {
Vector3D g = origin() + point[0] * u() + point[1] * v() ;
return g ;
}
Vector3D Surface::volumeOrigin() const {
double o_array[3] ;
_wtM->LocalToMaster ( Vector3D() , o_array ) ;
Vector3D o(o_array) ;
return o ;
}
double Surface::distance(const Vector3D& point ) const {
double pa[3] ;
_wtM->MasterToLocal( point , pa ) ;
Vector3D localPoint( pa ) ;
return _volSurf.distance( localPoint ) ;
//FG return ( _volSurf.type().isPlane() ? VolPlane(_volSurf).distance( localPoint ) : VolCylinder(_volSurf).distance( localPoint ) ) ;
}
bool Surface::insideBounds(const Vector3D& point, double epsilon) const {
double pa[3] ;
_wtM->MasterToLocal( point , pa ) ;
Vector3D localPoint( pa ) ;
return _volSurf.insideBounds( localPoint , epsilon) ;
//FG return ( _volSurf.type().isPlane() ? VolPlane(_volSurf).insideBounds( localPoint, epsilon ) : VolCylinder(_volSurf).insideBounds( localPoint , epsilon) ) ;
void Surface::initialize() {
// first we need to find the right volume for the local surface in the DetElement's volumes
std::list< PlacedVolume > pVList ;
PlacedVolume pv = _det.placement() ;
Volume theVol = _volSurf.volume() ;
if( ! findVolume( pv, theVol , pVList ) ){
theVol = _volSurf.volume() ;
std::stringstream sst ; sst << " ***** ERROR: Volume " << theVol.name() << " not found for DetElement " << _det.name() << " with surface " ;
throw std::runtime_error( sst.str() ) ;
// std::cout << " **** Surface::initialize() # placements for surface = " << pVList.size()
// << " worldTransform : "
// << std::endl ;
//=========== compute and cache world transform for surface ==========
const TGeoHMatrix& wm = _det.worldTransformation() ;
for( std::list<PlacedVolume>::iterator it= pVList.begin(), n = pVList.end() ; it != n ; ++it ){
PlacedVolume pv = *it ;
TGeoMatrix* m = pv->GetMatrix();
std::cout << " +++ matrix for placed volume : " << std::endl ;
m->Print() ;
}
#endif
// need to get the inverse transformation ( see Detector.cpp )
// std::auto_ptr<TGeoHMatrix> wtI( new TGeoHMatrix( wm.Inverse() ) ) ;
// has been fixed now, no need to get the inverse anymore:
dd4hep_ptr<TGeoHMatrix> wtI( new TGeoHMatrix( wm ) ) ;
//---- if the volSurface is not in the DetElement's volume, we need to mutliply the path to the volume to the
// DetElements world transform
for( std::list<PlacedVolume>::iterator it = ++( pVList.begin() ) , n = pVList.end() ; it != n ; ++it ){
PlacedVolume pvol = *it ;
TGeoMatrix* m = pvol->GetMatrix();
// std::cout << " +++ matrix for placed volume : " << std::endl ;
// m->Print() ;
// std::cout << " +++ new world transform matrix : " << std::endl ;
#if 0 //fixme: which convention to use here - the correct should be wtI, however it is the inverse of what is stored in DetElement ???
dd4hep_ptr<TGeoHMatrix> wt( new TGeoHMatrix( wtI->Inverse() ) ) ;
// cache the world transform for the surface
_wtM = wt.release() ;
#else
// wtI->Print() ;
// cache the world transform for the surface
_wtM = wtI.release() ;
#endif
// ============ now fill the global surface vectors ==========================
double ua[3], va[3], na[3], oa[3] ;
_wtM->LocalToMasterVect( _volSurf.u() , ua ) ;
_wtM->LocalToMasterVect( _volSurf.v() , va ) ;
_wtM->LocalToMasterVect( _volSurf.normal() , na ) ;
_wtM->LocalToMaster ( _volSurf.origin() , oa ) ;
_u.fill( ua ) ;
_v.fill( va ) ;
_n.fill( na ) ;
_o.fill( oa ) ;
// std::cout << " --- local and global surface vectors : ------- " << std::endl
// << " u : " << _volSurf.u() << " - " << _u << std::endl
// << " v : " << _volSurf.v() << " - " << _v << std::endl
// << " n : " << _volSurf.normal() << " - " << _n << std::endl
// << " o : " << _volSurf.origin() << " - " << _o << std::endl ;
// =========== check parallel and orthogonal to Z ===================
_type.checkOrthogonalToZ( *this ) ;
//======== set the unique surface ID from the DetElement ( and placements below ? )
// just use the DetElement ID for now ...
_id = _det.volumeID() ;
// typedef PlacedVolume::VolIDs IDV ;
// DetElement d = _det ;
// while( d.isValid() && d.parent().isValid() ){
// PlacedVolume pv = d.placement() ;
// if( pv.isValid() ){
// const IDV& idV = pv.volIDs() ;
// std::cout << " VolIDs : " << d.name() << std::endl ;
// for( unsigned i=0, n=idV.size() ; i<n ; ++i){
// std::cout << " " << idV[i].first << " - " << idV[i].second << std::endl ;
// }
// }
// d = d.parent() ;
// }
}
//===================================================================================================================
std::vector< std::pair<Vector3D, Vector3D> > Surface::getLines(unsigned nMax) {
const static double epsilon = 1e-6 ;
std::vector< std::pair<Vector3D, Vector3D> > lines ;
//--------------------------------------------
// check if there are lines defined in the VolSurface :
const std::vector< std::pair<Vector3D, Vector3D> >& local_lines = _volSurf.getLines() ;
if( local_lines.size() > 0 ) {
unsigned n=local_lines.size() ;
lines.reserve( n ) ;
for( unsigned i=0;i<n;++i){
DDSurfaces::Vector3D av,bv;
_wtM->LocalToMasterVect( local_lines[i].first , av.array() ) ;
_wtM->LocalToMasterVect( local_lines[i].second , bv.array() ) ;
lines.push_back( std::make_pair( av, bv ) );
}
return lines ;
}
//--------------------------------------------
// get local and global surface vectors
const DDSurfaces::Vector3D& lu = _volSurf.u() ;
const DDSurfaces::Vector3D& ln = _volSurf.normal() ;
const DDSurfaces::Vector3D& lo = _volSurf.origin() ;
Volume vol = volume() ;
const TGeoShape* shape = vol->GetShape() ;
if( type().isPlane() ) {
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. ) ;
ubl.fill( 1., 0., 0. ) ;
vbl.fill( 0., 0., 1. ) ;
uidx = 0 ;
vidx = 2 ;
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() ) ;
lines.push_back( std::make_pair( _o + boxDim[ uidx ] * ub + boxDim[ vidx ] * vb , _o - boxDim[ uidx ] * ub + boxDim[ vidx ] * vb ) ) ;
lines.push_back( std::make_pair( _o - boxDim[ uidx ] * ub + boxDim[ vidx ] * vb , _o - boxDim[ uidx ] * ub - boxDim[ vidx ] * vb ) ) ;
lines.push_back( std::make_pair( _o - boxDim[ uidx ] * ub - boxDim[ vidx ] * vb , _o + boxDim[ uidx ] * ub - boxDim[ vidx ] * vb ) ) ;
lines.push_back( std::make_pair( _o + boxDim[ uidx ] * ub - boxDim[ vidx ] * vb , _o + boxDim[ uidx ] * ub + boxDim[ vidx ] * vb ) ) ;
} else if( shape->IsA() == TGeoConeSeg::Class() ) {
// can only deal with special case of z-disk and origin in center of cone
if( type().isZDisk() && lo.rho() < epsilon ) {
double zhalf = cone->GetDZ() ;
double rmax1 = cone->GetRmax1() ;
double rmax2 = cone->GetRmax2() ;
double rmin1 = cone->GetRmin1() ;
double rmin2 = cone->GetRmin2() ;
// two circles around origin
// get radii at position of plane
double r0 = rmin1 + ( rmin2 - rmin1 ) / ( 2. * zhalf ) * ( zhalf + lo.z() ) ;
double r1 = rmax1 + ( rmax2 - rmax1 ) / ( 2. * zhalf ) * ( zhalf + lo.z() ) ;
unsigned n = nMax / 4 ;
double dPhi = 2.* ROOT::Math::Pi() / double( n ) ;
Vector3D rv00( r0*sin( i *dPhi ) , r0*cos( i *dPhi ) , 0. ) ;
Vector3D rv01( r0*sin( (i+1)*dPhi ) , r0*cos( (i+1)*dPhi ) , 0. ) ;
Vector3D rv10( r1*sin( i *dPhi ) , r1*cos( i *dPhi ) , 0. ) ;
Vector3D rv11( r1*sin( (i+1)*dPhi ) , r1*cos( (i+1)*dPhi ) , 0. ) ;
Vector3D pl0 = lo + rv00 ;
Vector3D pl1 = lo + rv01 ;
Vector3D pl2 = lo + rv10 ;
Vector3D pl3 = lo + rv11 ;
_wtM->LocalToMaster( pl0, pg0.array() ) ;
_wtM->LocalToMaster( pl1, pg1.array() ) ;
_wtM->LocalToMaster( pl2, pg2.array() ) ;
_wtM->LocalToMaster( pl3, pg3.array() ) ;
lines.push_back( std::make_pair( pg0, pg1 ) ) ;
lines.push_back( std::make_pair( pg2, pg3 ) ) ;
}
//add some vertical and horizontal lines so that the disc is seen in the rho-z projection
n = 4 ; dPhi = 2.* ROOT::Math::Pi() / double( n ) ;
Vector3D rv0( r0*sin( i * dPhi ) , r0*cos( i * dPhi ) , 0. ) ;
Vector3D rv1( r1*sin( i * dPhi ) , r1*cos( i * dPhi ) , 0. ) ;
Vector3D pl0 = lo + rv0 ;
Vector3D pl1 = lo + rv1 ;
_wtM->LocalToMaster( pl0, pg0.array() ) ;
_wtM->LocalToMaster( pl1, pg1.array() ) ;
lines.push_back( std::make_pair( pg0, pg1 ) ) ;
}
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
else if(shape->IsA() == TGeoTrap::Class()) {
TGeoTrap* trapezoid = ( TGeoTrap* ) shape;
double dx1 = trapezoid->GetBl1();
double dx2 = trapezoid->GetTl1();
double dz = trapezoid->GetH1();
//according to the TGeoTrap definition, the lengths are given such that the normal vector of the surface
//points in the e_z direction.
DDSurfaces::Vector3D ubl( 1., 0., 0. ) ;
DDSurfaces::Vector3D vbl( 0., 1., 0. ) ;
//the local span vectors are transformed into the main coordinate system (in LocalToMasterVect())
DDSurfaces::Vector3D ub ;
DDSurfaces::Vector3D vb ;
_wtM->LocalToMasterVect( ubl , ub.array() ) ;
_wtM->LocalToMasterVect( vbl , vb.array() ) ;
//the trapezoid is drawn as a set of four lines connecting its four corners
lines.reserve(4) ;
//_o is vector to the origin
lines.push_back( std::make_pair( _o + dx1 * ub - dz * vb , _o + dx2 * ub + dz * vb ) ) ;
lines.push_back( std::make_pair( _o + dx2 * ub + dz * vb , _o - dx2 * ub + dz * vb ) ) ;
lines.push_back( std::make_pair( _o - dx2 * ub + dz * vb , _o - dx1 * ub - dz * vb) ) ;
lines.push_back( std::make_pair( _o - dx1 * ub - dz * vb , _o + dx1 * ub - dz * vb) ) ;
return lines;
}
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
//added code by Thorben Quast for simplified set of lines for trapezoids with unequal lengths in x
else if(shape->IsA() == TGeoTrd1::Class()){
TGeoTrd1* trapezoid = ( TGeoTrd1* ) shape;
//all lengths are half length
double dx1 = trapezoid->GetDx1();
double dx2 = trapezoid->GetDx2();
//double dy = trapezoid->GetDy();
double dz = trapezoid->GetDz();
//the normal vector is parallel to e_y for all geometry cases in CLIC
//if that is at some point not the case anymore, then local plane vectors ubl, vbl
//must be initialized like it is done for the boxes (line 674 following)
DDSurfaces::Vector3D ubl( 1., 0., 0. ) ;
DDSurfaces::Vector3D vbl( 0., 0., 1. ) ;
//the local span vectors are transformed into the main coordinate system (in LocalToMasterVect())
DDSurfaces::Vector3D ub ;
DDSurfaces::Vector3D vb ;
_wtM->LocalToMasterVect( ubl , ub.array() ) ;
_wtM->LocalToMasterVect( vbl , vb.array() ) ;
//the trapezoid is drawn as a set of four lines connecting its four corners
lines.reserve(4) ;
//_o is vector to the origin
lines.push_back( std::make_pair( _o + dx1 * ub - dz * vb , _o + dx2 * ub + dz * vb ) ) ;
lines.push_back( std::make_pair( _o + dx2 * ub + dz * vb , _o - dx2 * ub + dz * vb ) ) ;
lines.push_back( std::make_pair( _o - dx2 * ub + dz * vb , _o - dx1 * ub - dz * vb) ) ;
lines.push_back( std::make_pair( _o - dx1 * ub - dz * vb , _o + dx1 * ub - dz * vb) ) ;
return lines;
}
//added code by Thorben Quast for simplified set of lines for trapezoids with unequal lengths in x AND y
else if(shape->IsA() == TGeoTrd2::Class()){
TGeoTrd2* trapezoid = ( TGeoTrd2* ) shape;
//all lengths are half length
double dx1 = trapezoid->GetDx1();
double dx2 = trapezoid->GetDx2();
//double dy1 = trapezoid->GetDy1();
//double dy2 = trapezoid->GetDy2();
double dz = trapezoid->GetDz();
//the normal vector is parallel to e_y for all geometry cases in CLIC
//if that is at some point not the case anymore, then local plane vectors ubl, vbl
//must be initialized like it is done for the boxes (line 674 following)
DDSurfaces::Vector3D ubl( 1., 0., 0. ) ;
DDSurfaces::Vector3D vbl( 0., 0., 1. ) ;
//the local span vectors are transformed into the main coordinate system (in LocalToMasterVect())
DDSurfaces::Vector3D ub ;
DDSurfaces::Vector3D vb ;
_wtM->LocalToMasterVect( ubl , ub.array() ) ;
_wtM->LocalToMasterVect( vbl , vb.array() ) ;
//the trapezoid is drawn as a set of four lines connecting its four corners
lines.reserve(4) ;
//_o is vector to the origin
lines.push_back( std::make_pair( _o + dx1 * ub - dz * vb , _o + dx2 * ub + dz * vb ) ) ;
lines.push_back( std::make_pair( _o + dx2 * ub + dz * vb , _o - dx2 * ub + dz * vb ) ) ;
lines.push_back( std::make_pair( _o - dx2 * ub + dz * vb , _o - dx1 * ub - dz * vb) ) ;
lines.push_back( std::make_pair( _o - dx1 * ub - dz * vb , _o + dx1 * ub - dz * vb) ) ;
return lines;
}
// ===== 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...
double dAlpha = 2.* ROOT::Math::Pi() / double( nMax ) ;
DDSurfaces::Vector3D luRot ;
luRot.fill( vecR ) ;
double dist = shape->DistFromInside( const_cast<double*> (lo.const_array()) , const_cast<double*> (luRot.const_array()) , 3, 0.1 ) ;
// local point at volume boundary
DDSurfaces::Vector3D lp = lo + dist * luRot ;
_wtM->LocalToMaster( lp , gp.array() ) ;
// std::cout << " **** normal:" << ln << " lu:" << lu << " alpha:" << alpha << " luRot:" << luRot << " lp :" << lp << " gp:" << gp << " dist : " << dist
// << " is point " << gp << " inside : " << shape->Contains( gp.const_array() )
// << " dist from outside for lo,lu " << shape->DistFromOutside( lo.const_array() , lu.const_array() , 3 )
// << " dist from inside for lo,ln " << shape->DistFromInside( lo.const_array() , ln.const_array() , 3 )
// << std::endl;
// shape->Dump() ;
if( i > 0 )
lines.push_back( std::make_pair( previous, gp ) ) ;
else
first = gp ;
previous = gp ;
}
lines.push_back( std::make_pair( previous, first ) ) ;
// if( shape->IsA() == TGeoTube::Class() ) {
if( shape->IsA() == TGeoConeSeg::Class() ) {
unsigned n = nMax / 4 ;
double dPhi = 2.* ROOT::Math::Pi() / double( n ) ;
Vector3D rv0( r*sin( i *dPhi ) , r*cos( i *dPhi ) , 0. ) ;
Vector3D rv1( r*sin( (i+1)*dPhi ) , r*cos( (i+1)*dPhi ) , 0. ) ;