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


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
                                                   >>  26 //
                                                   >>  27 // $Id: G4Trap.cc 69788 2013-05-15 12:06:57Z gcosmo $
                                                   >>  28 //
                                                   >>  29 // class G4Trap
                                                   >>  30 //
 26 // Implementation for G4Trap class                 31 // Implementation for G4Trap class
 27 //                                                 32 //
 28 // 21.03.95 P.Kent: Modified for `tolerant' ge <<  33 // History:
                                                   >>  34 //
                                                   >>  35 // 28.04.05 V.Grichine: new SurfaceNormal according to J. Apostolakis proposal 
                                                   >>  36 // 26.04.05 V.Grichine: new SurfaceNormal is default 
                                                   >>  37 // 19.04.05 V.Grichine: bug fixed in G4Trap("name",G4ThreeVector[8] vp)
                                                   >>  38 // 12.12.04 V.Grichine: SurfaceNormal with edges/vertices 
                                                   >>  39 // 15.11.04 V.Grichine: bug fixed in G4Trap("name",G4ThreeVector[8] vp)
                                                   >>  40 // 13.12.99 V.Grichine: bug fixed in DistanceToIn(p,v)
                                                   >>  41 // 19.11.99 V.Grichine: kUndef was added to Eside enum
                                                   >>  42 // 04.06.99 S.Giani: Fixed CalculateExtent in rotated case. 
                                                   >>  43 // 08.12.97 J.Allison: Added "nominal" constructor and method SetAllParameters.
                                                   >>  44 // 01.11.96 V.Grichine: Costructor for Right Angular Wedge from STEP, G4Trd/Para
 29 // 09.09.96 V.Grichine: Final modifications be     45 // 09.09.96 V.Grichine: Final modifications before to commit
 30 // 08.12.97 J.Allison: Added "nominal" constru <<  46 // 21.03.95 P.Kent: Modified for `tolerant' geometry
 31 // 28.04.05 V.Grichine: new SurfaceNormal acco <<  47 //
 32 // 18.04.17 E.Tcherniaev: complete revision, s <<  48 //////////////////////////////////////////////////////////////////////////////////// 
 33 // ------------------------------------------- << 
 34                                                    49 
 35 #include "G4Trap.hh"                               50 #include "G4Trap.hh"
 36                                                << 
 37 #if !defined(G4GEOM_USE_UTRAP)                 << 
 38                                                << 
 39 #include "globals.hh"                              51 #include "globals.hh"
 40 #include "G4GeomTools.hh"                      << 
 41                                                    52 
 42 #include "G4VoxelLimits.hh"                        53 #include "G4VoxelLimits.hh"
 43 #include "G4AffineTransform.hh"                    54 #include "G4AffineTransform.hh"
 44 #include "G4BoundingEnvelope.hh"               << 
 45                                                    55 
 46 #include "G4VPVParameterisation.hh"                56 #include "G4VPVParameterisation.hh"
 47                                                    57 
 48 #include "G4QuickRand.hh"                      <<  58 #include "Randomize.hh"
 49                                                    59 
 50 #include "G4VGraphicsScene.hh"                     60 #include "G4VGraphicsScene.hh"
 51 #include "G4Polyhedron.hh"                         61 #include "G4Polyhedron.hh"
                                                   >>  62 #include "G4NURBS.hh"
                                                   >>  63 #include "G4NURBSbox.hh"
 52                                                    64 
 53 using namespace CLHEP;                             65 using namespace CLHEP;
 54                                                    66 
                                                   >>  67 ////////////////////////////////////////////////////////////////////////
                                                   >>  68 //
                                                   >>  69 // Accuracy of coplanarity
                                                   >>  70 
                                                   >>  71 const G4double kCoplanar_Tolerance = 1E-4 ;
                                                   >>  72 
                                                   >>  73 //////////////////////////////////////////////////////////////////////////
                                                   >>  74 //
                                                   >>  75 // Private enum: Not for external use 
                                                   >>  76     
                                                   >>  77 enum Eside {kUndef,ks0,ks1,ks2,ks3,kPZ,kMZ};
                                                   >>  78 
 55 //////////////////////////////////////////////     79 //////////////////////////////////////////////////////////////////////////
 56 //                                                 80 //
 57 // Constructor - check and set half-widths as  <<  81 // Constructor - check and set half-widths as well as angles: 
 58 // final check of coplanarity                      82 // final check of coplanarity
 59                                                    83 
 60 G4Trap::G4Trap( const G4String& pName,             84 G4Trap::G4Trap( const G4String& pName,
 61                       G4double pDz,                85                       G4double pDz,
 62                       G4double pTheta, G4doubl     86                       G4double pTheta, G4double pPhi,
 63                       G4double pDy1, G4double      87                       G4double pDy1, G4double pDx1, G4double pDx2,
 64                       G4double pAlp1,              88                       G4double pAlp1,
 65                       G4double pDy2, G4double      89                       G4double pDy2, G4double pDx3, G4double pDx4,
 66                       G4double pAlp2 )         <<  90                       G4double pAlp2)
 67   : G4CSGSolid(pName), halfCarTolerance(0.5*kC <<  91   : G4CSGSolid(pName)
 68 {                                                  92 {
 69   fDz = pDz;                                   <<  93   if ( pDz <= 0 || pDy1 <= 0 || pDx1 <= 0 ||
 70   fTthetaCphi = std::tan(pTheta)*std::cos(pPhi <<  94        pDx2 <= 0 || pDy2 <= 0 || pDx3 <= 0 || pDx4 <= 0 )
 71   fTthetaSphi = std::tan(pTheta)*std::sin(pPhi <<  95   {
                                                   >>  96     std::ostringstream message;
                                                   >>  97     message << "Invalid length parameters for Solid: " << GetName() << G4endl
                                                   >>  98             << "        X - "
                                                   >>  99             << pDx1 << ", " << pDx2 << ", " << pDx3 << ", " << pDx4 << G4endl
                                                   >> 100             << "          Y - " << pDy1 << ", " << pDy2 << G4endl
                                                   >> 101             << "          Z - " << pDz;
                                                   >> 102     G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
                                                   >> 103                 FatalException, message);
                                                   >> 104   }
 72                                                   105 
 73   fDy1 = pDy1; fDx1 = pDx1; fDx2 = pDx2; fTalp << 106   fDz=pDz;
 74   fDy2 = pDy2; fDx3 = pDx3; fDx4 = pDx4; fTalp << 107   fTthetaCphi=std::tan(pTheta)*std::cos(pPhi);
                                                   >> 108   fTthetaSphi=std::tan(pTheta)*std::sin(pPhi);
                                                   >> 109       
                                                   >> 110   fDy1=pDy1;
                                                   >> 111   fDx1=pDx1;
                                                   >> 112   fDx2=pDx2;
                                                   >> 113   fTalpha1=std::tan(pAlp1);
                                                   >> 114      
                                                   >> 115   fDy2=pDy2;
                                                   >> 116   fDx3=pDx3;
                                                   >> 117   fDx4=pDx4;
                                                   >> 118   fTalpha2=std::tan(pAlp2);
 75                                                   119 
 76   CheckParameters();                           << 
 77   MakePlanes();                                   120   MakePlanes();
 78 }                                                 121 }
 79                                                   122 
 80 ////////////////////////////////////////////// << 123 ////////////////////////////////////////////////////////////////////////////
 81 //                                                124 //
 82 // Constructor - Design of trapezoid based on  << 125 // Constructor - Design of trapezoid based on 8 G4ThreeVector parameters, 
 83 // which are its vertices. Checking of planari << 126 // which are its vertices. Checking of planarity with preparation of 
 84 // fPlanes[] and than calculation of other mem    127 // fPlanes[] and than calculation of other members
 85                                                   128 
 86 G4Trap::G4Trap( const G4String& pName,            129 G4Trap::G4Trap( const G4String& pName,
 87                 const G4ThreeVector pt[8] )       130                 const G4ThreeVector pt[8] )
 88   : G4CSGSolid(pName), halfCarTolerance(0.5*kC << 131   : G4CSGSolid(pName)
 89 {                                                 132 {
                                                   >> 133   G4bool good;
                                                   >> 134 
 90   // Start with check of centering - the cente    135   // Start with check of centering - the center of gravity trap line
 91   // should cross the origin of frame             136   // should cross the origin of frame
 92   //                                           << 137 
 93   if (  pt[0].z() >= 0                         << 138   if (!(   pt[0].z() < 0 
 94         || pt[0].z() != pt[1].z()              << 139         && pt[0].z() == pt[1].z() && pt[0].z() == pt[2].z()
 95         || pt[0].z() != pt[2].z()              << 140         && pt[0].z() == pt[3].z()
 96         || pt[0].z() != pt[3].z()              << 141         && pt[4].z() > 0 
 97                                                << 142         && pt[4].z() == pt[5].z() && pt[4].z() == pt[6].z()
 98         || pt[4].z() <= 0                      << 143         && pt[4].z() == pt[7].z()
 99         || pt[4].z() != pt[5].z()              << 144         && std::fabs( pt[0].z() + pt[4].z() ) < kCarTolerance
100         || pt[4].z() != pt[6].z()              << 145         && pt[0].y() == pt[1].y() && pt[2].y() == pt[3].y()
101         || pt[4].z() != pt[7].z()              << 146         && pt[4].y() == pt[5].y() && pt[6].y() == pt[7].y()
102                                                << 147         && std::fabs( pt[0].y() + pt[2].y() + pt[4].y() + pt[6].y() ) < kCarTolerance 
103         || std::fabs( pt[0].z() + pt[4].z() )  << 148         && std::fabs( pt[0].x() + pt[1].x() + pt[4].x() + pt[5].x() + 
104                                                << 149            pt[2].x() + pt[3].x() + pt[6].x() + pt[7].x() ) < kCarTolerance ) )
105         || pt[0].y() != pt[1].y()              << 
106         || pt[2].y() != pt[3].y()              << 
107         || pt[4].y() != pt[5].y()              << 
108         || pt[6].y() != pt[7].y()              << 
109                                                << 
110         || std::fabs(pt[0].y()+pt[2].y()+pt[4] << 
111         || std::fabs(pt[0].x()+pt[1].x()+pt[4] << 
112                      pt[2].x()+pt[3].x()+pt[6] << 
113   {                                               150   {
114     std::ostringstream message;                   151     std::ostringstream message;
115     message << "Invalid vertice coordinates fo    152     message << "Invalid vertice coordinates for Solid: " << GetName();
116     G4Exception("G4Trap::G4Trap()", "GeomSolid    153     G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
117                 FatalException, message);         154                 FatalException, message);
118   }                                               155   }
                                                   >> 156     
                                                   >> 157   // Bottom side with normal approx. -Y
                                                   >> 158     
                                                   >> 159   good = MakePlane(pt[0],pt[4],pt[5],pt[1],fPlanes[0]);
119                                                   160 
120   // Set parameters                            << 161   if (!good)
121   //                                           << 162   {
122   fDz = (pt[7]).z();                           << 163     DumpInfo();
                                                   >> 164     G4Exception("G4Trap::G4Trap()", "GeomSolids0002", FatalException,
                                                   >> 165                 "Face at ~-Y not planar.");
                                                   >> 166   }
                                                   >> 167 
                                                   >> 168   // Top side with normal approx. +Y
                                                   >> 169     
                                                   >> 170   good = MakePlane(pt[2],pt[3],pt[7],pt[6],fPlanes[1]);
                                                   >> 171 
                                                   >> 172   if (!good)
                                                   >> 173   {
                                                   >> 174     std::ostringstream message;
                                                   >> 175     message << "Face at ~+Y not planar for Solid: " << GetName();
                                                   >> 176     G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
                                                   >> 177                 FatalException, message);
                                                   >> 178   }
                                                   >> 179 
                                                   >> 180   // Front side with normal approx. -X
                                                   >> 181     
                                                   >> 182   good = MakePlane(pt[0],pt[2],pt[6],pt[4],fPlanes[2]);
123                                                   183 
                                                   >> 184   if (!good)
                                                   >> 185   {
                                                   >> 186     std::ostringstream message;
                                                   >> 187     message << "Face at ~-X not planar for Solid: " << GetName();
                                                   >> 188     G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
                                                   >> 189                 FatalException, message);
                                                   >> 190   }
                                                   >> 191 
                                                   >> 192   // Back side iwth normal approx. +X
                                                   >> 193     
                                                   >> 194   good = MakePlane(pt[1],pt[5],pt[7],pt[3],fPlanes[3]);
                                                   >> 195   if (!good)
                                                   >> 196   {
                                                   >> 197     std::ostringstream message;
                                                   >> 198     message << "Face at ~+X not planar for Solid: " << GetName();
                                                   >> 199     G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
                                                   >> 200                 FatalException, message);
                                                   >> 201   }
                                                   >> 202   fDz = (pt[7]).z() ;
                                                   >> 203       
124   fDy1     = ((pt[2]).y()-(pt[1]).y())*0.5;       204   fDy1     = ((pt[2]).y()-(pt[1]).y())*0.5;
125   fDx1     = ((pt[1]).x()-(pt[0]).x())*0.5;       205   fDx1     = ((pt[1]).x()-(pt[0]).x())*0.5;
126   fDx2     = ((pt[3]).x()-(pt[2]).x())*0.5;       206   fDx2     = ((pt[3]).x()-(pt[2]).x())*0.5;
127   fTalpha1 = ((pt[2]).x()+(pt[3]).x()-(pt[1]).    207   fTalpha1 = ((pt[2]).x()+(pt[3]).x()-(pt[1]).x()-(pt[0]).x())*0.25/fDy1;
128                                                   208 
129   fDy2     = ((pt[6]).y()-(pt[5]).y())*0.5;       209   fDy2     = ((pt[6]).y()-(pt[5]).y())*0.5;
130   fDx3     = ((pt[5]).x()-(pt[4]).x())*0.5;       210   fDx3     = ((pt[5]).x()-(pt[4]).x())*0.5;
131   fDx4     = ((pt[7]).x()-(pt[6]).x())*0.5;       211   fDx4     = ((pt[7]).x()-(pt[6]).x())*0.5;
132   fTalpha2 = ((pt[6]).x()+(pt[7]).x()-(pt[5]).    212   fTalpha2 = ((pt[6]).x()+(pt[7]).x()-(pt[5]).x()-(pt[4]).x())*0.25/fDy2;
133                                                   213 
134   fTthetaCphi = ((pt[4]).x()+fDy2*fTalpha2+fDx    214   fTthetaCphi = ((pt[4]).x()+fDy2*fTalpha2+fDx3)/fDz;
135   fTthetaSphi = ((pt[4]).y()+fDy2)/fDz;           215   fTthetaSphi = ((pt[4]).y()+fDy2)/fDz;
136                                                << 
137   CheckParameters();                           << 
138   MakePlanes(pt);                              << 
139 }                                                 216 }
140                                                   217 
141 ////////////////////////////////////////////// << 218 //////////////////////////////////////////////////////////////////////////////
142 //                                                219 //
143 // Constructor for Right Angular Wedge from ST    220 // Constructor for Right Angular Wedge from STEP
144                                                   221 
145 G4Trap::G4Trap( const G4String& pName,            222 G4Trap::G4Trap( const G4String& pName,
146                       G4double pZ,                223                       G4double pZ,
147                       G4double pY,                224                       G4double pY,
148                       G4double pX, G4double pL    225                       G4double pX, G4double pLTX )
149   : G4CSGSolid(pName), halfCarTolerance(0.5*kC << 226   : G4CSGSolid(pName)
150 {                                                 227 {
151   fDz  = 0.5*pZ; fTthetaCphi = 0; fTthetaSphi  << 228   G4bool good;
152   fDy1 = 0.5*pY; fDx1 = 0.5*pX; fDx2 = 0.5*pLT << 
153   fDy2 = fDy1;   fDx3 = fDx1;   fDx4 = fDx2;   << 
154                                                   229 
155   CheckParameters();                           << 230   if ( pZ<=0 || pY<=0 || pX<=0 || pLTX<=0 || pLTX>pX )
156   MakePlanes();                                << 231   {
                                                   >> 232     std::ostringstream message;
                                                   >> 233     message << "Invalid length parameters for Solid: " << GetName();
                                                   >> 234     G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
                                                   >> 235                 FatalException, message);
                                                   >> 236   }
                                                   >> 237 
                                                   >> 238   fDz = 0.5*pZ ;
                                                   >> 239   fTthetaCphi = 0 ;
                                                   >> 240   fTthetaSphi = 0 ;
                                                   >> 241 
                                                   >> 242   fDy1 = 0.5*pY;
                                                   >> 243   fDx1 = 0.5*pX ;
                                                   >> 244   fDx2 = 0.5*pLTX;
                                                   >> 245   fTalpha1 =  0.5*(pLTX - pX)/pY;
                                                   >> 246 
                                                   >> 247   fDy2 = fDy1 ;
                                                   >> 248   fDx3 = fDx1;
                                                   >> 249   fDx4 = fDx2 ;
                                                   >> 250   fTalpha2 = fTalpha1 ;
                                                   >> 251 
                                                   >> 252   G4ThreeVector pt[8] ;
                                                   >> 253 
                                                   >> 254   pt[0]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
                                                   >> 255                       -fDz*fTthetaSphi-fDy1,-fDz);
                                                   >> 256   pt[1]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
                                                   >> 257                       -fDz*fTthetaSphi-fDy1,-fDz);
                                                   >> 258   pt[2]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
                                                   >> 259                       -fDz*fTthetaSphi+fDy1,-fDz);
                                                   >> 260   pt[3]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
                                                   >> 261                       -fDz*fTthetaSphi+fDy1,-fDz);
                                                   >> 262   pt[4]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
                                                   >> 263                       +fDz*fTthetaSphi-fDy2,+fDz);
                                                   >> 264   pt[5]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
                                                   >> 265                       +fDz*fTthetaSphi-fDy2,+fDz);
                                                   >> 266   pt[6]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
                                                   >> 267                       +fDz*fTthetaSphi+fDy2,+fDz);
                                                   >> 268   pt[7]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
                                                   >> 269                       +fDz*fTthetaSphi+fDy2,+fDz);
                                                   >> 270 
                                                   >> 271   // Bottom side with normal approx. -Y
                                                   >> 272   //
                                                   >> 273   good=MakePlane(pt[0],pt[4],pt[5],pt[1],fPlanes[0]);
                                                   >> 274   if (!good)
                                                   >> 275   {
                                                   >> 276     std::ostringstream message;
                                                   >> 277     message << "Face at ~-Y not planar for Solid: " << GetName();
                                                   >> 278     G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
                                                   >> 279                 FatalException, message);
                                                   >> 280   }
                                                   >> 281 
                                                   >> 282   // Top side with normal approx. +Y
                                                   >> 283   //
                                                   >> 284   good=MakePlane(pt[2],pt[3],pt[7],pt[6],fPlanes[1]);
                                                   >> 285   if (!good)
                                                   >> 286   {
                                                   >> 287     std::ostringstream message;
                                                   >> 288     message << "Face at ~+Y not planar for Solid: " << GetName();
                                                   >> 289     G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
                                                   >> 290                 FatalException, message);
                                                   >> 291   }
                                                   >> 292 
                                                   >> 293   // Front side with normal approx. -X
                                                   >> 294   //
                                                   >> 295   good=MakePlane(pt[0],pt[2],pt[6],pt[4],fPlanes[2]);
                                                   >> 296   if (!good)
                                                   >> 297   {
                                                   >> 298     std::ostringstream message;
                                                   >> 299     message << "Face at ~-X not planar for Solid: " << GetName();
                                                   >> 300     G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
                                                   >> 301                 FatalException, message);
                                                   >> 302   }
                                                   >> 303 
                                                   >> 304   // Back side iwth normal approx. +X
                                                   >> 305   //
                                                   >> 306   good=MakePlane(pt[1],pt[5],pt[7],pt[3],fPlanes[3]);
                                                   >> 307   if (!good)
                                                   >> 308   {
                                                   >> 309     std::ostringstream message;
                                                   >> 310     message << "Face at ~+X not planar for Solid: " << GetName();
                                                   >> 311     G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
                                                   >> 312                 FatalException, message);
                                                   >> 313   }
