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 9.4.p3)


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