Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/solids/CSG/src/G4Trap.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/G4Trap.cc (Version 11.3.0) and /geometry/solids/CSG/src/G4Trap.cc (Version 11.0.p4)


  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 // Implementation for G4Trap class                 26 // Implementation for G4Trap class
 27 //                                                 27 //
 28 // 21.03.95 P.Kent: Modified for `tolerant' ge     28 // 21.03.95 P.Kent: Modified for `tolerant' geometry
 29 // 09.09.96 V.Grichine: Final modifications be     29 // 09.09.96 V.Grichine: Final modifications before to commit
 30 // 08.12.97 J.Allison: Added "nominal" constru     30 // 08.12.97 J.Allison: Added "nominal" constructor and method SetAllParameters
 31 // 28.04.05 V.Grichine: new SurfaceNormal acco     31 // 28.04.05 V.Grichine: new SurfaceNormal according to J.Apostolakis proposal
 32 // 18.04.17 E.Tcherniaev: complete revision, s     32 // 18.04.17 E.Tcherniaev: complete revision, speed-up
 33 // -------------------------------------------     33 // --------------------------------------------------------------------
 34                                                    34 
 35 #include "G4Trap.hh"                               35 #include "G4Trap.hh"
 36                                                    36 
 37 #if !defined(G4GEOM_USE_UTRAP)                     37 #if !defined(G4GEOM_USE_UTRAP)
 38                                                    38 
 39 #include "globals.hh"                              39 #include "globals.hh"
 40 #include "G4GeomTools.hh"                          40 #include "G4GeomTools.hh"
 41                                                    41 
 42 #include "G4VoxelLimits.hh"                        42 #include "G4VoxelLimits.hh"
 43 #include "G4AffineTransform.hh"                    43 #include "G4AffineTransform.hh"
 44 #include "G4BoundingEnvelope.hh"                   44 #include "G4BoundingEnvelope.hh"
 45                                                    45 
 46 #include "G4VPVParameterisation.hh"                46 #include "G4VPVParameterisation.hh"
 47                                                    47 
 48 #include "G4QuickRand.hh"                          48 #include "G4QuickRand.hh"
 49                                                    49 
 50 #include "G4VGraphicsScene.hh"                     50 #include "G4VGraphicsScene.hh"
 51 #include "G4Polyhedron.hh"                         51 #include "G4Polyhedron.hh"
 52                                                    52 
 53 using namespace CLHEP;                             53 using namespace CLHEP;
 54                                                    54 
 55 //////////////////////////////////////////////     55 //////////////////////////////////////////////////////////////////////////
 56 //                                                 56 //
 57 // Constructor - check and set half-widths as      57 // Constructor - check and set half-widths as well as angles:
 58 // final check of coplanarity                      58 // final check of coplanarity
 59                                                    59 
 60 G4Trap::G4Trap( const G4String& pName,             60 G4Trap::G4Trap( const G4String& pName,
 61                       G4double pDz,                61                       G4double pDz,
 62                       G4double pTheta, G4doubl     62                       G4double pTheta, G4double pPhi,
 63                       G4double pDy1, G4double      63                       G4double pDy1, G4double pDx1, G4double pDx2,
 64                       G4double pAlp1,              64                       G4double pAlp1,
 65                       G4double pDy2, G4double      65                       G4double pDy2, G4double pDx3, G4double pDx4,
 66                       G4double pAlp2 )         <<  66                       G4double pAlp2)
 67   : G4CSGSolid(pName), halfCarTolerance(0.5*kC     67   : G4CSGSolid(pName), halfCarTolerance(0.5*kCarTolerance)
 68 {                                                  68 {
 69   fDz = pDz;                                       69   fDz = pDz;
 70   fTthetaCphi = std::tan(pTheta)*std::cos(pPhi     70   fTthetaCphi = std::tan(pTheta)*std::cos(pPhi);
 71   fTthetaSphi = std::tan(pTheta)*std::sin(pPhi     71   fTthetaSphi = std::tan(pTheta)*std::sin(pPhi);
 72                                                    72 
 73   fDy1 = pDy1; fDx1 = pDx1; fDx2 = pDx2; fTalp     73   fDy1 = pDy1; fDx1 = pDx1; fDx2 = pDx2; fTalpha1 = std::tan(pAlp1);
 74   fDy2 = pDy2; fDx3 = pDx3; fDx4 = pDx4; fTalp     74   fDy2 = pDy2; fDx3 = pDx3; fDx4 = pDx4; fTalpha2 = std::tan(pAlp2);
 75                                                    75 
 76   CheckParameters();                               76   CheckParameters();
 77   MakePlanes();                                    77   MakePlanes();
 78 }                                                  78 }
 79                                                    79 
 80 //////////////////////////////////////////////     80 //////////////////////////////////////////////////////////////////////////
 81 //                                                 81 //
 82 // Constructor - Design of trapezoid based on      82 // Constructor - Design of trapezoid based on 8 G4ThreeVector parameters,
 83 // which are its vertices. Checking of planari     83 // which are its vertices. Checking of planarity with preparation of
 84 // fPlanes[] and than calculation of other mem     84 // fPlanes[] and than calculation of other members
 85                                                    85 
 86 G4Trap::G4Trap( const G4String& pName,             86 G4Trap::G4Trap( const G4String& pName,
 87                 const G4ThreeVector pt[8] )        87                 const G4ThreeVector pt[8] )
 88   : G4CSGSolid(pName), halfCarTolerance(0.5*kC     88   : G4CSGSolid(pName), halfCarTolerance(0.5*kCarTolerance)
 89 {                                                  89 {
 90   // Start with check of centering - the cente     90   // Start with check of centering - the center of gravity trap line
 91   // should cross the origin of frame              91   // should cross the origin of frame
 92   //                                               92   //
 93   if (  pt[0].z() >= 0                         <<  93   if (!(   pt[0].z() < 0
 94         || pt[0].z() != pt[1].z()              <<  94         && pt[0].z() == pt[1].z()
 95         || pt[0].z() != pt[2].z()              <<  95         && pt[0].z() == pt[2].z()
 96         || pt[0].z() != pt[3].z()              <<  96         && pt[0].z() == pt[3].z()
 97                                                <<  97 
 98         || pt[4].z() <= 0                      <<  98         && pt[4].z() > 0
 99         || pt[4].z() != pt[5].z()              <<  99         && pt[4].z() == pt[5].z()
100         || pt[4].z() != pt[6].z()              << 100         && pt[4].z() == pt[6].z()
101         || pt[4].z() != pt[7].z()              << 101         && pt[4].z() == pt[7].z()
102                                                << 102 
103         || std::fabs( pt[0].z() + pt[4].z() )  << 103         && std::fabs( pt[0].z() + pt[4].z() ) < kCarTolerance
104                                                << 104 
105         || pt[0].y() != pt[1].y()              << 105         && pt[0].y() == pt[1].y()
106         || pt[2].y() != pt[3].y()              << 106         && pt[2].y() == pt[3].y()
107         || pt[4].y() != pt[5].y()              << 107         && pt[4].y() == pt[5].y()
108         || pt[6].y() != pt[7].y()              << 108         && pt[6].y() == pt[7].y()
109                                                << 109 
110         || std::fabs(pt[0].y()+pt[2].y()+pt[4] << 110         && std::fabs(pt[0].y()+pt[2].y()+pt[4].y()+pt[6].y()) < kCarTolerance
111         || std::fabs(pt[0].x()+pt[1].x()+pt[4] << 111         && std::fabs(pt[0].x()+pt[1].x()+pt[4].x()+pt[5].x() +
112                      pt[2].x()+pt[3].x()+pt[6] << 112                      pt[2].x()+pt[3].x()+pt[6].x()+pt[7].x()) < kCarTolerance ))
113   {                                               113   {
114     std::ostringstream message;                   114     std::ostringstream message;
115     message << "Invalid vertice coordinates fo    115     message << "Invalid vertice coordinates for Solid: " << GetName();
116     G4Exception("G4Trap::G4Trap()", "GeomSolid    116     G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
117                 FatalException, message);         117                 FatalException, message);
118   }                                               118   }
119                                                   119 
120   // Set parameters                               120   // Set parameters
121   //                                              121   //
122   fDz = (pt[7]).z();                              122   fDz = (pt[7]).z();
123                                                   123 
124   fDy1     = ((pt[2]).y()-(pt[1]).y())*0.5;       124   fDy1     = ((pt[2]).y()-(pt[1]).y())*0.5;
125   fDx1     = ((pt[1]).x()-(pt[0]).x())*0.5;       125   fDx1     = ((pt[1]).x()-(pt[0]).x())*0.5;
126   fDx2     = ((pt[3]).x()-(pt[2]).x())*0.5;       126   fDx2     = ((pt[3]).x()-(pt[2]).x())*0.5;
127   fTalpha1 = ((pt[2]).x()+(pt[3]).x()-(pt[1]).    127   fTalpha1 = ((pt[2]).x()+(pt[3]).x()-(pt[1]).x()-(pt[0]).x())*0.25/fDy1;
128                                                   128 
129   fDy2     = ((pt[6]).y()-(pt[5]).y())*0.5;       129   fDy2     = ((pt[6]).y()-(pt[5]).y())*0.5;
130   fDx3     = ((pt[5]).x()-(pt[4]).x())*0.5;       130   fDx3     = ((pt[5]).x()-(pt[4]).x())*0.5;
131   fDx4     = ((pt[7]).x()-(pt[6]).x())*0.5;       131   fDx4     = ((pt[7]).x()-(pt[6]).x())*0.5;
132   fTalpha2 = ((pt[6]).x()+(pt[7]).x()-(pt[5]).    132   fTalpha2 = ((pt[6]).x()+(pt[7]).x()-(pt[5]).x()-(pt[4]).x())*0.25/fDy2;
133                                                   133 
134   fTthetaCphi = ((pt[4]).x()+fDy2*fTalpha2+fDx    134   fTthetaCphi = ((pt[4]).x()+fDy2*fTalpha2+fDx3)/fDz;
135   fTthetaSphi = ((pt[4]).y()+fDy2)/fDz;           135   fTthetaSphi = ((pt[4]).y()+fDy2)/fDz;
136                                                   136 
137   CheckParameters();                              137   CheckParameters();
138   MakePlanes(pt);                                 138   MakePlanes(pt);
139 }                                                 139 }
140                                                   140 
141 //////////////////////////////////////////////    141 //////////////////////////////////////////////////////////////////////////
142 //                                                142 //
143 // Constructor for Right Angular Wedge from ST    143 // Constructor for Right Angular Wedge from STEP
144                                                   144 
145 G4Trap::G4Trap( const G4String& pName,            145 G4Trap::G4Trap( const G4String& pName,
146                       G4double pZ,                146                       G4double pZ,
147                       G4double pY,                147                       G4double pY,
148                       G4double pX, G4double pL    148                       G4double pX, G4double pLTX )
149   : G4CSGSolid(pName), halfCarTolerance(0.5*kC    149   : G4CSGSolid(pName), halfCarTolerance(0.5*kCarTolerance)
150 {                                                 150 {
151   fDz  = 0.5*pZ; fTthetaCphi = 0; fTthetaSphi     151   fDz  = 0.5*pZ; fTthetaCphi = 0; fTthetaSphi = 0;
152   fDy1 = 0.5*pY; fDx1 = 0.5*pX; fDx2 = 0.5*pLT    152   fDy1 = 0.5*pY; fDx1 = 0.5*pX; fDx2 = 0.5*pLTX; fTalpha1 = 0.5*(pLTX - pX)/pY;
153   fDy2 = fDy1;   fDx3 = fDx1;   fDx4 = fDx2;      153   fDy2 = fDy1;   fDx3 = fDx1;   fDx4 = fDx2;     fTalpha2 = fTalpha1;
154                                                   154 
155   CheckParameters();                              155   CheckParameters();
156   MakePlanes();                                   156   MakePlanes();
157 }                                                 157 }
158                                                   158 
159 //////////////////////////////////////////////    159 //////////////////////////////////////////////////////////////////////////
160 //                                                160 //
161 // Constructor for G4Trd                          161 // Constructor for G4Trd
162                                                   162 
163 G4Trap::G4Trap( const G4String& pName,            163 G4Trap::G4Trap( const G4String& pName,
164                       G4double pDx1,  G4double    164                       G4double pDx1,  G4double pDx2,
165                       G4double pDy1,  G4double    165                       G4double pDy1,  G4double pDy2,
166                       G4double pDz )              166                       G4double pDz )
167   : G4CSGSolid(pName), halfCarTolerance(0.5*kC    167   : G4CSGSolid(pName), halfCarTolerance(0.5*kCarTolerance), fTrapType(0)
168 {                                                 168 {
169   fDz  = pDz;  fTthetaCphi = 0; fTthetaSphi =     169   fDz  = pDz;  fTthetaCphi = 0; fTthetaSphi = 0;
170   fDy1 = pDy1; fDx1 = pDx1; fDx2 = pDx1; fTalp    170   fDy1 = pDy1; fDx1 = pDx1; fDx2 = pDx1; fTalpha1 = 0;
171   fDy2 = pDy2; fDx3 = pDx2; fDx4 = pDx2; fTalp    171   fDy2 = pDy2; fDx3 = pDx2; fDx4 = pDx2; fTalpha2 = 0;
172                                                   172 
173   CheckParameters();                              173   CheckParameters();
174   MakePlanes();                                   174   MakePlanes();
175 }                                                 175 }
176                                                   176 
177 //////////////////////////////////////////////    177 //////////////////////////////////////////////////////////////////////////
178 //                                                178 //
179 // Constructor for G4Para                         179 // Constructor for G4Para
180                                                   180 
181 G4Trap::G4Trap( const G4String& pName,            181 G4Trap::G4Trap( const G4String& pName,
182                       G4double pDx, G4double p    182                       G4double pDx, G4double pDy,
183                       G4double pDz,               183                       G4double pDz,
184                       G4double pAlpha,            184                       G4double pAlpha,
185                       G4double pTheta, G4doubl    185                       G4double pTheta, G4double pPhi )
186   : G4CSGSolid(pName), halfCarTolerance(0.5*kC    186   : G4CSGSolid(pName), halfCarTolerance(0.5*kCarTolerance)
187 {                                                 187 {
188   fDz = pDz;                                      188   fDz = pDz;
189   fTthetaCphi = std::tan(pTheta)*std::cos(pPhi    189   fTthetaCphi = std::tan(pTheta)*std::cos(pPhi);
190   fTthetaSphi = std::tan(pTheta)*std::sin(pPhi    190   fTthetaSphi = std::tan(pTheta)*std::sin(pPhi);
191                                                   191 
192   fDy1 = pDy; fDx1 = pDx; fDx2 = pDx; fTalpha1    192   fDy1 = pDy; fDx1 = pDx; fDx2 = pDx; fTalpha1 = std::tan(pAlpha);
193   fDy2 = pDy; fDx3 = pDx; fDx4 = pDx; fTalpha2    193   fDy2 = pDy; fDx3 = pDx; fDx4 = pDx; fTalpha2 = fTalpha1;
194                                                   194 
195   CheckParameters();                              195   CheckParameters();
196   MakePlanes();                                   196   MakePlanes();
197 }                                                 197 }
198                                                   198 
199 //////////////////////////////////////////////    199 //////////////////////////////////////////////////////////////////////////
200 //                                                200 //
201 // Nominal constructor for G4Trap whose parame    201 // Nominal constructor for G4Trap whose parameters are to be set by
202 // a G4VParamaterisation later.  Check and set    202 // a G4VParamaterisation later.  Check and set half-widths as well as
203 // angles: final check of coplanarity             203 // angles: final check of coplanarity
204                                                   204 
205 G4Trap::G4Trap( const G4String& pName )           205 G4Trap::G4Trap( const G4String& pName )
206   : G4CSGSolid (pName), halfCarTolerance(0.5*k    206   : G4CSGSolid (pName), halfCarTolerance(0.5*kCarTolerance),
207     fDz(1.), fTthetaCphi(0.), fTthetaSphi(0.),    207     fDz(1.), fTthetaCphi(0.), fTthetaSphi(0.),
208     fDy1(1.), fDx1(1.), fDx2(1.), fTalpha1(0.)    208     fDy1(1.), fDx1(1.), fDx2(1.), fTalpha1(0.),
209     fDy2(1.), fDx3(1.), fDx4(1.), fTalpha2(0.)    209     fDy2(1.), fDx3(1.), fDx4(1.), fTalpha2(0.)
210 {                                                 210 {
211   MakePlanes();                                   211   MakePlanes();
212 }                                                 212 }
213                                                   213 
214 //////////////////////////////////////////////    214 //////////////////////////////////////////////////////////////////////////
215 //                                                215 //
216 // Fake default constructor - sets only member    216 // Fake default constructor - sets only member data and allocates memory
217 //                            for usage restri    217 //                            for usage restricted to object persistency.
218 //                                                218 //
219 G4Trap::G4Trap( __void__& a )                     219 G4Trap::G4Trap( __void__& a )
220   : G4CSGSolid(a), halfCarTolerance(0.5*kCarTo    220   : G4CSGSolid(a), halfCarTolerance(0.5*kCarTolerance),
221     fDz(1.), fTthetaCphi(0.), fTthetaSphi(0.),    221     fDz(1.), fTthetaCphi(0.), fTthetaSphi(0.),
222     fDy1(1.), fDx1(1.), fDx2(1.), fTalpha1(0.)    222     fDy1(1.), fDx1(1.), fDx2(1.), fTalpha1(0.),
223     fDy2(1.), fDx3(1.), fDx4(1.), fTalpha2(0.)    223     fDy2(1.), fDx3(1.), fDx4(1.), fTalpha2(0.)
224 {                                                 224 {
225   MakePlanes();                                   225   MakePlanes();
226 }                                                 226 }
227                                                   227 
228 //////////////////////////////////////////////    228 //////////////////////////////////////////////////////////////////////////
229 //                                                229 //
230 // Destructor                                     230 // Destructor
231                                                   231 
232 G4Trap::~G4Trap() = default;                   << 232 G4Trap::~G4Trap()
                                                   >> 233 {
                                                   >> 234 }
233                                                   235 
234 //////////////////////////////////////////////    236 //////////////////////////////////////////////////////////////////////////
235 //                                                237 //
236 // Copy constructor                               238 // Copy constructor
237                                                   239 
238 G4Trap::G4Trap(const G4Trap& rhs)                 240 G4Trap::G4Trap(const G4Trap& rhs)
239   : G4CSGSolid(rhs), halfCarTolerance(rhs.half    241   : G4CSGSolid(rhs), halfCarTolerance(rhs.halfCarTolerance),
240     fDz(rhs.fDz), fTthetaCphi(rhs.fTthetaCphi)    242     fDz(rhs.fDz), fTthetaCphi(rhs.fTthetaCphi), fTthetaSphi(rhs.fTthetaSphi),
241     fDy1(rhs.fDy1), fDx1(rhs.fDx1), fDx2(rhs.f    243     fDy1(rhs.fDy1), fDx1(rhs.fDx1), fDx2(rhs.fDx2), fTalpha1(rhs.fTalpha1),
242     fDy2(rhs.fDy2), fDx3(rhs.fDx3), fDx4(rhs.f    244     fDy2(rhs.fDy2), fDx3(rhs.fDx3), fDx4(rhs.fDx4), fTalpha2(rhs.fTalpha2)
243 {                                                 245 {
244   for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs    246   for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs.fPlanes[i]; }
245   for (G4int i=0; i<6; ++i) { fAreas[i] = rhs.    247   for (G4int i=0; i<6; ++i) { fAreas[i] = rhs.fAreas[i]; }
246   fTrapType = rhs.fTrapType;                      248   fTrapType = rhs.fTrapType;
247 }                                                 249 }
248                                                   250 
249 //////////////////////////////////////////////    251 //////////////////////////////////////////////////////////////////////////
250 //                                                252 //
251 // Assignment operator                            253 // Assignment operator
252                                                   254 
253 G4Trap& G4Trap::operator = (const G4Trap& rhs)    255 G4Trap& G4Trap::operator = (const G4Trap& rhs)
254 {                                                 256 {
255   // Check assignment to self                     257   // Check assignment to self
256   //                                              258   //
257   if (this == &rhs)  { return *this; }            259   if (this == &rhs)  { return *this; }
258                                                   260 
259   // Copy base class data                         261   // Copy base class data
260   //                                              262   //
261   G4CSGSolid::operator=(rhs);                     263   G4CSGSolid::operator=(rhs);
262                                                   264 
263   // Copy data                                    265   // Copy data
264   //                                              266   //
265   halfCarTolerance = rhs.halfCarTolerance;        267   halfCarTolerance = rhs.halfCarTolerance;
266   fDz = rhs.fDz; fTthetaCphi = rhs.fTthetaCphi    268   fDz = rhs.fDz; fTthetaCphi = rhs.fTthetaCphi; fTthetaSphi = rhs.fTthetaSphi;
267   fDy1 = rhs.fDy1; fDx1 = rhs.fDx1; fDx2 = rhs    269   fDy1 = rhs.fDy1; fDx1 = rhs.fDx1; fDx2 = rhs.fDx2; fTalpha1 = rhs.fTalpha1;
268   fDy2 = rhs.fDy2; fDx3 = rhs.fDx3; fDx4 = rhs    270   fDy2 = rhs.fDy2; fDx3 = rhs.fDx3; fDx4 = rhs.fDx4; fTalpha2 = rhs.fTalpha2;
269   for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs    271   for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs.fPlanes[i]; }
270   for (G4int i=0; i<6; ++i) { fAreas[i] = rhs.    272   for (G4int i=0; i<6; ++i) { fAreas[i] = rhs.fAreas[i]; }
271   fTrapType = rhs.fTrapType;                      273   fTrapType = rhs.fTrapType;
272   return *this;                                   274   return *this;
273 }                                                 275 }
274                                                   276 
275 //////////////////////////////////////////////    277 //////////////////////////////////////////////////////////////////////////
276 //                                                278 //
277 // Set all parameters, as for constructor - ch    279 // Set all parameters, as for constructor - check and set half-widths
278 // as well as angles: final check of coplanari    280 // as well as angles: final check of coplanarity
279                                                   281 
280 void G4Trap::SetAllParameters ( G4double pDz,     282 void G4Trap::SetAllParameters ( G4double pDz,
281                                 G4double pThet    283                                 G4double pTheta,
282                                 G4double pPhi,    284                                 G4double pPhi,
283                                 G4double pDy1,    285                                 G4double pDy1,
284                                 G4double pDx1,    286                                 G4double pDx1,
285                                 G4double pDx2,    287                                 G4double pDx2,
286                                 G4double pAlp1    288                                 G4double pAlp1,
287                                 G4double pDy2,    289                                 G4double pDy2,
288                                 G4double pDx3,    290                                 G4double pDx3,
289                                 G4double pDx4,    291                                 G4double pDx4,
290                                 G4double pAlp2    292                                 G4double pAlp2 )
291 {                                                 293 {
292   // Reset data of the base class                 294   // Reset data of the base class
293   fCubicVolume = 0;                               295   fCubicVolume = 0;
294   fSurfaceArea = 0;                               296   fSurfaceArea = 0;
295   fRebuildPolyhedron = true;                      297   fRebuildPolyhedron = true;
296                                                   298 
297   // Set parameters                               299   // Set parameters
298   fDz = pDz;                                      300   fDz = pDz;
299   fTthetaCphi = std::tan(pTheta)*std::cos(pPhi    301   fTthetaCphi = std::tan(pTheta)*std::cos(pPhi);
300   fTthetaSphi = std::tan(pTheta)*std::sin(pPhi    302   fTthetaSphi = std::tan(pTheta)*std::sin(pPhi);
301                                                   303 
302   fDy1 = pDy1; fDx1 = pDx1; fDx2 = pDx2; fTalp    304   fDy1 = pDy1; fDx1 = pDx1; fDx2 = pDx2; fTalpha1 = std::tan(pAlp1);
303   fDy2 = pDy2; fDx3 = pDx3; fDx4 = pDx4; fTalp    305   fDy2 = pDy2; fDx3 = pDx3; fDx4 = pDx4; fTalpha2 = std::tan(pAlp2);
304                                                   306 
305   CheckParameters();                              307   CheckParameters();
306   MakePlanes();                                   308   MakePlanes();
307 }                                                 309 }
308                                                   310 
309 //////////////////////////////////////////////    311 //////////////////////////////////////////////////////////////////////////
310 //                                                312 //
311 // Check length parameters                        313 // Check length parameters
312                                                   314 
313 void G4Trap::CheckParameters()                    315 void G4Trap::CheckParameters()
314 {                                                 316 {
315   if (fDz<=0 ||                                   317   if (fDz<=0 ||
316       fDy1<=0 || fDx1<=0 || fDx2<=0 ||            318       fDy1<=0 || fDx1<=0 || fDx2<=0 ||
317       fDy2<=0 || fDx3<=0 || fDx4<=0)              319       fDy2<=0 || fDx3<=0 || fDx4<=0)
318   {                                               320   {
319     std::ostringstream message;                   321     std::ostringstream message;
320     message << "Invalid Length Parameters for     322     message << "Invalid Length Parameters for Solid: " << GetName()
321             << "\n  X - " <<fDx1<<", "<<fDx2<<    323             << "\n  X - " <<fDx1<<", "<<fDx2<<", "<<fDx3<<", "<<fDx4
322             << "\n  Y - " <<fDy1<<", "<<fDy2      324             << "\n  Y - " <<fDy1<<", "<<fDy2
323             << "\n  Z - " <<fDz;                  325             << "\n  Z - " <<fDz;
324     G4Exception("G4Trap::CheckParameters()", "    326     G4Exception("G4Trap::CheckParameters()", "GeomSolids0002",
325                 FatalException, message);         327                 FatalException, message);
326   }                                               328   }
327 }                                                 329 }
328                                                   330 
329 //////////////////////////////////////////////    331 //////////////////////////////////////////////////////////////////////////
330 //                                                332 //
331 // Compute vertices and set side planes           333 // Compute vertices and set side planes
332                                                   334 
333 void G4Trap::MakePlanes()                         335 void G4Trap::MakePlanes()
334 {                                                 336 {
335   G4double DzTthetaCphi = fDz*fTthetaCphi;        337   G4double DzTthetaCphi = fDz*fTthetaCphi;
336   G4double DzTthetaSphi = fDz*fTthetaSphi;        338   G4double DzTthetaSphi = fDz*fTthetaSphi;
337   G4double Dy1Talpha1   = fDy1*fTalpha1;          339   G4double Dy1Talpha1   = fDy1*fTalpha1;
338   G4double Dy2Talpha2   = fDy2*fTalpha2;          340   G4double Dy2Talpha2   = fDy2*fTalpha2;
339                                                   341 
340   G4ThreeVector pt[8] =                           342   G4ThreeVector pt[8] =
341   {                                               343   {
342     G4ThreeVector(-DzTthetaCphi-Dy1Talpha1-fDx    344     G4ThreeVector(-DzTthetaCphi-Dy1Talpha1-fDx1,-DzTthetaSphi-fDy1,-fDz),
343     G4ThreeVector(-DzTthetaCphi-Dy1Talpha1+fDx    345     G4ThreeVector(-DzTthetaCphi-Dy1Talpha1+fDx1,-DzTthetaSphi-fDy1,-fDz),
344     G4ThreeVector(-DzTthetaCphi+Dy1Talpha1-fDx    346     G4ThreeVector(-DzTthetaCphi+Dy1Talpha1-fDx2,-DzTthetaSphi+fDy1,-fDz),
345     G4ThreeVector(-DzTthetaCphi+Dy1Talpha1+fDx    347     G4ThreeVector(-DzTthetaCphi+Dy1Talpha1+fDx2,-DzTthetaSphi+fDy1,-fDz),
346     G4ThreeVector( DzTthetaCphi-Dy2Talpha2-fDx    348     G4ThreeVector( DzTthetaCphi-Dy2Talpha2-fDx3, DzTthetaSphi-fDy2, fDz),
347     G4ThreeVector( DzTthetaCphi-Dy2Talpha2+fDx    349     G4ThreeVector( DzTthetaCphi-Dy2Talpha2+fDx3, DzTthetaSphi-fDy2, fDz),
348     G4ThreeVector( DzTthetaCphi+Dy2Talpha2-fDx    350     G4ThreeVector( DzTthetaCphi+Dy2Talpha2-fDx4, DzTthetaSphi+fDy2, fDz),
349     G4ThreeVector( DzTthetaCphi+Dy2Talpha2+fDx    351     G4ThreeVector( DzTthetaCphi+Dy2Talpha2+fDx4, DzTthetaSphi+fDy2, fDz)
350   };                                              352   };
351                                                   353 
352   MakePlanes(pt);                                 354   MakePlanes(pt);
353 }                                                 355 }
354                                                   356 
355 //////////////////////////////////////////////    357 //////////////////////////////////////////////////////////////////////////
356 //                                                358 //
357 // Set side planes, check planarity               359 // Set side planes, check planarity
358                                                   360 
359 void G4Trap::MakePlanes(const G4ThreeVector pt    361 void G4Trap::MakePlanes(const G4ThreeVector pt[8])
360 {                                                 362 {
361   constexpr G4int iface[4][4] = { {0,4,5,1}, {    363   constexpr G4int iface[4][4] = { {0,4,5,1}, {2,3,7,6}, {0,2,6,4}, {1,5,7,3} };
362   const static G4String side[4] = { "~-Y", "~+    364   const static G4String side[4] = { "~-Y", "~+Y", "~-X", "~+X" };
363                                                   365 
364   for (G4int i=0; i<4; ++i)                       366   for (G4int i=0; i<4; ++i)
365   {                                               367   {
366     if (MakePlane(pt[iface[i][0]],                368     if (MakePlane(pt[iface[i][0]],
367                   pt[iface[i][1]],                369                   pt[iface[i][1]],
368                   pt[iface[i][2]],                370                   pt[iface[i][2]],
369                   pt[iface[i][3]],                371                   pt[iface[i][3]],
370                   fPlanes[i])) continue;          372                   fPlanes[i])) continue;
371                                                   373 
372     // Non planar side face                       374     // Non planar side face
373     G4ThreeVector normal(fPlanes[i].a,fPlanes[    375     G4ThreeVector normal(fPlanes[i].a,fPlanes[i].b,fPlanes[i].c);
374     G4double dmax = 0;                            376     G4double dmax = 0;
375     for (G4int k=0; k<4; ++k)                     377     for (G4int k=0; k<4; ++k)
376     {                                             378     {
377       G4double dist = normal.dot(pt[iface[i][k    379       G4double dist = normal.dot(pt[iface[i][k]]) + fPlanes[i].d;
378       if (std::abs(dist) > std::abs(dmax)) dma    380       if (std::abs(dist) > std::abs(dmax)) dmax = dist;
379     }                                             381     }
380     std::ostringstream message;                   382     std::ostringstream message;
381     message << "Side face " << side[i] << " is    383     message << "Side face " << side[i] << " is not planar for solid: "
382             << GetName() << "\nDiscrepancy: "     384             << GetName() << "\nDiscrepancy: " << dmax/mm << " mm\n";
383     StreamInfo(message);                          385     StreamInfo(message);
384     G4Exception("G4Trap::MakePlanes()", "GeomS    386     G4Exception("G4Trap::MakePlanes()", "GeomSolids0002",
385                 FatalException, message);         387                 FatalException, message);
386   }                                               388   }
387                                                   389 
388   // Re-compute parameters                        390   // Re-compute parameters
389   SetCachedValues();                              391   SetCachedValues();
390 }                                                 392 }
391                                                   393 
392 //////////////////////////////////////////////    394 //////////////////////////////////////////////////////////////////////////
393 //                                                395 //
394 // Calculate the coef's of the plane p1->p2->p    396 // Calculate the coef's of the plane p1->p2->p3->p4->p1
395 // where the ThreeVectors 1-4 are in anti-cloc    397 // where the ThreeVectors 1-4 are in anti-clockwise order when viewed
396 // from infront of the plane (i.e. from normal    398 // from infront of the plane (i.e. from normal direction).
397 //                                                399 //
398 // Return true if the points are coplanar, fal    400 // Return true if the points are coplanar, false otherwise
399                                                   401 
400 G4bool G4Trap::MakePlane( const G4ThreeVector&    402 G4bool G4Trap::MakePlane( const G4ThreeVector& p1,
401                           const G4ThreeVector&    403                           const G4ThreeVector& p2,
402                           const G4ThreeVector&    404                           const G4ThreeVector& p3,
403                           const G4ThreeVector&    405                           const G4ThreeVector& p4,
404                                 TrapSidePlane&    406                                 TrapSidePlane& plane )
405 {                                                 407 {
406   G4ThreeVector normal = ((p4 - p2).cross(p3 -    408   G4ThreeVector normal = ((p4 - p2).cross(p3 - p1)).unit();
407   if (std::abs(normal.x()) < DBL_EPSILON) norm    409   if (std::abs(normal.x()) < DBL_EPSILON) normal.setX(0);
408   if (std::abs(normal.y()) < DBL_EPSILON) norm    410   if (std::abs(normal.y()) < DBL_EPSILON) normal.setY(0);
409   if (std::abs(normal.z()) < DBL_EPSILON) norm    411   if (std::abs(normal.z()) < DBL_EPSILON) normal.setZ(0);
410   normal = normal.unit();                         412   normal = normal.unit();
411                                                   413 
412   G4ThreeVector centre = (p1 + p2 + p3 + p4)*0    414   G4ThreeVector centre = (p1 + p2 + p3 + p4)*0.25;
413   plane.a =  normal.x();                          415   plane.a =  normal.x();
414   plane.b =  normal.y();                          416   plane.b =  normal.y();
415   plane.c =  normal.z();                          417   plane.c =  normal.z();
416   plane.d = -normal.dot(centre);                  418   plane.d = -normal.dot(centre);
417                                                   419 
418   // compute distances and check planarity        420   // compute distances and check planarity
419   G4double d1 = std::abs(normal.dot(p1) + plan    421   G4double d1 = std::abs(normal.dot(p1) + plane.d);
420   G4double d2 = std::abs(normal.dot(p2) + plan    422   G4double d2 = std::abs(normal.dot(p2) + plane.d);
421   G4double d3 = std::abs(normal.dot(p3) + plan    423   G4double d3 = std::abs(normal.dot(p3) + plane.d);
422   G4double d4 = std::abs(normal.dot(p4) + plan    424   G4double d4 = std::abs(normal.dot(p4) + plane.d);
423   G4double dmax = std::max(std::max(std::max(d    425   G4double dmax = std::max(std::max(std::max(d1,d2),d3),d4);
424                                                   426 
425   return dmax <= 1000 * kCarTolerance;         << 427   return (dmax > 1000 * kCarTolerance) ? false : true;
426 }                                                 428 }
427                                                   429 
428 //////////////////////////////////////////////    430 //////////////////////////////////////////////////////////////////////////
429 //                                                431 //
430 // Recompute parameters using planes              432 // Recompute parameters using planes
431                                                   433 
432 void G4Trap::SetCachedValues()                    434 void G4Trap::SetCachedValues()
433 {                                                 435 {
434   // Set indeces                                  436   // Set indeces
435   constexpr  G4int iface[6][4] =                  437   constexpr  G4int iface[6][4] =
436       { {0,1,3,2}, {0,4,5,1}, {2,3,7,6}, {0,2,    438       { {0,1,3,2}, {0,4,5,1}, {2,3,7,6}, {0,2,6,4}, {1,5,7,3}, {4,6,7,5} };
437                                                   439 
438   // Get vertices                                 440   // Get vertices
439   G4ThreeVector pt[8];                            441   G4ThreeVector pt[8];
440   GetVertices(pt);                                442   GetVertices(pt);
441                                                   443 
442   // Set face areas                               444   // Set face areas
443   for (G4int i=0; i<6; ++i)                       445   for (G4int i=0; i<6; ++i)
444   {                                               446   {
445     fAreas[i] = G4GeomTools::QuadAreaNormal(pt    447     fAreas[i] = G4GeomTools::QuadAreaNormal(pt[iface[i][0]],
446                                             pt    448                                             pt[iface[i][1]],
447                                             pt    449                                             pt[iface[i][2]],
448                                             pt    450                                             pt[iface[i][3]]).mag();
449   }                                               451   }
450   for (G4int i=1; i<6; ++i) { fAreas[i] += fAr    452   for (G4int i=1; i<6; ++i) { fAreas[i] += fAreas[i - 1]; }
451                                                   453 
452   // Define type of trapezoid                     454   // Define type of trapezoid
453   fTrapType = 0;                                  455   fTrapType = 0;
454   if (fPlanes[0].b  == -1 && fPlanes[1].b == 1    456   if (fPlanes[0].b  == -1 && fPlanes[1].b == 1 &&
455       std::abs(fPlanes[0].a) < DBL_EPSILON &&     457       std::abs(fPlanes[0].a) < DBL_EPSILON &&
456       std::abs(fPlanes[0].c) < DBL_EPSILON &&     458       std::abs(fPlanes[0].c) < DBL_EPSILON &&
457       std::abs(fPlanes[1].a) < DBL_EPSILON &&     459       std::abs(fPlanes[1].a) < DBL_EPSILON &&
458       std::abs(fPlanes[1].c) < DBL_EPSILON)       460       std::abs(fPlanes[1].c) < DBL_EPSILON)
459   {                                               461   {
460     fTrapType = 1; // YZ section is a rectangl    462     fTrapType = 1; // YZ section is a rectangle ...
461     if (std::abs(fPlanes[2].a + fPlanes[3].a)     463     if (std::abs(fPlanes[2].a + fPlanes[3].a) < DBL_EPSILON &&
462         std::abs(fPlanes[2].c - fPlanes[3].c)     464         std::abs(fPlanes[2].c - fPlanes[3].c) < DBL_EPSILON &&
463         fPlanes[2].b == 0 &&                      465         fPlanes[2].b == 0 &&
464         fPlanes[3].b == 0)                        466         fPlanes[3].b == 0)
465     {                                             467     {
466       fTrapType = 2; // ... and XZ section is     468       fTrapType = 2; // ... and XZ section is a isosceles trapezoid
467       fPlanes[2].a = -fPlanes[3].a;               469       fPlanes[2].a = -fPlanes[3].a;
468       fPlanes[2].c =  fPlanes[3].c;               470       fPlanes[2].c =  fPlanes[3].c;
469     }                                             471     }
470     if (std::abs(fPlanes[2].a + fPlanes[3].a)     472     if (std::abs(fPlanes[2].a + fPlanes[3].a) < DBL_EPSILON &&
471         std::abs(fPlanes[2].b - fPlanes[3].b)     473         std::abs(fPlanes[2].b - fPlanes[3].b) < DBL_EPSILON &&
472         fPlanes[2].c == 0 &&                      474         fPlanes[2].c == 0 &&
473         fPlanes[3].c == 0)                        475         fPlanes[3].c == 0)
474     {                                             476     {
475       fTrapType = 3; // ... and XY section is     477       fTrapType = 3; // ... and XY section is a isosceles trapezoid
476       fPlanes[2].a = -fPlanes[3].a;               478       fPlanes[2].a = -fPlanes[3].a;
477       fPlanes[2].b =  fPlanes[3].b;               479       fPlanes[2].b =  fPlanes[3].b;
478     }                                             480     }
479   }                                               481   }
480 }                                                 482 }
481                                                   483 
482 //////////////////////////////////////////////    484 //////////////////////////////////////////////////////////////////////////
483 //                                                485 //
484 // Get volume                                     486 // Get volume
485                                                   487 
486 G4double G4Trap::GetCubicVolume()                 488 G4double G4Trap::GetCubicVolume()
487 {                                                 489 {
488   if (fCubicVolume == 0)                          490   if (fCubicVolume == 0)
489   {                                               491   {
490     G4ThreeVector pt[8];                          492     G4ThreeVector pt[8];
491     GetVertices(pt);                              493     GetVertices(pt);
492                                                   494 
493     G4double dz  = pt[4].z() - pt[0].z();         495     G4double dz  = pt[4].z() - pt[0].z();
494     G4double dy1 = pt[2].y() - pt[0].y();         496     G4double dy1 = pt[2].y() - pt[0].y();
495     G4double dx1 = pt[1].x() - pt[0].x();         497     G4double dx1 = pt[1].x() - pt[0].x();
496     G4double dx2 = pt[3].x() - pt[2].x();         498     G4double dx2 = pt[3].x() - pt[2].x();
497     G4double dy2 = pt[6].y() - pt[4].y();         499     G4double dy2 = pt[6].y() - pt[4].y();
498     G4double dx3 = pt[5].x() - pt[4].x();         500     G4double dx3 = pt[5].x() - pt[4].x();
499     G4double dx4 = pt[7].x() - pt[6].x();         501     G4double dx4 = pt[7].x() - pt[6].x();
500                                                   502 
501     fCubicVolume = ((dx1 + dx2 + dx3 + dx4)*(d    503     fCubicVolume = ((dx1 + dx2 + dx3 + dx4)*(dy1 + dy2) +
502                     (dx4 + dx3 - dx2 - dx1)*(d    504                     (dx4 + dx3 - dx2 - dx1)*(dy2 - dy1)/3)*dz*0.125;
503   }                                               505   }
504   return fCubicVolume;                            506   return fCubicVolume;
505 }                                                 507 }
506                                                   508 
507 //////////////////////////////////////////////    509 //////////////////////////////////////////////////////////////////////////
508 //                                                510 //
509 // Get surface area                               511 // Get surface area
510                                                   512 
511 G4double G4Trap::GetSurfaceArea()                 513 G4double G4Trap::GetSurfaceArea()
512 {                                                 514 {
513   if (fSurfaceArea == 0)                          515   if (fSurfaceArea == 0)
514   {                                               516   {
515     G4ThreeVector pt[8];                          517     G4ThreeVector pt[8];
516     G4int iface [6][4] =                          518     G4int iface [6][4] =
517       { {0,1,3,2}, {0,4,5,1}, {2,3,7,6}, {0,2,    519       { {0,1,3,2}, {0,4,5,1}, {2,3,7,6}, {0,2,6,4}, {1,5,7,3}, {4,6,7,5} };
518                                                   520 
519     GetVertices(pt);                              521     GetVertices(pt);
520     for (const auto & i : iface)               << 522     for (G4int i=0; i<6; ++i)
521     {                                             523     {
522       fSurfaceArea += G4GeomTools::QuadAreaNor << 524       fSurfaceArea += G4GeomTools::QuadAreaNormal(pt[iface[i][0]],
523                                                << 525                                                   pt[iface[i][1]],
524                                                << 526                                                   pt[iface[i][2]],
525                                                << 527                                                   pt[iface[i][3]]).mag();
526     }                                             528     }
527   }                                               529   }
528   return fSurfaceArea;                            530   return fSurfaceArea;
529 }                                                 531 }
530                                                   532 
531 //////////////////////////////////////////////    533 //////////////////////////////////////////////////////////////////////////
532 //                                                534 //
533 // Dispatch to parameterisation for replicatio    535 // Dispatch to parameterisation for replication mechanism dimension
534 // computation & modification.                    536 // computation & modification.
535                                                   537 
536 void G4Trap::ComputeDimensions(       G4VPVPar    538 void G4Trap::ComputeDimensions(       G4VPVParameterisation* p,
537                                 const G4int n,    539                                 const G4int n,
538                                 const G4VPhysi    540                                 const G4VPhysicalVolume* pRep )
539 {                                                 541 {
540   p->ComputeDimensions(*this,n,pRep);             542   p->ComputeDimensions(*this,n,pRep);
541 }                                                 543 }
542                                                   544 
543 //////////////////////////////////////////////    545 //////////////////////////////////////////////////////////////////////////
544 //                                                546 //
545 // Get bounding box                               547 // Get bounding box
546                                                   548 
547 void G4Trap::BoundingLimits(G4ThreeVector& pMi    549 void G4Trap::BoundingLimits(G4ThreeVector& pMin, G4ThreeVector& pMax) const
548 {                                                 550 {
549   G4ThreeVector pt[8];                            551   G4ThreeVector pt[8];
550   GetVertices(pt);                                552   GetVertices(pt);
551                                                   553 
552   G4double xmin = kInfinity, xmax = -kInfinity    554   G4double xmin = kInfinity, xmax = -kInfinity;
553   G4double ymin = kInfinity, ymax = -kInfinity    555   G4double ymin = kInfinity, ymax = -kInfinity;
554   for (const auto & i : pt)                    << 556   for (G4int i=0; i<8; ++i)
555   {                                               557   {
556     G4double x = i.x();                        << 558     G4double x = pt[i].x();
557     if (x < xmin) xmin = x;                       559     if (x < xmin) xmin = x;
558     if (x > xmax) xmax = x;                       560     if (x > xmax) xmax = x;
559     G4double y = i.y();                        << 561     G4double y = pt[i].y();
560     if (y < ymin) ymin = y;                       562     if (y < ymin) ymin = y;
561     if (y > ymax) ymax = y;                       563     if (y > ymax) ymax = y;
562   }                                               564   }
563                                                   565 
564   G4double dz   = GetZHalfLength();               566   G4double dz   = GetZHalfLength();
565   pMin.set(xmin,ymin,-dz);                        567   pMin.set(xmin,ymin,-dz);
566   pMax.set(xmax,ymax, dz);                        568   pMax.set(xmax,ymax, dz);
567                                                   569 
568   // Check correctness of the bounding box        570   // Check correctness of the bounding box
569   //                                              571   //
570   if (pMin.x() >= pMax.x() || pMin.y() >= pMax    572   if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
571   {                                               573   {
572     std::ostringstream message;                   574     std::ostringstream message;
573     message << "Bad bounding box (min >= max)     575     message << "Bad bounding box (min >= max) for solid: "
574             << GetName() << " !"                  576             << GetName() << " !"
575             << "\npMin = " << pMin                577             << "\npMin = " << pMin
576             << "\npMax = " << pMax;               578             << "\npMax = " << pMax;
577     G4Exception("G4Trap::BoundingLimits()", "G    579     G4Exception("G4Trap::BoundingLimits()", "GeomMgt0001",
578                 JustWarning, message);            580                 JustWarning, message);
579     DumpInfo();                                   581     DumpInfo();
580   }                                               582   }
581 }                                                 583 }
582                                                   584 
583 //////////////////////////////////////////////    585 //////////////////////////////////////////////////////////////////////////
584 //                                                586 //
585 // Calculate extent under transform and specif    587 // Calculate extent under transform and specified limit
586                                                   588 
587 G4bool G4Trap::CalculateExtent( const EAxis pA    589 G4bool G4Trap::CalculateExtent( const EAxis pAxis,
588                                 const G4VoxelL    590                                 const G4VoxelLimits& pVoxelLimit,
589                                 const G4Affine    591                                 const G4AffineTransform& pTransform,
590                                       G4double    592                                       G4double& pMin, G4double& pMax) const
591 {                                                 593 {
592   G4ThreeVector bmin, bmax;                       594   G4ThreeVector bmin, bmax;
593   G4bool exist;                                   595   G4bool exist;
594                                                   596 
595   // Check bounding box (bbox)                    597   // Check bounding box (bbox)
596   //                                              598   //
597   BoundingLimits(bmin,bmax);                      599   BoundingLimits(bmin,bmax);
598   G4BoundingEnvelope bbox(bmin,bmax);             600   G4BoundingEnvelope bbox(bmin,bmax);
599 #ifdef G4BBOX_EXTENT                              601 #ifdef G4BBOX_EXTENT
600   return bbox.CalculateExtent(pAxis,pVoxelLimi    602   return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
601 #endif                                            603 #endif
602   if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox    604   if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
603   {                                               605   {
604     return exist = pMin < pMax;                << 606     return exist = (pMin < pMax) ? true : false;
605   }                                               607   }
606                                                   608 
607   // Set bounding envelope (benv) and calculat    609   // Set bounding envelope (benv) and calculate extent
608   //                                              610   //
609   G4ThreeVector pt[8];                            611   G4ThreeVector pt[8];
610   GetVertices(pt);                                612   GetVertices(pt);
611                                                   613 
612   G4ThreeVectorList baseA(4), baseB(4);           614   G4ThreeVectorList baseA(4), baseB(4);
613   baseA[0] = pt[0];                               615   baseA[0] = pt[0];
614   baseA[1] = pt[1];                               616   baseA[1] = pt[1];
615   baseA[2] = pt[3];                               617   baseA[2] = pt[3];
616   baseA[3] = pt[2];                               618   baseA[3] = pt[2];
617                                                   619 
618   baseB[0] = pt[4];                               620   baseB[0] = pt[4];
619   baseB[1] = pt[5];                               621   baseB[1] = pt[5];
620   baseB[2] = pt[7];                               622   baseB[2] = pt[7];
621   baseB[3] = pt[6];                               623   baseB[3] = pt[6];
622                                                   624 
623   std::vector<const G4ThreeVectorList *> polyg    625   std::vector<const G4ThreeVectorList *> polygons(2);
624   polygons[0] = &baseA;                           626   polygons[0] = &baseA;
625   polygons[1] = &baseB;                           627   polygons[1] = &baseB;
626                                                   628 
627   G4BoundingEnvelope benv(bmin,bmax,polygons);    629   G4BoundingEnvelope benv(bmin,bmax,polygons);
628   exist = benv.CalculateExtent(pAxis,pVoxelLim    630   exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
629   return exist;                                   631   return exist;
630 }                                                 632 }
631                                                   633 
632 //////////////////////////////////////////////    634 //////////////////////////////////////////////////////////////////////////
633 //                                                635 //
634 // Return whether point is inside/outside/on_s    636 // Return whether point is inside/outside/on_surface
635                                                   637 
636 EInside G4Trap::Inside( const G4ThreeVector& p    638 EInside G4Trap::Inside( const G4ThreeVector& p ) const
637 {                                                 639 {
638   switch (fTrapType)                              640   switch (fTrapType)
639   {                                               641   {
640     case 0: // General case                       642     case 0: // General case
641     {                                             643     {
642       G4double dz = std::abs(p.z())-fDz;          644       G4double dz = std::abs(p.z())-fDz;
643       G4double dy1 = fPlanes[0].b*p.y()+fPlane    645       G4double dy1 = fPlanes[0].b*p.y()+fPlanes[0].c*p.z()+fPlanes[0].d;
644       G4double dy2 = fPlanes[1].b*p.y()+fPlane    646       G4double dy2 = fPlanes[1].b*p.y()+fPlanes[1].c*p.z()+fPlanes[1].d;
645       G4double dy = std::max(dz,std::max(dy1,d    647       G4double dy = std::max(dz,std::max(dy1,dy2));
646                                                   648 
647       G4double dx1 = fPlanes[2].a*p.x()+fPlane    649       G4double dx1 = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()
648                    + fPlanes[2].c*p.z()+fPlane    650                    + fPlanes[2].c*p.z()+fPlanes[2].d;
649       G4double dx2 = fPlanes[3].a*p.x()+fPlane    651       G4double dx2 = fPlanes[3].a*p.x()+fPlanes[3].b*p.y()
650                    + fPlanes[3].c*p.z()+fPlane    652                    + fPlanes[3].c*p.z()+fPlanes[3].d;
651       G4double dist = std::max(dy,std::max(dx1    653       G4double dist = std::max(dy,std::max(dx1,dx2));
652                                                   654 
653       return (dist > halfCarTolerance) ? kOuts    655       return (dist > halfCarTolerance) ? kOutside :
654         ((dist > -halfCarTolerance) ? kSurface    656         ((dist > -halfCarTolerance) ? kSurface : kInside);
655     }                                             657     }
656     case 1: // YZ section is a rectangle          658     case 1: // YZ section is a rectangle
657     {                                             659     {
658       G4double dz = std::abs(p.z())-fDz;          660       G4double dz = std::abs(p.z())-fDz;
659       G4double dy = std::max(dz,std::abs(p.y()    661       G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
660       G4double dx1 = fPlanes[2].a*p.x()+fPlane    662       G4double dx1 = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()
661                    + fPlanes[2].c*p.z()+fPlane    663                    + fPlanes[2].c*p.z()+fPlanes[2].d;
662       G4double dx2 = fPlanes[3].a*p.x()+fPlane    664       G4double dx2 = fPlanes[3].a*p.x()+fPlanes[3].b*p.y()
663                    + fPlanes[3].c*p.z()+fPlane    665                    + fPlanes[3].c*p.z()+fPlanes[3].d;
664       G4double dist = std::max(dy,std::max(dx1    666       G4double dist = std::max(dy,std::max(dx1,dx2));
665                                                   667 
666       return (dist > halfCarTolerance) ? kOuts    668       return (dist > halfCarTolerance) ? kOutside :
667         ((dist > -halfCarTolerance) ? kSurface    669         ((dist > -halfCarTolerance) ? kSurface : kInside);
668     }                                             670     }
669     case 2: // YZ section is a rectangle and      671     case 2: // YZ section is a rectangle and
670     {       // XZ section is an isosceles trap    672     {       // XZ section is an isosceles trapezoid
671       G4double dz = std::abs(p.z())-fDz;          673       G4double dz = std::abs(p.z())-fDz;
672       G4double dy = std::max(dz,std::abs(p.y()    674       G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
673       G4double dx = fPlanes[3].a*std::abs(p.x(    675       G4double dx = fPlanes[3].a*std::abs(p.x())
674                   + fPlanes[3].c*p.z()+fPlanes    676                   + fPlanes[3].c*p.z()+fPlanes[3].d;
675       G4double dist = std::max(dy,dx);            677       G4double dist = std::max(dy,dx);
676                                                   678 
677       return (dist > halfCarTolerance) ? kOuts    679       return (dist > halfCarTolerance) ? kOutside :
678         ((dist > -halfCarTolerance) ? kSurface    680         ((dist > -halfCarTolerance) ? kSurface : kInside);
679     }                                             681     }
680     case 3: // YZ section is a rectangle and      682     case 3: // YZ section is a rectangle and
681     {       // XY section is an isosceles trap    683     {       // XY section is an isosceles trapezoid
682       G4double dz = std::abs(p.z())-fDz;          684       G4double dz = std::abs(p.z())-fDz;
683       G4double dy = std::max(dz,std::abs(p.y()    685       G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
684       G4double dx = fPlanes[3].a*std::abs(p.x(    686       G4double dx = fPlanes[3].a*std::abs(p.x())
685                   + fPlanes[3].b*p.y()+fPlanes    687                   + fPlanes[3].b*p.y()+fPlanes[3].d;
686       G4double dist = std::max(dy,dx);            688       G4double dist = std::max(dy,dx);
687                                                   689 
688       return (dist > halfCarTolerance) ? kOuts    690       return (dist > halfCarTolerance) ? kOutside :
689         ((dist > -halfCarTolerance) ? kSurface    691         ((dist > -halfCarTolerance) ? kSurface : kInside);
690     }                                             692     }
691   }                                               693   }
692   return kOutside;                                694   return kOutside;
693 }                                                 695 }
694                                                   696 
695 //////////////////////////////////////////////    697 //////////////////////////////////////////////////////////////////////////
696 //                                                698 //
697 // Determine side, and return corresponding no    699 // Determine side, and return corresponding normal
698                                                   700 
699 G4ThreeVector G4Trap::SurfaceNormal( const G4T    701 G4ThreeVector G4Trap::SurfaceNormal( const G4ThreeVector& p ) const
700 {                                                 702 {
701   G4double nx = 0, ny = 0, nz = 0;                703   G4double nx = 0, ny = 0, nz = 0;
702   G4double dz = std::abs(p.z()) - fDz;            704   G4double dz = std::abs(p.z()) - fDz;
703   nz = std::copysign(G4double(std::abs(dz) <=     705   nz = std::copysign(G4double(std::abs(dz) <= halfCarTolerance), p.z());
704                                                   706 
705   switch (fTrapType)                              707   switch (fTrapType)
706   {                                               708   {
707     case 0: // General case                       709     case 0: // General case
708     {                                             710     {
709       for (G4int i=0; i<2; ++i)                   711       for (G4int i=0; i<2; ++i)
710       {                                           712       {
711         G4double dy = fPlanes[i].b*p.y() + fPl    713         G4double dy = fPlanes[i].b*p.y() + fPlanes[i].c*p.z() + fPlanes[i].d;
712         if (std::abs(dy) > halfCarTolerance) c    714         if (std::abs(dy) > halfCarTolerance) continue;
713         ny  = fPlanes[i].b;                       715         ny  = fPlanes[i].b;
714         nz += fPlanes[i].c;                       716         nz += fPlanes[i].c;
715         break;                                    717         break;
716       }                                           718       }
717       for (G4int i=2; i<4; ++i)                   719       for (G4int i=2; i<4; ++i)
718       {                                           720       {
719         G4double dx = fPlanes[i].a*p.x() +        721         G4double dx = fPlanes[i].a*p.x() +
720                       fPlanes[i].b*p.y() + fPl    722                       fPlanes[i].b*p.y() + fPlanes[i].c*p.z() + fPlanes[i].d;
721         if (std::abs(dx) > halfCarTolerance) c    723         if (std::abs(dx) > halfCarTolerance) continue;
722         nx  = fPlanes[i].a;                       724         nx  = fPlanes[i].a;
723         ny += fPlanes[i].b;                       725         ny += fPlanes[i].b;
724         nz += fPlanes[i].c;                       726         nz += fPlanes[i].c;
725         break;                                    727         break;
726       }                                           728       }
727       break;                                      729       break;
728     }                                             730     }
729     case 1: // YZ section - rectangle             731     case 1: // YZ section - rectangle
730     {                                             732     {
731       G4double dy = std::abs(p.y()) + fPlanes[    733       G4double dy = std::abs(p.y()) + fPlanes[1].d;
732       ny = std::copysign(G4double(std::abs(dy)    734       ny = std::copysign(G4double(std::abs(dy) <= halfCarTolerance), p.y());
733       for (G4int i=2; i<4; ++i)                   735       for (G4int i=2; i<4; ++i)
734       {                                           736       {
735         G4double dx = fPlanes[i].a*p.x() +        737         G4double dx = fPlanes[i].a*p.x() +
736                       fPlanes[i].b*p.y() + fPl    738                       fPlanes[i].b*p.y() + fPlanes[i].c*p.z() + fPlanes[i].d;
737         if (std::abs(dx) > halfCarTolerance) c    739         if (std::abs(dx) > halfCarTolerance) continue;
738         nx  = fPlanes[i].a;                       740         nx  = fPlanes[i].a;
739         ny += fPlanes[i].b;                       741         ny += fPlanes[i].b;
740         nz += fPlanes[i].c;                       742         nz += fPlanes[i].c;
741         break;                                    743         break;
742       }                                           744       }
743       break;                                      745       break;
744     }                                             746     }
745     case 2: // YZ section - rectangle, XZ sect    747     case 2: // YZ section - rectangle, XZ section - isosceles trapezoid
746     {                                             748     {
747       G4double dy = std::abs(p.y()) + fPlanes[    749       G4double dy = std::abs(p.y()) + fPlanes[1].d;
748       ny = std::copysign(G4double(std::abs(dy)    750       ny = std::copysign(G4double(std::abs(dy) <= halfCarTolerance), p.y());
749       G4double dx = fPlanes[3].a*std::abs(p.x(    751       G4double dx = fPlanes[3].a*std::abs(p.x()) +
750                     fPlanes[3].c*p.z() + fPlan    752                     fPlanes[3].c*p.z() + fPlanes[3].d;
751       G4double k = std::abs(dx) <= halfCarTole    753       G4double k = std::abs(dx) <= halfCarTolerance;
752       nx  = std::copysign(k, p.x())*fPlanes[3]    754       nx  = std::copysign(k, p.x())*fPlanes[3].a;
753       nz += k*fPlanes[3].c;                       755       nz += k*fPlanes[3].c;
754       break;                                      756       break;
755     }                                             757     }
756     case 3: // YZ section - rectangle, XY sect    758     case 3: // YZ section - rectangle, XY section - isosceles trapezoid
757     {                                             759     {
758       G4double dy = std::abs(p.y()) + fPlanes[    760       G4double dy = std::abs(p.y()) + fPlanes[1].d;
759       ny = std::copysign(G4double(std::abs(dy)    761       ny = std::copysign(G4double(std::abs(dy) <= halfCarTolerance), p.y());
760       G4double dx = fPlanes[3].a*std::abs(p.x(    762       G4double dx = fPlanes[3].a*std::abs(p.x()) +
761                     fPlanes[3].b*p.y() + fPlan    763                     fPlanes[3].b*p.y() + fPlanes[3].d;
762       G4double k = std::abs(dx) <= halfCarTole    764       G4double k = std::abs(dx) <= halfCarTolerance;
763       nx  = std::copysign(k, p.x())*fPlanes[3]    765       nx  = std::copysign(k, p.x())*fPlanes[3].a;
764       ny += k*fPlanes[3].b;                       766       ny += k*fPlanes[3].b;
765       break;                                      767       break;
766     }                                             768     }
767   }                                               769   }
768                                                   770 
769   // Return normal                                771   // Return normal
770   //                                              772   //
771   G4double mag2 = nx*nx + ny*ny + nz*nz;          773   G4double mag2 = nx*nx + ny*ny + nz*nz;
772   if (mag2 == 1)      return { nx,ny,nz };     << 774   if (mag2 == 1)      return G4ThreeVector(nx,ny,nz);
773   else if (mag2 != 0) return G4ThreeVector(nx,    775   else if (mag2 != 0) return G4ThreeVector(nx,ny,nz).unit(); // edge or corner
774   else                                            776   else
775   {                                               777   {
776     // Point is not on the surface                778     // Point is not on the surface
777     //                                            779     //
778 #ifdef G4CSGDEBUG                                 780 #ifdef G4CSGDEBUG
779     std::ostringstream message;                   781     std::ostringstream message;
780     G4long oldprc = message.precision(16);     << 782     G4int oldprc = message.precision(16);
781     message << "Point p is not on surface (!?)    783     message << "Point p is not on surface (!?) of solid: "
782             << GetName() << G4endl;               784             << GetName() << G4endl;
783     message << "Position:\n";                     785     message << "Position:\n";
784     message << "   p.x() = " << p.x()/mm << "     786     message << "   p.x() = " << p.x()/mm << " mm\n";
785     message << "   p.y() = " << p.y()/mm << "     787     message << "   p.y() = " << p.y()/mm << " mm\n";
786     message << "   p.z() = " << p.z()/mm << "     788     message << "   p.z() = " << p.z()/mm << " mm";
787     G4cout.precision(oldprc) ;                    789     G4cout.precision(oldprc) ;
788     G4Exception("G4Trap::SurfaceNormal(p)", "G    790     G4Exception("G4Trap::SurfaceNormal(p)", "GeomSolids1002",
789                 JustWarning, message );           791                 JustWarning, message );
790     DumpInfo();                                   792     DumpInfo();
791 #endif                                            793 #endif
792     return ApproxSurfaceNormal(p);                794     return ApproxSurfaceNormal(p);
793   }                                               795   }
794 }                                                 796 }
795                                                   797 
796 //////////////////////////////////////////////    798 //////////////////////////////////////////////////////////////////////////
797 //                                                799 //
798 // Algorithm for SurfaceNormal() following the    800 // Algorithm for SurfaceNormal() following the original specification
799 // for points not on the surface                  801 // for points not on the surface
800                                                   802 
801 G4ThreeVector G4Trap::ApproxSurfaceNormal( con    803 G4ThreeVector G4Trap::ApproxSurfaceNormal( const G4ThreeVector& p ) const
802 {                                                 804 {
803   G4double dist = -DBL_MAX;                       805   G4double dist = -DBL_MAX;
804   G4int iside = 0;                                806   G4int iside = 0;
805   for (G4int i=0; i<4; ++i)                       807   for (G4int i=0; i<4; ++i)
806   {                                               808   {
807     G4double d = fPlanes[i].a*p.x() +             809     G4double d = fPlanes[i].a*p.x() +
808                  fPlanes[i].b*p.y() +             810                  fPlanes[i].b*p.y() +
809                  fPlanes[i].c*p.z() + fPlanes[    811                  fPlanes[i].c*p.z() + fPlanes[i].d;
810     if (d > dist) { dist = d; iside = i; }        812     if (d > dist) { dist = d; iside = i; }
811   }                                               813   }
812                                                   814 
813   G4double distz = std::abs(p.z()) - fDz;         815   G4double distz = std::abs(p.z()) - fDz;
814   if (dist > distz)                               816   if (dist > distz)
815     return { fPlanes[iside].a, fPlanes[iside]. << 817     return G4ThreeVector(fPlanes[iside].a, fPlanes[iside].b, fPlanes[iside].c);
816   else                                            818   else
817     return { 0, 0, (G4double)((p.z() < 0) ? -1 << 819     return G4ThreeVector(0, 0, (p.z() < 0) ? -1 : 1);
818 }                                                 820 }
819                                                   821 
820 //////////////////////////////////////////////    822 //////////////////////////////////////////////////////////////////////////
821 //                                                823 //
822 // Calculate distance to shape from outside       824 // Calculate distance to shape from outside
823 //  - return kInfinity if no intersection         825 //  - return kInfinity if no intersection
824                                                   826 
825 G4double G4Trap::DistanceToIn(const G4ThreeVec    827 G4double G4Trap::DistanceToIn(const G4ThreeVector& p,
826                               const G4ThreeVec    828                               const G4ThreeVector& v ) const
827 {                                                 829 {
828   // Z intersections                              830   // Z intersections
829   //                                              831   //
830   if ((std::abs(p.z()) - fDz) >= -halfCarToler    832   if ((std::abs(p.z()) - fDz) >= -halfCarTolerance && p.z()*v.z() >= 0)
831     return kInfinity;                             833     return kInfinity;
832   G4double invz = (-v.z() == 0) ? DBL_MAX : -1    834   G4double invz = (-v.z() == 0) ? DBL_MAX : -1./v.z();
833   G4double dz = (invz < 0) ? fDz : -fDz;          835   G4double dz = (invz < 0) ? fDz : -fDz;
834   G4double tzmin = (p.z() + dz)*invz;             836   G4double tzmin = (p.z() + dz)*invz;
835   G4double tzmax = (p.z() - dz)*invz;             837   G4double tzmax = (p.z() - dz)*invz;
836                                                   838 
837   // Y intersections                              839   // Y intersections
838   //                                              840   //
839   G4double tymin = 0, tymax = DBL_MAX;            841   G4double tymin = 0, tymax = DBL_MAX;
840   G4int i = 0;                                    842   G4int i = 0;
841   for ( ; i<2; ++i)                               843   for ( ; i<2; ++i)
842   {                                               844   {
843     G4double cosa = fPlanes[i].b*v.y() + fPlan    845     G4double cosa = fPlanes[i].b*v.y() + fPlanes[i].c*v.z();
844     G4double dist = fPlanes[i].b*p.y() + fPlan    846     G4double dist = fPlanes[i].b*p.y() + fPlanes[i].c*p.z() + fPlanes[i].d;
845     if (dist >= -halfCarTolerance)                847     if (dist >= -halfCarTolerance)
846     {                                             848     {
847       if (cosa >= 0) return kInfinity;            849       if (cosa >= 0) return kInfinity;
848       G4double tmp  = -dist/cosa;                 850       G4double tmp  = -dist/cosa;
849       if (tymin < tmp) tymin = tmp;               851       if (tymin < tmp) tymin = tmp;
850     }                                             852     }
851     else if (cosa > 0)                            853     else if (cosa > 0)
852     {                                             854     {
853       G4double tmp  = -dist/cosa;                 855       G4double tmp  = -dist/cosa;
854       if (tymax > tmp) tymax = tmp;               856       if (tymax > tmp) tymax = tmp;
855     }                                             857     }
856   }                                               858   }
857                                                   859 
858   // Z intersections                              860   // Z intersections
859   //                                              861   //
860   G4double txmin = 0, txmax = DBL_MAX;            862   G4double txmin = 0, txmax = DBL_MAX;
861   for ( ; i<4; ++i)                               863   for ( ; i<4; ++i)
862   {                                               864   {
863     G4double cosa = fPlanes[i].a*v.x()+fPlanes    865     G4double cosa = fPlanes[i].a*v.x()+fPlanes[i].b*v.y()+fPlanes[i].c*v.z();
864     G4double dist = fPlanes[i].a*p.x()+fPlanes    866     G4double dist = fPlanes[i].a*p.x()+fPlanes[i].b*p.y()+fPlanes[i].c*p.z() +
865                     fPlanes[i].d;                 867                     fPlanes[i].d;
866     if (dist >= -halfCarTolerance)                868     if (dist >= -halfCarTolerance)
867     {                                             869     {
868       if (cosa >= 0) return kInfinity;            870       if (cosa >= 0) return kInfinity;
869       G4double tmp  = -dist/cosa;                 871       G4double tmp  = -dist/cosa;
870       if (txmin < tmp) txmin = tmp;               872       if (txmin < tmp) txmin = tmp;
871     }                                             873     }
872     else if (cosa > 0)                            874     else if (cosa > 0)
873     {                                             875     {
874       G4double tmp  = -dist/cosa;                 876       G4double tmp  = -dist/cosa;
875       if (txmax > tmp) txmax = tmp;               877       if (txmax > tmp) txmax = tmp;
876     }                                             878     }
877   }                                               879   }
878                                                   880 
879   // Find distance                                881   // Find distance
880   //                                              882   //
881   G4double tmin = std::max(std::max(txmin,tymi    883   G4double tmin = std::max(std::max(txmin,tymin),tzmin);
882   G4double tmax = std::min(std::min(txmax,tyma    884   G4double tmax = std::min(std::min(txmax,tymax),tzmax);
883                                                   885 
884   if (tmax <= tmin + halfCarTolerance) return     886   if (tmax <= tmin + halfCarTolerance) return kInfinity; // touch or no hit
885   return (tmin < halfCarTolerance ) ? 0. : tmi    887   return (tmin < halfCarTolerance ) ? 0. : tmin;
886 }                                                 888 }
887                                                   889 
888 //////////////////////////////////////////////    890 //////////////////////////////////////////////////////////////////////////
889 //                                                891 //
890 // Calculate exact shortest distance to any bo    892 // Calculate exact shortest distance to any boundary from outside
891 // This is the best fast estimation of the sho    893 // This is the best fast estimation of the shortest distance to trap
892 // - return 0 if point is inside                  894 // - return 0 if point is inside
893                                                   895 
894 G4double G4Trap::DistanceToIn( const G4ThreeVe    896 G4double G4Trap::DistanceToIn( const G4ThreeVector& p ) const
895 {                                                 897 {
896   switch (fTrapType)                              898   switch (fTrapType)
897   {                                               899   {
898     case 0: // General case                       900     case 0: // General case
899     {                                             901     {
900       G4double dz = std::abs(p.z())-fDz;          902       G4double dz = std::abs(p.z())-fDz;
901       G4double dy1 = fPlanes[0].b*p.y()+fPlane    903       G4double dy1 = fPlanes[0].b*p.y()+fPlanes[0].c*p.z()+fPlanes[0].d;
902       G4double dy2 = fPlanes[1].b*p.y()+fPlane    904       G4double dy2 = fPlanes[1].b*p.y()+fPlanes[1].c*p.z()+fPlanes[1].d;
903       G4double dy = std::max(dz,std::max(dy1,d    905       G4double dy = std::max(dz,std::max(dy1,dy2));
904                                                   906 
905       G4double dx1 = fPlanes[2].a*p.x()+fPlane    907       G4double dx1 = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()
906                    + fPlanes[2].c*p.z()+fPlane    908                    + fPlanes[2].c*p.z()+fPlanes[2].d;
907       G4double dx2 = fPlanes[3].a*p.x()+fPlane    909       G4double dx2 = fPlanes[3].a*p.x()+fPlanes[3].b*p.y()
908                    + fPlanes[3].c*p.z()+fPlane    910                    + fPlanes[3].c*p.z()+fPlanes[3].d;
909       G4double dist = std::max(dy,std::max(dx1    911       G4double dist = std::max(dy,std::max(dx1,dx2));
910       return (dist > 0) ? dist : 0.;              912       return (dist > 0) ? dist : 0.;
911     }                                             913     }
912     case 1: // YZ section is a rectangle          914     case 1: // YZ section is a rectangle
913     {                                             915     {
914       G4double dz = std::abs(p.z())-fDz;          916       G4double dz = std::abs(p.z())-fDz;
915       G4double dy = std::max(dz,std::abs(p.y()    917       G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
916       G4double dx1 = fPlanes[2].a*p.x()+fPlane    918       G4double dx1 = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()
917                    + fPlanes[2].c*p.z()+fPlane    919                    + fPlanes[2].c*p.z()+fPlanes[2].d;
918       G4double dx2 = fPlanes[3].a*p.x()+fPlane    920       G4double dx2 = fPlanes[3].a*p.x()+fPlanes[3].b*p.y()
919                    + fPlanes[3].c*p.z()+fPlane    921                    + fPlanes[3].c*p.z()+fPlanes[3].d;
920       G4double dist = std::max(dy,std::max(dx1    922       G4double dist = std::max(dy,std::max(dx1,dx2));
921       return (dist > 0) ? dist : 0.;              923       return (dist > 0) ? dist : 0.;
922     }                                             924     }
923     case 2: // YZ section is a rectangle and      925     case 2: // YZ section is a rectangle and
924     {       // XZ section is an isosceles trap    926     {       // XZ section is an isosceles trapezoid
925       G4double dz = std::abs(p.z())-fDz;          927       G4double dz = std::abs(p.z())-fDz;
926       G4double dy = std::max(dz,std::abs(p.y()    928       G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
927       G4double dx = fPlanes[3].a*std::abs(p.x(    929       G4double dx = fPlanes[3].a*std::abs(p.x())
928                   + fPlanes[3].c*p.z()+fPlanes    930                   + fPlanes[3].c*p.z()+fPlanes[3].d;
929       G4double dist = std::max(dy,dx);            931       G4double dist = std::max(dy,dx);
930       return (dist > 0) ? dist : 0.;              932       return (dist > 0) ? dist : 0.;
931     }                                             933     }
932     case 3: // YZ section is a rectangle and      934     case 3: // YZ section is a rectangle and
933     {       // XY section is an isosceles trap    935     {       // XY section is an isosceles trapezoid
934       G4double dz = std::abs(p.z())-fDz;          936       G4double dz = std::abs(p.z())-fDz;
935       G4double dy = std::max(dz,std::abs(p.y()    937       G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
936       G4double dx = fPlanes[3].a*std::abs(p.x(    938       G4double dx = fPlanes[3].a*std::abs(p.x())
937                   + fPlanes[3].b*p.y()+fPlanes    939                   + fPlanes[3].b*p.y()+fPlanes[3].d;
938       G4double dist = std::max(dy,dx);            940       G4double dist = std::max(dy,dx);
939       return (dist > 0) ? dist : 0.;              941       return (dist > 0) ? dist : 0.;
940     }                                             942     }
941   }                                               943   }
942   return 0.;                                      944   return 0.;
943 }                                                 945 }
944                                                   946 
945 //////////////////////////////////////////////    947 //////////////////////////////////////////////////////////////////////////
946 //                                                948 //
947 // Calculate distance to surface of shape from    949 // Calculate distance to surface of shape from inside and
948 // find normal at exit point, if required         950 // find normal at exit point, if required
949 // - when leaving the surface, return 0           951 // - when leaving the surface, return 0
950                                                   952 
951 G4double G4Trap::DistanceToOut(const G4ThreeVe    953 G4double G4Trap::DistanceToOut(const G4ThreeVector& p, const G4ThreeVector& v,
952                                const G4bool ca    954                                const G4bool calcNorm,
953                                      G4bool* v    955                                      G4bool* validNorm, G4ThreeVector* n) const
954 {                                                 956 {
955   // Z intersections                              957   // Z intersections
956   //                                              958   //
957   if ((std::abs(p.z()) - fDz) >= -halfCarToler    959   if ((std::abs(p.z()) - fDz) >= -halfCarTolerance && p.z()*v.z() > 0)
958   {                                               960   {
959     if (calcNorm)                                 961     if (calcNorm)
960     {                                             962     {
961       *validNorm = true;                          963       *validNorm = true;
962       n->set(0, 0, (p.z() < 0) ? -1 : 1);         964       n->set(0, 0, (p.z() < 0) ? -1 : 1);
963     }                                             965     }
964     return 0;                                     966     return 0;
965   }                                               967   }
966   G4double vz = v.z();                            968   G4double vz = v.z();
967   G4double tmax = (vz == 0) ? DBL_MAX : (std::    969   G4double tmax = (vz == 0) ? DBL_MAX : (std::copysign(fDz,vz) - p.z())/vz;
968   G4int iside = (vz < 0) ? -4 : -2; // little     970   G4int iside = (vz < 0) ? -4 : -2; // little trick: (-4+3)=-1, (-2+3)=+1
969                                                   971 
970   // Y intersections                              972   // Y intersections
971   //                                              973   //
972   G4int i = 0;                                    974   G4int i = 0;
973   for ( ; i<2; ++i)                               975   for ( ; i<2; ++i)
974   {                                               976   {
975     G4double cosa = fPlanes[i].b*v.y() + fPlan    977     G4double cosa = fPlanes[i].b*v.y() + fPlanes[i].c*v.z();
976     if (cosa > 0)                                 978     if (cosa > 0)
977     {                                             979     {
978       G4double dist = fPlanes[i].b*p.y() + fPl    980       G4double dist = fPlanes[i].b*p.y() + fPlanes[i].c*p.z() + fPlanes[i].d;
979       if (dist >= -halfCarTolerance)              981       if (dist >= -halfCarTolerance)
980       {                                           982       {
981         if (calcNorm)                             983         if (calcNorm)
982         {                                         984         {
983           *validNorm = true;                      985           *validNorm = true;
984           n->set(0, fPlanes[i].b, fPlanes[i].c    986           n->set(0, fPlanes[i].b, fPlanes[i].c);
985         }                                         987         }
986         return 0;                                 988         return 0;
987       }                                           989       }
988       G4double tmp = -dist/cosa;                  990       G4double tmp = -dist/cosa;
989       if (tmax > tmp) { tmax = tmp; iside = i;    991       if (tmax > tmp) { tmax = tmp; iside = i; }
990     }                                             992     }
991   }                                               993   }
992                                                   994 
993   // X intersections                              995   // X intersections
994   //                                              996   //
995   for ( ; i<4; ++i)                               997   for ( ; i<4; ++i)
996   {                                               998   {
997     G4double cosa = fPlanes[i].a*v.x()+fPlanes    999     G4double cosa = fPlanes[i].a*v.x()+fPlanes[i].b*v.y()+fPlanes[i].c*v.z();
998     if (cosa > 0)                                 1000     if (cosa > 0)
999     {                                             1001     {
1000       G4double dist = fPlanes[i].a*p.x() +       1002       G4double dist = fPlanes[i].a*p.x() +
1001                       fPlanes[i].b*p.y() + fP    1003                       fPlanes[i].b*p.y() + fPlanes[i].c*p.z() + fPlanes[i].d;
1002       if (dist >= -halfCarTolerance)             1004       if (dist >= -halfCarTolerance)
1003       {                                          1005       {
1004         if (calcNorm)                            1006         if (calcNorm)
1005         {                                        1007         {
1006            *validNorm = true;                    1008            *validNorm = true;
1007            n->set(fPlanes[i].a, fPlanes[i].b,    1009            n->set(fPlanes[i].a, fPlanes[i].b, fPlanes[i].c);
1008         }                                        1010         }
1009         return 0;                                1011         return 0;
1010       }                                          1012       }
1011       G4double tmp = -dist/cosa;                 1013       G4double tmp = -dist/cosa;
1012       if (tmax > tmp) { tmax = tmp; iside = i    1014       if (tmax > tmp) { tmax = tmp; iside = i; }
1013     }                                            1015     }
1014   }                                              1016   }
1015                                                  1017 
1016   // Set normal, if required, and return dist    1018   // Set normal, if required, and return distance
1017   //                                             1019   //
1018   if (calcNorm)                                  1020   if (calcNorm)
1019   {                                              1021   {
1020     *validNorm = true;                           1022     *validNorm = true;
1021     if (iside < 0)                               1023     if (iside < 0)
1022       n->set(0, 0, iside + 3); // (-4+3)=-1,     1024       n->set(0, 0, iside + 3); // (-4+3)=-1, (-2+3)=+1
1023     else                                         1025     else
1024       n->set(fPlanes[iside].a, fPlanes[iside]    1026       n->set(fPlanes[iside].a, fPlanes[iside].b, fPlanes[iside].c);
1025   }                                              1027   }
1026   return tmax;                                   1028   return tmax;
1027 }                                                1029 }
1028                                                  1030 
1029 /////////////////////////////////////////////    1031 //////////////////////////////////////////////////////////////////////////
1030 //                                               1032 //
1031 // Calculate exact shortest distance to any b    1033 // Calculate exact shortest distance to any boundary from inside
1032 // - Returns 0 is ThreeVector outside            1034 // - Returns 0 is ThreeVector outside
1033                                                  1035 
1034 G4double G4Trap::DistanceToOut( const G4Three    1036 G4double G4Trap::DistanceToOut( const G4ThreeVector& p ) const
1035 {                                                1037 {
1036 #ifdef G4CSGDEBUG                                1038 #ifdef G4CSGDEBUG
1037   if( Inside(p) == kOutside )                    1039   if( Inside(p) == kOutside )
1038   {                                              1040   {
1039     std::ostringstream message;                  1041     std::ostringstream message;
1040     G4long oldprc = message.precision(16);    << 1042     G4int oldprc = message.precision(16);
1041     message << "Point p is outside (!?) of so    1043     message << "Point p is outside (!?) of solid: " << GetName() << G4endl;
1042     message << "Position:\n";                    1044     message << "Position:\n";
1043     message << "   p.x() = " << p.x()/mm << "    1045     message << "   p.x() = " << p.x()/mm << " mm\n";
1044     message << "   p.y() = " << p.y()/mm << "    1046     message << "   p.y() = " << p.y()/mm << " mm\n";
1045     message << "   p.z() = " << p.z()/mm << "    1047     message << "   p.z() = " << p.z()/mm << " mm";
1046     G4cout.precision(oldprc);                    1048     G4cout.precision(oldprc);
1047     G4Exception("G4Trap::DistanceToOut(p)", "    1049     G4Exception("G4Trap::DistanceToOut(p)", "GeomSolids1002",
1048                 JustWarning, message );          1050                 JustWarning, message );
1049     DumpInfo();                                  1051     DumpInfo();
1050   }                                              1052   }
1051 #endif                                           1053 #endif
1052   switch (fTrapType)                             1054   switch (fTrapType)
1053   {                                              1055   {
1054     case 0: // General case                      1056     case 0: // General case
1055     {                                            1057     {
1056       G4double dz = std::abs(p.z())-fDz;         1058       G4double dz = std::abs(p.z())-fDz;
1057       G4double dy1 = fPlanes[0].b*p.y()+fPlan    1059       G4double dy1 = fPlanes[0].b*p.y()+fPlanes[0].c*p.z()+fPlanes[0].d;
1058       G4double dy2 = fPlanes[1].b*p.y()+fPlan    1060       G4double dy2 = fPlanes[1].b*p.y()+fPlanes[1].c*p.z()+fPlanes[1].d;
1059       G4double dy = std::max(dz,std::max(dy1,    1061       G4double dy = std::max(dz,std::max(dy1,dy2));
1060                                                  1062 
1061       G4double dx1 = fPlanes[2].a*p.x()+fPlan    1063       G4double dx1 = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()
1062                    + fPlanes[2].c*p.z()+fPlan    1064                    + fPlanes[2].c*p.z()+fPlanes[2].d;
1063       G4double dx2 = fPlanes[3].a*p.x()+fPlan    1065       G4double dx2 = fPlanes[3].a*p.x()+fPlanes[3].b*p.y()
1064                    + fPlanes[3].c*p.z()+fPlan    1066                    + fPlanes[3].c*p.z()+fPlanes[3].d;
1065       G4double dist = std::max(dy,std::max(dx    1067       G4double dist = std::max(dy,std::max(dx1,dx2));
1066       return (dist < 0) ? -dist : 0.;            1068       return (dist < 0) ? -dist : 0.;
1067     }                                            1069     }
1068     case 1: // YZ section is a rectangle         1070     case 1: // YZ section is a rectangle
1069     {                                            1071     {
1070       G4double dz = std::abs(p.z())-fDz;         1072       G4double dz = std::abs(p.z())-fDz;
1071       G4double dy = std::max(dz,std::abs(p.y(    1073       G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
1072       G4double dx1 = fPlanes[2].a*p.x()+fPlan    1074       G4double dx1 = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()
1073                    + fPlanes[2].c*p.z()+fPlan    1075                    + fPlanes[2].c*p.z()+fPlanes[2].d;
1074       G4double dx2 = fPlanes[3].a*p.x()+fPlan    1076       G4double dx2 = fPlanes[3].a*p.x()+fPlanes[3].b*p.y()
1075                    + fPlanes[3].c*p.z()+fPlan    1077                    + fPlanes[3].c*p.z()+fPlanes[3].d;
1076       G4double dist = std::max(dy,std::max(dx    1078       G4double dist = std::max(dy,std::max(dx1,dx2));
1077       return (dist < 0) ? -dist : 0.;            1079       return (dist < 0) ? -dist : 0.;
1078     }                                            1080     }
1079     case 2: // YZ section is a rectangle and     1081     case 2: // YZ section is a rectangle and
1080     {       // XZ section is an isosceles tra    1082     {       // XZ section is an isosceles trapezoid
1081       G4double dz = std::abs(p.z())-fDz;         1083       G4double dz = std::abs(p.z())-fDz;
1082       G4double dy = std::max(dz,std::abs(p.y(    1084       G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
1083       G4double dx = fPlanes[3].a*std::abs(p.x    1085       G4double dx = fPlanes[3].a*std::abs(p.x())
1084                   + fPlanes[3].c*p.z()+fPlane    1086                   + fPlanes[3].c*p.z()+fPlanes[3].d;
1085       G4double dist = std::max(dy,dx);           1087       G4double dist = std::max(dy,dx);
1086       return (dist < 0) ? -dist : 0.;            1088       return (dist < 0) ? -dist : 0.;
1087     }                                            1089     }
1088     case 3: // YZ section is a rectangle and     1090     case 3: // YZ section is a rectangle and
1089     {       // XY section is an isosceles tra    1091     {       // XY section is an isosceles trapezoid
1090       G4double dz = std::abs(p.z())-fDz;         1092       G4double dz = std::abs(p.z())-fDz;
1091       G4double dy = std::max(dz,std::abs(p.y(    1093       G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
1092       G4double dx = fPlanes[3].a*std::abs(p.x    1094       G4double dx = fPlanes[3].a*std::abs(p.x())
1093                   + fPlanes[3].b*p.y()+fPlane    1095                   + fPlanes[3].b*p.y()+fPlanes[3].d;
1094       G4double dist = std::max(dy,dx);           1096       G4double dist = std::max(dy,dx);
1095       return (dist < 0) ? -dist : 0.;            1097       return (dist < 0) ? -dist : 0.;
1096     }                                            1098     }
1097   }                                              1099   }
1098   return 0.;                                     1100   return 0.;
1099 }                                                1101 }
1100                                                  1102 
1101 /////////////////////////////////////////////    1103 //////////////////////////////////////////////////////////////////////////
1102 //                                               1104 //
1103 // GetEntityType                                 1105 // GetEntityType
1104                                                  1106 
1105 G4GeometryType G4Trap::GetEntityType() const     1107 G4GeometryType G4Trap::GetEntityType() const
1106 {                                                1108 {
1107   return {"G4Trap"};                          << 1109   return G4String("G4Trap");
1108 }                                             << 
1109                                               << 
1110 ///////////////////////////////////////////// << 
1111 //                                            << 
1112 // IsFaceted                                  << 
1113                                               << 
1114 G4bool G4Trap::IsFaceted() const              << 
1115 {                                             << 
1116   return true;                                << 
1117 }                                                1110 }
1118                                                  1111 
1119 /////////////////////////////////////////////    1112 //////////////////////////////////////////////////////////////////////////
1120 //                                               1113 //
1121 // Make a clone of the object                    1114 // Make a clone of the object
1122 //                                               1115 //
1123 G4VSolid* G4Trap::Clone() const                  1116 G4VSolid* G4Trap::Clone() const
1124 {                                                1117 {
1125   return new G4Trap(*this);                      1118   return new G4Trap(*this);
1126 }                                                1119 }
1127                                                  1120 
1128 /////////////////////////////////////////////    1121 //////////////////////////////////////////////////////////////////////////
1129 //                                               1122 //
1130 // Stream object contents to an output stream    1123 // Stream object contents to an output stream
1131                                                  1124 
1132 std::ostream& G4Trap::StreamInfo( std::ostrea    1125 std::ostream& G4Trap::StreamInfo( std::ostream& os ) const
1133 {                                                1126 {
1134   G4double phi    = GetPhi();                    1127   G4double phi    = GetPhi();   
1135   G4double theta  = GetTheta();                  1128   G4double theta  = GetTheta();
1136   G4double alpha1 = GetAlpha1();                 1129   G4double alpha1 = GetAlpha1();
1137   G4double alpha2 = GetAlpha2();                 1130   G4double alpha2 = GetAlpha2();
1138                                                  1131 
1139   G4long oldprc = os.precision(16);           << 1132   G4int oldprc = os.precision(16);
1140   os << "------------------------------------    1133   os << "-----------------------------------------------------------\n"
1141      << "    *** Dump for solid: " << GetName    1134      << "    *** Dump for solid: " << GetName() << " ***\n"
1142      << "    ================================    1135      << "    ===================================================\n"
1143      << " Solid type: G4Trap\n"                  1136      << " Solid type: G4Trap\n"
1144      << " Parameters:\n"                         1137      << " Parameters:\n"
1145      << "    half length Z: " << fDz/mm << "     1138      << "    half length Z: " << fDz/mm << " mm\n"
1146      << "    half length Y, face -Dz: " << fD    1139      << "    half length Y, face -Dz: " << fDy1/mm << " mm\n"
1147      << "    half length X, face -Dz, side -D    1140      << "    half length X, face -Dz, side -Dy1: " << fDx1/mm << " mm\n"
1148      << "    half length X, face -Dz, side +D    1141      << "    half length X, face -Dz, side +Dy1: " << fDx2/mm << " mm\n"
1149      << "    half length Y, face +Dz: " << fD    1142      << "    half length Y, face +Dz: " << fDy2/mm << " mm\n"
1150      << "    half length X, face +Dz, side -D    1143      << "    half length X, face +Dz, side -Dy2: " << fDx3/mm << " mm\n"
1151      << "    half length X, face +Dz, side +D    1144      << "    half length X, face +Dz, side +Dy2: " << fDx4/mm << " mm\n"
1152      << "    theta: " << theta/degree << " de    1145      << "    theta: " << theta/degree << " degrees\n"
1153      << "    phi:   " << phi/degree << " degr    1146      << "    phi:   " << phi/degree << " degrees\n"
1154      << "    alpha, face -Dz: " << alpha1/deg    1147      << "    alpha, face -Dz: " << alpha1/degree << " degrees\n"
1155      << "    alpha, face +Dz: " << alpha2/deg    1148      << "    alpha, face +Dz: " << alpha2/degree << " degrees\n"
1156      << "------------------------------------    1149      << "-----------------------------------------------------------\n";
1157   os.precision(oldprc);                          1150   os.precision(oldprc);
1158                                                  1151 
1159   return os;                                     1152   return os;
1160 }                                                1153 }
1161                                                  1154 
1162 /////////////////////////////////////////////    1155 //////////////////////////////////////////////////////////////////////////
1163 //                                               1156 //
1164 // Compute vertices from planes                  1157 // Compute vertices from planes
1165                                                  1158 
1166 void G4Trap::GetVertices(G4ThreeVector pt[8])    1159 void G4Trap::GetVertices(G4ThreeVector pt[8]) const
1167 {                                                1160 {
1168   for (G4int i=0; i<8; ++i)                      1161   for (G4int i=0; i<8; ++i)
1169   {                                              1162   {
1170     G4int iy = (i==0 || i==1 || i==4 || i==5)    1163     G4int iy = (i==0 || i==1 || i==4 || i==5) ? 0 : 1;
1171     G4int ix = (i==0 || i==2 || i==4 || i==6)    1164     G4int ix = (i==0 || i==2 || i==4 || i==6) ? 2 : 3;
1172     G4double z = (i < 4) ? -fDz : fDz;           1165     G4double z = (i < 4) ? -fDz : fDz;
1173     G4double y = -(fPlanes[iy].c*z + fPlanes[    1166     G4double y = -(fPlanes[iy].c*z + fPlanes[iy].d)/fPlanes[iy].b;
1174     G4double x = -(fPlanes[ix].b*y + fPlanes[    1167     G4double x = -(fPlanes[ix].b*y + fPlanes[ix].c*z
1175                    + fPlanes[ix].d)/fPlanes[i    1168                    + fPlanes[ix].d)/fPlanes[ix].a;
1176     pt[i].set(x,y,z);                            1169     pt[i].set(x,y,z);
1177   }                                              1170   }
1178 }                                                1171 }
1179                                                  1172 
1180 /////////////////////////////////////////////    1173 //////////////////////////////////////////////////////////////////////////
1181 //                                               1174 //
1182 // Generate random point on the surface          1175 // Generate random point on the surface
1183                                                  1176 
1184 G4ThreeVector G4Trap::GetPointOnSurface() con    1177 G4ThreeVector G4Trap::GetPointOnSurface() const
1185 {                                                1178 {
1186   // Set indeces                                 1179   // Set indeces
1187   constexpr G4int iface [6][4] =                 1180   constexpr G4int iface [6][4] =
1188     { {0,1,3,2}, {0,4,5,1}, {2,3,7,6}, {0,2,6    1181     { {0,1,3,2}, {0,4,5,1}, {2,3,7,6}, {0,2,6,4}, {1,5,7,3}, {4,6,7,5} };
1189                                                  1182 
1190   // Set vertices                                1183   // Set vertices
1191   G4ThreeVector pt[8];                           1184   G4ThreeVector pt[8];
1192   GetVertices(pt);                               1185   GetVertices(pt);
1193                                                  1186 
1194   // Select face                                 1187   // Select face
1195   //                                             1188   //
1196   G4double select = fAreas[5]*G4QuickRand();     1189   G4double select = fAreas[5]*G4QuickRand();
1197   G4int k = 5;                                   1190   G4int k = 5;
1198   k -= (G4int)(select <= fAreas[4]);          << 1191   k -= (select <= fAreas[4]);
1199   k -= (G4int)(select <= fAreas[3]);          << 1192   k -= (select <= fAreas[3]);
1200   k -= (G4int)(select <= fAreas[2]);          << 1193   k -= (select <= fAreas[2]);
1201   k -= (G4int)(select <= fAreas[1]);          << 1194   k -= (select <= fAreas[1]);
1202   k -= (G4int)(select <= fAreas[0]);          << 1195   k -= (select <= fAreas[0]);
1203                                                  1196 
1204   // Select sub-triangle                         1197   // Select sub-triangle
1205   //                                             1198   //
1206   G4int i0 = iface[k][0];                        1199   G4int i0 = iface[k][0];
1207   G4int i1 = iface[k][1];                        1200   G4int i1 = iface[k][1];
1208   G4int i2 = iface[k][2];                        1201   G4int i2 = iface[k][2];
1209   G4int i3 = iface[k][3];                        1202   G4int i3 = iface[k][3];
1210   G4double s2 = G4GeomTools::TriangleAreaNorm    1203   G4double s2 = G4GeomTools::TriangleAreaNormal(pt[i2],pt[i1],pt[i3]).mag();
1211   if (select > fAreas[k] - s2) i0 = i2;          1204   if (select > fAreas[k] - s2) i0 = i2;
1212                                                  1205 
1213   // Generate point                              1206   // Generate point
1214   //                                             1207   //
1215   G4double u = G4QuickRand();                    1208   G4double u = G4QuickRand();
1216   G4double v = G4QuickRand();                    1209   G4double v = G4QuickRand();
1217   if (u + v > 1.) { u = 1. - u; v = 1. - v; }    1210   if (u + v > 1.) { u = 1. - u; v = 1. - v; }
1218   return (1.-u-v)*pt[i0] + u*pt[i1] + v*pt[i3    1211   return (1.-u-v)*pt[i0] + u*pt[i1] + v*pt[i3];
1219 }                                                1212 }
1220                                                  1213 
1221 /////////////////////////////////////////////    1214 //////////////////////////////////////////////////////////////////////////
1222 //                                               1215 //
1223 // Methods for visualisation                     1216 // Methods for visualisation
1224                                                  1217 
1225 void G4Trap::DescribeYourselfTo ( G4VGraphics    1218 void G4Trap::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
1226 {                                                1219 {
1227   scene.AddSolid (*this);                        1220   scene.AddSolid (*this);
1228 }                                                1221 }
1229                                                  1222 
1230 G4Polyhedron* G4Trap::CreatePolyhedron () con    1223 G4Polyhedron* G4Trap::CreatePolyhedron () const
1231 {                                                1224 {
1232   G4double phi = std::atan2(fTthetaSphi, fTth    1225   G4double phi = std::atan2(fTthetaSphi, fTthetaCphi);
1233   G4double alpha1 = std::atan(fTalpha1);         1226   G4double alpha1 = std::atan(fTalpha1);
1234   G4double alpha2 = std::atan(fTalpha2);         1227   G4double alpha2 = std::atan(fTalpha2);
1235   G4double theta = std::atan(std::sqrt(fTthet    1228   G4double theta = std::atan(std::sqrt(fTthetaCphi*fTthetaCphi
1236                                       +fTthet    1229                                       +fTthetaSphi*fTthetaSphi));
1237                                                  1230 
1238   return new G4PolyhedronTrap(fDz, theta, phi    1231   return new G4PolyhedronTrap(fDz, theta, phi,
1239                               fDy1, fDx1, fDx    1232                               fDy1, fDx1, fDx2, alpha1,
1240                               fDy2, fDx3, fDx    1233                               fDy2, fDx3, fDx4, alpha2);
1241 }                                                1234 }
1242                                                  1235 
1243 #endif                                           1236 #endif
1244                                                  1237