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.2.p1)


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