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 10.4.p2)


  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 108482 2018-02-15 14:34:23Z gcosmo $
                                                   >>  28 //
                                                   >>  29 //
 26 // Implementation for G4Trd class                  30 // Implementation for G4Trd class
 27 //                                                 31 //
 28 // 12.01.95 P.Kent: First version              <<  32 // History:
 29 // 28.04.05 V.Grichine: new SurfaceNormal acco <<  33 //
 30 // 25.05.17 E.Tcherniaev: complete revision, s     34 // 25.05.17 E.Tcherniaev: complete revision, speed-up
 31 // ------------------------------------------- <<  35 // 23.09.16 E.Tcherniaev: use G4BoundingEnvelope for CalculateExtent(),
                                                   >>  36 //                      removed CreateRotatedVertices()
                                                   >>  37 // 28.04.05 V.Grichine: new SurfaceNormal according to J. Apostolakis proposal 
                                                   >>  38 // 26.04.05, V.Grichine, new SurfaceNoramal is default
                                                   >>  39 // 07.12.04, V.Grichine, SurfaceNoramal with edges/vertices.
                                                   >>  40 // 07.05.00, V.Grichine, in d = DistanceToIn(p,v), if d<0.5*kCarTolerance, d=0
                                                   >>  41 //    ~1996, V.Grichine, 1st implementation based on old code of P.Kent
                                                   >>  42 //
                                                   >>  43 //////////////////////////////////////////////////////////////////////////////
 32                                                    44 
 33 #include "G4Trd.hh"                                45 #include "G4Trd.hh"
 34                                                    46 
 35 #if !defined(G4GEOM_USE_UTRD)                      47 #if !defined(G4GEOM_USE_UTRD)
 36                                                    48 
 37 #include "G4GeomTools.hh"                          49 #include "G4GeomTools.hh"
 38                                                    50 
 39 #include "G4VoxelLimits.hh"                        51 #include "G4VoxelLimits.hh"
 40 #include "G4AffineTransform.hh"                    52 #include "G4AffineTransform.hh"
 41 #include "G4BoundingEnvelope.hh"                   53 #include "G4BoundingEnvelope.hh"
 42 #include "G4QuickRand.hh"                      <<  54 #include "Randomize.hh"
 43                                                    55 
 44 #include "G4VPVParameterisation.hh"                56 #include "G4VPVParameterisation.hh"
 45                                                    57 
 46 #include "G4VGraphicsScene.hh"                     58 #include "G4VGraphicsScene.hh"
 47                                                    59 
 48 using namespace CLHEP;                             60 using namespace CLHEP;
 49                                                    61 
 50 //////////////////////////////////////////////     62 //////////////////////////////////////////////////////////////////////////
 51 //                                                 63 //
 52 // Constructor - set & check half widths           64 // Constructor - set & check half widths
 53                                                    65 
 54 G4Trd::G4Trd(const G4String& pName,                66 G4Trd::G4Trd(const G4String& pName,
 55                    G4double pdx1, G4double pdx     67                    G4double pdx1, G4double pdx2,
 56                    G4double pdy1, G4double pdy     68                    G4double pdy1, G4double pdy2,
 57                    G4double pdz)                   69                    G4double pdz)
 58   : G4CSGSolid(pName), halfCarTolerance(0.5*kC     70   : G4CSGSolid(pName), halfCarTolerance(0.5*kCarTolerance),
 59     fDx1(pdx1), fDx2(pdx2), fDy1(pdy1), fDy2(p     71     fDx1(pdx1), fDx2(pdx2), fDy1(pdy1), fDy2(pdy2), fDz(pdz)
 60 {                                                  72 {
 61   CheckParameters();                               73   CheckParameters();
 62   MakePlanes();                                    74   MakePlanes();
 63 }                                                  75 }
 64                                                    76 
 65 //////////////////////////////////////////////     77 //////////////////////////////////////////////////////////////////////////
 66 //                                                 78 //
 67 // Fake default constructor - sets only member     79 // Fake default constructor - sets only member data and allocates memory
 68 //                            for usage restri     80 //                            for usage restricted to object persistency
 69 //                                                 81 //
 70 G4Trd::G4Trd( __void__& a )                        82 G4Trd::G4Trd( __void__& a )
 71   : G4CSGSolid(a), halfCarTolerance(0.5*kCarTo     83   : G4CSGSolid(a), halfCarTolerance(0.5*kCarTolerance),
 72     fDx1(1.), fDx2(1.), fDy1(1.), fDy2(1.), fD     84     fDx1(1.), fDx2(1.), fDy1(1.), fDy2(1.), fDz(1.)
 73 {                                                  85 {
 74   MakePlanes();                                    86   MakePlanes();
 75 }                                                  87 }
 76                                                    88 
 77 //////////////////////////////////////////////     89 //////////////////////////////////////////////////////////////////////////
 78 //                                                 90 //
 79 // Destructor                                      91 // Destructor
 80                                                    92 
 81 G4Trd::~G4Trd() = default;                     <<  93 G4Trd::~G4Trd()
                                                   >>  94 {
                                                   >>  95 }
 82                                                    96 
 83 //////////////////////////////////////////////     97 //////////////////////////////////////////////////////////////////////////
 84 //                                                 98 //
 85 // Copy constructor                                99 // Copy constructor
 86                                                   100 
 87 G4Trd::G4Trd(const G4Trd& rhs)                    101 G4Trd::G4Trd(const G4Trd& rhs)
 88   : G4CSGSolid(rhs), halfCarTolerance(rhs.half    102   : G4CSGSolid(rhs), halfCarTolerance(rhs.halfCarTolerance),
 89     fDx1(rhs.fDx1), fDx2(rhs.fDx2),               103     fDx1(rhs.fDx1), fDx2(rhs.fDx2),
 90     fDy1(rhs.fDy1), fDy2(rhs.fDy2), fDz(rhs.fD << 104     fDy1(rhs.fDy1), fDy2(rhs.fDy2), fDz(rhs.fDz)
 91     fHx(rhs.fHx), fHy(rhs.fHy)                 << 
 92 {                                                 105 {
 93   for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs    106   for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs.fPlanes[i]; }
 94 }                                                 107 }
 95                                                   108 
 96 //////////////////////////////////////////////    109 //////////////////////////////////////////////////////////////////////////
 97 //                                                110 //
 98 // Assignment operator                            111 // Assignment operator
 99                                                   112 
100 G4Trd& G4Trd::operator = (const G4Trd& rhs)       113 G4Trd& G4Trd::operator = (const G4Trd& rhs)
101 {                                                 114 {
102    // Check assignment to self                    115    // Check assignment to self
103    //                                             116    //
104    if (this == &rhs)  { return *this; }           117    if (this == &rhs)  { return *this; }
105                                                   118 
106    // Copy base class data                        119    // Copy base class data
107    //                                             120    //
108    G4CSGSolid::operator=(rhs);                    121    G4CSGSolid::operator=(rhs);
109                                                   122 
110    // Copy data                                   123    // Copy data
111    //                                             124    //
112    halfCarTolerance = rhs.halfCarTolerance;       125    halfCarTolerance = rhs.halfCarTolerance;
113    fDx1 = rhs.fDx1; fDx2 = rhs.fDx2;              126    fDx1 = rhs.fDx1; fDx2 = rhs.fDx2;
114    fDy1 = rhs.fDy1; fDy2 = rhs.fDy2;              127    fDy1 = rhs.fDy1; fDy2 = rhs.fDy2;
115    fDz = rhs.fDz;                                 128    fDz = rhs.fDz;
116    fHx = rhs.fHx; fHy = rhs.fHy;               << 
117    for (G4int i=0; i<4; ++i) { fPlanes[i] = rh    129    for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs.fPlanes[i]; }
118                                                   130 
119    return *this;                                  131    return *this;
120 }                                                 132 }
121                                                   133 
122 //////////////////////////////////////////////    134 //////////////////////////////////////////////////////////////////////////
123 //                                                135 //
124 // Set all parameters, as for constructor - se    136 // Set all parameters, as for constructor - set and check half-widths
125                                                   137 
126 void G4Trd::SetAllParameters(G4double pdx1, G4    138 void G4Trd::SetAllParameters(G4double pdx1, G4double pdx2,
127                              G4double pdy1, G4    139                              G4double pdy1, G4double pdy2, G4double pdz)
128 {                                                 140 {
129   // Reset data of the base class                 141   // Reset data of the base class
130   fCubicVolume = 0.;                           << 142   fCubicVolume = 0;
131   fSurfaceArea = 0.;                           << 143   fSurfaceArea = 0;
132   fRebuildPolyhedron = true;                      144   fRebuildPolyhedron = true;
133                                                   145 
134   // Set parameters                               146   // Set parameters
135   fDx1 = pdx1; fDx2 = pdx2;                       147   fDx1 = pdx1; fDx2 = pdx2;
136   fDy1 = pdy1; fDy2 = pdy2;                       148   fDy1 = pdy1; fDy2 = pdy2;
137   fDz  = pdz;                                     149   fDz  = pdz;
138                                                   150 
139   CheckParameters();                              151   CheckParameters();
140   MakePlanes();                                   152   MakePlanes();
141 }                                                 153 }
142                                                   154 
143 //////////////////////////////////////////////    155 //////////////////////////////////////////////////////////////////////////
144 //                                                156 //
145 // Check dimensions                               157 // Check dimensions
146                                                   158 
147 void G4Trd::CheckParameters()                     159 void G4Trd::CheckParameters()
148 {                                                 160 {
149   G4double dmin = 2*kCarTolerance;                161   G4double dmin = 2*kCarTolerance;
150   if ((fDx1 < 0 || fDx2 < 0 || fDy1 < 0 || fDy    162   if ((fDx1 < 0 || fDx2 < 0 || fDy1 < 0 || fDy2 < 0 || fDz < dmin) ||
151       (fDx1 < dmin && fDx2 < dmin) ||             163       (fDx1 < dmin && fDx2 < dmin) ||
152       (fDy1 < dmin && fDy2 < dmin))               164       (fDy1 < dmin && fDy2 < dmin))
153   {                                               165   {
154     std::ostringstream message;                   166     std::ostringstream message;
155     message << "Invalid (too small or negative    167     message << "Invalid (too small or negative) dimensions for Solid: "
156             << GetName()                          168             << GetName()
157             << "\n  X - " << fDx1 << ", " << f    169             << "\n  X - " << fDx1 << ", " << fDx2
158             << "\n  Y - " << fDy1 << ", " << f    170             << "\n  Y - " << fDy1 << ", " << fDy2
159             << "\n  Z - " << fDz;                 171             << "\n  Z - " << fDz;
160     G4Exception("G4Trd::CheckParameters()", "G    172     G4Exception("G4Trd::CheckParameters()", "GeomSolids0002",
161                 FatalException, message);         173                 FatalException, message);
162   }                                               174   }
163 }                                                 175 }
164                                                   176 
165 //////////////////////////////////////////////    177 //////////////////////////////////////////////////////////////////////////
166 //                                                178 //
167 // Set side planes                                179 // Set side planes
168                                                   180 
169 void G4Trd::MakePlanes()                          181 void G4Trd::MakePlanes()
170 {                                                 182 {
171   G4double dx = fDx1 - fDx2;                      183   G4double dx = fDx1 - fDx2;
172   G4double dy = fDy1 - fDy2;                      184   G4double dy = fDy1 - fDy2;
173   G4double dz = 2*fDz;                            185   G4double dz = 2*fDz;
174   fHx = std::sqrt(dy*dy + dz*dz);              << 186   G4double magx = std::sqrt(dx*dx + dz*dz);
175   fHy = std::sqrt(dx*dx + dz*dz);              << 187   G4double magy = std::sqrt(dy*dy + dz*dz);
176                                                   188 
177   // Set X planes at -Y & +Y                   << 189   // Set -Y & +Y planes
178   //                                              190   //
179   fPlanes[0].a =  0.;                             191   fPlanes[0].a =  0.;
180   fPlanes[0].b = -dz/fHx;                      << 192   fPlanes[0].b = -dz/magy;
181   fPlanes[0].c =  dy/fHx;                      << 193   fPlanes[0].c =  dy/magy;
182   fPlanes[0].d = fPlanes[0].b*fDy1 + fPlanes[0    194   fPlanes[0].d = fPlanes[0].b*fDy1 + fPlanes[0].c*fDz;
183                                                   195 
184   fPlanes[1].a =  fPlanes[0].a;                   196   fPlanes[1].a =  fPlanes[0].a;
185   fPlanes[1].b = -fPlanes[0].b;                   197   fPlanes[1].b = -fPlanes[0].b;
186   fPlanes[1].c =  fPlanes[0].c;                   198   fPlanes[1].c =  fPlanes[0].c;
187   fPlanes[1].d =  fPlanes[0].d;                   199   fPlanes[1].d =  fPlanes[0].d;
188                                                   200 
189   // Set Y planes at -X & +X                   << 201   // Set -X & +X planes
190   //                                              202   //
191   fPlanes[2].a = -dz/fHy;                      << 203   fPlanes[2].a = -dz/magx;
192   fPlanes[2].b =  0.;                             204   fPlanes[2].b =  0.;
193   fPlanes[2].c =  dx/fHy;                      << 205   fPlanes[2].c =  dx/magx;
194   fPlanes[2].d = fPlanes[2].a*fDx1 + fPlanes[2    206   fPlanes[2].d = fPlanes[2].a*fDx1 + fPlanes[2].c*fDz;
195                                                   207 
196   fPlanes[3].a = -fPlanes[2].a;                   208   fPlanes[3].a = -fPlanes[2].a;
197   fPlanes[3].b =  fPlanes[2].b;                   209   fPlanes[3].b =  fPlanes[2].b;
198   fPlanes[3].c =  fPlanes[2].c;                   210   fPlanes[3].c =  fPlanes[2].c;
199   fPlanes[3].d =  fPlanes[2].d;                   211   fPlanes[3].d =  fPlanes[2].d;
200 }                                                 212 }
201                                                   213 
202 //////////////////////////////////////////////    214 //////////////////////////////////////////////////////////////////////////
203 //                                                215 //
204 // Get volume                                     216 // Get volume
205                                                   217 
206 G4double G4Trd::GetCubicVolume()                  218 G4double G4Trd::GetCubicVolume()
207 {                                                 219 {
208   if (fCubicVolume == 0.)                      << 220   if (fCubicVolume == 0)
209   {                                               221   {
210     fCubicVolume = 2*fDz*( (fDx1+fDx2)*(fDy1+f    222     fCubicVolume = 2*fDz*( (fDx1+fDx2)*(fDy1+fDy2) +
211                            (fDx2-fDx1)*(fDy2-f    223                            (fDx2-fDx1)*(fDy2-fDy1)/3 );
212   }                                               224   }
213   return fCubicVolume;                            225   return fCubicVolume;
214 }                                                 226 }
215                                                   227 
216 //////////////////////////////////////////////    228 //////////////////////////////////////////////////////////////////////////
217 //                                                229 //
218 // Get surface area                               230 // Get surface area
219                                                   231 
220 G4double G4Trd::GetSurfaceArea()                  232 G4double G4Trd::GetSurfaceArea()
221 {                                                 233 {
222   if (fSurfaceArea == 0.)                      << 234   if (fSurfaceArea == 0)
223   {                                               235   {
224     fSurfaceArea =                                236     fSurfaceArea =
225       4*(fDx1*fDy1 + fDx2*fDy2) + 2*(fDx1+fDx2 << 237       4*(fDx1*fDy1+fDx2*fDy2) +
                                                   >> 238       2*(fDy1+fDy2)*std::hypot(fDx1-fDx2,2*fDz) +
                                                   >> 239       2*(fDx1+fDx2)*std::hypot(fDy1-fDy2,2*fDz);
226   }                                               240   }
227   return fSurfaceArea;                            241   return fSurfaceArea;
228 }                                                 242 }
229                                                   243 
230 //////////////////////////////////////////////    244 //////////////////////////////////////////////////////////////////////////
231 //                                                245 //
232 // Dispatch to parameterisation for replicatio    246 // Dispatch to parameterisation for replication mechanism dimension
233 // computation & modification                     247 // computation & modification
234                                                   248 
235 void G4Trd::ComputeDimensions(       G4VPVPara    249 void G4Trd::ComputeDimensions(       G4VPVParameterisation* p,
236                                const G4int n,     250                                const G4int n,
237                                const G4VPhysic    251                                const G4VPhysicalVolume* pRep )
238 {                                                 252 {
239   p->ComputeDimensions(*this,n,pRep);             253   p->ComputeDimensions(*this,n,pRep);
240 }                                                 254 }
241                                                   255 
242 //////////////////////////////////////////////    256 //////////////////////////////////////////////////////////////////////////
243 //                                                257 //
244 // Get bounding box                               258 // Get bounding box
245                                                   259 
246 void G4Trd::BoundingLimits(G4ThreeVector& pMin    260 void G4Trd::BoundingLimits(G4ThreeVector& pMin, G4ThreeVector& pMax) const
247 {                                                 261 {
248   G4double dx1 = GetXHalfLength1();               262   G4double dx1 = GetXHalfLength1();
249   G4double dx2 = GetXHalfLength2();               263   G4double dx2 = GetXHalfLength2();
250   G4double dy1 = GetYHalfLength1();               264   G4double dy1 = GetYHalfLength1();
251   G4double dy2 = GetYHalfLength2();               265   G4double dy2 = GetYHalfLength2();
252   G4double dz  = GetZHalfLength();                266   G4double dz  = GetZHalfLength();
253                                                   267 
254   G4double xmax = std::max(dx1,dx2);              268   G4double xmax = std::max(dx1,dx2);
255   G4double ymax = std::max(dy1,dy2);              269   G4double ymax = std::max(dy1,dy2);
256   pMin.set(-xmax,-ymax,-dz);                      270   pMin.set(-xmax,-ymax,-dz);
257   pMax.set( xmax, ymax, dz);                      271   pMax.set( xmax, ymax, dz);
258                                                   272 
259   // Check correctness of the bounding box        273   // Check correctness of the bounding box
260   //                                              274   //
261   if (pMin.x() >= pMax.x() || pMin.y() >= pMax    275   if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
262   {                                               276   {
263     std::ostringstream message;                   277     std::ostringstream message;
264     message << "Bad bounding box (min >= max)     278     message << "Bad bounding box (min >= max) for solid: "
265             << GetName() << " !"                  279             << GetName() << " !"
266             << "\npMin = " << pMin                280             << "\npMin = " << pMin
267             << "\npMax = " << pMax;               281             << "\npMax = " << pMax;
268     G4Exception("G4Trd::BoundingLimits()", "Ge    282     G4Exception("G4Trd::BoundingLimits()", "GeomMgt0001", JustWarning, message);
269     DumpInfo();                                   283     DumpInfo();
270   }                                               284   }
271 }                                                 285 }
272                                                   286 
273 //////////////////////////////////////////////    287 //////////////////////////////////////////////////////////////////////////
274 //                                                288 //
275 // Calculate extent under transform and specif    289 // Calculate extent under transform and specified limit
276                                                   290 
277 G4bool G4Trd::CalculateExtent( const EAxis pAx    291 G4bool G4Trd::CalculateExtent( const EAxis pAxis,
278                                const G4VoxelLi    292                                const G4VoxelLimits& pVoxelLimit,
279                                const G4AffineT    293                                const G4AffineTransform& pTransform,
280                                      G4double&    294                                      G4double& pMin, G4double& pMax ) const
281 {                                                 295 {
282   G4ThreeVector bmin, bmax;                       296   G4ThreeVector bmin, bmax;
283   G4bool exist;                                   297   G4bool exist;
284                                                   298 
285   // Check bounding box (bbox)                    299   // Check bounding box (bbox)
286   //                                              300   //
287   BoundingLimits(bmin,bmax);                      301   BoundingLimits(bmin,bmax);
288   G4BoundingEnvelope bbox(bmin,bmax);             302   G4BoundingEnvelope bbox(bmin,bmax);
289 #ifdef G4BBOX_EXTENT                              303 #ifdef G4BBOX_EXTENT
290   return bbox.CalculateExtent(pAxis,pVoxelLimi << 304   if (true) return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
291 #endif                                            305 #endif
292   if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox    306   if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
293   {                                               307   {
294     return exist = pMin < pMax;                << 308     return exist = (pMin < pMax) ? true : false;
295   }                                               309   }
296                                                   310 
297   // Set bounding envelope (benv) and calculat    311   // Set bounding envelope (benv) and calculate extent
298   //                                              312   //
299   G4double dx1 = GetXHalfLength1();               313   G4double dx1 = GetXHalfLength1();
300   G4double dx2 = GetXHalfLength2();               314   G4double dx2 = GetXHalfLength2();
301   G4double dy1 = GetYHalfLength1();               315   G4double dy1 = GetYHalfLength1();
302   G4double dy2 = GetYHalfLength2();               316   G4double dy2 = GetYHalfLength2();
303   G4double dz  = GetZHalfLength();                317   G4double dz  = GetZHalfLength();
304                                                   318 
305   G4ThreeVectorList baseA(4), baseB(4);           319   G4ThreeVectorList baseA(4), baseB(4);
306   baseA[0].set(-dx1,-dy1,-dz);                    320   baseA[0].set(-dx1,-dy1,-dz);
307   baseA[1].set( dx1,-dy1,-dz);                    321   baseA[1].set( dx1,-dy1,-dz);
308   baseA[2].set( dx1, dy1,-dz);                    322   baseA[2].set( dx1, dy1,-dz);
309   baseA[3].set(-dx1, dy1,-dz);                    323   baseA[3].set(-dx1, dy1,-dz);
310   baseB[0].set(-dx2,-dy2, dz);                    324   baseB[0].set(-dx2,-dy2, dz);
311   baseB[1].set( dx2,-dy2, dz);                    325   baseB[1].set( dx2,-dy2, dz);
312   baseB[2].set( dx2, dy2, dz);                    326   baseB[2].set( dx2, dy2, dz);
313   baseB[3].set(-dx2, dy2, dz);                    327   baseB[3].set(-dx2, dy2, dz);
314                                                   328 
315   std::vector<const G4ThreeVectorList *> polyg    329   std::vector<const G4ThreeVectorList *> polygons(2);
316   polygons[0] = &baseA;                           330   polygons[0] = &baseA;
317   polygons[1] = &baseB;                           331   polygons[1] = &baseB;
318                                                   332 
319   G4BoundingEnvelope benv(bmin,bmax,polygons);    333   G4BoundingEnvelope benv(bmin,bmax,polygons);
320   exist = benv.CalculateExtent(pAxis,pVoxelLim    334   exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
321   return exist;                                   335   return exist;
322 }                                                 336 }
323                                                   337 
324 //////////////////////////////////////////////    338 //////////////////////////////////////////////////////////////////////////
325 //                                                339 //
326 // Return whether point inside/outside/on surf    340 // Return whether point inside/outside/on surface, using tolerance
327                                                   341 
328 EInside G4Trd::Inside( const G4ThreeVector& p     342 EInside G4Trd::Inside( const G4ThreeVector& p ) const
329 {                                              << 343 {  
330   G4double dx = fPlanes[3].a*std::abs(p.x())+f    344   G4double dx = fPlanes[3].a*std::abs(p.x())+fPlanes[3].c*p.z()+fPlanes[3].d;
331   G4double dy = fPlanes[1].b*std::abs(p.y())+f    345   G4double dy = fPlanes[1].b*std::abs(p.y())+fPlanes[1].c*p.z()+fPlanes[1].d;
332   G4double dxy = std::max(dx,dy);                 346   G4double dxy = std::max(dx,dy);
333                                                   347 
334   G4double dz = std::abs(p.z())-fDz;              348   G4double dz = std::abs(p.z())-fDz;
335   G4double dist = std::max(dz,dxy);               349   G4double dist = std::max(dz,dxy);
336                                                   350 
337   return (dist > halfCarTolerance) ? kOutside  << 351   if (dist > halfCarTolerance) return kOutside;
338     ((dist > -halfCarTolerance) ? kSurface : k << 352   return (dist > -halfCarTolerance) ? kSurface : kInside;
339 }                                                 353 }
340                                                   354 
341 //////////////////////////////////////////////    355 //////////////////////////////////////////////////////////////////////////
342 //                                                356 //
343 // Determine side where point is, and return c    357 // Determine side where point is, and return corresponding normal
344                                                   358 
345 G4ThreeVector G4Trd::SurfaceNormal( const G4Th    359 G4ThreeVector G4Trd::SurfaceNormal( const G4ThreeVector& p ) const
346 {                                                 360 {
347   G4int nsurf = 0; // number of surfaces where    361   G4int nsurf = 0; // number of surfaces where p is placed
348                                                   362 
349   // Check Z faces                                363   // Check Z faces
350   //                                              364   //
351   G4double nz = 0;                                365   G4double nz = 0;
352   G4double dz = std::abs(p.z()) - fDz;            366   G4double dz = std::abs(p.z()) - fDz;
353   if (std::abs(dz) <= halfCarTolerance)           367   if (std::abs(dz) <= halfCarTolerance)
354   {                                               368   {
355     nz = (p.z() < 0) ? -1 : 1;                    369     nz = (p.z() < 0) ? -1 : 1;
356     ++nsurf;                                      370     ++nsurf;
357   }                                               371   }
358                                                   372 
359   // Check Y faces                                373   // Check Y faces
360   //                                              374   //
361   G4double ny = 0;                                375   G4double ny = 0;
362   G4double dy1 = fPlanes[0].b*p.y();              376   G4double dy1 = fPlanes[0].b*p.y();
363   G4double dy2 = fPlanes[0].c*p.z() + fPlanes[    377   G4double dy2 = fPlanes[0].c*p.z() + fPlanes[0].d;
364   if (std::abs(dy2 + dy1) <= halfCarTolerance)    378   if (std::abs(dy2 + dy1) <= halfCarTolerance)
365   {                                               379   {
366     ny += fPlanes[0].b;                           380     ny += fPlanes[0].b;
367     nz += fPlanes[0].c;                           381     nz += fPlanes[0].c;
368     ++nsurf;                                      382     ++nsurf;
369   }                                               383   }
370   if (std::abs(dy2 - dy1) <= halfCarTolerance)    384   if (std::abs(dy2 - dy1) <= halfCarTolerance)
371   {                                               385   {
372     ny += fPlanes[1].b;                           386     ny += fPlanes[1].b;
373     nz += fPlanes[1].c;                           387     nz += fPlanes[1].c;
374     ++nsurf;                                      388     ++nsurf;
375   }                                               389   }
376                                                   390 
377   // Check X faces                                391   // Check X faces
378   //                                              392   //
379   G4double nx = 0;                                393   G4double nx = 0;
380   G4double dx1 = fPlanes[2].a*p.x();              394   G4double dx1 = fPlanes[2].a*p.x();
381   G4double dx2 = fPlanes[2].c*p.z() + fPlanes[    395   G4double dx2 = fPlanes[2].c*p.z() + fPlanes[2].d;
382   if (std::abs(dx2 + dx1) <= halfCarTolerance)    396   if (std::abs(dx2 + dx1) <= halfCarTolerance)
383   {                                               397   {
384     nx += fPlanes[2].a;                           398     nx += fPlanes[2].a;
385     nz += fPlanes[2].c;                           399     nz += fPlanes[2].c;
386     ++nsurf;                                      400     ++nsurf;
387   }                                               401   }
388   if (std::abs(dx2 - dx1) <= halfCarTolerance)    402   if (std::abs(dx2 - dx1) <= halfCarTolerance)
389   {                                               403   {
390     nx += fPlanes[3].a;                           404     nx += fPlanes[3].a;
391     nz += fPlanes[3].c;                           405     nz += fPlanes[3].c;
392     ++nsurf;                                      406     ++nsurf;
393   }                                               407   }
394                                                   408 
395   // Return normal                                409   // Return normal
396   //                                              410   //
397   if (nsurf == 1)      return {nx,ny,nz};      << 411   if (nsurf == 1)      return G4ThreeVector(nx,ny,nz);
398   else if (nsurf != 0) return G4ThreeVector(nx    412   else if (nsurf != 0) return G4ThreeVector(nx,ny,nz).unit(); // edge or corner
399   else                                            413   else
400   {                                               414   {
401     // Point is not on the surface                415     // Point is not on the surface
402     //                                            416     //
403 #ifdef G4CSGDEBUG                                 417 #ifdef G4CSGDEBUG
404     std::ostringstream message;                   418     std::ostringstream message;
405     G4long oldprc = message.precision(16);     << 419     G4int oldprc = message.precision(16);
406     message << "Point p is not on surface (!?)    420     message << "Point p is not on surface (!?) of solid: "
407             << GetName() << G4endl;               421             << GetName() << G4endl;
408     message << "Position:\n";                     422     message << "Position:\n";
409     message << "   p.x() = " << p.x()/mm << "     423     message << "   p.x() = " << p.x()/mm << " mm\n";
410     message << "   p.y() = " << p.y()/mm << "     424     message << "   p.y() = " << p.y()/mm << " mm\n";
411     message << "   p.z() = " << p.z()/mm << "     425     message << "   p.z() = " << p.z()/mm << " mm";
412     G4cout.precision(oldprc) ;                    426     G4cout.precision(oldprc) ;
413     G4Exception("G4Trd::SurfaceNormal(p)", "Ge    427     G4Exception("G4Trd::SurfaceNormal(p)", "GeomSolids1002",
414                 JustWarning, message );           428                 JustWarning, message );
415     DumpInfo();                                   429     DumpInfo();
416 #endif                                            430 #endif
417     return ApproxSurfaceNormal(p);                431     return ApproxSurfaceNormal(p);
418   }                                               432   }
419 }                                                 433 }
420                                                   434 
421 //////////////////////////////////////////////    435 //////////////////////////////////////////////////////////////////////////
422 //                                                436 //
423 // Algorithm for SurfaceNormal() following the    437 // Algorithm for SurfaceNormal() following the original specification
424 // for points not on the surface                  438 // for points not on the surface
425                                                   439 
426 G4ThreeVector G4Trd::ApproxSurfaceNormal( cons    440 G4ThreeVector G4Trd::ApproxSurfaceNormal( const G4ThreeVector& p ) const
427 {                                                 441 {
428   G4double dist = -DBL_MAX;                       442   G4double dist = -DBL_MAX;
429   G4int iside = 0;                                443   G4int iside = 0;
430   for (G4int i=0; i<4; ++i)                       444   for (G4int i=0; i<4; ++i)
431   {                                               445   {
432     G4double d = fPlanes[i].a*p.x() +             446     G4double d = fPlanes[i].a*p.x() +
433                  fPlanes[i].b*p.y() +             447                  fPlanes[i].b*p.y() +
434                  fPlanes[i].c*p.z() + fPlanes[    448                  fPlanes[i].c*p.z() + fPlanes[i].d;
435     if (d > dist) { dist = d; iside = i; }        449     if (d > dist) { dist = d; iside = i; }
436   }                                               450   }
437                                                   451 
438   G4double distz = std::abs(p.z()) - fDz;         452   G4double distz = std::abs(p.z()) - fDz;
439   if (dist > distz)                               453   if (dist > distz)
440     return { fPlanes[iside].a, fPlanes[iside]. << 454     return G4ThreeVector(fPlanes[iside].a, fPlanes[iside].b, fPlanes[iside].c);
441   else                                            455   else
442     return { 0, 0, (G4double)((p.z() < 0) ? -1 << 456     return G4ThreeVector(0, 0, (p.z() < 0) ? -1 : 1);
443 }                                                 457 }
444                                                   458 
445 //////////////////////////////////////////////    459 //////////////////////////////////////////////////////////////////////////
446 //                                                460 //
447 // Calculate distance to shape from outside       461 // Calculate distance to shape from outside
448 //  - return kInfinity if no intersection         462 //  - return kInfinity if no intersection
449                                                   463 
450 G4double G4Trd::DistanceToIn(const G4ThreeVect    464 G4double G4Trd::DistanceToIn(const G4ThreeVector& p,
451                              const G4ThreeVect    465                              const G4ThreeVector& v ) const
452 {                                                 466 {
453   // Z intersections                              467   // Z intersections
454   //                                              468   //
455   if ((std::abs(p.z()) - fDz) >= -halfCarToler    469   if ((std::abs(p.z()) - fDz) >= -halfCarTolerance && p.z()*v.z() >= 0)
456     return kInfinity;                             470     return kInfinity;
457   G4double invz = (-v.z() == 0) ? DBL_MAX : -1    471   G4double invz = (-v.z() == 0) ? DBL_MAX : -1./v.z();
458   G4double dz = (invz < 0) ? fDz : -fDz;       << 472   G4double dz = (invz < 0) ? fDz : -fDz; 
459   G4double tzmin = (p.z() + dz)*invz;             473   G4double tzmin = (p.z() + dz)*invz;
460   G4double tzmax = (p.z() - dz)*invz;             474   G4double tzmax = (p.z() - dz)*invz;
461                                                   475 
462   // Y intersections                              476   // Y intersections
463   //                                              477   //
464   G4double tmin0 = tzmin, tmax0 = tzmax;          478   G4double tmin0 = tzmin, tmax0 = tzmax;
465   G4double ya = fPlanes[0].b*v.y(), yb = fPlan    479   G4double ya = fPlanes[0].b*v.y(), yb = fPlanes[0].c*v.z();
466   G4double yc = fPlanes[0].b*p.y(), yd = fPlan    480   G4double yc = fPlanes[0].b*p.y(), yd = fPlanes[0].c*p.z()+fPlanes[0].d;
467   G4double cos0 = yb + ya;                        481   G4double cos0 = yb + ya;
468   G4double dis0 = yd + yc;                        482   G4double dis0 = yd + yc;
469   if (dis0 >= -halfCarTolerance)                  483   if (dis0 >= -halfCarTolerance)
470   {                                               484   {
471     if (cos0 >= 0) return kInfinity;              485     if (cos0 >= 0) return kInfinity;
472     G4double tmp  = -dis0/cos0;                   486     G4double tmp  = -dis0/cos0;
473     if (tmin0 < tmp) tmin0 = tmp;                 487     if (tmin0 < tmp) tmin0 = tmp;
474   }                                               488   }
475   else if (cos0 > 0)                              489   else if (cos0 > 0)
476   {                                               490   {
477     G4double tmp  = -dis0/cos0;                   491     G4double tmp  = -dis0/cos0;
478     if (tmax0 > tmp) tmax0 = tmp;                 492     if (tmax0 > tmp) tmax0 = tmp;
479   }                                               493   }
480                                                   494 
481   G4double tmin1 = tmin0, tmax1 = tmax0;          495   G4double tmin1 = tmin0, tmax1 = tmax0;
482   G4double cos1 = yb - ya;                        496   G4double cos1 = yb - ya;
483   G4double dis1 = yd - yc;                        497   G4double dis1 = yd - yc;
484   if (dis1 >= -halfCarTolerance)                  498   if (dis1 >= -halfCarTolerance)
485   {                                               499   {
486     if (cos1 >= 0) return kInfinity;              500     if (cos1 >= 0) return kInfinity;
487     G4double tmp  = -dis1/cos1;                   501     G4double tmp  = -dis1/cos1;
488     if (tmin1 < tmp) tmin1 = tmp;                 502     if (tmin1 < tmp) tmin1 = tmp;
489   }                                               503   }
490   else if (cos1 > 0)                              504   else if (cos1 > 0)
491   {                                               505   {
492     G4double tmp  = -dis1/cos1;                   506     G4double tmp  = -dis1/cos1;
493     if (tmax1 > tmp) tmax1 = tmp;                 507     if (tmax1 > tmp) tmax1 = tmp;
494   }                                               508   }
495                                                   509 
496   // X intersections                              510   // X intersections
497   //                                              511   //
498   G4double tmin2 = tmin1, tmax2 = tmax1;          512   G4double tmin2 = tmin1, tmax2 = tmax1;
499   G4double xa = fPlanes[2].a*v.x(), xb = fPlan    513   G4double xa = fPlanes[2].a*v.x(), xb = fPlanes[2].c*v.z();
500   G4double xc = fPlanes[2].a*p.x(), xd = fPlan    514   G4double xc = fPlanes[2].a*p.x(), xd = fPlanes[2].c*p.z()+fPlanes[2].d;
501   G4double cos2 = xb + xa;                        515   G4double cos2 = xb + xa;
502   G4double dis2 = xd + xc;                        516   G4double dis2 = xd + xc;
503   if (dis2 >= -halfCarTolerance)                  517   if (dis2 >= -halfCarTolerance)
504   {                                               518   {
505     if (cos2 >= 0) return kInfinity;              519     if (cos2 >= 0) return kInfinity;
506     G4double tmp  = -dis2/cos2;                   520     G4double tmp  = -dis2/cos2;
507     if (tmin2 < tmp) tmin2 = tmp;                 521     if (tmin2 < tmp) tmin2 = tmp;
508   }                                               522   }
509   else if (cos2 > 0)                              523   else if (cos2 > 0)
510   {                                               524   {
511     G4double tmp  = -dis2/cos2;                   525     G4double tmp  = -dis2/cos2;
512     if (tmax2 > tmp) tmax2 = tmp;                 526     if (tmax2 > tmp) tmax2 = tmp;
513   }                                               527   }
514                                                   528 
515   G4double tmin3 = tmin2, tmax3 = tmax2;          529   G4double tmin3 = tmin2, tmax3 = tmax2;
516   G4double cos3 = xb - xa;                        530   G4double cos3 = xb - xa;
517   G4double dis3 = xd - xc;                        531   G4double dis3 = xd - xc;
518   if (dis3 >= -halfCarTolerance)                  532   if (dis3 >= -halfCarTolerance)
519   {                                               533   {
520     if (cos3 >= 0) return kInfinity;              534     if (cos3 >= 0) return kInfinity;
521     G4double tmp  = -dis3/cos3;                   535     G4double tmp  = -dis3/cos3;
522     if (tmin3 < tmp) tmin3 = tmp;                 536     if (tmin3 < tmp) tmin3 = tmp;
523   }                                               537   }
524   else if (cos3 > 0)                              538   else if (cos3 > 0)
525   {                                               539   {
526     G4double tmp  = -dis3/cos3;                   540     G4double tmp  = -dis3/cos3;
527     if (tmax3 > tmp) tmax3 = tmp;                 541     if (tmax3 > tmp) tmax3 = tmp;
528   }                                               542   }
529                                                   543 
530   // Find distance                                544   // Find distance
531   //                                              545   //
532   G4double tmin = tmin3, tmax = tmax3;            546   G4double tmin = tmin3, tmax = tmax3;
533   if (tmax <= tmin + halfCarTolerance) return     547   if (tmax <= tmin + halfCarTolerance) return kInfinity; // touch or no hit
534   return (tmin < halfCarTolerance ) ? 0. : tmi    548   return (tmin < halfCarTolerance ) ? 0. : tmin;
535 }                                                 549 }
536                                                   550 
537 //////////////////////////////////////////////    551 //////////////////////////////////////////////////////////////////////////
538 //                                                552 //
539 // Calculate exact shortest distance to any bo    553 // Calculate exact shortest distance to any boundary from outside
540 // This is the best fast estimation of the sho    554 // This is the best fast estimation of the shortest distance to trap
541 // - returns 0 if point is inside                 555 // - returns 0 if point is inside
542                                                   556 
543 G4double G4Trd::DistanceToIn( const G4ThreeVec    557 G4double G4Trd::DistanceToIn( const G4ThreeVector& p ) const
544 {                                                 558 {
545   G4double dx = fPlanes[3].a*std::abs(p.x())+f    559   G4double dx = fPlanes[3].a*std::abs(p.x())+fPlanes[3].c*p.z()+fPlanes[3].d;
546   G4double dy = fPlanes[1].b*std::abs(p.y())+f    560   G4double dy = fPlanes[1].b*std::abs(p.y())+fPlanes[1].c*p.z()+fPlanes[1].d;
547   G4double dxy = std::max(dx,dy);                 561   G4double dxy = std::max(dx,dy);
548                                                   562 
549   G4double dz = std::abs(p.z())-fDz;              563   G4double dz = std::abs(p.z())-fDz;
550   G4double dist = std::max(dz,dxy);               564   G4double dist = std::max(dz,dxy);
551                                                   565 
552   return (dist > 0) ? dist : 0.;                  566   return (dist > 0) ? dist : 0.;
553 }                                                 567 }
554                                                   568 
555 //////////////////////////////////////////////    569 //////////////////////////////////////////////////////////////////////////
556 //                                                570 //
557 // Calculate distance to surface of shape from    571 // Calculate distance to surface of shape from inside and
558 // find normal at exit point, if required         572 // find normal at exit point, if required
559 // - when leaving the surface, return 0           573 // - when leaving the surface, return 0
560                                                   574 
561 G4double G4Trd::DistanceToOut(const G4ThreeVec    575 G4double G4Trd::DistanceToOut(const G4ThreeVector& p, const G4ThreeVector& v,
562                               const G4bool cal    576                               const G4bool calcNorm,
563                                     G4bool* va << 577                                     G4bool *validNorm, G4ThreeVector *n) const
564 {                                                 578 {
565   // Z intersections                              579   // Z intersections
566   //                                              580   //
567   if ((std::abs(p.z()) - fDz) >= -halfCarToler    581   if ((std::abs(p.z()) - fDz) >= -halfCarTolerance && p.z()*v.z() > 0)
568   {                                               582   {
569     if (calcNorm)                                 583     if (calcNorm)
570     {                                             584     {
571       *validNorm = true;                          585       *validNorm = true;
572       n->set(0, 0, (p.z() < 0) ? -1 : 1);         586       n->set(0, 0, (p.z() < 0) ? -1 : 1);
573     }                                             587     }
574     return 0;                                     588     return 0;
575   }                                               589   }
576   G4double vz = v.z();                            590   G4double vz = v.z();
577   G4double tmax = (vz == 0) ? DBL_MAX : (std::    591   G4double tmax = (vz == 0) ? DBL_MAX : (std::copysign(fDz,vz) - p.z())/vz;
578   G4int iside = (vz < 0) ? -4 : -2; // little     592   G4int iside = (vz < 0) ? -4 : -2; // little trick: (-4+3)=-1, (-2+3)=+1
579                                                   593 
580   // Y intersections                              594   // Y intersections
581   //                                              595   //
582   G4int i = 0;                                    596   G4int i = 0;
583   for ( ; i<2; ++i)                               597   for ( ; i<2; ++i)
584   {                                               598   {
585     G4double cosa = fPlanes[i].b*v.y() + fPlan    599     G4double cosa = fPlanes[i].b*v.y() + fPlanes[i].c*v.z();
586     if (cosa > 0)                                 600     if (cosa > 0)
587     {                                             601     {
588       G4double dist = fPlanes[i].b*p.y()+fPlan    602       G4double dist = fPlanes[i].b*p.y()+fPlanes[i].c*p.z()+fPlanes[i].d;
589       if (dist >= -halfCarTolerance)              603       if (dist >= -halfCarTolerance)
590       {                                           604       {
591         if (calcNorm)                             605         if (calcNorm)
592         {                                         606         {
593           *validNorm = true;                      607           *validNorm = true;
594           n->set(0, fPlanes[i].b, fPlanes[i].c    608           n->set(0, fPlanes[i].b, fPlanes[i].c);
595         }                                         609         }
596         return 0;                                 610         return 0;
597       }                                           611       }
598       G4double tmp = -dist/cosa;                  612       G4double tmp = -dist/cosa;
599       if (tmax > tmp) { tmax = tmp; iside = i;    613       if (tmax > tmp) { tmax = tmp; iside = i; }
600     }                                             614     }
601   }                                               615   }
602                                                   616 
603   // X intersections                              617   // X intersections
604   //                                              618   //
605   for ( ; i<4; ++i)                               619   for ( ; i<4; ++i)
606   {                                               620   {
607     G4double cosa = fPlanes[i].a*v.x()+fPlanes    621     G4double cosa = fPlanes[i].a*v.x()+fPlanes[i].c*v.z();
608     if (cosa > 0)                                 622     if (cosa > 0)
609     {                                             623     {
610       G4double dist = fPlanes[i].a*p.x()+fPlan    624       G4double dist = fPlanes[i].a*p.x()+fPlanes[i].c*p.z()+fPlanes[i].d;
611       if (dist >= -halfCarTolerance)              625       if (dist >= -halfCarTolerance)
612       {                                           626       {
613         if (calcNorm)                             627         if (calcNorm)
614         {                                         628         {
615            *validNorm = true;                     629            *validNorm = true;
616            n->set(fPlanes[i].a, fPlanes[i].b,     630            n->set(fPlanes[i].a, fPlanes[i].b, fPlanes[i].c);
617         }                                         631         }
618         return 0;                                 632         return 0;
619       }                                           633       }
620       G4double tmp = -dist/cosa;                  634       G4double tmp = -dist/cosa;
621       if (tmax > tmp) { tmax = tmp; iside = i;    635       if (tmax > tmp) { tmax = tmp; iside = i; }
622     }                                             636     }
623   }                                               637   }
624                                                   638 
625   // Set normal, if required, and return dista    639   // Set normal, if required, and return distance
626   //                                              640   //
627   if (calcNorm)                                << 641   if (calcNorm) 
628   {                                               642   {
629     *validNorm = true;                            643     *validNorm = true;
630     if (iside < 0)                                644     if (iside < 0)
631       n->set(0, 0, iside + 3); // (-4+3)=-1, (    645       n->set(0, 0, iside + 3); // (-4+3)=-1, (-2+3)=+1
632     else                                          646     else
633       n->set(fPlanes[iside].a, fPlanes[iside].    647       n->set(fPlanes[iside].a, fPlanes[iside].b, fPlanes[iside].c);
634   }                                               648   }
635   return tmax;                                    649   return tmax;
636 }                                                 650 }
637                                                   651 
638 //////////////////////////////////////////////    652 //////////////////////////////////////////////////////////////////////////
639 //                                                653 //
640 // Calculate exact shortest distance to any bo    654 // Calculate exact shortest distance to any boundary from inside
641 // - returns 0 if point is outside                655 // - returns 0 if point is outside
642                                                   656 
643 G4double G4Trd::DistanceToOut( const G4ThreeVe    657 G4double G4Trd::DistanceToOut( const G4ThreeVector& p ) const
644 {                                                 658 {
645 #ifdef G4CSGDEBUG                                 659 #ifdef G4CSGDEBUG
646   if( Inside(p) == kOutside )                     660   if( Inside(p) == kOutside )
647   {                                               661   {
648     std::ostringstream message;                   662     std::ostringstream message;
649     G4long oldprc = message.precision(16);     << 663     G4int oldprc = message.precision(16);
650     message << "Point p is outside (!?) of sol    664     message << "Point p is outside (!?) of solid: " << GetName() << G4endl;
651     message << "Position:\n";                     665     message << "Position:\n";
652     message << "   p.x() = " << p.x()/mm << "     666     message << "   p.x() = " << p.x()/mm << " mm\n";
653     message << "   p.y() = " << p.y()/mm << "     667     message << "   p.y() = " << p.y()/mm << " mm\n";
654     message << "   p.z() = " << p.z()/mm << "     668     message << "   p.z() = " << p.z()/mm << " mm";
655     G4cout.precision(oldprc);                  << 669     G4cout.precision(oldprc) ;
656     G4Exception("G4Trd::DistanceToOut(p)", "Ge    670     G4Exception("G4Trd::DistanceToOut(p)", "GeomSolids1002",
657                 JustWarning, message );           671                 JustWarning, message );
658     DumpInfo();                                   672     DumpInfo();
659   }                                               673   }
660 #endif                                            674 #endif
661   G4double dx = fPlanes[3].a*std::abs(p.x())+f    675   G4double dx = fPlanes[3].a*std::abs(p.x())+fPlanes[3].c*p.z()+fPlanes[3].d;
662   G4double dy = fPlanes[1].b*std::abs(p.y())+f    676   G4double dy = fPlanes[1].b*std::abs(p.y())+fPlanes[1].c*p.z()+fPlanes[1].d;
663   G4double dxy = std::max(dx,dy);                 677   G4double dxy = std::max(dx,dy);
664                                                   678 
665   G4double dz = std::abs(p.z())-fDz;              679   G4double dz = std::abs(p.z())-fDz;
666   G4double dist = std::max(dz,dxy);               680   G4double dist = std::max(dz,dxy);
667                                                << 681   
668   return (dist < 0) ? -dist : 0.;                 682   return (dist < 0) ? -dist : 0.;
669 }                                                 683 }
670                                                   684 
671 //////////////////////////////////////////////    685 //////////////////////////////////////////////////////////////////////////
672 //                                                686 //
673 // GetEntityType                                  687 // GetEntityType
674                                                   688 
675 G4GeometryType G4Trd::GetEntityType() const       689 G4GeometryType G4Trd::GetEntityType() const
676 {                                                 690 {
677   return {"G4Trd"};                            << 691   return G4String("G4Trd");
678 }                                              << 
679                                                << 
680 ////////////////////////////////////////////// << 
681 //                                             << 
682 // IsFaceted                                   << 
683                                                << 
684 G4bool G4Trd::IsFaceted() const                << 
685 {                                              << 
686   return true;                                 << 
687 }                                                 692 }
688                                                   693 
689 //////////////////////////////////////////////    694 //////////////////////////////////////////////////////////////////////////
690 //                                                695 //
691 // Make a clone of the object                     696 // Make a clone of the object
692 //                                                697 //
693 G4VSolid* G4Trd::Clone() const                    698 G4VSolid* G4Trd::Clone() const
694 {                                                 699 {
695   return new G4Trd(*this);                        700   return new G4Trd(*this);
696 }                                                 701 }
697                                                   702 
698 //////////////////////////////////////////////    703 //////////////////////////////////////////////////////////////////////////
699 //                                                704 //
700 // Stream object contents to an output stream     705 // Stream object contents to an output stream
701                                                   706 
702 std::ostream& G4Trd::StreamInfo( std::ostream&    707 std::ostream& G4Trd::StreamInfo( std::ostream& os ) const
703 {                                                 708 {
704   G4long oldprc = os.precision(16);            << 709   G4int oldprc = os.precision(16);
705   os << "-------------------------------------    710   os << "-----------------------------------------------------------\n"
706      << "    *** Dump for solid - " << GetName    711      << "    *** Dump for solid - " << GetName() << " ***\n"
707      << "    =================================    712      << "    ===================================================\n"
708      << " Solid type: G4Trd\n"                    713      << " Solid type: G4Trd\n"
709      << " Parameters: \n"                         714      << " Parameters: \n"
710      << "    half length X, surface -dZ: " <<     715      << "    half length X, surface -dZ: " << fDx1/mm << " mm \n"
711      << "    half length X, surface +dZ: " <<     716      << "    half length X, surface +dZ: " << fDx2/mm << " mm \n"
712      << "    half length Y, surface -dZ: " <<     717      << "    half length Y, surface -dZ: " << fDy1/mm << " mm \n"
713      << "    half length Y, surface +dZ: " <<     718      << "    half length Y, surface +dZ: " << fDy2/mm << " mm \n"
714      << "    half length Z             : " <<     719      << "    half length Z             : " <<  fDz/mm << " mm \n"
715      << "-------------------------------------    720      << "-----------------------------------------------------------\n";
716   os.precision(oldprc);                           721   os.precision(oldprc);
717                                                   722 
718   return os;                                      723   return os;
719 }                                                 724 }
720                                                   725 
721 //////////////////////////////////////////////    726 //////////////////////////////////////////////////////////////////////////
722 //                                                727 //
723 // Return a point randomly and uniformly selec    728 // Return a point randomly and uniformly selected on the solid surface
724                                                   729 
725 G4ThreeVector G4Trd::GetPointOnSurface() const    730 G4ThreeVector G4Trd::GetPointOnSurface() const
726 {                                                 731 {
727   // Set areas                                 << 732   // Set vertices
728   //                                              733   //
729   G4double sxz = (fDx1 + fDx2)*fHx;            << 734   G4ThreeVector pt[8];
730   G4double syz = (fDy1 + fDy2)*fHy;            << 735   pt[0].set(-fDx1,-fDy1,-fDz);
731   G4double ssurf[6] = { 4.*fDx1*fDy1, sxz, sxz << 736   pt[1].set( fDx1,-fDy1,-fDz);
732   ssurf[1] += ssurf[0];                        << 737   pt[2].set(-fDx1, fDy1,-fDz);
733   ssurf[2] += ssurf[1];                        << 738   pt[3].set( fDx1, fDy1,-fDz);
734   ssurf[3] += ssurf[2];                        << 739   pt[4].set(-fDx2,-fDy2, fDz);
735   ssurf[4] += ssurf[3];                        << 740   pt[5].set( fDx2,-fDy2, fDz);
736   ssurf[5] += ssurf[4];                        << 741   pt[6].set(-fDx2, fDy2, fDz);
                                                   >> 742   pt[7].set( fDx2, fDy2, fDz);
                                                   >> 743 
                                                   >> 744   // Set faces  (-Z, -Y, +Y, -X, +X, +Z)
                                                   >> 745   //
                                                   >> 746   G4int iface [6][4] =
                                                   >> 747     { {0,1,3,2}, {0,4,5,1}, {2,3,7,6}, {0,2,6,4}, {1,5,7,3}, {4,6,7,5} };
                                                   >> 748 
                                                   >> 749   // Set areas
                                                   >> 750   // 
                                                   >> 751   G4double sxz = (fDy1 + fDy2)*std::hypot(fDx1 - fDx2, 2*fDz);
                                                   >> 752   G4double syz = (fDx1 + fDx2)*std::hypot(fDy1 - fDy2, 2*fDz);
                                                   >> 753   G4double sface[6] = { 4*fDx1*fDy1, syz, syz, sxz, sxz, 4*fDx2*fDy2 };
                                                   >> 754   for (G4int i=1; i<6; ++i) { sface[i] += sface[i-1]; }
737                                                   755 
738   // Select face                                  756   // Select face
739   //                                              757   //
740   G4double select = ssurf[5]*G4QuickRand();    << 758   G4double select = sface[5]*G4UniformRand();
741   G4int k = 5;                                    759   G4int k = 5;
742   k -= (G4int)(select <= ssurf[4]);            << 760   if (select <= sface[4]) k = 4;
743   k -= (G4int)(select <= ssurf[3]);            << 761   if (select <= sface[3]) k = 3;
744   k -= (G4int)(select <= ssurf[2]);            << 762   if (select <= sface[2]) k = 2;
745   k -= (G4int)(select <= ssurf[1]);            << 763   if (select <= sface[1]) k = 1;
746   k -= (G4int)(select <= ssurf[0]);            << 764   if (select <= sface[0]) k = 0;
747                                                << 765 
748   // Generate point on selected surface        << 766   // Select sub-triangle
749   //                                           << 767   //
750   G4double u = G4QuickRand();                  << 768   G4int i0 = iface[k][0];
751   G4double v = G4QuickRand();                  << 769   G4int i1 = iface[k][1];
752   switch(k)                                    << 770   G4int i2 = iface[k][2];
753   {                                            << 771   G4int i3 = iface[k][3];
754     case 0: // base at -Z                      << 772   G4double s1 = G4GeomTools::TriangleAreaNormal(pt[i0],pt[i1],pt[i3]).mag();
755     {                                          << 773   G4double s2 = G4GeomTools::TriangleAreaNormal(pt[i2],pt[i1],pt[i3]).mag();
756       return { (2.*u - 1.)*fDx1, (2.*v - 1.)*f << 774   if ((s1+s2)*G4UniformRand() > s1) i0 = i2;
757     }                                          << 775 
758     case 1: // X face at -Y                    << 776   // Generate point
759     {                                          << 777   //
760       if (u + v > 1.) { u = 1. - u; v = 1. - v << 778   G4double u = G4UniformRand();
761       G4ThreeVector p0(-fDx1,-fDy1,-fDz);      << 779   G4double v = G4UniformRand();
762       G4ThreeVector p1( fDx2,-fDy2, fDz);      << 780   if (u + v > 1.) { u = 1. - u; v = 1. - v; }
763       return (select <= ssurf[0] + fDx1*fHx) ? << 781   return (1.-u-v)*pt[i0] + u*pt[i1] + v*pt[i3];
764         (1. - u - v)*p0 + u*p1 + v*G4ThreeVect << 
765         (1. - u - v)*p0 + u*p1 + v*G4ThreeVect << 
766     }                                          << 
767     case 2: // X face at +Y                    << 
768     {                                          << 
769       if (u + v > 1.) { u = 1. - u; v = 1. - v << 
770       G4ThreeVector p0( fDx1, fDy1,-fDz);      << 
771       G4ThreeVector p1(-fDx2, fDy2, fDz);      << 
772       return (select <= ssurf[1] + fDx1*fHx) ? << 
773         (1. - u - v)*p0 + u*p1 + v*G4ThreeVect << 
774         (1. - u - v)*p0 + u*p1 + v*G4ThreeVect << 
775     }                                          << 
776     case 3: // Y face at -X                    << 
777     {                                          << 
778       if (u + v > 1.) { u = 1. - u; v = 1. - v << 
779       G4ThreeVector p0(-fDx1, fDy1,-fDz);      << 
780       G4ThreeVector p1(-fDx2,-fDy2, fDz);      << 
781       return (select <= ssurf[2] + fDy1*fHy) ? << 
782         (1. - u - v)*p0 + u*p1 + v*G4ThreeVect << 
783         (1. - u - v)*p0 + u*p1 + v*G4ThreeVect << 
784     }                                          << 
785     case 4: // Y face at +X                    << 
786     {                                          << 
787       if (u + v > 1.) { u = 1. - u; v = 1. - v << 
788       G4ThreeVector p0( fDx1,-fDy1,-fDz);      << 
789       G4ThreeVector p1( fDx2, fDy2, fDz);      << 
790       return (select <= ssurf[3] + fDy1*fHy) ? << 
791         (1. - u - v)*p0 + u*p1 + v*G4ThreeVect << 
792         (1. - u - v)*p0 + u*p1 + v*G4ThreeVect << 
793     }                                          << 
794     case 5: // base at +Z                      << 
795     {                                          << 
796       return { (2.*u - 1.)*fDx2, (2.*v - 1.)*f << 
797     }                                          << 
798   }                                            << 
799   return {0., 0., 0.};                         << 
800 }                                                 782 }
801                                                   783 
802 //////////////////////////////////////////////    784 //////////////////////////////////////////////////////////////////////////
803 //                                                785 //
804 // Methods for visualisation                      786 // Methods for visualisation
805                                                   787 
806 void G4Trd::DescribeYourselfTo ( G4VGraphicsSc    788 void G4Trd::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
807 {                                                 789 {
808   scene.AddSolid (*this);                         790   scene.AddSolid (*this);
809 }                                                 791 }
810                                                   792 
811 G4Polyhedron* G4Trd::CreatePolyhedron () const    793 G4Polyhedron* G4Trd::CreatePolyhedron () const
812 {                                                 794 {
813   return new G4PolyhedronTrd2 (fDx1, fDx2, fDy    795   return new G4PolyhedronTrd2 (fDx1, fDx2, fDy1, fDy2, fDz);
814 }                                                 796 }
815                                                   797 
816 #endif                                            798 #endif
817                                                   799