Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/solids/CSG/src/G4Trd.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/G4Trd.cc (Version 11.3.0) and /geometry/solids/CSG/src/G4Trd.cc (Version 5.2.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: G4Trd.cc,v 1.13 2003/06/16 16:53:43 gunter Exp $
                                                   >>  25 // GEANT4 tag $Name: geant4-05-02-patch-01 $
                                                   >>  26 //
                                                   >>  27 //
 26 // Implementation for G4Trd class                  28 // Implementation for G4Trd class
 27 //                                                 29 //
 28 // 12.01.95 P.Kent: First version              <<  30 // History:
 29 // 28.04.05 V.Grichine: new SurfaceNormal acco <<  31 //    ~1996, V.Grichine, 1st implementation based on old code of P.Kent
 30 // 25.05.17 E.Tcherniaev: complete revision, s <<  32 // 07.05.00, V.Grichine, in d = DistanceToIn(p,v), if d<0.5*kCarTolerance, d=0
 31 // ------------------------------------------- <<  33 //
                                                   >>  34 // ********************************************************************
 32                                                    35 
 33 #include "G4Trd.hh"                                36 #include "G4Trd.hh"
 34                                                    37 
 35 #if !defined(G4GEOM_USE_UTRD)                  <<  38 #include "G4VPVParameterisation.hh"
 36                                                << 
 37 #include "G4GeomTools.hh"                      << 
 38                                                << 
 39 #include "G4VoxelLimits.hh"                        39 #include "G4VoxelLimits.hh"
 40 #include "G4AffineTransform.hh"                    40 #include "G4AffineTransform.hh"
 41 #include "G4BoundingEnvelope.hh"               << 
 42 #include "G4QuickRand.hh"                      << 
 43                                                << 
 44 #include "G4VPVParameterisation.hh"            << 
 45                                                    41 
 46 #include "G4VGraphicsScene.hh"                     42 #include "G4VGraphicsScene.hh"
                                                   >>  43 #include "G4Polyhedron.hh"
                                                   >>  44 #include "G4NURBS.hh"
                                                   >>  45 #include "G4NURBSbox.hh"
 47                                                    46 
 48 using namespace CLHEP;                         <<  47 /////////////////////////////////////////////////////////////////////////
 49                                                << 
 50 ////////////////////////////////////////////// << 
 51 //                                             << 
 52 // Constructor - set & check half widths       << 
 53                                                << 
 54 G4Trd::G4Trd(const G4String& pName,            << 
 55                    G4double pdx1, G4double pdx << 
 56                    G4double pdy1, G4double pdy << 
 57                    G4double pdz)               << 
 58   : G4CSGSolid(pName), halfCarTolerance(0.5*kC << 
 59     fDx1(pdx1), fDx2(pdx2), fDy1(pdy1), fDy2(p << 
 60 {                                              << 
 61   CheckParameters();                           << 
 62   MakePlanes();                                << 
 63 }                                              << 
 64                                                << 
 65 ////////////////////////////////////////////// << 
 66 //                                             << 
 67 // Fake default constructor - sets only member << 
 68 //                            for usage restri << 
 69 //                                             << 
 70 G4Trd::G4Trd( __void__& a )                    << 
 71   : G4CSGSolid(a), halfCarTolerance(0.5*kCarTo << 
 72     fDx1(1.), fDx2(1.), fDy1(1.), fDy2(1.), fD << 
 73 {                                              << 
 74   MakePlanes();                                << 
 75 }                                              << 
 76                                                << 
 77 ////////////////////////////////////////////// << 
 78 //                                             << 
 79 // Destructor                                  << 
 80                                                << 
 81 G4Trd::~G4Trd() = default;                     << 
 82                                                << 
 83 ////////////////////////////////////////////// << 
 84 //                                             << 
 85 // Copy constructor                            << 
 86                                                << 
 87 G4Trd::G4Trd(const G4Trd& rhs)                 << 
 88   : G4CSGSolid(rhs), halfCarTolerance(rhs.half << 
 89     fDx1(rhs.fDx1), fDx2(rhs.fDx2),            << 
 90     fDy1(rhs.fDy1), fDy2(rhs.fDy2), fDz(rhs.fD << 
 91     fHx(rhs.fHx), fHy(rhs.fHy)                 << 
 92 {                                              << 
 93   for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs << 
 94 }                                              << 
 95                                                << 
 96 ////////////////////////////////////////////// << 
 97 //                                             << 
 98 // Assignment operator                         << 
 99                                                << 
100 G4Trd& G4Trd::operator = (const G4Trd& rhs)    << 
101 {                                              << 
102    // Check assignment to self                 << 
103    //                                          << 
104    if (this == &rhs)  { return *this; }        << 
105                                                << 
106    // Copy base class data                     << 
107    //                                          << 
108    G4CSGSolid::operator=(rhs);                 << 
109                                                << 
110    // Copy data                                << 
111    //                                          << 
112    halfCarTolerance = rhs.halfCarTolerance;    << 
113    fDx1 = rhs.fDx1; fDx2 = rhs.fDx2;           << 
114    fDy1 = rhs.fDy1; fDy2 = rhs.fDy2;           << 
115    fDz = rhs.fDz;                              << 
116    fHx = rhs.fHx; fHy = rhs.fHy;               << 
117    for (G4int i=0; i<4; ++i) { fPlanes[i] = rh << 
118                                                << 
119    return *this;                               << 
120 }                                              << 
121                                                << 
122 ////////////////////////////////////////////// << 
123 //                                                 48 //
124 // Set all parameters, as for constructor - se <<  49 // Constructor - check & set half widths
125                                                    50 
126 void G4Trd::SetAllParameters(G4double pdx1, G4 <<  51 G4Trd::G4Trd( const G4String& pName,
127                              G4double pdy1, G4 <<  52                     G4double pdx1,  G4double pdx2,
                                                   >>  53                     G4double pdy1,  G4double pdy2,
                                                   >>  54                     G4double pdz )
                                                   >>  55   : G4CSGSolid(pName)
128 {                                                  56 {
129   // Reset data of the base class              <<  57   CheckAndSetAllParameters (pdx1, pdx2, pdy1, pdy2, pdz);
130   fCubicVolume = 0.;                           << 
131   fSurfaceArea = 0.;                           << 
132   fRebuildPolyhedron = true;                   << 
133                                                << 
134   // Set parameters                            << 
135   fDx1 = pdx1; fDx2 = pdx2;                    << 
136   fDy1 = pdy1; fDy2 = pdy2;                    << 
137   fDz  = pdz;                                  << 
138                                                << 
139   CheckParameters();                           << 
140   MakePlanes();                                << 
141 }                                                  58 }
142                                                    59 
143 ////////////////////////////////////////////// <<  60 /////////////////////////////////////////////////////////////////////////
144 //                                                 61 //
145 // Check dimensions                            <<  62 // Set and check (coplanarity) of trd parameters
146                                                    63 
147 void G4Trd::CheckParameters()                  <<  64 void G4Trd::CheckAndSetAllParameters ( G4double pdx1,  G4double pdx2,
                                                   >>  65                                        G4double pdy1,  G4double pdy2,
                                                   >>  66                                        G4double pdz ) 
148 {                                                  67 {
149   G4double dmin = 2*kCarTolerance;             <<  68   if ( pdx1>0&&pdx2>0&&pdy1>0&&pdy2>0&&pdz>0 )
150   if ((fDx1 < 0 || fDx2 < 0 || fDy1 < 0 || fDy <<  69   {
151       (fDx1 < dmin && fDx2 < dmin) ||          <<  70     fDx1=pdx1; fDx2=pdx2;
152       (fDy1 < dmin && fDy2 < dmin))            <<  71     fDy1=pdy1; fDy2=pdy2;
153   {                                            <<  72     fDz=pdz;
154     std::ostringstream message;                <<  73   }
155     message << "Invalid (too small or negative <<  74   else
156             << GetName()                       <<  75   {
157             << "\n  X - " << fDx1 << ", " << f <<  76     if ( pdx1>=0 && pdx2>=0 && pdy1>=0 && pdy2>=0 && pdz>=0 )
158             << "\n  Y - " << fDy1 << ", " << f <<  77     {
159             << "\n  Z - " << fDz;              <<  78       // G4double  Minimum_length= (1+per_thousand) * kCarTolerance/2.;
160     G4Exception("G4Trd::CheckParameters()", "G <<  79       // FIX-ME : temporary solution for ZERO or very-small parameters
161                 FatalException, message);      <<  80       //
                                                   >>  81       G4double  Minimum_length= kCarTolerance/2.;
                                                   >>  82       fDx1=std::max(pdx1,Minimum_length); 
                                                   >>  83       fDx2=std::max(pdx2,Minimum_length); 
                                                   >>  84       fDy1=std::max(pdy1,Minimum_length); 
                                                   >>  85       fDy2=std::max(pdy2,Minimum_length); 
                                                   >>  86       fDz=std::max(pdz,Minimum_length);
                                                   >>  87     }
                                                   >>  88     else
                                                   >>  89     {
                                                   >>  90       G4cout << "ERROR - G4Trd()::CheckAndSetAllParameters(): " << GetName()
                                                   >>  91              << G4endl
                                                   >>  92              << "        Invalid dimensions, some are < 0 !" << G4endl
                                                   >>  93              << "          X - " << pdx1 << ", " << pdx2 << G4endl
                                                   >>  94              << "          Y - " << pdy1 << ", " << pdy2 << G4endl
                                                   >>  95              << "          Z - " << pdz << G4endl;
                                                   >>  96       G4cerr << "ERROR - G4Trd()::CheckAndSetAllParameters(): " << GetName()
                                                   >>  97              << G4endl
                                                   >>  98              << "        Invalid dimensions, some are < 0 !" << G4endl
                                                   >>  99              << "          X - " << pdx1 << ", " << pdx2 << G4endl
                                                   >> 100              << "          Y - " << pdy1 << ", " << pdy2 << G4endl
                                                   >> 101              << "          Z - " << pdz << G4endl;
                                                   >> 102       G4Exception("G4Trd::CheckAndSetAllParameters() - Invalid parameters");
                                                   >> 103     }
162   }                                               104   }
163 }                                                 105 }
164                                                   106 
165 //////////////////////////////////////////////    107 //////////////////////////////////////////////////////////////////////////
166 //                                                108 //
167 // Set side planes                             << 109 // Destructor
168                                                   110 
169 void G4Trd::MakePlanes()                       << 111 G4Trd::~G4Trd()
170 {                                                 112 {
171   G4double dx = fDx1 - fDx2;                   << 
172   G4double dy = fDy1 - fDy2;                   << 
173   G4double dz = 2*fDz;                         << 
174   fHx = std::sqrt(dy*dy + dz*dz);              << 
175   fHy = std::sqrt(dx*dx + dz*dz);              << 
176                                                << 
177   // Set X planes at -Y & +Y                   << 
178   //                                           << 
179   fPlanes[0].a =  0.;                          << 
180   fPlanes[0].b = -dz/fHx;                      << 
181   fPlanes[0].c =  dy/fHx;                      << 
182   fPlanes[0].d = fPlanes[0].b*fDy1 + fPlanes[0 << 
183                                                << 
184   fPlanes[1].a =  fPlanes[0].a;                << 
185   fPlanes[1].b = -fPlanes[0].b;                << 
186   fPlanes[1].c =  fPlanes[0].c;                << 
187   fPlanes[1].d =  fPlanes[0].d;                << 
188                                                << 
189   // Set Y planes at -X & +X                   << 
190   //                                           << 
191   fPlanes[2].a = -dz/fHy;                      << 
192   fPlanes[2].b =  0.;                          << 
193   fPlanes[2].c =  dx/fHy;                      << 
194   fPlanes[2].d = fPlanes[2].a*fDx1 + fPlanes[2 << 
195                                                << 
196   fPlanes[3].a = -fPlanes[2].a;                << 
197   fPlanes[3].b =  fPlanes[2].b;                << 
198   fPlanes[3].c =  fPlanes[2].c;                << 
199   fPlanes[3].d =  fPlanes[2].d;                << 
200 }                                                 113 }
201                                                   114 
202 ////////////////////////////////////////////// << 115 ////////////////////////////////////////////////////////////////////////////
203 //                                                116 //
204 // Get volume                                  << 
205                                                << 
206 G4double G4Trd::GetCubicVolume()               << 
207 {                                              << 
208   if (fCubicVolume == 0.)                      << 
209   {                                            << 
210     fCubicVolume = 2*fDz*( (fDx1+fDx2)*(fDy1+f << 
211                            (fDx2-fDx1)*(fDy2-f << 
212   }                                            << 
213   return fCubicVolume;                         << 
214 }                                              << 
215                                                << 
216 ////////////////////////////////////////////// << 
217 //                                                117 //
218 // Get surface area                            << 
219                                                   118 
220 G4double G4Trd::GetSurfaceArea()               << 119 void G4Trd::SetAllParameters ( G4double pdx1, G4double pdx2, G4double pdy1, 
                                                   >> 120                                G4double pdy2, G4double pdz ) 
221 {                                                 121 {
222   if (fSurfaceArea == 0.)                      << 122   CheckAndSetAllParameters (pdx1, pdx2, pdy1, pdy2, pdz);
223   {                                            << 
224     fSurfaceArea =                             << 
225       4*(fDx1*fDy1 + fDx2*fDy2) + 2*(fDx1+fDx2 << 
226   }                                            << 
227   return fSurfaceArea;                         << 
228 }                                                 123 }
229                                                   124 
230 ////////////////////////////////////////////// << 125 
                                                   >> 126 /////////////////////////////////////////////////////////////////////////
231 //                                                127 //
232 // Dispatch to parameterisation for replicatio    128 // Dispatch to parameterisation for replication mechanism dimension
233 // computation & modification                  << 129 // computation & modification.
234                                                   130 
235 void G4Trd::ComputeDimensions(       G4VPVPara    131 void G4Trd::ComputeDimensions(       G4VPVParameterisation* p,
236                                const G4int n,     132                                const G4int n,
237                                const G4VPhysic    133                                const G4VPhysicalVolume* pRep )
238 {                                                 134 {
239   p->ComputeDimensions(*this,n,pRep);             135   p->ComputeDimensions(*this,n,pRep);
240 }                                                 136 }
241                                                   137 
242 ////////////////////////////////////////////// << 
243 //                                             << 
244 // Get bounding box                            << 
245                                                << 
246 void G4Trd::BoundingLimits(G4ThreeVector& pMin << 
247 {                                              << 
248   G4double dx1 = GetXHalfLength1();            << 
249   G4double dx2 = GetXHalfLength2();            << 
250   G4double dy1 = GetYHalfLength1();            << 
251   G4double dy2 = GetYHalfLength2();            << 
252   G4double dz  = GetZHalfLength();             << 
253                                                << 
254   G4double xmax = std::max(dx1,dx2);           << 
255   G4double ymax = std::max(dy1,dy2);           << 
256   pMin.set(-xmax,-ymax,-dz);                   << 
257   pMax.set( xmax, ymax, dz);                   << 
258                                                << 
259   // Check correctness of the bounding box     << 
260   //                                           << 
261   if (pMin.x() >= pMax.x() || pMin.y() >= pMax << 
262   {                                            << 
263     std::ostringstream message;                << 
264     message << "Bad bounding box (min >= max)  << 
265             << GetName() << " !"               << 
266             << "\npMin = " << pMin             << 
267             << "\npMax = " << pMax;            << 
268     G4Exception("G4Trd::BoundingLimits()", "Ge << 
269     DumpInfo();                                << 
270   }                                            << 
271 }                                              << 
272                                                   138 
273 ////////////////////////////////////////////// << 139 ///////////////////////////////////////////////////////////////////////////
274 //                                                140 //
275 // Calculate extent under transform and specif    141 // Calculate extent under transform and specified limit
276                                                   142 
277 G4bool G4Trd::CalculateExtent( const EAxis pAx    143 G4bool G4Trd::CalculateExtent( const EAxis pAxis,
278                                const G4VoxelLi    144                                const G4VoxelLimits& pVoxelLimit,
279                                const G4AffineT    145                                const G4AffineTransform& pTransform,
280                                      G4double&    146                                      G4double& pMin, G4double& pMax ) const
281 {                                                 147 {
282   G4ThreeVector bmin, bmax;                    << 148   if (!pTransform.IsRotated())
283   G4bool exist;                                << 
284                                                << 
285   // Check bounding box (bbox)                 << 
286   //                                           << 
287   BoundingLimits(bmin,bmax);                   << 
288   G4BoundingEnvelope bbox(bmin,bmax);          << 
289 #ifdef G4BBOX_EXTENT                           << 
290   return bbox.CalculateExtent(pAxis,pVoxelLimi << 
291 #endif                                         << 
292   if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox << 
293   {                                               149   {
294     return exist = pMin < pMax;                << 150     // Special case handling for unrotated solids
                                                   >> 151     // Compute x/y/z mins and maxs respecting limits, with early returns
                                                   >> 152     // if outside limits. Then switch() on pAxis
                                                   >> 153 
                                                   >> 154     G4double xoffset,xMin,xMax;
                                                   >> 155     G4double yoffset,yMin,yMax;
                                                   >> 156     G4double zoffset,zMin,zMax;
                                                   >> 157 
                                                   >> 158     zoffset=pTransform.NetTranslation().z();
                                                   >> 159     zMin=zoffset-fDz;
                                                   >> 160     zMax=zoffset+fDz;
                                                   >> 161     if (pVoxelLimit.IsZLimited())
                                                   >> 162     {
                                                   >> 163       if ( (zMin>pVoxelLimit.GetMaxZExtent()+kCarTolerance)
                                                   >> 164         || (zMax<pVoxelLimit.GetMinZExtent()-kCarTolerance) )
                                                   >> 165       {
                                                   >> 166         return false;
                                                   >> 167       }
                                                   >> 168         else
                                                   >> 169       {
                                                   >> 170         if (zMin<pVoxelLimit.GetMinZExtent())
                                                   >> 171         {
                                                   >> 172           zMin=pVoxelLimit.GetMinZExtent();
                                                   >> 173         }
                                                   >> 174         if (zMax>pVoxelLimit.GetMaxZExtent())
                                                   >> 175         {
                                                   >> 176           zMax=pVoxelLimit.GetMaxZExtent();
                                                   >> 177         }
                                                   >> 178       }
                                                   >> 179     }
                                                   >> 180     xoffset=pTransform.NetTranslation().x();
                                                   >> 181     if (fDx2 >= fDx1)
                                                   >> 182     { 
                                                   >> 183       xMax =  xoffset+(fDx1+fDx2)/2+(zMax-zoffset)*(fDx2-fDx1)/(2*fDz) ;
                                                   >> 184       xMin = 2*xoffset - xMax ;
                                                   >> 185     }
                                                   >> 186     else
                                                   >> 187     {
                                                   >> 188       xMax =  xoffset+(fDx1+fDx2)/2+(zMin-zoffset)*(fDx2-fDx1)/(2*fDz) ;
                                                   >> 189       xMin =  2*xoffset - xMax ;
                                                   >> 190     }   
                                                   >> 191     if (pVoxelLimit.IsXLimited())
                                                   >> 192     {
                                                   >> 193       if ( (xMin>pVoxelLimit.GetMaxXExtent()+kCarTolerance)
                                                   >> 194         || (xMax<pVoxelLimit.GetMinXExtent()-kCarTolerance) )
                                                   >> 195       {
                                                   >> 196         return false;
                                                   >> 197       }
                                                   >> 198       else
                                                   >> 199       {
                                                   >> 200         if (xMin<pVoxelLimit.GetMinXExtent())
                                                   >> 201         {
                                                   >> 202           xMin=pVoxelLimit.GetMinXExtent();
                                                   >> 203         }
                                                   >> 204         if (xMax>pVoxelLimit.GetMaxXExtent())
                                                   >> 205         {
                                                   >> 206           xMax=pVoxelLimit.GetMaxXExtent();
                                                   >> 207         }
                                                   >> 208       }
                                                   >> 209     }
                                                   >> 210     yoffset= pTransform.NetTranslation().y() ;
                                                   >> 211     if(fDy2 >= fDy1)
                                                   >> 212     {
                                                   >> 213       yMax = yoffset+(fDy2+fDy1)/2+(zMax-zoffset)*(fDy2-fDy1)/(2*fDz) ;
                                                   >> 214       yMin = 2*yoffset - yMax ;
                                                   >> 215     }
                                                   >> 216     else
                                                   >> 217     {
                                                   >> 218       yMax = yoffset+(fDy2+fDy1)/2+(zMin-zoffset)*(fDy2-fDy1)/(2*fDz) ;
                                                   >> 219       yMin = 2*yoffset - yMax ;  
                                                   >> 220     }  
                                                   >> 221     if (pVoxelLimit.IsYLimited())
                                                   >> 222     {
                                                   >> 223       if ( (yMin>pVoxelLimit.GetMaxYExtent()+kCarTolerance)
                                                   >> 224         || (yMax<pVoxelLimit.GetMinYExtent()-kCarTolerance) )
                                                   >> 225       {
                                                   >> 226         return false;
                                                   >> 227       }
                                                   >> 228       else
                                                   >> 229       {
                                                   >> 230         if (yMin<pVoxelLimit.GetMinYExtent())
                                                   >> 231         {
                                                   >> 232           yMin=pVoxelLimit.GetMinYExtent();
                                                   >> 233         }
                                                   >> 234         if (yMax>pVoxelLimit.GetMaxYExtent())
                                                   >> 235         {
                                                   >> 236           yMax=pVoxelLimit.GetMaxYExtent();
                                                   >> 237         }
                                                   >> 238       }
                                                   >> 239     }
                                                   >> 240 
                                                   >> 241     switch (pAxis)
                                                   >> 242     {
                                                   >> 243       case kXAxis:
                                                   >> 244         pMin=xMin;
                                                   >> 245         pMax=xMax;
                                                   >> 246         break;
                                                   >> 247       case kYAxis:
                                                   >> 248         pMin=yMin;
                                                   >> 249         pMax=yMax;
                                                   >> 250         break;
                                                   >> 251       case kZAxis:
                                                   >> 252         pMin=zMin;
                                                   >> 253         pMax=zMax;
                                                   >> 254         break;
                                                   >> 255       default:
                                                   >> 256         break;
                                                   >> 257     }
                                                   >> 258 
                                                   >> 259     // Add 2*Tolerance to avoid precision troubles ?
                                                   >> 260     //
                                                   >> 261     pMin-=kCarTolerance;
                                                   >> 262     pMax+=kCarTolerance;
                                                   >> 263 
                                                   >> 264     return true;
295   }                                               265   }
                                                   >> 266   else
                                                   >> 267   {
                                                   >> 268     // General rotated case - create and clip mesh to boundaries
296                                                   269 
297   // Set bounding envelope (benv) and calculat << 270     G4bool existsAfterClip=false;
298   //                                           << 271     G4ThreeVectorList *vertices;
299   G4double dx1 = GetXHalfLength1();            << 272 
300   G4double dx2 = GetXHalfLength2();            << 273     pMin=+kInfinity;
301   G4double dy1 = GetYHalfLength1();            << 274     pMax=-kInfinity;
302   G4double dy2 = GetYHalfLength2();            << 275 
303   G4double dz  = GetZHalfLength();             << 276     // Calculate rotated vertex coordinates
304                                                << 277     //
305   G4ThreeVectorList baseA(4), baseB(4);        << 278     vertices=CreateRotatedVertices(pTransform);
306   baseA[0].set(-dx1,-dy1,-dz);                 << 279     ClipCrossSection(vertices,0,pVoxelLimit,pAxis,pMin,pMax);
307   baseA[1].set( dx1,-dy1,-dz);                 << 280     ClipCrossSection(vertices,4,pVoxelLimit,pAxis,pMin,pMax);
308   baseA[2].set( dx1, dy1,-dz);                 << 281     ClipBetweenSections(vertices,0,pVoxelLimit,pAxis,pMin,pMax);
309   baseA[3].set(-dx1, dy1,-dz);                 << 282       
310   baseB[0].set(-dx2,-dy2, dz);                 << 283     if (pMin!=kInfinity||pMax!=-kInfinity)
311   baseB[1].set( dx2,-dy2, dz);                 << 284     {
312   baseB[2].set( dx2, dy2, dz);                 << 285       existsAfterClip=true;
313   baseB[3].set(-dx2, dy2, dz);                 << 286 
314                                                << 287       // Add 2*tolerance to avoid precision troubles
315   std::vector<const G4ThreeVectorList *> polyg << 288       //
316   polygons[0] = &baseA;                        << 289       pMin-=kCarTolerance;
317   polygons[1] = &baseB;                        << 290       pMax+=kCarTolerance;
318                                                << 291         
319   G4BoundingEnvelope benv(bmin,bmax,polygons); << 292     }
320   exist = benv.CalculateExtent(pAxis,pVoxelLim << 293     else
321   return exist;                                << 294     {
                                                   >> 295       // Check for case where completely enveloping clipping volume
                                                   >> 296       // If point inside then we are confident that the solid completely
                                                   >> 297       // envelopes the clipping volume. Hence set min/max extents according
                                                   >> 298       // to clipping volume extents along the specified axis.
                                                   >> 299 
                                                   >> 300       G4ThreeVector clipCentre(
                                                   >> 301          (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
                                                   >> 302          (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
                                                   >> 303          (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5);
                                                   >> 304         
                                                   >> 305       if (Inside(pTransform.Inverse().TransformPoint(clipCentre))!=kOutside)
                                                   >> 306       {
                                                   >> 307         existsAfterClip=true;
                                                   >> 308         pMin=pVoxelLimit.GetMinExtent(pAxis);
                                                   >> 309         pMax=pVoxelLimit.GetMaxExtent(pAxis);
                                                   >> 310       }
                                                   >> 311     }
                                                   >> 312     delete vertices;
                                                   >> 313     return existsAfterClip;
                                                   >> 314   }
322 }                                                 315 }
323                                                   316 
324 ////////////////////////////////////////////// << 317 ///////////////////////////////////////////////////////////////////
325 //                                                318 //
326 // Return whether point inside/outside/on surf    319 // Return whether point inside/outside/on surface, using tolerance
327                                                   320 
328 EInside G4Trd::Inside( const G4ThreeVector& p     321 EInside G4Trd::Inside( const G4ThreeVector& p ) const
329 {                                              << 322 {  
330   G4double dx = fPlanes[3].a*std::abs(p.x())+f << 323   EInside in=kOutside;
331   G4double dy = fPlanes[1].b*std::abs(p.y())+f << 324   G4double x,y,zbase1,zbase2;
332   G4double dxy = std::max(dx,dy);              << 325     
                                                   >> 326   if (fabs(p.z())<=fDz-kCarTolerance/2)
                                                   >> 327   {
                                                   >> 328     zbase1=p.z()+fDz;  // Dist from -ve z plane
                                                   >> 329     zbase2=fDz-p.z();  // Dist from +ve z plane
333                                                   330 
334   G4double dz = std::abs(p.z())-fDz;           << 331     // Check whether inside x tolerance
335   G4double dist = std::max(dz,dxy);            << 332     //
                                                   >> 333     x=0.5*(fDx2*zbase1+fDx1*zbase2)/fDz - kCarTolerance/2;
                                                   >> 334     if (fabs(p.x())<=x)
                                                   >> 335     {
                                                   >> 336       y=0.5*((fDy2*zbase1+fDy1*zbase2))/fDz - kCarTolerance/2;
                                                   >> 337       if (fabs(p.y())<=y)
                                                   >> 338       {
                                                   >> 339         in=kInside;
                                                   >> 340       }
                                                   >> 341       else if (fabs(p.y())<=y+kCarTolerance)
                                                   >> 342       {
                                                   >> 343         in=kSurface;
                                                   >> 344       }
                                                   >> 345     }
                                                   >> 346     else if (fabs(p.x())<=x+kCarTolerance)
                                                   >> 347     {
                                                   >> 348       // y = y half width of shape at z of point + tolerant boundary
                                                   >> 349       //
                                                   >> 350       y=0.5*((fDy2*zbase1+fDy1*zbase2))/fDz + kCarTolerance/2;
                                                   >> 351       if (fabs(p.y())<=y)
                                                   >> 352       {
                                                   >> 353         in=kSurface;
                                                   >> 354       }
                                                   >> 355     }
                                                   >> 356   }
                                                   >> 357   else if (fabs(p.z())<=fDz+kCarTolerance/2)
                                                   >> 358   {
                                                   >> 359     // Only need to check outer tolerant boundaries
                                                   >> 360     //
                                                   >> 361     zbase1=p.z()+fDz;  // Dist from -ve z plane
                                                   >> 362     zbase2=fDz-p.z();   // Dist from +ve z plane
336                                                   363 
337   return (dist > halfCarTolerance) ? kOutside  << 364     // x = x half width of shape at z of point plus tolerance
338     ((dist > -halfCarTolerance) ? kSurface : k << 365     //
                                                   >> 366     x=0.5*(fDx2*zbase1+fDx1*zbase2)/fDz + kCarTolerance/2;
                                                   >> 367     if (fabs(p.x())<=x)
                                                   >> 368     {
                                                   >> 369       // y = y half width of shape at z of point
                                                   >> 370       //
                                                   >> 371       y=0.5*((fDy2*zbase1+fDy1*zbase2))/fDz + kCarTolerance/2;
                                                   >> 372       if (fabs(p.y())<=y) in=kSurface;
                                                   >> 373     }
                                                   >> 374   }  
                                                   >> 375   return in;
339 }                                                 376 }
340                                                   377 
341 //////////////////////////////////////////////    378 //////////////////////////////////////////////////////////////////////////
342 //                                                379 //
343 // Determine side where point is, and return c << 380 // Calculate side nearest to p, and return normal
                                                   >> 381 // If two sides are equidistant, normal of first side (x/y/z) 
                                                   >> 382 // encountered returned
344                                                   383 
345 G4ThreeVector G4Trd::SurfaceNormal( const G4Th    384 G4ThreeVector G4Trd::SurfaceNormal( const G4ThreeVector& p ) const
346 {                                                 385 {
347   G4int nsurf = 0; // number of surfaces where << 386   G4ThreeVector norm;
                                                   >> 387   G4double z,tanx,secx,newpx,widx;
                                                   >> 388   G4double tany,secy,newpy,widy;
                                                   >> 389   G4double distx,disty,distz,fcos;
                                                   >> 390 
                                                   >> 391   z=2.0*fDz;
                                                   >> 392 
                                                   >> 393   tanx=(fDx2-fDx1)/z;
                                                   >> 394   secx=sqrt(1.0+tanx*tanx);
                                                   >> 395   newpx=fabs(p.x())-p.z()*tanx;
                                                   >> 396   widx=fDx2-fDz*tanx;
                                                   >> 397 
                                                   >> 398   tany=(fDy2-fDy1)/z;
                                                   >> 399   secy=sqrt(1.0+tany*tany);
                                                   >> 400   newpy=fabs(p.y())-p.z()*tany;
                                                   >> 401   widy=fDy2-fDz*tany;
                                                   >> 402 
                                                   >> 403   distx=fabs(newpx-widx)/secx;  // perpendicular distance to x side
                                                   >> 404   disty=fabs(newpy-widy)/secy;  //                        to y side
                                                   >> 405   distz=fabs(fabs(p.z())-fDz);  //                        to z side
348                                                   406 
349   // Check Z faces                             << 407   // find closest side
350   //                                              408   //
351   G4double nz = 0;                             << 409   if (distx<=disty)
352   G4double dz = std::abs(p.z()) - fDz;         << 410   { 
353   if (std::abs(dz) <= halfCarTolerance)        << 411     if (distx<=distz) 
354   {                                            << 412     {
355     nz = (p.z() < 0) ? -1 : 1;                 << 413       // Closest to X
356     ++nsurf;                                   << 414       //
                                                   >> 415       fcos=1.0/secx;
                                                   >> 416       // normal=(+/-cos(ang),0,-sin(ang))
                                                   >> 417       if (p.x()>=0)
                                                   >> 418         norm=G4ThreeVector(fcos,0,-tanx*fcos);
                                                   >> 419       else
                                                   >> 420         norm=G4ThreeVector(-fcos,0,-tanx*fcos);
                                                   >> 421     }
                                                   >> 422     else
                                                   >> 423     {
                                                   >> 424       // Closest to Z
                                                   >> 425       //
                                                   >> 426       if (p.z()>=0)
                                                   >> 427         norm=G4ThreeVector(0,0,1);
                                                   >> 428       else
                                                   >> 429         norm=G4ThreeVector(0,0,-1);
                                                   >> 430     }
                                                   >> 431   }
                                                   >> 432   else
                                                   >> 433   {  
                                                   >> 434     if (disty<=distz)
                                                   >> 435     {
                                                   >> 436       // Closest to Y
                                                   >> 437       //
                                                   >> 438       fcos=1.0/secy;
                                                   >> 439       if (p.y()>=0)
                                                   >> 440         norm=G4ThreeVector(0,fcos,-tany*fcos);
                                                   >> 441       else
                                                   >> 442         norm=G4ThreeVector(0,-fcos,-tany*fcos);
                                                   >> 443     }
                                                   >> 444     else 
                                                   >> 445     {
                                                   >> 446       // Closest to Z
                                                   >> 447       //
                                                   >> 448       if (p.z()>=0)
                                                   >> 449         norm=G4ThreeVector(0,0,1);
                                                   >> 450       else
                                                   >> 451         norm=G4ThreeVector(0,0,-1);
                                                   >> 452     }
357   }                                               453   }
                                                   >> 454   return norm;   
                                                   >> 455 }
358                                                   456 
359   // Check Y faces                             << 457 ////////////////////////////////////////////////////////////////////////////
360   //                                           << 458 //
361   G4double ny = 0;                             << 459 // Calculate distance to shape from outside
362   G4double dy1 = fPlanes[0].b*p.y();           << 460 // - return kInfinity if no intersection
363   G4double dy2 = fPlanes[0].c*p.z() + fPlanes[ << 461 //
364   if (std::abs(dy2 + dy1) <= halfCarTolerance) << 462 // ALGORITHM:
                                                   >> 463 // For each component, calculate pair of minimum and maximum intersection
                                                   >> 464 // values for which the particle is in the extent of the shape
                                                   >> 465 // - The smallest (MAX minimum) allowed distance of the pairs is intersect
                                                   >> 466 // - Z plane intersectin uses tolerance
                                                   >> 467 // - XZ YZ planes use logic & *SLIGHTLY INCORRECT* tolerance
                                                   >> 468 //   (this saves at least 1 sqrt, 1 multiply and 1 divide... in applicable
                                                   >> 469 //    cases)
                                                   >> 470 // - Note: XZ and YZ planes each divide space into four regions,
                                                   >> 471 //   characterised by ss1 ss2
                                                   >> 472 // NOTE:
                                                   >> 473 //
                                                   >> 474 // `Inside' safe - meaningful answers given if point is inside the exact
                                                   >> 475 // shape.
                                                   >> 476 
                                                   >> 477 G4double G4Trd::DistanceToIn( const G4ThreeVector& p,
                                                   >> 478                               const G4ThreeVector& v ) const
                                                   >> 479 {  
                                                   >> 480   G4double snxt = kInfinity ;    // snxt = default return value
                                                   >> 481   G4double smin,smax;
                                                   >> 482   G4double s1,s2,tanxz,tanyz,ds1,ds2;
                                                   >> 483   G4double ss1,ss2,sn1=0.,sn2=0.,Dist;
                                                   >> 484 
                                                   >> 485   if ( v.z() )  // Calculate valid z intersect range
365   {                                               486   {
366     ny += fPlanes[0].b;                        << 487     if ( v.z() > 0 )   // Calculate smax: must be +ve or no intersection.
367     nz += fPlanes[0].c;                        << 488     {
368     ++nsurf;                                   << 489       Dist = fDz - p.z() ;  // to plane at +dz
                                                   >> 490 
                                                   >> 491       if (Dist >= 0.5*kCarTolerance)
                                                   >> 492       {
                                                   >> 493         smax = Dist/v.z() ;
                                                   >> 494         smin = -(fDz + p.z())/v.z() ;
                                                   >> 495       }
                                                   >> 496       else  return snxt ;
                                                   >> 497     }
                                                   >> 498     else // v.z <0
                                                   >> 499     {
                                                   >> 500       Dist=fDz+p.z();  // plane at -dz
                                                   >> 501 
                                                   >> 502       if ( Dist >= 0.5*kCarTolerance )
                                                   >> 503       {
                                                   >> 504         smax = -Dist/v.z() ;
                                                   >> 505         smin = (fDz - p.z())/v.z() ;
                                                   >> 506       }
                                                   >> 507       else return snxt ; 
                                                   >> 508     }
                                                   >> 509     if (smin < 0 ) smin = 0 ;
369   }                                               510   }
370   if (std::abs(dy2 - dy1) <= halfCarTolerance) << 511   else // v.z=0
371   {                                               512   {
372     ny += fPlanes[1].b;                        << 513     if (fabs(p.z()) >= fDz ) return snxt ;     // Outside & no intersect
373     nz += fPlanes[1].c;                        << 514     else
374     ++nsurf;                                   << 515     {
                                                   >> 516       smin = 0 ;    // Always inside z range
                                                   >> 517       smax = kInfinity;
                                                   >> 518     }
375   }                                               519   }
376                                                   520 
377   // Check X faces                             << 521   // Calculate x intersection range
378   //                                              522   //
379   G4double nx = 0;                             << 523   // Calc half width at p.z, and components towards planes
380   G4double dx1 = fPlanes[2].a*p.x();           << 524 
381   G4double dx2 = fPlanes[2].c*p.z() + fPlanes[ << 525   tanxz = (fDx2 - fDx1)*0.5/fDz ;
382   if (std::abs(dx2 + dx1) <= halfCarTolerance) << 526   s1    = 0.5*(fDx1+fDx2) + tanxz*p.z() ;  // x half width at p.z
                                                   >> 527   ds1   = v.x() - tanxz*v.z() ;       // Components of v towards faces at +-x
                                                   >> 528   ds2   = v.x() + tanxz*v.z() ;
                                                   >> 529   ss1   = s1 - p.x() ;         // -delta x to +ve plane
                                                   >> 530                                // -ve when outside
                                                   >> 531   ss2   = -s1 - p.x() ;        // -delta x to -ve plane
                                                   >> 532                                // +ve when outside
                                                   >> 533     
                                                   >> 534   if (ss1 < 0 && ss2 <= 0 )
383   {                                               535   {
384     nx += fPlanes[2].a;                        << 536     if (ds1 < 0)   // In +ve coord Area
385     nz += fPlanes[2].c;                        << 537     {
386     ++nsurf;                                   << 538       sn1 = ss1/ds1 ;
                                                   >> 539 
                                                   >> 540       if ( ds2 < 0 ) sn2 = ss2/ds2 ;           
                                                   >> 541       else           sn2 = kInfinity ;
                                                   >> 542     }
                                                   >> 543     else return snxt ;
387   }                                               544   }
388   if (std::abs(dx2 - dx1) <= halfCarTolerance) << 545   else if ( ss1 >= 0 && ss2 > 0 )
389   {                                               546   {
390     nx += fPlanes[3].a;                        << 547     if ( ds2 > 0 )  // In -ve coord Area
391     nz += fPlanes[3].c;                        << 548     {
392     ++nsurf;                                   << 549       sn1 = ss2/ds2 ;
                                                   >> 550 
                                                   >> 551       if (ds1 > 0)  sn2 = ss1/ds1 ;      
                                                   >> 552       else          sn2 = kInfinity;      
                                                   >> 553         
                                                   >> 554     }
                                                   >> 555     else   return snxt ;
393   }                                               556   }
                                                   >> 557   else if (ss1 >= 0 && ss2 <= 0 )
                                                   >> 558   {
                                                   >> 559     // Inside Area - calculate leaving distance
                                                   >> 560     // *Don't* use exact distance to side for tolerance
                                                   >> 561     //                                             = ss1*cos(ang xz)
                                                   >> 562     //                                             = ss1/sqrt(1.0+tanxz*tanxz)
                                                   >> 563     sn1 = 0 ;
394                                                   564 
395   // Return normal                             << 565     if ( ds1 > 0 )
396   //                                           << 566     {
397   if (nsurf == 1)      return {nx,ny,nz};      << 567       if (ss1 > 0.5*kCarTolerance) sn2 = ss1/ds1 ; // Leave +ve side extent
398   else if (nsurf != 0) return G4ThreeVector(nx << 568       else                         return snxt ;   // Leave immediately by +ve 
399   else                                         << 569     }
                                                   >> 570     else  sn2 = kInfinity ;
                                                   >> 571       
                                                   >> 572     if ( ds2 < 0 )
                                                   >> 573     {
                                                   >> 574       if ( ss2 < -0.5*kCarTolerance )
                                                   >> 575       {
                                                   >> 576         Dist = ss2/ds2 ;            // Leave -ve side extent
                                                   >> 577         if ( Dist < sn2 ) sn2 = Dist ;
                                                   >> 578       }
                                                   >> 579       else  return snxt ;
                                                   >> 580     }    
                                                   >> 581   }
                                                   >> 582   else if (ss1 < 0 && ss2 > 0 )
400   {                                               583   {
401     // Point is not on the surface             << 584     // Within +/- plane cross-over areas (not on boundaries ss1||ss2==0)
402     //                                         << 585 
403 #ifdef G4CSGDEBUG                              << 586     if ( ds1 >= 0 || ds2 <= 0 )
404     std::ostringstream message;                << 587     {   
405     G4long oldprc = message.precision(16);     << 588       return snxt ;
406     message << "Point p is not on surface (!?) << 589     }
407             << GetName() << G4endl;            << 590     else  // Will intersect & stay inside
408     message << "Position:\n";                  << 591     {
409     message << "   p.x() = " << p.x()/mm << "  << 592       sn1  = ss1/ds1 ;
410     message << "   p.y() = " << p.y()/mm << "  << 593       Dist = ss2/ds2 ;
411     message << "   p.z() = " << p.z()/mm << "  << 594       if (Dist > sn1 ) sn1 = Dist ;
412     G4cout.precision(oldprc) ;                 << 595       sn2 = kInfinity ;
413     G4Exception("G4Trd::SurfaceNormal(p)", "Ge << 596     }
414                 JustWarning, message );        << 
415     DumpInfo();                                << 
416 #endif                                         << 
417     return ApproxSurfaceNormal(p);             << 
418   }                                               597   }
419 }                                              << 
420                                                   598 
421 ////////////////////////////////////////////// << 599   // Reduce allowed range of distances as appropriate
422 //                                             << 
423 // Algorithm for SurfaceNormal() following the << 
424 // for points not on the surface               << 
425                                                   600 
426 G4ThreeVector G4Trd::ApproxSurfaceNormal( cons << 601   if ( sn1 > smin ) smin = sn1 ;
427 {                                              << 602   if ( sn2 < smax ) smax = sn2 ;
428   G4double dist = -DBL_MAX;                    << 
429   G4int iside = 0;                             << 
430   for (G4int i=0; i<4; ++i)                    << 
431   {                                            << 
432     G4double d = fPlanes[i].a*p.x() +          << 
433                  fPlanes[i].b*p.y() +          << 
434                  fPlanes[i].c*p.z() + fPlanes[ << 
435     if (d > dist) { dist = d; iside = i; }     << 
436   }                                            << 
437                                                   603 
438   G4double distz = std::abs(p.z()) - fDz;      << 604   // Check for incompatible ranges (eg z intersects between 50 ->100 and x
439   if (dist > distz)                            << 605   // only 10-40 -> no intersection)
440     return { fPlanes[iside].a, fPlanes[iside]. << 
441   else                                         << 
442     return { 0, 0, (G4double)((p.z() < 0) ? -1 << 
443 }                                              << 
444                                                   606 
445 ////////////////////////////////////////////// << 607   if ( smax < smin ) return snxt ;
446 //                                             << 
447 // Calculate distance to shape from outside    << 
448 //  - return kInfinity if no intersection      << 
449                                                   608 
450 G4double G4Trd::DistanceToIn(const G4ThreeVect << 609   // Calculate valid y intersection range 
451                              const G4ThreeVect << 610   // (repeat of x intersection code)
452 {                                              << 
453   // Z intersections                           << 
454   //                                           << 
455   if ((std::abs(p.z()) - fDz) >= -halfCarToler << 
456     return kInfinity;                          << 
457   G4double invz = (-v.z() == 0) ? DBL_MAX : -1 << 
458   G4double dz = (invz < 0) ? fDz : -fDz;       << 
459   G4double tzmin = (p.z() + dz)*invz;          << 
460   G4double tzmax = (p.z() - dz)*invz;          << 
461                                                   611 
462   // Y intersections                           << 612   tanyz = (fDy2-fDy1)*0.5/fDz ;
463   //                                           << 613   s2    = 0.5*(fDy1+fDy2) + tanyz*p.z() ;  // y half width at p.z
464   G4double tmin0 = tzmin, tmax0 = tzmax;       << 614   ds1   = v.y() - tanyz*v.z() ;       // Components of v towards faces at +-y
465   G4double ya = fPlanes[0].b*v.y(), yb = fPlan << 615   ds2   = v.y() + tanyz*v.z() ;
466   G4double yc = fPlanes[0].b*p.y(), yd = fPlan << 616   ss1   = s2 - p.y() ;         // -delta y to +ve plane
467   G4double cos0 = yb + ya;                     << 617   ss2   = -s2 - p.y() ;        // -delta y to -ve plane
468   G4double dis0 = yd + yc;                     << 618     
469   if (dis0 >= -halfCarTolerance)               << 619   if ( ss1 < 0 && ss2 <= 0 )
470   {                                               620   {
471     if (cos0 >= 0) return kInfinity;           << 621     if (ds1 < 0 ) // In +ve coord Area
472     G4double tmp  = -dis0/cos0;                << 622     {
473     if (tmin0 < tmp) tmin0 = tmp;              << 623       sn1 = ss1/ds1 ;
474   }                                            << 624       if ( ds2 < 0 )  sn2 = ss2/ds2 ;
475   else if (cos0 > 0)                           << 625       else            sn2 = kInfinity ;
476   {                                            << 626     }
477     G4double tmp  = -dis0/cos0;                << 627     else   return snxt ;
478     if (tmax0 > tmp) tmax0 = tmp;              << 
479   }                                               628   }
480                                                << 629   else if ( ss1 >= 0 && ss2 > 0 )
481   G4double tmin1 = tmin0, tmax1 = tmax0;       << 
482   G4double cos1 = yb - ya;                     << 
483   G4double dis1 = yd - yc;                     << 
484   if (dis1 >= -halfCarTolerance)               << 
485   {                                               630   {
486     if (cos1 >= 0) return kInfinity;           << 631     if ( ds2 > 0 )  // In -ve coord Area
487     G4double tmp  = -dis1/cos1;                << 632     {
488     if (tmin1 < tmp) tmin1 = tmp;              << 633       sn1 = ss2/ds2 ;
                                                   >> 634       if ( ds1 > 0 )  sn2 = ss1/ds1 ;
                                                   >> 635       else            sn2 = kInfinity ;      
                                                   >> 636     }
                                                   >> 637     else   return snxt ;
489   }                                               638   }
490   else if (cos1 > 0)                           << 639   else if (ss1 >= 0 && ss2 <= 0 )
491   {                                               640   {
492     G4double tmp  = -dis1/cos1;                << 641     // Inside Area - calculate leaving distance
493     if (tmax1 > tmp) tmax1 = tmp;              << 642     // *Don't* use exact distance to side for tolerance
494   }                                            << 643     //                                          = ss1*cos(ang yz)
                                                   >> 644     //                                          = ss1/sqrt(1.0+tanyz*tanyz)
                                                   >> 645     sn1 = 0 ;
495                                                   646 
496   // X intersections                           << 647     if ( ds1 > 0 )
497   //                                           << 648     {
498   G4double tmin2 = tmin1, tmax2 = tmax1;       << 649       if (ss1 > 0.5*kCarTolerance) sn2 = ss1/ds1 ; // Leave +ve side extent
499   G4double xa = fPlanes[2].a*v.x(), xb = fPlan << 650       else                         return snxt ;   // Leave immediately by +ve
500   G4double xc = fPlanes[2].a*p.x(), xd = fPlan << 651     }
501   G4double cos2 = xb + xa;                     << 652     else  sn2 = kInfinity ;
502   G4double dis2 = xd + xc;                     << 653       
503   if (dis2 >= -halfCarTolerance)               << 654     if ( ds2 < 0 )
504   {                                            << 655     {
505     if (cos2 >= 0) return kInfinity;           << 656       if ( ss2 < -0.5*kCarTolerance )
506     G4double tmp  = -dis2/cos2;                << 657       {
507     if (tmin2 < tmp) tmin2 = tmp;              << 658         Dist = ss2/ds2 ; // Leave -ve side extent
                                                   >> 659         if (Dist < sn2) sn2=Dist;
                                                   >> 660       }
                                                   >> 661       else  return snxt ;
                                                   >> 662     }    
508   }                                               663   }
509   else if (cos2 > 0)                           << 664   else if (ss1 < 0 && ss2 > 0 )
510   {                                               665   {
511     G4double tmp  = -dis2/cos2;                << 666     // Within +/- plane cross-over areas (not on boundaries ss1||ss2==0)
512     if (tmax2 > tmp) tmax2 = tmp;              << 
513   }                                            << 
514                                                   667 
515   G4double tmin3 = tmin2, tmax3 = tmax2;       << 668     if (ds1 >= 0 || ds2 <= 0 )  
516   G4double cos3 = xb - xa;                     << 669     {
517   G4double dis3 = xd - xc;                     << 670       return snxt ;
518   if (dis3 >= -halfCarTolerance)               << 671     }
519   {                                            << 672     else  // Will intersect & stay inside
520     if (cos3 >= 0) return kInfinity;           << 673     {
521     G4double tmp  = -dis3/cos3;                << 674       sn1 = ss1/ds1 ;
522     if (tmin3 < tmp) tmin3 = tmp;              << 675       Dist = ss2/ds2 ;
523   }                                            << 676       if (Dist > sn1 ) sn1 = Dist ;
524   else if (cos3 > 0)                           << 677       sn2 = kInfinity ;
525   {                                            << 678     }
526     G4double tmp  = -dis3/cos3;                << 
527     if (tmax3 > tmp) tmax3 = tmp;              << 
528   }                                               679   }
                                                   >> 680   
                                                   >> 681   // Reduce allowed range of distances as appropriate
529                                                   682 
530   // Find distance                             << 683   if ( sn1 > smin) smin = sn1 ;
531   //                                           << 684   if ( sn2 < smax) smax = sn2 ;
532   G4double tmin = tmin3, tmax = tmax3;         << 685 
533   if (tmax <= tmin + halfCarTolerance) return  << 686   // Check for incompatible ranges (eg x intersects between 50 ->100 and y
534   return (tmin < halfCarTolerance ) ? 0. : tmi << 687   // only 10-40 -> no intersection). Set snxt if ok
                                                   >> 688 
                                                   >> 689   if ( smax > smin ) snxt = smin ;
                                                   >> 690   if (snxt < 0.5*kCarTolerance ) snxt = 0.0 ;
                                                   >> 691 
                                                   >> 692   return snxt ;
535 }                                                 693 }
536                                                   694 
537 ////////////////////////////////////////////// << 695 /////////////////////////////////////////////////////////////////////////
538 //                                                696 //
539 // Calculate exact shortest distance to any bo << 697 // Approximate distance to shape
540 // This is the best fast estimation of the sho << 698 // Calculate perpendicular distances to z/x/y surfaces, return largest
541 // - returns 0 if point is inside              << 699 // which is the most fast estimation of shortest distance to Trd
                                                   >> 700 //  - Safe underestimate
                                                   >> 701 //  - If point within exact shape, return 0 
542                                                   702 
543 G4double G4Trd::DistanceToIn( const G4ThreeVec    703 G4double G4Trd::DistanceToIn( const G4ThreeVector& p ) const
544 {                                                 704 {
545   G4double dx = fPlanes[3].a*std::abs(p.x())+f << 705   G4double safe;
546   G4double dy = fPlanes[1].b*std::abs(p.y())+f << 706   G4double tanxz,distx,safx;
547   G4double dxy = std::max(dx,dy);              << 707   G4double tanyz,disty,safy;
                                                   >> 708   G4double zbase;
                                                   >> 709 
                                                   >> 710   safe=fabs(p.z())-fDz;
                                                   >> 711   if (safe<0) safe=0;      // Also used to ensure x/y distances
                                                   >> 712                            // POSITIVE 
                                                   >> 713 
                                                   >> 714   zbase=fDz+p.z();
548                                                   715 
549   G4double dz = std::abs(p.z())-fDz;           << 716   // Find distance along x direction to closest x plane
550   G4double dist = std::max(dz,dxy);            << 717   //
                                                   >> 718   tanxz=(fDx2-fDx1)*0.5/fDz;
                                                   >> 719   //    widx=fDx1+tanxz*(fDz+p.z()); // x width at p.z
                                                   >> 720   //    distx=fabs(p.x())-widx;      // distance to plane
                                                   >> 721   distx=fabs(p.x())-(fDx1+tanxz*zbase);
                                                   >> 722   if (distx>safe)
                                                   >> 723   {
                                                   >> 724     safx=distx/sqrt(1.0+tanxz*tanxz); // vector Dist=Dist*cos(ang)
                                                   >> 725     if (safx>safe) safe=safx;
                                                   >> 726   }
551                                                   727 
552   return (dist > 0) ? dist : 0.;               << 728   // Find distance along y direction to slanted wall
                                                   >> 729   tanyz=(fDy2-fDy1)*0.5/fDz;
                                                   >> 730   //    widy=fDy1+tanyz*(fDz+p.z()); // y width at p.z
                                                   >> 731   //    disty=fabs(p.y())-widy;      // distance to plane
                                                   >> 732   disty=fabs(p.y())-(fDy1+tanyz*zbase);
                                                   >> 733   if (disty>safe)    
                                                   >> 734   {
                                                   >> 735     safy=disty/sqrt(1.0+tanyz*tanyz); // distance along vector
                                                   >> 736     if (safy>safe) safe=safy;
                                                   >> 737   }
                                                   >> 738   return safe;
553 }                                                 739 }
554                                                   740 
555 ////////////////////////////////////////////// << 741 ////////////////////////////////////////////////////////////////////////
556 //                                                742 //
557 // Calculate distance to surface of shape from << 743 // Calcluate distance to surface of shape from inside
558 // find normal at exit point, if required      << 744 // Calculate distance to x/y/z planes - smallest is exiting distance
559 // - when leaving the surface, return 0        << 745 // - z planes have std. check for tolerance
560                                                << 746 // - xz yz planes have check based on distance || to x or y axis
561 G4double G4Trd::DistanceToOut(const G4ThreeVec << 747 //   (not corrected for slope of planes)
562                               const G4bool cal << 748 // ?BUG? If v.z==0 are there cases when snside not set????
563                                     G4bool* va << 749 
                                                   >> 750 G4double G4Trd::DistanceToOut( const G4ThreeVector& p,
                                                   >> 751                                const G4ThreeVector& v,
                                                   >> 752                                const G4bool calcNorm,
                                                   >> 753                                      G4bool *validNorm,
                                                   >> 754                                      G4ThreeVector *n ) const
564 {                                                 755 {
565   // Z intersections                           << 756   ESide side = kUndefined, snside = kUndefined;
566   //                                           << 757   G4double snxt,pdist;
567   if ((std::abs(p.z()) - fDz) >= -halfCarToler << 758   G4double central,ss1,ss2,ds1,ds2,sn=0.,sn2=0.;
                                                   >> 759   G4double tanxz=0.,cosxz=0.,tanyz=0.,cosyz=0.;
                                                   >> 760 
                                                   >> 761   if (calcNorm) *validNorm=true; // All normals are valid
                                                   >> 762 
                                                   >> 763   // Calculate z plane intersection
                                                   >> 764   if (v.z()>0)
568   {                                               765   {
569     if (calcNorm)                              << 766     pdist=fDz-p.z();
                                                   >> 767     if (pdist>kCarTolerance/2)
570     {                                             768     {
571       *validNorm = true;                       << 769       snxt=pdist/v.z();
572       n->set(0, 0, (p.z() < 0) ? -1 : 1);      << 770       side=kPZ;
                                                   >> 771     }
                                                   >> 772     else
                                                   >> 773     {
                                                   >> 774       if (calcNorm)
                                                   >> 775       {
                                                   >> 776         *n=G4ThreeVector(0,0,1);
                                                   >> 777       }
                                                   >> 778       return snxt=0;
573     }                                             779     }
574     return 0;                                  << 
575   }                                               780   }
576   G4double vz = v.z();                         << 781   else if (v.z()<0) 
577   G4double tmax = (vz == 0) ? DBL_MAX : (std:: << 782   {
578   G4int iside = (vz < 0) ? -4 : -2; // little  << 783     pdist=fDz+p.z();
                                                   >> 784     if (pdist>kCarTolerance/2)
                                                   >> 785     {
                                                   >> 786       snxt=-pdist/v.z();
                                                   >> 787       side=kMZ;
                                                   >> 788     }
                                                   >> 789     else
                                                   >> 790     {
                                                   >> 791       if (calcNorm)
                                                   >> 792       {
                                                   >> 793         *n=G4ThreeVector(0,0,-1);
                                                   >> 794       }
                                                   >> 795       return snxt=0;
                                                   >> 796     }
                                                   >> 797   }
                                                   >> 798   else
                                                   >> 799   {
                                                   >> 800     snxt=kInfinity;
                                                   >> 801   }
579                                                   802 
580   // Y intersections                           << 
581   //                                              803   //
582   G4int i = 0;                                 << 804   // Calculate x intersection
583   for ( ; i<2; ++i)                            << 805   //
584   {                                            << 806   tanxz=(fDx2-fDx1)*0.5/fDz;
585     G4double cosa = fPlanes[i].b*v.y() + fPlan << 807   central=0.5*(fDx1+fDx2);
586     if (cosa > 0)                              << 808 
                                                   >> 809   // +ve plane (1)
                                                   >> 810   //
                                                   >> 811   ss1=central+tanxz*p.z()-p.x();  // distance || x axis to plane
                                                   >> 812                                   // (+ve if point inside)
                                                   >> 813   ds1=v.x()-tanxz*v.z();    // component towards plane at +x
                                                   >> 814                             // (-ve if +ve -> -ve direction)
                                                   >> 815   // -ve plane (2)
                                                   >> 816   //
                                                   >> 817   ss2=-tanxz*p.z()-p.x()-central;  //distance || x axis to plane
                                                   >> 818                                    // (-ve if point inside)
                                                   >> 819   ds2=tanxz*v.z()+v.x();    // component towards plane at -x
                                                   >> 820 
                                                   >> 821   if (ss1>0&&ss2<0)
                                                   >> 822   {
                                                   >> 823     // Normal case - entirely inside region
                                                   >> 824     if (ds1<=0&&ds2<0)
                                                   >> 825     {   
                                                   >> 826       if (ss2<-kCarTolerance/2)
                                                   >> 827       {
                                                   >> 828         sn=ss2/ds2;  // Leave by -ve side
                                                   >> 829         snside=kMX;
                                                   >> 830       }
                                                   >> 831       else
                                                   >> 832       {
                                                   >> 833         sn=0; // Leave immediately by -ve side
                                                   >> 834         snside=kMX;
                                                   >> 835       }
                                                   >> 836     }
                                                   >> 837     else if (ds1>0&&ds2>=0)
                                                   >> 838     {
                                                   >> 839       if (ss1>kCarTolerance/2)
                                                   >> 840       {
                                                   >> 841         sn=ss1/ds1;  // Leave by +ve side
                                                   >> 842         snside=kPX;
                                                   >> 843       }
                                                   >> 844       else
                                                   >> 845       {
                                                   >> 846         sn=0; // Leave immediately by +ve side
                                                   >> 847         snside=kPX;
                                                   >> 848       }
                                                   >> 849     }
                                                   >> 850     else if (ds1>0&&ds2<0)
587     {                                             851     {
588       G4double dist = fPlanes[i].b*p.y()+fPlan << 852       if (ss1>kCarTolerance/2)
589       if (dist >= -halfCarTolerance)           << 
590       {                                           853       {
591         if (calcNorm)                          << 854         // sn=ss1/ds1;  // Leave by +ve side
                                                   >> 855         if (ss2<-kCarTolerance/2)
592         {                                         856         {
593           *validNorm = true;                   << 857           sn=ss1/ds1;  // Leave by +ve side
594           n->set(0, fPlanes[i].b, fPlanes[i].c << 858           sn2=ss2/ds2;
                                                   >> 859           if (sn2<sn)
                                                   >> 860           {
                                                   >> 861             sn=sn2;
                                                   >> 862             snside=kMX;
                                                   >> 863           }
                                                   >> 864           else
                                                   >> 865           {
                                                   >> 866             snside=kPX;
                                                   >> 867           }
595         }                                         868         }
596         return 0;                              << 869         else
                                                   >> 870         {
                                                   >> 871           sn=0; // Leave immediately by -ve
                                                   >> 872           snside=kMX;
                                                   >> 873         }      
597       }                                           874       }
598       G4double tmp = -dist/cosa;               << 875       else
599       if (tmax > tmp) { tmax = tmp; iside = i; << 876       {
                                                   >> 877         sn=0; // Leave immediately by +ve side
                                                   >> 878         snside=kPX;
                                                   >> 879       }
                                                   >> 880     }
                                                   >> 881     else
                                                   >> 882     {
                                                   >> 883       // Must be || to both
                                                   >> 884       //
                                                   >> 885       sn=kInfinity;    // Don't leave by either side
600     }                                             886     }
601   }                                               887   }
                                                   >> 888   else if (ss1<=0&&ss2<0)
                                                   >> 889   {
                                                   >> 890     // Outside, in +ve Area
                                                   >> 891     
                                                   >> 892     if (ds1>0)
                                                   >> 893     {
                                                   >> 894       sn=0;       // Away from shape
                                                   >> 895                   // Left by +ve side
                                                   >> 896       snside=kPX;
                                                   >> 897     }
                                                   >> 898     else
                                                   >> 899     {
                                                   >> 900       if (ds2<0)
                                                   >> 901       {
                                                   >> 902         // Ignore +ve plane and use -ve plane intersect
                                                   >> 903         //
                                                   >> 904         sn=ss2/ds2; // Leave by -ve side
                                                   >> 905         snside=kMX;
                                                   >> 906       }
                                                   >> 907       else
                                                   >> 908       {
                                                   >> 909         // Must be || to both -> exit determined by other axes
                                                   >> 910         //
                                                   >> 911         sn=kInfinity; // Don't leave by either side
                                                   >> 912       }
                                                   >> 913     }
                                                   >> 914   }
                                                   >> 915   else if (ss1>0&&ss2>=0)
                                                   >> 916   {
                                                   >> 917     // Outside, in -ve Area
602                                                   918 
603   // X intersections                           << 919     if (ds2<0)
604   //                                           << 920     {
605   for ( ; i<4; ++i)                            << 921       sn=0;       // away from shape
                                                   >> 922                   // Left by -ve side
                                                   >> 923       snside=kMX;
                                                   >> 924     }
                                                   >> 925     else
                                                   >> 926     {
                                                   >> 927       if (ds1>0)
                                                   >> 928       {
                                                   >> 929         // Ignore +ve plane and use -ve plane intersect
                                                   >> 930         //
                                                   >> 931         sn=ss1/ds1; // Leave by +ve side
                                                   >> 932         snside=kPX;
                                                   >> 933       }
                                                   >> 934       else
                                                   >> 935       {
                                                   >> 936         // Must be || to both -> exit determined by other axes
                                                   >> 937         //
                                                   >> 938         sn=kInfinity; // Don't leave by either side
                                                   >> 939       }
                                                   >> 940     }
                                                   >> 941   }
                                                   >> 942 
                                                   >> 943   // Update minimum exit distance
                                                   >> 944 
                                                   >> 945   if (sn<snxt)
                                                   >> 946   {
                                                   >> 947     snxt=sn;
                                                   >> 948     side=snside;
                                                   >> 949   }
                                                   >> 950   if (snxt>0)
606   {                                               951   {
607     G4double cosa = fPlanes[i].a*v.x()+fPlanes << 952     // Calculate y intersection
608     if (cosa > 0)                              << 953 
                                                   >> 954     tanyz=(fDy2-fDy1)*0.5/fDz;
                                                   >> 955     central=0.5*(fDy1+fDy2);
                                                   >> 956 
                                                   >> 957     // +ve plane (1)
                                                   >> 958     //
                                                   >> 959     ss1=central+tanyz*p.z()-p.y(); // distance || y axis to plane
                                                   >> 960                                    // (+ve if point inside)
                                                   >> 961     ds1=v.y()-tanyz*v.z();  // component towards +ve plane
                                                   >> 962                             // (-ve if +ve -> -ve direction)
                                                   >> 963     // -ve plane (2)
                                                   >> 964     //
                                                   >> 965     ss2=-tanyz*p.z()-p.y()-central; // distance || y axis to plane
                                                   >> 966                                     // (-ve if point inside)
                                                   >> 967     ds2=tanyz*v.z()+v.y();  // component towards -ve plane
                                                   >> 968 
                                                   >> 969     if (ss1>0&&ss2<0)
609     {                                             970     {
610       G4double dist = fPlanes[i].a*p.x()+fPlan << 971       // Normal case - entirely inside region
611       if (dist >= -halfCarTolerance)           << 972 
                                                   >> 973       if (ds1<=0&&ds2<0)
                                                   >> 974       {   
                                                   >> 975         if (ss2<-kCarTolerance/2)
                                                   >> 976         {
                                                   >> 977           sn=ss2/ds2;  // Leave by -ve side
                                                   >> 978           snside=kMY;
                                                   >> 979         }
                                                   >> 980         else
                                                   >> 981         {
                                                   >> 982           sn=0; // Leave immediately by -ve side
                                                   >> 983           snside=kMY;
                                                   >> 984         }
                                                   >> 985       }
                                                   >> 986       else if (ds1>0&&ds2>=0)
                                                   >> 987       {
                                                   >> 988         if (ss1>kCarTolerance/2)
                                                   >> 989         {
                                                   >> 990           sn=ss1/ds1;  // Leave by +ve side
                                                   >> 991           snside=kPY;
                                                   >> 992         }
                                                   >> 993         else
                                                   >> 994         {
                                                   >> 995           sn=0; // Leave immediately by +ve side
                                                   >> 996           snside=kPY;
                                                   >> 997         }
                                                   >> 998       }
                                                   >> 999       else if (ds1>0&&ds2<0)
612       {                                           1000       {
613         if (calcNorm)                          << 1001         if (ss1>kCarTolerance/2)
614         {                                         1002         {
615            *validNorm = true;                  << 1003           // sn=ss1/ds1;  // Leave by +ve side
616            n->set(fPlanes[i].a, fPlanes[i].b,  << 1004           if (ss2<-kCarTolerance/2)
                                                   >> 1005           {
                                                   >> 1006             sn=ss1/ds1;  // Leave by +ve side
                                                   >> 1007             sn2=ss2/ds2;
                                                   >> 1008             if (sn2<sn)
                                                   >> 1009             {
                                                   >> 1010               sn=sn2;
                                                   >> 1011               snside=kMY;
                                                   >> 1012             }
                                                   >> 1013             else
                                                   >> 1014             {
                                                   >> 1015               snside=kPY;
                                                   >> 1016             }
                                                   >> 1017           }
                                                   >> 1018           else
                                                   >> 1019           {
                                                   >> 1020             sn=0; // Leave immediately by -ve
                                                   >> 1021             snside=kMY;
                                                   >> 1022           }
617         }                                         1023         }
618         return 0;                              << 1024         else
                                                   >> 1025         {
                                                   >> 1026           sn=0; // Leave immediately by +ve side
                                                   >> 1027           snside=kPY;
                                                   >> 1028         }
                                                   >> 1029       }
                                                   >> 1030       else
                                                   >> 1031       {
                                                   >> 1032         // Must be || to both
                                                   >> 1033         //
                                                   >> 1034         sn=kInfinity;    // Don't leave by either side
619       }                                           1035       }
620       G4double tmp = -dist/cosa;               << 1036     }
621       if (tmax > tmp) { tmax = tmp; iside = i; << 1037     else if (ss1<=0&&ss2<0)
                                                   >> 1038     {
                                                   >> 1039       // Outside, in +ve Area
                                                   >> 1040 
                                                   >> 1041       if (ds1>0)
                                                   >> 1042       {
                                                   >> 1043         sn=0;       // Away from shape
                                                   >> 1044                     // Left by +ve side
                                                   >> 1045         snside=kPY;
                                                   >> 1046       }
                                                   >> 1047       else
                                                   >> 1048       {
                                                   >> 1049         if (ds2<0)
                                                   >> 1050         {
                                                   >> 1051           // Ignore +ve plane and use -ve plane intersect
                                                   >> 1052           //
                                                   >> 1053           sn=ss2/ds2; // Leave by -ve side
                                                   >> 1054           snside=kMY;
                                                   >> 1055         }
                                                   >> 1056         else
                                                   >> 1057         {
                                                   >> 1058           // Must be || to both -> exit determined by other axes
                                                   >> 1059           //
                                                   >> 1060           sn=kInfinity; // Don't leave by either side
                                                   >> 1061         }
                                                   >> 1062       }
                                                   >> 1063     }
                                                   >> 1064     else if (ss1>0&&ss2>=0)
                                                   >> 1065     {
                                                   >> 1066       // Outside, in -ve Area
                                                   >> 1067       if (ds2<0)
                                                   >> 1068       {
                                                   >> 1069         sn=0;       // away from shape
                                                   >> 1070                     // Left by -ve side
                                                   >> 1071         snside=kMY;
                                                   >> 1072       }
                                                   >> 1073       else
                                                   >> 1074       {
                                                   >> 1075         if (ds1>0)
                                                   >> 1076         {
                                                   >> 1077           // Ignore +ve plane and use -ve plane intersect
                                                   >> 1078           //
                                                   >> 1079           sn=ss1/ds1; // Leave by +ve side
                                                   >> 1080           snside=kPY;
                                                   >> 1081         }
                                                   >> 1082         else
                                                   >> 1083         {
                                                   >> 1084           // Must be || to both -> exit determined by other axes
                                                   >> 1085           //
                                                   >> 1086           sn=kInfinity; // Don't leave by either side
                                                   >> 1087         }
                                                   >> 1088       }
                                                   >> 1089     }
                                                   >> 1090 
                                                   >> 1091     // Update minimum exit distance
                                                   >> 1092 
                                                   >> 1093     if (sn<snxt)
                                                   >> 1094     {
                                                   >> 1095       snxt=sn;
                                                   >> 1096       side=snside;
622     }                                             1097     }
623   }                                               1098   }
624                                                   1099 
625   // Set normal, if required, and return dista << 
626   //                                           << 
627   if (calcNorm)                                   1100   if (calcNorm)
628   {                                               1101   {
629     *validNorm = true;                         << 1102     switch (side)
630     if (iside < 0)                             << 1103     {
631       n->set(0, 0, iside + 3); // (-4+3)=-1, ( << 1104       case kPX:
632     else                                       << 1105         cosxz=1.0/sqrt(1.0+tanxz*tanxz);
633       n->set(fPlanes[iside].a, fPlanes[iside]. << 1106         *n=G4ThreeVector(cosxz,0,-tanxz*cosxz);
                                                   >> 1107         break;
                                                   >> 1108       case kMX:
                                                   >> 1109         cosxz=-1.0/sqrt(1.0+tanxz*tanxz);
                                                   >> 1110         *n=G4ThreeVector(cosxz,0,tanxz*cosxz);
                                                   >> 1111         break;
                                                   >> 1112       case kPY:
                                                   >> 1113         cosyz=1.0/sqrt(1.0+tanyz*tanyz);
                                                   >> 1114         *n=G4ThreeVector(0,cosyz,-tanyz*cosyz);
                                                   >> 1115         break;
                                                   >> 1116       case kMY:
                                                   >> 1117         cosyz=-1.0/sqrt(1.0+tanyz*tanyz);
                                                   >> 1118         *n=G4ThreeVector(0,cosyz,tanyz*cosyz);
                                                   >> 1119         break;
                                                   >> 1120       case kPZ:
                                                   >> 1121         *n=G4ThreeVector(0,0,1);
                                                   >> 1122         break;
                                                   >> 1123       case kMZ:
                                                   >> 1124         *n=G4ThreeVector(0,0,-1);
                                                   >> 1125         break;
                                                   >> 1126       default:
                                                   >> 1127         DumpInfo();
                                                   >> 1128         G4Exception("G4Trd::DistanceToOut() - Invalid enum");
                                                   >> 1129         break;
                                                   >> 1130     }
634   }                                               1131   }
635   return tmax;                                 << 1132   return snxt; 
636 }                                                 1133 }
637                                                   1134 
638 ////////////////////////////////////////////// << 1135 ///////////////////////////////////////////////////////////////////////////
639 //                                                1136 //
640 // Calculate exact shortest distance to any bo    1137 // Calculate exact shortest distance to any boundary from inside
641 // - returns 0 if point is outside             << 1138 // - Returns 0 is point outside
642                                                   1139 
643 G4double G4Trd::DistanceToOut( const G4ThreeVe    1140 G4double G4Trd::DistanceToOut( const G4ThreeVector& p ) const
644 {                                                 1141 {
                                                   >> 1142   G4double safe;
                                                   >> 1143   G4double tanxz,xdist,saf1;
                                                   >> 1144   G4double tanyz,ydist,saf2;
                                                   >> 1145   G4double zbase;
                                                   >> 1146 
645 #ifdef G4CSGDEBUG                                 1147 #ifdef G4CSGDEBUG
646   if( Inside(p) == kOutside )                     1148   if( Inside(p) == kOutside )
647   {                                               1149   {
648     std::ostringstream message;                << 1150      G4cout.precision(16) ;
649     G4long oldprc = message.precision(16);     << 1151      G4cout << G4endl ;
650     message << "Point p is outside (!?) of sol << 1152      DumpInfo();
651     message << "Position:\n";                  << 1153      G4cout << "Position:"  << G4endl << G4endl ;
652     message << "   p.x() = " << p.x()/mm << "  << 1154      G4cout << "p.x() = "   << p.x()/mm << " mm" << G4endl ;
653     message << "   p.y() = " << p.y()/mm << "  << 1155      G4cout << "p.y() = "   << p.y()/mm << " mm" << G4endl ;
654     message << "   p.z() = " << p.z()/mm << "  << 1156      G4cout << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl ;
655     G4cout.precision(oldprc);                  << 1157      G4cout << "G4Trd::DistanceToOut(p) - point p is outside ?!" << G4endl ;
656     G4Exception("G4Trd::DistanceToOut(p)", "Ge << 1158      G4cerr << "G4Trd::DistanceToOut(p) - point p is outside ?!" << G4endl ;
657                 JustWarning, message );        << 
658     DumpInfo();                                << 
659   }                                               1159   }
660 #endif                                            1160 #endif
661   G4double dx = fPlanes[3].a*std::abs(p.x())+f << 
662   G4double dy = fPlanes[1].b*std::abs(p.y())+f << 
663   G4double dxy = std::max(dx,dy);              << 
664                                                   1161 
665   G4double dz = std::abs(p.z())-fDz;           << 1162   safe=fDz-fabs(p.z());  // z perpendicular Dist
666   G4double dist = std::max(dz,dxy);            << 
667                                                   1163 
668   return (dist < 0) ? -dist : 0.;              << 1164   zbase=fDz+p.z();
669 }                                              << 
670                                                   1165 
671 ////////////////////////////////////////////// << 1166   // xdist = distance perpendicular to z axis to closest x plane from p
672 //                                             << 1167   //       = (x half width of shape at p.z) - fabs(p.x)
673 // GetEntityType                               << 1168   //
                                                   >> 1169   tanxz=(fDx2-fDx1)*0.5/fDz;
                                                   >> 1170   xdist=fDx1+tanxz*zbase-fabs(p.x());
                                                   >> 1171   saf1=xdist/sqrt(1.0+tanxz*tanxz); // x*cos(ang_xz) =
                                                   >> 1172                                     // shortest (perpendicular)
                                                   >> 1173                                     // distance to plane
                                                   >> 1174   tanyz=(fDy2-fDy1)*0.5/fDz;
                                                   >> 1175   ydist=fDy1+tanyz*zbase-fabs(p.y());
                                                   >> 1176   saf2=ydist/sqrt(1.0+tanyz*tanyz);
674                                                   1177 
675 G4GeometryType G4Trd::GetEntityType() const    << 1178   // Return minimum x/y/z distance
676 {                                              << 1179   //
677   return {"G4Trd"};                            << 1180   if (safe>saf1) safe=saf1;
                                                   >> 1181   if (safe>saf2) safe=saf2;
                                                   >> 1182 
                                                   >> 1183   if (safe<0) safe=0;
                                                   >> 1184   return safe;     
678 }                                                 1185 }
679                                                   1186 
680 ////////////////////////////////////////////// << 1187 ////////////////////////////////////////////////////////////////////////////
681 //                                                1188 //
682 // IsFaceted                                   << 1189 // Create a List containing the transformed vertices
683                                                << 1190 // Ordering [0-3] -fDz cross section
684 G4bool G4Trd::IsFaceted() const                << 1191 //          [4-7] +fDz cross section such that [0] is below [4],
685 {                                              << 1192 //                                             [1] below [5] etc.
686   return true;                                 << 1193 // Note:
                                                   >> 1194 //  Caller has deletion resposibility
                                                   >> 1195 
                                                   >> 1196 G4ThreeVectorList*
                                                   >> 1197 G4Trd::CreateRotatedVertices( const G4AffineTransform& pTransform ) const
                                                   >> 1198 {
                                                   >> 1199   G4ThreeVectorList *vertices;
                                                   >> 1200   vertices=new G4ThreeVectorList();
                                                   >> 1201   vertices->reserve(8);
                                                   >> 1202   if (vertices)
                                                   >> 1203   {
                                                   >> 1204     G4ThreeVector vertex0(-fDx1,-fDy1,-fDz);
                                                   >> 1205     G4ThreeVector vertex1(fDx1,-fDy1,-fDz);
                                                   >> 1206     G4ThreeVector vertex2(fDx1,fDy1,-fDz);
                                                   >> 1207     G4ThreeVector vertex3(-fDx1,fDy1,-fDz);
                                                   >> 1208     G4ThreeVector vertex4(-fDx2,-fDy2,fDz);
                                                   >> 1209     G4ThreeVector vertex5(fDx2,-fDy2,fDz);
                                                   >> 1210     G4ThreeVector vertex6(fDx2,fDy2,fDz);
                                                   >> 1211     G4ThreeVector vertex7(-fDx2,fDy2,fDz);
                                                   >> 1212 
                                                   >> 1213     vertices->push_back(pTransform.TransformPoint(vertex0));
                                                   >> 1214     vertices->push_back(pTransform.TransformPoint(vertex1));
                                                   >> 1215     vertices->push_back(pTransform.TransformPoint(vertex2));
                                                   >> 1216     vertices->push_back(pTransform.TransformPoint(vertex3));
                                                   >> 1217     vertices->push_back(pTransform.TransformPoint(vertex4));
                                                   >> 1218     vertices->push_back(pTransform.TransformPoint(vertex5));
                                                   >> 1219     vertices->push_back(pTransform.TransformPoint(vertex6));
                                                   >> 1220     vertices->push_back(pTransform.TransformPoint(vertex7));
                                                   >> 1221   }
                                                   >> 1222   else
                                                   >> 1223   {
                                                   >> 1224     DumpInfo();
                                                   >> 1225     G4Exception("G4Trd::CreateRotatedVertices() - Out of memory");
                                                   >> 1226   }
                                                   >> 1227   return vertices;
687 }                                                 1228 }
688                                                   1229 
689 //////////////////////////////////////////////    1230 //////////////////////////////////////////////////////////////////////////
690 //                                                1231 //
691 // Make a clone of the object                  << 1232 // GetEntityType
692 //                                             << 1233 
693 G4VSolid* G4Trd::Clone() const                 << 1234 G4GeometryType G4Trd::GetEntityType() const
694 {                                                 1235 {
695   return new G4Trd(*this);                     << 1236   return G4String("G4Trd");
696 }                                                 1237 }
697                                                   1238 
698 //////////////////////////////////////////////    1239 //////////////////////////////////////////////////////////////////////////
699 //                                                1240 //
700 // Stream object contents to an output stream     1241 // Stream object contents to an output stream
701                                                   1242 
702 std::ostream& G4Trd::StreamInfo( std::ostream&    1243 std::ostream& G4Trd::StreamInfo( std::ostream& os ) const
703 {                                                 1244 {
704   G4long oldprc = os.precision(16);            << 
705   os << "-------------------------------------    1245   os << "-----------------------------------------------------------\n"
706      << "    *** Dump for solid - " << GetName    1246      << "    *** Dump for solid - " << GetName() << " ***\n"
707      << "    =================================    1247      << "    ===================================================\n"
708      << " Solid type: G4Trd\n"                    1248      << " Solid type: G4Trd\n"
709      << " Parameters: \n"                         1249      << " Parameters: \n"
710      << "    half length X, surface -dZ: " <<     1250      << "    half length X, surface -dZ: " << fDx1/mm << " mm \n"
711      << "    half length X, surface +dZ: " <<     1251      << "    half length X, surface +dZ: " << fDx2/mm << " mm \n"
712      << "    half length Y, surface -dZ: " <<     1252      << "    half length Y, surface -dZ: " << fDy1/mm << " mm \n"
713      << "    half length Y, surface +dZ: " <<     1253      << "    half length Y, surface +dZ: " << fDy2/mm << " mm \n"
714      << "    half length Z             : " <<  << 1254      << "    half length Z             : " << fDz/mm << " mm \n"
715      << "-------------------------------------    1255      << "-----------------------------------------------------------\n";
716   os.precision(oldprc);                        << 
717                                                   1256 
718   return os;                                      1257   return os;
719 }                                                 1258 }
720                                                   1259 
721 ////////////////////////////////////////////// << 1260 ///////////////////////////////////////////////////////////////////////
722 //                                             << 
723 // Return a point randomly and uniformly selec << 
724                                                << 
725 G4ThreeVector G4Trd::GetPointOnSurface() const << 
726 {                                              << 
727   // Set areas                                 << 
728   //                                           << 
729   G4double sxz = (fDx1 + fDx2)*fHx;            << 
730   G4double syz = (fDy1 + fDy2)*fHy;            << 
731   G4double ssurf[6] = { 4.*fDx1*fDy1, sxz, sxz << 
732   ssurf[1] += ssurf[0];                        << 
733   ssurf[2] += ssurf[1];                        << 
734   ssurf[3] += ssurf[2];                        << 
735   ssurf[4] += ssurf[3];                        << 
736   ssurf[5] += ssurf[4];                        << 
737                                                << 
738   // Select face                               << 
739   //                                           << 
740   G4double select = ssurf[5]*G4QuickRand();    << 
741   G4int k = 5;                                 << 
742   k -= (G4int)(select <= ssurf[4]);            << 
743   k -= (G4int)(select <= ssurf[3]);            << 
744   k -= (G4int)(select <= ssurf[2]);            << 
745   k -= (G4int)(select <= ssurf[1]);            << 
746   k -= (G4int)(select <= ssurf[0]);            << 
747                                                << 
748   // Generate point on selected surface        << 
749   //                                           << 
750   G4double u = G4QuickRand();                  << 
751   G4double v = G4QuickRand();                  << 
752   switch(k)                                    << 
753   {                                            << 
754     case 0: // base at -Z                      << 
755     {                                          << 
756       return { (2.*u - 1.)*fDx1, (2.*v - 1.)*f << 
757     }                                          << 
758     case 1: // X face at -Y                    << 
759     {                                          << 
760       if (u + v > 1.) { u = 1. - u; v = 1. - v << 
761       G4ThreeVector p0(-fDx1,-fDy1,-fDz);      << 
762       G4ThreeVector p1( fDx2,-fDy2, fDz);      << 
763       return (select <= ssurf[0] + fDx1*fHx) ? << 
764         (1. - u - v)*p0 + u*p1 + v*G4ThreeVect << 
765         (1. - u - v)*p0 + u*p1 + v*G4ThreeVect << 
766     }                                          << 
767     case 2: // X face at +Y                    << 
768     {                                          << 
769       if (u + v > 1.) { u = 1. - u; v = 1. - v << 
770       G4ThreeVector p0( fDx1, fDy1,-fDz);      << 
771       G4ThreeVector p1(-fDx2, fDy2, fDz);      << 
772       return (select <= ssurf[1] + fDx1*fHx) ? << 
773         (1. - u - v)*p0 + u*p1 + v*G4ThreeVect << 
774         (1. - u - v)*p0 + u*p1 + v*G4ThreeVect << 
775     }                                          << 
776     case 3: // Y face at -X                    << 
777     {                                          << 
778       if (u + v > 1.) { u = 1. - u; v = 1. - v << 
779       G4ThreeVector p0(-fDx1, fDy1,-fDz);      << 
780       G4ThreeVector p1(-fDx2,-fDy2, fDz);      << 
781       return (select <= ssurf[2] + fDy1*fHy) ? << 
782         (1. - u - v)*p0 + u*p1 + v*G4ThreeVect << 
783         (1. - u - v)*p0 + u*p1 + v*G4ThreeVect << 
784     }                                          << 
785     case 4: // Y face at +X                    << 
786     {                                          << 
787       if (u + v > 1.) { u = 1. - u; v = 1. - v << 
788       G4ThreeVector p0( fDx1,-fDy1,-fDz);      << 
789       G4ThreeVector p1( fDx2, fDy2, fDz);      << 
790       return (select <= ssurf[3] + fDy1*fHy) ? << 
791         (1. - u - v)*p0 + u*p1 + v*G4ThreeVect << 
792         (1. - u - v)*p0 + u*p1 + v*G4ThreeVect << 
793     }                                          << 
794     case 5: // base at +Z                      << 
795     {                                          << 
796       return { (2.*u - 1.)*fDx2, (2.*v - 1.)*f << 
797     }                                          << 
798   }                                            << 
799   return {0., 0., 0.};                         << 
800 }                                              << 
801                                                << 
802 ////////////////////////////////////////////// << 
803 //                                                1261 //
804 // Methods for visualisation                      1262 // Methods for visualisation
805                                                   1263 
806 void G4Trd::DescribeYourselfTo ( G4VGraphicsSc    1264 void G4Trd::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
807 {                                                 1265 {
808   scene.AddSolid (*this);                      << 1266   scene.AddThis (*this);
809 }                                                 1267 }
810                                                   1268 
811 G4Polyhedron* G4Trd::CreatePolyhedron () const    1269 G4Polyhedron* G4Trd::CreatePolyhedron () const
812 {                                                 1270 {
813   return new G4PolyhedronTrd2 (fDx1, fDx2, fDy    1271   return new G4PolyhedronTrd2 (fDx1, fDx2, fDy1, fDy2, fDz);
814 }                                                 1272 }
815                                                   1273 
816 #endif                                         << 1274 G4NURBS* G4Trd::CreateNURBS () const
                                                   >> 1275 {
                                                   >> 1276   //  return new G4NURBSbox (fDx, fDy, fDz);
                                                   >> 1277   return 0;
                                                   >> 1278 }
                                                   >> 1279 
                                                   >> 1280 //
                                                   >> 1281 //
                                                   >> 1282 ///////////////////////////////////////////////////////////////////////////
817                                                   1283