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 8.3)


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