Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/solids/CSG/src/G4Para.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /geometry/solids/CSG/src/G4Para.cc (Version 11.3.0) and /geometry/solids/CSG/src/G4Para.cc (Version 4.1.p1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                    <<   3 // * DISCLAIMER                                                       *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th <<   5 // * The following disclaimer summarizes all the specific disclaimers *
  6 // * the Geant4 Collaboration.  It is provided <<   6 // * of contributors to this software. The specific disclaimers,which *
  7 // * conditions of the Geant4 Software License <<   7 // * govern, are listed with their locations in:                      *
  8 // * LICENSE and available at  http://cern.ch/ <<   8 // *   http://cern.ch/geant4/license                                  *
  9 // * include a list of copyright holders.      << 
 10 // *                                                9 // *                                                                  *
 11 // * Neither the authors of this software syst     10 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     11 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     12 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     13 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  <<  14 // * use.                                                             *
 16 // * for the full disclaimer and the limitatio << 
 17 // *                                               15 // *                                                                  *
 18 // * This  code  implementation is the result  <<  16 // * This  code  implementation is the  intellectual property  of the *
 19 // * technical work of the GEANT4 collaboratio <<  17 // * GEANT4 collaboration.                                            *
 20 // * By using,  copying,  modifying or  distri <<  18 // * By copying,  distributing  or modifying the Program (or any work *
 21 // * any work based  on the software)  you  ag <<  19 // * based  on  the Program)  you indicate  your  acceptance of  this *
 22 // * use  in  resulting  scientific  publicati <<  20 // * statement, and all its terms.                                    *
 23 // * acceptance of all terms of the Geant4 Sof << 
 24 // *******************************************     21 // ********************************************************************
 25 //                                                 22 //
                                                   >>  23 //
                                                   >>  24 // $Id: G4Para.cc,v 1.11 2002/05/15 09:37:33 gcosmo Exp $
                                                   >>  25 // GEANT4 tag $Name: geant4-04-01 $
                                                   >>  26 //
                                                   >>  27 // class G4Para
                                                   >>  28 //
 26 // Implementation for G4Para class                 29 // Implementation for G4Para class
 27 //                                                 30 //
 28 // 21.03.95 P.Kent: Modified for `tolerant' ge <<  31 // History:
 29 // 31.10.96 V.Grichine: Modifications accordin <<  32 // 21.03.95 P.Kent Modified for `tolerant' geom
 30 // 28.04.05 V.Grichine: new SurfaceNormal acco <<  33 // 31.10.96 V.Grichine Modifications according G4Box/Tubs before to commit
 31 // 29.05.17 E.Tcherniaev: complete revision, s <<  34 // 18.11.99 V. Grichine , kUndef was added to ESide
 32 ////////////////////////////////////////////// <<  35 // 14.02.02 V. Grichine , bug fixed in Inside according to proposal of D. Wright
 33                                                    36 
 34 #include "G4Para.hh"                               37 #include "G4Para.hh"
 35                                                << 
 36 #if !defined(G4GEOM_USE_UPARA)                 << 
 37                                                << 
 38 #include "G4VoxelLimits.hh"                        38 #include "G4VoxelLimits.hh"
 39 #include "G4AffineTransform.hh"                    39 #include "G4AffineTransform.hh"
 40 #include "G4BoundingEnvelope.hh"               << 
 41 #include "Randomize.hh"                        << 
 42                                                    40 
 43 #include "G4VPVParameterisation.hh"                41 #include "G4VPVParameterisation.hh"
 44                                                    42 
 45 #include "G4VGraphicsScene.hh"                     43 #include "G4VGraphicsScene.hh"
                                                   >>  44 #include "G4Polyhedron.hh"
                                                   >>  45 #include "G4NURBS.hh"
                                                   >>  46 #include "G4NURBSbox.hh"
 46                                                    47 
 47 using namespace CLHEP;                         <<  48 // Private enum: Not for external use 
                                                   >>  49     
                                                   >>  50 enum ESide {kUndef,kPX,kMX,kPY,kMY,kPZ,kMZ};
 48                                                    51 
 49 ////////////////////////////////////////////// <<  52 // used internally for normal routine
 50 //                                             << 
 51 //  Constructor - set & check half widths      << 
 52                                                    53 
 53 G4Para::G4Para(const G4String& pName,          <<  54 enum ENSide {kNZ,kNX,kNY};
 54                      G4double pDx, G4double pD << 
 55                      G4double pAlpha, G4double << 
 56   : G4CSGSolid(pName), halfCarTolerance(0.5*kC << 
 57 {                                              << 
 58   SetAllParameters(pDx, pDy, pDz, pAlpha, pThe << 
 59   fRebuildPolyhedron = false;  // default valu << 
 60 }                                              << 
 61                                                    55 
 62 ////////////////////////////////////////////// <<  56 /////////////////////////////////////////////////////////////////////
 63 //                                                 57 //
 64 // Constructor - design of trapezoid based on  <<  58 // Constructor - check and set half-widths
 65                                                    59 
 66 G4Para::G4Para( const G4String& pName,         <<  60 void G4Para::SetAllParameters(G4double pDx, G4double pDy, G4double pDz, 
 67                 const G4ThreeVector pt[8] )    <<  61             G4double pAlpha, G4double pTheta, G4double pPhi)
 68   : G4CSGSolid(pName), halfCarTolerance(0.5*kC << 
 69 {                                                  62 {
 70   // Find dimensions and trigonometric values  <<  63   if (pDx>0&&pDy>0&&pDz>0)
 71   //                                           << 
 72   fDx = (pt[3].x() - pt[2].x())*0.5;           << 
 73   fDy = (pt[2].y() - pt[1].y())*0.5;           << 
 74   fDz = pt[7].z();                             << 
 75   CheckParameters(); // check dimensions       << 
 76                                                << 
 77   fTalpha = (pt[2].x() + pt[3].x() - pt[1].x() << 
 78   fTthetaCphi = (pt[4].x() + fDy*fTalpha + fDx << 
 79   fTthetaSphi = (pt[4].y() + fDy)/fDz;         << 
 80   MakePlanes();                                << 
 81                                                << 
 82   // Recompute vertices                        << 
 83   //                                           << 
 84   G4ThreeVector v[8];                          << 
 85   G4double DyTalpha = fDy*fTalpha;             << 
 86   G4double DzTthetaSphi = fDz*fTthetaSphi;     << 
 87   G4double DzTthetaCphi = fDz*fTthetaCphi;     << 
 88   v[0].set(-DzTthetaCphi-DyTalpha-fDx, -DzTthe << 
 89   v[1].set(-DzTthetaCphi-DyTalpha+fDx, -DzTthe << 
 90   v[2].set(-DzTthetaCphi+DyTalpha-fDx, -DzTthe << 
 91   v[3].set(-DzTthetaCphi+DyTalpha+fDx, -DzTthe << 
 92   v[4].set( DzTthetaCphi-DyTalpha-fDx,  DzTthe << 
 93   v[5].set( DzTthetaCphi-DyTalpha+fDx,  DzTthe << 
 94   v[6].set( DzTthetaCphi+DyTalpha-fDx,  DzTthe << 
 95   v[7].set( DzTthetaCphi+DyTalpha+fDx,  DzTthe << 
 96                                                << 
 97   // Compare with original vertices            << 
 98   //                                           << 
 99   for (G4int i=0; i<8; ++i)                    << 
100   {                                                64   {
101     G4double delx = std::abs(pt[i].x() - v[i]. <<  65      fDx=pDx;
102     G4double dely = std::abs(pt[i].y() - v[i]. <<  66      fDy=pDy;
103     G4double delz = std::abs(pt[i].z() - v[i]. <<  67      fDz=pDz;
104     G4double discrepancy = std::max(std::max(d <<  68      fTalpha=tan(pAlpha);
105     if (discrepancy > 0.1*kCarTolerance)       <<  69      fTthetaCphi=tan(pTheta)*cos(pPhi);
106     {                                          <<  70      fTthetaSphi=tan(pTheta)*sin(pPhi);
107       std::ostringstream message;              <<  71   }
108       G4long oldprc = message.precision(16);   <<  72   else
109       message << "Invalid vertice coordinates  <<  73   {
110               << "\nVertix #" << i << ", discr <<  74     G4Exception("Error in G4Para::SetAllParameters - Invalid Length Parameters");
111               << "\n  original   : " << pt[i]  << 
112               << "\n  recomputed : " << v[i];  << 
113       G4cout.precision(oldprc);                << 
114       G4Exception("G4Para::G4Para()", "GeomSol << 
115                   FatalException, message);    << 
116                                                << 
117     }                                          << 
118   }                                                75   }
119 }                                                  76 }
120                                                    77 
121 ////////////////////////////////////////////// <<  78 ///////////////////////////////////////////////////////////////////////////
122 //                                             << 
123 // Fake default constructor - sets only member << 
124 //                            for usage restri << 
125                                                << 
126 G4Para::G4Para( __void__& a )                  << 
127   : G4CSGSolid(a), halfCarTolerance(0.5*kCarTo << 
128 {                                              << 
129   SetAllParameters(1., 1., 1., 0., 0., 0.);    << 
130   fRebuildPolyhedron = false; // default value << 
131 }                                              << 
132                                                << 
133 ////////////////////////////////////////////// << 
134 //                                             << 
135 // Destructor                                  << 
136                                                << 
137 G4Para::~G4Para() = default;                   << 
138                                                << 
139 ////////////////////////////////////////////// << 
140 //                                             << 
141 // Copy constructor                            << 
142                                                << 
143 G4Para::G4Para(const G4Para& rhs)              << 
144   : G4CSGSolid(rhs), halfCarTolerance(rhs.half << 
145     fDx(rhs.fDx), fDy(rhs.fDy), fDz(rhs.fDz),  << 
146     fTthetaCphi(rhs.fTthetaCphi),fTthetaSphi(r << 
147 {                                              << 
148   for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs << 
149 }                                              << 
150                                                << 
151 ////////////////////////////////////////////// << 
152 //                                                 79 //
153 // Assignment operator                         << 
154                                                    80 
155 G4Para& G4Para::operator = (const G4Para& rhs) <<  81 G4Para::G4Para(const G4String& pName,G4double pDx, G4double pDy, G4double pDz,
                                                   >>  82            G4double pAlpha, G4double pTheta,G4double pPhi) : G4CSGSolid(pName)
156 {                                                  83 {
157    // Check assignment to self                 <<  84     if (pDx>0&&pDy>0&&pDz>0)
158    //                                          <<  85     {
159    if (this == &rhs)  { return *this; }        <<  86            SetAllParameters( pDx, pDy, pDz, pAlpha, pTheta, pPhi);
160                                                <<  87     }
161    // Copy base class data                     <<  88     else
162    //                                          <<  89     {
163    G4CSGSolid::operator=(rhs);                 <<  90        G4Exception("Error in G4Para::G4Para - Invalid Length Parameters");
164                                                <<  91     }
165    // Copy data                                << 
166    //                                          << 
167    halfCarTolerance = rhs.halfCarTolerance;    << 
168    fDx = rhs.fDx;                              << 
169    fDy = rhs.fDy;                              << 
170    fDz = rhs.fDz;                              << 
171    fTalpha = rhs.fTalpha;                      << 
172    fTthetaCphi = rhs.fTthetaCphi;              << 
173    fTthetaSphi = rhs.fTthetaSphi;              << 
174    for (G4int i=0; i<4; ++i) { fPlanes[i] = rh << 
175                                                << 
176    return *this;                               << 
177 }                                                  92 }
178                                                    93 
179 ////////////////////////////////////////////// <<  94 ////////////////////////////////////////////////////////////////////////
180 //                                                 95 //
181 // Set all parameters, as for constructor - se <<  96 // Constructor - Design of trapezoid based on 8 G4ThreeVector parameters, 
182                                                <<  97 // which are its vertices. Checking of planarity with preparation of 
183 void G4Para::SetAllParameters(G4double pDx, G4 <<  98 // fPlanes[] and than calculation of other members
184                               G4double pAlpha, << 
185 {                                              << 
186   // Reset data of the base class              << 
187   fCubicVolume = 0;                            << 
188   fSurfaceArea = 0;                            << 
189   fRebuildPolyhedron = true;                   << 
190                                                << 
191   // Set parameters                            << 
192   fDx = pDx;                                   << 
193   fDy = pDy;                                   << 
194   fDz = pDz;                                   << 
195   fTalpha = std::tan(pAlpha);                  << 
196   fTthetaCphi = std::tan(pTheta)*std::cos(pPhi << 
197   fTthetaSphi = std::tan(pTheta)*std::sin(pPhi << 
198                                                    99 
199   CheckParameters();                           << 
200   MakePlanes();                                << 
201 }                                              << 
202                                                   100 
203 ////////////////////////////////////////////// << 101 G4Para::G4Para( const G4String& pName,
204 //                                             << 102           const G4ThreeVector pt[8]) : G4CSGSolid(pName)
205 // Check dimensions                            << 
206                                                << 
207 void G4Para::CheckParameters()                 << 
208 {                                                 103 {
209   if (fDx < 2*kCarTolerance ||                 << 104     if (   pt[0].z()<0 && pt[0].z()==pt[1].z() && pt[0].z()==pt[2].z() && pt[0].z()==pt[3].z()
210       fDy < 2*kCarTolerance ||                 << 105   && pt[4].z()>0 && pt[4].z()==pt[5].z() && pt[4].z()==pt[6].z() && pt[4].z()==pt[7].z()
211       fDz < 2*kCarTolerance)                   << 106   && (pt[0].z()+pt[4].z())== 0
212   {                                            << 107   && pt[0].y()==pt[1].y() && pt[2].y()==pt[3].y()
213     std::ostringstream message;                << 108   && pt[4].y()==pt[5].y() && pt[6].y()==pt[7].y()
214     message << "Invalid (too small or negative << 109   && (pt[0].y()+pt[2].y()+pt[4].y()+pt[6].y())==0      )
215             << GetName()                       << 110      {
216             << "\n  X - " << fDx               << 111      fDz = (pt[7]).z() ;
217             << "\n  Y - " << fDy               << 112       
218             << "\n  Z - " << fDz;              << 113      fDy = ((pt[2]).y()-(pt[1]).y())*0.5 ;
219     G4Exception("G4Para::CheckParameters()", " << 114      fDx = ((pt[1]).x()-(pt[0]).x())*0.5 ;
220                 FatalException, message);      << 115      fDx = ((pt[3]).x()-(pt[2]).x())*0.5 ;
221   }                                            << 116      fTalpha = ((pt[2]).x()+(pt[3]).x()-(pt[1]).x()-(pt[0]).x())*0.25/fDy ;
222 }                                              << 117 
                                                   >> 118      //    fDy = ((pt[6]).y()-(pt[5]).y())*0.5 ;
                                                   >> 119      //    fDx = ((pt[5]).x()-(pt[4]).x())*0.5 ;
                                                   >> 120      //    fDx = ((pt[7]).x()-(pt[6]).x())*0.5 ;
                                                   >> 121      //    fTalpha = ((pt[6]).x()+(pt[7]).x()-(pt[5]).x()-(pt[4]).x())*0.25/fDy ;
223                                                   122 
224 ////////////////////////////////////////////// << 123      fTthetaCphi = ((pt[4]).x()+fDy*fTalpha+fDx)/fDz ;
225 //                                             << 124      fTthetaSphi = ((pt[4]).y()+fDy)/fDz ;
226 // Set side planes                             << 
227                                                << 
228 void G4Para::MakePlanes()                      << 
229 {                                              << 
230   G4ThreeVector vx(1, 0, 0);                   << 
231   G4ThreeVector vy(fTalpha, 1, 0);             << 
232   G4ThreeVector vz(fTthetaCphi, fTthetaSphi, 1 << 
233                                                << 
234   // Set -Y & +Y planes                        << 
235   //                                           << 
236   G4ThreeVector ynorm = (vx.cross(vz)).unit(); << 
237                                                << 
238   fPlanes[0].a = 0.;                           << 
239   fPlanes[0].b = ynorm.y();                    << 
240   fPlanes[0].c = ynorm.z();                    << 
241   fPlanes[0].d = fPlanes[0].b*fDy; // point (0 << 
242                                                << 
243   fPlanes[1].a =  0.;                          << 
244   fPlanes[1].b = -fPlanes[0].b;                << 
245   fPlanes[1].c = -fPlanes[0].c;                << 
246   fPlanes[1].d =  fPlanes[0].d;                << 
247                                                << 
248   // Set -X & +X planes                        << 
249   //                                           << 
250   G4ThreeVector xnorm = (vz.cross(vy)).unit(); << 
251                                                << 
252   fPlanes[2].a = xnorm.x();                    << 
253   fPlanes[2].b = xnorm.y();                    << 
254   fPlanes[2].c = xnorm.z();                    << 
255   fPlanes[2].d = fPlanes[2].a*fDx; // point (f << 
256                                                << 
257   fPlanes[3].a = -fPlanes[2].a;                << 
258   fPlanes[3].b = -fPlanes[2].b;                << 
259   fPlanes[3].c = -fPlanes[2].c;                << 
260   fPlanes[3].d =  fPlanes[2].d;                << 
261 }                                              << 
262                                                   125 
263 ////////////////////////////////////////////// << 126   }
264 //                                             << 127     else
265 // Get volume                                  << 128   {
                                                   >> 129       G4Exception("Error in G4Para::G4Para - Invalid vertice coordinates");
                                                   >> 130   }
266                                                   131 
267 G4double G4Para::GetCubicVolume()              << 132     
268 {                                              << 
269   // It is like G4Box, since para transformati << 
270   if (fCubicVolume == 0)                       << 
271   {                                            << 
272     fCubicVolume = 8*fDx*fDy*fDz;              << 
273   }                                            << 
274   return fCubicVolume;                         << 
275 }                                                 133 }
276                                                   134 
277 //////////////////////////////////////////////    135 //////////////////////////////////////////////////////////////////////////
278 //                                                136 //
279 // Get surface area                            << 
280                                                   137 
281 G4double G4Para::GetSurfaceArea()              << 138 G4Para::~G4Para()
282 {                                                 139 {
283   if(fSurfaceArea == 0)                        << 140    ;
284   {                                            << 
285     G4ThreeVector vx(fDx, 0, 0);               << 
286     G4ThreeVector vy(fDy*fTalpha, fDy, 0);     << 
287     G4ThreeVector vz(fDz*fTthetaCphi, fDz*fTth << 
288                                                << 
289     G4double sxy = fDx*fDy; // (vx.cross(vy)). << 
290     G4double sxz = (vx.cross(vz)).mag();       << 
291     G4double syz = (vy.cross(vz)).mag();       << 
292                                                << 
293     fSurfaceArea = 8*(sxy+sxz+syz);            << 
294   }                                            << 
295   return fSurfaceArea;                         << 
296 }                                                 141 }
297                                                   142 
298 //////////////////////////////////////////////    143 //////////////////////////////////////////////////////////////////////////
299 //                                                144 //
300 // Dispatch to parameterisation for replicatio    145 // Dispatch to parameterisation for replication mechanism dimension
301 // computation & modification                  << 146 // computation & modification.
302                                                   147 
303 void G4Para::ComputeDimensions(      G4VPVPara << 148  void G4Para::ComputeDimensions(G4VPVParameterisation* p,
304                                 const G4int n,    149                                 const G4int n,
305                                 const G4VPhysi << 150                                 const G4VPhysicalVolume* pRep)
306 {                                                 151 {
307   p->ComputeDimensions(*this,n,pRep);          << 152     p->ComputeDimensions(*this,n,pRep);
308 }                                                 153 }
309                                                   154 
310 ////////////////////////////////////////////// << 
311 //                                             << 
312 // Get bounding box                            << 
313                                                << 
314 void G4Para::BoundingLimits(G4ThreeVector& pMi << 
315 {                                              << 
316   G4double dz = GetZHalfLength();              << 
317   G4double dx = GetXHalfLength();              << 
318   G4double dy = GetYHalfLength();              << 
319                                                << 
320   G4double x0 = dz*fTthetaCphi;                << 
321   G4double x1 = dy*GetTanAlpha();              << 
322   G4double xmin =                              << 
323     std::min(                                  << 
324     std::min(                                  << 
325     std::min(-x0-x1-dx,-x0+x1-dx),x0-x1-dx),x0 << 
326   G4double xmax =                              << 
327     std::max(                                  << 
328     std::max(                                  << 
329     std::max(-x0-x1+dx,-x0+x1+dx),x0-x1+dx),x0 << 
330                                                << 
331   G4double y0 = dz*fTthetaSphi;                << 
332   G4double ymin = std::min(-y0-dy,y0-dy);      << 
333   G4double ymax = std::max(-y0+dy,y0+dy);      << 
334                                                << 
335   pMin.set(xmin,ymin,-dz);                     << 
336   pMax.set(xmax,ymax, dz);                     << 
337                                                << 
338   // Check correctness of the bounding box     << 
339   //                                           << 
340   if (pMin.x() >= pMax.x() || pMin.y() >= pMax << 
341   {                                            << 
342     std::ostringstream message;                << 
343     message << "Bad bounding box (min >= max)  << 
344             << GetName() << " !"               << 
345             << "\npMin = " << pMin             << 
346             << "\npMax = " << pMax;            << 
347     G4Exception("G4Para::BoundingLimits()", "G << 
348                 JustWarning, message);         << 
349     DumpInfo();                                << 
350   }                                            << 
351 }                                              << 
352                                                   155 
353 ////////////////////////////////////////////// << 156 //////////////////////////////////////////////////////////////
354 //                                                157 //
355 // Calculate extent under transform and specif    158 // Calculate extent under transform and specified limit
356                                                   159 
357 G4bool G4Para::CalculateExtent( const EAxis pA << 160 G4bool G4Para::CalculateExtent(const EAxis pAxis,
358                                 const G4VoxelL << 161                                const G4VoxelLimits& pVoxelLimit,
359                                 const G4Affine << 162                                const G4AffineTransform& pTransform,
360                                      G4double& << 163                                G4double& pMin, G4double& pMax) const
361 {                                              << 164 {
362   G4ThreeVector bmin, bmax;                    << 165     G4bool flag;
363   G4bool exist;                                << 166 
364                                                << 167     if (!pTransform.IsRotated())
365   // Check bounding box (bbox)                 << 168    {  
366   //                                           << 169                       // Special case handling for unrotated trapezoids
367   BoundingLimits(bmin,bmax);                   << 170                       // Compute z/x/y/ mins and maxs respecting limits, with early returns
368   G4BoundingEnvelope bbox(bmin,bmax);          << 171                       // if outside limits. Then switch() on pAxis
369 #ifdef G4BBOX_EXTENT                           << 172     G4int i ; 
370   return bbox.CalculateExtent(pAxis,pVoxelLimi << 173           G4double xoffset,xMin,xMax;
371 #endif                                         << 174           G4double yoffset,yMin,yMax;
372   if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox << 175           G4double zoffset,zMin,zMax;
373   {                                            << 176     G4double temp[8] ;          // some points for intersection with zMin/zMax
374     return exist = pMin < pMax;                << 177     
375   }                                            << 178           xoffset=pTransform.NetTranslation().x();      
376                                                << 179           yoffset=pTransform.NetTranslation().y();
377   // Set bounding envelope (benv) and calculat << 180           zoffset=pTransform.NetTranslation().z();
378   //                                           << 181  
379   G4double dz = GetZHalfLength();              << 182           G4ThreeVector pt[8];   // vertices after translation
380   G4double dx = GetXHalfLength();              << 183           pt[0]=G4ThreeVector(xoffset-fDz*fTthetaCphi-fDy*fTalpha-fDx,
381   G4double dy = GetYHalfLength();              << 184                         yoffset-fDz*fTthetaSphi-fDy,zoffset-fDz);
382                                                << 185           pt[1]=G4ThreeVector(xoffset-fDz*fTthetaCphi-fDy*fTalpha+fDx,
383   G4double x0 = dz*fTthetaCphi;                << 186                   yoffset-fDz*fTthetaSphi-fDy,zoffset-fDz);
384   G4double x1 = dy*GetTanAlpha();              << 187           pt[2]=G4ThreeVector(xoffset-fDz*fTthetaCphi+fDy*fTalpha-fDx,
385   G4double y0 = dz*fTthetaSphi;                << 188                   yoffset-fDz*fTthetaSphi+fDy,zoffset-fDz);
386                                                << 189           pt[3]=G4ThreeVector(xoffset-fDz*fTthetaCphi+fDy*fTalpha+fDx,
387   G4ThreeVectorList baseA(4), baseB(4);        << 190                   yoffset-fDz*fTthetaSphi+fDy,zoffset-fDz);
388   baseA[0].set(-x0-x1-dx,-y0-dy,-dz);          << 191           pt[4]=G4ThreeVector(xoffset+fDz*fTthetaCphi-fDy*fTalpha-fDx,
389   baseA[1].set(-x0-x1+dx,-y0-dy,-dz);          << 192                   yoffset+fDz*fTthetaSphi-fDy,zoffset+fDz);
390   baseA[2].set(-x0+x1+dx,-y0+dy,-dz);          << 193           pt[5]=G4ThreeVector(xoffset+fDz*fTthetaCphi-fDy*fTalpha+fDx,
391   baseA[3].set(-x0+x1-dx,-y0+dy,-dz);          << 194                   yoffset+fDz*fTthetaSphi-fDy,zoffset+fDz);
392                                                << 195           pt[6]=G4ThreeVector(xoffset+fDz*fTthetaCphi+fDy*fTalpha-fDx,
393   baseB[0].set(+x0-x1-dx, y0-dy, dz);          << 196                   yoffset+fDz*fTthetaSphi+fDy,zoffset+fDz);
394   baseB[1].set(+x0-x1+dx, y0-dy, dz);          << 197           pt[7]=G4ThreeVector(xoffset+fDz*fTthetaCphi+fDy*fTalpha+fDx,
395   baseB[2].set(+x0+x1+dx, y0+dy, dz);          << 198                   yoffset+fDz*fTthetaSphi+fDy,zoffset+fDz);
396   baseB[3].set(+x0+x1-dx, y0+dy, dz);          << 199       zMin=zoffset-fDz;
397                                                << 200       zMax=zoffset+fDz;
398   std::vector<const G4ThreeVectorList *> polyg << 201       if (pVoxelLimit.IsZLimited())
399   polygons[0] = &baseA;                        << 202     {
400   polygons[1] = &baseB;                        << 203         if (zMin>pVoxelLimit.GetMaxZExtent()+kCarTolerance
401                                                << 204       ||zMax<pVoxelLimit.GetMinZExtent()-kCarTolerance)
402   G4BoundingEnvelope benv(bmin,bmax,polygons); << 205       {
403   exist = benv.CalculateExtent(pAxis,pVoxelLim << 206           return false;
404   return exist;                                << 207       }
405 }                                              << 208         else
406                                                << 209       {
407 ////////////////////////////////////////////// << 210           if (zMin<pVoxelLimit.GetMinZExtent())
408 //                                             << 211         {
409 // Determine where is point p, inside/on_surfa << 212             zMin=pVoxelLimit.GetMinZExtent();
410 //                                             << 213         }
411                                                << 214           if (zMax>pVoxelLimit.GetMaxZExtent())
412 EInside G4Para::Inside( const G4ThreeVector& p << 215         {
413 {                                              << 216             zMax=pVoxelLimit.GetMaxZExtent();
414   G4double xx = fPlanes[2].a*p.x()+fPlanes[2]. << 217         }
415   G4double dx = std::abs(xx) + fPlanes[2].d;   << 218       }
                                                   >> 219     }
                                                   >> 220 
                                                   >> 221             temp[0] = pt[0].y()+(pt[4].y()-pt[0].y())*(zMin-pt[0].z())/(pt[4].z()-pt[0].z()) ;
                                                   >> 222             temp[1] = pt[0].y()+(pt[4].y()-pt[0].y())*(zMax-pt[0].z())/(pt[4].z()-pt[0].z()) ;
                                                   >> 223       temp[2] = pt[2].y()+(pt[6].y()-pt[2].y())*(zMin-pt[2].z())/(pt[6].z()-pt[2].z()) ;
                                                   >> 224       temp[3] = pt[2].y()+(pt[6].y()-pt[2].y())*(zMax-pt[2].z())/(pt[6].z()-pt[2].z()) ;        
                                                   >> 225       yMax = yoffset - fabs(fDz*fTthetaSphi) - fDy - fDy ;
                                                   >> 226       yMin = -yMax ;
                                                   >> 227       for(i=0;i<4;i++)
                                                   >> 228       {
                                                   >> 229          if(temp[i] > yMax) yMax = temp[i] ;
                                                   >> 230          if(temp[i] < yMin) yMin = temp[i] ;
                                                   >> 231       }
                                                   >> 232       
                                                   >> 233       if (pVoxelLimit.IsYLimited())
                                                   >> 234     {
                                                   >> 235         if (yMin>pVoxelLimit.GetMaxYExtent()+kCarTolerance
                                                   >> 236       ||yMax<pVoxelLimit.GetMinYExtent()-kCarTolerance)
                                                   >> 237       {
                                                   >> 238           return false;
                                                   >> 239       }
                                                   >> 240         else
                                                   >> 241       {
                                                   >> 242           if (yMin<pVoxelLimit.GetMinYExtent())
                                                   >> 243         {
                                                   >> 244             yMin=pVoxelLimit.GetMinYExtent();
                                                   >> 245         }
                                                   >> 246           if (yMax>pVoxelLimit.GetMaxYExtent())
                                                   >> 247         {
                                                   >> 248             yMax=pVoxelLimit.GetMaxYExtent();
                                                   >> 249         }
                                                   >> 250       }
                                                   >> 251     }
                                                   >> 252 
                                                   >> 253             temp[0] = pt[0].x()+(pt[4].x()-pt[0].x())*(zMin-pt[0].z())/(pt[4].z()-pt[0].z()) ;
                                                   >> 254             temp[1] = pt[0].x()+(pt[4].x()-pt[0].x())*(zMax-pt[0].z())/(pt[4].z()-pt[0].z()) ;
                                                   >> 255       temp[2] = pt[2].x()+(pt[6].x()-pt[2].x())*(zMin-pt[2].z())/(pt[6].z()-pt[2].z()) ;
                                                   >> 256       temp[3] = pt[2].x()+(pt[6].x()-pt[2].x())*(zMax-pt[2].z())/(pt[6].z()-pt[2].z()) ;
                                                   >> 257             temp[4] = pt[3].x()+(pt[7].x()-pt[3].x())*(zMin-pt[3].z())/(pt[7].z()-pt[3].z()) ;
                                                   >> 258             temp[5] = pt[3].x()+(pt[7].x()-pt[3].x())*(zMax-pt[3].z())/(pt[7].z()-pt[3].z()) ;
                                                   >> 259       temp[6] = pt[1].x()+(pt[5].x()-pt[1].x())*(zMin-pt[1].z())/(pt[5].z()-pt[1].z()) ;
                                                   >> 260       temp[7] = pt[1].x()+(pt[5].x()-pt[1].x())*(zMax-pt[1].z())/(pt[5].z()-pt[1].z()) ;
                                                   >> 261       
                                                   >> 262       xMax = xoffset - fabs(fDz*fTthetaCphi) - fDx - fDx -fDx - fDx;
                                                   >> 263       xMin = -xMax ;
                                                   >> 264       for(i=0;i<8;i++)
                                                   >> 265       {
                                                   >> 266          if(temp[i] > xMax) xMax = temp[i] ;
                                                   >> 267          if(temp[i] < xMin) xMin = temp[i] ;
                                                   >> 268       }
                                                   >> 269                                             // xMax/Min = f(yMax/Min) ?
                                                   >> 270       if (pVoxelLimit.IsXLimited())
                                                   >> 271     {
                                                   >> 272         if (xMin>pVoxelLimit.GetMaxXExtent()+kCarTolerance
                                                   >> 273       ||xMax<pVoxelLimit.GetMinXExtent()-kCarTolerance)
                                                   >> 274       {
                                                   >> 275           return false;
                                                   >> 276       }
                                                   >> 277         else
                                                   >> 278       {
                                                   >> 279           if (xMin<pVoxelLimit.GetMinXExtent())
                                                   >> 280         {
                                                   >> 281             xMin=pVoxelLimit.GetMinXExtent();
                                                   >> 282         }
                                                   >> 283           if (xMax>pVoxelLimit.GetMaxXExtent())
                                                   >> 284         {
                                                   >> 285             xMax=pVoxelLimit.GetMaxXExtent();
                                                   >> 286         }
                                                   >> 287       }
                                                   >> 288     }
                                                   >> 289 
                                                   >> 290       switch (pAxis)
                                                   >> 291     {
                                                   >> 292     case kXAxis:
                                                   >> 293         pMin=xMin;
                                                   >> 294         pMax=xMax;
                                                   >> 295         break;
                                                   >> 296     case kYAxis:
                                                   >> 297         pMin=yMin;
                                                   >> 298         pMax=yMax;
                                                   >> 299         break;
                                                   >> 300     case kZAxis:
                                                   >> 301         pMin=zMin;
                                                   >> 302         pMax=zMax;
                                                   >> 303         break;
                                                   >> 304     default:
                                                   >> 305         break;
                                                   >> 306     }
416                                                   307 
417   G4double yy = fPlanes[0].b*p.y()+fPlanes[0]. << 308       pMin-=kCarTolerance;
418   G4double dy = std::abs(yy) + fPlanes[0].d;   << 309       pMax+=kCarTolerance;
419   G4double dxy = std::max(dx,dy);              << 
420                                                   310 
421   G4double dz = std::abs(p.z())-fDz;           << 311       flag = true;
422   G4double dist = std::max(dxy,dz);            << 312    }
                                                   >> 313     else
                                                   >> 314    {
                                                   >> 315 // General rotated case - create and clip mesh to boundaries
423                                                   316 
424   if (dist > halfCarTolerance) return kOutside << 317       G4bool existsAfterClip=false;
425   return (dist > -halfCarTolerance) ? kSurface << 318       G4ThreeVectorList *vertices;
426 }                                              << 
427                                                   319 
428 ////////////////////////////////////////////// << 320       pMin=+kInfinity;
429 //                                             << 321       pMax=-kInfinity;
430 // Determine side where point is, and return c << 322 // Calculate rotated vertex coordinates
431                                                << 323 
432 G4ThreeVector G4Para::SurfaceNormal( const G4T << 324       vertices=CreateRotatedVertices(pTransform);
433 {                                              << 325       ClipCrossSection(vertices,0,pVoxelLimit,pAxis,pMin,pMax);
434   G4int nsurf = 0; // number of surfaces where << 326       ClipCrossSection(vertices,4,pVoxelLimit,pAxis,pMin,pMax);
                                                   >> 327       ClipBetweenSections(vertices,0,pVoxelLimit,pAxis,pMin,pMax);
                                                   >> 328       
                                                   >> 329       if (pMin!=kInfinity||pMax!=-kInfinity)
                                                   >> 330     {
                                                   >> 331         existsAfterClip=true;
                                                   >> 332         
                                                   >> 333 // Add 2*tolerance to avoid precision troubles
                                                   >> 334         pMin-=kCarTolerance;
                                                   >> 335         pMax+=kCarTolerance;
                                                   >> 336         
                                                   >> 337     }
                                                   >> 338       else
                                                   >> 339     {
                                                   >> 340 // Check for case where completely enveloping clipping volume
                                                   >> 341 // If point inside then we are confident that the solid completely
                                                   >> 342 // envelopes the clipping volume. Hence set min/max extents according
                                                   >> 343 // to clipping volume extents along the specified axis.
                                                   >> 344        
                                                   >> 345            G4ThreeVector clipCentre(
                                                   >> 346       (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
                                                   >> 347       (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
                                                   >> 348       (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5);
                                                   >> 349         
                                                   >> 350         if (Inside(pTransform.Inverse().TransformPoint(clipCentre))!=kOutside)
                                                   >> 351       {
                                                   >> 352           existsAfterClip=true;
                                                   >> 353           pMin=pVoxelLimit.GetMinExtent(pAxis);
                                                   >> 354           pMax=pVoxelLimit.GetMaxExtent(pAxis);
                                                   >> 355       }
                                                   >> 356     }
                                                   >> 357       delete vertices ;          //  'new' in the function called
                                                   >> 358       flag = existsAfterClip ;
                                                   >> 359    }
                                                   >> 360    return flag;
                                                   >> 361 }
                                                   >> 362 
                                                   >> 363 /////////////////////////////////////////////////////////////////////////////
                                                   >> 364 //
                                                   >> 365 // Check in p is inside/on surface/outside solid
                                                   >> 366 
                                                   >> 367 EInside G4Para::Inside(const G4ThreeVector& p) const
                                                   >> 368 {
                                                   >> 369   G4double xt, yt, yt1;
                                                   >> 370   EInside  in = kOutside;
                                                   >> 371 
                                                   >> 372   yt1 = p.y() - fTthetaSphi*p.z();
                                                   >> 373   yt  = fabs(yt1) ;
435                                                   374 
436   // Check Z faces                             << 375   // xt = fabs( p.x() - fTthetaCphi*p.z() - fTalpha*yt );
437   //                                           << 
438   G4double nz = 0;                             << 
439   G4double dz = std::abs(p.z()) - fDz;         << 
440   if (std::abs(dz) <= halfCarTolerance)        << 
441   {                                            << 
442     nz = (p.z() < 0) ? -1 : 1;                 << 
443     ++nsurf;                                   << 
444   }                                            << 
445                                                   376 
446   // Check Y faces                             << 377   xt = fabs( p.x() - fTthetaCphi*p.z() - fTalpha*yt1 );
447   //                                           << 
448   G4double ny = 0;                             << 
449   G4double yy = fPlanes[0].b*p.y()+fPlanes[0]. << 
450   if (std::abs(fPlanes[0].d + yy) <= halfCarTo << 
451   {                                            << 
452     ny  = fPlanes[0].b;                        << 
453     nz += fPlanes[0].c;                        << 
454     ++nsurf;                                   << 
455   }                                            << 
456   else if (std::abs(fPlanes[1].d - yy) <= half << 
457   {                                            << 
458     ny  = fPlanes[1].b;                        << 
459     nz += fPlanes[1].c;                        << 
460     ++nsurf;                                   << 
461   }                                            << 
462                                                   378 
463   // Check X faces                             << 379   if ( fabs( p.z() ) <= fDz - kCarTolerance*0.5)
464   //                                           << 
465   G4double nx = 0;                             << 
466   G4double xx = fPlanes[2].a*p.x()+fPlanes[2]. << 
467   if (std::abs(fPlanes[2].d + xx) <= halfCarTo << 
468   {                                               380   {
469     nx  = fPlanes[2].a;                        << 381     if (yt <= fDy - kCarTolerance*0.5)
470     ny += fPlanes[2].b;                        << 382     {
471     nz += fPlanes[2].c;                        << 383       if      ( xt <= fDx - kCarTolerance*0.5 ) in = kInside;
472     ++nsurf;                                   << 384       else if ( xt <= fDx + kCarTolerance*0.5 ) in = kSurface;
                                                   >> 385     }
                                                   >> 386     else if ( yt <= fDy + kCarTolerance*0.5)
                                                   >> 387     {
                                                   >> 388       if ( xt <= fDx + kCarTolerance*0.5 ) in = kSurface; 
                                                   >> 389     }
473   }                                               390   }
474   else if (std::abs(fPlanes[3].d - xx) <= half << 391   else  if ( fabs(p.z()) <= fDz + kCarTolerance*0.5 )
475   {                                               392   {
476     nx  = fPlanes[3].a;                        << 
477     ny += fPlanes[3].b;                        << 
478     nz += fPlanes[3].c;                        << 
479     ++nsurf;                                   << 
480   }                                            << 
481                                                   393 
482   // Return normal                             << 394     if ( yt <= fDy + kCarTolerance*0.5)
483   //                                           << 395     {
484   if (nsurf == 1)      return {nx,ny,nz};      << 396       if ( xt <= fDx + kCarTolerance*0.5 ) in = kSurface; 
485   else if (nsurf != 0) return G4ThreeVector(nx << 397     }
486   else                                         << 
487   {                                            << 
488     // Point is not on the surface             << 
489     //                                         << 
490 #ifdef G4CSGDEBUG                              << 
491     std::ostringstream message;                << 
492     G4int oldprc = message.precision(16);      << 
493     message << "Point p is not on surface (!?) << 
494             << GetName() << G4endl;            << 
495     message << "Position:\n";                  << 
496     message << "   p.x() = " << p.x()/mm << "  << 
497     message << "   p.y() = " << p.y()/mm << "  << 
498     message << "   p.z() = " << p.z()/mm << "  << 
499     G4cout.precision(oldprc) ;                 << 
500     G4Exception("G4Para::SurfaceNormal(p)", "G << 
501                 JustWarning, message );        << 
502     DumpInfo();                                << 
503 #endif                                         << 
504     return ApproxSurfaceNormal(p);             << 
505   }                                               398   }
                                                   >> 399   return in;
506 }                                                 400 }
507                                                   401 
508 ////////////////////////////////////////////// << 402 ///////////////////////////////////////////////////////////////////////////
509 //                                                403 //
510 // Algorithm for SurfaceNormal() following the << 404 // Calculate side nearest to p, and return normal
511 // for points not on the surface               << 405 // If 2+ sides equidistant, first side's normal returned (arbitrarily)
512                                                   406 
513 G4ThreeVector G4Para::ApproxSurfaceNormal( con << 407 G4ThreeVector G4Para::SurfaceNormal( const G4ThreeVector& p) const
514 {                                                 408 {
515   G4double dist = -DBL_MAX;                    << 409     ENSide  side;
516   G4int iside = 0;                             << 410     G4ThreeVector norm;
517   for (G4int i=0; i<4; ++i)                    << 411     G4double distx,disty,distz;
518   {                                            << 412     G4double newpx,newpy,xshift;
519     G4double d = fPlanes[i].a*p.x() +          << 413     G4double calpha,salpha; // Sin/Cos(alpha) - needed to recalc G4Parameter 
520                  fPlanes[i].b*p.y() +          << 414     G4double tntheta,cosntheta; // tan and cos of normal's theta component
521                  fPlanes[i].c*p.z() + fPlanes[ << 415     G4double ycomp;
522     if (d > dist) { dist = d; iside = i; }     << 
523   }                                            << 
524                                                   416 
525   G4double distz = std::abs(p.z()) - fDz;      << 417     newpx=p.x()-fTthetaCphi*p.z();
526   if (dist > distz)                            << 418     newpy=p.y()-fTthetaSphi*p.z();
527     return { fPlanes[iside].a, fPlanes[iside]. << 
528   else                                         << 
529     return { 0, 0, (G4double)((p.z() < 0) ? -1 << 
530 }                                              << 
531                                                   419 
532 ////////////////////////////////////////////// << 420     calpha=1/sqrt(1+fTalpha*fTalpha);
                                                   >> 421     if (fTalpha)
                                                   >> 422   {
                                                   >> 423       salpha=-calpha/fTalpha; // NOTE: actually use MINUS sin(alpha)
                                                   >> 424   }
                                                   >> 425     else
                                                   >> 426   {
                                                   >> 427       salpha=0;
                                                   >> 428   }
                                                   >> 429 
                                                   >> 430     xshift=newpx*calpha+newpy*salpha;
                                                   >> 431 
                                                   >> 432     distx=fabs(fabs(xshift)-fDx*calpha);
                                                   >> 433     disty=fabs(fabs(newpy)-fDy);
                                                   >> 434     distz=fabs(fabs(p.z())-fDz);
                                                   >> 435     
                                                   >> 436     if (distx<disty)
                                                   >> 437   {
                                                   >> 438       if (distx<distz) side=kNX;
                                                   >> 439       else side=kNZ;
                                                   >> 440   }
                                                   >> 441     else
                                                   >> 442   {
                                                   >> 443       if (disty<distz) side=kNY;
                                                   >> 444       else side=kNZ;
                                                   >> 445   }
                                                   >> 446 
                                                   >> 447     switch (side)
                                                   >> 448   {
                                                   >> 449   case kNX:
                                                   >> 450       tntheta=fTthetaCphi*calpha+fTthetaSphi*salpha;
                                                   >> 451       if (xshift<0)
                                                   >> 452     {
                                                   >> 453         cosntheta=-1/sqrt(1+tntheta*tntheta);
                                                   >> 454     }
                                                   >> 455       else
                                                   >> 456     {
                                                   >> 457         cosntheta=1/sqrt(1+tntheta*tntheta);
                                                   >> 458     }
                                                   >> 459       norm=G4ThreeVector(calpha*cosntheta,salpha*cosntheta,-tntheta*cosntheta);
                                                   >> 460       break;
                                                   >> 461   case kNY:
                                                   >> 462       if (newpy<0)
                                                   >> 463     {
                                                   >> 464         ycomp=-1/sqrt(1+fTthetaSphi*fTthetaSphi);
                                                   >> 465     }
                                                   >> 466       else
                                                   >> 467     {
                                                   >> 468         ycomp=1/sqrt(1+fTthetaSphi*fTthetaSphi);
                                                   >> 469     }
                                                   >> 470       norm=G4ThreeVector(0,ycomp,-fTthetaSphi*ycomp);
                                                   >> 471       break;
                                                   >> 472   case kNZ:
                                                   >> 473 // Closest to Z
                                                   >> 474       if (p.z()>=0)
                                                   >> 475     {
                                                   >> 476         norm=G4ThreeVector(0,0,1);
                                                   >> 477     }
                                                   >> 478       else
                                                   >> 479     {
                                                   >> 480         norm=G4ThreeVector(0,0,-1);
                                                   >> 481     }
                                                   >> 482       break;
                                                   >> 483   }
                                                   >> 484     return norm;
                                                   >> 485 }
                                                   >> 486 
                                                   >> 487 //////////////////////////////////////////////////////////////////////////////
                                                   >> 488 //
                                                   >> 489 // Calculate distance to shape from outside - return kInfinity if no intersection
                                                   >> 490 //
                                                   >> 491 // ALGORITHM:
                                                   >> 492 // For each component, calculate pair of minimum and maximum intersection
                                                   >> 493 // values for which the particle is in the extent of the shape
                                                   >> 494 // - The smallest (MAX minimum) allowed distance of the pairs is intersect
                                                   >> 495 // - Z plane intersectin uses tolerance
                                                   >> 496 // - XZ YZ planes use logic & *SLIGHTLY INCORRECT* tolerance
                                                   >> 497 //   (this saves at least 1 sqrt, 1 multiply and 1 divide... in applicable
                                                   >> 498 //    cases)
                                                   >> 499 // - Note: XZ and YZ planes each divide space into four regions,
                                                   >> 500 //   characterised by ss1 ss2
                                                   >> 501 
                                                   >> 502 G4double G4Para::DistanceToIn(const G4ThreeVector& p,const G4ThreeVector& v) const
                                                   >> 503 {
                                                   >> 504     G4double snxt;    // snxt = default return value
                                                   >> 505     G4double smin,smax;
                                                   >> 506     G4double tmin,tmax;
                                                   >> 507     G4double yt,vy,xt,vx;
                                                   >> 508     G4double max;
                                                   >> 509 //
                                                   >> 510 // Z Intersection range
                                                   >> 511 //
                                                   >> 512     if (v.z()>0)
                                                   >> 513   {
                                                   >> 514       max=fDz-p.z();
                                                   >> 515       if (max>kCarTolerance*0.5)
                                                   >> 516     {
                                                   >> 517         smax=max/v.z();
                                                   >> 518         smin=(-fDz-p.z())/v.z();
                                                   >> 519     }
                                                   >> 520       else
                                                   >> 521     {
                                                   >> 522         return snxt=kInfinity;
                                                   >> 523     }
                                                   >> 524   }
                                                   >> 525     else if (v.z()<0)
                                                   >> 526   {
                                                   >> 527       max=-fDz-p.z();
                                                   >> 528       if (max<-kCarTolerance*0.5)
                                                   >> 529     {
                                                   >> 530         smax=max/v.z();
                                                   >> 531         smin=(fDz-p.z())/v.z();
                                                   >> 532     }
                                                   >> 533       else
                                                   >> 534     {
                                                   >> 535         return snxt=kInfinity;
                                                   >> 536     }
                                                   >> 537   }
                                                   >> 538     else
                                                   >> 539   {
                                                   >> 540       if (fabs(p.z())<=fDz) // Inside
                                                   >> 541     {
                                                   >> 542         smin=0;
                                                   >> 543         smax=kInfinity;
                                                   >> 544     }
                                                   >> 545       else
                                                   >> 546     {
                                                   >> 547         return snxt=kInfinity;
                                                   >> 548     }
                                                   >> 549   }
                                                   >> 550     
533 //                                                551 //
534 // Calculate distance to shape from outside    << 552 // Y G4Parallel planes intersection
535 //  - return kInfinity if no intersection      << 553 //
536                                                << 554     yt=p.y()-fTthetaSphi*p.z();
537 G4double G4Para::DistanceToIn(const G4ThreeVec << 555     vy=v.y()-fTthetaSphi*v.z();
538                               const G4ThreeVec << 
539 {                                              << 
540   // Z intersections                           << 
541   //                                           << 
542   if ((std::abs(p.z()) - fDz) >= -halfCarToler << 
543     return kInfinity;                          << 
544   G4double invz = (-v.z() == 0) ? DBL_MAX : -1 << 
545   G4double dz = (invz < 0) ? fDz : -fDz;       << 
546   G4double tzmin = (p.z() + dz)*invz;          << 
547   G4double tzmax = (p.z() - dz)*invz;          << 
548                                                << 
549   // Y intersections                           << 
550   //                                           << 
551   G4double tmin0 = tzmin, tmax0 = tzmax;       << 
552   G4double cos0 = fPlanes[0].b*v.y() + fPlanes << 
553   G4double disy = fPlanes[0].b*p.y() + fPlanes << 
554   G4double dis0 = fPlanes[0].d + disy;         << 
555   if (dis0 >= -halfCarTolerance)               << 
556   {                                            << 
557     if (cos0 >= 0) return kInfinity;           << 
558     G4double tmp  = -dis0/cos0;                << 
559     if (tmin0 < tmp) tmin0 = tmp;              << 
560   }                                            << 
561   else if (cos0 > 0)                           << 
562   {                                            << 
563     G4double tmp  = -dis0/cos0;                << 
564     if (tmax0 > tmp) tmax0 = tmp;              << 
565   }                                            << 
566                                                << 
567   G4double tmin1 = tmin0, tmax1 = tmax0;       << 
568   G4double cos1 = -cos0;                       << 
569   G4double dis1 = fPlanes[1].d - disy;         << 
570   if (dis1 >= -halfCarTolerance)               << 
571   {                                            << 
572     if (cos1 >= 0) return kInfinity;           << 
573     G4double tmp  = -dis1/cos1;                << 
574     if (tmin1 < tmp) tmin1 = tmp;              << 
575   }                                            << 
576   else if (cos1 > 0)                           << 
577   {                                            << 
578     G4double tmp  = -dis1/cos1;                << 
579     if (tmax1 > tmp) tmax1 = tmp;              << 
580   }                                            << 
581                                                << 
582   // X intersections                           << 
583   //                                           << 
584   G4double tmin2 = tmin1, tmax2 = tmax1;       << 
585   G4double cos2 = fPlanes[2].a*v.x() + fPlanes << 
586   G4double disx = fPlanes[2].a*p.x() + fPlanes << 
587   G4double dis2 = fPlanes[2].d + disx;         << 
588   if (dis2 >= -halfCarTolerance)               << 
589   {                                            << 
590     if (cos2 >= 0) return kInfinity;           << 
591     G4double tmp  = -dis2/cos2;                << 
592     if (tmin2 < tmp) tmin2 = tmp;              << 
593   }                                            << 
594   else if (cos2 > 0)                           << 
595   {                                            << 
596     G4double tmp  = -dis2/cos2;                << 
597     if (tmax2 > tmp) tmax2 = tmp;              << 
598   }                                            << 
599                                                << 
600   G4double tmin3 = tmin2, tmax3 = tmax2;       << 
601   G4double cos3 = -cos2;                       << 
602   G4double dis3 = fPlanes[3].d - disx;         << 
603   if (dis3 >= -halfCarTolerance)               << 
604   {                                            << 
605     if (cos3 >= 0) return kInfinity;           << 
606     G4double tmp  = -dis3/cos3;                << 
607     if (tmin3 < tmp) tmin3 = tmp;              << 
608   }                                            << 
609   else if (cos3 > 0)                           << 
610   {                                            << 
611     G4double tmp  = -dis3/cos3;                << 
612     if (tmax3 > tmp) tmax3 = tmp;              << 
613   }                                            << 
614                                                   556 
615   // Find distance                             << 557     if (vy>0)
616   //                                           << 558   {
617   G4double tmin = tmin3, tmax = tmax3;         << 559       max=fDy-yt;
618   if (tmax <= tmin + halfCarTolerance) return  << 560       if (max>kCarTolerance*0.5)
619   return (tmin < halfCarTolerance ) ? 0. : tmi << 561     {
                                                   >> 562         tmax=max/vy;
                                                   >> 563         tmin=(-fDy-yt)/vy;
                                                   >> 564     }
                                                   >> 565       else
                                                   >> 566     {
                                                   >> 567         return snxt=kInfinity;
                                                   >> 568     }
                                                   >> 569   }
                                                   >> 570     else if (vy<0)
                                                   >> 571   {
                                                   >> 572       max=-fDy-yt;
                                                   >> 573       if (max<-kCarTolerance*0.5)
                                                   >> 574     {
                                                   >> 575         tmax=max/vy;
                                                   >> 576         tmin=(fDy-yt)/vy;
                                                   >> 577     }
                                                   >> 578       else
                                                   >> 579     {
                                                   >> 580         return snxt=kInfinity;
                                                   >> 581     }
                                                   >> 582   }
                                                   >> 583     else
                                                   >> 584   {
                                                   >> 585       if (fabs(yt)<=fDy)
                                                   >> 586     {
                                                   >> 587         tmin=0;
                                                   >> 588         tmax=kInfinity;
                                                   >> 589     }
                                                   >> 590       else
                                                   >> 591     {
                                                   >> 592         return snxt=kInfinity;
                                                   >> 593     }
                                                   >> 594   }
                                                   >> 595 
                                                   >> 596 // Re-Calc valid intersection range
                                                   >> 597     if (tmin>smin) smin=tmin;
                                                   >> 598     if (tmax<smax) smax=tmax;
                                                   >> 599     if (smax<=smin)
                                                   >> 600   {
                                                   >> 601       return snxt=kInfinity;
                                                   >> 602   }
                                                   >> 603     else
                                                   >> 604   {
                                                   >> 605 //
                                                   >> 606 // X G4Parallel planes intersection
                                                   >> 607 //
                                                   >> 608       xt=p.x()-fTthetaCphi*p.z()-fTalpha*yt;
                                                   >> 609       vx=v.x()-fTthetaCphi*v.z()-fTalpha*vy;
                                                   >> 610       if (vx>0)
                                                   >> 611     {
                                                   >> 612         max=fDx-xt;
                                                   >> 613         if (max>kCarTolerance*0.5)
                                                   >> 614       {
                                                   >> 615           tmax=max/vx;
                                                   >> 616           tmin=(-fDx-xt)/vx;
                                                   >> 617       }
                                                   >> 618         else
                                                   >> 619       {
                                                   >> 620           return snxt=kInfinity;
                                                   >> 621       }
                                                   >> 622     }
                                                   >> 623       else if (vx<0)
                                                   >> 624     {
                                                   >> 625         max=-fDx-xt;
                                                   >> 626         if (max<-kCarTolerance*0.5)
                                                   >> 627       {
                                                   >> 628           tmax=max/vx;
                                                   >> 629           tmin=(fDx-xt)/vx;
                                                   >> 630       }
                                                   >> 631         else
                                                   >> 632       {
                                                   >> 633           return snxt=kInfinity;
                                                   >> 634       }
                                                   >> 635     }
                                                   >> 636       else
                                                   >> 637     {
                                                   >> 638         if (fabs(xt)<=fDx)
                                                   >> 639       {
                                                   >> 640           tmin=0;
                                                   >> 641           tmax=kInfinity;
                                                   >> 642       }
                                                   >> 643         else
                                                   >> 644       {
                                                   >> 645           return snxt=kInfinity;
                                                   >> 646       }
                                                   >> 647     }
                                                   >> 648       if (tmin>smin) smin=tmin;
                                                   >> 649       if (tmax<smax) smax=tmax;
                                                   >> 650   }
                                                   >> 651    
                                                   >> 652     if (smax>0&&smin<smax)
                                                   >> 653   {
                                                   >> 654       if (smin>0)
                                                   >> 655     {
                                                   >> 656         snxt=smin;
                                                   >> 657     }
                                                   >> 658       else
                                                   >> 659     {
                                                   >> 660         snxt=0;
                                                   >> 661     }
                                                   >> 662   }
                                                   >> 663     else
                                                   >> 664   {
                                                   >> 665       snxt=kInfinity;
                                                   >> 666   }
                                                   >> 667     return snxt;
620 }                                                 668 }
621                                                   669 
622 ////////////////////////////////////////////// << 670 ////////////////////////////////////////////////////////////////////////////
623 //                                                671 //
624 // Calculate exact shortest distance to any bo    672 // Calculate exact shortest distance to any boundary from outside
625 // - returns 0 is point inside                 << 673 // - Returns 0 is point inside
626                                                   674 
627 G4double G4Para::DistanceToIn( const G4ThreeVe << 675 G4double G4Para::DistanceToIn(const G4ThreeVector& p) const
628 {                                                 676 {
629   G4double xx = fPlanes[2].a*p.x()+fPlanes[2]. << 677     G4double safe;
630   G4double dx = std::abs(xx) + fPlanes[2].d;   << 678     G4double distz1,distz2,disty1,disty2,distx1,distx2;
631                                                << 679     G4double trany,cosy,tranx,cosx;
632   G4double yy = fPlanes[0].b*p.y()+fPlanes[0]. << 680 
633   G4double dy = std::abs(yy) + fPlanes[0].d;   << 681 // Z planes
634   G4double dxy = std::max(dx,dy);              << 682     distz1=p.z()-fDz;
635                                                << 683     distz2=-fDz-p.z();
636   G4double dz = std::abs(p.z())-fDz;           << 684     if (distz1>distz2)
637   G4double dist = std::max(dxy,dz);            << 685   {
                                                   >> 686       safe=distz1;
                                                   >> 687   }
                                                   >> 688     else
                                                   >> 689   {
                                                   >> 690       safe=distz2;
                                                   >> 691   }
                                                   >> 692 
                                                   >> 693     trany=p.y()-fTthetaSphi*p.z(); // Transformed y into `box' system
                                                   >> 694         // Transformed x into `box' system
                                                   >> 695     cosy=1.0/sqrt(1.0+fTthetaSphi*fTthetaSphi);
                                                   >> 696     disty1=(trany-fDy)*cosy;
                                                   >> 697     disty2=(-fDy-trany)*cosy;
                                                   >> 698     
                                                   >> 699     if (disty1>safe) safe=disty1;
                                                   >> 700     if (disty2>safe) safe=disty2;
638                                                   701 
639   return (dist > 0) ? dist : 0.;               << 702     tranx=p.x()-fTthetaCphi*p.z()-fTalpha*trany;
                                                   >> 703     cosx=1.0/sqrt(1.0+fTalpha*fTalpha+fTthetaCphi*fTthetaCphi);
                                                   >> 704     distx1=(tranx-fDx)*cosx;
                                                   >> 705     distx2=(-fDx-tranx)*cosx;
                                                   >> 706     
                                                   >> 707     if (distx1>safe) safe=distx1;
                                                   >> 708     if (distx2>safe) safe=distx2;
                                                   >> 709     
                                                   >> 710     if (safe<0) safe=0;
                                                   >> 711     return safe;  
640 }                                                 712 }
641                                                   713 
642 //////////////////////////////////////////////    714 //////////////////////////////////////////////////////////////////////////
643 //                                                715 //
644 // Calculate distance to surface of shape from << 716 // Calcluate distance to surface of shape from inside
645 // find normal at exit point                   << 717 // Calculate distance to x/y/z planes - smallest is exiting distance
646 // - when leaving the surface, return 0        << 
647                                                << 
648 G4double G4Para::DistanceToOut(const G4ThreeVe << 
649                                const G4bool ca << 
650                                      G4bool* v << 
651 {                                              << 
652   // Z intersections                           << 
653   //                                           << 
654   if ((std::abs(p.z()) - fDz) >= -halfCarToler << 
655   {                                            << 
656     if (calcNorm)                              << 
657     {                                          << 
658       *validNorm = true;                       << 
659       n->set(0, 0, (p.z() < 0) ? -1 : 1);      << 
660     }                                          << 
661     return 0.;                                 << 
662   }                                            << 
663   G4double vz = v.z();                         << 
664   G4double tmax = (vz == 0) ? DBL_MAX : (std:: << 
665   G4int iside = (vz < 0) ? -4 : -2; // little  << 
666                                                << 
667   // Y intersections                           << 
668   //                                           << 
669   G4double cos0 = fPlanes[0].b*v.y() + fPlanes << 
670   if (cos0 > 0)                                << 
671   {                                            << 
672     G4double dis0 = fPlanes[0].b*p.y() + fPlan << 
673     if (dis0 >= -halfCarTolerance)             << 
674     {                                          << 
675       if (calcNorm)                            << 
676       {                                        << 
677         *validNorm = true;                     << 
678         n->set(0, fPlanes[0].b, fPlanes[0].c); << 
679       }                                        << 
680       return 0.;                               << 
681     }                                          << 
682     G4double tmp = -dis0/cos0;                 << 
683     if (tmax > tmp) { tmax = tmp; iside = 0; } << 
684   }                                            << 
685                                                   718 
686   G4double cos1 = -cos0;                       << 719 G4double G4Para::DistanceToOut(const G4ThreeVector& p,const G4ThreeVector& v,
687   if (cos1 > 0)                                << 720              const G4bool calcNorm,
688   {                                            << 721              G4bool *validNorm,G4ThreeVector *n) const
689     G4double dis1 = fPlanes[1].b*p.y() + fPlan << 722 {
690     if (dis1 >= -halfCarTolerance)             << 723     ESide side = kUndef;
691     {                                          << 724     G4double snxt;    // snxt = return value
692       if (calcNorm)                            << 725     G4double max,tmax;
693       {                                        << 726     G4double yt,vy,xt,vx;
694         *validNorm = true;                     << 727 
695         n->set(0, fPlanes[1].b, fPlanes[1].c); << 728     G4double ycomp,calpha,salpha,tntheta,cosntheta;
696       }                                        << 729 //
697       return 0.;                               << 730 // Z Intersections
698     }                                          << 731 //
699     G4double tmp = -dis1/cos1;                 << 732     if (v.z()>0)
700     if (tmax > tmp) { tmax = tmp; iside = 1; } << 733   {
701   }                                            << 734       max=fDz-p.z();
                                                   >> 735       if (max>kCarTolerance*0.5)
                                                   >> 736     {
                                                   >> 737         snxt=max/v.z();
                                                   >> 738         side=kPZ;
                                                   >> 739     }
                                                   >> 740       else
                                                   >> 741     {
                                                   >> 742         if (calcNorm)
                                                   >> 743       {
                                                   >> 744           *validNorm=true;
                                                   >> 745           *n=G4ThreeVector(0,0,1);
                                                   >> 746       }
                                                   >> 747         return snxt=0;
                                                   >> 748     }
                                                   >> 749   }
                                                   >> 750     else if (v.z()<0)
                                                   >> 751   {
                                                   >> 752       max=-fDz-p.z();
                                                   >> 753       if (max<-kCarTolerance*0.5)
                                                   >> 754     {
                                                   >> 755         snxt=max/v.z();
                                                   >> 756         side=kMZ;
                                                   >> 757     }
                                                   >> 758       else
                                                   >> 759     {
                                                   >> 760         if (calcNorm)
                                                   >> 761       {
                                                   >> 762           *validNorm=true;
                                                   >> 763           *n=G4ThreeVector(0,0,-1);
                                                   >> 764       }
                                                   >> 765         return snxt=0;
                                                   >> 766     }
                                                   >> 767   }
                                                   >> 768     else
                                                   >> 769   {
                                                   >> 770       snxt=kInfinity;
                                                   >> 771   }
702                                                   772 
703   // X intersections                           << 773     
704   //                                           << 774 //
705   G4double cos2 = fPlanes[2].a*v.x() + fPlanes << 775 // Y plane intersection
706   if (cos2 > 0)                                << 776 //
707   {                                            << 777     yt=p.y()-fTthetaSphi*p.z();
708     G4double dis2 = fPlanes[2].a*p.x()+fPlanes << 778     vy=v.y()-fTthetaSphi*v.z();
709     if (dis2 >= -halfCarTolerance)             << 
710     {                                          << 
711       if (calcNorm)                            << 
712       {                                        << 
713          *validNorm = true;                    << 
714          n->set(fPlanes[2].a, fPlanes[2].b, fP << 
715       }                                        << 
716       return 0.;                               << 
717     }                                          << 
718     G4double tmp = -dis2/cos2;                 << 
719     if (tmax > tmp) { tmax = tmp; iside = 2; } << 
720   }                                            << 
721                                                   779 
722   G4double cos3 = -cos2;                       << 780     if (vy>0)
723   if (cos3 > 0)                                << 781   {
724   {                                            << 782       max=fDy-yt;
725     G4double dis3 = fPlanes[3].a*p.x()+fPlanes << 783       if (max>kCarTolerance*0.5)
726     if (dis3 >= -halfCarTolerance)             << 784     {
727     {                                          << 785         tmax=max/vy;
728       if (calcNorm)                            << 786         if (tmax<snxt)
729       {                                        << 787       {
730          *validNorm = true;                    << 788           snxt=tmax;
731          n->set(fPlanes[3].a, fPlanes[3].b, fP << 789           side=kPY;
732       }                                        << 790       }
733       return 0.;                               << 791     }
734     }                                          << 792       else
735     G4double tmp = -dis3/cos3;                 << 793     {
736     if (tmax > tmp) { tmax = tmp; iside = 3; } << 794         if (calcNorm)
737   }                                            << 795       {      
                                                   >> 796           *validNorm=true; // Leaving via plus Y
                                                   >> 797           ycomp=1/sqrt(1+fTthetaSphi*fTthetaSphi);
                                                   >> 798           *n=G4ThreeVector(0,ycomp,-fTthetaSphi*ycomp);
                                                   >> 799       }
                                                   >> 800         return snxt=0;
                                                   >> 801     }
                                                   >> 802   }
                                                   >> 803     else if (vy<0)
                                                   >> 804   {
                                                   >> 805       max=-fDy-yt;
                                                   >> 806       if (max<-kCarTolerance*0.5)
                                                   >> 807     {
                                                   >> 808         tmax=max/vy;
                                                   >> 809         if (tmax<snxt)
                                                   >> 810       {
                                                   >> 811           snxt=tmax;
                                                   >> 812           side=kMY;
                                                   >> 813       }
                                                   >> 814     }
                                                   >> 815       else
                                                   >> 816     {
                                                   >> 817         if (calcNorm)
                                                   >> 818       {
                                                   >> 819           *validNorm=true; // Leaving via minus Y
                                                   >> 820           ycomp=-1/sqrt(1+fTthetaSphi*fTthetaSphi);
                                                   >> 821           *n=G4ThreeVector(0,ycomp,-fTthetaSphi*ycomp);
                                                   >> 822 
                                                   >> 823       }
                                                   >> 824         return snxt=0;
                                                   >> 825     }
                                                   >> 826   }
                                                   >> 827 //
                                                   >> 828 // X plane intersection
                                                   >> 829 //
                                                   >> 830     xt=p.x()-fTthetaCphi*p.z()-fTalpha*yt;
                                                   >> 831     vx=v.x()-fTthetaCphi*v.z()-fTalpha*vy;
                                                   >> 832     if (vx>0)
                                                   >> 833   {
                                                   >> 834       max=fDx-xt;
                                                   >> 835       if (max>kCarTolerance*0.5)
                                                   >> 836     {
                                                   >> 837         tmax=max/vx;
                                                   >> 838         if (tmax<snxt)
                                                   >> 839       {
                                                   >> 840           snxt=tmax;
                                                   >> 841           side=kPX;
                                                   >> 842       }
                                                   >> 843     }
                                                   >> 844       else
                                                   >> 845     {
                                                   >> 846         if (calcNorm)
                                                   >> 847       {
                                                   >> 848           *validNorm=true; // Leaving via plus X
                                                   >> 849           calpha=1/sqrt(1+fTalpha*fTalpha);
                                                   >> 850           if (fTalpha)
                                                   >> 851         {
                                                   >> 852             salpha=-calpha/fTalpha; // NOTE: actually use MINUS sin(alpha)
                                                   >> 853         }
                                                   >> 854           else
                                                   >> 855         {
                                                   >> 856             salpha=0;
                                                   >> 857         }
                                                   >> 858           tntheta=fTthetaCphi*calpha+fTthetaSphi*salpha;
                                                   >> 859           cosntheta=1/sqrt(1+tntheta*tntheta);
                                                   >> 860           *n=G4ThreeVector(calpha*cosntheta,salpha*cosntheta,-tntheta*cosntheta);
                                                   >> 861       }
                                                   >> 862         return snxt=0;
                                                   >> 863     }
                                                   >> 864   }
                                                   >> 865     else if (vx<0)
                                                   >> 866   {
                                                   >> 867       max=-fDx-xt;
                                                   >> 868       if (max<-kCarTolerance*0.5)
                                                   >> 869     {
                                                   >> 870         tmax=max/vx;
                                                   >> 871         if (tmax<snxt)
                                                   >> 872       {
                                                   >> 873           snxt=tmax;
                                                   >> 874           side=kMX;
                                                   >> 875       }
                                                   >> 876     }
                                                   >> 877       else
                                                   >> 878     {
                                                   >> 879         if (calcNorm)
                                                   >> 880       {
                                                   >> 881           *validNorm=true; // Leaving via minus X
                                                   >> 882           calpha=1/sqrt(1+fTalpha*fTalpha);
                                                   >> 883           if (fTalpha)
                                                   >> 884         {
                                                   >> 885             salpha=-calpha/fTalpha; // NOTE: actually use MINUS sin(alpha)
                                                   >> 886         }
                                                   >> 887           else
                                                   >> 888         {
                                                   >> 889             salpha=0;
                                                   >> 890         }
                                                   >> 891           tntheta=fTthetaCphi*calpha+fTthetaSphi*salpha;
                                                   >> 892           cosntheta=-1/sqrt(1+tntheta*tntheta);
                                                   >> 893           *n=G4ThreeVector(calpha*cosntheta,salpha*cosntheta,-tntheta*cosntheta);          
                                                   >> 894           return snxt=0;
                                                   >> 895       }
                                                   >> 896     }
                                                   >> 897   }
738                                                   898 
739   // Set normal, if required, and return dista << 899     if (calcNorm)
740   //                                           << 900   {
741   if (calcNorm)                                << 901       *validNorm=true;
742   {                                            << 902       switch (side)
743     *validNorm = true;                         << 903     {
744     if (iside < 0)                             << 904     case kMZ:
745       n->set(0, 0, iside + 3); // (-4+3)=-1, ( << 905         *n=G4ThreeVector(0,0,-1);
746     else                                       << 906         break;
747       n->set(fPlanes[iside].a, fPlanes[iside]. << 907     case kPZ:
748   }                                            << 908         *n=G4ThreeVector(0,0,1);
749   return tmax;                                 << 909         break;
                                                   >> 910     case kMY:
                                                   >> 911         ycomp=-1/sqrt(1+fTthetaSphi*fTthetaSphi);
                                                   >> 912         *n=G4ThreeVector(0,ycomp,-fTthetaSphi*ycomp);
                                                   >> 913         break;        
                                                   >> 914     case kPY:
                                                   >> 915         ycomp=1/sqrt(1+fTthetaSphi*fTthetaSphi);
                                                   >> 916         *n=G4ThreeVector(0,ycomp,-fTthetaSphi*ycomp);
                                                   >> 917         break;        
                                                   >> 918     case kMX:
                                                   >> 919         calpha=1/sqrt(1+fTalpha*fTalpha);
                                                   >> 920         if (fTalpha)
                                                   >> 921       {
                                                   >> 922           salpha=-calpha/fTalpha; // NOTE: actually use MINUS sin(alpha)
                                                   >> 923       }
                                                   >> 924         else
                                                   >> 925       {
                                                   >> 926           salpha=0;
                                                   >> 927       }
                                                   >> 928         tntheta=fTthetaCphi*calpha+fTthetaSphi*salpha;
                                                   >> 929         cosntheta=-1/sqrt(1+tntheta*tntheta);
                                                   >> 930         *n=G4ThreeVector(calpha*cosntheta,salpha*cosntheta,-tntheta*cosntheta);
                                                   >> 931         break;
                                                   >> 932     case kPX:
                                                   >> 933         calpha=1/sqrt(1+fTalpha*fTalpha);
                                                   >> 934         if (fTalpha)
                                                   >> 935       {
                                                   >> 936           salpha=-calpha/fTalpha; // NOTE: actually use MINUS sin(alpha)
                                                   >> 937       }
                                                   >> 938         else
                                                   >> 939       {
                                                   >> 940           salpha=0;
                                                   >> 941       }
                                                   >> 942         tntheta=fTthetaCphi*calpha+fTthetaSphi*salpha;
                                                   >> 943         cosntheta=1/sqrt(1+tntheta*tntheta);
                                                   >> 944         *n=G4ThreeVector(calpha*cosntheta,salpha*cosntheta,-tntheta*cosntheta);
                                                   >> 945         break;
                                                   >> 946     default:
                                                   >> 947         G4Exception("Invalid enum in G4Para::DistanceToOut");
                                                   >> 948         break;
                                                   >> 949         
                                                   >> 950     }
                                                   >> 951   }
                                                   >> 952     return snxt;
750 }                                                 953 }
751                                                   954 
752 ////////////////////////////////////////////// << 955 /////////////////////////////////////////////////////////////////////////////
753 //                                                956 //
754 // Calculate exact shortest distance to any bo    957 // Calculate exact shortest distance to any boundary from inside
755 // - returns 0 is point outside                << 958 // - Returns 0 is point outside
756                                                   959 
757 G4double G4Para::DistanceToOut( const G4ThreeV << 960 G4double G4Para::DistanceToOut(const G4ThreeVector& p) const
758 {                                                 961 {
                                                   >> 962     G4double safe;
                                                   >> 963     G4double distz1,distz2,disty1,disty2,distx1,distx2;
                                                   >> 964     G4double trany,cosy,tranx,cosx;
                                                   >> 965 
759 #ifdef G4CSGDEBUG                                 966 #ifdef G4CSGDEBUG
760   if( Inside(p) == kOutside )                     967   if( Inside(p) == kOutside )
761   {                                               968   {
762     std::ostringstream message;                << 969      G4cout.precision(16) ;
763     G4int oldprc = message.precision(16);      << 970      G4cout << G4endl ;
764     message << "Point p is outside (!?) of sol << 971      G4cout << "Parallelpiped parameters:" << G4endl << G4endl ;
765     message << "Position:\n";                  << 972      G4cout << "fDx = "     << fDx/mm << " mm" << G4endl ;
766     message << "   p.x() = " << p.x()/mm << "  << 973      G4cout << "fDy = "     << fDy/mm << " mm" << G4endl ;
767     message << "   p.y() = " << p.y()/mm << "  << 974      G4cout << "fDz = "     << fDz/mm << " mm" << G4endl ;
768     message << "   p.z() = " << p.z()/mm << "  << 975      G4cout << "fTalpha =     " << fTalpha/degree << " degree" << G4endl;
769     G4cout.precision(oldprc) ;                 << 976      G4cout << "fTthetaCphi = " << fTthetaCphi/degree << " degree" << G4endl;
770     G4Exception("G4Para::DistanceToOut(p)", "G << 977      G4cout << "fTthetaSphi = " << fTthetaSphi/degree << " degree" << G4endl << G4endl ;
771                 JustWarning, message );        << 978      G4cout << "Position:"  << G4endl << G4endl ;
772     DumpInfo();                                << 979      G4cout << "p.x() = "   << p.x()/mm << " mm" << G4endl ;
773     }                                          << 980      G4cout << "p.y() = "   << p.y()/mm << " mm" << G4endl ;
                                                   >> 981      G4cout << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl ;
                                                   >> 982      G4cout << "G4Para::DistanceToOut(p) - point p is outside ?!" << G4endl ;
                                                   >> 983      // G4Exception("Invalid call in G4Para::DistanceToOut(p), point p is outside") ;
                                                   >> 984   }
774 #endif                                            985 #endif
775   G4double xx = fPlanes[2].a*p.x()+fPlanes[2]. << 
776   G4double dx = std::abs(xx) + fPlanes[2].d;   << 
777                                                << 
778   G4double yy = fPlanes[0].b*p.y()+fPlanes[0]. << 
779   G4double dy = std::abs(yy) + fPlanes[0].d;   << 
780   G4double dxy = std::max(dx,dy);              << 
781                                                   986 
782   G4double dz = std::abs(p.z())-fDz;           << 987 // Z planes
783   G4double dist = std::max(dxy,dz);            << 988     distz1=fDz-p.z();
                                                   >> 989     distz2=fDz+p.z();
                                                   >> 990     if (distz1<distz2)
                                                   >> 991   {
                                                   >> 992       safe=distz1;
                                                   >> 993   }
                                                   >> 994     else
                                                   >> 995   {
                                                   >> 996       safe=distz2;
                                                   >> 997   }
                                                   >> 998 
                                                   >> 999     trany=p.y()-fTthetaSphi*p.z(); // Transformed y into `box' system
                                                   >> 1000         // Transformed x into `box' system
                                                   >> 1001     cosy=1.0/sqrt(1.0+fTthetaSphi*fTthetaSphi);
                                                   >> 1002     disty1=(fDy-trany)*cosy;
                                                   >> 1003     disty2=(fDy+trany)*cosy;
                                                   >> 1004     
                                                   >> 1005     if (disty1<safe) safe=disty1;
                                                   >> 1006     if (disty2<safe) safe=disty2;
784                                                   1007 
785   return (dist < 0) ? -dist : 0.;              << 1008     tranx=p.x()-fTthetaCphi*p.z()-fTalpha*trany;
                                                   >> 1009     cosx=1.0/sqrt(1.0+fTalpha*fTalpha+fTthetaCphi*fTthetaCphi);
                                                   >> 1010     distx1=(fDx-tranx)*cosx;
                                                   >> 1011     distx2=(fDx+tranx)*cosx;
                                                   >> 1012     
                                                   >> 1013     if (distx1<safe) safe=distx1;
                                                   >> 1014     if (distx2<safe) safe=distx2;
                                                   >> 1015     
                                                   >> 1016     if (safe<0) safe=0;
                                                   >> 1017     return safe;  
786 }                                                 1018 }
787                                                   1019 
788 ////////////////////////////////////////////// << 1020 ////////////////////////////////////////////////////////////////////////////////
789 //                                                1021 //
790 // GetEntityType                               << 1022 // Create a List containing the transformed vertices
791                                                << 1023 // Ordering [0-3] -fDz cross section
792 G4GeometryType G4Para::GetEntityType() const   << 1024 //          [4-7] +fDz cross section such that [0] is below [4],
793 {                                              << 1025 //                                             [1] below [5] etc.
794   return {"G4Para"};                           << 1026 // Note:
                                                   >> 1027 //  Caller has deletion resposibility
                                                   >> 1028 
                                                   >> 1029 G4ThreeVectorList*
                                                   >> 1030    G4Para::CreateRotatedVertices(const G4AffineTransform& pTransform) const
                                                   >> 1031 {
                                                   >> 1032     G4ThreeVectorList *vertices;
                                                   >> 1033     vertices=new G4ThreeVectorList();
                                                   >> 1034     vertices->reserve(8);
                                                   >> 1035     if (vertices)
                                                   >> 1036   {
                                                   >> 1037       G4ThreeVector vertex0(-fDz*fTthetaCphi-fDy*fTalpha-fDx,-fDz*fTthetaSphi-fDy,-fDz);
                                                   >> 1038       G4ThreeVector vertex1(-fDz*fTthetaCphi-fDy*fTalpha+fDx,-fDz*fTthetaSphi-fDy,-fDz);
                                                   >> 1039       G4ThreeVector vertex2(-fDz*fTthetaCphi+fDy*fTalpha-fDx,-fDz*fTthetaSphi+fDy,-fDz);
                                                   >> 1040       G4ThreeVector vertex3(-fDz*fTthetaCphi+fDy*fTalpha+fDx,-fDz*fTthetaSphi+fDy,-fDz);
                                                   >> 1041       G4ThreeVector vertex4(+fDz*fTthetaCphi-fDy*fTalpha-fDx,+fDz*fTthetaSphi-fDy,+fDz);
                                                   >> 1042       G4ThreeVector vertex5(+fDz*fTthetaCphi-fDy*fTalpha+fDx,+fDz*fTthetaSphi-fDy,+fDz);
                                                   >> 1043       G4ThreeVector vertex6(+fDz*fTthetaCphi+fDy*fTalpha-fDx,+fDz*fTthetaSphi+fDy,+fDz);
                                                   >> 1044       G4ThreeVector vertex7(+fDz*fTthetaCphi+fDy*fTalpha+fDx,+fDz*fTthetaSphi+fDy,+fDz);
                                                   >> 1045 
                                                   >> 1046       vertices->push_back(pTransform.TransformPoint(vertex0));
                                                   >> 1047       vertices->push_back(pTransform.TransformPoint(vertex1));
                                                   >> 1048       vertices->push_back(pTransform.TransformPoint(vertex2));
                                                   >> 1049       vertices->push_back(pTransform.TransformPoint(vertex3));
                                                   >> 1050       vertices->push_back(pTransform.TransformPoint(vertex4));
                                                   >> 1051       vertices->push_back(pTransform.TransformPoint(vertex5));
                                                   >> 1052       vertices->push_back(pTransform.TransformPoint(vertex6));
                                                   >> 1053       vertices->push_back(pTransform.TransformPoint(vertex7));
                                                   >> 1054   }
                                                   >> 1055     else
                                                   >> 1056   {
                                                   >> 1057       G4Exception("G4Para::CreateRotatedVertices - Out of memory !");
                                                   >> 1058   }
                                                   >> 1059     return vertices;
795 }                                                 1060 }
796                                                   1061 
797 ////////////////////////////////////////////// << 1062 ////////////////////////////////////////////////////////////////////////////
798 //                                                1063 //
799 // IsFaceted                                   << 1064 // Methods for visualisation
800                                                   1065 
801 G4bool G4Para::IsFaceted() const               << 1066 void G4Para::DescribeYourselfTo (G4VGraphicsScene& scene) const
802 {                                                 1067 {
803   return true;                                 << 1068    scene.AddThis (*this);
804 }                                                 1069 }
805                                                   1070 
806 ////////////////////////////////////////////// << 1071 G4Polyhedron* G4Para::CreatePolyhedron () const
807 //                                             << 
808 // Make a clone of the object                  << 
809 //                                             << 
810 G4VSolid* G4Para::Clone() const                << 
811 {                                                 1072 {
812   return new G4Para(*this);                    << 1073     G4double phi = atan2(fTthetaSphi, fTthetaCphi);
                                                   >> 1074     G4double alpha = atan(fTalpha);
                                                   >> 1075     G4double theta = atan(sqrt(fTthetaCphi*fTthetaCphi
                                                   >> 1076                                +fTthetaSphi*fTthetaSphi));
                                                   >> 1077     
                                                   >> 1078     return new G4PolyhedronPara(fDx, fDy, fDz, alpha, theta, phi);
813 }                                                 1079 }
814                                                   1080 
815 ////////////////////////////////////////////// << 1081 G4NURBS* G4Para::CreateNURBS () const
816 //                                             << 
817 // Stream object contents to an output stream  << 
818                                                << 
819 std::ostream& G4Para::StreamInfo( std::ostream << 
820 {                                                 1082 {
821   G4double alpha = std::atan(fTalpha);         << 1083    // return new G4NURBSbox (fDx, fDy, fDz);
822   G4double theta = std::atan(std::sqrt(fTtheta << 1084    return 0 ;
823                                        fTtheta << 
824   G4double phi   = std::atan2(fTthetaSphi,fTth << 
825                                                << 
826   G4long oldprc = os.precision(16);            << 
827   os << "------------------------------------- << 
828      << "    *** Dump for solid - " << GetName << 
829      << "    ================================= << 
830      << " Solid type: G4Para\n"                << 
831      << " Parameters:\n"                       << 
832      << "    half length X: " << fDx/mm << " m << 
833      << "    half length Y: " << fDy/mm << " m << 
834      << "    half length Z: " << fDz/mm << " m << 
835      << "    alpha: " << alpha/degree << "degr << 
836      << "    theta: " << theta/degree << "degr << 
837      << "    phi: " << phi/degree << "degrees\ << 
838      << "------------------------------------- << 
839   os.precision(oldprc);                        << 
840                                                << 
841   return os;                                   << 
842 }                                                 1085 }
843                                                   1086 
844 ////////////////////////////////////////////// << 
845 //                                                1087 //
846 // Return a point randomly and uniformly selec << 
847                                                << 
848 G4ThreeVector G4Para::GetPointOnSurface() cons << 
849 {                                              << 
850   G4double DyTalpha = fDy*fTalpha;             << 
851   G4double DzTthetaSphi = fDz*fTthetaSphi;     << 
852   G4double DzTthetaCphi = fDz*fTthetaCphi;     << 
853                                                << 
854   // Set vertices                              << 
855   //                                           << 
856   G4ThreeVector pt[8];                         << 
857   pt[0].set(-DzTthetaCphi-DyTalpha-fDx, -DzTth << 
858   pt[1].set(-DzTthetaCphi-DyTalpha+fDx, -DzTth << 
859   pt[2].set(-DzTthetaCphi+DyTalpha-fDx, -DzTth << 
860   pt[3].set(-DzTthetaCphi+DyTalpha+fDx, -DzTth << 
861   pt[4].set( DzTthetaCphi-DyTalpha-fDx,  DzTth << 
862   pt[5].set( DzTthetaCphi-DyTalpha+fDx,  DzTth << 
863   pt[6].set( DzTthetaCphi+DyTalpha-fDx,  DzTth << 
864   pt[7].set( DzTthetaCphi+DyTalpha+fDx,  DzTth << 
865                                                << 
866   // Set areas (-Z, -Y, +Y, -X, +X, +Z)        << 
867   //                                           << 
868   G4ThreeVector vx(fDx, 0, 0);                 << 
869   G4ThreeVector vy(DyTalpha, fDy, 0);          << 
870   G4ThreeVector vz(DzTthetaCphi, DzTthetaSphi, << 
871                                                << 
872   G4double sxy = fDx*fDy; // (vx.cross(vy)).ma << 
873   G4double sxz = (vx.cross(vz)).mag();         << 
874   G4double syz = (vy.cross(vz)).mag();         << 
875                                                << 
876   G4double sface[6] = { sxy, syz, syz, sxz, sx << 
877   for (G4int i=1; i<6; ++i) { sface[i] += sfac << 
878                                                << 
879   // Select face                               << 
880   //                                           << 
881   G4double select = sface[5]*G4UniformRand();  << 
882   G4int k = 5;                                 << 
883   if (select <= sface[4]) k = 4;               << 
884   if (select <= sface[3]) k = 3;               << 
885   if (select <= sface[2]) k = 2;               << 
886   if (select <= sface[1]) k = 1;               << 
887   if (select <= sface[0]) k = 0;               << 
888                                                << 
889   // Generate point                            << 
890   //                                           << 
891   G4int ip[6][3] = {{0,1,2}, {0,4,1}, {2,3,6}, << 
892   G4double u = G4UniformRand();                << 
893   G4double v = G4UniformRand();                << 
894   return (1.-u-v)*pt[ip[k][0]] + u*pt[ip[k][1] << 
895 }                                              << 
896                                                << 
897 ////////////////////////////////////////////// << 
898 //                                                1088 //
899 // Methods for visualisation                   << 1089 ///////////////////////////  End of G4Para.cc ///////////////////////////
900                                                << 
901 void G4Para::DescribeYourselfTo ( G4VGraphicsS << 
902 {                                              << 
903   scene.AddSolid (*this);                      << 
904 }                                              << 
905                                                << 
906 G4Polyhedron* G4Para::CreatePolyhedron () cons << 
907 {                                              << 
908   G4double phi = std::atan2(fTthetaSphi, fTthe << 
909   G4double alpha = std::atan(fTalpha);         << 
910   G4double theta = std::atan(std::sqrt(fTtheta << 
911                                        fTtheta << 
912                                                << 
913   return new G4PolyhedronPara(fDx, fDy, fDz, a << 
914 }                                              << 
915 #endif                                         << 
916                                                   1090