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


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