157 }                                                 314 }
158                                                   315 
159 ////////////////////////////////////////////// << 316 ///////////////////////////////////////////////////////////////////////////////
160 //                                                317 //
161 // Constructor for G4Trd                          318 // Constructor for G4Trd
162                                                   319 
163 G4Trap::G4Trap( const G4String& pName,            320 G4Trap::G4Trap( const G4String& pName,
164                       G4double pDx1,  G4double    321                       G4double pDx1,  G4double pDx2,
165                       G4double pDy1,  G4double    322                       G4double pDy1,  G4double pDy2,
166                       G4double pDz )              323                       G4double pDz )
167   : G4CSGSolid(pName), halfCarTolerance(0.5*kC << 324   : G4CSGSolid(pName)
168 {                                                 325 {
169   fDz  = pDz;  fTthetaCphi = 0; fTthetaSphi =  << 326   G4bool good;
170   fDy1 = pDy1; fDx1 = pDx1; fDx2 = pDx1; fTalp << 
171   fDy2 = pDy2; fDx3 = pDx2; fDx4 = pDx2; fTalp << 
172                                                   327 
173   CheckParameters();                           << 328   if ( pDz<=0 || pDy1<=0 || pDx1<=0 || pDx2<=0 || pDy2<=0 )
174   MakePlanes();                                << 329   {
                                                   >> 330     std::ostringstream message;
                                                   >> 331     message << "Invalid length parameters for Solid: " << GetName();
                                                   >> 332     G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
                                                   >> 333                 FatalException, message);
                                                   >> 334   }
                                                   >> 335 
                                                   >> 336   fDz = pDz;
                                                   >> 337   fTthetaCphi = 0 ;
                                                   >> 338   fTthetaSphi = 0 ;
                                                   >> 339       
                                                   >> 340   fDy1 = pDy1 ;
                                                   >> 341   fDx1 = pDx1 ;
                                                   >> 342   fDx2 = pDx1 ;
                                                   >> 343   fTalpha1 = 0 ;
                                                   >> 344      
                                                   >> 345   fDy2 = pDy2 ;
                                                   >> 346   fDx3 = pDx2 ;
                                                   >> 347   fDx4 = pDx2 ;
                                                   >> 348   fTalpha2 = 0 ;
                                                   >> 349 
                                                   >> 350   G4ThreeVector pt[8] ;
                                                   >> 351      
                                                   >> 352   pt[0]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
                                                   >> 353                       -fDz*fTthetaSphi-fDy1,-fDz);
                                                   >> 354   pt[1]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
                                                   >> 355                       -fDz*fTthetaSphi-fDy1,-fDz);
                                                   >> 356   pt[2]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
                                                   >> 357                       -fDz*fTthetaSphi+fDy1,-fDz);
                                                   >> 358   pt[3]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
                                                   >> 359                       -fDz*fTthetaSphi+fDy1,-fDz);
                                                   >> 360   pt[4]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
                                                   >> 361                       +fDz*fTthetaSphi-fDy2,+fDz);
                                                   >> 362   pt[5]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
                                                   >> 363                       +fDz*fTthetaSphi-fDy2,+fDz);
                                                   >> 364   pt[6]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
                                                   >> 365                       +fDz*fTthetaSphi+fDy2,+fDz);
                                                   >> 366   pt[7]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
                                                   >> 367                       +fDz*fTthetaSphi+fDy2,+fDz);
                                                   >> 368 
                                                   >> 369   // Bottom side with normal approx. -Y
                                                   >> 370   //
                                                   >> 371   good=MakePlane(pt[0],pt[4],pt[5],pt[1],fPlanes[0]);
                                                   >> 372   if (!good)
                                                   >> 373   {
                                                   >> 374     std::ostringstream message;
                                                   >> 375     message << "Face at ~-Y not planar for Solid: " << GetName();
                                                   >> 376     G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
                                                   >> 377                 FatalException, message);
                                                   >> 378   }
                                                   >> 379 
                                                   >> 380   // Top side with normal approx. +Y
                                                   >> 381   //
                                                   >> 382   good=MakePlane(pt[2],pt[3],pt[7],pt[6],fPlanes[1]);
                                                   >> 383   if (!good)
                                                   >> 384   {
                                                   >> 385     std::ostringstream message;
                                                   >> 386     message << "Face at ~+Y not planar for Solid: " << GetName();
                                                   >> 387     G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
                                                   >> 388                 FatalException, message);
                                                   >> 389   }
                                                   >> 390 
                                                   >> 391   // Front side with normal approx. -X
                                                   >> 392   //
                                                   >> 393   good=MakePlane(pt[0],pt[2],pt[6],pt[4],fPlanes[2]);
                                                   >> 394   if (!good)
                                                   >> 395   {
                                                   >> 396     std::ostringstream message;
                                                   >> 397     message << "Face at ~-X not planar for Solid: " << GetName();
                                                   >> 398     G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
                                                   >> 399                 FatalException, message);
                                                   >> 400   }
                                                   >> 401 
                                                   >> 402   // Back side iwth normal approx. +X
                                                   >> 403   //
                                                   >> 404   good=MakePlane(pt[1],pt[5],pt[7],pt[3],fPlanes[3]);
                                                   >> 405   if (!good)
                                                   >> 406   {
                                                   >> 407     std::ostringstream message;
                                                   >> 408     message << "Face at ~+X not planar for Solid: " << GetName();
                                                   >> 409     G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
                                                   >> 410                 FatalException, message);
                                                   >> 411   }
175 }                                                 412 }
176                                                   413 
177 ////////////////////////////////////////////// << 414 ////////////////////////////////////////////////////////////////////////////
178 //                                                415 //
179 // Constructor for G4Para                         416 // Constructor for G4Para
180                                                   417 
181 G4Trap::G4Trap( const G4String& pName,            418 G4Trap::G4Trap( const G4String& pName,
182                       G4double pDx, G4double p    419                       G4double pDx, G4double pDy,
183                       G4double pDz,               420                       G4double pDz,
184                       G4double pAlpha,            421                       G4double pAlpha,
185                       G4double pTheta, G4doubl    422                       G4double pTheta, G4double pPhi )
186   : G4CSGSolid(pName), halfCarTolerance(0.5*kC << 423   : G4CSGSolid(pName)       
187 {                                                 424 {
188   fDz = pDz;                                   << 425   G4bool good;
189   fTthetaCphi = std::tan(pTheta)*std::cos(pPhi << 426 
190   fTthetaSphi = std::tan(pTheta)*std::sin(pPhi << 427   if ( pDz<=0 || pDy<=0 || pDx<=0 )
                                                   >> 428   {
                                                   >> 429     std::ostringstream message;
                                                   >> 430     message << "Invalid length parameters for Solid: " << GetName();
                                                   >> 431     G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
                                                   >> 432                 FatalException, message);
                                                   >> 433   }
191                                                   434 
192   fDy1 = pDy; fDx1 = pDx; fDx2 = pDx; fTalpha1 << 435   fDz = pDz ;
193   fDy2 = pDy; fDx3 = pDx; fDx4 = pDx; fTalpha2 << 436   fTthetaCphi = std::tan(pTheta)*std::cos(pPhi) ;
                                                   >> 437   fTthetaSphi = std::tan(pTheta)*std::sin(pPhi) ;
                                                   >> 438      
                                                   >> 439   fDy1 = pDy ;
                                                   >> 440   fDx1 = pDx ;
                                                   >> 441   fDx2 = pDx ;
                                                   >> 442   fTalpha1 = std::tan(pAlpha) ;
                                                   >> 443     
                                                   >> 444   fDy2 = pDy ;
                                                   >> 445   fDx3 = pDx ;
                                                   >> 446   fDx4 = pDx ;
                                                   >> 447   fTalpha2 = fTalpha1 ;
                                                   >> 448 
                                                   >> 449   G4ThreeVector pt[8] ;
                                                   >> 450      
                                                   >> 451   pt[0]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
                                                   >> 452                       -fDz*fTthetaSphi-fDy1,-fDz);
                                                   >> 453   pt[1]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
                                                   >> 454                       -fDz*fTthetaSphi-fDy1,-fDz);
                                                   >> 455   pt[2]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
                                                   >> 456                       -fDz*fTthetaSphi+fDy1,-fDz);
                                                   >> 457   pt[3]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
                                                   >> 458                       -fDz*fTthetaSphi+fDy1,-fDz);
                                                   >> 459   pt[4]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
                                                   >> 460                       +fDz*fTthetaSphi-fDy2,+fDz);
                                                   >> 461   pt[5]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
                                                   >> 462                       +fDz*fTthetaSphi-fDy2,+fDz);
                                                   >> 463   pt[6]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
                                                   >> 464                       +fDz*fTthetaSphi+fDy2,+fDz);
                                                   >> 465   pt[7]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
                                                   >> 466                       +fDz*fTthetaSphi+fDy2,+fDz);
194                                                   467 
195   CheckParameters();                           << 468   // Bottom side with normal approx. -Y
196   MakePlanes();                                << 469   //
                                                   >> 470   good=MakePlane(pt[0],pt[4],pt[5],pt[1],fPlanes[0]);
                                                   >> 471   if (!good)
                                                   >> 472   {
                                                   >> 473     std::ostringstream message;
                                                   >> 474     message << "Face at ~-Y not planar for Solid: " << GetName();
                                                   >> 475     G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
                                                   >> 476                 FatalException, message);
                                                   >> 477   }
                                                   >> 478 
                                                   >> 479   // Top side with normal approx. +Y
                                                   >> 480   //
                                                   >> 481   good=MakePlane(pt[2],pt[3],pt[7],pt[6],fPlanes[1]);
                                                   >> 482   if (!good)
                                                   >> 483   {
                                                   >> 484     std::ostringstream message;
                                                   >> 485     message << "Face at ~+Y not planar for Solid: " << GetName();
                                                   >> 486     G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
                                                   >> 487                 FatalException, message);
                                                   >> 488   }
                                                   >> 489 
                                                   >> 490   // Front side with normal approx. -X
                                                   >> 491   //
                                                   >> 492   good=MakePlane(pt[0],pt[2],pt[6],pt[4],fPlanes[2]);
                                                   >> 493   if (!good)
                                                   >> 494   {
                                                   >> 495     std::ostringstream message;
                                                   >> 496     message << "Face at ~-X not planar for Solid: " << GetName();
                                                   >> 497     G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
                                                   >> 498                 FatalException, message);
                                                   >> 499   }
                                                   >> 500 
                                                   >> 501   // Back side iwth normal approx. +X
                                                   >> 502   //
                                                   >> 503   good=MakePlane(pt[1],pt[5],pt[7],pt[3],fPlanes[3]);
                                                   >> 504   if (!good)
                                                   >> 505   {
                                                   >> 506     std::ostringstream message;
                                                   >> 507     message << "Face at ~+X not planar for Solid: " << GetName();
                                                   >> 508     G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
                                                   >> 509                 FatalException, message);
                                                   >> 510   }
