Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/solids/specific/src/G4Tet.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/specific/src/G4Tet.cc (Version 11.3.0) and /geometry/solids/specific/src/G4Tet.cc (Version 9.3.p2)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the  intell     18 // * This  code  implementation is the  intellectual property  of the *
 19 // * Vanderbilt University Free Electron Laser     19 // * Vanderbilt University Free Electron Laser Center                 *
 20 // * Vanderbilt University, Nashville, TN, USA     20 // * Vanderbilt University, Nashville, TN, USA                        *
 21 // * Development supported by:                     21 // * Development supported by:                                        *
 22 // * United States MFEL program  under grant F     22 // * United States MFEL program  under grant FA9550-04-1-0045         *
 23 // * and NASA under contract number NNG04CT05P     23 // * and NASA under contract number NNG04CT05P                        *
 24 // * Written by Marcus H. Mendenhall and Rober     24 // * Written by Marcus H. Mendenhall and Robert A. Weller.            *
 25 // *                                               25 // *                                                                  *
 26 // * Contributed to the Geant4 Core, January,      26 // * Contributed to the Geant4 Core, January, 2005.                   *
 27 // *                                               27 // *                                                                  *
 28 // *******************************************     28 // ********************************************************************
 29 //                                                 29 //
                                                   >>  30 // $Id: G4Tet.cc,v 1.11.12.1 2010/09/08 15:54:59 gcosmo Exp $
                                                   >>  31 // GEANT4 tag $Name: geant4-09-03-patch-02 $
                                                   >>  32 //
                                                   >>  33 // class G4Tet
                                                   >>  34 //
 30 // Implementation for G4Tet class                  35 // Implementation for G4Tet class
 31 //                                                 36 //
 32 // 03.09.2004 - Marcus Mendenhall, created     <<  37 // History:
 33 // 08.01.2020 - Evgueni Tcherniaev, complete r <<  38 //
                                                   >>  39 //  20040903 - Marcus Mendenhall, created G4Tet
                                                   >>  40 //  20041101 - Marcus Mendenhall, optimized constant dot products with
                                                   >>  41 //             fCdotNijk values
                                                   >>  42 //  20041101 - MHM removed tracking error by clipping DistanceToOut to 0
                                                   >>  43 //             for surface cases
                                                   >>  44 //  20041101 - MHM many speed optimizations in if statements
                                                   >>  45 //  20041101 - MHM changed vdotn comparisons to 1e-12 instead of 0.0 to
                                                   >>  46 //             avoid nearly-parallel problems
                                                   >>  47 //  20041102 - MHM Added extra distance into solid to DistanceToIn(p,v)
                                                   >>  48 //             hit testing
                                                   >>  49 //  20041102 - MHM added ability to check for degeneracy without throwing
                                                   >>  50 //             G4Exception
                                                   >>  51 //  20041103 - MHM removed many unused variables from class
                                                   >>  52 //  20040803 - Dionysios Anninos, added GetPointOnSurface() method
                                                   >>  53 //  20061112 - MHM added code for G4VSolid GetSurfaceArea()
                                                   >>  54 //
 34 // -------------------------------------------     55 // --------------------------------------------------------------------
 35                                                    56 
 36 #include "G4Tet.hh"                                57 #include "G4Tet.hh"
 37                                                    58 
 38 #if !defined(G4GEOM_USE_UTET)                  <<  59 const char G4Tet::CVSVers[]="$Id: G4Tet.cc,v 1.11.12.1 2010/09/08 15:54:59 gcosmo Exp $";
 39                                                    60 
 40 #include "G4VoxelLimits.hh"                        61 #include "G4VoxelLimits.hh"
 41 #include "G4AffineTransform.hh"                    62 #include "G4AffineTransform.hh"
 42 #include "G4BoundingEnvelope.hh"               << 
 43                                                    63 
 44 #include "G4VPVParameterisation.hh"                64 #include "G4VPVParameterisation.hh"
 45                                                    65 
 46 #include "G4QuickRand.hh"                      <<  66 #include "Randomize.hh"
 47                                                    67 
 48 #include "G4VGraphicsScene.hh"                     68 #include "G4VGraphicsScene.hh"
 49 #include "G4Polyhedron.hh"                         69 #include "G4Polyhedron.hh"
                                                   >>  70 #include "G4NURBS.hh"
                                                   >>  71 #include "G4NURBSbox.hh"
 50 #include "G4VisExtent.hh"                          72 #include "G4VisExtent.hh"
 51                                                    73 
 52 #include "G4AutoLock.hh"                       <<  74 #include "G4ThreeVector.hh"
 53                                                    75 
 54 namespace                                      <<  76 #include <cmath>
 55 {                                              << 
 56   G4Mutex polyhedronMutex = G4MUTEX_INITIALIZE << 
 57 }                                              << 
 58                                                    77 
 59 using namespace CLHEP;                             78 using namespace CLHEP;
 60                                                    79 
 61 //////////////////////////////////////////////     80 ////////////////////////////////////////////////////////////////////////
 62 //                                                 81 //
 63 // Constructor - create a tetrahedron              82 // Constructor - create a tetrahedron
                                                   >>  83 // This class is implemented separately from general polyhedra,
                                                   >>  84 // because the simplex geometry can be computed very quickly,
                                                   >>  85 // which may become important in situations imported from mesh generators,
                                                   >>  86 // in which a very large number of G4Tets are created.
 64 // A Tet has all of its geometrical informatio     87 // A Tet has all of its geometrical information precomputed
 65 //                                             <<  88 
 66 G4Tet::G4Tet(const G4String& pName,                89 G4Tet::G4Tet(const G4String& pName,
 67              const G4ThreeVector& p0,          <<  90                    G4ThreeVector anchor,
 68              const G4ThreeVector& p1,          <<  91                    G4ThreeVector p2,
 69              const G4ThreeVector& p2,          <<  92                    G4ThreeVector p3,
 70              const G4ThreeVector& p3, G4bool*  <<  93                    G4ThreeVector p4, G4bool *degeneracyFlag)
 71   : G4VSolid(pName)                            <<  94   : G4VSolid(pName), fpPolyhedron(0), warningFlag(0)
 72 {                                              <<  95 {
 73   // Check for degeneracy                      <<  96   // fV<x><y> is vector from vertex <y> to vertex <x>
 74   G4bool degenerate = CheckDegeneracy(p0, p1,  <<  97   //
 75   if (degeneracyFlag != nullptr)               <<  98   G4ThreeVector fV21=p2-anchor;
                                                   >>  99   G4ThreeVector fV31=p3-anchor;
                                                   >> 100   G4ThreeVector fV41=p4-anchor;
                                                   >> 101 
                                                   >> 102   // make sure this is a correctly oriented set of points for the tetrahedron
                                                   >> 103   //
                                                   >> 104   G4double signed_vol=fV21.cross(fV31).dot(fV41);
                                                   >> 105   if(signed_vol<0.0)
 76   {                                               106   {
 77     *degeneracyFlag = degenerate;              << 107     G4ThreeVector temp(p4);
 78   }                                            << 108     p4=p3;
                                                   >> 109     p3=temp;
                                                   >> 110     temp=fV41;
                                                   >> 111     fV41=fV31;
                                                   >> 112     fV31=temp; 
                                                   >> 113   }
                                                   >> 114   fCubicVolume = std::fabs(signed_vol) / 6.;
                                                   >> 115 
                                                   >> 116   G4ThreeVector fV24=p2-p4;
                                                   >> 117   G4ThreeVector fV43=p4-p3;
                                                   >> 118   G4ThreeVector fV32=p3-p2;
                                                   >> 119 
                                                   >> 120   fXMin=std::min(std::min(std::min(anchor.x(), p2.x()),p3.x()),p4.x());
                                                   >> 121   fXMax=std::max(std::max(std::max(anchor.x(), p2.x()),p3.x()),p4.x());
                                                   >> 122   fYMin=std::min(std::min(std::min(anchor.y(), p2.y()),p3.y()),p4.y());
                                                   >> 123   fYMax=std::max(std::max(std::max(anchor.y(), p2.y()),p3.y()),p4.y());
                                                   >> 124   fZMin=std::min(std::min(std::min(anchor.z(), p2.z()),p3.z()),p4.z());
                                                   >> 125   fZMax=std::max(std::max(std::max(anchor.z(), p2.z()),p3.z()),p4.z());
                                                   >> 126 
                                                   >> 127   fDx=(fXMax-fXMin)*0.5; fDy=(fYMax-fYMin)*0.5; fDz=(fZMax-fZMin)*0.5;
                                                   >> 128 
                                                   >> 129   fMiddle=G4ThreeVector(fXMax+fXMin, fYMax+fYMin, fZMax+fZMin)*0.5;
                                                   >> 130   fMaxSize=std::max(std::max(std::max((anchor-fMiddle).mag(),
                                                   >> 131                                       (p2-fMiddle).mag()),
                                                   >> 132                              (p3-fMiddle).mag()),
                                                   >> 133                     (p4-fMiddle).mag());
                                                   >> 134 
                                                   >> 135   G4bool degenerate=std::fabs(signed_vol) < 1e-9*fMaxSize*fMaxSize*fMaxSize;
                                                   >> 136 
                                                   >> 137   if(degeneracyFlag) *degeneracyFlag=degenerate;
 79   else if (degenerate)                            138   else if (degenerate)
 80   {                                               139   {
 81     std::ostringstream message;                << 140     G4Exception("G4Tet::G4Tet()", "InvalidSetup", FatalException,
 82     message << "Degenerate tetrahedron: " << G << 141                 "Degenerate tetrahedron not allowed.");
 83             << "  anchor: " << p0 << "\n"      << 
 84             << "  p1    : " << p1 << "\n"      << 
 85             << "  p2    : " << p2 << "\n"      << 
 86             << "  p3    : " << p3 << "\n"      << 
 87             << "  volume: "                    << 
 88             << std::abs((p1 - p0).cross(p2 - p << 
 89     G4Exception("G4Tet::G4Tet()", "GeomSolids0 << 
 90   }                                               142   }
 91                                                   143 
 92   // Define surface thickness                  << 144   fTol=1e-9*(std::fabs(fXMin)+std::fabs(fXMax)+std::fabs(fYMin)
 93   halfTolerance = 0.5 * kCarTolerance;         << 145             +std::fabs(fYMax)+std::fabs(fZMin)+std::fabs(fZMax));
 94                                                << 146   //fTol=kCarTolerance;
 95   // Set data members                          << 147 
 96   Initialize(p0, p1, p2, p3);                  << 148   fAnchor=anchor;
                                                   >> 149   fP2=p2;
                                                   >> 150   fP3=p3;
                                                   >> 151   fP4=p4;
                                                   >> 152 
                                                   >> 153   G4ThreeVector fCenter123=(anchor+p2+p3)*(1.0/3.0); // face center
                                                   >> 154   G4ThreeVector fCenter134=(anchor+p4+p3)*(1.0/3.0);
                                                   >> 155   G4ThreeVector fCenter142=(anchor+p4+p2)*(1.0/3.0);
                                                   >> 156   G4ThreeVector fCenter234=(p2+p3+p4)*(1.0/3.0);
                                                   >> 157 
                                                   >> 158   // compute area of each triangular face by cross product
                                                   >> 159   // and sum for total surface area
                                                   >> 160 
                                                   >> 161   G4ThreeVector normal123=fV31.cross(fV21);
                                                   >> 162   G4ThreeVector normal134=fV41.cross(fV31);
                                                   >> 163   G4ThreeVector normal142=fV21.cross(fV41);
                                                   >> 164   G4ThreeVector normal234=fV32.cross(fV43);
                                                   >> 165 
                                                   >> 166   fSurfaceArea=(
                                                   >> 167       normal123.mag()+
                                                   >> 168       normal134.mag()+
                                                   >> 169       normal142.mag()+
                                                   >> 170       normal234.mag()
                                                   >> 171   )/2.0;
                                                   >> 172 
                                                   >> 173   fNormal123=normal123.unit();
                                                   >> 174   fNormal134=normal134.unit();
                                                   >> 175   fNormal142=normal142.unit();
                                                   >> 176   fNormal234=normal234.unit();
                                                   >> 177 
                                                   >> 178   fCdotN123=fCenter123.dot(fNormal123);
                                                   >> 179   fCdotN134=fCenter134.dot(fNormal134);
                                                   >> 180   fCdotN142=fCenter142.dot(fNormal142);
                                                   >> 181   fCdotN234=fCenter234.dot(fNormal234);
 97 }                                                 182 }
 98                                                   183 
 99 ////////////////////////////////////////////// << 184 //////////////////////////////////////////////////////////////////////////
100 //                                                185 //
101 // Fake default constructor - sets only member    186 // Fake default constructor - sets only member data and allocates memory
102 //                            for usage restri    187 //                            for usage restricted to object persistency.
103 //                                                188 //
104 G4Tet::G4Tet( __void__& a )                       189 G4Tet::G4Tet( __void__& a )
105   : G4VSolid(a)                                << 190   : G4VSolid(a), fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0),
                                                   >> 191     fAnchor(0,0,0), fP2(0,0,0), fP3(0,0,0), fP4(0,0,0), fMiddle(0,0,0),
                                                   >> 192     fNormal123(0,0,0), fNormal142(0,0,0), fNormal134(0,0,0),
                                                   >> 193     fNormal234(0,0,0), warningFlag(0),
                                                   >> 194     fCdotN123(0.), fCdotN142(0.), fCdotN134(0.), fCdotN234(0.),
                                                   >> 195     fXMin(0.), fXMax(0.), fYMin(0.), fYMax(0.), fZMin(0.), fZMax(0.),
                                                   >> 196     fDx(0.), fDy(0.), fDz(0.), fTol(0.), fMaxSize(0.)
106 {                                                 197 {
107 }                                                 198 }
108                                                   199 
109 ////////////////////////////////////////////// << 200 //////////////////////////////////////////////////////////////////////////
110 //                                                201 //
111 // Destructor                                     202 // Destructor
112 //                                             << 203 
113 G4Tet::~G4Tet()                                   204 G4Tet::~G4Tet()
114 {                                                 205 {
115   delete fpPolyhedron; fpPolyhedron = nullptr; << 206   delete fpPolyhedron;
116 }                                                 207 }
117                                                   208 
118 ////////////////////////////////////////////// << 209 //////////////////////////////////////////////////////////////////////////
119 //                                             << 
120 // Copy constructor                            << 
121 //                                                210 //
122 G4Tet::G4Tet(const G4Tet& rhs)                 << 211 // CheckDegeneracy
123   : G4VSolid(rhs)                              << 
124 {                                              << 
125    halfTolerance = rhs.halfTolerance;          << 
126    fCubicVolume = rhs.fCubicVolume;            << 
127    fSurfaceArea = rhs.fSurfaceArea;            << 
128    for (G4int i = 0; i < 4; ++i) { fVertex[i]  << 
129    for (G4int i = 0; i < 4; ++i) { fNormal[i]  << 
130    for (G4int i = 0; i < 4; ++i) { fDist[i] =  << 
131    for (G4int i = 0; i < 4; ++i) { fArea[i] =  << 
132    fBmin = rhs.fBmin;                          << 
133    fBmax = rhs.fBmax;                          << 
134 }                                              << 
135                                                   212 
136 ////////////////////////////////////////////// << 213 G4bool G4Tet::CheckDegeneracy( G4ThreeVector anchor,
137 //                                             << 214                                G4ThreeVector p2,
138 // Assignment operator                         << 215                                G4ThreeVector p3,
139 //                                             << 216                                G4ThreeVector p4 )
140 G4Tet& G4Tet::operator = (const G4Tet& rhs)    << 
141 {                                                 217 {
142    // Check assignment to self                 << 218   G4bool result;
143    //                                          << 219   G4Tet *object=new G4Tet("temp",anchor,p2,p3,p4,&result);
144    if (this == &rhs)  { return *this; }        << 220   delete object;
145                                                << 221   return result;
146    // Copy base class data                     << 
147    //                                          << 
148    G4VSolid::operator=(rhs);                   << 
149                                                << 
150    // Copy data                                << 
151    //                                          << 
152    halfTolerance = rhs.halfTolerance;          << 
153    fCubicVolume = rhs.fCubicVolume;            << 
154    fSurfaceArea = rhs.fSurfaceArea;            << 
155    for (G4int i = 0; i < 4; ++i) { fVertex[i]  << 
156    for (G4int i = 0; i < 4; ++i) { fNormal[i]  << 
157    for (G4int i = 0; i < 4; ++i) { fDist[i] =  << 
158    for (G4int i = 0; i < 4; ++i) { fArea[i] =  << 
159    fBmin = rhs.fBmin;                          << 
160    fBmax = rhs.fBmax;                          << 
161    fRebuildPolyhedron = false;                 << 
162    delete fpPolyhedron; fpPolyhedron = nullptr << 
163                                                << 
164    return *this;                               << 
165 }                                                 222 }
166                                                   223 
167 ////////////////////////////////////////////// << 224 //////////////////////////////////////////////////////////////////////////
168 //                                             << 
169 // Return true if tetrahedron is degenerate    << 
170 // Tetrahedron is concidered as degenerate in  << 
171 // height is less than degeneracy tolerance    << 
172 //                                                225 //
173 G4bool G4Tet::CheckDegeneracy(const G4ThreeVec << 226 // Dispatch to parameterisation for replication mechanism dimension
174                               const G4ThreeVec << 227 // computation & modification.
175                               const G4ThreeVec << 
176                               const G4ThreeVec << 
177 {                                              << 
178   G4double hmin = 4. * kCarTolerance; // degen << 
179                                                << 
180   // Calculate volume                          << 
181   G4double vol = std::abs((p1 - p0).cross(p2 - << 
182                                                << 
183   // Calculate face areas squared              << 
184   G4double ss[4];                              << 
185   ss[0] = ((p1 - p0).cross(p2 - p0)).mag2();   << 
186   ss[1] = ((p2 - p0).cross(p3 - p0)).mag2();   << 
187   ss[2] = ((p3 - p0).cross(p1 - p0)).mag2();   << 
188   ss[3] = ((p2 - p1).cross(p3 - p1)).mag2();   << 
189                                                << 
190   // Find face with max area                   << 
191   G4int k = 0;                                 << 
192   for (G4int i = 1; i < 4; ++i) { if (ss[i] >  << 
193                                                   228 
194   // Check: vol^2 / s^2 <= hmin^2              << 229 void G4Tet::ComputeDimensions(G4VPVParameterisation* ,
195   return (vol*vol <= ss[k]*hmin*hmin);         << 230                               const G4int ,
                                                   >> 231                               const G4VPhysicalVolume* )
                                                   >> 232 {
196 }                                                 233 }
197                                                   234 
198 ////////////////////////////////////////////// << 235 //////////////////////////////////////////////////////////////////////////
199 //                                             << 
200 // Set data members                            << 
201 //                                                236 //
202 void G4Tet::Initialize(const G4ThreeVector& p0 << 237 // Calculate extent under transform and specified limit
203                        const G4ThreeVector& p1 << 
204                        const G4ThreeVector& p2 << 
205                        const G4ThreeVector& p3 << 
206 {                                              << 
207   // Set vertices                              << 
208   fVertex[0] = p0;                             << 
209   fVertex[1] = p1;                             << 
210   fVertex[2] = p2;                             << 
211   fVertex[3] = p3;                             << 
212                                                   238 
213   G4ThreeVector norm[4];                       << 239 G4bool G4Tet::CalculateExtent(const EAxis pAxis,
214   norm[0] = (p2 - p0).cross(p1 - p0);          << 240                               const G4VoxelLimits& pVoxelLimit,
215   norm[1] = (p3 - p0).cross(p2 - p0);          << 241                               const G4AffineTransform& pTransform,
216   norm[2] = (p1 - p0).cross(p3 - p0);          << 242                                     G4double& pMin, G4double& pMax) const
217   norm[3] = (p2 - p1).cross(p3 - p1);          << 243 {
218   G4double volume = norm[0].dot(p3 - p0);      << 244   G4double xMin,xMax;
219   if (volume > 0.)                             << 245   G4double yMin,yMax;
220   {                                            << 246   G4double zMin,zMax;
221     for (auto & i : norm) { i = -i; }          << 247 
                                                   >> 248   if (pTransform.IsRotated())
                                                   >> 249   {
                                                   >> 250     G4ThreeVector pp0=pTransform.TransformPoint(fAnchor);
                                                   >> 251     G4ThreeVector pp1=pTransform.TransformPoint(fP2);
                                                   >> 252     G4ThreeVector pp2=pTransform.TransformPoint(fP3);
                                                   >> 253     G4ThreeVector pp3=pTransform.TransformPoint(fP4);
                                                   >> 254 
                                                   >> 255     xMin    = std::min(std::min(std::min(pp0.x(), pp1.x()),pp2.x()),pp3.x());
                                                   >> 256     xMax    = std::max(std::max(std::max(pp0.x(), pp1.x()),pp2.x()),pp3.x());
                                                   >> 257     yMin    = std::min(std::min(std::min(pp0.y(), pp1.y()),pp2.y()),pp3.y());
                                                   >> 258     yMax    = std::max(std::max(std::max(pp0.y(), pp1.y()),pp2.y()),pp3.y());
                                                   >> 259     zMin    = std::min(std::min(std::min(pp0.z(), pp1.z()),pp2.z()),pp3.z());
                                                   >> 260     zMax    = std::max(std::max(std::max(pp0.z(), pp1.z()),pp2.z()),pp3.z());
                                                   >> 261 
                                                   >> 262   }
                                                   >> 263   else
                                                   >> 264   {
                                                   >> 265     G4double xoffset = pTransform.NetTranslation().x() ;
                                                   >> 266     xMin    = xoffset + fXMin;
                                                   >> 267     xMax    = xoffset + fXMax;
                                                   >> 268     G4double yoffset = pTransform.NetTranslation().y() ;
                                                   >> 269     yMin    = yoffset + fYMin;
                                                   >> 270     yMax    = yoffset + fYMax;
                                                   >> 271     G4double zoffset = pTransform.NetTranslation().z() ;
                                                   >> 272     zMin    = zoffset + fZMin;
                                                   >> 273     zMax    = zoffset + fZMax;
                                                   >> 274   }
                                                   >> 275 
                                                   >> 276   if (pVoxelLimit.IsXLimited())
                                                   >> 277   {
                                                   >> 278     if ( (xMin > pVoxelLimit.GetMaxXExtent()+fTol) || 
                                                   >> 279          (xMax < pVoxelLimit.GetMinXExtent()-fTol)  )  { return false; }
                                                   >> 280     else
                                                   >> 281     {
                                                   >> 282       xMin = std::max(xMin, pVoxelLimit.GetMinXExtent());
                                                   >> 283       xMax = std::min(xMax, pVoxelLimit.GetMaxXExtent());
                                                   >> 284     }
222   }                                               285   }
223                                                   286 
224   // Set normals to face planes                << 287   if (pVoxelLimit.IsYLimited())
225   for (G4int i = 0; i < 4; ++i) { fNormal[i] = << 
226                                                << 
227   // Set distances to planes                   << 
228   for (G4int i = 0; i < 3; ++i) { fDist[i] = f << 
229   fDist[3] = fNormal[3].dot(p1);               << 
230                                                << 
231   // Set face areas                            << 
232   for (G4int i = 0; i < 4; ++i) { fArea[i] = 0 << 
233                                                << 
234   // Set bounding box                          << 
235   for (G4int i = 0; i < 3; ++i)                << 
236   {                                               288   {
237     fBmin[i] = std::min(std::min(std::min(p0[i << 289     if ( (yMin > pVoxelLimit.GetMaxYExtent()+fTol) ||
238     fBmax[i] = std::max(std::max(std::max(p0[i << 290          (yMax < pVoxelLimit.GetMinYExtent()-fTol)  )  { return false; }
239   }                                            << 291     else
240                                                << 292     {
241   // Set volume and surface area               << 293       yMin = std::max(yMin, pVoxelLimit.GetMinYExtent());
242   fCubicVolume = std::abs(volume)/6.;          << 294       yMax = std::min(yMax, pVoxelLimit.GetMaxYExtent());
243   fSurfaceArea = fArea[0] + fArea[1] + fArea[2 << 295     }
244 }                                              << 296     }
245                                                   297 
246 ////////////////////////////////////////////// << 298     if (pVoxelLimit.IsZLimited())
247 //                                             << 299     {
248 // Set vertices                                << 300       if ( (zMin > pVoxelLimit.GetMaxZExtent()+fTol) ||
249 //                                             << 301            (zMax < pVoxelLimit.GetMinZExtent()-fTol)  )  { return false; }
250 void G4Tet::SetVertices(const G4ThreeVector& p << 302     else
251                         const G4ThreeVector& p << 303     {
252                         const G4ThreeVector& p << 304       zMin = std::max(zMin, pVoxelLimit.GetMinZExtent());
253                         const G4ThreeVector& p << 305       zMax = std::min(zMax, pVoxelLimit.GetMaxZExtent());
254 {                                              << 306     }
255   // Check for degeneracy                      << 
256   G4bool degenerate = CheckDegeneracy(p0, p1,  << 
257   if (degeneracyFlag != nullptr)               << 
258   {                                            << 
259     *degeneracyFlag = degenerate;              << 
260   }                                               307   }
261   else if (degenerate)                         << 308 
                                                   >> 309   switch (pAxis)
262   {                                               310   {
263     std::ostringstream message;                << 311     case kXAxis:
264     message << "Degenerate tetrahedron is not  << 312       pMin=xMin;
265             << "  anchor: " << p0 << "\n"      << 313       pMax=xMax;
266             << "  p1    : " << p1 << "\n"      << 314       break;
267             << "  p2    : " << p2 << "\n"      << 315     case kYAxis:
268             << "  p3    : " << p3 << "\n"      << 316       pMin=yMin;
269             << "  volume: "                    << 317       pMax=yMax;
270             << std::abs((p1 - p0).cross(p2 - p << 318       break;
271     G4Exception("G4Tet::SetVertices()", "GeomS << 319     case kZAxis:
272                 FatalException, message);      << 320       pMin=zMin;
                                                   >> 321       pMax=zMax;
                                                   >> 322       break;
                                                   >> 323     default:
                                                   >> 324       break;
273   }                                               325   }
274                                                   326 
275   // Set data members                          << 327   return true;
276   Initialize(p0, p1, p2, p3);                  << 328 } 
277                                                << 
278   // Set flag to rebuild polyhedron            << 
279   fRebuildPolyhedron = true;                   << 
280 }                                              << 
281                                                   329 
282 ////////////////////////////////////////////// << 330 /////////////////////////////////////////////////////////////////////////
283 //                                             << 
284 // Return four vertices                        << 
285 //                                                331 //
286 void G4Tet::GetVertices(G4ThreeVector& p0,     << 332 // Return whether point inside/outside/on surface, using tolerance
287                         G4ThreeVector& p1,     << 
288                         G4ThreeVector& p2,     << 
289                         G4ThreeVector& p3) con << 
290 {                                              << 
291   p0 = fVertex[0];                             << 
292   p1 = fVertex[1];                             << 
293   p2 = fVertex[2];                             << 
294   p3 = fVertex[3];                             << 
295 }                                              << 
296                                                   333 
297 ////////////////////////////////////////////// << 334 EInside G4Tet::Inside(const G4ThreeVector& p) const
298 //                                             << 
299 // Return std::vector of vertices              << 
300 //                                             << 
301 std::vector<G4ThreeVector> G4Tet::GetVertices( << 
302 {                                                 335 {
303   std::vector<G4ThreeVector> vertices(4);      << 336   G4double r123, r134, r142, r234;
304   for (G4int i = 0; i < 4; ++i) { vertices[i]  << 
305   return vertices;                             << 
306 }                                              << 
307                                                   337 
308 ////////////////////////////////////////////// << 338   // this is written to allow if-statement truncation so the outside test
309 //                                             << 339   // (where most of the world is) can fail very quickly and efficiently
310 // Dispatch to parameterisation for replicatio << 
311 // computation & modification.                 << 
312 //                                             << 
313 void G4Tet::ComputeDimensions(G4VPVParameteris << 
314                               const G4int ,    << 
315                               const G4VPhysica << 
316 {                                              << 
317 }                                              << 
318                                                   340 
319 ////////////////////////////////////////////// << 341   if ( (r123=p.dot(fNormal123)-fCdotN123) > fTol ||
320 //                                             << 342        (r134=p.dot(fNormal134)-fCdotN134) > fTol ||
321 // Set bounding box                            << 343        (r142=p.dot(fNormal142)-fCdotN142) > fTol ||
322 //                                             << 344        (r234=p.dot(fNormal234)-fCdotN234) > fTol )
323 void G4Tet::SetBoundingLimits(const G4ThreeVec << 345   {
324                               const G4ThreeVec << 346     return kOutside; // at least one is out!
325 {                                              << 347   }
326   G4int iout[4] = { 0, 0, 0, 0 };              << 348   else if( (r123 < -fTol)&&(r134 < -fTol)&&(r142 < -fTol)&&(r234 < -fTol) )
327   for (G4int i = 0; i < 4; ++i)                << 
328   {                                               349   {
329     iout[i] = (G4int)(fVertex[i].x() < pMin.x( << 350     return kInside; // all are definitively inside
330                       fVertex[i].y() < pMin.y( << 
331                       fVertex[i].z() < pMin.z( << 
332                       fVertex[i].x() > pMax.x( << 
333                       fVertex[i].y() > pMax.y( << 
334                       fVertex[i].z() > pMax.z( << 
335   }                                               351   }
336   if (iout[0] + iout[1] + iout[2] + iout[3] != << 352   else
337   {                                               353   {
338     std::ostringstream message;                << 354     return kSurface; // too close to tell
339     message << "Attempt to set bounding box th << 
340             << GetName() << " !\n"             << 
341             << "  Specified bounding box limit << 
342             << "    pmin: " << pMin << "\n"    << 
343             << "    pmax: " << pMax << "\n"    << 
344             << "  Tetrahedron vertices:\n"     << 
345             << "    anchor " << fVertex[0] <<  << 
346             << "    p1 "     << fVertex[1] <<  << 
347             << "    p2 "     << fVertex[2] <<  << 
348             << "    p3 "     << fVertex[3] <<  << 
349     G4Exception("G4Tet::SetBoundingLimits()",  << 
350                 FatalException, message);      << 
351   }                                               355   }
352   fBmin = pMin;                                << 
353   fBmax = pMax;                                << 
354 }                                                 356 }
355                                                   357 
356 ////////////////////////////////////////////// << 358 ///////////////////////////////////////////////////////////////////////
357 //                                                359 //
358 // Return bounding box                         << 360 // Calculate side nearest to p, and return normal
359 //                                             << 361 // If two sides are equidistant, normal of first side (x/y/z) 
360 void G4Tet::BoundingLimits(G4ThreeVector& pMin << 362 // encountered returned.
                                                   >> 363 // This assumes that we are looking from the inside!
                                                   >> 364 
                                                   >> 365 G4ThreeVector G4Tet::SurfaceNormal( const G4ThreeVector& p) const
361 {                                                 366 {
362   pMin = fBmin;                                << 367   G4double r123=std::fabs(p.dot(fNormal123)-fCdotN123);
363   pMax = fBmax;                                << 368   G4double r134=std::fabs(p.dot(fNormal134)-fCdotN134);
364 }                                              << 369   G4double r142=std::fabs(p.dot(fNormal142)-fCdotN142);
                                                   >> 370   G4double r234=std::fabs(p.dot(fNormal234)-fCdotN234);
                                                   >> 371 
                                                   >> 372   if( (r123<=r134) && (r123<=r142) && (r123<=r234) )  { return fNormal123; }
                                                   >> 373   else if ( (r134<=r142) && (r134<=r234) )  { return fNormal134; }
                                                   >> 374   else if (r142 <= r234)  { return fNormal142; }
                                                   >> 375   return fNormal234;
                                                   >> 376 }
                                                   >> 377 
                                                   >> 378 ///////////////////////////////////////////////////////////////////////////
                                                   >> 379 //
                                                   >> 380 // Calculate distance to box from an outside point
                                                   >> 381 // - return kInfinity if no intersection.
                                                   >> 382 // All this is very unrolled, for speed.
365                                                   383 
366 ////////////////////////////////////////////// << 384 G4double G4Tet::DistanceToIn(const G4ThreeVector& p,
367 //                                             << 385                              const G4ThreeVector& v) const
368 // Calculate extent under transform and specif << 
369 //                                             << 
370 G4bool G4Tet::CalculateExtent(const EAxis pAxi << 
371                               const G4VoxelLim << 
372                               const G4AffineTr << 
373                                     G4double&  << 
374 {                                                 386 {
375   G4ThreeVector bmin, bmax;                    << 387     G4ThreeVector vu(v.unit()), hp;
                                                   >> 388     G4double vdotn, t, tmin=kInfinity;
376                                                   389 
377   // Check bounding box (bbox)                 << 390     G4double extraDistance=10.0*fTol; // a little ways into the solid
378   //                                           << 
379   BoundingLimits(bmin,bmax);                   << 
380   G4BoundingEnvelope bbox(bmin,bmax);          << 
381                                                   391 
382   // Use simple bounding-box to help in the ca << 392     vdotn=-vu.dot(fNormal123);
383   //                                           << 393     if(vdotn > 1e-12)
384   return bbox.CalculateExtent(pAxis,pVoxelLimi << 394     { // this is a candidate face, since it is pointing at us
                                                   >> 395       t=(p.dot(fNormal123)-fCdotN123)/vdotn; // #  distance to intersection
                                                   >> 396       if( (t>=-fTol) && (t<tmin) )
                                                   >> 397       { // if not true, we're going away from this face or it's not close
                                                   >> 398         hp=p+vu*(t+extraDistance); // a little beyond point of intersection
                                                   >> 399         if ( ( hp.dot(fNormal134)-fCdotN134 < 0.0 ) &&
                                                   >> 400              ( hp.dot(fNormal142)-fCdotN142 < 0.0 ) &&
                                                   >> 401              ( hp.dot(fNormal234)-fCdotN234 < 0.0 ) )
                                                   >> 402         {
                                                   >> 403           tmin=t;
                                                   >> 404         }
                                                   >> 405       }
                                                   >> 406     }
385                                                   407 
386 #if 0                                          << 408     vdotn=-vu.dot(fNormal134);
387   // Precise extent computation (disabled by d << 409     if(vdotn > 1e-12)
388   //                                           << 410     { // # this is a candidate face, since it is pointing at us
389   G4bool exist;                                << 411       t=(p.dot(fNormal134)-fCdotN134)/vdotn; // #  distance to intersection
390   if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox << 412       if( (t>=-fTol) && (t<tmin) )
391   {                                            << 413       { // if not true, we're going away from this face
392     return exist = (pMin < pMax) ? true : fals << 414         hp=p+vu*(t+extraDistance); // a little beyond point of intersection
393   }                                            << 415         if ( ( hp.dot(fNormal123)-fCdotN123 < 0.0 ) && 
                                                   >> 416              ( hp.dot(fNormal142)-fCdotN142 < 0.0 ) &&
                                                   >> 417              ( hp.dot(fNormal234)-fCdotN234 < 0.0 ) )
                                                   >> 418         {
                                                   >> 419           tmin=t;
                                                   >> 420         }
                                                   >> 421       }
                                                   >> 422     }
394                                                   423 
395   // Set bounding envelope (benv) and calculat << 424     vdotn=-vu.dot(fNormal142);
396   //                                           << 425     if(vdotn > 1e-12)
397   std::vector<G4ThreeVector> vec = GetVertices << 426     { // # this is a candidate face, since it is pointing at us
                                                   >> 427       t=(p.dot(fNormal142)-fCdotN142)/vdotn; // #  distance to intersection
                                                   >> 428       if( (t>=-fTol) && (t<tmin) )
                                                   >> 429       { // if not true, we're going away from this face
                                                   >> 430         hp=p+vu*(t+extraDistance); // a little beyond point of intersection
                                                   >> 431         if ( ( hp.dot(fNormal123)-fCdotN123 < 0.0 ) &&
                                                   >> 432              ( hp.dot(fNormal134)-fCdotN134 < 0.0 ) &&
                                                   >> 433              ( hp.dot(fNormal234)-fCdotN234 < 0.0 ) )
                                                   >> 434         {
                                                   >> 435           tmin=t;
                                                   >> 436         }
                                                   >> 437       }
                                                   >> 438     }
398                                                   439 
399   G4ThreeVectorList anchor(1);                 << 440     vdotn=-vu.dot(fNormal234);
400   anchor[0].set(vec[0].x(),vec[0].y(),vec[0].z << 441     if(vdotn > 1e-12)
                                                   >> 442     { // # this is a candidate face, since it is pointing at us
                                                   >> 443       t=(p.dot(fNormal234)-fCdotN234)/vdotn; // #  distance to intersection
                                                   >> 444       if( (t>=-fTol) && (t<tmin) )
                                                   >> 445       { // if not true, we're going away from this face
                                                   >> 446         hp=p+vu*(t+extraDistance); // a little beyond point of intersection
                                                   >> 447         if ( ( hp.dot(fNormal123)-fCdotN123 < 0.0 ) &&
                                                   >> 448              ( hp.dot(fNormal134)-fCdotN134 < 0.0 ) &&
                                                   >> 449              ( hp.dot(fNormal142)-fCdotN142 < 0.0 ) )
                                                   >> 450         {
                                                   >> 451           tmin=t;
                                                   >> 452         }
                                                   >> 453       }
                                                   >> 454     }
401                                                   455 
402   G4ThreeVectorList base(3);                   << 456   return std::max(0.0,tmin);
403   base[0].set(vec[1].x(),vec[1].y(),vec[1].z() << 457 }
404   base[1].set(vec[2].x(),vec[2].y(),vec[2].z() << 
405   base[2].set(vec[3].x(),vec[3].y(),vec[3].z() << 
406                                                   458 
407   std::vector<const G4ThreeVectorList *> polyg << 459 //////////////////////////////////////////////////////////////////////////
408   polygons[0] = &anchor;                       << 460 // 
409   polygons[1] = &base;                         << 461 // Approximate distance to tet.
                                                   >> 462 // returns distance to sphere centered on bounding box
                                                   >> 463 // - If inside return 0
410                                                   464 
411   G4BoundingEnvelope benv(bmin,bmax,polygons); << 465 G4double G4Tet::DistanceToIn(const G4ThreeVector& p) const
412   return exists = benv.CalculateExtent(pAxis,p << 466 {
413 #endif                                         << 467   G4double dd=(p-fMiddle).mag() - fMaxSize - fTol;
                                                   >> 468   return std::max(0.0, dd);
414 }                                                 469 }
415                                                   470 
416 ////////////////////////////////////////////// << 471 /////////////////////////////////////////////////////////////////////////
417 //                                                472 //
418 // Return whether point inside/outside/on surf << 473 // Calcluate distance to surface of box from inside
419 //                                             << 474 // by calculating distances to box's x/y/z planes.
420 EInside G4Tet::Inside(const G4ThreeVector& p)  << 475 // Smallest distance is exact distance to exiting.
                                                   >> 476 
                                                   >> 477 G4double G4Tet::DistanceToOut( const G4ThreeVector& p,const G4ThreeVector& v,
                                                   >> 478                                const G4bool calcNorm,
                                                   >> 479                                      G4bool *validNorm, G4ThreeVector *n) const
421 {                                                 480 {
422   G4double dd[4];                              << 481     G4ThreeVector vu(v.unit());
423   for (G4int i = 0; i < 4; ++i) { dd[i] = fNor << 482     G4double t1=kInfinity,t2=kInfinity,t3=kInfinity,t4=kInfinity, vdotn, tt;
424                                                   483 
425   G4double dist = std::max(std::max(std::max(d << 484     vdotn=vu.dot(fNormal123);
426   return (dist <= -halfTolerance) ?            << 485     if(vdotn > 1e-12)  // #we're heading towards this face, so it is a candidate
427     kInside : ((dist <= halfTolerance) ? kSurf << 486     {
428 }                                              << 487       t1=(fCdotN123-p.dot(fNormal123))/vdotn; // #  distance to intersection
                                                   >> 488     }
429                                                   489 
430 ////////////////////////////////////////////// << 490     vdotn=vu.dot(fNormal134);
431 //                                             << 491     if(vdotn > 1e-12) // #we're heading towards this face, so it is a candidate
432 // Return unit normal to surface at p          << 492     {
433 //                                             << 493       t2=(fCdotN134-p.dot(fNormal134))/vdotn; // #  distance to intersection
434 G4ThreeVector G4Tet::SurfaceNormal( const G4Th << 494     }
435 {                                              << 
436   G4double k[4];                               << 
437   for (G4int i = 0; i < 4; ++i)                << 
438   {                                            << 
439     k[i] = (G4double)(std::abs(fNormal[i].dot( << 
440   }                                            << 
441   G4double nsurf = k[0] + k[1] + k[2] + k[3];  << 
442   G4ThreeVector norm =                         << 
443     k[0]*fNormal[0] + k[1]*fNormal[1] + k[2]*f << 
444                                                << 
445   if (nsurf == 1.) return norm;                << 
446   else if (nsurf > 1.) return norm.unit(); //  << 
447   {                                            << 
448 #ifdef G4SPECSDEBUG                            << 
449     std::ostringstream message;                << 
450     G4long oldprc = message.precision(16);     << 
451     message << "Point p is not on surface (!?) << 
452             << GetName() << "\n";              << 
453     message << "Position:\n";                  << 
454     message << "   p.x() = " << p.x()/mm << "  << 
455     message << "   p.y() = " << p.y()/mm << "  << 
456     message << "   p.z() = " << p.z()/mm << "  << 
457     G4cout.precision(oldprc);                  << 
458     G4Exception("G4Tet::SurfaceNormal(p)", "Ge << 
459                 JustWarning, message );        << 
460     DumpInfo();                                << 
461 #endif                                         << 
462     return ApproxSurfaceNormal(p);             << 
463   }                                            << 
464 }                                              << 
465                                                   495 
466 ////////////////////////////////////////////// << 496     vdotn=vu.dot(fNormal142);
467 //                                             << 497     if(vdotn > 1e-12) // #we're heading towards this face, so it is a candidate
468 // Find surface nearest to point and return co << 498     {
469 // This method normally should not be called   << 499       t3=(fCdotN142-p.dot(fNormal142))/vdotn; // #  distance to intersection
470 //                                             << 500     }
471 G4ThreeVector G4Tet::ApproxSurfaceNormal(const << 
472 {                                              << 
473   G4double dist = -DBL_MAX;                    << 
474   G4int iside = 0;                             << 
475   for (G4int i = 0; i < 4; ++i)                << 
476   {                                            << 
477     G4double d = fNormal[i].dot(p) - fDist[i]; << 
478     if (d > dist) { dist = d; iside = i; }     << 
479   }                                            << 
480   return fNormal[iside];                       << 
481 }                                              << 
482                                                   501 
483 ////////////////////////////////////////////// << 502     vdotn=vu.dot(fNormal234);
484 //                                             << 503     if(vdotn > 1e-12) // #we're heading towards this face, so it is a candidate
485 // Calculate distance to surface from outside, << 504     {
486 // return kInfinity if no intersection         << 505       t4=(fCdotN234-p.dot(fNormal234))/vdotn; // #  distance to intersection
487 //                                             << 506     }
488 G4double G4Tet::DistanceToIn(const G4ThreeVect << 507 
489                              const G4ThreeVect << 508     tt=std::min(std::min(std::min(t1,t2),t3),t4);
490 {                                              << 509 
491   G4double tin = -DBL_MAX, tout = DBL_MAX;     << 510     if (warningFlag && (tt == kInfinity || tt < -fTol))
492   for (G4int i = 0; i < 4; ++i)                << 
493   {                                            << 
494     G4double cosa = fNormal[i].dot(v);         << 
495     G4double dist = fNormal[i].dot(p) - fDist[ << 
496     if (dist >= -halfTolerance)                << 
497     {                                             511     {
498       if (cosa >= 0.) { return kInfinity; }    << 512       DumpInfo();
499       tin = std::max(tin, -dist/cosa);         << 513       G4cout << "p = " << p / mm << "mm" << G4endl;
                                                   >> 514       G4cout << "v = " << v  << G4endl;
                                                   >> 515       G4cout << "t1, t2, t3, t4 (mm) "
                                                   >> 516              << t1/mm << ", " << t2/mm << ", " << t3/mm << ", " << t4/mm
                                                   >> 517              << G4endl << G4endl;
                                                   >> 518       G4Exception("G4Tet::DistanceToOut(p,v,...)", "Notification", JustWarning,
                                                   >> 519                   "No good intersection found or already outside!?" );
                                                   >> 520       if(validNorm)
                                                   >> 521       {
                                                   >> 522         *validNorm=false; // flag normal as meaningless
                                                   >> 523       }
500     }                                             524     }
501     else if (cosa > 0.)                        << 525     else if(calcNorm && n)
502     {                                             526     {
503       tout = std::min(tout, -dist/cosa);       << 527       static G4ThreeVector normal;
                                                   >> 528       if(tt==t1)        { normal=fNormal123; }
                                                   >> 529       else if (tt==t2)  { normal=fNormal134; }
                                                   >> 530       else if (tt==t3)  { normal=fNormal142; }
                                                   >> 531       else if (tt==t4)  { normal=fNormal234; }
                                                   >> 532       n=&normal;
                                                   >> 533       if(validNorm) { *validNorm=true; }
504     }                                             534     }
505   }                                            << 
506                                                   535 
507   return (tout - tin <= halfTolerance) ?       << 536     return std::max(tt,0.0); // avoid tt<0.0 by a tiny bit
508     kInfinity : ((tin < halfTolerance) ? 0. :  << 537                              // if we are right on a face
509 }                                                 538 }
510                                                   539 
511 ////////////////////////////////////////////// << 540 ////////////////////////////////////////////////////////////////////////////
512 //                                                541 //
513 // Estimate safety distance to surface from ou << 542 // Calculate exact shortest distance to any boundary from inside
514 //                                             << 543 // - If outside return 0
515 G4double G4Tet::DistanceToIn(const G4ThreeVect << 544 
                                                   >> 545 G4double G4Tet::DistanceToOut(const G4ThreeVector& p) const
516 {                                                 546 {
517   G4double dd[4];                              << 547   G4double t1,t2,t3,t4;
518   for (G4int i = 0; i < 4; ++i) { dd[i] = fNor << 548   t1=fCdotN123-p.dot(fNormal123); //  distance to plane, positive if inside
                                                   >> 549   t2=fCdotN134-p.dot(fNormal134); //  distance to plane
                                                   >> 550   t3=fCdotN142-p.dot(fNormal142); //  distance to plane
                                                   >> 551   t4=fCdotN234-p.dot(fNormal234); //  distance to plane
519                                                   552 
520   G4double dist = std::max(std::max(std::max(d << 553   // if any one of these is negative, we are outside,
521   return (dist > 0.) ? dist : 0.;              << 554   // so return zero in that case
                                                   >> 555 
                                                   >> 556   G4double tmin=std::min(std::min(std::min(t1,t2),t3),t4);
                                                   >> 557   return (tmin < fTol)? 0:tmin;
522 }                                                 558 }
523                                                   559 
524 //////////////////////////////////////////////    560 ////////////////////////////////////////////////////////////////////////
525 //                                                561 //
526 // Calcluate distance to surface from inside   << 562 // Create a List containing the transformed vertices
527 //                                             << 563 // Note: Caller has deletion responsibility
528 G4double G4Tet::DistanceToOut(const G4ThreeVec << 564 
529                               const G4ThreeVec << 565 G4ThreeVectorList*
530                               const G4bool cal << 566 G4Tet::CreateRotatedVertices(const G4AffineTransform& pTransform) const
531                                     G4bool* va << 
532                                     G4ThreeVec << 
533 {                                                 567 {
534   // Calculate distances and cosines           << 568   G4ThreeVectorList* vertices = new G4ThreeVectorList();
535   G4double cosa[4], dist[4];                   << 
536   G4int ind[4] = {0}, nside = 0;               << 
537   for (G4int i = 0; i < 4; ++i)                << 
538   {                                            << 
539     G4double tmp = fNormal[i].dot(v);          << 
540     cosa[i] = tmp;                             << 
541     ind[nside] = (G4int)(tmp > 0) * i;         << 
542     nside += (G4int)(tmp > 0);                 << 
543     dist[i] = fNormal[i].dot(p) - fDist[i];    << 
544   }                                            << 
545                                                   569 
546   // Find intersection (in most of cases nside << 570   if (vertices)
547   G4double tout = DBL_MAX;                     << 
548   G4int iside = 0;                             << 
549   for (G4int i = 0; i < nside; ++i)            << 
550   {                                               571   {
551     G4int k = ind[i];                          << 572     vertices->reserve(4);
552     // Check: leaving the surface              << 573     G4ThreeVector vertex0(fAnchor);
553     if (dist[k] >= -halfTolerance) { tout = 0. << 574     G4ThreeVector vertex1(fP2);
554     // Compute distance to intersection        << 575     G4ThreeVector vertex2(fP3);
555     G4double tmp = -dist[k]/cosa[k];           << 576     G4ThreeVector vertex3(fP4);
556     if (tmp < tout) { tout = tmp; iside = k; } << 
557   }                                            << 
558                                                   577 
559   // Set normal, if required, and return dista << 578     vertices->push_back(pTransform.TransformPoint(vertex0));
560   if (calcNorm)                                << 579     vertices->push_back(pTransform.TransformPoint(vertex1));
                                                   >> 580     vertices->push_back(pTransform.TransformPoint(vertex2));
                                                   >> 581     vertices->push_back(pTransform.TransformPoint(vertex3));
                                                   >> 582   }
                                                   >> 583   else
561   {                                               584   {
562     *validNorm = true;                         << 585     DumpInfo();
563     *n = fNormal[iside];                       << 586     G4Exception("G4Tet::CreateRotatedVertices()",
                                                   >> 587                 "FatalError", FatalException,
                                                   >> 588                 "Error in allocation of vertices. Out of memory !");
564   }                                               589   }
565   return tout;                                 << 590   return vertices;
566 }                                              << 
567                                                << 
568 ////////////////////////////////////////////// << 
569 //                                             << 
570 // Calculate safety distance to surface from i << 
571 //                                             << 
572 G4double G4Tet::DistanceToOut(const G4ThreeVec << 
573 {                                              << 
574   G4double dd[4];                              << 
575   for (G4int i = 0; i < 4; ++i) { dd[i] = fDis << 
576                                                << 
577   G4double dist = std::min(std::min(std::min(d << 
578   return (dist > 0.) ? dist : 0.;              << 
579 }                                                 591 }
580                                                   592 
581 ////////////////////////////////////////////// << 593 //////////////////////////////////////////////////////////////////////////
582 //                                                594 //
583 // GetEntityType                                  595 // GetEntityType
584 //                                             << 596 
585 G4GeometryType G4Tet::GetEntityType() const       597 G4GeometryType G4Tet::GetEntityType() const
586 {                                                 598 {
587   return {"G4Tet"};                            << 599   return G4String("G4Tet");
588 }                                                 600 }
589                                                   601 
590 ////////////////////////////////////////////// << 602 //////////////////////////////////////////////////////////////////////////
591 //                                             << 
592 // IsFaceted                                   << 
593 //                                                603 //
594 G4bool G4Tet::IsFaceted() const                << 604 // Stream object contents to an output stream
                                                   >> 605 
                                                   >> 606 std::ostream& G4Tet::StreamInfo(std::ostream& os) const
595 {                                                 607 {
596   return true;                                 << 608   os << "-----------------------------------------------------------\n"
                                                   >> 609   << "    *** Dump for solid - " << GetName() << " ***\n"
                                                   >> 610   << "    ===================================================\n"
                                                   >> 611   << " Solid type: G4Tet\n"
                                                   >> 612   << " Parameters: \n"
                                                   >> 613   << "    anchor: " << fAnchor/mm << " mm \n"
                                                   >> 614   << "    p2: " << fP2/mm << " mm \n"
                                                   >> 615   << "    p3: " << fP3/mm << " mm \n"
                                                   >> 616   << "    p4: " << fP4/mm << " mm \n"
                                                   >> 617   << "    normal123: " << fNormal123 << " \n"
                                                   >> 618   << "    normal134: " << fNormal134 << " \n"
                                                   >> 619   << "    normal142: " << fNormal142 << " \n"
                                                   >> 620   << "    normal234: " << fNormal234 << " \n"
                                                   >> 621   << "-----------------------------------------------------------\n";
                                                   >> 622 
                                                   >> 623   return os;
597 }                                                 624 }
598                                                   625 
                                                   >> 626 
599 //////////////////////////////////////////////    627 ////////////////////////////////////////////////////////////////////////
600 //                                                628 //
601 // Make a clone of the object                  << 629 // GetPointOnFace
602 //                                                630 //
603 G4VSolid* G4Tet::Clone() const                 << 631 // Auxiliary method for get point on surface
                                                   >> 632 
                                                   >> 633 G4ThreeVector G4Tet::GetPointOnFace(G4ThreeVector p1, G4ThreeVector p2,
                                                   >> 634                                     G4ThreeVector p3, G4double& area) const
604 {                                                 635 {
605   return new G4Tet(*this);                     << 636   G4double lambda1,lambda2;
                                                   >> 637   G4ThreeVector v, w;
                                                   >> 638 
                                                   >> 639   v = p3 - p1;
                                                   >> 640   w = p1 - p2;
                                                   >> 641 
                                                   >> 642   lambda1 = RandFlat::shoot(0.,1.);
                                                   >> 643   lambda2 = RandFlat::shoot(0.,lambda1);
                                                   >> 644 
                                                   >> 645   area = 0.5*(v.cross(w)).mag();
                                                   >> 646 
                                                   >> 647   return (p2 + lambda1*w + lambda2*v);
606 }                                                 648 }
607                                                   649 
608 ////////////////////////////////////////////// << 650 ////////////////////////////////////////////////////////////////////////////
609 //                                                651 //
610 // Stream object contents to an output stream  << 652 // GetPointOnSurface
611 //                                             << 653 
612 std::ostream& G4Tet::StreamInfo(std::ostream&  << 654 G4ThreeVector G4Tet::GetPointOnSurface() const
613 {                                                 655 {
614   G4long oldprc = os.precision(16);            << 656   G4double chose,aOne,aTwo,aThree,aFour;
615   os << "------------------------------------- << 657   G4ThreeVector p1, p2, p3, p4;
616      << "    *** Dump for solid - " << GetName << 658   
617      << "    ================================= << 659   p1 = GetPointOnFace(fAnchor,fP2,fP3,aOne);
618      << " Solid type: " << GetEntityType() <<  << 660   p2 = GetPointOnFace(fAnchor,fP4,fP3,aTwo);
619      << " Parameters: \n"                      << 661   p3 = GetPointOnFace(fAnchor,fP4,fP2,aThree);
620      << "    anchor: " << fVertex[0]/mm << " m << 662   p4 = GetPointOnFace(fP4,fP3,fP2,aFour);
621      << "    p1    : " << fVertex[1]/mm << " m << 663   
622      << "    p2    : " << fVertex[2]/mm << " m << 664   chose = RandFlat::shoot(0.,aOne+aTwo+aThree+aFour);
623      << "    p3    : " << fVertex[3]/mm << " m << 665   if( (chose>=0.) && (chose <aOne) ) {return p1;}
624      << "------------------------------------- << 666   else if( (chose>=aOne) && (chose < aOne+aTwo) ) {return p2;}
625   os.precision(oldprc);                        << 667   else if( (chose>=aOne+aTwo) && (chose<aOne+aTwo+aThree) ) {return p3;}
626   return os;                                   << 668   return p4;
627 }                                                 669 }
628                                                   670 
629 //////////////////////////////////////////////    671 ////////////////////////////////////////////////////////////////////////
630 //                                                672 //
631 // Return random point on the surface          << 673 // GetVertices
632 //                                             << 
633 G4ThreeVector G4Tet::GetPointOnSurface() const << 
634 {                                              << 
635   constexpr G4int iface[4][3] = { {0,1,2}, {0, << 
636                                                   674 
637   // Select face                               << 675 std::vector<G4ThreeVector> G4Tet::GetVertices() const 
638   G4double select = fSurfaceArea*G4QuickRand() << 676 {
639   G4int i = 0;                                 << 677   std::vector<G4ThreeVector> vertices(4);
640   i += (G4int)(select > fArea[0]);             << 678   vertices[0] = fAnchor;
641   i += (G4int)(select > fArea[0] + fArea[1]);  << 679   vertices[1] = fP2;
642   i += (G4int)(select > fArea[0] + fArea[1] +  << 680   vertices[2] = fP3;
643                                                << 681   vertices[3] = fP4;
644   // Set selected triangle                     << 
645   G4ThreeVector p0 = fVertex[iface[i][0]];     << 
646   G4ThreeVector e1 = fVertex[iface[i][1]] - p0 << 
647   G4ThreeVector e2 = fVertex[iface[i][2]] - p0 << 
648                                                   682 
649   // Return random point                       << 683   return vertices;
650   G4double r1 = G4QuickRand();                 << 
651   G4double r2 = G4QuickRand();                 << 
652   return (r1 + r2 > 1.) ?                      << 
653     p0 + e1*(1. - r1) + e2*(1. - r2) : p0 + e1 << 
654 }                                                 684 }
655                                                   685 
656 //////////////////////////////////////////////    686 ////////////////////////////////////////////////////////////////////////
657 //                                                687 //
658 // Return volume of the tetrahedron            << 688 // GetCubicVolume
659 //                                             << 689 
660 G4double G4Tet::GetCubicVolume()                  690 G4double G4Tet::GetCubicVolume()
661 {                                                 691 {
662   return fCubicVolume;                            692   return fCubicVolume;
663 }                                                 693 }
664                                                   694 
665 //////////////////////////////////////////////    695 ////////////////////////////////////////////////////////////////////////
666 //                                                696 //
667 // Return surface area of the tetrahedron      << 697 // GetSurfaceArea
668 //                                             << 698 
669 G4double G4Tet::GetSurfaceArea()                  699 G4double G4Tet::GetSurfaceArea()
670 {                                                 700 {
671   return fSurfaceArea;                            701   return fSurfaceArea;
672 }                                                 702 }
673                                                   703 
674 ////////////////////////////////////////////// << 
675 //                                             << 
676 // Methods for visualisation                      704 // Methods for visualisation
                                                   >> 705 
                                                   >> 706 ////////////////////////////////////////////////////////////////////////
677 //                                                707 //
678 void G4Tet::DescribeYourselfTo (G4VGraphicsSce << 708 // DescribeYourselfTo
                                                   >> 709 
                                                   >> 710 void G4Tet::DescribeYourselfTo (G4VGraphicsScene& scene) const 
679 {                                                 711 {
680   scene.AddSolid (*this);                         712   scene.AddSolid (*this);
681 }                                                 713 }
682                                                   714 
683 //////////////////////////////////////////////    715 ////////////////////////////////////////////////////////////////////////
684 //                                                716 //
685 // Return VisExtent                            << 717 // GetExtent
686 //                                             << 718 
687 G4VisExtent G4Tet::GetExtent() const           << 719 G4VisExtent G4Tet::GetExtent() const 
688 {                                                 720 {
689   return { fBmin.x(), fBmax.x(),               << 721   return G4VisExtent (fXMin, fXMax, fYMin, fYMax, fZMin, fZMax);
690            fBmin.y(), fBmax.y(),               << 
691            fBmin.z(), fBmax.z() };             << 
692 }                                                 722 }
693                                                   723 
694 //////////////////////////////////////////////    724 ////////////////////////////////////////////////////////////////////////
695 //                                                725 //
696 // CreatePolyhedron                               726 // CreatePolyhedron
697 //                                             << 
698 G4Polyhedron* G4Tet::CreatePolyhedron() const  << 
699 {                                              << 
700   // Check orientation of vertices             << 
701   G4ThreeVector v1 = fVertex[1] - fVertex[0];  << 
702   G4ThreeVector v2 = fVertex[2] - fVertex[0];  << 
703   G4ThreeVector v3 = fVertex[3] - fVertex[0];  << 
704   G4bool invert = v1.cross(v2).dot(v3) < 0.;   << 
705   G4int k2 = (invert) ? 3 : 2;                 << 
706   G4int k3 = (invert) ? 2 : 3;                 << 
707                                                   727 
708   // Set coordinates of vertices               << 728 G4Polyhedron* G4Tet::CreatePolyhedron () const 
                                                   >> 729 {
                                                   >> 730   G4Polyhedron *ph=new G4Polyhedron;
709   G4double xyz[4][3];                             731   G4double xyz[4][3];
710   for (G4int i = 0; i < 3; ++i)                << 732   static G4int faces[4][4]={{1,3,2,0},{1,4,3,0},{1,2,4,0},{2,3,4,0}};
711   {                                            << 733   xyz[0][0]=fAnchor.x(); xyz[0][1]=fAnchor.y(); xyz[0][2]=fAnchor.z();
712     xyz[0][i] = fVertex[0][i];                 << 734   xyz[1][0]=fP2.x(); xyz[1][1]=fP2.y(); xyz[1][2]=fP2.z();
713     xyz[1][i] = fVertex[1][i];                 << 735   xyz[2][0]=fP3.x(); xyz[2][1]=fP3.y(); xyz[2][2]=fP3.z();
714     xyz[2][i] = fVertex[k2][i];                << 736   xyz[3][0]=fP4.x(); xyz[3][1]=fP4.y(); xyz[3][2]=fP4.z();
715     xyz[3][i] = fVertex[k3][i];                << 
716   }                                            << 
717                                                   737 
718   // Create polyhedron                         << 
719   G4int faces[4][4] = { {1,3,2,0}, {1,4,3,0},  << 
720   auto  ph = new G4Polyhedron;                 << 
721   ph->createPolyhedron(4,4,xyz,faces);            738   ph->createPolyhedron(4,4,xyz,faces);
722                                                   739 
723   return ph;                                      740   return ph;
724 }                                                 741 }
725                                                   742 
726 //////////////////////////////////////////////    743 ////////////////////////////////////////////////////////////////////////
727 //                                                744 //
728 // GetPolyhedron                               << 745 // CreateNURBS
                                                   >> 746 
                                                   >> 747 G4NURBS* G4Tet::CreateNURBS () const 
                                                   >> 748 {
                                                   >> 749   return new G4NURBSbox (fDx, fDy, fDz);
                                                   >> 750 }
                                                   >> 751 
                                                   >> 752 ////////////////////////////////////////////////////////////////////////
729 //                                                753 //
730 G4Polyhedron* G4Tet::GetPolyhedron() const     << 754 // GetPolyhedron
                                                   >> 755 
                                                   >> 756 G4Polyhedron* G4Tet::GetPolyhedron () const
731 {                                                 757 {
732   if (fpPolyhedron == nullptr ||               << 758   if (!fpPolyhedron ||
733       fRebuildPolyhedron ||                    << 
734       fpPolyhedron->GetNumberOfRotationStepsAt    759       fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
735       fpPolyhedron->GetNumberOfRotationSteps()    760       fpPolyhedron->GetNumberOfRotationSteps())
736   {                                            << 761     {
737     G4AutoLock l(&polyhedronMutex);            << 762       delete fpPolyhedron;
738     delete fpPolyhedron;                       << 763       fpPolyhedron = CreatePolyhedron();
739     fpPolyhedron = CreatePolyhedron();         << 764     }
740     fRebuildPolyhedron = false;                << 
741     l.unlock();                                << 
742   }                                            << 
743   return fpPolyhedron;                            765   return fpPolyhedron;
744 }                                                 766 }
745                                                << 
746 #endif                                         << 
747                                                   767