eic-opticks/u4/U4Polycone.h

0001 #pragma once
0002 /**
0003 U4Polycone.h
0004 ===============
0005
0006 Polycone Z-Nudging
0007 -------------------
0008
0009 Polycone Z-nudging is simpler than the general case because:
0010
0011 1. no transforms
0012 2. no need to look for coincidences, as every Z-joint is coincident
0013    and every subtracted inner end face is coincident with the outer that
0014    it is subtracted from.
0015
0016
0017 Polycone example
0018 ----------------
0019
0020 Initial::
0021
0022     U4Polycone::desc num 4 rz 4 R_inner 2 R_outer 2 Z 3
0023       0 RZ     43.000    195.000      0.000
0024       1 RZ     43.000    195.000    -15.000
0025       2 RZ     55.500     70.000    -15.000
0026       3 RZ     55.500     70.000   -101.000
0027
0028 After reversal::
0029
0030     U4Polycone::desc num 4 rz 4 R_inner 2 R_outer 2 Z 3
0031       0 RZ     55.500     70.000   -101.000
0032       1 RZ     55.500     70.000    -15.000
0033       2 RZ     43.000    195.000    -15.000
0034       3 RZ     43.000    195.000      0.000
0035                rmin      rmax          z
0036
0037
0038
0039        -195                 -43    :    43                  195
0040
0041          0___________________0     :     0___________________0                 z = 0      (0)
0042          |                   |     :     |                   |
0043          |                   |     :     |                   |
0044          1___________________1     :     1___________________1                 z = -15    (1)
0045                      2    2        :        2    2                             z = -15    (2)
0046                      |    |        :        |    |
0047                      |    |        :        |    |
0048                      |    |        :        |    |
0049                      3____3        :        3____3                             z = -101   (3)
0050
0051         -195       -70  -55.5      0       55.5  70         195
0052
0053
0054 **/
0055
0056
0057 #include <vector>
0058 #include <set>
0059 #include <algorithm>
0060
0061 #include "G4Polycone.hh"
0062
0063 #include "sn.h"
0064 #include "ssys.h"
0065
0066 struct RZ
0067 {
0068     double rmin ;
0069     double rmax ;
0070     double z ;
0071
0072     std::string desc() const ;
0073 };
0074
0075 inline std::string RZ::desc() const
0076 {
0077     std::stringstream ss ;
0078     ss << "RRZ"
0079        << " " << std::setw(10) << std::fixed << std::setprecision(3) << rmin
0080        << " " << std::setw(10) << std::fixed << std::setprecision(3) << rmax
0081        << " " << std::setw(10) << std::fixed << std::setprecision(3) << z
0082        ;
0083     std::string str = ss.str();
0084     return str ;
0085 }
0086
0087
0088 struct U4Polycone
0089 {
0090     static sn*  Convert( const G4Polycone* polycone, int lvid, int depth, int level );
0091 private:
0092
0093     std::string desc() const ;
0094     static void GetMinMax( double& mn, double& mx, const std::set<double>& vv );
0095
0096
0097     U4Polycone(const G4Polycone* poly, int lvid, int depth, int level );
0098     bool checkZOrder( bool z_ascending );
0099     void init();
0100     void init_RZ();
0101     void init_phicut();
0102     void init_outer();
0103     void init_inner();
0104
0105     void collectPrims(std::vector<sn*>& prims, bool outside );
0106
0107
0108     // MEMBERS
0109
0110     int lvid ;
0111     int depth ;
0112     int level ;
0113
0114     bool enable_nudge ;
0115     const G4Polycone* polycone ;
0116     const G4PolyconeHistorical* ph ;
0117
0118     int num ;
0119     std::vector<RZ> rz ;
0120
0121
0122     static constexpr const char* U4Polycone__ENABLE_PHICUT = "U4Polycone__ENABLE_PHICUT" ;
0123     bool   ENABLE_PHICUT ;
0124     double phi_start ;
0125     double phi_end ;
0126     sn*    phicut ;
0127
0128     std::set<double> R_inner ;
0129     std::set<double> R_outer ;
0130     std::set<double> Z ;
0131
0132     int num_R_inner ;
0133     double R_inner_min ;
0134     double R_inner_max ;
0135
0136     int num_R_outer ;
0137     double R_outer_min ;
0138     double R_outer_max ;
0139
0140     int num_Z ;
0141     double Z_min ;
0142     double Z_max ;
0143
0144     bool   has_inner ;
0145
0146     std::vector<sn*> outer_prims ;
0147     std::vector<sn*> inner_prims ;
0148     sn*    inner ;
0149     sn*    outer ;
0150     sn*    root ;
0151
0152     const char* label ;
0153
0154 };
0155
0156
0157 inline sn* U4Polycone::Convert( const G4Polycone* polycone, int lvid, int depth, int level )
0158 {
0159     U4Polycone upoly(polycone, lvid, depth, level ) ;
0160
0161     if(level > 0)
0162     {
0163        std::cerr << "U4Polycone::Convert" << std::endl ;
0164        std::cerr << upoly.root->render(5) ;
0165     }
0166     return upoly.root ;
0167 }
0168
0169
0170
0171
0172 inline std::string U4Polycone::desc() const
0173 {
0174     std::stringstream ss ;
0175     ss << "U4Polycone::desc"
0176        << " lvid " << lvid
0177        << " depth " << depth
0178        << " level " << level
0179        << " enable_nudge " << ( enable_nudge ? "YES" : "NO " )
0180        << " num " << num
0181        << " ENABLE_PHICUT " << std::setw(10) << ( ENABLE_PHICUT ? "YES" : "NO " )
0182        << " phi_start " << std::setw(10) << phi_start
0183        << " phi_end " << std::setw(10) << phi_end
0184        << " phicut " << std::setw(10) << ( phicut ? "YES" : "NO " )
0185        << " rz " << rz.size()
0186        << std::endl
0187        << " num_R_inner " << std::setw(3) << num_R_inner
0188        << " R_inner_min " << std::setw(10) << R_inner_min
0189        << " R_inner_max " << std::setw(10) << R_inner_max
0190        << std::endl
0191        << " num_R_outer " << std::setw(3) << num_R_outer
0192        << " R_outer_min " << std::setw(10) << R_outer_min
0193        << " R_outer_max " << std::setw(10) << R_outer_max
0194        << std::endl
0195        << " num_Z       " << std::setw(3) << num_Z
0196        << " Z_min       " << std::setw(10) << Z_min
0197        << " Z_max       " << std::setw(10) << Z_max
0198        << std::endl
0199        << " has_inner " << ( has_inner ? "YES" : "NO" )
0200        << " root " << std::setw(3) << ( root ? root->index() : -1 )
0201        << " label " << ( label ? label : "-" )
0202        << std::endl
0203        ;
0204
0205     for(int i=0 ; i < num ; i++) ss << std::setw(3) << i << " " << rz[i].desc() << std::endl ;
0206     std::string str = ss.str();
0207     return str ;
0208 }
0209
0210
0211
0212
0213 inline void U4Polycone::GetMinMax( double& mn, double& mx, const std::set<double>& vv )
0214 {
0215     mn = *vv.begin() ;
0216     mx = *vv.begin() ;
0217     for(std::set<double>::const_iterator it = vv.begin() ; it != vv.end() ; it++ )
0218     {
0219         mn = std::min( mn, *it );
0220         mx = std::max( mx, *it );
0221     }
0222 }
0223
0224 inline U4Polycone::U4Polycone(const G4Polycone* polycone_, int lvid_, int depth_, int level_ )
0225     :
0226     lvid(lvid_),
0227     depth(depth_),
0228     level(level_),
0229     enable_nudge(!ssys::getenvbool("U4Polycone__DISABLE_NUDGE")),
0230     polycone(polycone_),
0231     ph(polycone->GetOriginalParameters()),
0232     num(ph->Num_z_planes),
0233     ENABLE_PHICUT(ssys::getenvbool(U4Polycone__ENABLE_PHICUT)),
0234     phi_start(0.),
0235     phi_end(2.*M_PI),
0236     phicut(nullptr),
0237     num_R_inner(0),
0238     R_inner_min(0),
0239     R_inner_max(0),
0240     num_R_outer(0),
0241     R_outer_min(0),
0242     R_outer_max(0),
0243     num_Z(0),
0244     Z_min(0),
0245     Z_max(0),
0246     has_inner(false),
0247     inner(nullptr),
0248     outer(nullptr),
0249     root(nullptr),
0250     label("NOT-WITH-SND")
0251 {
0252     init();
0253     if(level > 0 ) std::cerr
0254         << "U4Polycone::U4Polycone "
0255         << ( label ? label : "-" )
0256         << std::endl
0257         << desc()
0258         << std::endl
0259         ;
0260
0261 }
0262
0263 inline bool U4Polycone::checkZOrder( bool z_ascending )
0264 {
0265     int count_z_order = 0 ;
0266     for( int i=1 ; i < num ; i++)
0267     {
0268         bool z_order = z_ascending ? rz[i-1].z <= rz[i].z : rz[i].z <= rz[i-1].z ;
0269         if(z_order) count_z_order += 1 ;
0270     }
0271     bool all_z_order = count_z_order == num - 1 ;
0272     return all_z_order ;
0273 }
0274
0275
0276 /**
0277 U4Polycone::init
0278 -----------------
0279
0280 inner and outer can be sn::Cylinder OR sn::Collection (eg sn::UnionTree)
0281 and root can be those also OR sn::Boolean difference of them.
0282 So the technique used to support phi range needs to be applicable
0283 to many types of node including trees and singles.
0284
0285 Ideas how to do that:
0286
0287 1. extend s_pa from 6 to 8 elem and include phi range there,
0288    this has disadvantage of adding empty params to almost all s_pa
0289    that are only used for those with phi segment (and a subsequent
0290    theta segment would push from 8 to 10)
0291
0292 2. add s_au for auxiliary params analogous to how s_pa is related to sn,
0293    that is quite a lot of code change
0294
0295 3. HMM: use boolean intersection with a special segment (CSG_PHICUT)
0296    type that the phi (and in future theta perhaps) range in its param
0297
0298    * THIS IS ADVANTAGEOUS FROM POINT OF VIEW OF EXPRESSING THE INFO
0299      WITHOUT MUCH CHANGE TO EXISTING CODE
0300
0301 **/
0302
0303
0304 inline void U4Polycone::init()
0305 {
0306     init_RZ();
0307     init_phicut();
0308     init_outer();
0309
0310     sn* _root = nullptr ;
0311     {
0312         if(has_inner == false)
0313         {
0314             _root = outer ;
0315         }
0316         else
0317         {
0318             init_inner();
0319             assert( inner );
0320             _root = sn::Boolean(CSG_DIFFERENCE, outer, inner );
0321         }
0322     }
0323
0324
0325     if(phicut == nullptr)
0326     {
0327         root = _root ;
0328     }
0329     else
0330     {
0331         if(ENABLE_PHICUT)
0332         {
0333             root = sn::Boolean(CSG_INTERSECTION, phicut, _root );
0334         }
0335         else
0336         {
0337             std::cerr
0338                << "U4Polycone::init FATAL geometry with unsupported phicut : "
0339                << " enable experimental support with envvar "
0340                << " [" <<  U4Polycone__ENABLE_PHICUT << "]"
0341                << "\n"
0342                ;
0343
0344             assert(0);
0345             std::raise(SIGINT);
0346         }
0347     }
0348
0349 }
0350
0351 /**
0352 U4Polycone::init_RZ
0353 ---------------------
0354
0355 1. fill (RZ)rz vector and insert values into std::set
0356    to give number of unique rmin, rmax, z
0357
0358 2. if necessary reverse the rz vector to make all z ascending
0359
0360 3. get the min/max ranges of rmin, rmax, z and determine if there
0361    is an inner based on the rmin range
0362
0363
0364 **/
0365
0366 inline void U4Polycone::init_RZ()
0367 {
0368     rz.resize(num);
0369
0370     for (int i=0; i < num ; i++)
0371     {
0372         RZ& rzi = rz[i] ;
0373         rzi.rmin = ph->Rmin[i] ;
0374         rzi.rmax = ph->Rmax[i] ;
0375         rzi.z    = ph->Z_values[i] ;
0376
0377         R_inner.insert(rzi.rmin);
0378         R_outer.insert(rzi.rmax);
0379         Z.insert(rzi.z);
0380     }
0381
0382     num_R_inner = R_inner.size();
0383     num_R_outer = R_outer.size();
0384     num_Z = Z.size();
0385
0386
0387     bool all_z_descending = checkZOrder(false);
0388     if(all_z_descending)
0389     {
0390         if(level > 0) std::cerr
0391            << "U4Polycone::init_RZ"
0392            << label
0393            << " all_z_descending detected, reversing "
0394            << std::endl
0395            ;
0396         std::reverse( std::begin(rz), std::end(rz) ) ;
0397     }
0398     bool all_z_ascending  = checkZOrder(true  );
0399     assert( all_z_ascending );
0400     if(!all_z_ascending) std::raise(SIGINT);
0401
0402     GetMinMax(R_inner_min, R_inner_max, R_inner);
0403     GetMinMax(R_outer_min, R_outer_max, R_outer);
0404     GetMinMax(Z_min, Z_max, Z);
0405
0406     assert( Z_max > Z_min );
0407     bool no_inner = R_inner_min == 0. && R_inner_max == 0. ;
0408     has_inner = !no_inner ;
0409 }
0410
0411
0412 /**
0413 U4Polycone::init_phicut
0414 -------------------------
0415
0416
0417
0418            Y
0419      .     |      n
0420        .   |    .
0421          . | .
0422            +------X
0423             .
0424               .
0425
0426
0427 Special cased "hemi" for simpler testing::
0428
0429                              +Y
0430                               .
0431                     + + + + + . + + + + +
0432                     + + + + + . + + + + +
0433                     + + + + + . + + + + +
0434                     + + + + + . + + + + +
0435                     + + + + + . + + + + +
0436   phi_end = pi  -X +----------+---------+ +X phi_start = 0.
0437                               :
0438                               :
0439                               :
0440                               n [0,-1,0]
0441                              -Y
0442
0443 Normal in -Y direction corresponds to the halfspace Y>0
0444
0445 TODO: handle all half-space situations not just one
0446
0447
0448 **/
0449
0450 inline void U4Polycone::init_phicut()
0451 {
0452     double eps = 1e-6 ;
0453     double phi_delta = ph->Opening_angle/CLHEP::radian ;
0454     phi_start = ph->Start_angle/CLHEP::radian ;
0455     phi_end   = phi_start + phi_delta ;
0456
0457     bool has_phicut = phi_start > 0. || phi_end < 2.0*CLHEP::pi  ;
0458     bool has_half = std::abs( phi_start - 0.f ) < eps  && std::abs(phi_end-CLHEP::pi) < eps ;
0459
0460     if( has_phicut )
0461     {
0462         if( has_half )
0463         {
0464             phicut = sn::HalfSpace( 0., -1., 0., 0. );
0465         }
0466         else
0467         {
0468             phicut = sn::PhiCut( phi_start, phi_end );
0469         }
0470     }
0471
0472     if(has_phicut) std::cerr
0473        << "U4Polycone::init_phicut"
0474        << " phi_start   " << std::setw(10) << std::fixed << std::setprecision(4) << phi_start
0475        << " phi_end     " << std::setw(10) << std::fixed << std::setprecision(4) << phi_end
0476        << " phi_delta   " << std::setw(10) << std::fixed << std::setprecision(4) << phi_delta
0477        << " has_phicut " << ( has_phicut ? "YES" : "NO " )
0478        << " has_half " << ( has_half ? "YES" : "NO " )
0479        << " ENABLE_PHICUT " << ( ENABLE_PHICUT ? "YES" : "NO " )
0480        << "\n"
0481        ;
0482
0483 }
0484
0485 /**
0486 U4Polycone::init_outer
0487 ------------------------
0488
0489 1. populate outer_prims vector with cones and cylinders
0490 2. when more than one outer prim invoke sn::ZNudgeOverlapJoints.
0491    This does not change geometry as it just changes internal joints
0492    between prims to avoid coincident faces (assuming sane prim sizes).
0493 3. collect the vector of prims into a binary union tree of nodes
0494
0495 **/
0496
0497 inline void U4Polycone::init_outer()
0498 {
0499     if( num_R_outer == 1 )
0500     {
0501         assert( R_outer_min == R_outer_max );
0502         outer = sn::Cylinder( R_outer_max, Z_min, Z_max );
0503     }
0504     else
0505     {
0506         collectPrims( outer_prims, true  ); // outside:true
0507         int num_outer_prim = outer_prims.size() ;
0508
0509         if(level > 0) std::cerr
0510             << "U4Polycone::init_outer."
0511             << " num_outer_prim " << num_outer_prim
0512             << std::endl
0513             ;
0514
0515         if(num_outer_prim > 1) sn::ZNudgeOverlapJoints(lvid, outer_prims, enable_nudge);
0516         outer = sn::Collection(outer_prims) ;
0517     }
0518
0519 }
0520
0521 /**
0522 U4Polycone::init_inner
0523 -----------------------
0524
0525 1. if there is only a single inner create a single cylinder
0526
0527    * HMM : what about a single cone inner ?
0528
0529
0530 **/
0531
0532 inline void U4Polycone::init_inner()
0533 {
0534     assert( has_inner ) ;
0535
0536     if(level > 0) std::cerr
0537        << "U4Polycone::init_inner "
0538        << std::endl
0539        ;
0540
0541     if( num_R_inner == 1 )  // cylinder inner
0542     {
0543         assert( R_inner_min == R_inner_max );
0544         inner = sn::Cylinder(R_inner_min, Z_min, Z_max);
0545     }
0546     else
0547     {
0548         collectPrims( inner_prims, false  ); // outside:false
0549         int num_inner_prim = inner_prims.size() ;
0550
0551         if(level > 0) std::cerr
0552             << "U4Polycone::init."
0553             << label
0554             << " num_inner_prim " << num_inner_prim
0555             << std::endl
0556             ;
0557
0558         sn::ZNudgeExpandEnds(lvid, inner_prims, enable_nudge);    // only for inner
0559         if(num_inner_prim > 1) sn::ZNudgeOverlapJoints(lvid, inner_prims, enable_nudge);
0560         inner = sn::Collection( inner_prims );
0561     }
0562 }
0563
0564
0565 /**
0566 U4Polycone::collectPrims
0567 --------------------------
0568
0569 Populate prims vectors with snd indices or sn pointers
0570 to cylinder or cone nodes created using values
0571 from the (RZ)rz vector. For outside:true use Rmax values
0572 otherwise use Rmin values for the inner.
0573
0574 **/
0575
0576 void U4Polycone::collectPrims(std::vector<sn*>& prims,  bool outside  )
0577 {
0578     bool inner = !outside ;
0579     // loop over pairs of planes
0580     int num_rz = rz.size();
0581
0582     for( int i=1 ; i < num_rz ; i++ )
0583     {
0584         const RZ& rz1 = rz[i-1] ;   // zplane struct rmin, rmax, z
0585         const RZ& rz2 = rz[i] ;
0586         double r1 = outside ? rz1.rmax : rz1.rmin ;
0587         double r2 = outside ? rz2.rmax : rz2.rmin ;
0588         double z1 = rz1.z ;
0589         double z2 = rz2.z ;
0590
0591         if( z1 == z2 )
0592         {
0593             if(level > 0) std::cerr << "U4Polycone::collectPrims skipping prim as z2 == z1  " << std::endl ;
0594             continue ;
0595         }
0596
0597
0598         bool z_ascending = z2 > z1 ;
0599         assert(z_ascending);
0600         if(!z_ascending) std::raise(SIGINT);
0601
0602         bool is_cylinder = r1 == r2 ;
0603         int idx = -1 ;
0604
0605         if( inner && r1 == 0. && r2 == 0. )
0606         {
0607             if(level > 0) std::cerr << "U4Polycone::collectPrims skipping inner as r1 == r2 == 0.  " << std::endl ;
0608             continue ;
0609         }
0610
0611         /*
0612         if(is_cylinder)
0613         {
0614              bool r2_nonzero = r2 > 0. ;
0615              if(!r2_nonzero) std::cerr
0616                   << "U4Polycone::collectPrims"
0617                   << " ERROR - CYLINDER PAIR UNEXPECTED ZERO r2 "
0618                   << " lvid " << lvid
0619                   << " r2_nonzero " << ( r2_nonzero ? "YES" : "NO ")
0620                   << " num_rz " << num_rz
0621                   << " i " << i
0622                   << "\n"
0623                   << desc()
0624                   << "\n"
0625                   ;
0626              assert(r2_nonzero);
0627         }
0628         */
0629
0630         sn* pr = is_cylinder ? sn::Cylinder(r2, z1, z2 ) : sn::Cone( r1, z1, r2, z2 ) ;
0631         pr->lvid = lvid ;  // so this before setting root for debug purposes
0632         prims.push_back(pr);
0633         idx = pr->index() ;
0634
0635         if( level > 0 ) std::cerr
0636              << "U4Polycone::collectPrims"
0637              << " outside " << ( outside ? "YES" : "NO " )
0638              << " idx " << idx
0639              << " is_cylinder " << ( is_cylinder ? "YES" : "NO " )
0640              << std::endl
0641              ;
0642     }
0643 }
0644
0645