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


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