197 }                                                 511 }
198                                                   512 
199 ////////////////////////////////////////////// << 513 ///////////////////////////////////////////////////////////////////////////
200 //                                                514 //
201 // Nominal constructor for G4Trap whose parame    515 // Nominal constructor for G4Trap whose parameters are to be set by
202 // a G4VParamaterisation later.  Check and set    516 // a G4VParamaterisation later.  Check and set half-widths as well as
203 // angles: final check of coplanarity             517 // angles: final check of coplanarity
204                                                   518 
205 G4Trap::G4Trap( const G4String& pName )           519 G4Trap::G4Trap( const G4String& pName )
206   : G4CSGSolid (pName), halfCarTolerance(0.5*k << 520   : G4CSGSolid (pName), fDz(1.), fTthetaCphi(0.), fTthetaSphi(0.),
207     fDz(1.), fTthetaCphi(0.), fTthetaSphi(0.), << 
208     fDy1(1.), fDx1(1.), fDx2(1.), fTalpha1(0.)    521     fDy1(1.), fDx1(1.), fDx2(1.), fTalpha1(0.),
209     fDy2(1.), fDx3(1.), fDx4(1.), fTalpha2(0.)    522     fDy2(1.), fDx3(1.), fDx4(1.), fTalpha2(0.)
210 {                                                 523 {
211   MakePlanes();                                   524   MakePlanes();
212 }                                                 525 }
213                                                   526 
214 ////////////////////////////////////////////// << 527 ///////////////////////////////////////////////////////////////////////
215 //                                                528 //
216 // Fake default constructor - sets only member    529 // Fake default constructor - sets only member data and allocates memory
217 //                            for usage restri    530 //                            for usage restricted to object persistency.
218 //                                                531 //
219 G4Trap::G4Trap( __void__& a )                     532 G4Trap::G4Trap( __void__& a )
220   : G4CSGSolid(a), halfCarTolerance(0.5*kCarTo << 533   : G4CSGSolid(a), fDz(1.), fTthetaCphi(0.), fTthetaSphi(0.),
221     fDz(1.), fTthetaCphi(0.), fTthetaSphi(0.), << 
222     fDy1(1.), fDx1(1.), fDx2(1.), fTalpha1(0.)    534     fDy1(1.), fDx1(1.), fDx2(1.), fTalpha1(0.),
223     fDy2(1.), fDx3(1.), fDx4(1.), fTalpha2(0.)    535     fDy2(1.), fDx3(1.), fDx4(1.), fTalpha2(0.)
224 {                                                 536 {
225   MakePlanes();                                   537   MakePlanes();
226 }                                                 538 }
227                                                   539 
228 ////////////////////////////////////////////// << 540 ////////////////////////////////////////////////////////////////////////
229 //                                                541 //
230 // Destructor                                     542 // Destructor
231                                                   543 
232 G4Trap::~G4Trap() = default;                   << 544 G4Trap::~G4Trap()
                                                   >> 545 {
                                                   >> 546 }
233                                                   547 
234 //////////////////////////////////////////////    548 //////////////////////////////////////////////////////////////////////////
235 //                                                549 //
236 // Copy constructor                               550 // Copy constructor
237                                                   551 
238 G4Trap::G4Trap(const G4Trap& rhs)                 552 G4Trap::G4Trap(const G4Trap& rhs)
239   : G4CSGSolid(rhs), halfCarTolerance(rhs.half << 553   : G4CSGSolid(rhs), fDz(rhs.fDz),
240     fDz(rhs.fDz), fTthetaCphi(rhs.fTthetaCphi) << 554     fTthetaCphi(rhs.fTthetaCphi), fTthetaSphi(rhs.fTthetaSphi),
241     fDy1(rhs.fDy1), fDx1(rhs.fDx1), fDx2(rhs.f    555     fDy1(rhs.fDy1), fDx1(rhs.fDx1), fDx2(rhs.fDx2), fTalpha1(rhs.fTalpha1),
242     fDy2(rhs.fDy2), fDx3(rhs.fDx3), fDx4(rhs.f    556     fDy2(rhs.fDy2), fDx3(rhs.fDx3), fDx4(rhs.fDx4), fTalpha2(rhs.fTalpha2)
243 {                                                 557 {
244   for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs << 558   for (size_t i=0; i<4; ++i)
245   for (G4int i=0; i<6; ++i) { fAreas[i] = rhs. << 559   {
246   fTrapType = rhs.fTrapType;                   << 560     fPlanes[i].a = rhs.fPlanes[i].a;
                                                   >> 561     fPlanes[i].b = rhs.fPlanes[i].b;
                                                   >> 562     fPlanes[i].c = rhs.fPlanes[i].c;
                                                   >> 563     fPlanes[i].d = rhs.fPlanes[i].d;
                                                   >> 564   }
247 }                                                 565 }
248                                                   566 
249 //////////////////////////////////////////////    567 //////////////////////////////////////////////////////////////////////////
250 //                                                568 //
251 // Assignment operator                            569 // Assignment operator
252                                                   570 
253 G4Trap& G4Trap::operator = (const G4Trap& rhs) << 571 G4Trap& G4Trap::operator = (const G4Trap& rhs) 
254 {                                                 572 {
255   // Check assignment to self                     573   // Check assignment to self
256   //                                              574   //
257   if (this == &rhs)  { return *this; }            575   if (this == &rhs)  { return *this; }
258                                                   576 
259   // Copy base class data                         577   // Copy base class data
260   //                                              578   //
261   G4CSGSolid::operator=(rhs);                     579   G4CSGSolid::operator=(rhs);
262                                                   580 
263   // Copy data                                    581   // Copy data
264   //                                              582   //
265   halfCarTolerance = rhs.halfCarTolerance;     << 583   fDz = rhs.fDz;
266   fDz = rhs.fDz; fTthetaCphi = rhs.fTthetaCphi << 584   fTthetaCphi = rhs.fTthetaCphi; fTthetaSphi = rhs.fTthetaSphi;
267   fDy1 = rhs.fDy1; fDx1 = rhs.fDx1; fDx2 = rhs    585   fDy1 = rhs.fDy1; fDx1 = rhs.fDx1; fDx2 = rhs.fDx2; fTalpha1 = rhs.fTalpha1;
268   fDy2 = rhs.fDy2; fDx3 = rhs.fDx3; fDx4 = rhs    586   fDy2 = rhs.fDy2; fDx3 = rhs.fDx3; fDx4 = rhs.fDx4; fTalpha2 = rhs.fTalpha2;
269   for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs << 587   for (size_t i=0; i<4; ++i)
270   for (G4int i=0; i<6; ++i) { fAreas[i] = rhs. << 588   {
271   fTrapType = rhs.fTrapType;                   << 589     fPlanes[i].a = rhs.fPlanes[i].a;
                                                   >> 590     fPlanes[i].b = rhs.fPlanes[i].b;
                                                   >> 591     fPlanes[i].c = rhs.fPlanes[i].c;
                                                   >> 592     fPlanes[i].d = rhs.fPlanes[i].d;
                                                   >> 593   }
                                                   >> 594 
272   return *this;                                   595   return *this;
273 }                                                 596 }
274                                                   597 
275 ////////////////////////////////////////////// << 598 ///////////////////////////////////////////////////////////////////////
276 //                                                599 //
277 // Set all parameters, as for constructor - ch    600 // Set all parameters, as for constructor - check and set half-widths
278 // as well as angles: final check of coplanari    601 // as well as angles: final check of coplanarity
279                                                   602 
280 void G4Trap::SetAllParameters ( G4double pDz,     603 void G4Trap::SetAllParameters ( G4double pDz,
281                                 G4double pThet    604                                 G4double pTheta,
282                                 G4double pPhi,    605                                 G4double pPhi,
283                                 G4double pDy1,    606                                 G4double pDy1,
284                                 G4double pDx1,    607                                 G4double pDx1,
285                                 G4double pDx2,    608                                 G4double pDx2,
286                                 G4double pAlp1    609                                 G4double pAlp1,
287                                 G4double pDy2,    610                                 G4double pDy2,
288                                 G4double pDx3,    611                                 G4double pDx3,
289                                 G4double pDx4,    612                                 G4double pDx4,
290                                 G4double pAlp2    613                                 G4double pAlp2 )
291 {                                                 614 {
292   // Reset data of the base class              << 615   if ( pDz<=0 || pDy1<=0 || pDx1<=0 || pDx2<=0 || pDy2<=0 || pDx3<=0 || pDx4<=0 )
293   fCubicVolume = 0;                            << 
294   fSurfaceArea = 0;                            << 
295   fRebuildPolyhedron = true;                   << 
296                                                << 
297   // Set parameters                            << 
298   fDz = pDz;                                   << 
299   fTthetaCphi = std::tan(pTheta)*std::cos(pPhi << 
300   fTthetaSphi = std::tan(pTheta)*std::sin(pPhi << 
301                                                << 
302   fDy1 = pDy1; fDx1 = pDx1; fDx2 = pDx2; fTalp << 
303   fDy2 = pDy2; fDx3 = pDx3; fDx4 = pDx4; fTalp << 
304                                                << 
305   CheckParameters();                           << 
306   MakePlanes();                                << 
307 }                                              << 
308                                                << 
309 ////////////////////////////////////////////// << 
310 //                                             << 
311 // Check length parameters                     << 
312                                                << 
313 void G4Trap::CheckParameters()                 << 
314 {                                              << 
315   if (fDz<=0 ||                                << 
316       fDy1<=0 || fDx1<=0 || fDx2<=0 ||         << 
317       fDy2<=0 || fDx3<=0 || fDx4<=0)           << 
318   {                                               616   {
319     std::ostringstream message;                   617     std::ostringstream message;
320     message << "Invalid Length Parameters for  << 618     message << "Invalid Length Parameters for Solid: " << GetName() << G4endl
321             << "\n  X - " <<fDx1<<", "<<fDx2<< << 619             << "        X - "
322             << "\n  Y - " <<fDy1<<", "<<fDy2   << 620             << pDx1 << ", " << pDx2 << ", " << pDx3 << ", " << pDx4 << G4endl
323             << "\n  Z - " <<fDz;               << 621             << "          Y - " << pDy1 << ", " << pDy2 << G4endl
324     G4Exception("G4Trap::CheckParameters()", " << 622             << "          Z - " << pDz;
                                                   >> 623     G4Exception("G4Trap::SetAllParameters()", "GeomSolids0002",
325                 FatalException, message);         624                 FatalException, message);
326   }                                               625   }
                                                   >> 626   fCubicVolume= 0.;
                                                   >> 627   fSurfaceArea= 0.;
                                                   >> 628   fpPolyhedron = 0;
                                                   >> 629   fDz=pDz;
                                                   >> 630   fTthetaCphi=std::tan(pTheta)*std::cos(pPhi);
                                                   >> 631   fTthetaSphi=std::tan(pTheta)*std::sin(pPhi);
                                                   >> 632      
                                                   >> 633   fDy1=pDy1;
                                                   >> 634   fDx1=pDx1;
                                                   >> 635   fDx2=pDx2;
                                                   >> 636   fTalpha1=std::tan(pAlp1);
                                                   >> 637     
                                                   >> 638   fDy2=pDy2;
                                                   >> 639   fDx3=pDx3;
                                                   >> 640   fDx4=pDx4;
                                                   >> 641   fTalpha2=std::tan(pAlp2);
                                                   >> 642 
                                                   >> 643   MakePlanes();
327 }                                                 644 }
328                                                   645 
329 //////////////////////////////////////////////    646 //////////////////////////////////////////////////////////////////////////
330 //                                                647 //
331 // Compute vertices and set side planes        << 648 // Checking of coplanarity
332                                                   649 
333 void G4Trap::MakePlanes()                      << 650 G4bool G4Trap::MakePlanes()
334 {                                                 651 {
335   G4double DzTthetaCphi = fDz*fTthetaCphi;     << 652   G4bool good = true;
336   G4double DzTthetaSphi = fDz*fTthetaSphi;     << 
337   G4double Dy1Talpha1   = fDy1*fTalpha1;       << 
338   G4double Dy2Talpha2   = fDy2*fTalpha2;       << 
339                                                << 
340   G4ThreeVector pt[8] =                        << 
341   {                                            << 
342     G4ThreeVector(-DzTthetaCphi-Dy1Talpha1-fDx << 
343     G4ThreeVector(-DzTthetaCphi-Dy1Talpha1+fDx << 
344     G4ThreeVector(-DzTthetaCphi+Dy1Talpha1-fDx << 
345     G4ThreeVector(-DzTthetaCphi+Dy1Talpha1+fDx << 
346     G4ThreeVector( DzTthetaCphi-Dy2Talpha2-fDx << 
347     G4ThreeVector( DzTthetaCphi-Dy2Talpha2+fDx << 
348     G4ThreeVector( DzTthetaCphi+Dy2Talpha2-fDx << 
349     G4ThreeVector( DzTthetaCphi+Dy2Talpha2+fDx << 
350   };                                           << 
351                                                   653 
352   MakePlanes(pt);                              << 654   G4ThreeVector pt[8] ;
353 }                                              << 655      
                                                   >> 656   pt[0]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
                                                   >> 657                       -fDz*fTthetaSphi-fDy1,-fDz);
                                                   >> 658   pt[1]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
                                                   >> 659                       -fDz*fTthetaSphi-fDy1,-fDz);
                                                   >> 660   pt[2]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
                                                   >> 661                       -fDz*fTthetaSphi+fDy1,-fDz);
                                                   >> 662   pt[3]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
                                                   >> 663                       -fDz*fTthetaSphi+fDy1,-fDz);
                                                   >> 664   pt[4]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
                                                   >> 665                       +fDz*fTthetaSphi-fDy2,+fDz);
                                                   >> 666   pt[5]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
                                                   >> 667                       +fDz*fTthetaSphi-fDy2,+fDz);
                                                   >> 668   pt[6]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
                                                   >> 669                       +fDz*fTthetaSphi+fDy2,+fDz);
                                                   >> 670   pt[7]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
                                                   >> 671                       +fDz*fTthetaSphi+fDy2,+fDz);
354                                                   672 
355 ////////////////////////////////////////////// << 673   // Bottom side with normal approx. -Y
356 //                                             << 674   //
357 // Set side planes, check planarity            << 675   good=MakePlane(pt[0],pt[4],pt[5],pt[1],fPlanes[0]) ;
358                                                << 676   if (!good)
359 void G4Trap::MakePlanes(const G4ThreeVector pt << 677   {
360 {                                              << 678     std::ostringstream message;
361   constexpr G4int iface[4][4] = { {0,4,5,1}, { << 679     message << "Face at ~-Y not planar for Solid: " << GetName();
362   const static G4String side[4] = { "~-Y", "~+ << 680     G4Exception("G4Trap::MakePlanes()", "GeomSolids0002",
                                                   >> 681                 FatalException, message);
                                                   >> 682   }
363                                                   683 
364   for (G4int i=0; i<4; ++i)                    << 684   // Top side with normal approx. +Y
                                                   >> 685   //
                                                   >> 686   good=MakePlane(pt[2],pt[3],pt[7],pt[6],fPlanes[1]);
                                                   >> 687   if (!good)
365   {                                               688   {
366     if (MakePlane(pt[iface[i][0]],             << 689     std::ostringstream message;
367                   pt[iface[i][1]],             << 690     message << "Face at ~+Y not planar for Solid: " << GetName();
368                   pt[iface[i][2]],             << 691     G4Exception("G4Trap::MakePlanes()", "GeomSolids0002",
369                   pt[iface[i][3]],             << 692                 FatalException, message);
370                   fPlanes[i])) continue;       << 693   }
371                                                   694 
372     // Non planar side face                    << 695   // Front side with normal approx. -X
373     G4ThreeVector normal(fPlanes[i].a,fPlanes[ << 696   //
374     G4double dmax = 0;                         << 697   good=MakePlane(pt[0],pt[2],pt[6],pt[4],fPlanes[2]);
375     for (G4int k=0; k<4; ++k)                  << 698   if (!good)
376     {                                          << 699   {
377       G4double dist = normal.dot(pt[iface[i][k << 
378       if (std::abs(dist) > std::abs(dmax)) dma << 
379     }                                          << 
380     std::ostringstream message;                   700     std::ostringstream message;
381     message << "Side face " << side[i] << " is << 701     message << "Face at ~-X not planar for Solid: " << GetName();
382             << GetName() << "\nDiscrepancy: "  << 702     G4Exception("G4Trap::MakePlanes()", "GeomSolids0002",
383     StreamInfo(message);                       << 703                 FatalException, message);
                                                   >> 704   }
                                                   >> 705    
                                                   >> 706   // Back side iwth normal approx. +X
                                                   >> 707   //
                                                   >> 708   good = MakePlane(pt[1],pt[5],pt[7],pt[3],fPlanes[3]);
                                                   >> 709   if ( !good )
                                                   >> 710   {
                                                   >> 711     std::ostringstream message;
                                                   >> 712     message << "Face at ~+X not planar for Solid: " << GetName();
384     G4Exception("G4Trap::MakePlanes()", "GeomS    713     G4Exception("G4Trap::MakePlanes()", "GeomSolids0002",
385                 FatalException, message);         714                 FatalException, message);
386   }                                               715   }
387                                                   716 
388   // Re-compute parameters                     << 717   return good;
389   SetCachedValues();                           << 
390 }                                                 718 }
391                                                   719 
392 ////////////////////////////////////////////// << 720 //////////////////////////////////////////////////////////////////////////////
393 //                                                721 //
394 // Calculate the coef's of the plane p1->p2->p    722 // Calculate the coef's of the plane p1->p2->p3->p4->p1
395 // where the ThreeVectors 1-4 are in anti-cloc << 723 // where the ThreeVectors 1-4 are in anti-clockwise order when viewed from
396 // from infront of the plane (i.e. from normal << 724 // infront of the plane (i.e. from normal direction).
397 //                                                725 //
398 // Return true if the points are coplanar, fal << 726 // Return true if the ThreeVectors are coplanar + set coef;s
                                                   >> 727 //        false if ThreeVectors are not coplanar
399                                                   728 
400 G4bool G4Trap::MakePlane( const G4ThreeVector&    729 G4bool G4Trap::MakePlane( const G4ThreeVector& p1,
401                           const G4ThreeVector&    730                           const G4ThreeVector& p2,
402                           const G4ThreeVector&    731                           const G4ThreeVector& p3,
403                           const G4ThreeVector&    732                           const G4ThreeVector& p4,
404                                 TrapSidePlane&    733                                 TrapSidePlane& plane )
405 {                                                 734 {
406   G4ThreeVector normal = ((p4 - p2).cross(p3 - << 735   G4double a, b, c, sd;
407   if (std::abs(normal.x()) < DBL_EPSILON) norm << 736   G4ThreeVector v12, v13, v14, Vcross;
408   if (std::abs(normal.y()) < DBL_EPSILON) norm << 
409   if (std::abs(normal.z()) < DBL_EPSILON) norm << 
410   normal = normal.unit();                      << 
411                                                << 
412   G4ThreeVector centre = (p1 + p2 + p3 + p4)*0 << 
413   plane.a =  normal.x();                       << 
414   plane.b =  normal.y();                       << 
415   plane.c =  normal.z();                       << 
416   plane.d = -normal.dot(centre);               << 
417                                                << 
418   // compute distances and check planarity     << 
419   G4double d1 = std::abs(normal.dot(p1) + plan << 
420   G4double d2 = std::abs(normal.dot(p2) + plan << 
421   G4double d3 = std::abs(normal.dot(p3) + plan << 
422   G4double d4 = std::abs(normal.dot(p4) + plan << 
423   G4double dmax = std::max(std::max(std::max(d << 
424                                                   737 
425   return dmax <= 1000 * kCarTolerance;         << 738   G4bool good;
426 }                                              << 
427                                                   739 
428 ////////////////////////////////////////////// << 740   v12    = p2 - p1;
429 //                                             << 741   v13    = p3 - p1;
430 // Recompute parameters using planes           << 742   v14    = p4 - p1;
431                                                << 743   Vcross = v12.cross(v13);
432 void G4Trap::SetCachedValues()                 << 
433 {                                              << 
434   // Set indeces                               << 
435   constexpr  G4int iface[6][4] =               << 
436       { {0,1,3,2}, {0,4,5,1}, {2,3,7,6}, {0,2, << 
437                                                   744 
438   // Get vertices                              << 745   if (std::fabs(Vcross.dot(v14)/(Vcross.mag()*v14.mag())) > kCoplanar_Tolerance)
439   G4ThreeVector pt[8];                         << 
440   GetVertices(pt);                             << 
441                                                << 
442   // Set face areas                            << 
443   for (G4int i=0; i<6; ++i)                    << 
444   {                                               746   {
445     fAreas[i] = G4GeomTools::QuadAreaNormal(pt << 747     good = false;
446                                             pt << 
447                                             pt << 
448                                             pt << 
449   }                                            << 
450   for (G4int i=1; i<6; ++i) { fAreas[i] += fAr << 
451                                                << 
452   // Define type of trapezoid                  << 
453   fTrapType = 0;                               << 
454   if (fPlanes[0].b  == -1 && fPlanes[1].b == 1 << 
455       std::abs(fPlanes[0].a) < DBL_EPSILON &&  << 
456       std::abs(fPlanes[0].c) < DBL_EPSILON &&  << 
457       std::abs(fPlanes[1].a) < DBL_EPSILON &&  << 
458       std::abs(fPlanes[1].c) < DBL_EPSILON)    << 
459   {                                            << 
460     fTrapType = 1; // YZ section is a rectangl << 
461     if (std::abs(fPlanes[2].a + fPlanes[3].a)  << 
462         std::abs(fPlanes[2].c - fPlanes[3].c)  << 
463         fPlanes[2].b == 0 &&                   << 
464         fPlanes[3].b == 0)                     << 
465     {                                          << 
466       fTrapType = 2; // ... and XZ section is  << 
467       fPlanes[2].a = -fPlanes[3].a;            << 
468       fPlanes[2].c =  fPlanes[3].c;            << 
469     }                                          << 
470     if (std::abs(fPlanes[2].a + fPlanes[3].a)  << 
471         std::abs(fPlanes[2].b - fPlanes[3].b)  << 
472         fPlanes[2].c == 0 &&                   << 
473         fPlanes[3].c == 0)                     << 
474     {                                          << 
475       fTrapType = 3; // ... and XY section is  << 
476       fPlanes[2].a = -fPlanes[3].a;            << 
477       fPlanes[2].b =  fPlanes[3].b;            << 
478     }                                          << 
479   }                                               748   }
480 }                                              << 749   else
481                                                << 
482 ////////////////////////////////////////////// << 
483 //                                             << 
484 // Get volume                                  << 
485                                                << 
486 G4double G4Trap::GetCubicVolume()              << 
487 {                                              << 
488   if (fCubicVolume == 0)                       << 
489   {                                            << 
490     G4ThreeVector pt[8];                       << 
491     GetVertices(pt);                           << 
492                                                << 
493     G4double dz  = pt[4].z() - pt[0].z();      << 
494     G4double dy1 = pt[2].y() - pt[0].y();      << 
495     G4double dx1 = pt[1].x() - pt[0].x();      << 
496     G4double dx2 = pt[3].x() - pt[2].x();      << 
497     G4double dy2 = pt[6].y() - pt[4].y();      << 
498     G4double dx3 = pt[5].x() - pt[4].x();      << 
499     G4double dx4 = pt[7].x() - pt[6].x();      << 
500                                                << 
501     fCubicVolume = ((dx1 + dx2 + dx3 + dx4)*(d << 
502                     (dx4 + dx3 - dx2 - dx1)*(d << 
503   }                                            << 
504   return fCubicVolume;                         << 
505 }                                              << 
506                                                << 
507 ////////////////////////////////////////////// << 
508 //                                             << 
509 // Get surface area                            << 
510                                                << 
511 G4double G4Trap::GetSurfaceArea()              << 
512 {                                              << 
513   if (fSurfaceArea == 0)                       << 
514   {                                               750   {
515     G4ThreeVector pt[8];                       << 751     // a,b,c correspond to the x/y/z components of the
516     G4int iface [6][4] =                       << 752     // normal vector to the plane
517       { {0,1,3,2}, {0,4,5,1}, {2,3,7,6}, {0,2, << 753      
518                                                << 754     // a  = (p2.y()-p1.y())*(p1.z()+p2.z())+(p3.y()-p2.y())*(p2.z()+p3.z());
519     GetVertices(pt);                           << 755     // a += (p4.y()-p3.y())*(p3.z()+p4.z())+(p1.y()-p4.y())*(p4.z()+p1.z()); // ?   
520     for (const auto & i : iface)               << 756     // b  = (p2.z()-p1.z())*(p1.x()+p2.x())+(p3.z()-p2.z())*(p2.x()+p3.x());
                                                   >> 757     // b += (p4.z()-p3.z())*(p3.x()+p4.x())+(p1.z()-p4.z())*(p4.x()+p1.x()); // ?      
                                                   >> 758     // c  = (p2.x()-p1.x())*(p1.y()+p2.y())+(p3.x()-p2.x())*(p2.y()+p3.y());
                                                   >> 759     // c += (p4.x()-p3.x())*(p3.y()+p4.y())+(p1.x()-p4.x())*(p4.y()+p1.y()); // ?
                                                   >> 760 
                                                   >> 761     // Let create diagonals 4-2 and 3-1 than (4-2)x(3-1) provides
                                                   >> 762     // vector perpendicular to the plane directed to outside !!!
                                                   >> 763     // and a,b,c, = f(1,2,3,4) external relative to trap normal
                                                   >> 764 
                                                   >> 765     a = +(p4.y() - p2.y())*(p3.z() - p1.z())
                                                   >> 766         - (p3.y() - p1.y())*(p4.z() - p2.z());
                                                   >> 767 
                                                   >> 768     b = -(p4.x() - p2.x())*(p3.z() - p1.z())
                                                   >> 769         + (p3.x() - p1.x())*(p4.z() - p2.z());
                                                   >> 770  
                                                   >> 771     c = +(p4.x() - p2.x())*(p3.y() - p1.y())
                                                   >> 772         - (p3.x() - p1.x())*(p4.y() - p2.y());
                                                   >> 773 
                                                   >> 774     sd = std::sqrt( a*a + b*b + c*c ); // so now vector plane.(a,b,c) is unit 
                                                   >> 775 
                                                   >> 776     if( sd > 0 )
                                                   >> 777     {
                                                   >> 778       plane.a = a/sd;
                                                   >> 779       plane.b = b/sd;
                                                   >> 780       plane.c = c/sd;
                                                   >> 781     }
                                                   >> 782     else
521     {                                             783     {
522       fSurfaceArea += G4GeomTools::QuadAreaNor << 784       std::ostringstream message;
523                                                << 785       message << "Invalid parameters: norm.mod() <= 0, for Solid: "
524                                                << 786               << GetName();
525                                                << 787       G4Exception("G4Trap::MakePlanes()", "GeomSolids0002",
                                                   >> 788                   FatalException, message) ;
526     }                                             789     }
                                                   >> 790     // Calculate D: p1 in in plane so D=-n.p1.Vect()
                                                   >> 791     
                                                   >> 792     plane.d = -( plane.a*p1.x() + plane.b*p1.y() + plane.c*p1.z() );
                                                   >> 793 
                                                   >> 794     good = true;
527   }                                               795   }
528   return fSurfaceArea;                         << 796   return good;
529 }                                                 797 }
530                                                   798 
531 ////////////////////////////////////////////// << 799 //////////////////////////////////////////////////////////////////////////////
532 //                                                800 //
533 // Dispatch to parameterisation for replicatio    801 // Dispatch to parameterisation for replication mechanism dimension
534 // computation & modification.                    802 // computation & modification.
535                                                   803 
536 void G4Trap::ComputeDimensions(       G4VPVPar    804 void G4Trap::ComputeDimensions(       G4VPVParameterisation* p,
537                                 const G4int n,    805                                 const G4int n,
538                                 const G4VPhysi    806                                 const G4VPhysicalVolume* pRep )
539 {                                                 807 {
540   p->ComputeDimensions(*this,n,pRep);             808   p->ComputeDimensions(*this,n,pRep);
541 }                                                 809 }
542                                                   810 
543 ////////////////////////////////////////////// << 
544 //                                             << 
545 // Get bounding box                            << 
546                                                << 
547 void G4Trap::BoundingLimits(G4ThreeVector& pMi << 
548 {                                              << 
549   G4ThreeVector pt[8];                         << 
550   GetVertices(pt);                             << 
551                                                << 
552   G4double xmin = kInfinity, xmax = -kInfinity << 
553   G4double ymin = kInfinity, ymax = -kInfinity << 
554   for (const auto & i : pt)                    << 
555   {                                            << 
556     G4double x = i.x();                        << 
557     if (x < xmin) xmin = x;                    << 
558     if (x > xmax) xmax = x;                    << 
559     G4double y = i.y();                        << 
560     if (y < ymin) ymin = y;                    << 
561     if (y > ymax) ymax = y;                    << 
562   }                                            << 
563                                                << 
564   G4double dz   = GetZHalfLength();            << 
565   pMin.set(xmin,ymin,-dz);                     << 
566   pMax.set(xmax,ymax, dz);                     << 
567                                                << 
568   // Check correctness of the bounding box     << 
569   //                                           << 
570   if (pMin.x() >= pMax.x() || pMin.y() >= pMax << 
571   {                                            << 
572     std::ostringstream message;                << 
573     message << "Bad bounding box (min >= max)  << 
574             << GetName() << " !"               << 
575             << "\npMin = " << pMin             << 
576             << "\npMax = " << pMax;            << 
577     G4Exception("G4Trap::BoundingLimits()", "G << 
578                 JustWarning, message);         << 
579     DumpInfo();                                << 
580   }                                            << 
581 }                                              << 
582                                                   811 
583 ////////////////////////////////////////////// << 812 ////////////////////////////////////////////////////////////////////////
584 //                                                813 //
585 // Calculate extent under transform and specif    814 // Calculate extent under transform and specified limit
586                                                   815 
587 G4bool G4Trap::CalculateExtent( const EAxis pA    816 G4bool G4Trap::CalculateExtent( const EAxis pAxis,
588                                 const G4VoxelL    817                                 const G4VoxelLimits& pVoxelLimit,
589                                 const G4Affine    818                                 const G4AffineTransform& pTransform,
590                                       G4double    819                                       G4double& pMin, G4double& pMax) const
591 {                                                 820 {
592   G4ThreeVector bmin, bmax;                    << 821   G4double xMin, xMax, yMin, yMax, zMin, zMax;
593   G4bool exist;                                << 822   G4bool flag;
594                                                   823 
595   // Check bounding box (bbox)                 << 824   if (!pTransform.IsRotated())
596   //                                           << 825   {  
597   BoundingLimits(bmin,bmax);                   << 826     // Special case handling for unrotated trapezoids
598   G4BoundingEnvelope bbox(bmin,bmax);          << 827     // Compute z/x/y/ mins and maxs respecting limits, with early returns
599 #ifdef G4BBOX_EXTENT                           << 828     // if outside limits. Then switch() on pAxis
600   return bbox.CalculateExtent(pAxis,pVoxelLimi << 829 
601 #endif                                         << 830     G4int i ; 
602   if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox << 831     G4double xoffset;
603   {                                            << 832     G4double yoffset;
604     return exist = pMin < pMax;                << 833     G4double zoffset;
                                                   >> 834     G4double temp[8] ;     // some points for intersection with zMin/zMax
                                                   >> 835     G4ThreeVector pt[8];   // vertices after translation
                                                   >> 836     
                                                   >> 837     xoffset=pTransform.NetTranslation().x();      
                                                   >> 838     yoffset=pTransform.NetTranslation().y();
                                                   >> 839     zoffset=pTransform.NetTranslation().z();
                                                   >> 840  
                                                   >> 841     pt[0]=G4ThreeVector(xoffset-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
                                                   >> 842                         yoffset-fDz*fTthetaSphi-fDy1,zoffset-fDz);
                                                   >> 843     pt[1]=G4ThreeVector(xoffset-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
                                                   >> 844                         yoffset-fDz*fTthetaSphi-fDy1,zoffset-fDz);
                                                   >> 845     pt[2]=G4ThreeVector(xoffset-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
                                                   >> 846                         yoffset-fDz*fTthetaSphi+fDy1,zoffset-fDz);
                                                   >> 847     pt[3]=G4ThreeVector(xoffset-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
                                                   >> 848                         yoffset-fDz*fTthetaSphi+fDy1,zoffset-fDz);
                                                   >> 849     pt[4]=G4ThreeVector(xoffset+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
                                                   >> 850                         yoffset+fDz*fTthetaSphi-fDy2,zoffset+fDz);
                                                   >> 851     pt[5]=G4ThreeVector(xoffset+fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
                                                   >> 852                         yoffset+fDz*fTthetaSphi-fDy2,zoffset+fDz);
                                                   >> 853     pt[6]=G4ThreeVector(xoffset+fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
                                                   >> 854                         yoffset+fDz*fTthetaSphi+fDy2,zoffset+fDz);
                                                   >> 855     pt[7]=G4ThreeVector(xoffset+fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
                                                   >> 856                         yoffset+fDz*fTthetaSphi+fDy2,zoffset+fDz);
                                                   >> 857     zMin=zoffset-fDz;
                                                   >> 858     zMax=zoffset+fDz;
                                                   >> 859 
                                                   >> 860     if ( pVoxelLimit.IsZLimited() )
                                                   >> 861     {
                                                   >> 862       if ( (zMin > pVoxelLimit.GetMaxZExtent() + kCarTolerance)
                                                   >> 863         || (zMax < pVoxelLimit.GetMinZExtent() - kCarTolerance) )
                                                   >> 864       {
                                                   >> 865         return false;
                                                   >> 866       }
                                                   >> 867       else
                                                   >> 868       {
                                                   >> 869         if ( zMin < pVoxelLimit.GetMinZExtent() )
                                                   >> 870         {
                                                   >> 871           zMin = pVoxelLimit.GetMinZExtent() ;
                                                   >> 872         }
                                                   >> 873         if ( zMax > pVoxelLimit.GetMaxZExtent() )
                                                   >> 874         {
                                                   >> 875           zMax = pVoxelLimit.GetMaxZExtent() ;
                                                   >> 876         }
                                                   >> 877       }
                                                   >> 878     }
                                                   >> 879     temp[0] = pt[0].y()+(pt[4].y()-pt[0].y())*(zMin-pt[0].z())
                                                   >> 880                        /(pt[4].z()-pt[0].z()) ;
                                                   >> 881     temp[1] = pt[0].y()+(pt[4].y()-pt[0].y())*(zMax-pt[0].z())
                                                   >> 882                        /(pt[4].z()-pt[0].z()) ;
                                                   >> 883     temp[2] = pt[2].y()+(pt[6].y()-pt[2].y())*(zMin-pt[2].z())
                                                   >> 884                        /(pt[6].z()-pt[2].z()) ;
                                                   >> 885     temp[3] = pt[2].y()+(pt[6].y()-pt[2].y())*(zMax-pt[2].z())
                                                   >> 886                        /(pt[6].z()-pt[2].z()) ;
                                                   >> 887 
                                                   >> 888     yMax = yoffset - std::fabs(fDz*fTthetaSphi) - fDy1 - fDy2 ;
                                                   >> 889     yMin = -yMax ;
                                                   >> 890 
                                                   >> 891     for( i = 0 ; i < 4 ; i++ )
                                                   >> 892     {
                                                   >> 893       if( temp[i] > yMax ) yMax = temp[i] ;
                                                   >> 894       if( temp[i] < yMin ) yMin = temp[i] ;
                                                   >> 895     }    
                                                   >> 896     if ( pVoxelLimit.IsYLimited() )
                                                   >> 897     {
                                                   >> 898       if ( (yMin > pVoxelLimit.GetMaxYExtent() + kCarTolerance)
                                                   >> 899         || (yMax < pVoxelLimit.GetMinYExtent() - kCarTolerance) )
                                                   >> 900       {
                                                   >> 901         return false;
                                                   >> 902       }
                                                   >> 903       else
                                                   >> 904       {
                                                   >> 905         if ( yMin < pVoxelLimit.GetMinYExtent() )
                                                   >> 906         {
                                                   >> 907           yMin = pVoxelLimit.GetMinYExtent() ;
                                                   >> 908         }
                                                   >> 909         if ( yMax > pVoxelLimit.GetMaxYExtent() )
                                                   >> 910         {
                                                   >> 911           yMax = pVoxelLimit.GetMaxYExtent() ;
                                                   >> 912         }
                                                   >> 913       }
                                                   >> 914     }
                                                   >> 915     temp[0] = pt[0].x()+(pt[4].x()-pt[0].x())
                                                   >> 916                        *(zMin-pt[0].z())/(pt[4].z()-pt[0].z()) ;
                                                   >> 917     temp[1] = pt[0].x()+(pt[4].x()-pt[0].x())
                                                   >> 918                        *(zMax-pt[0].z())/(pt[4].z()-pt[0].z()) ;
                                                   >> 919     temp[2] = pt[2].x()+(pt[6].x()-pt[2].x())
                                                   >> 920                        *(zMin-pt[2].z())/(pt[6].z()-pt[2].z()) ;
                                                   >> 921     temp[3] = pt[2].x()+(pt[6].x()-pt[2].x())
                                                   >> 922                        *(zMax-pt[2].z())/(pt[6].z()-pt[2].z()) ;
                                                   >> 923     temp[4] = pt[3].x()+(pt[7].x()-pt[3].x())
                                                   >> 924                        *(zMin-pt[3].z())/(pt[7].z()-pt[3].z()) ;
                                                   >> 925     temp[5] = pt[3].x()+(pt[7].x()-pt[3].x())
                                                   >> 926                        *(zMax-pt[3].z())/(pt[7].z()-pt[3].z()) ;
                                                   >> 927     temp[6] = pt[1].x()+(pt[5].x()-pt[1].x())
                                                   >> 928                        *(zMin-pt[1].z())/(pt[5].z()-pt[1].z()) ;
                                                   >> 929     temp[7] = pt[1].x()+(pt[5].x()-pt[1].x())
                                                   >> 930                        *(zMax-pt[1].z())/(pt[5].z()-pt[1].z()) ;
                                                   >> 931       
                                                   >> 932     xMax = xoffset - std::fabs(fDz*fTthetaCphi) - fDx1 - fDx2 -fDx3 - fDx4 ;
                                                   >> 933     xMin = -xMax ;
                                                   >> 934 
                                                   >> 935     for( i = 0 ; i < 8 ; i++ )
                                                   >> 936     {
                                                   >> 937       if( temp[i] > xMax) xMax = temp[i] ;
                                                   >> 938       if( temp[i] < xMin) xMin = temp[i] ;
                                                   >> 939     }                                            
                                                   >> 940     if (pVoxelLimit.IsXLimited())   // xMax/Min = f(yMax/Min) ?
                                                   >> 941     {
                                                   >> 942       if ( (xMin > pVoxelLimit.GetMaxXExtent() + kCarTolerance)
                                                   >> 943         || (xMax < pVoxelLimit.GetMinXExtent() - kCarTolerance) )
                                                   >> 944       {
                                                   >> 945         return false;
                                                   >> 946       }
                                                   >> 947       else
                                                   >> 948       {
                                                   >> 949         if ( xMin < pVoxelLimit.GetMinXExtent() )
                                                   >> 950         {
                                                   >> 951           xMin = pVoxelLimit.GetMinXExtent() ;
                                                   >> 952         }
                                                   >> 953         if ( xMax > pVoxelLimit.GetMaxXExtent() )
                                                   >> 954         {
                                                   >> 955           xMax = pVoxelLimit.GetMaxXExtent() ;
                                                   >> 956         }
                                                   >> 957       }
                                                   >> 958     }
                                                   >> 959     switch (pAxis)
                                                   >> 960     {
                                                   >> 961       case kXAxis:
                                                   >> 962         pMin=xMin;
                                                   >> 963         pMax=xMax;
                                                   >> 964         break;
                                                   >> 965 
                                                   >> 966       case kYAxis:
                                                   >> 967         pMin=yMin;
                                                   >> 968         pMax=yMax;
                                                   >> 969         break;
                                                   >> 970 
                                                   >> 971       case kZAxis:
                                                   >> 972         pMin=zMin;
                                                   >> 973         pMax=zMax;
                                                   >> 974         break;
                                                   >> 975 
                                                   >> 976       default:
                                                   >> 977         break;
                                                   >> 978     }
                                                   >> 979     pMin -= kCarTolerance;
                                                   >> 980     pMax += kCarTolerance;
                                                   >> 981 
                                                   >> 982     flag = true;
605   }                                               983   }
                                                   >> 984   else    // General rotated case -
                                                   >> 985   {
                                                   >> 986     G4bool existsAfterClip = false ;
                                                   >> 987     G4ThreeVectorList*       vertices;
                                                   >> 988     pMin                   = +kInfinity;
                                                   >> 989     pMax                   = -kInfinity;
                                                   >> 990       
                                                   >> 991     // Calculate rotated vertex coordinates. Operator 'new' is called
                                                   >> 992 
                                                   >> 993     vertices = CreateRotatedVertices(pTransform);
                                                   >> 994       
                                                   >> 995     xMin = +kInfinity; yMin = +kInfinity; zMin = +kInfinity;
                                                   >> 996     xMax = -kInfinity; yMax = -kInfinity; zMax = -kInfinity;
                                                   >> 997       
                                                   >> 998     for( G4int nv = 0 ; nv < 8 ; nv++ )
                                                   >> 999     { 
                                                   >> 1000       if( (*vertices)[nv].x() > xMax ) xMax = (*vertices)[nv].x();
                                                   >> 1001       if( (*vertices)[nv].y() > yMax ) yMax = (*vertices)[nv].y();
                                                   >> 1002       if( (*vertices)[nv].z() > zMax ) zMax = (*vertices)[nv].z();
                                                   >> 1003       
                                                   >> 1004       if( (*vertices)[nv].x() < xMin ) xMin = (*vertices)[nv].x();
                                                   >> 1005       if( (*vertices)[nv].y() < yMin ) yMin = (*vertices)[nv].y();
                                                   >> 1006       if( (*vertices)[nv].z() < zMin ) zMin = (*vertices)[nv].z();
                                                   >> 1007     }
                                                   >> 1008     if ( pVoxelLimit.IsZLimited() )
                                                   >> 1009     {
                                                   >> 1010       if ( (zMin > pVoxelLimit.GetMaxZExtent() + kCarTolerance)
                                                   >> 1011         || (zMax < pVoxelLimit.GetMinZExtent() - kCarTolerance) )
                                                   >> 1012       {
                                                   >> 1013         delete vertices ;    //  'new' in the function called
                                                   >> 1014         return false;
                                                   >> 1015       }
                                                   >> 1016       else
                                                   >> 1017       {
                                                   >> 1018         if ( zMin < pVoxelLimit.GetMinZExtent() )
                                                   >> 1019         {
                                                   >> 1020           zMin = pVoxelLimit.GetMinZExtent() ;
                                                   >> 1021         }
                                                   >> 1022         if ( zMax > pVoxelLimit.GetMaxZExtent() )
                                                   >> 1023         {
                                                   >> 1024           zMax = pVoxelLimit.GetMaxZExtent() ;
                                                   >> 1025         }
                                                   >> 1026       }
                                                   >> 1027     } 
                                                   >> 1028     if ( pVoxelLimit.IsYLimited() )
                                                   >> 1029     {
                                                   >> 1030       if ( (yMin > pVoxelLimit.GetMaxYExtent() + kCarTolerance)
                                                   >> 1031         || (yMax < pVoxelLimit.GetMinYExtent() - kCarTolerance) )
                                                   >> 1032       {
                                                   >> 1033         delete vertices ;    //  'new' in the function called
                                                   >> 1034         return false;
                                                   >> 1035       }
                                                   >> 1036       else
                                                   >> 1037       {
                                                   >> 1038         if ( yMin < pVoxelLimit.GetMinYExtent() )
                                                   >> 1039         {
                                                   >> 1040           yMin = pVoxelLimit.GetMinYExtent() ;
                                                   >> 1041         }
                                                   >> 1042         if ( yMax > pVoxelLimit.GetMaxYExtent() )
                                                   >> 1043         {
                                                   >> 1044           yMax = pVoxelLimit.GetMaxYExtent() ;
                                                   >> 1045         }
                                                   >> 1046       }
                                                   >> 1047     }
                                                   >> 1048     if ( pVoxelLimit.IsXLimited() )
                                                   >> 1049     {
                                                   >> 1050       if ( (xMin > pVoxelLimit.GetMaxXExtent() + kCarTolerance)
                                                   >> 1051         || (xMax < pVoxelLimit.GetMinXExtent() - kCarTolerance) )
                                                   >> 1052       {
                                                   >> 1053         delete vertices ;    //  'new' in the function called
                                                   >> 1054         return false ;
                                                   >> 1055       } 
                                                   >> 1056       else
                                                   >> 1057       {
                                                   >> 1058         if ( xMin < pVoxelLimit.GetMinXExtent() )
                                                   >> 1059         {
                                                   >> 1060           xMin = pVoxelLimit.GetMinXExtent() ;
                                                   >> 1061         }
                                                   >> 1062         if ( xMax > pVoxelLimit.GetMaxXExtent() )
                                                   >> 1063         {
                                                   >> 1064           xMax = pVoxelLimit.GetMaxXExtent() ;
                                                   >> 1065         }
                                                   >> 1066       }
                                                   >> 1067     }
                                                   >> 1068     switch (pAxis)
                                                   >> 1069     {
                                                   >> 1070       case kXAxis:
                                                   >> 1071         pMin=xMin;
                                                   >> 1072         pMax=xMax;
                                                   >> 1073         break;
606                                                   1074 
607   // Set bounding envelope (benv) and calculat << 1075       case kYAxis:
608   //                                           << 1076         pMin=yMin;
609   G4ThreeVector pt[8];                         << 1077         pMax=yMax;
610   GetVertices(pt);                             << 1078         break;
                                                   >> 1079 
                                                   >> 1080       case kZAxis:
                                                   >> 1081         pMin=zMin;
                                                   >> 1082         pMax=zMax;
                                                   >> 1083         break;
611                                                   1084 
612   G4ThreeVectorList baseA(4), baseB(4);        << 1085       default:
613   baseA[0] = pt[0];                            << 1086         break;
614   baseA[1] = pt[1];                            << 1087     }
615   baseA[2] = pt[3];                            << 1088     if ( (pMin != kInfinity) || (pMax != -kInfinity) )
616   baseA[3] = pt[2];                            << 1089     {
617                                                << 1090       existsAfterClip=true;
618   baseB[0] = pt[4];                            << 1091         
619   baseB[1] = pt[5];                            << 1092       // Add tolerance to avoid precision troubles
620   baseB[2] = pt[7];                            << 1093       //
621   baseB[3] = pt[6];                            << 1094       pMin -= kCarTolerance ;
622                                                << 1095       pMax += kCarTolerance ;      
623   std::vector<const G4ThreeVectorList *> polyg << 1096     }
624   polygons[0] = &baseA;                        << 1097     delete vertices ;          //  'new' in the function called
625   polygons[1] = &baseB;                        << 1098     flag = existsAfterClip ;
626                                                << 1099   }
627   G4BoundingEnvelope benv(bmin,bmax,polygons); << 1100   return flag;
628   exist = benv.CalculateExtent(pAxis,pVoxelLim << 
629   return exist;                                << 
630 }                                                 1101 }
631                                                   1102 
632 ////////////////////////////////////////////// << 1103 
                                                   >> 1104 ////////////////////////////////////////////////////////////////////////
633 //                                                1105 //
634 // Return whether point is inside/outside/on_s << 1106 // Return whether point inside/outside/on surface, using tolerance
635                                                   1107 
636 EInside G4Trap::Inside( const G4ThreeVector& p    1108 EInside G4Trap::Inside( const G4ThreeVector& p ) const
637 {                                                 1109 {
638   switch (fTrapType)                           << 1110   EInside in;
                                                   >> 1111   G4double Dist;
                                                   >> 1112   G4int i;
                                                   >> 1113   if ( std::fabs(p.z()) <= fDz-kCarTolerance*0.5)
                                                   >> 1114   {
                                                   >> 1115     in = kInside;
                                                   >> 1116 
                                                   >> 1117     for ( i = 0;i < 4;i++ )
                                                   >> 1118     {
                                                   >> 1119       Dist = fPlanes[i].a*p.x() + fPlanes[i].b*p.y()
                                                   >> 1120             +fPlanes[i].c*p.z() + fPlanes[i].d;
                                                   >> 1121 
                                                   >> 1122       if      (Dist >  kCarTolerance*0.5)  return in = kOutside;
                                                   >> 1123       else if (Dist > -kCarTolerance*0.5)         in = kSurface;
                                                   >> 1124        
                                                   >> 1125     }
                                                   >> 1126   }
                                                   >> 1127   else if (std::fabs(p.z()) <= fDz+kCarTolerance*0.5)
639   {                                               1128   {
640     case 0: // General case                    << 1129     in = kSurface;
                                                   >> 1130 
                                                   >> 1131     for ( i = 0; i < 4; i++ )
641     {                                             1132     {
642       G4double dz = std::abs(p.z())-fDz;       << 1133       Dist =  fPlanes[i].a*p.x() + fPlanes[i].b*p.y()
643       G4double dy1 = fPlanes[0].b*p.y()+fPlane << 1134              +fPlanes[i].c*p.z() + fPlanes[i].d;
644       G4double dy2 = fPlanes[1].b*p.y()+fPlane << 
645       G4double dy = std::max(dz,std::max(dy1,d << 
646                                                << 
647       G4double dx1 = fPlanes[2].a*p.x()+fPlane << 
648                    + fPlanes[2].c*p.z()+fPlane << 
649       G4double dx2 = fPlanes[3].a*p.x()+fPlane << 
650                    + fPlanes[3].c*p.z()+fPlane << 
651       G4double dist = std::max(dy,std::max(dx1 << 
652                                                << 
653       return (dist > halfCarTolerance) ? kOuts << 
654         ((dist > -halfCarTolerance) ? kSurface << 
655     }                                          << 
656     case 1: // YZ section is a rectangle       << 
657     {                                          << 
658       G4double dz = std::abs(p.z())-fDz;       << 
659       G4double dy = std::max(dz,std::abs(p.y() << 
660       G4double dx1 = fPlanes[2].a*p.x()+fPlane << 
661                    + fPlanes[2].c*p.z()+fPlane << 
662       G4double dx2 = fPlanes[3].a*p.x()+fPlane << 
663                    + fPlanes[3].c*p.z()+fPlane << 
664       G4double dist = std::max(dy,std::max(dx1 << 
665                                                << 
666       return (dist > halfCarTolerance) ? kOuts << 
667         ((dist > -halfCarTolerance) ? kSurface << 
668     }                                          << 
669     case 2: // YZ section is a rectangle and   << 
670     {       // XZ section is an isosceles trap << 
671       G4double dz = std::abs(p.z())-fDz;       << 
672       G4double dy = std::max(dz,std::abs(p.y() << 
673       G4double dx = fPlanes[3].a*std::abs(p.x( << 
674                   + fPlanes[3].c*p.z()+fPlanes << 
675       G4double dist = std::max(dy,dx);         << 
676                                                << 
677       return (dist > halfCarTolerance) ? kOuts << 
678         ((dist > -halfCarTolerance) ? kSurface << 
679     }                                          << 
680     case 3: // YZ section is a rectangle and   << 
681     {       // XY section is an isosceles trap << 
682       G4double dz = std::abs(p.z())-fDz;       << 
683       G4double dy = std::max(dz,std::abs(p.y() << 
684       G4double dx = fPlanes[3].a*std::abs(p.x( << 
685                   + fPlanes[3].b*p.y()+fPlanes << 
686       G4double dist = std::max(dy,dx);         << 
687                                                   1135 
688       return (dist > halfCarTolerance) ? kOuts << 1136       if (Dist > kCarTolerance*0.5)        return in = kOutside;      
689         ((dist > -halfCarTolerance) ? kSurface << 
690     }                                             1137     }
691   }                                               1138   }
692   return kOutside;                             << 1139   else  in = kOutside;
                                                   >> 1140   
                                                   >> 1141   return in;
693 }                                                 1142 }
694                                                   1143 
695 ////////////////////////////////////////////// << 1144 /////////////////////////////////////////////////////////////////////////////
696 //                                                1145 //
697 // Determine side, and return corresponding no << 1146 // Calculate side nearest to p, and return normal
                                                   >> 1147 // If 2+ sides equidistant, first side's normal returned (arbitrarily)
698                                                   1148 
699 G4ThreeVector G4Trap::SurfaceNormal( const G4T    1149 G4ThreeVector G4Trap::SurfaceNormal( const G4ThreeVector& p ) const
700 {                                                 1150 {
701   G4double nx = 0, ny = 0, nz = 0;             << 1151   G4int i, noSurfaces = 0;
702   G4double dz = std::abs(p.z()) - fDz;         << 1152   G4double dist, distz, distx, disty, distmx, distmy, safe = kInfinity;
703   nz = std::copysign(G4double(std::abs(dz) <=  << 1153   G4double delta    = 0.5*kCarTolerance;
                                                   >> 1154   G4ThreeVector norm, sumnorm(0.,0.,0.);
704                                                   1155 
705   switch (fTrapType)                           << 1156   for (i = 0; i < 4; i++)
706   {                                               1157   {
707     case 0: // General case                    << 1158     dist =  std::fabs(fPlanes[i].a*p.x() + fPlanes[i].b*p.y()
708     {                                          << 1159           + fPlanes[i].c*p.z() + fPlanes[i].d);
709       for (G4int i=0; i<2; ++i)                << 1160     if ( dist < safe )
710       {                                        << 
711         G4double dy = fPlanes[i].b*p.y() + fPl << 
712         if (std::abs(dy) > halfCarTolerance) c << 
713         ny  = fPlanes[i].b;                    << 
714         nz += fPlanes[i].c;                    << 
715         break;                                 << 
716       }                                        << 
717       for (G4int i=2; i<4; ++i)                << 
718       {                                        << 
719         G4double dx = fPlanes[i].a*p.x() +     << 
720                       fPlanes[i].b*p.y() + fPl << 
721         if (std::abs(dx) > halfCarTolerance) c << 
722         nx  = fPlanes[i].a;                    << 
723         ny += fPlanes[i].b;                    << 
724         nz += fPlanes[i].c;                    << 
725         break;                                 << 
726       }                                        << 
727       break;                                   << 
728     }                                          << 
729     case 1: // YZ section - rectangle          << 
730     {                                             1161     {
731       G4double dy = std::abs(p.y()) + fPlanes[ << 1162       safe = dist;
732       ny = std::copysign(G4double(std::abs(dy) << 
733       for (G4int i=2; i<4; ++i)                << 
734       {                                        << 
735         G4double dx = fPlanes[i].a*p.x() +     << 
736                       fPlanes[i].b*p.y() + fPl << 
737         if (std::abs(dx) > halfCarTolerance) c << 
738         nx  = fPlanes[i].a;                    << 
739         ny += fPlanes[i].b;                    << 
740         nz += fPlanes[i].c;                    << 
741         break;                                 << 
742       }                                        << 
743       break;                                   << 
744     }                                          << 
745     case 2: // YZ section - rectangle, XZ sect << 
746     {                                          << 
747       G4double dy = std::abs(p.y()) + fPlanes[ << 
748       ny = std::copysign(G4double(std::abs(dy) << 
749       G4double dx = fPlanes[3].a*std::abs(p.x( << 
750                     fPlanes[3].c*p.z() + fPlan << 
751       G4double k = std::abs(dx) <= halfCarTole << 
752       nx  = std::copysign(k, p.x())*fPlanes[3] << 
753       nz += k*fPlanes[3].c;                    << 
754       break;                                   << 
755     }                                          << 
756     case 3: // YZ section - rectangle, XY sect << 
757     {                                          << 
758       G4double dy = std::abs(p.y()) + fPlanes[ << 
759       ny = std::copysign(G4double(std::abs(dy) << 
760       G4double dx = fPlanes[3].a*std::abs(p.x( << 
761                     fPlanes[3].b*p.y() + fPlan << 
762       G4double k = std::abs(dx) <= halfCarTole << 
763       nx  = std::copysign(k, p.x())*fPlanes[3] << 
764       ny += k*fPlanes[3].b;                    << 
765       break;                                   << 
766     }                                             1163     }
767   }                                               1164   }
                                                   >> 1165   distz  = std::fabs( std::fabs( p.z() ) - fDz );
768                                                   1166 
769   // Return normal                             << 1167   distmy = std::fabs( fPlanes[0].a*p.x() + fPlanes[0].b*p.y()
770   //                                           << 1168                     + fPlanes[0].c*p.z() + fPlanes[0].d      );
771   G4double mag2 = nx*nx + ny*ny + nz*nz;       << 1169 
772   if (mag2 == 1)      return { nx,ny,nz };     << 1170   disty  = std::fabs( fPlanes[1].a*p.x() + fPlanes[1].b*p.y()
773   else if (mag2 != 0) return G4ThreeVector(nx, << 1171                     + fPlanes[1].c*p.z() + fPlanes[1].d      );
774   else                                         << 1172 
                                                   >> 1173   distmx = std::fabs( fPlanes[2].a*p.x() + fPlanes[2].b*p.y()
                                                   >> 1174                     + fPlanes[2].c*p.z() + fPlanes[2].d      );
                                                   >> 1175 
                                                   >> 1176   distx  = std::fabs( fPlanes[3].a*p.x() + fPlanes[3].b*p.y()
                                                   >> 1177                     + fPlanes[3].c*p.z() + fPlanes[3].d      );
                                                   >> 1178 
                                                   >> 1179   G4ThreeVector nX  = G4ThreeVector(fPlanes[3].a,fPlanes[3].b,fPlanes[3].c);
                                                   >> 1180   G4ThreeVector nmX = G4ThreeVector(fPlanes[2].a,fPlanes[2].b,fPlanes[2].c);
                                                   >> 1181   G4ThreeVector nY  = G4ThreeVector(fPlanes[1].a,fPlanes[1].b,fPlanes[1].c);
                                                   >> 1182   G4ThreeVector nmY = G4ThreeVector(fPlanes[0].a,fPlanes[0].b,fPlanes[0].c);
                                                   >> 1183   G4ThreeVector nZ  = G4ThreeVector(0.,0.,1.0);
                                                   >> 1184 
                                                   >> 1185   if (distx <= delta)      
                                                   >> 1186   {
                                                   >> 1187     noSurfaces ++;
                                                   >> 1188     sumnorm += nX;     
                                                   >> 1189   }
                                                   >> 1190   if (distmx <= delta)      
                                                   >> 1191   {
                                                   >> 1192     noSurfaces ++;
                                                   >> 1193     sumnorm += nmX;      
                                                   >> 1194   }
                                                   >> 1195   if (disty <= delta)
                                                   >> 1196   {
                                                   >> 1197     noSurfaces ++;
                                                   >> 1198     sumnorm += nY;  
                                                   >> 1199   }
                                                   >> 1200   if (distmy <= delta)
                                                   >> 1201   {
                                                   >> 1202     noSurfaces ++;
                                                   >> 1203     sumnorm += nmY;  
                                                   >> 1204   }
                                                   >> 1205   if (distz <= delta)  
                                                   >> 1206   {
                                                   >> 1207     noSurfaces ++;
                                                   >> 1208     if ( p.z() >= 0.)  sumnorm += nZ;
                                                   >> 1209     else               sumnorm -= nZ; 
                                                   >> 1210   }
                                                   >> 1211   if ( noSurfaces == 0 )
775   {                                               1212   {
776     // Point is not on the surface             << 
777     //                                         << 
778 #ifdef G4CSGDEBUG                                 1213 #ifdef G4CSGDEBUG
779     std::ostringstream message;                << 
780     G4long oldprc = message.precision(16);     << 
781     message << "Point p is not on surface (!?) << 
782             << GetName() << G4endl;            << 
783     message << "Position:\n";                  << 
784     message << "   p.x() = " << p.x()/mm << "  << 
785     message << "   p.y() = " << p.y()/mm << "  << 
786     message << "   p.z() = " << p.z()/mm << "  << 
787     G4cout.precision(oldprc) ;                 << 
788     G4Exception("G4Trap::SurfaceNormal(p)", "G    1214     G4Exception("G4Trap::SurfaceNormal(p)", "GeomSolids1002",
789                 JustWarning, message );        << 1215                 JustWarning, "Point p is not on surface !?" );
790     DumpInfo();                                << 1216 #endif 
791 #endif                                         << 1217      norm = ApproxSurfaceNormal(p);
792     return ApproxSurfaceNormal(p);             << 
793   }                                               1218   }
                                                   >> 1219   else if ( noSurfaces == 1 ) norm = sumnorm;
                                                   >> 1220   else                        norm = sumnorm.unit();
                                                   >> 1221   return norm;
794 }                                                 1222 }
795                                                   1223 
796 ////////////////////////////////////////////// << 1224 ////////////////////////////////////////////////////////////////////////////////////
797 //                                                1225 //
798 // Algorithm for SurfaceNormal() following the    1226 // Algorithm for SurfaceNormal() following the original specification
799 // for points not on the surface                  1227 // for points not on the surface
800                                                   1228 
801 G4ThreeVector G4Trap::ApproxSurfaceNormal( con    1229 G4ThreeVector G4Trap::ApproxSurfaceNormal( const G4ThreeVector& p ) const
802 {                                                 1230 {
803   G4double dist = -DBL_MAX;                    << 1231   G4double safe=kInfinity,Dist,safez;
804   G4int iside = 0;                             << 1232   G4int i,imin=0;
805   for (G4int i=0; i<4; ++i)                    << 1233   for (i=0;i<4;i++)
806   {                                            << 1234   {
807     G4double d = fPlanes[i].a*p.x() +          << 1235     Dist=std::fabs(fPlanes[i].a*p.x()+fPlanes[i].b*p.y()
808                  fPlanes[i].b*p.y() +          << 1236         +fPlanes[i].c*p.z()+fPlanes[i].d);
809                  fPlanes[i].c*p.z() + fPlanes[ << 1237     if (Dist<safe)
810     if (d > dist) { dist = d; iside = i; }     << 1238     {
                                                   >> 1239       safe=Dist;
                                                   >> 1240       imin=i;
                                                   >> 1241     }
                                                   >> 1242   }
                                                   >> 1243   safez=std::fabs(std::fabs(p.z())-fDz);
                                                   >> 1244   if (safe<safez)
                                                   >> 1245   {
                                                   >> 1246     return G4ThreeVector(fPlanes[imin].a,fPlanes[imin].b,fPlanes[imin].c);
811   }                                               1247   }
812                                                << 
813   G4double distz = std::abs(p.z()) - fDz;      << 
814   if (dist > distz)                            << 
815     return { fPlanes[iside].a, fPlanes[iside]. << 
816   else                                            1248   else
817     return { 0, 0, (G4double)((p.z() < 0) ? -1 << 1249   {
                                                   >> 1250     if (p.z()>0)
                                                   >> 1251     {
                                                   >> 1252       return G4ThreeVector(0,0,1);
                                                   >> 1253     }
                                                   >> 1254     else
                                                   >> 1255     {
                                                   >> 1256       return G4ThreeVector(0,0,-1);
                                                   >> 1257     }
                                                   >> 1258   }
818 }                                                 1259 }
819                                                   1260 
820 ////////////////////////////////////////////// << 1261 ////////////////////////////////////////////////////////////////////////////
                                                   >> 1262 //
                                                   >> 1263 // Calculate distance to shape from outside - return kInfinity if no intersection
821 //                                                1264 //
822 // Calculate distance to shape from outside    << 1265 // ALGORITHM:
823 //  - return kInfinity if no intersection      << 1266 // For each component, calculate pair of minimum and maximum intersection
                                                   >> 1267 // values for which the particle is in the extent of the shape
                                                   >> 1268 // - The smallest (MAX minimum) allowed distance of the pairs is intersect
824                                                   1269 
825 G4double G4Trap::DistanceToIn(const G4ThreeVec << 1270 G4double G4Trap::DistanceToIn( const G4ThreeVector& p,
826                               const G4ThreeVec << 1271                                const G4ThreeVector& v ) const
827 {                                                 1272 {
828   // Z intersections                           << 
829   //                                           << 
830   if ((std::abs(p.z()) - fDz) >= -halfCarToler << 
831     return kInfinity;                          << 
832   G4double invz = (-v.z() == 0) ? DBL_MAX : -1 << 
833   G4double dz = (invz < 0) ? fDz : -fDz;       << 
834   G4double tzmin = (p.z() + dz)*invz;          << 
835   G4double tzmax = (p.z() - dz)*invz;          << 
836                                                   1273 
837   // Y intersections                           << 1274   G4double snxt;    // snxt = default return value
                                                   >> 1275   G4double max,smax,smin;
                                                   >> 1276   G4double pdist,Comp,vdist;
                                                   >> 1277   G4int i;
                                                   >> 1278   //
                                                   >> 1279   // Z Intersection range
838   //                                              1280   //
839   G4double tymin = 0, tymax = DBL_MAX;         << 1281   if ( v.z() > 0 )
840   G4int i = 0;                                 << 
841   for ( ; i<2; ++i)                            << 
842   {                                               1282   {
843     G4double cosa = fPlanes[i].b*v.y() + fPlan << 1283     max = fDz - p.z() ;
844     G4double dist = fPlanes[i].b*p.y() + fPlan << 1284     if (max > 0.5*kCarTolerance)
845     if (dist >= -halfCarTolerance)             << 
846     {                                             1285     {
847       if (cosa >= 0) return kInfinity;         << 1286       smax = max/v.z();
848       G4double tmp  = -dist/cosa;              << 1287       smin = (-fDz-p.z())/v.z();
849       if (tymin < tmp) tymin = tmp;            << 
850     }                                             1288     }
851     else if (cosa > 0)                         << 1289     else
852     {                                             1290     {
853       G4double tmp  = -dist/cosa;              << 1291       return snxt=kInfinity;
854       if (tymax > tmp) tymax = tmp;            << 
855     }                                             1292     }
856   }                                               1293   }
857                                                << 1294   else if (v.z() < 0 )
858   // Z intersections                           << 1295   {
859   //                                           << 1296     max = - fDz - p.z() ;
860   G4double txmin = 0, txmax = DBL_MAX;         << 1297     if (max < -0.5*kCarTolerance )
861   for ( ; i<4; ++i)                            << 1298     {
                                                   >> 1299       smax=max/v.z();
                                                   >> 1300       smin=(fDz-p.z())/v.z();
                                                   >> 1301     }
                                                   >> 1302     else
                                                   >> 1303     {
                                                   >> 1304       return snxt=kInfinity;
                                                   >> 1305     }
                                                   >> 1306   }
                                                   >> 1307   else
862   {                                               1308   {
863     G4double cosa = fPlanes[i].a*v.x()+fPlanes << 1309     if (std::fabs(p.z())<fDz - 0.5*kCarTolerance) // Inside was <=fDz
864     G4double dist = fPlanes[i].a*p.x()+fPlanes << 
865                     fPlanes[i].d;              << 
866     if (dist >= -halfCarTolerance)             << 
867     {                                             1310     {
868       if (cosa >= 0) return kInfinity;         << 1311       smin=0;
869       G4double tmp  = -dist/cosa;              << 1312       smax=kInfinity;
870       if (txmin < tmp) txmin = tmp;            << 
871     }                                             1313     }
872     else if (cosa > 0)                         << 1314     else
873     {                                             1315     {
874       G4double tmp  = -dist/cosa;              << 1316       return snxt=kInfinity;
875       if (txmax > tmp) txmax = tmp;            << 
876     }                                             1317     }
877   }                                               1318   }
878                                                   1319 
879   // Find distance                             << 1320   for (i=0;i<4;i++)
                                                   >> 1321   {
                                                   >> 1322     pdist=fPlanes[i].a*p.x()+fPlanes[i].b*p.y()
                                                   >> 1323          +fPlanes[i].c*p.z()+fPlanes[i].d;
                                                   >> 1324     Comp=fPlanes[i].a*v.x()+fPlanes[i].b*v.y()+fPlanes[i].c*v.z();
                                                   >> 1325     if ( pdist >= -0.5*kCarTolerance )      // was >0
                                                   >> 1326     {
                                                   >> 1327       //
                                                   >> 1328       // Outside the plane -> this is an extent entry distance
                                                   >> 1329       //
                                                   >> 1330       if (Comp >= 0)   // was >0
                                                   >> 1331       {
                                                   >> 1332         return snxt=kInfinity ;
                                                   >> 1333       }
                                                   >> 1334       else 
                                                   >> 1335       {
                                                   >> 1336         vdist=-pdist/Comp;
                                                   >> 1337         if (vdist>smin)
                                                   >> 1338         {
                                                   >> 1339           if (vdist<smax)
                                                   >> 1340           {
                                                   >> 1341             smin = vdist;
                                                   >> 1342           }
                                                   >> 1343           else
                                                   >> 1344           {
                                                   >> 1345             return snxt=kInfinity;
                                                   >> 1346           }
                                                   >> 1347         }
                                                   >> 1348       }
                                                   >> 1349     }
                                                   >> 1350     else
                                                   >> 1351     {
                                                   >> 1352       //
                                                   >> 1353       // Inside the plane -> couble  be an extent exit distance (smax)
                                                   >> 1354       //
                                                   >> 1355       if (Comp>0)  // Will leave extent
                                                   >> 1356       {
                                                   >> 1357         vdist=-pdist/Comp;
                                                   >> 1358         if (vdist<smax)
                                                   >> 1359         {
                                                   >> 1360           if (vdist>smin)
                                                   >> 1361           {
                                                   >> 1362             smax=vdist;
                                                   >> 1363           }
                                                   >> 1364           else
                                                   >> 1365           {
                                                   >> 1366             return snxt=kInfinity;
                                                   >> 1367           }
                                                   >> 1368         }  
                                                   >> 1369       }
                                                   >> 1370     }
                                                   >> 1371   }
880   //                                              1372   //
881   G4double tmin = std::max(std::max(txmin,tymi << 1373   // Checks in non z plane intersections ensure smin<smax
882   G4double tmax = std::min(std::min(txmax,tyma << 1374   //
883                                                << 1375   if (smin >=0 )
884   if (tmax <= tmin + halfCarTolerance) return  << 1376   {
885   return (tmin < halfCarTolerance ) ? 0. : tmi << 1377     snxt = smin ;
                                                   >> 1378   }
                                                   >> 1379   else
                                                   >> 1380   {
                                                   >> 1381     snxt = 0 ;
                                                   >> 1382   }
                                                   >> 1383   return snxt;
886 }                                                 1384 }
887                                                   1385 
888 ////////////////////////////////////////////// << 1386 ///////////////////////////////////////////////////////////////////////////
889 //                                                1387 //
890 // Calculate exact shortest distance to any bo    1388 // Calculate exact shortest distance to any boundary from outside
891 // This is the best fast estimation of the sho    1389 // This is the best fast estimation of the shortest distance to trap
892 // - return 0 if point is inside               << 1390 // - Returns 0 is ThreeVector inside
893                                                   1391 
894 G4double G4Trap::DistanceToIn( const G4ThreeVe    1392 G4double G4Trap::DistanceToIn( const G4ThreeVector& p ) const
895 {                                                 1393 {
896   switch (fTrapType)                           << 1394   G4double safe=0.0,Dist;
                                                   >> 1395   G4int i;
                                                   >> 1396   safe=std::fabs(p.z())-fDz;
                                                   >> 1397   for (i=0;i<4;i++)
897   {                                               1398   {
898     case 0: // General case                    << 1399     Dist=fPlanes[i].a*p.x()+fPlanes[i].b*p.y()
899     {                                          << 1400         +fPlanes[i].c*p.z()+fPlanes[i].d;
900       G4double dz = std::abs(p.z())-fDz;       << 1401     if (Dist > safe) safe=Dist;
901       G4double dy1 = fPlanes[0].b*p.y()+fPlane << 
902       G4double dy2 = fPlanes[1].b*p.y()+fPlane << 
903       G4double dy = std::max(dz,std::max(dy1,d << 
904                                                << 
905       G4double dx1 = fPlanes[2].a*p.x()+fPlane << 
906                    + fPlanes[2].c*p.z()+fPlane << 
907       G4double dx2 = fPlanes[3].a*p.x()+fPlane << 
908                    + fPlanes[3].c*p.z()+fPlane << 
909       G4double dist = std::max(dy,std::max(dx1 << 
910       return (dist > 0) ? dist : 0.;           << 
911     }                                          << 
912     case 1: // YZ section is a rectangle       << 
913     {                                          << 
914       G4double dz = std::abs(p.z())-fDz;       << 
915       G4double dy = std::max(dz,std::abs(p.y() << 
916       G4double dx1 = fPlanes[2].a*p.x()+fPlane << 
917                    + fPlanes[2].c*p.z()+fPlane << 
918       G4double dx2 = fPlanes[3].a*p.x()+fPlane << 
919                    + fPlanes[3].c*p.z()+fPlane << 
920       G4double dist = std::max(dy,std::max(dx1 << 
921       return (dist > 0) ? dist : 0.;           << 
922     }                                          << 
923     case 2: // YZ section is a rectangle and   << 
924     {       // XZ section is an isosceles trap << 
925       G4double dz = std::abs(p.z())-fDz;       << 
926       G4double dy = std::max(dz,std::abs(p.y() << 
927       G4double dx = fPlanes[3].a*std::abs(p.x( << 
928                   + fPlanes[3].c*p.z()+fPlanes << 
929       G4double dist = std::max(dy,dx);         << 
930       return (dist > 0) ? dist : 0.;           << 
931     }                                          << 
932     case 3: // YZ section is a rectangle and   << 
933     {       // XY section is an isosceles trap << 
934       G4double dz = std::abs(p.z())-fDz;       << 
935       G4double dy = std::max(dz,std::abs(p.y() << 
936       G4double dx = fPlanes[3].a*std::abs(p.x( << 
937                   + fPlanes[3].b*p.y()+fPlanes << 
938       G4double dist = std::max(dy,dx);         << 
939       return (dist > 0) ? dist : 0.;           << 
940     }                                          << 
941   }                                               1402   }
942   return 0.;                                   << 1403   if (safe<0) safe=0;
                                                   >> 1404   return safe;  
943 }                                                 1405 }
944                                                   1406 
945 ////////////////////////////////////////////// << 1407 /////////////////////////////////////////////////////////////////////////////////
946 //                                                1408 //
947 // Calculate distance to surface of shape from << 1409 // Calculate distance to surface of shape from inside
948 // find normal at exit point, if required      << 1410 // Calculate distance to x/y/z planes - smallest is exiting distance
949 // - when leaving the surface, return 0        << 
950                                                   1411 
951 G4double G4Trap::DistanceToOut(const G4ThreeVe    1412 G4double G4Trap::DistanceToOut(const G4ThreeVector& p, const G4ThreeVector& v,
952                                const G4bool ca    1413                                const G4bool calcNorm,
953                                      G4bool* v << 1414                                      G4bool *validNorm, G4ThreeVector *n) const
954 {                                                 1415 {
955   // Z intersections                           << 1416   Eside side = kUndef;
                                                   >> 1417   G4double snxt;    // snxt = return value
                                                   >> 1418   G4double pdist,Comp,vdist,max;
                                                   >> 1419   //
                                                   >> 1420   // Z Intersections
956   //                                              1421   //
957   if ((std::abs(p.z()) - fDz) >= -halfCarToler << 1422   if (v.z()>0)
                                                   >> 1423   {
                                                   >> 1424     max=fDz-p.z();
                                                   >> 1425     if (max>kCarTolerance/2)
                                                   >> 1426     {
                                                   >> 1427       snxt=max/v.z();
                                                   >> 1428       side=kPZ;
                                                   >> 1429     }
                                                   >> 1430     else
                                                   >> 1431     {
                                                   >> 1432       if (calcNorm)
                                                   >> 1433       {
                                                   >> 1434         *validNorm=true;
                                                   >> 1435         *n=G4ThreeVector(0,0,1);
                                                   >> 1436       }
                                                   >> 1437       return snxt=0;
                                                   >> 1438     }
                                                   >> 1439   }
                                                   >> 1440   else if (v.z()<0)
958   {                                               1441   {
959     if (calcNorm)                              << 1442     max=-fDz-p.z();
                                                   >> 1443     if (max<-kCarTolerance/2)
960     {                                             1444     {
961       *validNorm = true;                       << 1445       snxt=max/v.z();
962       n->set(0, 0, (p.z() < 0) ? -1 : 1);      << 1446       side=kMZ;
                                                   >> 1447     }
                                                   >> 1448     else
                                                   >> 1449     {
                                                   >> 1450       if (calcNorm)
                                                   >> 1451       {
                                                   >> 1452         *validNorm=true;
                                                   >> 1453         *n=G4ThreeVector(0,0,-1);
                                                   >> 1454       }
                                                   >> 1455       return snxt=0;
963     }                                             1456     }
964     return 0;                                  << 
965   }                                               1457   }
966   G4double vz = v.z();                         << 1458   else
967   G4double tmax = (vz == 0) ? DBL_MAX : (std:: << 1459   {
968   G4int iside = (vz < 0) ? -4 : -2; // little  << 1460     snxt=kInfinity;
                                                   >> 1461   }
969                                                   1462 
970   // Y intersections                           << 
971   //                                              1463   //
972   G4int i = 0;                                 << 1464   // Intersections with planes[0] (expanded because of setting enum)
973   for ( ; i<2; ++i)                            << 1465   //
                                                   >> 1466   pdist=fPlanes[0].a*p.x()+fPlanes[0].b*p.y()+fPlanes[0].c*p.z()+fPlanes[0].d;
                                                   >> 1467   Comp=fPlanes[0].a*v.x()+fPlanes[0].b*v.y()+fPlanes[0].c*v.z();
                                                   >> 1468   if (pdist>0)
974   {                                               1469   {
975     G4double cosa = fPlanes[i].b*v.y() + fPlan << 1470     // Outside the plane
976     if (cosa > 0)                              << 1471     if (Comp>0)
977     {                                             1472     {
978       G4double dist = fPlanes[i].b*p.y() + fPl << 1473       // Leaving immediately
979       if (dist >= -halfCarTolerance)           << 1474       if (calcNorm)
980       {                                           1475       {
981         if (calcNorm)                          << 1476         *validNorm=true;
982         {                                      << 1477         *n=G4ThreeVector(fPlanes[0].a,fPlanes[0].b,fPlanes[0].c);
983           *validNorm = true;                   << 1478       }
984           n->set(0, fPlanes[i].b, fPlanes[i].c << 1479       return snxt=0;
985         }                                      << 1480     }
986         return 0;                              << 1481   }
                                                   >> 1482   else if (pdist<-kCarTolerance/2)
                                                   >> 1483   {
                                                   >> 1484     // Inside the plane
                                                   >> 1485     if (Comp>0)
                                                   >> 1486     {
                                                   >> 1487       // Will leave extent
                                                   >> 1488       vdist=-pdist/Comp;
                                                   >> 1489       if (vdist<snxt)
                                                   >> 1490       {
                                                   >> 1491         snxt=vdist;
                                                   >> 1492         side=ks0;
                                                   >> 1493       }
                                                   >> 1494     }
                                                   >> 1495   }
                                                   >> 1496   else
                                                   >> 1497   {
                                                   >> 1498     // On surface
                                                   >> 1499     if (Comp>0)
                                                   >> 1500     {
                                                   >> 1501       if (calcNorm)
                                                   >> 1502       {
                                                   >> 1503         *validNorm=true;
                                                   >> 1504         *n=G4ThreeVector(fPlanes[0].a,fPlanes[0].b,fPlanes[0].c);
987       }                                           1505       }
988       G4double tmp = -dist/cosa;               << 1506       return snxt=0;
989       if (tmax > tmp) { tmax = tmp; iside = i; << 
990     }                                             1507     }
991   }                                               1508   }
992                                                   1509 
993   // X intersections                           << 
994   //                                              1510   //
995   for ( ; i<4; ++i)                            << 1511   // Intersections with planes[1] (expanded because of setting enum)
                                                   >> 1512   //
                                                   >> 1513   pdist=fPlanes[1].a*p.x()+fPlanes[1].b*p.y()+fPlanes[1].c*p.z()+fPlanes[1].d;
                                                   >> 1514   Comp=fPlanes[1].a*v.x()+fPlanes[1].b*v.y()+fPlanes[1].c*v.z();
                                                   >> 1515   if (pdist>0)
996   {                                               1516   {
997     G4double cosa = fPlanes[i].a*v.x()+fPlanes << 1517     // Outside the plane
998     if (cosa > 0)                              << 1518     if (Comp>0)
999     {                                             1519     {
1000       G4double dist = fPlanes[i].a*p.x() +    << 1520       // Leaving immediately
1001                       fPlanes[i].b*p.y() + fP << 1521       if (calcNorm)
1002       if (dist >= -halfCarTolerance)          << 
1003       {                                          1522       {
1004         if (calcNorm)                         << 1523         *validNorm=true;
1005         {                                     << 1524         *n=G4ThreeVector(fPlanes[1].a,fPlanes[1].b,fPlanes[1].c);
1006            *validNorm = true;                 << 1525       }
1007            n->set(fPlanes[i].a, fPlanes[i].b, << 1526       return snxt=0;
1008         }                                     << 1527     }
1009         return 0;                             << 1528   }
                                                   >> 1529   else if (pdist<-kCarTolerance/2)
                                                   >> 1530   {
                                                   >> 1531     // Inside the plane
                                                   >> 1532     if (Comp>0)
                                                   >> 1533     {
                                                   >> 1534       // Will leave extent
                                                   >> 1535       vdist=-pdist/Comp;
                                                   >> 1536       if (vdist<snxt)
                                                   >> 1537       {
                                                   >> 1538         snxt=vdist;
                                                   >> 1539         side=ks1;
                                                   >> 1540       }
                                                   >> 1541     }
                                                   >> 1542   }
                                                   >> 1543   else
                                                   >> 1544   {
                                                   >> 1545     // On surface
                                                   >> 1546     if (Comp>0)
                                                   >> 1547     {
                                                   >> 1548       if (calcNorm)
                                                   >> 1549       {
                                                   >> 1550         *validNorm=true;
                                                   >> 1551         *n=G4ThreeVector(fPlanes[1].a,fPlanes[1].b,fPlanes[1].c);
                                                   >> 1552       }
                                                   >> 1553       return snxt=0;
                                                   >> 1554     }
                                                   >> 1555   }
                                                   >> 1556 
                                                   >> 1557   //
                                                   >> 1558   // Intersections with planes[2] (expanded because of setting enum)
                                                   >> 1559   //
                                                   >> 1560   pdist=fPlanes[2].a*p.x()+fPlanes[2].b*p.y()+fPlanes[2].c*p.z()+fPlanes[2].d;
                                                   >> 1561   Comp=fPlanes[2].a*v.x()+fPlanes[2].b*v.y()+fPlanes[2].c*v.z();
                                                   >> 1562   if (pdist>0)
                                                   >> 1563   {
                                                   >> 1564     // Outside the plane
                                                   >> 1565     if (Comp>0)
                                                   >> 1566     {
                                                   >> 1567       // Leaving immediately
                                                   >> 1568       if (calcNorm)
                                                   >> 1569       {
                                                   >> 1570         *validNorm=true;
                                                   >> 1571         *n=G4ThreeVector(fPlanes[2].a,fPlanes[2].b,fPlanes[2].c);
                                                   >> 1572       }
                                                   >> 1573       return snxt=0;
                                                   >> 1574     }
                                                   >> 1575   }
                                                   >> 1576   else if (pdist<-kCarTolerance/2)
                                                   >> 1577   {
                                                   >> 1578     // Inside the plane
                                                   >> 1579     if (Comp>0)
                                                   >> 1580     {
                                                   >> 1581       // Will leave extent
                                                   >> 1582       vdist=-pdist/Comp;
                                                   >> 1583       if (vdist<snxt)
                                                   >> 1584       {
                                                   >> 1585         snxt=vdist;
                                                   >> 1586         side=ks2;
1010       }                                          1587       }
1011       G4double tmp = -dist/cosa;              << 1588     }
1012       if (tmax > tmp) { tmax = tmp; iside = i << 1589   }
                                                   >> 1590   else
                                                   >> 1591   {
                                                   >> 1592     // On surface
                                                   >> 1593     if (Comp>0)
                                                   >> 1594     {
                                                   >> 1595       if (calcNorm)
                                                   >> 1596       {
                                                   >> 1597         *validNorm=true;
                                                   >> 1598         *n=G4ThreeVector(fPlanes[2].a,fPlanes[2].b,fPlanes[2].c);
                                                   >> 1599       }
                                                   >> 1600       return snxt=0;
1013     }                                            1601     }
1014   }                                              1602   }
1015                                                  1603 
1016   // Set normal, if required, and return dist << 
1017   //                                             1604   //
                                                   >> 1605   // Intersections with planes[3] (expanded because of setting enum)
                                                   >> 1606   //
                                                   >> 1607   pdist=fPlanes[3].a*p.x()+fPlanes[3].b*p.y()+fPlanes[3].c*p.z()+fPlanes[3].d;
                                                   >> 1608   Comp=fPlanes[3].a*v.x()+fPlanes[3].b*v.y()+fPlanes[3].c*v.z();
                                                   >> 1609   if (pdist>0)
                                                   >> 1610   {
                                                   >> 1611     // Outside the plane
                                                   >> 1612     if (Comp>0)
                                                   >> 1613     {
                                                   >> 1614       // Leaving immediately
                                                   >> 1615       if (calcNorm)
                                                   >> 1616       {
                                                   >> 1617         *validNorm=true;
                                                   >> 1618         *n=G4ThreeVector(fPlanes[3].a,fPlanes[3].b,fPlanes[3].c);
                                                   >> 1619       }
                                                   >> 1620       return snxt=0;
                                                   >> 1621     }
                                                   >> 1622   }
                                                   >> 1623   else if (pdist<-kCarTolerance/2)
                                                   >> 1624   {
                                                   >> 1625     // Inside the plane
                                                   >> 1626     if (Comp>0)
                                                   >> 1627     {
                                                   >> 1628       // Will leave extent
                                                   >> 1629       vdist=-pdist/Comp;
                                                   >> 1630       if (vdist<snxt)
                                                   >> 1631       {
                                                   >> 1632         snxt=vdist;
                                                   >> 1633         side=ks3;
                                                   >> 1634       }
                                                   >> 1635     }
                                                   >> 1636   }
                                                   >> 1637   else
                                                   >> 1638   {
                                                   >> 1639     // On surface
                                                   >> 1640     if (Comp>0)
                                                   >> 1641     {
                                                   >> 1642       if (calcNorm)
                                                   >> 1643       {
                                                   >> 1644         *validNorm=true;
                                                   >> 1645         *n=G4ThreeVector(fPlanes[3].a,fPlanes[3].b,fPlanes[3].c);
                                                   >> 1646       }
                                                   >> 1647       return snxt=0;
                                                   >> 1648     }
                                                   >> 1649   }
                                                   >> 1650 
                                                   >> 1651   // set normal
1018   if (calcNorm)                                  1652   if (calcNorm)
1019   {                                              1653   {
1020     *validNorm = true;                        << 1654     *validNorm=true;
1021     if (iside < 0)                            << 1655     switch(side)
1022       n->set(0, 0, iside + 3); // (-4+3)=-1,  << 1656     {
1023     else                                      << 1657       case ks0:
1024       n->set(fPlanes[iside].a, fPlanes[iside] << 1658         *n=G4ThreeVector(fPlanes[0].a,fPlanes[0].b,fPlanes[0].c);
                                                   >> 1659         break;
                                                   >> 1660       case ks1:
                                                   >> 1661         *n=G4ThreeVector(fPlanes[1].a,fPlanes[1].b,fPlanes[1].c);
                                                   >> 1662         break;
                                                   >> 1663       case ks2:
                                                   >> 1664         *n=G4ThreeVector(fPlanes[2].a,fPlanes[2].b,fPlanes[2].c);
                                                   >> 1665         break;
                                                   >> 1666       case ks3:
                                                   >> 1667         *n=G4ThreeVector(fPlanes[3].a,fPlanes[3].b,fPlanes[3].c);
                                                   >> 1668         break;
                                                   >> 1669       case kMZ:
                                                   >> 1670         *n=G4ThreeVector(0,0,-1);
                                                   >> 1671         break;
                                                   >> 1672       case kPZ:
                                                   >> 1673         *n=G4ThreeVector(0,0,1);
                                                   >> 1674         break;
                                                   >> 1675       default:
                                                   >> 1676         G4cout << G4endl;
                                                   >> 1677         DumpInfo();
                                                   >> 1678         std::ostringstream message;
                                                   >> 1679         G4int oldprc = message.precision(16);
                                                   >> 1680         message << "Undefined side for valid surface normal to solid."
                                                   >> 1681                 << G4endl
                                                   >> 1682                 << "Position:"  << G4endl << G4endl
                                                   >> 1683                 << "p.x() = "   << p.x()/mm << " mm" << G4endl
                                                   >> 1684                 << "p.y() = "   << p.y()/mm << " mm" << G4endl
                                                   >> 1685                 << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl
                                                   >> 1686                 << "Direction:" << G4endl << G4endl
                                                   >> 1687                 << "v.x() = "   << v.x() << G4endl
                                                   >> 1688                 << "v.y() = "   << v.y() << G4endl
                                                   >> 1689                 << "v.z() = "   << v.z() << G4endl << G4endl
                                                   >> 1690                 << "Proposed distance :" << G4endl << G4endl
                                                   >> 1691                 << "snxt = "    << snxt/mm << " mm" << G4endl;
                                                   >> 1692         message.precision(oldprc);
                                                   >> 1693         G4Exception("G4Trap::DistanceToOut(p,v,..)","GeomSolids1002",
                                                   >> 1694                     JustWarning, message);
                                                   >> 1695         break;
                                                   >> 1696     }
1025   }                                              1697   }
1026   return tmax;                                << 1698   return snxt;
1027 }                                                1699 }
1028                                                  1700 
1029 ///////////////////////////////////////////// << 1701 //////////////////////////////////////////////////////////////////////////////
1030 //                                               1702 //
1031 // Calculate exact shortest distance to any b    1703 // Calculate exact shortest distance to any boundary from inside
1032 // - Returns 0 is ThreeVector outside            1704 // - Returns 0 is ThreeVector outside
1033                                                  1705 
1034 G4double G4Trap::DistanceToOut( const G4Three    1706 G4double G4Trap::DistanceToOut( const G4ThreeVector& p ) const
1035 {                                                1707 {
                                                   >> 1708   G4double safe=0.0,Dist;
                                                   >> 1709   G4int i;
                                                   >> 1710 
1036 #ifdef G4CSGDEBUG                                1711 #ifdef G4CSGDEBUG
1037   if( Inside(p) == kOutside )                    1712   if( Inside(p) == kOutside )
1038   {                                              1713   {
1039     std::ostringstream message;               << 1714      G4int oldprc = G4cout.precision(16) ;
1040     G4long oldprc = message.precision(16);    << 1715      G4cout << G4endl ;
1041     message << "Point p is outside (!?) of so << 1716      DumpInfo();
1042     message << "Position:\n";                 << 1717      G4cout << "Position:"  << G4endl << G4endl ;
1043     message << "   p.x() = " << p.x()/mm << " << 1718      G4cout << "p.x() = "   << p.x()/mm << " mm" << G4endl ;
1044     message << "   p.y() = " << p.y()/mm << " << 1719      G4cout << "p.y() = "   << p.y()/mm << " mm" << G4endl ;
1045     message << "   p.z() = " << p.z()/mm << " << 1720      G4cout << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl ;
1046     G4cout.precision(oldprc);                 << 1721      G4cout.precision(oldprc) ;
1047     G4Exception("G4Trap::DistanceToOut(p)", " << 1722      G4Exception("G4Trap::DistanceToOut(p)",
1048                 JustWarning, message );       << 1723                  "GeomSolids1002", JustWarning, "Point p is outside !?" );
1049     DumpInfo();                               << 
1050   }                                              1724   }
1051 #endif                                           1725 #endif
1052   switch (fTrapType)                          << 1726 
                                                   >> 1727   safe=fDz-std::fabs(p.z());
                                                   >> 1728   if (safe<0) safe=0;
                                                   >> 1729   else
1053   {                                              1730   {
1054     case 0: // General case                   << 1731     for (i=0;i<4;i++)
1055     {                                            1732     {
1056       G4double dz = std::abs(p.z())-fDz;      << 1733       Dist=-(fPlanes[i].a*p.x()+fPlanes[i].b*p.y()
1057       G4double dy1 = fPlanes[0].b*p.y()+fPlan << 1734             +fPlanes[i].c*p.z()+fPlanes[i].d);
1058       G4double dy2 = fPlanes[1].b*p.y()+fPlan << 1735       if (Dist<safe) safe=Dist;
1059       G4double dy = std::max(dz,std::max(dy1, << 1736     }
1060                                               << 1737     if (safe<0) safe=0;
1061       G4double dx1 = fPlanes[2].a*p.x()+fPlan << 1738   }
1062                    + fPlanes[2].c*p.z()+fPlan << 1739   return safe;  
1063       G4double dx2 = fPlanes[3].a*p.x()+fPlan << 
1064                    + fPlanes[3].c*p.z()+fPlan << 
1065       G4double dist = std::max(dy,std::max(dx << 
1066       return (dist < 0) ? -dist : 0.;         << 
1067     }                                         << 
1068     case 1: // YZ section is a rectangle      << 
1069     {                                         << 
1070       G4double dz = std::abs(p.z())-fDz;      << 
1071       G4double dy = std::max(dz,std::abs(p.y( << 
1072       G4double dx1 = fPlanes[2].a*p.x()+fPlan << 
1073                    + fPlanes[2].c*p.z()+fPlan << 
1074       G4double dx2 = fPlanes[3].a*p.x()+fPlan << 
1075                    + fPlanes[3].c*p.z()+fPlan << 
1076       G4double dist = std::max(dy,std::max(dx << 
1077       return (dist < 0) ? -dist : 0.;         << 
1078     }                                         << 
1079     case 2: // YZ section is a rectangle and  << 
1080     {       // XZ section is an isosceles tra << 
1081       G4double dz = std::abs(p.z())-fDz;      << 
1082       G4double dy = std::max(dz,std::abs(p.y( << 
1083       G4double dx = fPlanes[3].a*std::abs(p.x << 
1084                   + fPlanes[3].c*p.z()+fPlane << 
1085       G4double dist = std::max(dy,dx);        << 
1086       return (dist < 0) ? -dist : 0.;         << 
1087     }                                         << 
1088     case 3: // YZ section is a rectangle and  << 
1089     {       // XY section is an isosceles tra << 
1090       G4double dz = std::abs(p.z())-fDz;      << 
1091       G4double dy = std::max(dz,std::abs(p.y( << 
1092       G4double dx = fPlanes[3].a*std::abs(p.x << 
1093                   + fPlanes[3].b*p.y()+fPlane << 
1094       G4double dist = std::max(dy,dx);        << 
1095       return (dist < 0) ? -dist : 0.;         << 
1096     }                                         << 
1097   }                                           << 
1098   return 0.;                                  << 
1099 }                                                1740 }
1100                                                  1741 
1101 /////////////////////////////////////////////    1742 //////////////////////////////////////////////////////////////////////////
1102 //                                               1743 //
1103 // GetEntityType                              << 1744 // Create a List containing the transformed vertices
1104                                               << 1745 // Ordering [0-3] -fDz cross section
1105 G4GeometryType G4Trap::GetEntityType() const  << 1746 //          [4-7] +fDz cross section such that [0] is below [4],
1106 {                                             << 1747 //                                             [1] below [5] etc.
1107   return {"G4Trap"};                          << 1748 // Note:
                                                   >> 1749 //  Caller has deletion resposibility
                                                   >> 1750 
                                                   >> 1751 G4ThreeVectorList*
                                                   >> 1752 G4Trap::CreateRotatedVertices( const G4AffineTransform& pTransform ) const
                                                   >> 1753 {
                                                   >> 1754   G4ThreeVectorList *vertices;
                                                   >> 1755   vertices=new G4ThreeVectorList();
                                                   >> 1756   if (vertices)
                                                   >> 1757   {
                                                   >> 1758     vertices->reserve(8);
                                                   >> 1759     G4ThreeVector vertex0(-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
                                                   >> 1760                           -fDz*fTthetaSphi-fDy1,-fDz);
                                                   >> 1761     G4ThreeVector vertex1(-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
                                                   >> 1762                           -fDz*fTthetaSphi-fDy1,-fDz);
                                                   >> 1763     G4ThreeVector vertex2(-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
                                                   >> 1764                           -fDz*fTthetaSphi+fDy1,-fDz);
                                                   >> 1765     G4ThreeVector vertex3(-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
                                                   >> 1766                           -fDz*fTthetaSphi+fDy1,-fDz);
                                                   >> 1767     G4ThreeVector vertex4(+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
                                                   >> 1768                           +fDz*fTthetaSphi-fDy2,+fDz);
                                                   >> 1769     G4ThreeVector vertex5(+fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
                                                   >> 1770                           +fDz*fTthetaSphi-fDy2,+fDz);
                                                   >> 1771     G4ThreeVector vertex6(+fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
                                                   >> 1772                           +fDz*fTthetaSphi+fDy2,+fDz);
                                                   >> 1773     G4ThreeVector vertex7(+fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
                                                   >> 1774                           +fDz*fTthetaSphi+fDy2,+fDz);
                                                   >> 1775 
                                                   >> 1776     vertices->push_back(pTransform.TransformPoint(vertex0));
                                                   >> 1777     vertices->push_back(pTransform.TransformPoint(vertex1));
                                                   >> 1778     vertices->push_back(pTransform.TransformPoint(vertex2));
                                                   >> 1779     vertices->push_back(pTransform.TransformPoint(vertex3));
                                                   >> 1780     vertices->push_back(pTransform.TransformPoint(vertex4));
                                                   >> 1781     vertices->push_back(pTransform.TransformPoint(vertex5));
                                                   >> 1782     vertices->push_back(pTransform.TransformPoint(vertex6));
                                                   >> 1783     vertices->push_back(pTransform.TransformPoint(vertex7));
                                                   >> 1784   }
                                                   >> 1785   else
                                                   >> 1786   {
                                                   >> 1787     DumpInfo();
                                                   >> 1788     G4Exception("G4Trap::CreateRotatedVertices()",
                                                   >> 1789                 "GeomSolids0003", FatalException,
                                                   >> 1790                 "Error in allocation of vertices. Out of memory !");
                                                   >> 1791   }
                                                   >> 1792   return vertices;
1108 }                                                1793 }
1109                                                  1794 
1110 /////////////////////////////////////////////    1795 //////////////////////////////////////////////////////////////////////////
1111 //                                               1796 //
1112 // IsFaceted                                  << 1797 // GetEntityType
1113                                                  1798 
1114 G4bool G4Trap::IsFaceted() const              << 1799 G4GeometryType G4Trap::GetEntityType() const
1115 {                                                1800 {
1116   return true;                                << 1801   return G4String("G4Trap");
1117 }                                                1802 }
1118                                                  1803 
1119 /////////////////////////////////////////////    1804 //////////////////////////////////////////////////////////////////////////
1120 //                                               1805 //
1121 // Make a clone of the object                    1806 // Make a clone of the object
1122 //                                               1807 //
1123 G4VSolid* G4Trap::Clone() const                  1808 G4VSolid* G4Trap::Clone() const
1124 {                                                1809 {
1125   return new G4Trap(*this);                      1810   return new G4Trap(*this);
1126 }                                                1811 }
1127                                                  1812 
1128 /////////////////////////////////////////////    1813 //////////////////////////////////////////////////////////////////////////
1129 //                                               1814 //
1130 // Stream object contents to an output stream    1815 // Stream object contents to an output stream
1131                                                  1816 
1132 std::ostream& G4Trap::StreamInfo( std::ostrea    1817 std::ostream& G4Trap::StreamInfo( std::ostream& os ) const
1133 {                                                1818 {
1134   G4double phi    = GetPhi();                 << 1819   G4int oldprc = os.precision(16);
1135   G4double theta  = GetTheta();               << 
1136   G4double alpha1 = GetAlpha1();              << 
1137   G4double alpha2 = GetAlpha2();              << 
1138                                               << 
1139   G4long oldprc = os.precision(16);           << 
1140   os << "------------------------------------    1820   os << "-----------------------------------------------------------\n"
1141      << "    *** Dump for solid: " << GetName << 1821      << "    *** Dump for solid - " << GetName() << " ***\n"
1142      << "    ================================    1822      << "    ===================================================\n"
1143      << " Solid type: G4Trap\n"                  1823      << " Solid type: G4Trap\n"
1144      << " Parameters:\n"                      << 1824      << " Parameters: \n"
1145      << "    half length Z: " << fDz/mm << "  << 1825      << "    half length Z: " << fDz/mm << " mm \n"
1146      << "    half length Y, face -Dz: " << fD << 1826      << "    half length Y of face -fDz: " << fDy1/mm << " mm \n"
1147      << "    half length X, face -Dz, side -D << 1827      << "    half length X of side -fDy1, face -fDz: " << fDx1/mm << " mm \n"
1148      << "    half length X, face -Dz, side +D << 1828      << "    half length X of side +fDy1, face -fDz: " << fDx2/mm << " mm \n"
1149      << "    half length Y, face +Dz: " << fD << 1829      << "    half length Y of face +fDz: " << fDy2/mm << " mm \n"
1150      << "    half length X, face +Dz, side -D << 1830      << "    half length X of side -fDy2, face +fDz: " << fDx3/mm << " mm \n"
1151      << "    half length X, face +Dz, side +D << 1831      << "    half length X of side +fDy2, face +fDz: " << fDx4/mm << " mm \n"
1152      << "    theta: " << theta/degree << " de << 1832      << "    std::tan(theta)*std::cos(phi): " << fTthetaCphi/degree << " degrees \n"
1153      << "    phi:   " << phi/degree << " degr << 1833      << "    std::tan(theta)*std::sin(phi): " << fTthetaSphi/degree << " degrees \n"
1154      << "    alpha, face -Dz: " << alpha1/deg << 1834      << "    std::tan(alpha), -fDz: " << fTalpha1/degree << " degrees \n"
1155      << "    alpha, face +Dz: " << alpha2/deg << 1835      << "    std::tan(alpha), +fDz: " << fTalpha2/degree << " degrees \n"
                                                   >> 1836      << "    trap side plane equations:\n"
                                                   >> 1837      << "        " << fPlanes[0].a << " X + " << fPlanes[0].b << " Y + "
                                                   >> 1838                    << fPlanes[0].c << " Z + " << fPlanes[0].d << " = 0\n"
                                                   >> 1839      << "        " << fPlanes[1].a << " X + " << fPlanes[1].b << " Y + "
                                                   >> 1840                    << fPlanes[1].c << " Z + " << fPlanes[1].d << " = 0\n"
                                                   >> 1841      << "        " << fPlanes[2].a << " X + " << fPlanes[2].b << " Y + "
                                                   >> 1842                    << fPlanes[2].c << " Z + " << fPlanes[2].d << " = 0\n"
                                                   >> 1843      << "        " << fPlanes[3].a << " X + " << fPlanes[3].b << " Y + "
                                                   >> 1844                    << fPlanes[3].c << " Z + " << fPlanes[3].d << " = 0\n"
1156      << "------------------------------------    1845      << "-----------------------------------------------------------\n";
1157   os.precision(oldprc);                          1846   os.precision(oldprc);
1158                                                  1847 
1159   return os;                                     1848   return os;
1160 }                                                1849 }
1161                                                  1850 
1162 ///////////////////////////////////////////// << 1851 /////////////////////////////////////////////////////////////////////////
                                                   >> 1852 //
                                                   >> 1853 // GetPointOnPlane
1163 //                                               1854 //
1164 // Compute vertices from planes               << 1855 // Auxiliary method for Get Point on Surface
1165                                                  1856 
1166 void G4Trap::GetVertices(G4ThreeVector pt[8]) << 1857 G4ThreeVector G4Trap::GetPointOnPlane(G4ThreeVector p0, G4ThreeVector p1, 
1167 {                                             << 1858                                       G4ThreeVector p2, G4ThreeVector p3,
1168   for (G4int i=0; i<8; ++i)                   << 1859                                       G4double& area) const
1169   {                                           << 1860 {
1170     G4int iy = (i==0 || i==1 || i==4 || i==5) << 1861   G4double lambda1, lambda2, chose, aOne, aTwo;
1171     G4int ix = (i==0 || i==2 || i==4 || i==6) << 1862   G4ThreeVector t, u, v, w, Area, normal;
1172     G4double z = (i < 4) ? -fDz : fDz;        << 1863   
1173     G4double y = -(fPlanes[iy].c*z + fPlanes[ << 1864   t = p1 - p0;
1174     G4double x = -(fPlanes[ix].b*y + fPlanes[ << 1865   u = p2 - p1;
1175                    + fPlanes[ix].d)/fPlanes[i << 1866   v = p3 - p2;
1176     pt[i].set(x,y,z);                         << 1867   w = p0 - p3;
1177   }                                           << 1868 
                                                   >> 1869   Area = G4ThreeVector(w.y()*v.z() - w.z()*v.y(),
                                                   >> 1870                        w.z()*v.x() - w.x()*v.z(),
                                                   >> 1871                        w.x()*v.y() - w.y()*v.x());
                                                   >> 1872   
                                                   >> 1873   aOne = 0.5*Area.mag();
                                                   >> 1874   
                                                   >> 1875   Area = G4ThreeVector(t.y()*u.z() - t.z()*u.y(),
                                                   >> 1876                        t.z()*u.x() - t.x()*u.z(),
                                                   >> 1877                        t.x()*u.y() - t.y()*u.x());
                                                   >> 1878   
                                                   >> 1879   aTwo = 0.5*Area.mag();
                                                   >> 1880   
                                                   >> 1881   area = aOne + aTwo;
                                                   >> 1882   
                                                   >> 1883   chose = RandFlat::shoot(0.,aOne+aTwo);
                                                   >> 1884 
                                                   >> 1885   if( (chose>=0.) && (chose < aOne) )
                                                   >> 1886   {
                                                   >> 1887     lambda1 = RandFlat::shoot(0.,1.);
                                                   >> 1888     lambda2 = RandFlat::shoot(0.,lambda1);
                                                   >> 1889     return (p2+lambda1*v+lambda2*w);    
                                                   >> 1890   }
                                                   >> 1891   
                                                   >> 1892   // else
                                                   >> 1893 
                                                   >> 1894   lambda1 = RandFlat::shoot(0.,1.);
                                                   >> 1895   lambda2 = RandFlat::shoot(0.,lambda1);
                                                   >> 1896 
                                                   >> 1897   return (p0+lambda1*t+lambda2*u);    
1178 }                                                1898 }
1179                                                  1899 
1180 ///////////////////////////////////////////// << 1900 ///////////////////////////////////////////////////////////////
1181 //                                               1901 //
1182 // Generate random point on the surface       << 1902 // GetPointOnSurface
1183                                                  1903 
1184 G4ThreeVector G4Trap::GetPointOnSurface() con    1904 G4ThreeVector G4Trap::GetPointOnSurface() const
1185 {                                                1905 {
1186   // Set indeces                              << 1906   G4double aOne, aTwo, aThree, aFour, aFive, aSix, chose;
1187   constexpr G4int iface [6][4] =              << 1907   G4ThreeVector One, Two, Three, Four, Five, Six, test;
1188     { {0,1,3,2}, {0,4,5,1}, {2,3,7,6}, {0,2,6 << 
1189                                               << 
1190   // Set vertices                             << 
1191   G4ThreeVector pt[8];                           1908   G4ThreeVector pt[8];
1192   GetVertices(pt);                            << 1909      
1193                                               << 1910   pt[0] = G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
1194   // Select face                              << 1911                         -fDz*fTthetaSphi-fDy1,-fDz);
1195   //                                          << 1912   pt[1] = G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
1196   G4double select = fAreas[5]*G4QuickRand();  << 1913                         -fDz*fTthetaSphi-fDy1,-fDz);
1197   G4int k = 5;                                << 1914   pt[2] = G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
1198   k -= (G4int)(select <= fAreas[4]);          << 1915                         -fDz*fTthetaSphi+fDy1,-fDz);
1199   k -= (G4int)(select <= fAreas[3]);          << 1916   pt[3] = G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
1200   k -= (G4int)(select <= fAreas[2]);          << 1917                         -fDz*fTthetaSphi+fDy1,-fDz);
1201   k -= (G4int)(select <= fAreas[1]);          << 1918   pt[4] = G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
1202   k -= (G4int)(select <= fAreas[0]);          << 1919                         +fDz*fTthetaSphi-fDy2,+fDz);
1203                                               << 1920   pt[5] = G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
1204   // Select sub-triangle                      << 1921                         +fDz*fTthetaSphi-fDy2,+fDz);
1205   //                                          << 1922   pt[6] = G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
1206   G4int i0 = iface[k][0];                     << 1923                         +fDz*fTthetaSphi+fDy2,+fDz);
1207   G4int i1 = iface[k][1];                     << 1924   pt[7] = G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
1208   G4int i2 = iface[k][2];                     << 1925                         +fDz*fTthetaSphi+fDy2,+fDz);
1209   G4int i3 = iface[k][3];                     << 1926   
1210   G4double s2 = G4GeomTools::TriangleAreaNorm << 1927   // make sure we provide the points in a clockwise fashion
1211   if (select > fAreas[k] - s2) i0 = i2;       << 1928 
1212                                               << 1929   One   = GetPointOnPlane(pt[0],pt[1],pt[3],pt[2], aOne);
1213   // Generate point                           << 1930   Two   = GetPointOnPlane(pt[4],pt[5],pt[7],pt[6], aTwo);
1214   //                                          << 1931   Three = GetPointOnPlane(pt[6],pt[7],pt[3],pt[2], aThree);
1215   G4double u = G4QuickRand();                 << 1932   Four  = GetPointOnPlane(pt[4],pt[5],pt[1],pt[0], aFour); 
1216   G4double v = G4QuickRand();                 << 1933   Five  = GetPointOnPlane(pt[0],pt[2],pt[6],pt[4], aFive);
1217   if (u + v > 1.) { u = 1. - u; v = 1. - v; } << 1934   Six   = GetPointOnPlane(pt[1],pt[3],pt[7],pt[5], aSix);
1218   return (1.-u-v)*pt[i0] + u*pt[i1] + v*pt[i3 << 1935  
                                                   >> 1936   chose = RandFlat::shoot(0.,aOne+aTwo+aThree+aFour+aFive+aSix);
                                                   >> 1937   if( (chose>=0.) && (chose<aOne) )                    
                                                   >> 1938     { return One; }
                                                   >> 1939   else if( (chose>=aOne) && (chose<aOne+aTwo) )  
                                                   >> 1940     { return Two; }
                                                   >> 1941   else if( (chose>=aOne+aTwo) && (chose<aOne+aTwo+aThree) )
                                                   >> 1942     { return Three; }
                                                   >> 1943   else if( (chose>=aOne+aTwo+aThree) && (chose<aOne+aTwo+aThree+aFour) )
                                                   >> 1944     { return Four; }
                                                   >> 1945   else if( (chose>=aOne+aTwo+aThree+aFour)
                                                   >> 1946         && (chose<aOne+aTwo+aThree+aFour+aFive) )
                                                   >> 1947     { return Five; }
                                                   >> 1948   return Six;
1219 }                                                1949 }
1220                                                  1950 
1221 /////////////////////////////////////////////    1951 //////////////////////////////////////////////////////////////////////////
1222 //                                               1952 //
1223 // Methods for visualisation                     1953 // Methods for visualisation
1224                                                  1954 
1225 void G4Trap::DescribeYourselfTo ( G4VGraphics    1955 void G4Trap::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
1226 {                                                1956 {
1227   scene.AddSolid (*this);                        1957   scene.AddSolid (*this);
1228 }                                                1958 }
1229                                                  1959 
1230 G4Polyhedron* G4Trap::CreatePolyhedron () con    1960 G4Polyhedron* G4Trap::CreatePolyhedron () const
1231 {                                                1961 {
1232   G4double phi = std::atan2(fTthetaSphi, fTth    1962   G4double phi = std::atan2(fTthetaSphi, fTthetaCphi);
1233   G4double alpha1 = std::atan(fTalpha1);         1963   G4double alpha1 = std::atan(fTalpha1);
1234   G4double alpha2 = std::atan(fTalpha2);         1964   G4double alpha2 = std::atan(fTalpha2);
1235   G4double theta = std::atan(std::sqrt(fTthet << 1965   G4double theta = std::atan(std::sqrt(fTthetaCphi*fTthetaCphi+fTthetaSphi*fTthetaSphi));
1236                                       +fTthet << 
1237                                                  1966 
1238   return new G4PolyhedronTrap(fDz, theta, phi    1967   return new G4PolyhedronTrap(fDz, theta, phi,
1239                               fDy1, fDx1, fDx    1968                               fDy1, fDx1, fDx2, alpha1,
1240                               fDy2, fDx3, fDx    1969                               fDy2, fDx3, fDx4, alpha2);
1241 }                                                1970 }
1242                                                  1971 
1243 #endif                                        << 1972 G4NURBS* G4Trap::CreateNURBS () const
                                                   >> 1973 {
                                                   >> 1974    // return new G4NURBSbox (fDx, fDy, fDz);
                                                   >> 1975    return 0 ;
                                                   >> 1976 }
1244                                                  1977