eic-opticks/sysrap/sgeomtools.h

0001 #pragma once
0002 /**
0003 sgeomtools.h
0004 =============
0005
0006 Adapt some extracts from G4GeomTools
0007
0008 **/
0009
0010 #include <cmath>
0011 #include "scuda.h"
0012
0013 struct sgeomtools
0014 {
0015     static constexpr const double kCarTolerance = 1e-5 ;
0016
0017     static void Set( double2& p, double x, double y );
0018
0019     static bool DiskExtent(
0020             double rmin,
0021             double rmax,
0022             double startPhi,
0023             double deltaPhi,
0024             double2& pmin,
0025             double2& pmax);
0026
0027     static void DiskExtent(
0028             double rmin,
0029             double rmax,
0030             double sinPhiStart,
0031             double cosPhiStart,
0032             double sinPhiEnd,
0033             double cosPhiEnd,
0034             double2& pmin,
0035             double2& pmax);
0036 };
0037
0038
0039 inline void sgeomtools::Set( double2& p, double x, double y )
0040 {
0041     p.x = x ;
0042     p.y = y ;
0043 }
0044
0045 inline bool sgeomtools::DiskExtent(
0046             double rmin,
0047             double rmax,
0048             double startPhi,
0049             double deltaPhi,
0050             double2& pmin,
0051             double2& pmax)   // static
0052 {
0053     Set(pmin,0.,0.);
0054     Set(pmax,0.,0.);
0055
0056     if(rmin < 0.)    return false ;
0057     if(rmax <= rmin) return false ;
0058     if(deltaPhi <= 0.) return false ;
0059
0060     Set(pmin,-rmax,-rmax);
0061     Set(pmax, rmax, rmax);
0062
0063     if(deltaPhi >= M_PI*2. ) return true ;
0064     double endPhi = startPhi + deltaPhi ;
0065
0066     DiskExtent( rmin,
0067                 rmax,
0068                 std::sin(startPhi),
0069                 std::cos(startPhi),
0070                 std::sin(endPhi),
0071                 std::cos(endPhi),
0072                 pmin,
0073                 pmax);
0074
0075     return true ;
0076 }
0077
0078
0079 inline void sgeomtools::DiskExtent(
0080             double rmin,
0081             double rmax,
0082             double sinStart,
0083             double cosStart,
0084             double sinEnd,
0085             double cosEnd,
0086             double2& pmin,
0087             double2& pmax)   // static
0088 {
0089     Set(pmin,-rmax,-rmax);
0090     Set(pmax, rmax, rmax);
0091
0092   if (std::abs(sinEnd-sinStart) < kCarTolerance &&
0093       std::abs(cosEnd-cosStart) < kCarTolerance) return;
0094
0095   // get start and end quadrants
0096   //
0097   //      1 | 0
0098   //     ---+---
0099   //      3 | 2
0100   //
0101   int icase = (cosEnd < 0) ? 1 : 0;
0102   if (sinEnd   < 0) icase += 2;
0103   if (cosStart < 0) icase += 4;
0104   if (sinStart < 0) icase += 8;
0105
0106   switch (icase)
0107   {
0108   // start quadrant 0
0109   case  0:                                 // start->end : 0->0
0110     if (sinEnd < sinStart) break;
0111     Set(pmin,rmin*cosEnd,rmin*sinStart);
0112     Set(pmax,rmax*cosStart,rmax*sinEnd  );
0113     break;
0114   case  1:                                 // start->end : 0->1
0115     Set(pmin,rmax*cosEnd,std::min(rmin*sinStart,rmin*sinEnd));
0116     Set(pmax,rmax*cosStart,rmax  );
0117     break;
0118   case  2:                                 // start->end : 0->2
0119     Set(pmin,-rmax,-rmax);
0120     Set(pmax,std::max(rmax*cosStart,rmax*cosEnd),rmax);
0121     break;
0122   case  3:                                 // start->end : 0->3
0123     Set(pmin,-rmax,rmax*sinEnd);
0124     Set(pmax,rmax*cosStart,rmax);
0125     break;
0126   // start quadrant 1
0127   case  4:                                 // start->end : 1->0
0128     Set(pmin,-rmax,-rmax);
0129     Set(pmax,rmax,std::max(rmax*sinStart,rmax*sinEnd));
0130     break;
0131   case  5:                                 // start->end : 1->1
0132     if (sinEnd > sinStart) break;
0133     Set(pmin,rmax*cosEnd,rmin*sinEnd  );
0134     Set(pmax,rmin*cosStart,rmax*sinStart);
0135     break;
0136   case  6:                                 // start->end : 1->2
0137     Set(pmin,-rmax,-rmax);
0138     Set(pmax,rmax*cosEnd,rmax*sinStart);
0139     break;
0140   case  7:                                 // start->end : 1->3
0141     Set(pmin,-rmax,rmax*sinEnd);
0142     Set(pmax,std::max(rmin*cosStart,rmin*cosEnd),rmax*sinStart);
0143     break;
0144   // start quadrant 2
0145   case  8:                                 // start->end : 2->0
0146     Set(pmin,std::min(rmin*cosStart,rmin*cosEnd),rmax*sinStart);
0147     Set(pmax,rmax,rmax*sinEnd);
0148     break;
0149   case  9:                                 // start->end : 2->1
0150     Set(pmin,rmax*cosEnd,rmax*sinStart);
0151     Set(pmax,rmax,rmax);
0152     break;
0153   case 10:                                 // start->end : 2->2
0154     if (sinEnd < sinStart) break;
0155     Set(pmin,rmin*cosStart,rmax*sinStart);
0156     Set(pmax,rmax*cosEnd,rmin*sinEnd  );
0157     break;
0158   case 11:                                 // start->end : 2->3
0159     Set(pmin,-rmax,std::min(rmax*sinStart,rmax*sinEnd));
0160     Set(pmax,rmax,rmax);
0161     break;
0162   // start quadrant 3
0163   case 12:                                 // start->end : 3->0
0164     Set(pmin,rmax*cosStart,-rmax);
0165     Set(pmax,rmax,rmax*sinEnd);
0166     break;
0167   case 13:                                 // start->end : 3->1
0168     Set(pmin,std::min(rmax*cosStart,rmax*cosEnd),-rmax);
0169     Set(pmax,rmax,rmax);
0170     break;
0171   case 14:                                 // start->end : 3->2
0172     Set(pmin,rmax*cosStart,-rmax);
0173     Set(pmax,rmax*cosEnd,std::max(rmin*sinStart,rmin*sinEnd));
0174     break;
0175   case 15:                                 // start->end : 3->3
0176     if (sinEnd > sinStart) break;
0177     Set(pmin,rmax*cosStart,rmax*sinEnd);
0178     Set(pmax,rmin*cosEnd,rmin*sinStart);
0179     break;
0180   }
0181   return;
0182 }
0183
0184