Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-04-10 07:50:33 

0001 #pragma once
0002 /**
0003 U4Solid.h : Convert G4VSolid CSG trees into sn.h trees
0004 ==========================================================
0005 
0006 Canonical usage from U4Tree::initSolid
0007 
0008 
0009 Priors:
0010 
0011 x4/X4Solid
0012    G4VSolid->NNode
0013 
0014 u4/U4SolidTree
0015    some overlap with X4Entity, tree mechanics might be better that X4Solid
0016 
0017 npy/NNode
0018    organic mess that need to leave behind
0019 
0020 npy/NNodeUncoincide npy/NNodeNudger
0021    some of this rather complex code will
0022    need to be moved to an "sndUncoincide"
0023 
0024 **/
0025 
0026 
0027 #include <set>
0028 #include <csignal>
0029 
0030 #include "ssys.h"
0031 #include "scuda.h"
0032 
0033 #include "sn.h"
0034 #include "stran.h"
0035 #include "stra.h"
0036 #include "OpticksCSG.h"
0037 
0038 #include "G4VSolid.hh"
0039 
0040 #include "G4Orb.hh"
0041 #include "G4Sphere.hh"
0042 #include "G4Ellipsoid.hh"
0043 #include "G4Box.hh"
0044 #include "G4Tubs.hh"
0045 #include "G4CutTubs.hh"
0046 #include "G4Polycone.hh"
0047 #include "G4Cons.hh"
0048 #include "G4Hype.hh"
0049 #include "G4MultiUnion.hh"
0050 #include "G4Torus.hh"
0051 #include "G4UnionSolid.hh"
0052 #include "G4IntersectionSolid.hh"
0053 #include "G4SubtractionSolid.hh"
0054 
0055 #include "G4BooleanSolid.hh"
0056 #include "G4RotationMatrix.hh"
0057 #include <CLHEP/Units/SystemOfUnits.h>
0058 
0059 
0060 #include "U4Polycone.h"
0061 #include "U4Transform.h"
0062 #include "U4MultiUnion.h"
0063 
0064 enum {
0065     _G4Other,
0066     _G4Orb,
0067     _G4Sphere,
0068     _G4Ellipsoid,
0069     _G4Box,
0070     _G4Tubs,
0071     _G4Polycone,
0072     _G4Cons,
0073     _G4Hype,
0074     _G4MultiUnion,
0075     _G4Torus,
0076     _G4UnionSolid,
0077     _G4IntersectionSolid,
0078     _G4SubtractionSolid,
0079     _G4DisplacedSolid,
0080     _G4CutTubs
0081  };
0082 
0083 struct U4Solid
0084 {
0085     static constexpr const char* G4Orb_               = "Orb" ;
0086     static constexpr const char* G4Sphere_            = "Sph" ;
0087     static constexpr const char* G4Ellipsoid_         = "Ell" ;
0088     static constexpr const char* G4Box_               = "Box" ;
0089     static constexpr const char* G4Tubs_              = "Tub" ;
0090     static constexpr const char* G4Polycone_          = "Pol" ;
0091     static constexpr const char* G4Cons_              = "Con" ;
0092     static constexpr const char* G4Hype_              = "Hyp" ;
0093     static constexpr const char* G4MultiUnion_        = "Mul" ;
0094     static constexpr const char* G4Torus_             = "Tor" ;
0095     static constexpr const char* G4UnionSolid_        = "Uni" ;
0096     static constexpr const char* G4IntersectionSolid_ = "Int" ;
0097     static constexpr const char* G4SubtractionSolid_  = "Sub" ;
0098     static constexpr const char* G4DisplacedSolid_    = "Dis" ;
0099     static constexpr const char* G4CutTubs_           = "TuC" ;
0100 
0101     static constexpr const char* _U4Solid__IsFlaggedLVID = "U4Solid__IsFlaggedLVID" ;
0102     static const int   IsFlaggedLVID_ ;
0103     static const char* IsFlaggedName_ ;
0104     static const char* IsFlaggedType_ ;
0105 
0106     static bool IsFlaggedLVID(int q_lvid);
0107     static bool IsFlaggedName(const char* name);
0108     static bool IsFlaggedType(const char* type);
0109 
0110     static const char*  brief_KEY ;
0111     std::string         brief() const ;
0112     std::string         desc() const ;
0113     static const char*  Name( const G4VSolid* solid) ;
0114     static unsigned     Hint( const char* name );
0115     static const char*  EType(const G4VSolid* solid) ;
0116     static int          Level(int level_, const char* name, const char* entityType );
0117 
0118     static int          Type(const char* entityType ) ;
0119     static const char*  Tag( int type ) ;
0120     const char* tag() const ;
0121 
0122 
0123     static sn* Convert(const G4VSolid* solid, int lvid, int depth, int level=-1 ) ;
0124 
0125 private:
0126     U4Solid( const G4VSolid* solid, int lvid, int depth, int level  );
0127 
0128     void init();
0129     void init_Constituents();
0130     void init_Tree();
0131     void init_Tree_Shrink();
0132 
0133     void init_Orb();
0134     void init_Sphere();
0135     void init_Ellipsoid();
0136     void init_Box();
0137     void init_Tubs();
0138     void init_Polycone();
0139     void init_Cons();
0140     void init_Hype();
0141     void init_MultiUnion();
0142     void init_Torus();
0143     void init_UnionSolid();
0144     void init_IntersectionSolid();
0145     void init_SubtractionSolid();
0146     void init_CutTubs();
0147 
0148     sn* init_Sphere_(char layer);
0149     sn* init_Cons_(char layer);
0150 
0151 
0152     static OpticksCSG_t BooleanOperator(const G4BooleanSolid* solid );
0153     static bool  IsBoolean(   const G4VSolid* solid) ;
0154     static bool  IsDisplaced( const G4VSolid* solid) ;
0155 
0156     // TODO: clean up these statics, most not needed
0157     static const G4VSolid* Left(  const G4VSolid* node) ;
0158     static const G4VSolid* Right( const G4VSolid* node) ;  // NB the G4VSolid returned might be a G4DisplacedSolid wrapping the G4VSolid
0159     static       G4VSolid* Left_(       G4VSolid* node) ;
0160     static       G4VSolid* Right_(      G4VSolid* node) ;
0161 
0162     static const G4VSolid* Moved(  G4RotationMatrix* rot, G4ThreeVector* tla, const G4VSolid* node) ;
0163     static       G4VSolid* Moved_( G4RotationMatrix* rot, G4ThreeVector* tla,       G4VSolid* node) ;
0164     static const G4VSolid* Moved( const G4VSolid* node) ;
0165     static       G4VSolid* Moved_(      G4VSolid* node) ;
0166 
0167     void init_BooleanSolid();
0168     void init_DisplacedSolid();
0169 
0170     // members
0171     const G4VSolid* solid ;
0172     const char*     name ;
0173     unsigned        hint ;
0174     const char*     entityType ;
0175     int             level ;
0176     int             type ;
0177 
0178     int             lvid ;
0179     int             depth ;   // recursion depth across different G4VSolid
0180 
0181     sn*             root ;
0182 
0183 };
0184 
0185 inline const int U4Solid::IsFlaggedLVID_ = ssys::getenvint(_U4Solid__IsFlaggedLVID, -1) ;
0186 inline bool U4Solid::IsFlaggedLVID(int q_lvid)
0187 {
0188     return IsFlaggedLVID_ == q_lvid ;
0189 }
0190 
0191 
0192 
0193 inline const char* U4Solid::IsFlaggedName_ = ssys::getenvvar("U4Solid__IsFlaggedName", "PLACEHOLDER") ;
0194 inline bool U4Solid::IsFlaggedName(const char* name)
0195 {
0196     return name && IsFlaggedName_ && strstr(name, IsFlaggedName_ ) ;
0197 }
0198 
0199 inline const char* U4Solid::IsFlaggedType_ = ssys::getenvvar("U4Solid__IsFlaggedType", "G4Dummy") ;
0200 inline bool U4Solid::IsFlaggedType(const char* type)
0201 {
0202     return type && IsFlaggedType_ && strstr(type, IsFlaggedType_ ) ;
0203 }
0204 
0205 inline const char* U4Solid::brief_KEY = "lvid/depth/tag/root/level" ;
0206 
0207 inline std::string U4Solid::brief() const
0208 {
0209     std::stringstream ss ;
0210     ss << std::setw(3) << lvid
0211        << "/" << depth
0212        << "/" << tag()
0213        << "/"  << std::setw(3) << ( root ? root->index() : -1 )
0214        << "/" << std::setw(1) << level
0215        ;
0216     std::string str = ss.str();
0217     return str ;
0218 }
0219 
0220 inline std::string U4Solid::desc() const
0221 {
0222     std::stringstream ss ;
0223     ss
0224        << brief()
0225        << " level " << level
0226        << " entityType " << std::setw(16) << ( entityType ? entityType : "-" )
0227        << " name " << ( name ? name : "-" )
0228        ;
0229     std::string str = ss.str();
0230     return str ;
0231 }
0232 
0233 inline const char* U4Solid::Name(const G4VSolid* solid)  // static
0234 {
0235     G4String _name = solid->GetName() ; // bizarre: G4VSolid::GetName returns by value, not reference
0236     const char* name = _name.c_str();
0237     return strdup(name) ;
0238 }
0239 
0240 inline unsigned U4Solid::Hint(const char* name)  // static
0241 {
0242     return sn::NameHint(name);
0243 }
0244 
0245 
0246 inline const char* U4Solid::EType(const G4VSolid* solid)  // static
0247 {
0248     G4GeometryType _etype = solid->GetEntityType();  // G4GeometryType typedef for G4String
0249     return strdup(_etype.c_str()) ;
0250 }
0251 
0252 
0253 /**
0254 U4Solid::Level
0255 ----------------
0256 
0257 IsFlaggedName
0258    returns true for solid name strings that are present within the U4Solid__IsFlaggedName envvar value
0259 
0260 IsFlaggedType
0261    returns true for solid entityType strings that are present within the U4Solid__IsFlaggedType envvar value
0262    for example::
0263 
0264         export U4Solid__IsFlaggedType=G4MultiUnion
0265 
0266 
0267 
0268 **/
0269 
0270 
0271 inline int U4Solid::Level(int level_, const char* name, const char* entityType )   // static
0272 {
0273     if(level_ > 0 ) return level_ ;
0274     if(IsFlaggedName(name) || IsFlaggedType(entityType)) return 1 ;
0275     return -1 ;
0276 }
0277 
0278 inline int U4Solid::Type(const char* name)   // static
0279 {
0280     int type = _G4Other ;
0281     if( strcmp(name, "G4Orb") == 0 )               type = _G4Orb ;
0282     if( strcmp(name, "G4Sphere") == 0 )            type = _G4Sphere ;
0283     if( strcmp(name, "G4Ellipsoid") == 0 )         type = _G4Ellipsoid ;
0284     if( strcmp(name, "G4Box") == 0 )               type = _G4Box ;
0285     if( strcmp(name, "G4Tubs") == 0 )              type = _G4Tubs ;
0286     if( strcmp(name, "G4Polycone") == 0 )          type = _G4Polycone ;
0287     if( strcmp(name, "G4Cons") == 0 )              type = _G4Cons ;
0288     if( strcmp(name, "G4Hype") == 0 )              type = _G4Hype ;
0289     if( strcmp(name, "G4MultiUnion") == 0 )        type = _G4MultiUnion ;
0290     if( strcmp(name, "G4Torus") == 0 )             type = _G4Torus ;
0291     if( strcmp(name, "G4UnionSolid") == 0 )        type = _G4UnionSolid ;
0292     if( strcmp(name, "G4SubtractionSolid") == 0 )  type = _G4SubtractionSolid ;
0293     if( strcmp(name, "G4IntersectionSolid") == 0 ) type = _G4IntersectionSolid ;
0294     if( strcmp(name, "G4DisplacedSolid") == 0 )    type = _G4DisplacedSolid ;
0295     if( strcmp(name, "G4CutTubs") == 0 )           type = _G4CutTubs ;
0296     return type ;
0297 }
0298 
0299 inline const char* U4Solid::Tag(int type)   // static
0300 {
0301     const char* tag = nullptr ;
0302     switch(type)
0303     {
0304         case _G4Orb:               tag = G4Orb_               ; break ;
0305         case _G4Sphere:            tag = G4Sphere_            ; break ;
0306         case _G4Ellipsoid:         tag = G4Ellipsoid_         ; break ;
0307         case _G4Box:               tag = G4Box_               ; break ;
0308         case _G4Tubs:              tag = G4Tubs_              ; break ;
0309         case _G4Polycone:          tag = G4Polycone_          ; break ;
0310         case _G4Cons:              tag = G4Cons_              ; break ;
0311         case _G4Hype:              tag = G4Hype_              ; break ;
0312         case _G4MultiUnion:        tag = G4MultiUnion_        ; break ;
0313         case _G4Torus:             tag = G4Torus_             ; break ;
0314         case _G4UnionSolid:        tag = G4UnionSolid_        ; break ;
0315         case _G4SubtractionSolid:  tag = G4SubtractionSolid_  ; break ;
0316         case _G4IntersectionSolid: tag = G4IntersectionSolid_ ; break ;
0317         case _G4DisplacedSolid:    tag = G4DisplacedSolid_    ; break ;
0318         case _G4CutTubs:           tag = G4CutTubs_           ; break ;
0319     }
0320     return tag ;
0321 }
0322 
0323 inline const char* U4Solid::tag() const { return Tag(type) ; }
0324 
0325 
0326 /**
0327 U4Solid::Convert
0328 -----------------
0329 
0330 Canonically invoked from U4Tree::initSolid
0331 
0332 **/
0333 
0334 
0335 inline sn* U4Solid::Convert(const G4VSolid* solid, int lvid, int depth, int level ) // static
0336 {
0337     bool flagged_LVID = IsFlaggedLVID(lvid);
0338     if(flagged_LVID) std::cout
0339         << "U4Solid::Convert"
0340         << " solid.GetName[" << solid->GetName() << "]"
0341         << " lvid " << lvid
0342         << " " << _U4Solid__IsFlaggedLVID << " : " << ( flagged_LVID ? "YES" : "NO " )
0343         << " depth " << depth
0344         << " level " << level
0345         << "\n"
0346         ;
0347 
0348 
0349     U4Solid so(solid, lvid, depth, level );
0350     return so.root ;
0351 }
0352 
0353 
0354 
0355 
0356 inline U4Solid::U4Solid(const G4VSolid* solid_, int lvid_, int depth_, int level_ )
0357     :
0358     solid(solid_),
0359     name(Name(solid_)),
0360     hint(Hint(name)),
0361     entityType(EType(solid)),
0362     level(Level(level_,name,entityType)),
0363     type(Type(entityType)),
0364     lvid(lvid_),
0365     depth(depth_),
0366     root(nullptr)
0367 {
0368     init() ;
0369 }
0370 
0371 inline void U4Solid::init()
0372 {
0373     if(level > 0) std::cerr
0374         <<  ( depth == 0 ? "\n" : "" )
0375         << "[U4Solid::init " << brief()
0376         << " " << ( depth == 0 ? brief_KEY : "" )
0377         << " " << ( depth == 0 ? name : "" )
0378         << std::endl
0379         ;
0380 
0381     init_Constituents();
0382     init_Tree() ;
0383 
0384     if(level > 0) std::cerr << "]U4Solid::init " << brief() << std::endl ;
0385 }
0386 
0387 
0388 /**
0389 U4Solid::init_Constituents
0390 --------------------------
0391 
0392 Creates sn nodes that capture the G4VSolid type and parameters.
0393 This is recursive for compound nodes like booleans.
0394 
0395 **/
0396 
0397 inline void U4Solid::init_Constituents()
0398 {
0399     switch(type)
0400     {
0401         case _G4Orb               : init_Orb()                   ; break ;
0402         case _G4Sphere            : init_Sphere()                ; break ;
0403         case _G4Ellipsoid         : init_Ellipsoid()             ; break ;
0404         case _G4Box               : init_Box()                   ; break ;
0405         case _G4Tubs              : init_Tubs()                  ; break ;
0406         case _G4Polycone          : init_Polycone()              ; break ;
0407         case _G4Cons              : init_Cons()                  ; break ;
0408         case _G4Hype              : init_Hype()                  ; break ;
0409         case _G4MultiUnion        : init_MultiUnion()            ; break ;
0410         case _G4Torus             : init_Torus()                 ; break ;
0411         case _G4UnionSolid        : init_UnionSolid()            ; break ;
0412         case _G4IntersectionSolid : init_IntersectionSolid()     ; break ;
0413         case _G4SubtractionSolid  : init_SubtractionSolid()      ; break ;
0414         case _G4DisplacedSolid    : init_DisplacedSolid()        ; break ;
0415         case _G4CutTubs           : init_CutTubs()               ; break ;
0416     }
0417 
0418     if(!root) std::cerr << "U4Solid::init_Constituents UNHANDLED SOLID TYPE " << type << "\n" << desc() << "\n" ;
0419     assert( root);
0420     root->set_hint_note(hint);
0421 }
0422 
0423 
0424 
0425 /**
0426 U4Solid::init_Tree
0427 ------------------
0428 
0429 Only invoked from the top solid after init_Constituents so the
0430 entire tree of G4VSolid has been recursively traversed
0431 when the sn::postconvert is called
0432 
0433 **/
0434 
0435 inline void U4Solid::init_Tree()
0436 {
0437     if( depth != 0 )  return ;
0438 
0439     init_Tree_Shrink();
0440 
0441     root->postconvert(lvid);
0442 }
0443 
0444 /**
0445 U4Solid::init_Tree_Shrink
0446 --------------------------
0447 
0448 try to shrink tree if prim are listnode hinted
0449 
0450 **/
0451 
0452 inline void U4Solid::init_Tree_Shrink()
0453 {
0454     if( depth != 0 )  return ;
0455 
0456     sn* root0 = root ;
0457 
0458     if(root0->has_candidate_listnode_discontiguous())
0459     {
0460         root = sn::CreateSmallerTreeWithListNode_discontiguous(root0);
0461         root->check_idx("U4Solid::init_Tree_Shrink.discontiguous");
0462     }
0463     else if(root0->has_candidate_listnode_contiguous())
0464     {
0465         root = sn::CreateSmallerTreeWithListNode_contiguous(root0);
0466         root->check_idx("U4Solid::init_Tree_Shrink.contiguous");
0467     }
0468 
0469     if(root != root0)
0470     {
0471 
0472         std::cerr << "U4Solid::init_Tree_Shrink CHANGED root with sn::CreateSmallerTreeWithListNode_discontiguous/contiguous\n" ;
0473         //std::cerr << "U4Solid::init_Tree_Shrink NOT DELETING \n" ;
0474         delete root0 ;
0475     }
0476 }
0477 
0478 
0479 
0480 
0481 
0482 
0483 
0484 
0485 
0486 
0487 
0488 
0489 
0490 
0491 
0492 
0493 
0494 inline void U4Solid::init_Orb()
0495 {
0496     const G4Orb* orb = dynamic_cast<const G4Orb*>(solid);
0497     assert(orb);
0498     double radius = orb->GetRadius()/CLHEP::mm ;
0499     root = sn::Sphere(radius) ;
0500 }
0501 
0502 inline void U4Solid::init_Sphere()
0503 {
0504     sn* outer = init_Sphere_('O') ; assert( outer ) ;
0505     sn* inner = init_Sphere_('I');
0506     root = inner == nullptr ? outer : sn::Boolean( CSG_DIFFERENCE, outer, inner ) ;
0507 }
0508 
0509 /**
0510 U4Solid::init_Sphere_
0511 -----------------------
0512 
0513 TODO: bring over phicut handling from X4Solid
0514 
0515 ::
0516 
0517       .       +z
0518                :
0519             ___:___       __ z = r,   theta = 0
0520            /   :   \
0521           /    :    \
0522       -x + - - : - - + +x
0523           \    :    /
0524            \___:___/      __ z = -r,  theta = pi
0525                :
0526               -z
0527 
0528 **/
0529 
0530 inline sn* U4Solid::init_Sphere_(char layer)
0531 {
0532     assert( layer == 'I' || layer == 'O' );
0533     const G4Sphere* sphere = dynamic_cast<const G4Sphere*>(solid);
0534     assert(sphere);
0535 
0536     double rmax = sphere->GetOuterRadius()/CLHEP::mm ;
0537     double rmin = sphere->GetInnerRadius()/CLHEP::mm ;
0538     double radius = layer == 'I' ? rmin : rmax ;
0539     if(radius == 0.) return nullptr ;
0540 
0541     double startThetaAngle = sphere->GetStartThetaAngle()/CLHEP::radian ;
0542     double deltaThetaAngle = sphere->GetDeltaThetaAngle()/CLHEP::radian ;
0543 
0544     double rTheta = startThetaAngle ;
0545     double lTheta = startThetaAngle + deltaThetaAngle ;
0546     assert( rTheta >= 0. && rTheta <= CLHEP::pi ) ;
0547     assert( lTheta >= 0. && lTheta <= CLHEP::pi ) ;
0548 
0549     double zmin = radius*std::cos(lTheta) ;
0550     double zmax = radius*std::cos(rTheta) ;
0551     assert( zmax > zmin ) ;
0552 
0553     bool z_slice = startThetaAngle > 0. || deltaThetaAngle < CLHEP::pi ;
0554 
0555     double startPhi = sphere->GetStartPhiAngle()/CLHEP::radian ;
0556     double deltaPhi = sphere->GetDeltaPhiAngle()/CLHEP::radian ;
0557     bool has_deltaPhi = startPhi != 0. || deltaPhi != 2.*CLHEP::pi  ;
0558 
0559     bool has_deltaPhi_expect = has_deltaPhi == false ;
0560     assert( has_deltaPhi_expect );
0561     if(!has_deltaPhi_expect) std::raise(SIGINT);
0562 
0563     return z_slice ? sn::ZSphere( radius, zmin, zmax ) : sn::Sphere(radius ) ;
0564 }
0565 
0566 
0567 /**
0568 U4Solid::init_Ellipsoid
0569 --------------------------
0570 
0571 Implemented using a sphere with an associated scale transform.
0572 The sphere radius is picked to be the ZSemiAxis with scale factors being::
0573 
0574     ( XSemiAxis/ZSemiAxis, YSemiAxis/ZSemiAxis, ZSemiAxis/ZSemiAxis )
0575     ( XSemiAxis/ZSemiAxis, YSemiAxis/ZSemiAxis, 1. )
0576 
0577 This means that for a PMT bulb like shape flattened in Z the scale factors might be::
0578 
0579      ( 1.346, 1.346, 1.00 )
0580 
0581 So the underlying sphere is being expanded equally in x and y directions
0582 with z unscaled.
0583 
0584 **/
0585 
0586 inline void U4Solid::init_Ellipsoid()
0587 {
0588     const G4Ellipsoid* ellipsoid = static_cast<const G4Ellipsoid*>(solid);
0589     assert(ellipsoid);
0590 
0591     double sx = ellipsoid->GetSemiAxisMax(0)/CLHEP::mm ;
0592     double sy = ellipsoid->GetSemiAxisMax(1)/CLHEP::mm ;
0593     double sz = ellipsoid->GetSemiAxisMax(2)/CLHEP::mm ;
0594 
0595     glm::tvec3<double> sc( sx/sz, sy/sz, 1. ) ;   // unity scaling in z, so z-coords are unaffected
0596     glm::tmat4x4<double> scale(1.);
0597     U4Transform::GetScaleTransform(scale, sc.x, sc.y, sc.z );
0598 
0599 
0600     double zcut1 = ellipsoid->GetZBottomCut()/CLHEP::mm ;
0601     double zcut2 = ellipsoid->GetZTopCut()/CLHEP::mm ;
0602 
0603     double zmin = zcut1 > -sz ? zcut1 : -sz ;
0604     double zmax = zcut2 <  sz ? zcut2 :  sz ;
0605     assert( zmax > zmin ) ;
0606 
0607     bool upper_cut = zmax <  sz ;
0608     bool lower_cut = zmin > -sz ;
0609     bool zslice = lower_cut || upper_cut ;
0610 
0611     if(level > 0) std::cerr
0612         << "U4Solid::init_Ellipsoid"
0613         << " upper_cut " << upper_cut
0614         << " lower_cut " << lower_cut
0615         << " zcut1 " << zcut1
0616         << " zcut2 " << zcut2
0617         << " zmin " << zmin
0618         << " zmax " << zmax
0619         << " sx " << sx
0620         << " sy " << sy
0621         << " sz " << sz
0622         << " zslice " << ( zslice ? "YES" : "NO" )
0623         << std::endl
0624         ;
0625 
0626 
0627     if( upper_cut == false && lower_cut == false )
0628     {
0629         root = sn::Sphere(sz) ;
0630     }
0631     else if( upper_cut == true && lower_cut == true )
0632     {
0633         root = sn::ZSphere(sz, zmin, zmax) ;
0634     }
0635     else if ( upper_cut == false && lower_cut == true )   // PMT mask uses this
0636     {
0637         double zmax_safe = zmax + 0.1 ;
0638         root = sn::ZSphere( sz, zmin, zmax_safe )  ;
0639 
0640     }
0641     else if ( upper_cut == true && lower_cut == false )
0642     {
0643         double zmin_safe = zmin - 0.1 ;
0644         root = sn::ZSphere( sz, zmin_safe, zmax )  ;
0645     }
0646 
0647     // zmin_safe/zmax_safe use safety offset when there is no cut
0648     // this avoids rare apex(nadir) hole bug
0649     // see notes/issues/unexpected_zsphere_miss_from_inside_for_rays_that_would_be_expected_to_intersect_close_to_apex.rst
0650     root->setXF(scale);
0651 
0652 }
0653 
0654 
0655 inline void U4Solid::init_Hype()
0656 {
0657     assert(0);
0658 }
0659 
0660 
0661 /**
0662 U4Solid::init_MultiUnion
0663 --------------------------
0664 
0665 * non-sn-primitive subs like polycone get placeholdered
0666 
0667 
0668 **/
0669 
0670 
0671 
0672 inline void U4Solid::init_MultiUnion()
0673 {
0674     const G4MultiUnion* const muni = static_cast<const G4MultiUnion*>(solid);
0675     assert(muni);
0676 
0677     int hint = CSG::HintCode(name);
0678     int type = ( hint == CSG_DISCONTIGUOUS || hint == CSG_CONTIGUOUS ) ? hint : CSG_CONTIGUOUS ;
0679 
0680     unsigned sub_num = muni->GetNumberOfSolids() ;
0681     if(level > 0 ) std::cout
0682         << "[U4Solid::init_MultiUnion"
0683         << " {" << brief() << "} "
0684         << " level " << level
0685         << " name " << name
0686         << " hint " << hint
0687         << " CSG::Name(hint) " << CSG::Name(hint)
0688         << " type " << type
0689         << " CSG::Name(type) " << CSG::Name(type)
0690         << " sub_num " << sub_num
0691         << "\n"
0692         << U4MultiUnion::Desc( muni )
0693         << "\n"
0694         ;
0695 
0696     std::vector<sn*> prims ;
0697     for( unsigned i=0 ; i < sub_num ; i++)
0698     {
0699         const G4VSolid* sub = muni->GetSolid(i);
0700         G4String sub_name = sub->GetName() ;
0701 
0702         glm::tmat4x4<double> xf(1.) ;
0703         U4Transform::GetMultiUnionItemTransform( xf, muni, i );
0704 
0705         sn* p = Convert( sub,  lvid, depth+1, level );
0706         assert(p);
0707         p->combineXF(xf);
0708         bool p_expect = p->is_primitive();
0709 
0710         if(level > 0 || !p_expect ) std::cout
0711             << "-U4Solid::init_MultiUnion"
0712             << " i " << std::setw(5) << i
0713             << " sub_name " << sub_name
0714             << " p_expect " << ( p_expect ? "YES" : "NO " )
0715             << "\n"
0716             << stra<double>::Desc(xf)
0717             << "\n"
0718             ;
0719 
0720         //assert( p_expect );
0721         if(!p_expect) p = sn::Notsupported();
0722         // sn::Notsupported is a placeholder "prim" signalling that the geometry cannot
0723         // be analytically intersected and needs to be triangulated
0724 
0725         prims.push_back(p) ;
0726     }
0727     root = sn::Compound( prims, type );
0728 
0729     if(level > 0 ) std::cout
0730         << "]U4Solid::init_MultiUnion"
0731         << " level " << level
0732         << "\n"
0733         ;
0734 
0735 
0736 
0737 }
0738 
0739 
0740 
0741 /**
0742 U4Solid::init_Torus
0743 --------------------
0744 
0745 Not planning to get analytic GPU intersect working anytime soon,
0746 just need analytic placeholder whilst use triangulated geom.
0747 
0748 **/
0749 
0750 inline void U4Solid::init_Torus()
0751 {
0752     const G4Torus* const torus = static_cast<const G4Torus*>(solid);
0753     assert(torus);
0754 
0755     double rmin_mm = torus->GetRmin()/CLHEP::mm ;
0756     double rmax_mm = torus->GetRmax()/CLHEP::mm ;
0757     double rtor_mm = torus->GetRtor()/CLHEP::mm ;
0758     double startPhi_deg = torus->GetSPhi()/CLHEP::degree ;
0759     double deltaPhi_deg = torus->GetDPhi()/CLHEP::degree ;
0760 
0761     if(level > 0) std::cout
0762        << " U4Solid::init_Torus "
0763        << " rmin_mm " << rmin_mm
0764        << " rmax_mm " << rmax_mm
0765        << " rtor_mm " << rtor_mm
0766        << " startPhi_deg " << startPhi_deg
0767        << " deltaPhi_deg " << deltaPhi_deg
0768        << "\n"
0769        ;
0770 
0771     assert( rmax_mm > rmin_mm );
0772     assert( rtor_mm > rmax_mm );
0773 
0774     root = sn::Torus(rmin_mm, rmax_mm, rtor_mm, startPhi_deg, deltaPhi_deg ) ;
0775 }
0776 
0777 
0778 
0779 
0780 inline void U4Solid::init_Box()
0781 {
0782     const G4Box* box = dynamic_cast<const G4Box*>(solid);
0783     assert(box);
0784 
0785     double fx = 2.0*box->GetXHalfLength()/CLHEP::mm ;
0786     double fy = 2.0*box->GetYHalfLength()/CLHEP::mm ;
0787     double fz = 2.0*box->GetZHalfLength()/CLHEP::mm ;
0788 
0789     root = sn::Box3(fx, fy, fz) ;
0790 }
0791 
0792 
0793 inline void U4Solid::init_Tubs()
0794 {
0795     const G4Tubs* tubs = dynamic_cast<const G4Tubs*>(solid);
0796     assert(tubs);
0797 
0798     double hz = tubs->GetZHalfLength()/CLHEP::mm ;
0799     double rmax = tubs->GetOuterRadius()/CLHEP::mm ;
0800     double rmin = tubs->GetInnerRadius()/CLHEP::mm ;
0801     bool has_inner = rmin > 0. ;
0802 
0803     sn* outer = sn::Cylinder(rmax, -hz, hz );
0804 
0805 
0806     if(has_inner == false)
0807     {
0808         root = outer ;
0809     }
0810     else
0811     {
0812         bool   do_nudge_inner = true ;
0813         double nudge_inner = 0.01 ;
0814         double dz = do_nudge_inner ? hz*nudge_inner : 0. ;
0815 
0816         sn* inner = sn::Cylinder(rmin, -(hz+dz), hz+dz );
0817         root = sn::Boolean( CSG_DIFFERENCE, outer, inner );
0818     }
0819 
0820 }
0821 
0822 /**
0823 U4Solid::init_CutTubs
0824 -----------------------
0825 
0826 ParameterHolder : Not planning to get analytic GPU intersect
0827 working anytime soon just need analytic placeholder
0828 whilst use triangulated geom.
0829 
0830 **/
0831 
0832 inline void U4Solid::init_CutTubs()
0833 {
0834     const G4CutTubs* tubs = dynamic_cast<const G4CutTubs*>(solid);
0835     assert(tubs);
0836 
0837     double hz = tubs->GetZHalfLength()/CLHEP::mm ;
0838     double rmax = tubs->GetOuterRadius()/CLHEP::mm ;
0839     double rmin = tubs->GetInnerRadius()/CLHEP::mm ;
0840     bool has_inner = rmin > 0. ;
0841 
0842     G4ThreeVector lo = tubs->GetLowNorm();
0843     G4ThreeVector hi = tubs->GetHighNorm();
0844 
0845     double pz_nrm_x = hi.x();
0846     double pz_nrm_y = hi.y();
0847     double pz_nrm_z = hi.z();
0848     assert( pz_nrm_z > 0. );  // expect outward normal away from pz top edge
0849 
0850     double nz_nrm_x = lo.x();
0851     double nz_nrm_y = lo.y();
0852     double nz_nrm_z = lo.z();
0853     assert( nz_nrm_z < 0. );  // expect outward normal away from nz bot edge
0854 
0855     assert( pz_nrm_y == 0. );
0856     assert( nz_nrm_y == 0. );
0857 
0858     sn* outer = sn::CutCylinder(rmax, hz,
0859           pz_nrm_x, pz_nrm_y, pz_nrm_z,
0860           nz_nrm_x, nz_nrm_y, nz_nrm_z );
0861 
0862     if(has_inner == false)
0863     {
0864         root = outer ;
0865     }
0866     else
0867     {
0868         bool   do_nudge_inner = true ;
0869         double nudge_inner = 0.01 ;  // HMM: MAYBE NEED MORE WITH CUTS ?
0870         double dz = do_nudge_inner ? hz*nudge_inner : 0. ;
0871 
0872         sn* inner = sn::CutCylinder(rmin, hz+dz,
0873           pz_nrm_x, pz_nrm_y, pz_nrm_z,
0874           nz_nrm_x, nz_nrm_y, nz_nrm_z );
0875 
0876         root = sn::Boolean( CSG_DIFFERENCE, outer, inner );
0877     }
0878 }
0879 
0880 
0881 
0882 
0883 inline void U4Solid::init_Polycone()
0884 {
0885     const G4Polycone* polycone = dynamic_cast<const G4Polycone*>(solid);
0886     assert(polycone);
0887     root = U4Polycone::Convert(polycone, lvid, depth, level );
0888 
0889     if(level > 0 ) std::cerr
0890         << "U4Solid::init_Polycone"
0891         << " level " << level
0892         << std::endl
0893         << desc()
0894         << std::endl
0895         ;
0896 }
0897 
0898 
0899 
0900 inline sn* U4Solid::init_Cons_(char layer)
0901 {
0902     assert( layer == 'I' || layer == 'O' );
0903     const G4Cons* cone = dynamic_cast<const G4Cons*>(solid);
0904     assert(cone);
0905 
0906     double rmax1    = cone->GetOuterRadiusMinusZ()/CLHEP::mm ;
0907     double rmax2    = cone->GetOuterRadiusPlusZ()/CLHEP::mm  ;
0908 
0909     double rmin1    = cone->GetInnerRadiusMinusZ()/CLHEP::mm ;
0910     double rmin2    = cone->GetInnerRadiusPlusZ()/CLHEP::mm  ;
0911 
0912     double hz       = cone->GetZHalfLength()/CLHEP::mm   ;
0913 
0914     //double startPhi = cone->GetStartPhiAngle()/CLHEP::degree ;
0915     //double deltaPhi = cone->GetDeltaPhiAngle()/CLHEP::degree ;
0916 
0917     double r1 = layer == 'I' ? rmin1 : rmax1 ;
0918     double r2 = layer == 'I' ? rmin2 : rmax2 ;
0919     double z1 = -hz ;
0920     double z2 = hz ;
0921 
0922     bool invalid =  r1 == 0. && r2 == 0. ;
0923 
0924     return invalid ? nullptr : sn::Cone( r1, z1, r2, z2 ) ;
0925 }
0926 
0927 inline void U4Solid::init_Cons()
0928 {
0929     sn* outer = init_Cons_('O');  assert( outer );
0930     sn* inner = init_Cons_('I');
0931     root = inner == nullptr ? outer : sn::Boolean( CSG_DIFFERENCE, outer, inner ) ;
0932 
0933     if(level > 0 ) std::cerr
0934         << "U4Solid::init_Cons"
0935         << " level " << level
0936         << desc()
0937         ;
0938 
0939 }
0940 
0941 
0942 
0943 
0944 inline void U4Solid::init_UnionSolid()
0945 {
0946     init_BooleanSolid();
0947 }
0948 inline void U4Solid::init_IntersectionSolid()
0949 {
0950     init_BooleanSolid();
0951 }
0952 inline void U4Solid::init_SubtractionSolid()
0953 {
0954     init_BooleanSolid();
0955 }
0956 
0957 inline OpticksCSG_t U4Solid::BooleanOperator(const G4BooleanSolid* solid ) // static
0958 {
0959     OpticksCSG_t _operator = CSG_ZERO ;
0960     if      (dynamic_cast<const G4IntersectionSolid*>(solid)) _operator = CSG_INTERSECTION ;
0961     else if (dynamic_cast<const G4SubtractionSolid*>(solid))  _operator = CSG_DIFFERENCE ;
0962     else if (dynamic_cast<const G4UnionSolid*>(solid))        _operator = CSG_UNION ;
0963     assert( _operator != CSG_ZERO ) ;
0964     return _operator ;
0965 }
0966 
0967 
0968 inline bool U4Solid::IsBoolean(const G4VSolid* solid) // static
0969 {
0970     return dynamic_cast<const G4BooleanSolid*>(solid) != nullptr ;
0971 }
0972 inline bool U4Solid::IsDisplaced(const G4VSolid* solid) // static
0973 {
0974     return dynamic_cast<const G4DisplacedSolid*>(solid) != nullptr ;
0975 }
0976 inline const G4VSolid* U4Solid::Left(const G4VSolid* solid ) // static
0977 {
0978     return IsBoolean(solid) ? solid->GetConstituentSolid(0) : nullptr ;
0979 }
0980 inline const G4VSolid* U4Solid::Right(const G4VSolid* solid ) // static
0981 {
0982     return IsBoolean(solid) ? solid->GetConstituentSolid(1) : nullptr ;
0983 }
0984 inline G4VSolid* U4Solid::Left_(G4VSolid* solid ) // static
0985 {
0986     return IsBoolean(solid) ? solid->GetConstituentSolid(0) : nullptr ;
0987 }
0988 inline G4VSolid* U4Solid::Right_(G4VSolid* solid ) // static
0989 {
0990     return IsBoolean(solid) ? solid->GetConstituentSolid(1) : nullptr ;
0991 }
0992 
0993 /**
0994 U4Solid::Moved
0995 ---------------
0996 
0997 When node isa G4DisplacedSolid sets the rotation and translation and returns the constituentMovedSolid
0998 otherwise returns the input node.
0999 
1000 cf X4Transform3D::GetDisplacementTransform does the same
1001 
1002 **/
1003 inline const G4VSolid* U4Solid::Moved( G4RotationMatrix* rot, G4ThreeVector* tla, const G4VSolid* node )  // static
1004 {
1005     const G4DisplacedSolid* disp = dynamic_cast<const G4DisplacedSolid*>(node) ;
1006     if(disp)
1007     {
1008         if(rot) *rot = disp->GetFrameRotation();
1009         if(tla) *tla = disp->GetObjectTranslation();
1010         // HMM: looks a bit fishy being inconsistent here : until look at g4-cls G4DisplacedSolid
1011     }
1012     return disp ? disp->GetConstituentMovedSolid() : node  ;
1013 }
1014 inline G4VSolid* U4Solid::Moved_( G4RotationMatrix* rot, G4ThreeVector* tla, G4VSolid* node )  // static
1015 {
1016     G4DisplacedSolid* disp = dynamic_cast<G4DisplacedSolid*>(node) ;
1017     if(disp)
1018     {
1019         if(rot) *rot = disp->GetFrameRotation();
1020         if(tla) *tla = disp->GetObjectTranslation();
1021         // HMM: looks a bit fishy being inconsistent here : until look at g4-cls G4DisplacedSolid
1022     }
1023     return disp ? disp->GetConstituentMovedSolid() : node  ;
1024 }
1025 
1026 inline const G4VSolid* U4Solid::Moved( const G4VSolid* node )  // static
1027 {
1028     const G4DisplacedSolid* disp = dynamic_cast<const G4DisplacedSolid*>(node) ;
1029     return disp ? disp->GetConstituentMovedSolid() : node  ;
1030 }
1031 
1032 inline G4VSolid* U4Solid::Moved_( G4VSolid* node )  // static
1033 {
1034     G4DisplacedSolid* disp = dynamic_cast<G4DisplacedSolid*>(node) ;
1035     return disp ? disp->GetConstituentMovedSolid() : node  ;
1036 }
1037 
1038 
1039 inline void U4Solid::init_BooleanSolid()
1040 {
1041     const G4BooleanSolid* boo = static_cast<const G4BooleanSolid*>(solid);
1042     assert(boo);
1043 
1044     OpticksCSG_t op = BooleanOperator(boo);
1045     G4VSolid*    left  = const_cast<G4VSolid*>(boo->GetConstituentSolid(0));
1046     G4VSolid*    right = const_cast<G4VSolid*>(boo->GetConstituentSolid(1));
1047 
1048     bool is_left_displaced = dynamic_cast<G4DisplacedSolid*>(left) != nullptr ;
1049     bool is_right_displaced = dynamic_cast<G4DisplacedSolid*>(right) != nullptr ;
1050 
1051     bool is_left_displaced_expect = is_left_displaced == false ;
1052     assert( is_left_displaced_expect && "not expecting left displacement " );
1053     if(!is_left_displaced_expect) std::raise(SIGINT);
1054 
1055     sn* l = Convert( left,  lvid, depth+1, level );
1056     sn* r = Convert( right, lvid, depth+1, level );
1057 
1058     if(l->xform && level > 0) std::cout
1059         << "U4Solid::init_BooleanSolid "
1060         << " observe transform on left node "
1061         << " l.xform " << l->xform->desc()
1062         << std::endl
1063         ;
1064 
1065     if(is_right_displaced) assert( r->xform && "expecting transform on right displaced node " );
1066 
1067     root = sn::Boolean( op, l, r );
1068 }
1069 
1070 
1071 /**
1072 U4Solid::init_DisplacedSolid
1073 ------------------------------
1074 
1075 When booleans are created with transforms the right hand side
1076 gets wrapped into a DisplacedSolid.::
1077 
1078     U4Solid::init SUCCEEDED desc: U4Solid::desc solid Y type Y U4Solid::Tag(type) Box root 57
1079     U4Solid::init SUCCEEDED desc: U4Solid::desc solid Y type Y U4Solid::Tag(type) Box root 58
1080     U4Solid::init SUCCEEDED desc: U4Solid::desc solid Y type Y U4Solid::Tag(type) Dis root 58
1081     U4Solid::init SUCCEEDED desc: U4Solid::desc solid Y type Y U4Solid::Tag(type) Uni root 59
1082     U4Solid::init SUCCEEDED desc: U4Solid::desc solid Y type Y U4Solid::Tag(type) Box root 60
1083     U4Solid::init SUCCEEDED desc: U4Solid::desc solid Y type Y U4Solid::Tag(type) Dis root 60
1084 
1085 
1086 * "Dis" root not incremented ?
1087 
1088    * thats because Dis are "internal" nodes for the Geant4 model that just carry the transform
1089      they dont actually correspond to a separate constituent
1090 
1091 Dis ARE FUNNY NODES THAT JUST ACT TO HOLD THE TRANSFORM
1092 **/
1093 inline void U4Solid::init_DisplacedSolid()
1094 {
1095     const G4DisplacedSolid* disp = static_cast<const G4DisplacedSolid*>(solid);
1096     assert(disp);
1097 
1098     glm::tmat4x4<double> xf(1.) ;
1099     U4Transform::GetDispTransform( xf, disp );
1100 
1101     const G4VSolid* moved = Moved(solid);
1102     assert(moved);
1103 
1104     bool single_disp = dynamic_cast<const G4DisplacedSolid*>(moved) == nullptr ;
1105     assert(single_disp && "only single disp is expected" );
1106     if(!single_disp) std::raise(SIGINT);
1107 
1108     root = Convert( moved, lvid, depth+1, level );
1109     root->combineXF(xf);
1110 }
1111 
1112
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!