From 27938def09e2ffadf7dacd2d6050e59f145fbb5b Mon Sep 17 00:00:00 2001 From: Frank Gaede <frank.gaede@desy.de> Date: Tue, 22 Apr 2014 17:11:58 +0000 Subject: [PATCH] - added some cylinders to ILDExTPC - added drawing for cylinders in Surface.ccp/teve_display - added copy ctor and assignment to Vector3D -... --- DDRec/include/DDRec/Surface.h | 7 +- DDRec/src/Surface.cpp | 286 +++++++++++++++++++++-- DDSurfaces/include/DDSurfaces/ISurface.h | 10 +- DDSurfaces/include/DDSurfaces/Vector3D.h | 10 + UtilityApps/src/teve_display.cpp | 22 +- examples/ILDExDet/src/ILDExTPC_geo.cpp | 57 ++++- 6 files changed, 354 insertions(+), 38 deletions(-) diff --git a/DDRec/include/DDRec/Surface.h b/DDRec/include/DDRec/Surface.h index 605947e0b..0d6a30cd2 100644 --- a/DDRec/include/DDRec/Surface.h +++ b/DDRec/include/DDRec/Surface.h @@ -354,7 +354,12 @@ namespace DD4hep { /** 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 ) ; + // std::vector< Vector3D > getVertices( unsigned nMax=360 ) ; + + /** Get lines constraining the surface for drawing ( might not be exact boundaries) - + * at most nMax lines are returned. + */ + std::vector< std::pair< Vector3D, Vector3D> > getLines(unsigned nMax=100) ; protected: void initialize() ; diff --git a/DDRec/src/Surface.cpp b/DDRec/src/Surface.cpp index 1ce04b846..5d2d00c3b 100644 --- a/DDRec/src/Surface.cpp +++ b/DDRec/src/Surface.cpp @@ -150,12 +150,12 @@ namespace DD4hep { // unsigned count = volList.size() ; // for(unsigned i=0 ; i < count ; ++i) { - // std::cout << " **" ; + // std::cout << " **" ; // } - // std::cout << " searching for volume: " << std::hex << theVol.ptr() << " <-> pv.volume : " << pv.volume().ptr() - // << " pv.volume().ptr() == theVol.ptr() " << (pv.volume().ptr() == theVol.ptr() ) - // << std::endl ; - + // 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() ) { @@ -169,7 +169,9 @@ namespace DD4hep { if ( !node ) { - throw std::runtime_error("*** findVolume: Invalid placement: - node pointer Null ! " ); + // 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(); @@ -377,12 +379,12 @@ namespace DD4hep { } //=================================================================================================================== - std::vector< Vector3D > Surface::getVertices(unsigned nMax) { + std::vector< std::pair<Vector3D, Vector3D> > Surface::getLines(unsigned nMax) { const static double epsilon = 1e-6 ; - std::vector< Vector3D > _vert ; + std::vector< std::pair<Vector3D, Vector3D> > _vert ; // get local and global surface vectors @@ -438,12 +440,12 @@ namespace DD4hep { _wtM->LocalToMasterVect( ubl , ub.array() ) ; _wtM->LocalToMasterVect( vbl , vb.array() ) ; - _vert.resize(4) ; + _vert.reserve(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 ; + _vert.push_back( std::make_pair( _o + boxDim[ uidx ] * ub + boxDim[ vidx ] * vb , _o - boxDim[ uidx ] * ub + boxDim[ vidx ] * vb ) ) ; + _vert.push_back( std::make_pair( _o - boxDim[ uidx ] * ub + boxDim[ vidx ] * vb , _o - boxDim[ uidx ] * ub - boxDim[ vidx ] * vb ) ) ; + _vert.push_back( std::make_pair( _o - boxDim[ uidx ] * ub - boxDim[ vidx ] * vb , _o + boxDim[ uidx ] * ub - boxDim[ vidx ] * vb ) ) ; + _vert.push_back( std::make_pair( _o + boxDim[ uidx ] * ub - boxDim[ vidx ] * vb , _o + boxDim[ uidx ] * ub + boxDim[ vidx ] * vb ) ) ; return _vert ; } @@ -458,12 +460,15 @@ namespace DD4hep { // 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 ) ; + _vert.reserve( nMax ) ; double dAlpha = 2.* ROOT::Math::Pi() / double( nMax ) ; TVector3 normal( ln.x() , ln.y() , ln.z() ) ; + + DDSurfaces::Vector3D first, previous ; + for(unsigned i=0 ; i< nMax ; ++i ){ double alpha = double(i) * dAlpha ; @@ -489,19 +494,264 @@ namespace DD4hep { // std::cout << " **** normal:" << ln << " lu:" << lu << " alpha:" << alpha << " luRot:" << luRot << " lp :" << lp << " gp:" << gp << std::endl; - _vert[i] = gp ; + if( i > 0 ) + _vert.push_back( std::make_pair( previous, gp ) ) ; + else + first = gp ; + + previous = gp ; } + _vert.push_back( std::make_pair( previous, first ) ) ; - } else { // cylinder + } else if( type().isCylinder() ) { - // TODO - } + // if( shape->IsA() == TGeoTube::Class() ) { + if( shape->IsA() == TGeoConeSeg::Class() ) { + + _vert.reserve( nMax ) ; + + TGeoTube* tube = ( TGeoTube* ) shape ; + + double zHalf = tube->GetDZ() ; + + Vector3D zv( 0. , 0. , zHalf ) ; + + double r = lo.rho() ; + + + + unsigned n = nMax / 4 ; + 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.push_back( std::make_pair( pg0, pg1 ) ) ; + _vert.push_back( std::make_pair( pg1, pg2 ) ) ; + _vert.push_back( std::make_pair( pg2, pg3 ) ) ; + _vert.push_back( std::make_pair( pg3, pg0 ) ) ; + } + + // unsigned n = nMax / 4 ; + // 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.push_back( std::make_pair( pg0, pg1 ) ) ; + // _vert.push_back( std::make_pair( pg1, pg2 ) ) ; + // _vert.push_back( std::make_pair( pg2, pg3 ) ) ; + // _vert.push_back( std::make_pair( pg3, pg0 ) ) ; + // } + } + } return _vert ; } + // //=================================================================================================================== + + // 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 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 ; + + // } + //=================================================================================================================== diff --git a/DDSurfaces/include/DDSurfaces/ISurface.h b/DDSurfaces/include/DDSurfaces/ISurface.h index db030abbd..779e301c5 100644 --- a/DDSurfaces/include/DDSurfaces/ISurface.h +++ b/DDSurfaces/include/DDSurfaces/ISurface.h @@ -9,7 +9,7 @@ namespace DDSurfaces { - struct SurfaceType ; + class SurfaceType ; typedef long long int long64 ; @@ -78,12 +78,13 @@ namespace DDSurfaces { public: /// enum for defining the bits used to decode the properties enum{ - Cylinder = 1, + Cylinder = 0, Plane, Sensitive, Helper, ParallelToZ, - OrthogonalToZ + OrthogonalToZ, + Invisible } ; ///default c'tor @@ -136,6 +137,9 @@ namespace DDSurfaces { /// true if surface is orthogonal to Z bool isOrthogonalToZ() const { return _bits[ SurfaceType::OrthogonalToZ ] ; } + /// true if surface is not invisble - for drawing only + bool isVisible() const { return ! _bits[ SurfaceType::Invisible ] ; } + /// true if this is a cylinder parallel to Z bool isZCylinder() const { return ( _bits[ SurfaceType::Cylinder ] && _bits[ SurfaceType::ParallelToZ ] ) ; } diff --git a/DDSurfaces/include/DDSurfaces/Vector3D.h b/DDSurfaces/include/DDSurfaces/Vector3D.h index dae2cc3f4..97e63a754 100644 --- a/DDSurfaces/include/DDSurfaces/Vector3D.h +++ b/DDSurfaces/include/DDSurfaces/Vector3D.h @@ -25,6 +25,9 @@ namespace DDSurfaces { Vector3D() : _x(0.0),_y(0.0),_z(0.0) {} + /** Copy constructor*/ + Vector3D(const Vector3D& v) : _x(v[0]),_y(v[1]),_z(v[2]) {} + /** Constructor for float array.*/ Vector3D(const float* v) : _x(v[0]),_y(v[1]),_z(v[2]) {} @@ -54,6 +57,13 @@ namespace DDSurfaces { // _z( t.z() ){ // } + //assignment operator + Vector3D& operator=(const Vector3D& v) { + _x = v[0] ; + _y = v[1] ; + _z = v[2] ; + return *this ; + } /// fill vector from arbitrary class that defines operator[] template <class T> diff --git a/UtilityApps/src/teve_display.cpp b/UtilityApps/src/teve_display.cpp index 2c000d3e0..fc7d6c831 100644 --- a/UtilityApps/src/teve_display.cpp +++ b/UtilityApps/src/teve_display.cpp @@ -196,18 +196,21 @@ TEveStraightLineSet* getSurfaces() { Surface* surf = *it ; - const std::vector< Vector3D > vert = surf->getVertices() ; + if( ! surf->type().isVisible() ) + continue ; - if( vert.empty() ) - std::cout << " **** drawSurfaces() : empty vertices for surface " << *surf << std::endl ; + const std::vector< std::pair<Vector3D,Vector3D> > lines = surf->getLines() ; - unsigned nV = vert.size() ; + if( lines.empty() ) + std::cout << " **** drawSurfaces() : empty lines vector for surface " << *surf << std::endl ; - for( unsigned i=0 ; i<nV ; ++i){ + unsigned nL = lines.size() ; - unsigned j = ( i < nV-1 ? i+1 : 0 ) ; + for( unsigned i=0 ; i<nL ; ++i){ - ls->AddLine( vert[i].x(), vert[i].y(), vert[i].z(), vert[j].x() , vert[j].y() , vert[j].z() ) ; + + ls->AddLine( lines[i].first.x(), lines[i].first.y(), lines[i].first.z(), + lines[i].second.x(), lines[i].second.y(), lines[i].second.z() ) ; } ls->SetLineColor( kRed ) ; @@ -215,12 +218,11 @@ TEveStraightLineSet* getSurfaces() { ls->SetMarkerSize(.1); ls->SetMarkerStyle(4); - } // don't know what to do ... - - // gEve->AddElement(ls); + } return ls; } + //===================================================================================== void make_gui() { diff --git a/examples/ILDExDet/src/ILDExTPC_geo.cpp b/examples/ILDExDet/src/ILDExTPC_geo.cpp index 8ee073c38..d6aed607f 100644 --- a/examples/ILDExDet/src/ILDExTPC_geo.cpp +++ b/examples/ILDExDet/src/ILDExTPC_geo.cpp @@ -11,6 +11,7 @@ #include "DD4hep/DetFactoryHelper.h" #include "DD4hep/Detector.h" #include "DD4hep/TGeoUnits.h" +#include "DDRec/Surface.h" #include "TPCData.h" @@ -23,6 +24,8 @@ using namespace std; //using namespace tgeo ; using namespace DD4hep; using namespace DD4hep::Geometry; +using namespace DD4hep::DDRec ; +using namespace DDSurfaces ; static Ref_t create_element(LCDD& lcdd, xml_h e, SensitiveDetector sens) { xml_det_t x_det = e; @@ -76,23 +79,47 @@ static Ref_t create_element(LCDD& lcdd, xml_h e, SensitiveDetector sens) { Volume part_vol(part_nam,part_tub,part_mat); Position part_pos(px_pos.x(),px_pos.y(),px_pos.z()); - RotationZYX part_rot(px_rot.z(),px_rot.y(),px_rot.x()); - bool reflect = px_det.reflect(); + RotationZYX part_rot(px_rot.z(),px_rot.y(),px_rot.x()); + bool reflect = px_det.reflect(); sens.setType("tracker"); part_vol.setVisAttributes(lcdd,px_det.visStr()); + + // vectors for measurement plane + Vector3D u( 0. , 1. , 0. ) ; + Vector3D v( 0. , 0. , 1. ) ; + Vector3D n( 1. , 0. , 0. ) ; + + + double rcyl = ( px_tube.rmax() + px_tube.rmin() ) * 0.5 ; + double drcyl = px_tube.rmax() - px_tube.rmin() ; + Vector3D ocyl( rcyl , 0. , 0. ) ; + + //cache the important volumes in TPCData for later use without having to know their name switch(part_det.id()) { - case 2: + case 2: { tpcData->innerWall=part_det; + // add a surface to the det element + VolCylinder surf( part_vol , SurfaceType( SurfaceType::Helper ) , drcyl*.5 , drcyl*.5 , u,v,n , ocyl ) ; + volSurfaceList( part_det )->push_back( surf ) ; + } break; - case 3: + + case 3: { tpcData->outerWall=part_det; + // add a surface to the det element + VolCylinder surf( part_vol , SurfaceType( SurfaceType::Helper ) , drcyl*.5 , drcyl*.5 , u,v,n , ocyl ) ; + volSurfaceList( part_det )->push_back( surf ) ; + } break; + case 4: { + tpcData->gasVolume=part_det; + xml_comp_t px_lay(px_det.child(_U(layer))); int nTPClayer( px_lay.attr<int>(_U(number)) ); @@ -112,6 +139,7 @@ static Ref_t create_element(LCDD& lcdd, xml_h e, SensitiveDetector sens) { g_phiMax = 6.283185307e+00 ; // FIXME: where to define ? is it allways 2PI ? //------------------------------------------------------ + for(int i=0 ; i < nTPClayer ; ++i){ Tube gas_tubL( r0 + (2*i) * dR , r0 + (2*i+1) * dR , zh ); @@ -122,11 +150,21 @@ static Ref_t create_element(LCDD& lcdd, xml_h e, SensitiveDetector sens) { Volume gas_volU( _toString( i, "tpc_row_upper_%03d") , gas_tubU, part_mat); gas_volU.setSensitiveDetector( sens ); - part_vol.placeVolume( gas_volU, RotationZYX(0,0,0) ).addPhysVolID("layer",i) ; + + PlacedVolume pv = part_vol.placeVolume( gas_volU, RotationZYX(0,0,0) ) ; + pv.addPhysVolID("layer",i) ; + + DetElement layerDE( tpc , _toString( i, "tpc_row_upper_%03d") , x_det.id() ); + layerDE.setPlacement( pv ) ; + + Vector3D o( r0 + (2*i+1) * dR , 0. , 0. ) ; + + VolCylinder surf( gas_volU , SurfaceType(SurfaceType::Sensitive, SurfaceType::Invisible ) , dR , dR , u,v,n ,o ) ; + + volSurfaceList( layerDE )->push_back( surf ) ; } - tpcData->gasVolume=part_det; break; } case 5: @@ -136,6 +174,12 @@ static Ref_t create_element(LCDD& lcdd, xml_h e, SensitiveDetector sens) { //Endplate if(part_det.id()== 0){ tpcData->endplate=part_det; + + // add a plane to the endcap volume + // note: u and v are exchanged: normal is along z ... + VolPlane surf( part_vol , SurfaceType( SurfaceType::Helper ) , px_tube.zhalf() , x_tube.zhalf(), u , n , v ) ; + volSurfaceList( part_det )->push_back( surf ) ; + //modules int mdcount=0; for(xml_coll_t m(px_det,_U(modules)); m; ++m) { @@ -180,6 +224,7 @@ static Ref_t create_element(LCDD& lcdd, xml_h e, SensitiveDetector sens) { part_phv.addPhysVolID("side",0); part_det.setPlacement(part_phv); tpc.add(part_det); + //now reflect it if(reflect){ Position r_pos(px_pos.x(),px_pos.y(),-px_pos.z()); -- GitLab