Back to home page

EIC code displayed by LXR
 
 

    

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

0001 #include <cassert>
0002 #include <iostream>
0003 #include <sstream>
0004 #include <iomanip>
0005 #include <cstdlib>
0006 
0007 #include "SLOG.hh"
0008 
0009 #include "G4SolidStore.hh"
0010 #include "G4UnionSolid.hh"
0011 #include "G4IntersectionSolid.hh"
0012 #include "G4SubtractionSolid.hh"
0013 #include "G4Ellipsoid.hh"
0014 #include "G4Tubs.hh"
0015 #include "G4Polycone.hh"
0016 #include "G4Torus.hh"
0017 #include "G4Box.hh"
0018 #include "G4Orb.hh"
0019 #include "G4MultiUnion.hh"
0020 #include "G4Sphere.hh"
0021 #include "G4Version.hh"
0022 
0023 #include "G4RotationMatrix.hh"
0024 
0025 #include "scanvas.h"
0026 #include "U4SolidTree.hh"
0027 
0028 
0029 const bool U4SolidTree::verbose = getenv("U4SolidTree_verbose") != nullptr ; 
0030 
0031 const plog::Severity U4SolidTree::LEVEL = SLOG::EnvLevel("U4SolidTree", "DEBUG") ; 
0032 
0033 G4VSolid* U4SolidTree::ApplyZCutTree( const G4VSolid* original, double zcut ) // static
0034 {
0035     if(verbose)
0036     std::cout << "[ U4SolidTree::ApplyZCutTree zcut " << zcut << " original.GetName " << original->GetName() << std::endl ; 
0037 
0038     U4SolidTree* zs = new U4SolidTree(original); 
0039     zs->apply_cut( zcut );  
0040     //zs->dump("U4SolidTree::ApplyZCutTree"); 
0041 
0042     if(verbose)
0043     std::cout << "] U4SolidTree::ApplyZCutTree  zs.root.GetName " << zs->root->GetName()  << std::endl ; 
0044     return zs->root ; 
0045 }
0046 
0047 void U4SolidTree::Draw(const G4VSolid* original, const char* msg ) // static
0048 {
0049     LOG(LEVEL) << "[" ; 
0050     if(original == nullptr )
0051     {
0052         std::cout << "ERROR U4SolidTree::Draw got nullptr original : msg " << msg << std::endl ; 
0053         return ; 
0054     }    
0055 
0056 
0057     U4SolidTree* zs = new U4SolidTree(original); 
0058     zs->draw(msg);
0059     zs->dumpNames(msg); 
0060     zs->zdump(msg); 
0061 
0062     LOG(LEVEL) << "]" ; 
0063 }
0064 
0065 U4SolidTree::U4SolidTree(const G4VSolid* original_ ) 
0066     :
0067     original(original_),
0068     root(DeepClone(original_)),
0069     edited(false),
0070     parent_map(new std::map<const G4VSolid*, const G4VSolid*>),
0071 
0072     in_map(    new std::map<const G4VSolid*, int>),
0073     rin_map(   new std::map<const G4VSolid*, int>),
0074     pre_map(   new std::map<const G4VSolid*, int>),
0075     rpre_map(  new std::map<const G4VSolid*, int>),
0076     post_map(  new std::map<const G4VSolid*, int>),
0077     rpost_map( new std::map<const G4VSolid*, int>),
0078 
0079     zcls_map(  new std::map<const G4VSolid*, int>),
0080     depth_map( new std::map<const G4VSolid*, int>),
0081     mkr_map(   new std::map<const G4VSolid*, char>),
0082 
0083     width(0),
0084     height(0), 
0085     extra_width(1),    // +1 for annotations to the right
0086     extra_height(1+1), // +1 as height zero tree is still one node, +1 for annotation  
0087     canvas(  new scanvas(width+extra_width, height+extra_height, 8, 5) ),
0088     names( new std::vector<std::string>),
0089     nameprefix(nullptr),
0090     inorder( new std::vector<const G4VSolid*> ),
0091     rinorder( new std::vector<const G4VSolid*> ),
0092     preorder( new std::vector<const G4VSolid*> ),
0093     rpreorder( new std::vector<const G4VSolid*> ),
0094     postorder( new std::vector<const G4VSolid*> ),
0095     rpostorder( new std::vector<const G4VSolid*> ),
0096     crux( new std::vector<G4VSolid*> )
0097 {
0098     init(); 
0099 }
0100 
0101 
0102 
0103 
0104 
0105 void U4SolidTree::init()
0106 {
0107     if(verbose) std::cout << "U4SolidTree::init" << std::endl ; 
0108     instrumentTree(); 
0109 }
0110 
0111 void U4SolidTree::dump(const char* msg) const 
0112 {
0113     dumpTree(msg); 
0114     dumpUp(msg); 
0115 }
0116 
0117 void U4SolidTree::instrumentTree()
0118 {
0119     if(verbose) std::cout << "U4SolidTree::instrumentTree parent " << std::endl ; 
0120     parent_map->clear(); 
0121     (*parent_map)[root] = nullptr ;  // root has no parent by definition
0122     parent_r( root, 0 ); 
0123 
0124     if(verbose) std::cout << "U4SolidTree::instrumentTree depth " << std::endl ; 
0125     depth_map->clear(); 
0126     depth_r(root, 0 ); 
0127 
0128     if(verbose) std::cout << "U4SolidTree::instrumentTree inorder " << std::endl ; 
0129     inorder->clear();
0130     inorder_r(root, 0 ); 
0131 
0132     if(verbose) std::cout << "U4SolidTree::instrumentTree rinorder " << std::endl ; 
0133     rinorder->clear();
0134     rinorder_r(root, 0 ); 
0135 
0136     if(verbose) std::cout << "U4SolidTree::instrumentTree preorder " << std::endl ; 
0137     preorder->clear();
0138     preorder_r(root, 0 ); 
0139 
0140     if(verbose) std::cout << "U4SolidTree::instrumentTree rpreorder " << std::endl ; 
0141     rpreorder->clear();
0142     rpreorder_r(root, 0 ); 
0143 
0144     if(verbose) std::cout << "U4SolidTree::instrumentTree postorder " << std::endl ; 
0145     postorder->clear(); 
0146     postorder_r(root, 0 ); 
0147 
0148     if(verbose) std::cout << "U4SolidTree::instrumentTree rpostorder " << std::endl ; 
0149     rpostorder->clear(); 
0150     rpostorder_r(root, 0 ); 
0151 
0152     if(verbose) std::cout << "U4SolidTree::instrumentTree names" << std::endl ; 
0153     names->clear(); 
0154     collectNames_inorder_r( root, 0 ); 
0155 
0156     if(verbose)
0157     {
0158         std::cout << "U4SolidTree::instrumentTree names.size " << names->size() << std::endl ; 
0159         for(unsigned i=0 ; i < names->size() ; i++) std::cout << "[" << (*names)[i] << "]" << std::endl ;  
0160     }
0161     if(verbose) std::cout << "U4SolidTree::instrumentTree nameprefix" << std::endl ; 
0162     nameprefix = CommonPrefix(names); 
0163     if(verbose) std::cout << "U4SolidTree::instrumentTree nameprefix [" << nameprefix << "]" << std::endl ; 
0164 
0165     if(verbose) std::cout << "U4SolidTree::instrumentTree [ original_num_node " << std::endl ; 
0166     int original_num_node = NumNode_r(original, 0); 
0167     if(verbose) std::cout << "U4SolidTree::instrumentTree ] original_num_node : " << original_num_node  << std::endl ; 
0168 
0169     if(verbose) std::cout << "U4SolidTree::instrumentTree [ root_num_node " << std::endl ; 
0170     int root_num_node = NumNode_r(root, 0); 
0171     if(verbose) std::cout << "U4SolidTree::instrumentTree ] root_num_node : " << root_num_node << std::endl ; 
0172 
0173     int depth_size = depth_map->size(); 
0174     int inorder_size = inorder->size(); 
0175     int rinorder_size = rinorder->size(); 
0176     int preorder_size = preorder->size(); 
0177     int rpreorder_size = rpreorder->size(); 
0178     int postorder_size = postorder->size(); 
0179     int rpostorder_size = rpostorder->size(); 
0180 
0181     if(verbose) std::cout 
0182         << "U4SolidTree::instrumentTree"
0183         << " depth_size " << depth_size
0184         << " inorder_size " << inorder_size
0185         << " rinorder_size " << rinorder_size
0186         << " preorder_size " << preorder_size
0187         << " rpreorder_size " << rpreorder_size
0188         << " postorder_size " << postorder_size
0189         << " rpostorder_size " << rpostorder_size
0190     //    << " original_num_node " << original_num_node
0191         << " root_num_node " << root_num_node
0192         << std::endl 
0193         ;
0194 
0195     assert( depth_size == root_num_node );   
0196     assert( inorder_size == root_num_node );   
0197     assert( rinorder_size == root_num_node );   
0198     assert( preorder_size == root_num_node );   
0199     assert( rpreorder_size == root_num_node );   
0200     assert( postorder_size == root_num_node );   
0201     assert( rpostorder_size == root_num_node );   
0202 
0203     if( edited == false )
0204     {
0205         assert( original_num_node == root_num_node ); 
0206     }
0207 
0208     width = root_num_node ; 
0209     height = maxdepth() ; 
0210 
0211     if(verbose) printf("U4SolidTree::instrumentTree width %d height %d \n", width, height );     
0212     canvas->resize( width+extra_width, height+extra_height );    
0213 }
0214 
0215 
0216 const char* U4SolidTree::CommonPrefix(const std::vector<std::string>* a) // static
0217 {
0218     std::vector<std::string> aa(*a); 
0219 
0220     if(verbose)
0221     {
0222         std::cout << "U4SolidTree::CommonPrefix aa.size (before sort) " << aa.size() << std::endl ;  
0223         for(unsigned i=0 ; i < aa.size() ; i++) std::cout << "[" << aa[i] << "]" << std::endl ;  
0224     }
0225 
0226     std::sort( aa.begin(), aa.end() );  
0227 
0228     if(verbose)
0229     {
0230         std::cout << "U4SolidTree::CommonPrefix aa.size (after sort) " << aa.size() << std::endl ;  
0231         for(unsigned i=0 ; i < aa.size() ; i++) std::cout << "[" << aa[i] << "]" << std::endl ;  
0232     }
0233 
0234     const std::string& s1 = aa[0] ; 
0235     const std::string& s2 = aa[aa.size()-1] ; 
0236 
0237     std::string prefix(s1) ; 
0238     for(unsigned i=0 ; i < s1.size() ; i++) 
0239     {
0240         if( s1[i] != s2[i] ) 
0241         {
0242             prefix = s1.substr(0,i) ; 
0243             break ; 
0244         }
0245     }
0246     const char* prefix_ = strdup(prefix.c_str()) ; 
0247 
0248     if(verbose)
0249     {
0250         std::cout << "U4SolidTree::CommonPrefix s1 " << s1 << std::endl ;  
0251         std::cout << "U4SolidTree::CommonPrefix s2 " << s2 << std::endl ;  
0252         std::cout << "U4SolidTree::CommonPrefix prefix " << prefix << std::endl ;  
0253         std::cout << "U4SolidTree::CommonPrefix prefix_ " << prefix_ << std::endl ;  
0254     }
0255 
0256     return prefix_  ; 
0257 } 
0258 
0259 /**
0260 U4SolidTree::parent_r
0261 -----------------------
0262 
0263 Note that the parent_map uses the raw constituent G4DisplacedSolid 
0264 rather than the moved G4VSolid that it points to in order to have 
0265 treewise access to the transform up the lineage. 
0266 
0267 **/
0268 void U4SolidTree::parent_r( const G4VSolid* node_, int depth)
0269 {
0270     if( node_ == nullptr ) return ; 
0271     const G4VSolid* node = Moved(node_ ); 
0272 
0273     if(IsBoolean(node))
0274     {
0275         const G4VSolid* l = Left(node) ; 
0276         const G4VSolid* r = Right(node) ; 
0277 
0278         parent_r( l, depth+1 );  
0279         parent_r( r, depth+1 );  
0280 
0281         // postorder visit 
0282         (*parent_map)[l] = node_ ; 
0283         (*parent_map)[r] = node_ ; 
0284     }
0285 }
0286 
0287 void U4SolidTree::depth_r(const G4VSolid* node_, int depth)
0288 {
0289     if( node_ == nullptr ) return ; 
0290     const G4VSolid* node = Moved(node_ ); 
0291 
0292     depth_r(Left(node), depth+1) ; 
0293     depth_r(Right(node), depth+1) ; 
0294 
0295     (*depth_map)[node_] = depth ;   
0296 }
0297 
0298 void U4SolidTree::inorder_r(const G4VSolid* node_, int depth)
0299 {
0300     if( node_ == nullptr ) return ; 
0301     const G4VSolid* node = Moved(node_ ); 
0302 
0303     inorder_r(Left(node), depth+1) ; 
0304 
0305     (*in_map)[node_] = inorder->size();   
0306     inorder->push_back(node_) ; 
0307 
0308     inorder_r(Right(node), depth+1) ; 
0309 } 
0310 
0311 void U4SolidTree::rinorder_r(const G4VSolid* node_, int depth)
0312 {
0313     if( node_ == nullptr ) return ; 
0314     const G4VSolid* node = Moved(node_ ); 
0315 
0316     rinorder_r(Right(node), depth+1) ; 
0317 
0318     (*rin_map)[node_] = rinorder->size();   
0319     rinorder->push_back(node_) ; 
0320 
0321     rinorder_r(Left(node), depth+1) ; 
0322 } 
0323 
0324 void U4SolidTree::preorder_r(const G4VSolid* node_, int depth)
0325 {
0326     if( node_ == nullptr ) return ; 
0327     const G4VSolid* node = Moved(node_ ); 
0328 
0329     (*pre_map)[node_] = preorder->size();   
0330     preorder->push_back(node_) ; 
0331 
0332     preorder_r(Left(node), depth+1) ; 
0333     preorder_r(Right(node), depth+1) ; 
0334 } 
0335 
0336 void U4SolidTree::rpreorder_r(const G4VSolid* node_, int depth)
0337 {
0338     if( node_ == nullptr ) return ; 
0339     const G4VSolid* node = Moved(node_ ); 
0340 
0341     (*rpre_map)[node_] = rpreorder->size();   
0342     rpreorder->push_back(node_) ; 
0343 
0344     rpreorder_r(Right(node), depth+1) ; 
0345     rpreorder_r(Left(node), depth+1) ; 
0346 } 
0347 
0348 void U4SolidTree::postorder_r(const G4VSolid* node_, int depth)
0349 {
0350     if( node_ == nullptr ) return ; 
0351     const G4VSolid* node = Moved(node_ ); 
0352 
0353     postorder_r(Left(node), depth+1) ; 
0354     postorder_r(Right(node), depth+1) ; 
0355 
0356     (*post_map)[node_] = postorder->size();   
0357     postorder->push_back(node_) ; 
0358 }
0359 
0360 void U4SolidTree::rpostorder_r(const G4VSolid* node_, int depth)
0361 {
0362     if( node_ == nullptr ) return ; 
0363     const G4VSolid* node = Moved(node_ ); 
0364 
0365     rpostorder_r(Right(node), depth+1) ; 
0366     rpostorder_r(Left(node), depth+1) ; 
0367 
0368     (*rpost_map)[node_] = rpostorder->size();   
0369     rpostorder->push_back(node_) ; 
0370 }
0371 
0372 const G4VSolid* U4SolidTree::parent(  const G4VSolid* node ) const { return parent_map->count(node) == 1 ? (*parent_map)[node] : nullptr ; }
0373 G4VSolid*       U4SolidTree::parent_( const G4VSolid* node ) const { const G4VSolid* p = parent(node); return const_cast<G4VSolid*>(p) ; }
0374 
0375 int U4SolidTree::depth( const G4VSolid* node_) const { return (*depth_map)[node_] ; }
0376 int U4SolidTree::in(    const G4VSolid* node_) const { return (*in_map)[node_] ; }
0377 int U4SolidTree::rin(   const G4VSolid* node_) const { return (*rin_map)[node_] ; }
0378 int U4SolidTree::pre(   const G4VSolid* node_) const { return (*pre_map)[node_] ; }
0379 int U4SolidTree::rpre(  const G4VSolid* node_) const { return (*rpre_map)[node_] ; }
0380 int U4SolidTree::post(  const G4VSolid* node_) const { return (*post_map)[node_] ; }
0381 int U4SolidTree::rpost( const G4VSolid* node_) const { return (*rpost_map)[node_] ; }
0382 
0383 /**
0384 U4SolidTree::index
0385 ---------------
0386 
0387 Returns the index of a node within various traversal orders, 
0388 obtained by lookups on the maps collected by instrumentTree.
0389 
0390 
0391 IN inorder
0392    left-to-right index (aka side) 
0393 
0394 RIN reverse inorder
0395    right-to-left index 
0396 
0397 PRE preorder
0398    with left unbalanced tree this does a anti-clockwise rotation starting from root
0399    [by observation is opposite to RPOST]
0400 
0401 RPRE reverse preorder
0402    with left unbalanced does curios zig-zag that looks to be exact opposite of 
0403    the familiar postorder traversal, kinda undo of the postorder.
0404    [by observation is opposite to POST]
0405 
0406    **COULD BE USEFUL FOR TREE EDITING FROM RIGHT** 
0407 
0408 POST postorder
0409    old familar traversal starting from bottom left 
0410    [by observation is opposite to RPRE]
0411 
0412 RPOST reverse postorder
0413    with left unbalanced does a clockwise cycle starting
0414    from rightmost visiting all prim before getting 
0415    to any operators
0416    [by observation is opposite to PRE]
0417 
0418 **/
0419 
0420 int U4SolidTree::index( const G4VSolid* n, int mode ) const  
0421 {
0422     int idx = -1 ; 
0423     switch(mode)
0424     {
0425         case IN:    idx = in(n)    ; break ; 
0426         case RIN:   idx = rin(n)   ; break ; 
0427         case PRE:   idx = pre(n)   ; break ; 
0428         case RPRE:  idx = rpre(n)  ; break ; 
0429         case POST:  idx = post(n)  ; break ; 
0430         case RPOST: idx = rpost(n) ; break ;  
0431     }
0432     return idx ; 
0433 }
0434 
0435 const char* U4SolidTree::IN_ = "IN" ; 
0436 const char* U4SolidTree::RIN_ = "RIN" ; 
0437 const char* U4SolidTree::PRE_ = "PRE" ; 
0438 const char* U4SolidTree::RPRE_ = "RPRE" ; 
0439 const char* U4SolidTree::POST_ = "POST" ; 
0440 const char* U4SolidTree::RPOST_ = "RPOST" ; 
0441 
0442 const char* U4SolidTree::OrderName(int mode) // static
0443 {
0444     const char* s = nullptr ; 
0445     switch(mode)
0446     {
0447         case IN:    s = IN_    ; break ; 
0448         case RIN:   s = RIN_   ; break ; 
0449         case PRE:   s = PRE_   ; break ; 
0450         case RPRE:  s = RPRE_  ; break ; 
0451         case POST:  s = POST_  ; break ; 
0452         case RPOST: s = RPOST_ ; break ;  
0453     }
0454     return s ; 
0455 }
0456 
0457 /**
0458 U4SolidTree::zcls
0459 --------------
0460 
0461 Returns zcls value associated with a node, when *move* is true any 
0462 G4DisplacedSolid *node_* are dereferenced to get to the moved solid *node* 
0463 within.
0464 
0465 Am veering towards using move:false as standard because can 
0466 always get the moved solid from the G4DisplacedSolid but not vice versa.
0467 Note that it is a technicality : either could work its just a case of
0468 which is most convenient. 
0469 
0470 **/
0471 
0472 int U4SolidTree::zcls( const G4VSolid* node_ ) const 
0473 { 
0474     return zcls_map->count(node_) == 1 ? (*zcls_map)[node_] : UNDEFINED ; 
0475 }
0476 
0477 void U4SolidTree::set_zcls( const G4VSolid* node_, int zc )
0478 {
0479     (*zcls_map)[node_] = zc  ;
0480 } 
0481 
0482 char U4SolidTree::mkr( const G4VSolid* node_) const 
0483 { 
0484     return mkr_map->count(node_) == 1 ? (*mkr_map)[node_] : ' ' ; 
0485 }
0486 
0487 void U4SolidTree::set_mkr( const G4VSolid* node_, char mk )
0488 {
0489     (*mkr_map)[node_] = mk  ;
0490 } 
0491 
0492 
0493 bool U4SolidTree::is_include( const G4VSolid* node_) const 
0494 {
0495     return zcls(node_) == INCLUDE ; 
0496 }
0497 bool U4SolidTree::is_exclude( const G4VSolid* node_) const 
0498 {
0499     return zcls(node_) == EXCLUDE ; 
0500 }
0501 
0502 bool U4SolidTree::is_exclude_include( const G4VSolid* node_) const 
0503 {
0504     if(!IsBoolean(node_)) return false ; 
0505     const G4VSolid* left_  = Left(node_); 
0506     const G4VSolid* right_ = Right(node_); 
0507     return is_exclude(left_) &&  is_include(right_) ; 
0508 }
0509 
0510 bool U4SolidTree::is_include_exclude( const G4VSolid* node_) const 
0511 {
0512     if(!IsBoolean(node_)) return false ; 
0513     const G4VSolid* left_  = Left(node_); 
0514     const G4VSolid* right_ = Right(node_); 
0515     return is_include(left_) &&  is_exclude(right_) ; 
0516 }
0517 
0518 bool U4SolidTree::is_crux( const G4VSolid* node_ ) const 
0519 {
0520     if(!IsBoolean(node_)) return false ; 
0521     const G4VSolid* left_  = Left(node_); 
0522     const G4VSolid* right_ = Right(node_); 
0523     bool exclude_include = zcls(left_) == EXCLUDE &&  zcls(right_) == INCLUDE ;
0524     bool include_exclude = zcls(left_) == INCLUDE &&  zcls(right_) == EXCLUDE ; 
0525     return exclude_include ^ include_exclude ;   // XOR one or other, not both 
0526 }
0527 
0528 int U4SolidTree::num_prim_r(const G4VSolid* node_) const
0529 {   
0530     if(!node_) return 0 ; 
0531     const G4VSolid* n = Moved(node_); 
0532     const G4VSolid* l = Left(n); 
0533     const G4VSolid* r = Right(n); 
0534     return ( l && r ) ? num_prim_r(l) + num_prim_r(r) : 1 ;
0535 }
0536 
0537 int U4SolidTree::num_prim() const 
0538 {
0539     return num_prim_r(root); 
0540 }
0541 
0542 int U4SolidTree::num_node() const
0543 {
0544     return NumNode_r(root, 0) ; 
0545 }
0546 
0547 /**
0548 U4SolidTree::NumNode_r
0549 ------------------
0550 
0551 In order to visit all nodes it is essential to deref potentially 
0552 G4DisplacedSolid with Moved otherwise may for example miss displaced booleans. 
0553 Noticed this with the JUNO PMT with tubs subtract torus neck. 
0554 
0555 **/
0556 
0557 int U4SolidTree::NumNode_r(const G4VSolid* node_, int depth) // static 
0558 {
0559     int num = node_ ? 1 : 0 ;
0560     if( node_ )
0561     {
0562         const G4VSolid* node = Moved(node_ ); 
0563         const G4VSolid* left = Left(node); 
0564         const G4VSolid* right = Right(node); 
0565 
0566         if( left && right )
0567         {   
0568             num += NumNode_r( left,  depth+1 );  
0569             num += NumNode_r( right, depth+1 );  
0570         }   
0571 
0572         if(verbose) std::cout 
0573             << "U4SolidTree::NumNode_r " 
0574             << " depth " << std::setw(3) << depth
0575             << " type " << std::setw(20) << EntityTypeName(node)
0576             << " name " << std::setw(20) << node->GetName()
0577             << " left " << std::setw(20) << ( left ? left->GetName() : "-" ) 
0578             << " right " << std::setw(20) << ( right ? right->GetName() : "-" )
0579             << " num " << std::setw(3) << num
0580             << std::endl 
0581             ;       
0582     }
0583     return num ; 
0584 } 
0585 
0586 int U4SolidTree::num_node_select(int qcls) const
0587 {
0588     return num_node_select_r(root, qcls) ; 
0589 }
0590 int U4SolidTree::num_node_select_r(const G4VSolid* node_, int qcls) const 
0591 {
0592     int zcl = zcls(node_) ;
0593     const G4VSolid* node = Moved(node_ ); 
0594  
0595 
0596     int num = ( node_ && zcl == qcls ) ? 1 : 0 ;
0597     if( node )
0598     {
0599         if(verbose) std::cout 
0600             << "U4SolidTree::num_node_select_r"
0601             << " zcl " << zcl
0602             << " zcn " << ClassifyMaskName(zcl)
0603             << " Desc " << Desc(node)
0604             << std::endl
0605             ;  
0606 
0607         const G4VSolid* l = Left(node); 
0608         const G4VSolid* r = Right(node); 
0609         if( l && r )
0610         {   
0611             num += num_node_select_r( l, qcls );  
0612             num += num_node_select_r( r, qcls );  
0613         }   
0614     }
0615     return num ; 
0616 } 
0617 
0618 const char* U4SolidTree::desc() const 
0619 {
0620     int num_node_     = num_node(); 
0621     int num_prim_     = num_prim(); 
0622 
0623     int num_undefined = num_node_select(UNDEFINED);  
0624     int num_exclude   = num_node_select(EXCLUDE);  
0625     int num_include   = num_node_select(INCLUDE);  
0626     int num_mixed     = num_node_select(MIXED);  
0627 
0628     std::stringstream ss ; 
0629     ss 
0630        << " NODE:" << num_node_
0631        << " PRIM:" << num_prim_
0632        << " UNDEFINED:" << num_undefined
0633        << " EXCLUDE:" << num_exclude
0634        << " INCLUDE:" << num_include
0635        << " MIXED:" << num_mixed
0636        ; 
0637      
0638     std::string s = ss.str(); 
0639     return strdup(s.c_str());
0640 }
0641 
0642 void U4SolidTree::prune(bool act, int pass)
0643 {
0644     int num_include = num_node_select(INCLUDE) ;
0645 
0646     if(verbose)
0647     printf("U4SolidTree::prune act:%d pass:%d num_include %d \n", act, pass, num_include);
0648 
0649     if( num_include == 0)
0650     {
0651         if(verbose)
0652         printf("U4SolidTree::prune act:%d pass:%d find zero remaining nodes : num_include %d, will set root to nullptr \n", act, pass, num_include);
0653 
0654         if(act)
0655         {
0656             if(verbose)
0657             printf("U4SolidTree:::prune act:%d pass:%d setting root to nullptr \n", act, pass);
0658 
0659             root = nullptr ;
0660             edited = true ; 
0661         }
0662     }
0663 
0664     if(crux->size() == 0 ) return ;
0665     bool expect= crux->size() == 1 ; 
0666     if(!expect) exit(EXIT_FAILURE); 
0667     assert(expect) ;       // more than one crux node not expected
0668 
0669     G4VSolid* x = (*crux)[0] ;
0670     prune_crux(x, act, pass);
0671 }
0672 
0673 /**
0674 U4SolidTree::prune_crux
0675 -------------------
0676 
0677 **/
0678 
0679 void U4SolidTree::prune_crux(G4VSolid* x, bool act, int pass)
0680 {
0681     assert( x );
0682     bool ie = is_include_exclude(x) ;  // include left, exclude right
0683     bool ei = is_exclude_include(x) ;  // exclude left, include right 
0684     bool expect_0 =  ie ^ ei ;     // XOR definition of crux node 
0685     if(!expect_0) exit(EXIT_FAILURE) ; 
0686 
0687     G4VSolid* survivor = ie ? Left_(x) : Right_(x) ;
0688     bool expect_1 = survivor != nullptr ; 
0689     if(!expect_1) exit(EXIT_FAILURE); 
0690     assert( survivor );
0691 
0692     set_mkr(survivor, 'S') ;
0693 
0694     G4VSolid* p = parent_(x); 
0695 
0696     if(!is_include(p)) // hmm problems with this will not be apparent with the maximally unbalanced trees 
0697     {
0698         if(verbose) 
0699         printf("U4SolidTree::prune_crux act:%d pass:%d parent disqualified as not INCLUDE : handle as root of the new tree\n", act, pass) ; 
0700         p = nullptr ;    // as *p* is just local variable are changing it even when act=false 
0701     }
0702 
0703     if( p != nullptr )   // non-root prune
0704     {
0705         G4VSolid* p_left = Left_(p); 
0706         G4VSolid* p_right = Right_(p); 
0707 
0708         bool x_is_p_left  = x == p_left ; 
0709         bool x_is_p_right = x == p_right ; 
0710 
0711         set_mkr(p, 'P') ;
0712 
0713         if(x_is_p_left)
0714         {
0715             if(act)
0716             {
0717                 if(verbose)
0718                 printf("U4SolidTree:::prune_crux act:%d pass:%d SetLeft changing p.left to survivor \n", act, pass);
0719                 SetLeft(p, survivor) ; 
0720                 edited = true ; 
0721             }
0722         }
0723         else if(x_is_p_right)
0724         {
0725             if(act)
0726             {
0727                 if(verbose)
0728                 printf("U4SolidTree:prune_crux act:%d pass:%d SetRight changing p.right to survivor  \n", act, pass);
0729                 SetRight(p, survivor) ; 
0730                 edited = true ; 
0731             }
0732         }
0733     }
0734     else           // root prune
0735     {
0736         if( act )
0737         {
0738 
0739             G4String root_name = root->GetName(); 
0740             G4String survivor_name = survivor->GetName() ; 
0741             if(verbose)
0742             printf("U4SolidTree::prune_crux act:%d pass:%d changing root %s to survivor %s \n", act, pass, root_name.c_str(), survivor_name.c_str() );
0743             root = survivor ;
0744             edited = true ; 
0745         }
0746     }
0747 
0748     if(act)
0749     {
0750         instrumentTree();
0751     }
0752 }
0753 
0754 void U4SolidTree::draw(const char* msg, int pass) 
0755 {
0756     canvas->clear();
0757 
0758     //int mode = RPRE ; 
0759     int mode = IN ; 
0760     draw_r(root, mode);
0761 
0762     canvas->draw(   -1, -1, 0,0,  "zdelta" ); 
0763     canvas->draw(   -1, -1, 0,2,  "az1" ); 
0764     canvas->draw(   -1, -1, 0,3,  "az0" ); 
0765 
0766     std::cout 
0767         << msg 
0768         << " [" << std::setw(2) << pass << "]" 
0769         << " nameprefix " << ( nameprefix ? nameprefix : "-" ) 
0770         << " " << desc() 
0771         << " Order:" << OrderName(mode) 
0772         << std::endl
0773         ;
0774  
0775     for(unsigned i=0 ; i < names->size() ; i++ ) 
0776     {
0777         const std::string& name = (*names)[i] ; 
0778         std::string nam = nameprefix ? name.substr(strlen(nameprefix)) : "" ; 
0779         std::string snam = nam.substr(0,6) ;  
0780         if(false) std::cout 
0781             << std::setw(3) << i 
0782             << " : " 
0783             << std::setw(20) << name 
0784             << " : " 
0785             << std::setw(20) << nam 
0786             << " : [" 
0787             << std::setw(6) << snam << "]" 
0788             << std::endl
0789             ;
0790         canvas->draw(  i,  -1, 0,4, snam.c_str() ); 
0791     } 
0792     canvas->print(); 
0793 }
0794 
0795 void U4SolidTree::dumpNames(const char* msg) const 
0796 {
0797     std::cout << msg << std::endl ; 
0798     for(unsigned i=0 ; i < names->size() ; i++ ) 
0799     {
0800         const std::string& name = (*names)[i] ; 
0801         std::string nam = nameprefix ? name.substr(strlen(nameprefix)) : "" ; 
0802         std::string snam = nam.substr(0,6) ;  
0803         std::cout 
0804             << std::setw(3) << i 
0805             << " : " 
0806             << std::setw(35) << name 
0807             << " : " 
0808             << std::setw(35) << nam 
0809             << " : [" 
0810             << std::setw(10) << snam << "]" 
0811             << std::endl
0812             ;
0813     } 
0814 }
0815 
0816 void U4SolidTree::zdump(const char* msg) const 
0817 {
0818     std::cout << msg << std::endl ; 
0819     int mode = IN ; 
0820     zdump_r(root, mode); 
0821 }
0822 
0823 void U4SolidTree::zdump_r( const G4VSolid* node_, int mode ) const 
0824 {
0825     if( node_ == nullptr ) return ;
0826 
0827     const G4VSolid* node = Moved(node_); 
0828     zdump_r( Left(node),  mode  );
0829     zdump_r( Right(node), mode );
0830 
0831     int ix = in(node_) ;            // inorder index, aka "side", increasing from left to right 
0832     int iy = depth(node_) ;         // increasing downwards
0833     int idx = index(node_, mode);  // index for presentation 
0834 
0835     const char* tag = EntityTag(node, true) ; 
0836     int zcl = zcls(node_);                 
0837     const char* zcn = ClassifyMaskName(zcl) ; 
0838     G4String name = node->GetName(); 
0839 
0840 
0841     bool can_z = U4SolidTree::CanZ(node) ;
0842 
0843     if(can_z)
0844     {
0845         double zdelta = getZ(node_) ;  
0846         double z0, z1 ; 
0847         ZRange(z0, z1, node);  
0848 
0849         double az0 = z0 + zdelta ;  
0850         double az1 = z1 + zdelta ;  
0851 
0852         std::cout 
0853             << " ix " << std::setw(2) << ix
0854             << " iy " << std::setw(2) << iy
0855             << " idx " << std::setw(2) << idx
0856             << " tag " << std::setw(10) << tag 
0857             << " zcn " << std::setw(10) << zcn 
0858             << " zdelta " << std::setw(10) << std::fixed << std::setprecision(3) << zdelta 
0859             << " az0 " << std::setw(10) << std::fixed << std::setprecision(3) << az0
0860             << " az1 " << std::setw(10) << std::fixed << std::setprecision(3) << az1
0861             << " name " << name
0862             << std::endl 
0863             ; 
0864     }
0865 }
0866 
0867 
0868 
0869 /**
0870 U4SolidTree::draw_r
0871 ----------------
0872 
0873 Recursively paints nodes of the tree onto the canvas
0874 using the *mode* traversal order to label the nodes
0875 
0876 **/
0877 
0878 void U4SolidTree::draw_r( const G4VSolid* node_, int mode )
0879 {
0880     if( node_ == nullptr ) return ;
0881 
0882     const G4VSolid* node = Moved(node_); 
0883     draw_r( Left(node),  mode );
0884     draw_r( Right(node), mode );
0885 
0886     int ix = in(node_) ;            // inorder index, aka "side", increasing from left to right 
0887     int iy = depth(node_) ;         // increasing downwards
0888     int idx = index(node_, mode);  // index for presentation 
0889 
0890     const char* tag = EntityTag(node, true) ; 
0891     int zcl = zcls(node_);                 
0892     const char* zcn = ClassifyMaskName(zcl) ; 
0893     char mk = mkr(node); 
0894 
0895     canvas->draw(   ix, iy, 0,0,  tag); 
0896     canvas->draw(   ix, iy, 0,1,  zcn); 
0897     canvas->draw(   ix, iy, 0,2,  idx); 
0898     canvas->drawch( ix, iy, 0,3,  mk ); 
0899 
0900     bool can_y = U4SolidTree::CanY(node) ;
0901     bool can_z = U4SolidTree::CanZ(node) ;
0902     if(verbose) std::cout << "U4SolidTree::draw_r can_z " << can_z << std::endl ;  
0903 
0904     if(can_y)
0905     {
0906         double y_delta = getY(node_) ;  
0907         double y0, y1 ; 
0908         YRange(y0, y1, node);  
0909 
0910         const char* fmt = "%7.2f" ;  
0911         canvas->drawf(   ix, -1, 0,0,  y_delta   , fmt); 
0912         canvas->drawf(   ix, -1, 0,2,  y1+y_delta, fmt  ); 
0913         canvas->drawf(   ix, -1, 0,3,  y0+y_delta, fmt ); 
0914     }
0915     else if(can_z)
0916     {
0917         double z_delta = getZ(node_) ;  
0918         double z0, z1 ; 
0919         ZRange(z0, z1, node);  
0920 
0921         // below is coercing double into int val 
0922         canvas->draw(   ix, -1, 0,0,  z_delta ); 
0923         canvas->draw(   ix, -1, 0,2,  z1+z_delta ); 
0924         canvas->draw(   ix, -1, 0,3,  z0+z_delta ); 
0925     }
0926 
0927 
0928 
0929 
0930 
0931 
0932 
0933 }
0934 
0935 
0936 /**
0937 U4SolidTree::EntityTypeName
0938 -----------------------
0939 
0940 Unexpectedly G4 returns EntityType by value rather than by reference
0941 so have to strdup to avoid corruption when the G4String goes out of scope. 
0942 
0943 **/
0944 
0945 const char* U4SolidTree::EntityTypeName(const G4VSolid* solid)   // static
0946 {
0947     G4GeometryType type = solid->GetEntityType();  // G4GeometryType typedef for G4String
0948     return strdup(type.c_str()); 
0949 }
0950 
0951 const char* U4SolidTree::EntityTag(const G4VSolid* node_, bool move)   // static
0952 {
0953     const G4VSolid*  node = move ? Moved(nullptr, nullptr, node_ ) : node_ ; 
0954     return EntityTag_(node); 
0955 }
0956 
0957 
0958 const char* U4SolidTree::G4Ellipsoid_          = "Ell" ; 
0959 const char* U4SolidTree::G4Tubs_               = "Tub" ;
0960 const char* U4SolidTree::G4Polycone_           = "Pol" ;
0961 const char* U4SolidTree::G4Torus_              = "Tor" ;
0962 const char* U4SolidTree::G4Box_                = "Box" ;
0963 const char* U4SolidTree::G4Orb_                = "Orb" ;
0964 const char* U4SolidTree::G4MultiUnion_         = "MUN" ;
0965 const char* U4SolidTree::G4Sphere_             = "Sph" ;
0966 
0967 const char* U4SolidTree::G4UnionSolid_         = "Uni" ;
0968 const char* U4SolidTree::G4SubtractionSolid_   = "Sub" ;
0969 const char* U4SolidTree::G4IntersectionSolid_  = "Int" ;
0970 const char* U4SolidTree::G4DisplacedSolid_     = "Dis" ;
0971 
0972 // TODO: get this functionality from U4Solid 
0973 
0974 const char* U4SolidTree::EntityTag_( const G4VSolid* solid )
0975 {
0976     int etype = EntityType(solid); 
0977     const char* s = nullptr ;
0978     switch(etype)
0979     {
0980        case _G4Ellipsoid:         s = G4Ellipsoid_         ; break ; 
0981        case _G4Tubs:              s = G4Tubs_              ; break ; 
0982        case _G4Polycone:          s = G4Polycone_          ; break ; 
0983        case _G4Torus:             s = G4Torus_             ; break ; 
0984        case _G4Box:               s = G4Box_               ; break ; 
0985        case _G4Orb:               s = G4Orb_               ; break ; 
0986        case _G4MultiUnion:        s = G4MultiUnion_        ; break ; 
0987        case _G4Sphere:            s = G4Sphere_            ; break ; 
0988  
0989        case _G4UnionSolid:        s = G4UnionSolid_        ; break ; 
0990        case _G4SubtractionSolid:  s = G4SubtractionSolid_  ; break ; 
0991        case _G4IntersectionSolid: s = G4IntersectionSolid_ ; break ; 
0992        case _G4DisplacedSolid:    s = G4DisplacedSolid_    ; break ; 
0993     }
0994     return s ; 
0995 }
0996  
0997 
0998 const char* U4SolidTree::DirtyEntityTag_( const G4VSolid* node )
0999 {
1000     G4GeometryType type = node->GetEntityType();  // G4GeometryType typedef for G4String
1001     char* tag = strdup(type.c_str() + 2);  // +2 skip "G4"
1002     assert( strlen(tag) > 3 ); 
1003     tag[3] = '\0' ;
1004     return tag ;  
1005 }
1006  
1007 
1008 
1009 int U4SolidTree::EntityType(const G4VSolid* solid)   // static 
1010 {
1011     G4GeometryType etype = solid->GetEntityType();  // G4GeometryType typedef for G4String
1012     const char* name = etype.c_str(); 
1013 
1014     int type = _G4Other ; 
1015     if( strcmp(name, "G4Ellipsoid") == 0 )         type = _G4Ellipsoid ; 
1016     if( strcmp(name, "G4Tubs") == 0 )              type = _G4Tubs ; 
1017     if( strcmp(name, "G4Polycone") == 0 )          type = _G4Polycone ; 
1018     if( strcmp(name, "G4Torus") == 0 )             type = _G4Torus ; 
1019     if( strcmp(name, "G4Box") == 0 )               type = _G4Box ; 
1020     if( strcmp(name, "G4Orb") == 0 )               type = _G4Orb ; 
1021     if( strcmp(name, "G4MultiUnion") == 0 )        type = _G4MultiUnion ; 
1022     if( strcmp(name, "G4Sphere") == 0 )            type = _G4Sphere ; 
1023 
1024     if( strcmp(name, "G4UnionSolid") == 0 )        type = _G4UnionSolid ; 
1025     if( strcmp(name, "G4SubtractionSolid") == 0 )  type = _G4SubtractionSolid ; 
1026     if( strcmp(name, "G4IntersectionSolid") == 0 ) type = _G4IntersectionSolid ; 
1027     if( strcmp(name, "G4DisplacedSolid") == 0 )    type = _G4DisplacedSolid ; 
1028     return type ; 
1029 }
1030 
1031 std::string U4SolidTree::Desc(const G4VSolid* solid) // static
1032 {
1033     std::stringstream ss ; 
1034 
1035     ss << EntityTypeName(solid)
1036        << " name " << solid->GetName()
1037        << " bool " << IsBoolean(solid)
1038        << " disp " << IsDisplaced(solid)
1039        ; 
1040 
1041     std::string s = ss.str(); 
1042     return s ; 
1043 }
1044 
1045 bool U4SolidTree::IsBoolean(const G4VSolid* solid) // static
1046 {
1047     return dynamic_cast<const G4BooleanSolid*>(solid) != nullptr ; 
1048 }
1049 bool U4SolidTree::IsDisplaced(const G4VSolid* solid) // static
1050 {
1051     return dynamic_cast<const G4DisplacedSolid*>(solid) != nullptr ; 
1052 }
1053 const G4VSolid* U4SolidTree::Left(const G4VSolid* solid ) // static
1054 {
1055     return IsBoolean(solid) ? solid->GetConstituentSolid(0) : nullptr ; 
1056 }
1057 const G4VSolid* U4SolidTree::Right(const G4VSolid* solid ) // static
1058 {
1059     return IsBoolean(solid) ? solid->GetConstituentSolid(1) : nullptr ; 
1060 }
1061 G4VSolid* U4SolidTree::Left_(G4VSolid* solid ) // static
1062 {
1063     return IsBoolean(solid) ? solid->GetConstituentSolid(0) : nullptr ; 
1064 }
1065 G4VSolid* U4SolidTree::Right_(G4VSolid* solid ) // static
1066 {
1067     return IsBoolean(solid) ? solid->GetConstituentSolid(1) : nullptr ; 
1068 }
1069 
1070 /**
1071 U4SolidTree::Moved
1072 ---------------
1073 
1074 When node isa G4DisplacedSolid sets the rotation and translation and returns the constituentMovedSolid
1075 otherwise returns the input node.
1076 
1077 **/
1078 const G4VSolid* U4SolidTree::Moved( G4RotationMatrix* rot, G4ThreeVector* tla, const G4VSolid* node )  // static
1079 {
1080     const G4DisplacedSolid* disp = dynamic_cast<const G4DisplacedSolid*>(node) ; 
1081     if(disp)
1082     {
1083         if(rot) *rot = disp->GetFrameRotation();
1084         if(tla) *tla = disp->GetObjectTranslation();  
1085     }
1086     return disp ? disp->GetConstituentMovedSolid() : node  ;
1087 }
1088 
1089 const G4VSolid* U4SolidTree::Moved( const G4VSolid* node )  // static
1090 {
1091     const G4DisplacedSolid* disp = dynamic_cast<const G4DisplacedSolid*>(node) ; 
1092     return disp ? disp->GetConstituentMovedSolid() : node  ;
1093 }
1094 
1095 G4VSolid* U4SolidTree::Moved_( G4RotationMatrix* rot, G4ThreeVector* tla, G4VSolid* node )  // static
1096 {
1097     G4DisplacedSolid* disp = dynamic_cast<G4DisplacedSolid*>(node) ; 
1098     if(disp)
1099     {
1100         if(rot) *rot = disp->GetFrameRotation();
1101         if(tla) *tla = disp->GetObjectTranslation();  
1102     }
1103     return disp ? disp->GetConstituentMovedSolid() : node  ;
1104 }
1105 
1106 G4VSolid* U4SolidTree::Moved_( G4VSolid* node )  // static
1107 {
1108     G4DisplacedSolid* disp = dynamic_cast<G4DisplacedSolid*>(node) ; 
1109     return disp ? disp->GetConstituentMovedSolid() : node  ;
1110 }
1111 
1112 
1113 
1114 
1115 int U4SolidTree::maxdepth() const  
1116 {
1117     return Maxdepth_r( root, 0 );  
1118 }
1119 
1120 int U4SolidTree::Maxdepth_r( const G4VSolid* node_, int depth)  // static 
1121 {
1122     return IsBoolean(node_) ? std::max( Maxdepth_r(U4SolidTree::Left(node_), depth+1), Maxdepth_r(U4SolidTree::Right(node_), depth+1)) : depth ; 
1123 }
1124 
1125 /**
1126 U4SolidTree::dumpTree
1127 ------------------
1128  
1129 Postorder traversal of CSG tree
1130 **/
1131 void U4SolidTree::dumpTree(const char* msg ) const 
1132 {
1133     std::cout << msg << " maxdepth(aka height) " << height << std::endl ; 
1134     dumpTree_r(root, 0 ); 
1135 }
1136 
1137 void U4SolidTree::dumpTree_r( const G4VSolid* node_, int depth  ) const
1138 {
1139     if(IsBoolean(node_))
1140     {
1141         dumpTree_r(U4SolidTree::Left(node_) , depth+1) ; 
1142         dumpTree_r(U4SolidTree::Right(node_), depth+1) ; 
1143     }
1144 
1145     // postorder visit 
1146 
1147     G4RotationMatrix node_rot ;   // node (not tree) transforms
1148     G4ThreeVector    node_tla(0., 0., 0. ); 
1149     const G4VSolid*  node = Moved(&node_rot, &node_tla, node_ ); 
1150     assert( node ); 
1151 
1152     double zdelta_always_zero = getZ(node); 
1153 
1154     bool expect = zdelta_always_zero == 0. ; 
1155     if(!expect) exit(EXIT_FAILURE); 
1156 
1157     assert(expect ); 
1158     // Hmm thats tricky, using node arg always gives zero.
1159     // Must use node_ which is the one that might be G4DisplacedSolid. 
1160 
1161     double zdelta = getZ(node_) ;  
1162 
1163     std::cout 
1164         << " type " << std::setw(20) << EntityTypeName(node) 
1165         << " name " << std::setw(20) << node->GetName() 
1166         << " depth " << depth 
1167         << " zdelta "
1168         << std::fixed << std::setw(7) << std::setprecision(2) << zdelta
1169         << " node_tla (" 
1170         << std::fixed << std::setw(7) << std::setprecision(2) << node_tla.x() << " " 
1171         << std::fixed << std::setw(7) << std::setprecision(2) << node_tla.y() << " "
1172         << std::fixed << std::setw(7) << std::setprecision(2) << node_tla.z() << ")"
1173         ;
1174 
1175     if(U4SolidTree::CanZ(node))
1176     {
1177         double z0, z1 ; 
1178         ZRange(z0, z1, node);  
1179         std::cout 
1180             << " z1 " << std::fixed << std::setw(7) << std::setprecision(2) << z1
1181             << " z0 " << std::fixed << std::setw(7) << std::setprecision(2) << z0 
1182             << " az1 " << std::fixed << std::setw(7) << std::setprecision(2) << ( z1 + zdelta )
1183             << " az0 " << std::fixed << std::setw(7) << std::setprecision(2) << ( z0 + zdelta )
1184             ;
1185     }
1186     std::cout 
1187         << std::endl
1188         ; 
1189 }
1190 
1191 /**
1192 U4SolidTree::classifyTree
1193 ----------------------
1194 
1195 postorder traveral classifies every node doing bitwise-OR
1196 combination of the child classifications.
1197 
1198 **/
1199 
1200 int U4SolidTree::classifyTree(double zcut)  
1201 {
1202     crux->clear(); 
1203     if(verbose) std::cout << "U4SolidTree::classifyTree against zcut " << zcut  << std::endl ; 
1204     int zc = classifyTree_r(root, 0, zcut); 
1205     return zc ; 
1206 }
1207 
1208 
1209 int U4SolidTree::classifyTree_r(G4VSolid* node_, int depth, double zcut )
1210 {
1211     int zcl = 0 ; 
1212     int sid = in(node_);    // inorder 
1213     int pos = post(node_); 
1214 
1215     if(IsBoolean(node_))
1216     {
1217         int left_zcl = classifyTree_r(Left_(node_) , depth+1, zcut) ; 
1218         int right_zcl = classifyTree_r(Right_(node_), depth+1, zcut) ; 
1219 
1220         zcl |= left_zcl ; 
1221         zcl |= right_zcl ; 
1222 
1223 
1224         if(left_zcl == INCLUDE && right_zcl == EXCLUDE )
1225         {
1226             crux->push_back(node_); 
1227             set_mkr( node_, 'X' ); 
1228         }
1229         else if(left_zcl == EXCLUDE && right_zcl == INCLUDE )
1230         {
1231             crux->push_back(node_); 
1232             set_mkr( node_, 'Y' ); 
1233         }
1234 
1235         if(false) std::cout 
1236             << "U4SolidTree::classifyTree_r" 
1237             << " sid " << std::setw(2) << sid
1238             << " pos " << std::setw(2) << pos
1239             << " left_zcl " << ClassifyMaskName(left_zcl)
1240             << " right_zcl " << ClassifyMaskName(right_zcl)
1241             << " zcl " << ClassifyMaskName(zcl)
1242             << std::endl
1243             ;
1244 
1245         set_zcls( node_, zcl ); 
1246     }
1247     else
1248     {
1249         // node_ is the raw one which may be G4DisplacedSolid
1250         double zd = getZ(node_) ;  
1251         G4VSolid* node = Moved_(node_); 
1252         bool can_z = CanZ(node) ; 
1253         assert(can_z); 
1254         if(can_z)
1255         {
1256             double z0, z1 ; 
1257             ZRange(z0, z1, node);  
1258 
1259             double az0 =  z0 + zd ; 
1260             double az1 =  z1 + zd ; 
1261 
1262             zcl = ClassifyZCut( az0, az1, zcut ); 
1263 
1264             if(verbose) std::cout 
1265                 << "U4SolidTree::classifyTree_r"
1266                 << " sid " << std::setw(2) << sid
1267                 << " pos " << std::setw(2) << pos
1268                 << " zd " << std::fixed << std::setw(7) << std::setprecision(2) << zd
1269                 << " z1 " << std::fixed << std::setw(7) << std::setprecision(2) << z1
1270                 << " z0 " << std::fixed << std::setw(7) << std::setprecision(2) << z0 
1271                 << " az1 " << std::fixed << std::setw(7) << std::setprecision(2) << az1
1272                 << " az0 " << std::fixed << std::setw(7) << std::setprecision(2) << az0 
1273                 << " zcl " << ClassifyName(zcl)
1274                 << std::endl
1275                 ;
1276 
1277             set_zcls( node_, zcl );   // HMM node_ or node ?
1278         }
1279     }
1280     return zcl ; 
1281 }
1282 
1283 int U4SolidTree::classifyMask(const G4VSolid* top) const   // NOT USED ?
1284 {
1285     return classifyMask_r(top, 0); 
1286 }
1287 int U4SolidTree::classifyMask_r( const G4VSolid* node_, int depth ) const 
1288 {
1289     int mask = 0 ; 
1290     if(IsBoolean(node_))
1291     {
1292         mask |= classifyMask_r( Left(node_) , depth+1 ) ; 
1293         mask |= classifyMask_r( Right(node_), depth+1 ) ; 
1294     }
1295     else
1296     {
1297         mask |= zcls(node_) ; 
1298     }
1299     return mask ; 
1300 }
1301 
1302 /**
1303 U4SolidTree::ClassifyZCut
1304 ------------------------
1305 
1306 Inclusion status of solid with regard to a particular zcut::
1307 
1308                        --- 
1309                         .
1310                         .   EXCLUDE  : zcut entirely above the solid
1311                         .
1312                         .
1313       +---zd+z1----+   --- 
1314       |            |    .   
1315       | . zd . . . |    .   STRADDLE : zcut within z range of solid
1316       |            |    .
1317       +---zd+z0 ---+   ---
1318                         .
1319                         .   INCLUDE  : zcut fully below the solid 
1320                         .
1321                         .
1322                        ---  
1323 
1324 **/
1325 
1326 int U4SolidTree::ClassifyZCut( double az0, double az1, double zcut ) // static
1327 {
1328     assert( az1 > az0 ); 
1329     int cls = UNDEFINED ; 
1330     if(       zcut <= az0 )              cls = INCLUDE ; 
1331     else if ( zcut < az1 && zcut > az0 ) cls = STRADDLE ; 
1332     else if ( zcut >= az1              ) cls = EXCLUDE ; 
1333     return cls ; 
1334 }
1335 
1336 const char* U4SolidTree::UNDEFINED_ = "UNDEFINED" ; 
1337 const char* U4SolidTree::INCLUDE_   = "INCLUDE" ; 
1338 const char* U4SolidTree::STRADDLE_  = "STRADDLE" ; 
1339 const char* U4SolidTree::EXCLUDE_   = "EXCLUDE" ; 
1340 
1341 const char* U4SolidTree::ClassifyName( int zcl ) // static 
1342 {
1343     const char* s = nullptr ; 
1344     switch( zcl )
1345     {
1346         case UNDEFINED: s = UNDEFINED_ ; break ; 
1347         case INCLUDE  : s = INCLUDE_ ; break ; 
1348         case STRADDLE : s = STRADDLE_ ; break ; 
1349         case EXCLUDE  : s = EXCLUDE_ ; break ; 
1350     }
1351     return s ; 
1352 }
1353 
1354 const char* U4SolidTree::ClassifyMaskName( int zcl ) // static
1355 {
1356     std::stringstream ss ; 
1357 
1358     if( zcl == UNDEFINED ) ss << UNDEFINED_[0] ; 
1359     if( zcl & INCLUDE )  ss << INCLUDE_[0] ; 
1360     if( zcl & STRADDLE ) ss << STRADDLE_[0] ; 
1361     if( zcl & EXCLUDE )  ss << EXCLUDE_[0] ; 
1362     std::string s = ss.str(); 
1363     return strdup(s.c_str()); 
1364 }
1365 
1366 /**
1367 U4SolidTree::apply_cut
1368 ------------------
1369 
1370 See opticks/sysrap/tests/TreePruneTest.cc and opticks/sysrap/tests/tree.sh
1371 for development of tree cutting and pruning.
1372 
1373 Steps:
1374 
1375 1. classify the nodes of the tree against the zcut 
1376 2. change STRADDLE nodes params and transforms according to the zcut
1377 3. reclassify against the zcut, so STRADDLE nodes should become INCLUDE nodes
1378 4. edit the tree to remove the EXCLUDE nodes
1379 
1380 **/
1381 
1382 void U4SolidTree::apply_cut(double zcut)
1383 {
1384     if(verbose)
1385     printf("U4SolidTree::apply_cut %7.2f \n", zcut );
1386 
1387     if(verbose)
1388     Draw(root, "U4SolidTree::apply_cut root before cut"); 
1389 
1390     unsigned pass = 0 ;
1391     unsigned maxpass = 10 ;
1392 
1393     while( root != nullptr && zcls(root) != INCLUDE && pass < maxpass )
1394     {
1395         classifyTree(zcut);   // set n.cls n.mkr
1396 
1397         cutTree_r(root, 0, zcut); 
1398 
1399         classifyTree(zcut);   // set n.cls n.mkr
1400 
1401         prune(false, pass);
1402 
1403         if(verbose)
1404         draw("U4SolidTree::apply_cut before prune", pass );
1405 
1406         prune(true, pass);
1407 
1408         classifyTree(zcut);
1409         instrumentTree();
1410 
1411         if(verbose)
1412         draw("tree::apply_cut after prune and re-classify", pass );
1413 
1414         pass++ ;
1415     }
1416 }
1417 
1418 void U4SolidTree::collectNames_inorder_r( const G4VSolid* n_, int depth )
1419 {
1420     if(n_ == nullptr) return ; 
1421     const G4VSolid* n = Moved(nullptr, nullptr, n_) ;  
1422     collectNames_inorder_r(Left(n),  depth+1); 
1423 
1424     G4String name = n->GetName(); 
1425     names->push_back(name);  
1426 
1427     collectNames_inorder_r(Right(n), depth+1); 
1428 }
1429 
1430 void U4SolidTree::cutTree_r( const G4VSolid* node_, int depth, double zcut )
1431 {
1432     if(IsBoolean(node_))
1433     {
1434         cutTree_r( U4SolidTree::Left(node_) , depth+1, zcut ) ; 
1435         cutTree_r( U4SolidTree::Right(node_), depth+1, zcut ) ; 
1436     }
1437     else
1438     {
1439         // get tree frame zcut into local frame of the node 
1440         double zdelta = getZ(node_) ; 
1441         double local_zcut = zcut - zdelta ; 
1442 
1443         int zcl = zcls(node_);   
1444         if( zcl == STRADDLE )
1445         {
1446             if(verbose) std::cout 
1447                 << "U4SolidTree::cutTree_r"
1448                 << " depth " << std::setw(2) << depth
1449                 << " zdelta " << std::fixed << std::setw(10) << std::setprecision(4) << zdelta 
1450                 << " zcut " << std::fixed << std::setw(10) << std::setprecision(4) << zcut
1451                 << " local_zcut " << std::fixed << std::setw(10) << std::setprecision(4) << local_zcut
1452                 << std::endl
1453                 ;
1454 
1455             ApplyZCut( const_cast<G4VSolid*>(node_), local_zcut ); 
1456         } 
1457     }
1458 }
1459 
1460 void U4SolidTree::collectNodes( std::vector<const G4VSolid*>& nodes, const G4VSolid* top, int query_zcls  )
1461 {
1462     collectNodes_r(nodes, top, 0, query_zcls);  
1463 }
1464 
1465 void U4SolidTree::collectNodes_r( std::vector<const G4VSolid*>& nodes, const G4VSolid* node_, int query_zcls, int depth  )
1466 {
1467     if(IsBoolean(node_))
1468     {
1469         collectNodes_r( nodes, Left(node_) , query_zcls, depth+1 ) ; 
1470         collectNodes_r( nodes, Right(node_), query_zcls, depth+1 ) ; 
1471     }
1472     else
1473     {
1474         int zcl = zcls(node_) ; 
1475         if( zcl == query_zcls )
1476         {
1477             nodes.push_back(node_) ; // node_ or node ?
1478         } 
1479     }
1480 } 
1481 
1482 void U4SolidTree::ApplyZCut( G4VSolid* node_, double local_zcut ) // static
1483 {
1484     G4VSolid* node = Moved_(node_ ); 
1485     if(verbose) std::cout << "U4SolidTree::ApplyZCut " << EntityTypeName(node) << std::endl ; 
1486     switch(EntityType(node))
1487     {
1488         case _G4Ellipsoid: ApplyZCut_G4Ellipsoid( node  , local_zcut);  break ; 
1489         case _G4Tubs:      ApplyZCut_G4Tubs(      node_ , local_zcut);  break ; // cutting tubs requires changing transform, hence node_
1490         case _G4Polycone:  ApplyZCut_G4Polycone(  node  , local_zcut);  break ; 
1491         default: 
1492         { 
1493             std::cout 
1494                 << "U4SolidTree::ApplyZCut FATAL : not implemented for entityType " 
1495                 << EntityTypeName(node) 
1496                 << std::endl 
1497                 ; 
1498             assert(0) ; 
1499         } ;
1500     }
1501 }
1502 
1503 /**
1504 U4SolidTree::ApplyZCut_G4Ellipsoid
1505 --------------------------------
1506      
1507 ::
1508 
1509     local                                             absolute 
1510     frame                                             frame    
1511 
1512     z1  +-----------------------------------------+    zd   
1513          \                                       /
1514            
1515             .                              .
1516     _________________________________________________ zcut 
1517                 .                     .
1518                                  
1519     z0                 .      .   
1520                          
1521                                      
1522 **/
1523 
1524 
1525 void U4SolidTree::ApplyZCut_G4Ellipsoid( G4VSolid* node, double local_zcut)
1526 {  
1527     G4Ellipsoid* ellipsoid =  dynamic_cast<G4Ellipsoid*>(node) ;  
1528     assert(ellipsoid); 
1529 
1530     double z0 = ellipsoid->GetZBottomCut() ; 
1531     double z1 = ellipsoid->GetZTopCut() ;
1532     
1533     double new_z0 = local_zcut ; 
1534 
1535     bool expect = new_z0 >= z0 && new_z0 < z1 && z1 > z0 ; 
1536     if(!expect) exit(EXIT_FAILURE); 
1537     assert( expect ); 
1538 
1539     double new_z1 = z1 ; 
1540 
1541     ellipsoid->SetZCuts( new_z0, new_z1 ); 
1542 }
1543 
1544 /**
1545 U4SolidTree::ApplyZCut_G4Polycone_NotWorking
1546 -----------------------------------------
1547 
1548 Currently limited to only 2 Z-planes, 
1549 to support more that 2 would need to delve 
1550 into the r-z details which should be straightforward, 
1551 just it is not yet implemented. 
1552 
1553 SMOKING GUN : CHANGING OriginalParameters looks like
1554 its working when just look at Opticks results but the 
1555 change does not get thru to Geant4 
1556 
1557 
1558 **/
1559 
1560 void U4SolidTree::ApplyZCut_G4Polycone_NotWorking( G4VSolid* node, double local_zcut)
1561 {  
1562     G4Polycone* polycone = dynamic_cast<G4Polycone*>(node) ;  
1563     assert(polycone); 
1564     G4PolyconeHistorical* pars = polycone->GetOriginalParameters(); 
1565 
1566     unsigned num_z = pars->Num_z_planes ; 
1567     for(unsigned i=1 ; i < num_z ; i++)
1568     {
1569         double z0 = pars->Z_values[i-1] ; 
1570         double z1 = pars->Z_values[i] ; 
1571         bool expect = z1 > z0 ; 
1572         assert(expect);   
1573         if(!expect) exit(EXIT_FAILURE); 
1574     }
1575 
1576     assert( num_z == 2 );    // simplifying assumption 
1577     pars->Z_values[0] = local_zcut ;   // get into polycone local frame
1578 
1579     polycone->SetOriginalParameters(pars); // OOPS : SEEMS Geant4 IGNORES
1580 }
1581 
1582 
1583 /**
1584 U4SolidTree::ApplyZCut_G4Polycone
1585 ------------------------------
1586 
1587 Use placement new to replace the polycone using the ctor again 
1588 with different params so Geant4 cannot ignore 
1589 
1590 Note that the radii is modified according to the cut, in 
1591 order to cut without changing shape other than the cut.::
1592 
1593                                      
1594                                       .  (z1,r1)
1595                                .
1596                          .   
1597                    .
1598              .      (zc,rc)
1599        +    
1600      (z0,r0)
1601 
1602 **/
1603 
1604 void U4SolidTree::ApplyZCut_G4Polycone( G4VSolid* node, double local_zcut)
1605 {  
1606     G4Polycone* polycone = dynamic_cast<G4Polycone*>(node) ;  
1607     assert(polycone); 
1608 
1609     G4PolyconeHistorical* pars = polycone->GetOriginalParameters(); 
1610 
1611     bool expect ; 
1612   
1613     unsigned num_z = pars->Num_z_planes ; 
1614     G4double* zp = new G4double[num_z] ; 
1615     G4double* ri = new G4double[num_z] ; 
1616     G4double* ro = new G4double[num_z] ; 
1617 
1618     for(unsigned i=0 ; i < num_z ; i++)
1619     {
1620         zp[i] = pars->Z_values[i] ; 
1621         ri[i] = pars->Rmin[i]; 
1622         ro[i] = pars->Rmax[i]; 
1623     }
1624 
1625     for(unsigned i=1 ; i < num_z ; i++)
1626     {
1627         double z0 = zp[i-1] ; 
1628         double z1 = zp[i] ; 
1629         expect = z1 > z0 ; 
1630         assert(expect);   
1631         if(!expect) exit(EXIT_FAILURE); 
1632     }
1633 
1634     assert( num_z == 2 );  // simplifying assumption 
1635 
1636     double zfrac = (local_zcut - zp[0])/(zp[1] - zp[0]) ;
1637  
1638     double ri_zcut = ri[0] + zfrac*( ri[1] - ri[0] ) ;   
1639     double ro_zcut = ro[0] + zfrac*( ro[1] - ro[0] ) ; 
1640 
1641     zp[0] = local_zcut ;   
1642     ri[0] = ri_zcut ; 
1643     ro[0] = ro_zcut ; 
1644    
1645 
1646     G4String name = polycone->GetName() ; 
1647     G4double startPhi = polycone->GetStartPhi() ;  
1648     G4double endPhi = polycone->GetEndPhi() ; 
1649 
1650     G4SolidStore::GetInstance()->DeRegister(polycone);  // avoids G4SolidStore segv at cleanup
1651 
1652     // placement new 
1653     G4Polycone* polycone_replacement = new (polycone) G4Polycone( name, startPhi, endPhi, num_z, zp, ri, ro ); 
1654     expect = polycone_replacement == polycone ;  
1655     assert( expect ); 
1656     if(!expect) exit(EXIT_FAILURE); 
1657 
1658 }
1659 
1660 
1661 
1662 
1663 /**
1664 U4SolidTree::ApplyZCut_G4Tubs
1665 --------------------------
1666 
1667 * SEE THE TAIL OF THIS FILE FOR A DERIVATION OF new_hz AND zoffset 
1668 
1669 G4Tubs is more difficult to cut than G4Polycone because it is 
1670 symmetrically defined, so cutting requires shifting too. 
1671 Replacing all G4Tubs with G4Polycone at the initial clone stage
1672 avoids using this method at all. 
1673 
1674 BUT if you prefer not to promote all G4Tubs to G4Polycone 
1675 it is also possible to cut G4Tubs by changing the 
1676 G4DisplacedSolid transforms appropriately. 
1677 
1678 Initially tried using G4DisplacedSolid::SetObjectTranslation
1679 but looking at the implementation that only changes one
1680 of the two transforms held by the object. Hence used 
1681 placement new yet again to rerun the ctor, after deregistering 
1682 the disp solid from G4SolidStore.
1683 
1684 NOTE THAT THE G4Tubs MUST HAVE AN ASSOCIATED TRANSFORM (EVEN IDENTITY MATRIX)
1685 FOR THIS TO WORK : SO ALWAYS USE THE G4BooleanSolid ctor with rot and tla
1686   
1687 **/
1688 
1689 void U4SolidTree::ApplyZCut_G4Tubs( G4VSolid* node_ , double local_zcut )
1690 { 
1691     if( PROMOTE_TUBS_TO_POLYCONE )
1692         assert(0 && "All G4Tubs should have been promoted to G4Polycone at clone stage, how did you get here ?") ; 
1693 
1694     bool expect ; 
1695     G4RotationMatrix node_rot ; 
1696     G4ThreeVector    node_tla(0., 0., 0. ); 
1697     G4VSolid*  node = Moved_(&node_rot, &node_tla, node_ ); 
1698 
1699     G4Tubs* tubs = dynamic_cast<G4Tubs*>(node) ;  
1700     expect = tubs != nullptr ; 
1701     assert(expect); 
1702     if(!expect) exit(EXIT_FAILURE); 
1703 
1704     G4DisplacedSolid* disp = dynamic_cast<G4DisplacedSolid*>(node_) ; 
1705     expect = disp != nullptr ; // transform must be associated as must change offset to cut G4Tubs
1706     assert( expect ); 
1707     if(!expect) exit(EXIT_FAILURE); 
1708 
1709     double hz = tubs->GetZHalfLength() ; 
1710     double new_hz  = (hz - local_zcut)/2. ;  
1711     double zoffset = (hz + local_zcut)/2. ; 
1712 
1713     tubs->SetZHalfLength(new_hz);  
1714     node_tla.setZ( node_tla.z() + zoffset ); 
1715 
1716     std::cout 
1717         << "U4SolidTree::ApplyZCut_G4Tubs"
1718         << " hz " << hz
1719         << " local_zcut " << local_zcut
1720         << " new_hz " << new_hz
1721         << " zoffset " << zoffset 
1722         << std::endl 
1723         ;
1724 
1725     G4String disp_name = disp->GetName() ; 
1726     //disp->CleanTransformations() ; // needed ? perhaps avoids small leak 
1727     G4SolidStore::GetInstance()->DeRegister(disp);  // avoids G4SolidStore segv at cleanup
1728 
1729     G4DisplacedSolid* disp2 = new (disp) G4DisplacedSolid(disp_name, node, &node_rot, node_tla );  
1730 
1731     expect =  disp2 == disp ; 
1732     assert(expect); 
1733     if(!expect) exit(EXIT_FAILURE); 
1734 }
1735 
1736 
1737 
1738 
1739 /**
1740 U4SolidTree::dumpUp
1741 ----------------
1742 
1743 Ordinary postorder recursive traverse in order to get to all nodes. 
1744 This approach should allow to obtain combination transforms in
1745 complex trees. 
1746 
1747 **/
1748 
1749 void U4SolidTree::dumpUp(const char* msg) const 
1750 {
1751     assert( parent_map ); 
1752     std::cout << msg << std::endl ; 
1753     dumpUp_r(root, 0); 
1754 }
1755 
1756 void U4SolidTree::dumpUp_r(const G4VSolid* node, int depth) const  
1757 {
1758     if(IsBoolean(node))
1759     {
1760         dumpUp_r(U4SolidTree::Left(  node ), depth+1 ); 
1761         dumpUp_r(U4SolidTree::Right( node ), depth+1 ); 
1762     }
1763     else
1764     {
1765         G4RotationMatrix* tree_rot = nullptr ; 
1766         G4ThreeVector    tree_tla(0., 0., 0. ); 
1767         getTreeTransform(tree_rot, &tree_tla, node ); 
1768 
1769         const G4VSolid* nd = Moved(nullptr, nullptr, node); 
1770         std::cout 
1771             << "U4SolidTree::dumpUp_r" 
1772             << " depth " << depth 
1773             << " type " << std::setw(20) << EntityTypeName(nd) 
1774             << " name " << std::setw(20) << nd->GetName() 
1775             << " tree_tla (" 
1776             << std::fixed << std::setw(7) << std::setprecision(2) << tree_tla.x() << " " 
1777             << std::fixed << std::setw(7) << std::setprecision(2) << tree_tla.y() << " "
1778             << std::fixed << std::setw(7) << std::setprecision(2) << tree_tla.z() 
1779             << ")"
1780             << std::endl
1781             ; 
1782    }
1783 }
1784 
1785 /**
1786 U4SolidTree::getTreeTransform
1787 -------------------------------
1788 
1789 Would normally use parent links to determine all transforms relevant to a node, 
1790 but Geant4 boolean trees do not have parent links. 
1791 Hence use an external parent_map to provide uplinks enabling iteration 
1792 up the tree from any node up to the root. 
1793 
1794 **/
1795 
1796 void U4SolidTree::getTreeTransform( G4RotationMatrix* rot, G4ThreeVector* tla, const G4VSolid* node ) const 
1797 {
1798     bool expect = rot == nullptr ; 
1799     assert(expect && "non null rotation not implemented"); 
1800     if(!expect) exit(EXIT_FAILURE); 
1801 
1802     const G4VSolid* nd = node ; 
1803 
1804     unsigned count = 0 ; 
1805     while(nd)
1806     {
1807         G4RotationMatrix r ; 
1808         G4ThreeVector    t(0., 0., 0. ); 
1809         const G4VSolid* dn = Moved( &r, &t, nd ); 
1810         assert( r.isIdentity() );  // simplifying assumption 
1811 
1812         *tla += t ;    // add up the translations 
1813 
1814         if(false) std::cout
1815             << "U4SolidTree::getTreeTransform" 
1816             << " count " << std::setw(2) << count 
1817             << " dn.name " << std::setw(20) << dn->GetName()
1818             << " dn.type " << std::setw(20) << dn->GetEntityType()
1819             << " dn.t (" 
1820             << std::fixed << std::setw(7) << std::setprecision(2) << t.x() << " " 
1821             << std::fixed << std::setw(7) << std::setprecision(2) << t.y() << " "
1822             << std::fixed << std::setw(7) << std::setprecision(2) << t.z() 
1823             << ")"
1824             << " tla (" 
1825             << std::fixed << std::setw(7) << std::setprecision(2) << tla->x() << " " 
1826             << std::fixed << std::setw(7) << std::setprecision(2) << tla->y() << " "
1827             << std::fixed << std::setw(7) << std::setprecision(2) << tla->z() 
1828             << ")"
1829             << std::endl 
1830             ; 
1831 
1832         nd = (*parent_map)[nd] ; // parentmap lineage uses G4DisplacedSolid so not using *dn* here
1833         count += 1 ; 
1834     }     
1835 }
1836 
1837 /**
1838 U4SolidTree::DeepClone
1839 -------------------
1840 
1841 Clones a CSG tree of solids, assuming that the tree is
1842 composed only of the limited set of primitives that are supported. 
1843 
1844 G4BooleanSolid copy ctor just steals constituent pointers so 
1845 it does not make an independent copy.  
1846 Unlike the primitive copy ctors (at least those looked at: G4Polycone, G4Tubs) 
1847 which appear to make properly independent copies 
1848 
1849 **/
1850 
1851 G4VSolid* U4SolidTree::DeepClone( const  G4VSolid* solid )  // static 
1852 {
1853     LOG(LEVEL) << "[" ; 
1854     G4RotationMatrix* rot = nullptr ; 
1855     G4ThreeVector* tla = nullptr ; 
1856     int depth = 0 ; 
1857     G4VSolid* clone = DeepClone_r(solid, depth, rot, tla );  
1858     LOG(LEVEL) << "]" ; 
1859     return clone ; 
1860 }
1861 
1862 /**
1863 U4SolidTree::DeepClone_r
1864 --------------------------
1865 
1866 G4DisplacedSolid is a wrapper for the right hand side boolean constituent 
1867 which serves the purpose of holding the transform. The G4DisplacedSolid 
1868 is automatically created by the G4BooleanSolid ctor when there is an associated transform.  
1869 
1870 The below *rot* and *tla* look at first glance like they are not used. 
1871 But look more closely, the recursive DeepClone_r calls within BooleanClone are using them 
1872 across the generations. **This structure is necessary for BooleanClone because the 
1873 transform from the child is needed when cloning the parent**.
1874 
1875 The transform is wrapped around the right hand child (either leaf or operator) 
1876 but for the clone it is needed at construction of the parent.  
1877 
1878 **/
1879 
1880 G4VSolid* U4SolidTree::DeepClone_r( const G4VSolid* node_, int depth, G4RotationMatrix* rot, G4ThreeVector* tla )  // static 
1881 {
1882     const G4VSolid* node = Moved( rot, tla, node_ );  // <-- checking this node for a wrapper, not expecting one for the root node
1883     DeepCloneDump(depth, node_, node, tla ); 
1884     G4VSolid* clone = IsBoolean(node) ? BooleanClone(node, depth, rot, tla ) : PrimitiveClone(node) ; 
1885     ExpectNonNull( clone, "U4SolidTree::DeepClone_r GOT null clone " ); 
1886     return clone ; 
1887 }    
1888 
1889 /**
1890 U4SolidTree::BooleanClone
1891 -----------------------
1892 
1893 The left and right G4VSolid outputs from DeepClone_r will not be G4DisplacedSolid because
1894 those get "dereferenced" by Moved and the rot/tla set.  This means that the information 
1895 from the G4DisplacedSolid is available. This approach is necessary as the G4DisplacedSolid
1896 is an "internal" object that the G4BooleanSolid ctor creates from the rot and tla ctor arguments. 
1897 
1898 HMM : rot and tla arguments were not used and they are not always null ... 
1899 that suggests there is a lack of generality here for booleans of booleans etc... 
1900 with translations that apply on top of translations. 
1901 However for JUNO PMTs are not expecting any rotations
1902 and are not expecting more than one level of translation. 
1903 
1904 **/
1905 
1906 G4VSolid* U4SolidTree::BooleanClone( const  G4VSolid* solid, int depth, G4RotationMatrix* rot, G4ThreeVector* tla ) // static
1907 {
1908     if(verbose) std::cout << "U4SolidTree::BooleanClone" << std::endl ; 
1909 
1910     ExpectNoRotation(    rot, "U4SolidTree::BooleanClone expect_rot ERROR " ); 
1911     ExpectNoTranslation( tla, "U4SolidTree::BooleanClone expect_tla ERROR (not expecting more than one level of translation) ", true ); 
1912 
1913     G4String name = solid->GetName() ; 
1914     G4RotationMatrix lrot, rrot ;  
1915     G4ThreeVector    ltra, rtra ; 
1916 
1917     const G4BooleanSolid* src_boolean = dynamic_cast<const G4BooleanSolid*>(solid) ; 
1918     G4VSolid* left  = DeepClone_r( src_boolean->GetConstituentSolid(0), depth+1, &lrot, &ltra ) ; 
1919     G4VSolid* right = DeepClone_r( src_boolean->GetConstituentSolid(1), depth+1, &rrot, &rtra ) ; 
1920 
1921     ExpectNoDisplaced(     left, "U4SolidTree::BooleanClone expect_left ERROR "); 
1922     ExpectNoRotation(     &lrot, "U4SolidTree::BooleanClone expect_lrot ERROR : G4 should never have left transforms  " ); 
1923     ExpectNoTranslation(  &ltra, "U4SolidTree::BooleanClone expect_ltra ERROR : not expectsing left translation", false ); 
1924     ExpectNoRotation(     &rrot, "U4SolidTree::BooleanClone expect_rrot ERROR : simplifying assumptions  " ); 
1925 
1926     G4VSolid* clone = nullptr ; 
1927     switch(EntityType(solid))
1928     {
1929         case _G4UnionSolid        : clone = new G4UnionSolid(       name, left, right, &rrot, rtra ) ; break ; 
1930         case _G4SubtractionSolid  : clone = new G4SubtractionSolid( name, left, right, &rrot, rtra ) ; break ;
1931         case _G4IntersectionSolid : clone = new G4IntersectionSolid(name, left, right, &rrot, rtra ) ; break ; 
1932     } 
1933     CheckBooleanClone( clone, left, right ); 
1934     return clone ; 
1935 }
1936 
1937 
1938 
1939 
1940 
1941 void U4SolidTree::DeepCloneDump( int depth, const G4VSolid* node_ , const G4VSolid* node , G4ThreeVector* tla )
1942 {
1943     LOG(LEVEL) 
1944         << "[ depth " << depth 
1945         << " EntityTypeName(node_) " << EntityTypeName(node_)  
1946         << " EntityTypeName(node) "  << EntityTypeName(node)  
1947         ; 
1948      
1949     // if node_ isa G4DisplacedSolid node will not be and rot/tla will be set 
1950 
1951     if(verbose) 
1952     {
1953         std::cout 
1954             << "U4SolidTree::DeepClone_r(preorder visit)"
1955             << " type " << std::setw(20) << EntityTypeName(node)
1956             << " name " << std::setw(20) << node->GetName()
1957             << " depth " << std::setw(2) << depth
1958             ; 
1959         if(tla) std::cout 
1960             << " tla (" 
1961             << tla->x() 
1962             << " " 
1963             << tla->y() 
1964             << " " 
1965             << tla->z() 
1966             << ")"  
1967             ; 
1968         std::cout 
1969             << std::endl 
1970             ; 
1971     }
1972 } 
1973 
1974 void U4SolidTree::ExpectNoRotation( const G4RotationMatrix* rot, const char* msg )
1975 {
1976     bool expect_rot = rot == nullptr || rot->isIdentity() ;   
1977     if(!expect_rot) std::cout << msg  << std::endl ; 
1978     assert( expect_rot ); 
1979     if(!expect_rot) exit(EXIT_FAILURE); 
1980 }
1981 
1982 void U4SolidTree::ExpectNoTranslation( const G4ThreeVector* tla, const char* msg, bool skip_assert )
1983 {
1984     if(!tla) return ; 
1985     bool expect_tla = tla->x() == 0. && tla->y() == 0. && tla->z() == 0. ;
1986     if(!expect_tla) 
1987     {
1988         std::cout 
1989             << msg 
1990             << " tla( " 
1991             << tla->x() 
1992             << " " 
1993             << tla->y() 
1994             << " " 
1995             << tla->z() 
1996             << ") " 
1997             << std::endl
1998             ; 
1999 
2000         if( skip_assert )
2001         {
2002             LOG(fatal) << " TMP SKIP ASSERT " ; 
2003         }
2004         else
2005         {
2006             assert( expect_tla );
2007             if(!expect_tla) exit(EXIT_FAILURE); 
2008         }
2009     }
2010 }
2011 
2012 
2013 void U4SolidTree::ExpectNoDisplaced( const G4VSolid* solid, const char* msg )
2014 {
2015     bool expect = dynamic_cast<const G4DisplacedSolid*>(solid) == nullptr ; 
2016     if(!expect) std::cout << msg << std::endl ; 
2017     if(!expect) exit(EXIT_FAILURE); 
2018     assert( expect );
2019 }
2020 
2021 void U4SolidTree::ExpectNonNull( const G4VSolid* clone, const char* msg )
2022 {
2023     bool expect = clone != nullptr ; 
2024     if(!expect) std::cout << msg << std::endl ; 
2025 
2026     assert(expect);
2027     if(!expect) exit(EXIT_FAILURE); 
2028 }
2029 
2030 
2031 
2032 
2033 
2034 void U4SolidTree::CheckBooleanClone( const G4VSolid* clone, const G4VSolid* left, const G4VSolid* right ) // static
2035 {
2036     bool expect ; 
2037 
2038     if(!clone) std::cout << "U4SolidTree::CheckBooleanClone FATAL " << std::endl ; 
2039 
2040     expect = clone != nullptr ; 
2041     assert(expect);
2042     if(!expect) exit(EXIT_FAILURE); 
2043  
2044     const G4BooleanSolid* boolean = dynamic_cast<const G4BooleanSolid*>(clone) ; 
2045 
2046     // lhs is never wrapped in G4DisplacedSolid 
2047     const G4VSolid* lhs = boolean->GetConstituentSolid(0) ; 
2048     const G4DisplacedSolid* lhs_disp = dynamic_cast<const G4DisplacedSolid*>(lhs) ; 
2049     expect = lhs_disp == nullptr && lhs == left ; 
2050     assert(expect) ;      
2051     if(!expect) exit(EXIT_FAILURE); 
2052 
2053     // rhs will be wrapped in G4DisplacedSolid as above G4BooleanSolid ctor has transform rrot/rtla
2054     const G4VSolid* rhs = boolean->GetConstituentSolid(1) ; 
2055     const G4DisplacedSolid* rhs_disp = dynamic_cast<const G4DisplacedSolid*>(rhs) ; 
2056 
2057     expect = rhs_disp != nullptr && rhs != right ;  
2058     assert(expect);    
2059     if(!expect) exit(EXIT_FAILURE); 
2060 
2061     const G4VSolid* right_check = rhs_disp->GetConstituentMovedSolid() ;
2062     expect = right_check == right ; 
2063 
2064     assert(expect);    
2065     if(!expect) exit(EXIT_FAILURE); 
2066 }
2067 
2068 void U4SolidTree::GetBooleanBytes(char** bytes, int& num_bytes, const G4VSolid* solid ) // static
2069 {
2070     int type = EntityType(solid); 
2071     switch(type)
2072     {
2073         case _G4UnionSolid        : num_bytes = sizeof(G4UnionSolid)        ; break ; 
2074         case _G4SubtractionSolid  : num_bytes = sizeof(G4SubtractionSolid)  ; break ; 
2075         case _G4IntersectionSolid : num_bytes = sizeof(G4IntersectionSolid) ; break ; 
2076     } 
2077 
2078     *bytes = new char[num_bytes]; 
2079     memcpy( *bytes, (char*) solid, num_bytes ); 
2080 }
2081 
2082 int U4SolidTree::CompareBytes( char* bytes0, char* bytes1, int num_bytes ) // static
2083 {
2084     int mismatch = 0 ; 
2085     for(int i=0 ; i < num_bytes ; i++ ) 
2086     {   
2087         if( bytes0[i] != bytes1[i] ) 
2088         {   
2089             printf("mismatch %5d : %3d : %3d  \n", i, int(bytes0[i]), int(bytes1[i]) ); 
2090             mismatch++ ;   
2091         }   
2092     }   
2093     printf("mismatch %d\n", mismatch); 
2094     return mismatch ; 
2095 }
2096 
2097 /**
2098 U4SolidTree::PlacementNewDupe
2099 --------------------------
2100 
2101 Technical test that should do nothing and in the process
2102 demonstrate that can use placement new to replace a boolean solid 
2103 object with a bit perfect duplicate in the same memory address. 
2104 
2105 1. access and hold innards of the solid in local scope
2106 2. de-registers the solid from G4SolidStore
2107 3. placement new instantiate replacement object at the same memory address 
2108    putting the pieces back together again 
2109 
2110 The motivation for this is that tree pruning needs to be able 
2111 to SetRight/SetLeft but G4BooleanSolid has no such methods. 
2112 So in order to workaround this can use placement new to construct 
2113 a duplicate object or one with some inputs changed.
2114 
2115 G4VSolid registers/deregisters with G4SolidStore in ctor/dtor.
2116 This means that using placement new to replace a solid causes memory 
2117 corruption with::
2118 
2119     malloc: *** error for object 0x7f888c5084b0: pointer being freed was not allocated
2120 
2121 Avoided this problem by manually deregistering from G4SolidStore
2122 Although G4VSolid dtor deregisters it would be wrong to delete 
2123 as do not want to deallocate the memory, want to reuse it.
2124 
2125 When using the no transform ctor the dupe is a bit perfect match. 
2126 With the transform ctor see ~2 bytes different from fPtrSolidB.
2127 Revealed using dirty offsetof(G4UnionSolid, member) checks using:: 
2128 
2129    #define private   public
2130    #define protected public
2131 
2132 That is not a problem, just a little leak. 
2133 The reason is that the transform ctor allocates, so will normally get a different fPtrSolidB::
2134 
2135      77 G4BooleanSolid::G4BooleanSolid( const G4String& pName,
2136      78                                       G4VSolid* pSolidA ,
2137      79                                       G4VSolid* pSolidB ,
2138      80                                       G4RotationMatrix* rotMatrix,
2139      81                                 const G4ThreeVector& transVector    ) :
2140      82   G4VSolid(pName), fStatistics(1000000), fCubVolEpsilon(0.001),
2141      83   fAreaAccuracy(-1.), fCubicVolume(-1.), fSurfaceArea(-1.),
2142      84   fRebuildPolyhedron(false), fpPolyhedron(0), fPrimitivesSurfaceArea(0.),
2143      85   createdDisplacedSolid(true)
2144      86 {
2145      87   fPtrSolidA = pSolidA ;
2146      88   fPtrSolidB = new G4DisplacedSolid("placedB",pSolidB,rotMatrix,transVector) ;
2147      89 }
2148 
2149 **/
2150 
2151 void U4SolidTree::PlacementNewDupe( G4VSolid* solid) // static
2152 {
2153     G4BooleanSolid* src = dynamic_cast<G4BooleanSolid*>(solid) ; 
2154     assert( src ); 
2155 
2156     G4String name = src->GetName() ; 
2157     G4VSolid* left = src->GetConstituentSolid(0) ; 
2158     G4VSolid* disp_right = src->GetConstituentSolid(1) ;  // may be G4DisplacedSolid
2159 
2160     G4RotationMatrix rrot ; 
2161     G4ThreeVector    rtla ; 
2162     G4VSolid* right = Moved_( &rrot, &rtla, disp_right ); 
2163     int type = EntityType(solid); 
2164 
2165     G4SolidStore::GetInstance()->DeRegister(solid);
2166 
2167     G4VSolid* dupe = nullptr ; 
2168     switch(type)
2169     {
2170         case _G4UnionSolid        : dupe = new (solid) G4UnionSolid(        name, left, right, &rrot, rtla ) ; break ; 
2171         case _G4SubtractionSolid  : dupe = new (solid) G4SubtractionSolid(  name, left, right, &rrot, rtla ) ; break ;
2172         case _G4IntersectionSolid : dupe = new (solid) G4IntersectionSolid( name, left, right, &rrot, rtla ) ; break ; 
2173     } 
2174     bool expect = dupe == solid ; 
2175     assert(expect); 
2176     if(!expect) exit(EXIT_FAILURE); 
2177 
2178 } 
2179 
2180 /**
2181 U4SolidTree::SetRight
2182 -----------------
2183 
2184 The right hand side CSG constituent of *node* is changed to the provided *right* G4VSolid, 
2185 with the transform rotation/translation applied to the *right* solid.
2186 Use nullptr for rrot or rtla when no rotation or translation is required.
2187 
2188 As there is no G4BooleanSolid::SetConstituentSolid method this sneakily replaces *node* with 
2189 another at the same memory addess (using placement new) with the *right* changed by re-construction.
2190 
2191 Note that it is not appropriate to use a G4DisplacedSolid argument for either
2192 *solid* or *right* as G4DisplacedSolid objects are internal implementation details of 
2193 G4BooleanSolid that automatically gets instanciated by the G4BooleanSolid ctors 
2194 with transform arguments. So that means that G4DisplacedSolid are not suitable inputs
2195 to G4BooleanSolid ctors. 
2196 
2197 **/
2198 
2199 void U4SolidTree::SetRight(  G4VSolid* node, G4VSolid* right, G4RotationMatrix* rrot, G4ThreeVector* rtla )
2200 {
2201     assert( dynamic_cast<G4DisplacedSolid*>(node) == nullptr ) ; 
2202     assert( dynamic_cast<G4DisplacedSolid*>(right) == nullptr ) ; 
2203 
2204     G4BooleanSolid* src = dynamic_cast<G4BooleanSolid*>(node) ; 
2205     assert( src ); 
2206 
2207     int type = EntityType(src); 
2208     G4String name = src->GetName() ; 
2209     G4VSolid* left = src->GetConstituentSolid(0) ; 
2210 
2211     G4SolidStore::GetInstance()->DeRegister(src);
2212 
2213     G4VSolid* replacement = nullptr ; 
2214 
2215     G4ThreeVector tlate(0.,0.,0.); 
2216     if( rtla == nullptr ) rtla = &tlate ; 
2217     switch(type)
2218     {
2219         case _G4UnionSolid        : replacement = new (src) G4UnionSolid(        name, left, right, rrot, *rtla ) ; break ; 
2220         case _G4SubtractionSolid  : replacement = new (src) G4SubtractionSolid(  name, left, right, rrot, *rtla ) ; break ;
2221         case _G4IntersectionSolid : replacement = new (src) G4IntersectionSolid( name, left, right, rrot, *rtla ) ; break ; 
2222     }
2223 
2224     bool expect = replacement == src ;
2225     assert(expect); 
2226     if(!expect) exit(EXIT_FAILURE); 
2227 }
2228 
2229 /**
2230 U4SolidTree::SetLeft
2231 -------------------
2232 
2233 The lefthand side constituent of *node* is changed to *left* 
2234 This is implemented using placement new trickery to replace the 
2235 *node* with another at the same memory address with changed *left*.
2236 
2237 **/
2238 
2239 void U4SolidTree::SetLeft(  G4VSolid* node, G4VSolid* left)  // static 
2240 {
2241     assert( dynamic_cast<G4DisplacedSolid*>(node) == nullptr ) ; 
2242     assert( dynamic_cast<G4DisplacedSolid*>(left) == nullptr ) ; 
2243 
2244     G4BooleanSolid* src = dynamic_cast<G4BooleanSolid*>(node) ; 
2245     assert( src ); 
2246 
2247     int type = EntityType(src); 
2248     G4String name = src->GetName() ; 
2249     G4VSolid* disp_right = src->GetConstituentSolid(1) ;  // may be G4DisplacedSolid
2250 
2251     G4RotationMatrix rrot ; 
2252     G4ThreeVector    rtla(0.,0.,0.) ; 
2253     G4VSolid* right = Moved_( &rrot, &rtla, disp_right );  // if *right* isa G4DisplacedSolid get the moved solid inside and transforms
2254  
2255     G4SolidStore::GetInstance()->DeRegister(src);
2256 
2257     G4VSolid* replacement = nullptr ; 
2258     switch(type)
2259     {
2260         case _G4UnionSolid        : replacement = new (src) G4UnionSolid(        name, left, right, &rrot, rtla ) ; break ; 
2261         case _G4SubtractionSolid  : replacement = new (src) G4SubtractionSolid(  name, left, right, &rrot, rtla ) ; break ;
2262         case _G4IntersectionSolid : replacement = new (src) G4IntersectionSolid( name, left, right, &rrot, rtla ) ; break ; 
2263     } 
2264     bool expect = replacement == src ; 
2265     assert(expect); 
2266     if(!expect) exit(EXIT_FAILURE); 
2267 }
2268 
2269 
2270 /**
2271 U4SolidTree::PrimitiveClone
2272 -----------------------------
2273 
2274 Create G4VSolid using copy ctor of appropriate EntityType
2275 
2276 **/
2277 
2278 G4VSolid* U4SolidTree::PrimitiveClone( const  G4VSolid* solid )  // static 
2279 {
2280     LOG(LEVEL) << EntityTypeName(solid) ; 
2281 
2282     G4VSolid* clone = nullptr ; 
2283     int type = EntityType(solid); 
2284     if( type == _G4Ellipsoid )
2285     {
2286         const G4Ellipsoid* ellipsoid = dynamic_cast<const G4Ellipsoid*>(solid) ; 
2287         clone = new G4Ellipsoid(*ellipsoid) ;
2288     }
2289     else if( type == _G4Tubs && PROMOTE_TUBS_TO_POLYCONE == false)
2290     {
2291         const G4Tubs* tubs = dynamic_cast<const G4Tubs*>(solid) ; 
2292         clone = new G4Tubs(*tubs) ;  
2293     }
2294     else if( type == _G4Tubs && PROMOTE_TUBS_TO_POLYCONE == true)
2295     {
2296         clone = PromoteTubsToPolycone( solid );  
2297     }
2298     else if( type == _G4Polycone )
2299     {
2300         const G4Polycone* polycone = dynamic_cast<const G4Polycone*>(solid) ; 
2301         clone = new G4Polycone(*polycone) ;  
2302     }
2303     else if( type == _G4Torus )
2304     {
2305         const G4Torus* torus = dynamic_cast<const G4Torus*>(solid) ; 
2306         clone = new G4Torus(*torus) ;  
2307     }
2308     else if( type == _G4Box )
2309     {
2310         const G4Box* box = dynamic_cast<const G4Box*>(solid) ; 
2311         clone = new G4Box(*box) ;  
2312     }
2313     else if( type == _G4Orb )
2314     {
2315         const G4Orb* orb = dynamic_cast<const G4Orb*>(solid) ; 
2316         clone = new G4Orb(*orb) ;  
2317     }
2318     else if( type == _G4MultiUnion )
2319     {
2320         const G4MultiUnion* mun = dynamic_cast<const G4MultiUnion*>(solid) ; 
2321         clone = new G4MultiUnion(*mun) ;  
2322     }
2323     else if( type == _G4Sphere )
2324     {
2325         const G4Sphere* sph = dynamic_cast<const G4Sphere*>(solid) ; 
2326         clone = new G4Sphere(*sph) ;  
2327     }
2328     else
2329     {
2330         std::cout 
2331             << "U4SolidTree::PrimitiveClone FATAL unimplemented prim type " << EntityTypeName(solid) 
2332             << std::endl 
2333             ;
2334         assert(0); 
2335     } 
2336     return clone ; 
2337 }
2338 
2339 
2340 G4VSolid* U4SolidTree::PrimitiveClone( const  G4VSolid* solid, const char* name )  // static 
2341 {
2342     G4VSolid* clone = PrimitiveClone(solid); 
2343     clone->SetName(name); 
2344     return clone ; 
2345 }
2346 
2347 
2348 //const bool U4SolidTree::PROMOTE_TUBS_TO_POLYCONE = true ; 
2349 const bool U4SolidTree::PROMOTE_TUBS_TO_POLYCONE = false ; 
2350 
2351 G4VSolid* U4SolidTree::PromoteTubsToPolycone( const G4VSolid* solid ) // static
2352 {
2353     const G4Tubs* tubs = dynamic_cast<const G4Tubs*>(solid) ; 
2354     assert(tubs); 
2355 
2356     G4String name = tubs->GetName();
2357 #if G4VERSION_NUMBER < 1100
2358     double dz = tubs->GetDz(); 
2359     double rmin = tubs->GetRMin(); 
2360     double rmax = tubs->GetRMax(); 
2361     double sphi = tubs->GetSPhi(); 
2362     double dphi = tubs->GetDPhi(); 
2363 #else
2364     double dz   = tubs->GetZHalfLength();
2365     double rmin = tubs->GetInnerRadius();
2366     double rmax = tubs->GetOuterRadius();
2367     double sphi = tubs->GetStartPhiAngle();
2368     double dphi = tubs->GetDeltaPhiAngle();
2369 #endif
2370     
2371     G4int numZPlanes = 2 ; 
2372     G4double zPlane[] = { -dz    , dz } ;   
2373     G4double rInner[] = {  rmin  , rmin   } ;   
2374     G4double rOuter[] = {  rmax  , rmax } ;    
2375 
2376     G4Polycone* polycone = new G4Polycone(
2377                                name,
2378                                sphi,
2379                                dphi,
2380                                numZPlanes,
2381                                zPlane,
2382                                rInner,
2383                                rOuter
2384                                );  
2385 
2386     G4VSolid* clone = polycone ; 
2387     return clone ; 
2388 }
2389 
2390 
2391 double U4SolidTree::getX( const G4VSolid* node ) const
2392 {
2393     G4RotationMatrix* tree_rot = nullptr ; 
2394     G4ThreeVector    tree_tla(0., 0., 0. ); 
2395     getTreeTransform(tree_rot, &tree_tla, node ); 
2396 
2397     double x_delta = tree_tla.x() ; 
2398     return x_delta ; 
2399 }
2400 
2401 double U4SolidTree::getY( const G4VSolid* node ) const
2402 {
2403     G4RotationMatrix* tree_rot = nullptr ; 
2404     G4ThreeVector    tree_tla(0., 0., 0. ); 
2405     getTreeTransform(tree_rot, &tree_tla, node ); 
2406 
2407     double y_delta = tree_tla.y() ; 
2408     return y_delta ; 
2409 }
2410 
2411 double U4SolidTree::getZ( const G4VSolid* node ) const
2412 {
2413     G4RotationMatrix* tree_rot = nullptr ; 
2414     G4ThreeVector    tree_tla(0., 0., 0. ); 
2415     getTreeTransform(tree_rot, &tree_tla, node ); 
2416 
2417     double z_delta = tree_tla.z() ; 
2418     return z_delta ; 
2419 }
2420 
2421 
2422 
2423 
2424 
2425 
2426 bool U4SolidTree::CanX(  const G4VSolid* solid ) 
2427 {
2428     int type = EntityType(solid) ; 
2429     return type == _G4Box ; 
2430 }
2431 bool U4SolidTree::CanY(  const G4VSolid* solid ) 
2432 {
2433     int type = EntityType(solid) ; 
2434     return type == _G4Box ; 
2435 }
2436 bool U4SolidTree::CanZ( const G4VSolid* solid ) // static
2437 {
2438     int type = EntityType(solid) ; 
2439     bool can = type == _G4Ellipsoid || type == _G4Tubs || type == _G4Polycone || type == _G4Torus || type == _G4Box ; 
2440     G4String name = solid->GetName(); 
2441 
2442     if( can == false && verbose )
2443     {
2444         std::cout 
2445             << "U4SolidTree::CanZ ERROR "
2446             << " false for solid.EntityTypeName " << EntityTypeName(solid)
2447             << " solid.name " << name.c_str() 
2448             << std::endl ; 
2449             ;
2450     }
2451 
2452     return can ; 
2453 }
2454 
2455 
2456 void U4SolidTree::XRange( double& x0, double& x1, const G4VSolid* solid) // static  
2457 {
2458     switch(EntityType(solid))
2459     {
2460         case _G4Box:       GetXRange( dynamic_cast<const G4Box*>(solid),       x0, x1 );  break ; 
2461         case _G4Other:    { std::cout << "U4SolidTree::GetX FATAL : not implemented for entityType " << EntityTypeName(solid) << std::endl ; assert(0) ; } ; break ;  
2462     }
2463 }
2464 void U4SolidTree::YRange( double& y0, double& y1, const G4VSolid* solid) // static  
2465 {
2466     switch(EntityType(solid))
2467     {
2468         case _G4Box:       GetYRange( dynamic_cast<const G4Box*>(solid),       y0, y1 );  break ; 
2469         case _G4Other:    { std::cout << "U4SolidTree::GetY FATAL : not implemented for entityType " << EntityTypeName(solid) << std::endl ; assert(0) ; } ; break ;  
2470     }
2471 }
2472 void U4SolidTree::ZRange( double& z0, double& z1, const G4VSolid* solid) // static  
2473 {
2474     switch(EntityType(solid))
2475     {
2476         case _G4Ellipsoid: GetZRange( dynamic_cast<const G4Ellipsoid*>(solid), z0, z1 );  break ; 
2477         case _G4Tubs:      GetZRange( dynamic_cast<const G4Tubs*>(solid)    ,  z0, z1 );  break ; 
2478         case _G4Polycone:  GetZRange( dynamic_cast<const G4Polycone*>(solid),  z0, z1 );  break ; 
2479         case _G4Torus:     GetZRange( dynamic_cast<const G4Torus*>(solid),     z0, z1 );  break ; 
2480         case _G4Box:       GetZRange( dynamic_cast<const G4Box*>(solid),       z0, z1 );  break ; 
2481         case _G4Other:    { std::cout << "U4SolidTree::GetZ FATAL : not implemented for entityType " << EntityTypeName(solid) << std::endl ; assert(0) ; } ; break ;  
2482     }
2483 }
2484 
2485 
2486 
2487 void U4SolidTree::GetXRange( const G4Box* const box, double& _x0, double& _x1 )  // static 
2488 {
2489     _x1 = box->GetXHalfLength() ; 
2490     _x0 = -_x1 ; 
2491     assert( _x1 > 0. ); 
2492 }
2493 void U4SolidTree::GetYRange( const G4Box* const box, double& _y0, double& _y1 )  // static 
2494 {
2495     _y1 = box->GetYHalfLength() ; 
2496     _y0 = -_y1 ; 
2497     assert( _y1 > 0. ); 
2498 }
2499 void U4SolidTree::GetZRange( const G4Box* const box, double& _z0, double& _z1 )  // static 
2500 {
2501     _z1 = box->GetZHalfLength() ; 
2502     _z0 = -_z1 ; 
2503     assert( _z1 > 0. ); 
2504 }
2505 
2506 
2507 
2508 
2509 void U4SolidTree::GetZRange( const G4Ellipsoid* const ellipsoid, double& _z0, double& _z1 )  // static 
2510 {
2511     _z1 = ellipsoid->GetZTopCut() ; 
2512     _z0 = ellipsoid->GetZBottomCut() ;  
2513 }
2514 void U4SolidTree::GetZRange( const G4Tubs* const tubs, double& _z0, double& _z1 )  // static 
2515 {
2516     _z1 = tubs->GetZHalfLength() ;  
2517     _z0 = -_z1 ;  
2518     assert( _z1 > 0. ); 
2519 }
2520 void U4SolidTree::GetZRange( const G4Polycone* const polycone, double& _z0, double& _z1 )  // static 
2521 {
2522     G4PolyconeHistorical* pars = polycone->GetOriginalParameters(); 
2523     unsigned num_z = pars->Num_z_planes ; 
2524     for(unsigned i=1 ; i < num_z ; i++)
2525     {
2526         double z0 = pars->Z_values[i-1] ; 
2527         double z1 = pars->Z_values[i] ; 
2528         bool expect = z1 > z0 ; 
2529         assert(expect);   
2530         if(!expect) exit(EXIT_FAILURE); 
2531     }
2532     _z1 = pars->Z_values[num_z-1] ; 
2533     _z0 = pars->Z_values[0] ;  
2534 }
2535 void U4SolidTree::GetZRange( const G4Torus* const torus, double& _z0, double& _z1 )  // static 
2536 {
2537     G4double rmax = torus->GetRmax() ; 
2538     _z1 = rmax ; 
2539     _z0 = -rmax ;  
2540 }
2541 
2542 
2543 
2544 /**
2545 U4SolidTree::ApplyZCut_G4Tubs
2546 ----------------------------
2547 
2548 An alternative to using this is to promote all G4Tubs to G4Polycone at the clone stage,  
2549 avoiding the need to change boolean displacements as G4Polycone is not symmetrically defined. 
2550 
2551 * INITIALLY HAD PROBLEM WITH THE IMPLEMENTATION : CHANGING OFFSETS DID NOT WORK
2552 * SOLVED WITH PLACEMENT NEW APPLIED TO THE G4DisplacedSolid 
2553 
2554 Cutting G4Tubs::
2555 
2556 
2557      zd+hz  +---------+               +---------+     new_zd + new_hz
2558             |         |               |         |  
2559             |         |               |         |
2560             |         |               |         |
2561             |         |             __|_________|__   new_zd
2562             |         |               |         |
2563      zd   --|---------|--             |         |
2564             |         |               |         |
2565             |         |               |         |
2566          .  | . . . . | . .zcut . . . +---------+ . . new_zd - new_hz  . . . . . .
2567             |         | 
2568             |         |
2569     zd-hz   +---------+ 
2570 
2571 
2572      cut position
2573 
2574           zcut = new_zd - new_hz 
2575           new_zd = zcut + new_hz  
2576 
2577 
2578      original height:  2*hz                         
2579       
2580      cut height :     
2581 
2582           loc_zcut = zcut - zd 
2583 
2584           2*new_hz = 2*hz - (zcut-(zd-hz)) 
2585 
2586                    = 2*hz - ( zcut - zd + hz )
2587 
2588                    = 2*hz -  zcut + zd - hz 
2589 
2590                    = hz + zd - zcut 
2591 
2592                    = hz - (zcut - zd)             
2593 
2594                    = hz - loc_zcut 
2595 
2596 
2597                                     hz + zd - zcut
2598                         new_hz =  -----------------
2599                                          2
2600 
2601                                     hz - loc_zcut 
2602                         new_hz =  --------------------      new_hz( loc_zcut:-hz ) = hz     unchanged
2603                                           2                 new_hz( loc_zcut:0   ) = hz/2   halved
2604                                                             new_hz( loc_zcut:+hz ) =  0     made to disappear 
2605 
2606 
2607 
2608 Simpler way to derive the same thing, is to use the initial local frame::
2609 
2610 
2611        +hz  +---------+               +---------+     zoff + new_hz
2612             |         |               |         |  
2613             |         |               |         |
2614             |         |               |         |
2615             |         |             __|_________|__   zoff 
2616             |         |               |         |
2617         0 --|---------|- . . . . . . . . . . . . . . . 0 
2618             |         |               |         |
2619             |         |               |         |
2620    loc_zcut | . . . . | . . . . .  .  +---------+ . . zoff - new_hz  . . . . . .
2621             |         | 
2622             |         |
2623       -hz   +---------+. . . . . . . . . . . . 
2624 
2625 
2626 
2627             loc_zcut = zoff - new_hz
2628 
2629                 zoff = loc_zcut + new_hz 
2630 
2631                      = loc_zcut +  (hz - loc_zcut)/2
2632 
2633                      =  2*loc_zcut + (hz - loc_zcut)
2634                         ------------------------------
2635                                    2
2636 
2637                zoff  =  loc_zcut + hz                   zoff( loc_zcut:-hz ) = 0      unchanged
2638                         ---------------                 zoff( loc_zcut:0   ) = hz/2   
2639                               2                         zoff( loc_zcut:+hz ) = hz     makes sense, 
2640                                                                                       think about just before it disappears
2641 
2642 
2643 Simpler way,  notice the top and bottom line equations, add or subtract and rearrange::
2644 
2645             hz       = zoff + new_hz
2646 
2647             loc_zcut = zoff - new_hz
2648 
2649 
2650       Add them:
2651 
2652             hz + loc_zcut = 2*zoff 
2653 
2654                    zoff =  hz + loc_zcut
2655                            --------------
2656                                  2 
2657       Subtract them:
2658 
2659                hz - loc_zcut = 2*new_hz
2660 
2661                    new_hz = hz - loc_zcut 
2662                             ---------------
2663                                  2
2664 
2665 **/
2666
[ 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!