Geant4 Cross Reference

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


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration and of QinetiQ Ltd,   *
 20 // * subject to DEFCON 705 IPR conditions.         20 // * subject to DEFCON 705 IPR conditions.                            *
 21 // * By using,  copying,  modifying or  distri     21 // * By using,  copying,  modifying or  distributing the software (or *
 22 // * any work based  on the software)  you  ag     22 // * any work based  on the software)  you  agree  to acknowledge its *
 23 // * use  in  resulting  scientific  publicati     23 // * use  in  resulting  scientific  publications,  and indicate your *
 24 // * acceptance of all terms of the Geant4 Sof     24 // * acceptance of all terms of the Geant4 Software license.          *
 25 // *******************************************     25 // ********************************************************************
 26 //                                                 26 //
 27 // G4QuadrangularFacet class implementation.       27 // G4QuadrangularFacet class implementation.
 28 //                                                 28 //
 29 // 31 October 2004 P R Truscott, QinetiQ Ltd,      29 // 31 October 2004 P R Truscott, QinetiQ Ltd, UK - Created.
 30 // 12 October 2012 M Gayer, CERN                   30 // 12 October 2012 M Gayer, CERN
 31 //                 New implementation reducing     31 //                 New implementation reducing memory requirements by 50%,
 32 //                 and considerable CPU speedu     32 //                 and considerable CPU speedup together with the new
 33 //                 implementation of G4Tessell     33 //                 implementation of G4TessellatedSolid.
 34 // 29 February 2016 E Tcherniaev, CERN             34 // 29 February 2016 E Tcherniaev, CERN
 35 //                 Added exhaustive tests to c     35 //                 Added exhaustive tests to catch various problems with a
 36 //                 quadrangular facet: colline     36 //                 quadrangular facet: collinear vertices, non planar surface,
 37 //                 degenerate, concave or self     37 //                 degenerate, concave or self intersecting quadrilateral.
 38 // -------------------------------------------     38 // --------------------------------------------------------------------
 39                                                    39 
 40 #include "G4QuadrangularFacet.hh"                  40 #include "G4QuadrangularFacet.hh"
 41 #include "geomdefs.hh"                             41 #include "geomdefs.hh"
 42 #include "Randomize.hh"                            42 #include "Randomize.hh"
 43                                                    43  
 44 using namespace std;                               44 using namespace std;
 45                                                    45 
 46 //////////////////////////////////////////////     46 ///////////////////////////////////////////////////////////////////////////////
 47 //                                                 47 //
 48 // Constructing two adjacent G4TriangularFacet     48 // Constructing two adjacent G4TriangularFacet
 49 // Not efficient, but practical...                 49 // Not efficient, but practical...
 50 //                                                 50 //
 51 G4QuadrangularFacet::G4QuadrangularFacet (cons     51 G4QuadrangularFacet::G4QuadrangularFacet (const G4ThreeVector& vt0,
 52                                           cons     52                                           const G4ThreeVector& vt1,
 53                                           cons     53                                           const G4ThreeVector& vt2,
 54                                           cons     54                                           const G4ThreeVector& vt3,
 55                                                    55                                                 G4FacetVertexType vertexType)
                                                   >>  56   : G4VFacet()
 56 {                                                  57 {
 57   G4double delta   =  1.0 * kCarTolerance; //      58   G4double delta   =  1.0 * kCarTolerance; // dimension tolerance
 58   G4double epsilon = 0.01 * kCarTolerance; //      59   G4double epsilon = 0.01 * kCarTolerance; // planarity tolerance
 59                                                    60 
 60   G4ThreeVector e1, e2, e3;                        61   G4ThreeVector e1, e2, e3;
 61   SetVertex(0, vt0);                               62   SetVertex(0, vt0);
 62   if (vertexType == ABSOLUTE)                      63   if (vertexType == ABSOLUTE)
 63   {                                                64   {
 64     SetVertex(1, vt1);                             65     SetVertex(1, vt1);
 65     SetVertex(2, vt2);                             66     SetVertex(2, vt2);
 66     SetVertex(3, vt3);                             67     SetVertex(3, vt3);
 67                                                    68 
 68     e1 = vt1 - vt0;                                69     e1 = vt1 - vt0;
 69     e2 = vt2 - vt0;                                70     e2 = vt2 - vt0;
 70     e3 = vt3 - vt0;                                71     e3 = vt3 - vt0;
 71   }                                                72   }
 72   else                                             73   else
 73   {                                                74   {
 74     SetVertex(1, vt0 + vt1);                       75     SetVertex(1, vt0 + vt1);
 75     SetVertex(2, vt0 + vt2);                       76     SetVertex(2, vt0 + vt2);
 76     SetVertex(3, vt0 + vt3);                       77     SetVertex(3, vt0 + vt3);
 77                                                    78 
 78     e1 = vt1;                                      79     e1 = vt1;
 79     e2 = vt2;                                      80     e2 = vt2;
 80     e3 = vt3;                                      81     e3 = vt3;
 81   }                                                82   }
 82                                                    83 
 83   // Check length of sides and diagonals           84   // Check length of sides and diagonals
 84   //                                               85   //
 85   G4double leng1 = e1.mag();                       86   G4double leng1 = e1.mag();
 86   G4double leng2 = (e2-e1).mag();                  87   G4double leng2 = (e2-e1).mag();
 87   G4double leng3 = (e3-e2).mag();                  88   G4double leng3 = (e3-e2).mag();
 88   G4double leng4 = e3.mag();                       89   G4double leng4 = e3.mag();
 89                                                    90 
 90   G4double diag1 = e2.mag();                       91   G4double diag1 = e2.mag();
 91   G4double diag2 = (e3-e1).mag();                  92   G4double diag2 = (e3-e1).mag();
 92                                                    93 
 93   if (leng1 <= delta || leng2 <= delta || leng     94   if (leng1 <= delta || leng2 <= delta || leng3 <= delta || leng4 <= delta ||
 94       diag1 <= delta || diag2 <= delta)            95       diag1 <= delta || diag2 <= delta)
 95   {                                                96   {
 96     ostringstream message;                         97     ostringstream message;
 97     message << "Sides/diagonals of facet are t     98     message << "Sides/diagonals of facet are too small." << G4endl
 98             << "P0 = " << GetVertex(0) << G4en     99             << "P0 = " << GetVertex(0) << G4endl
 99             << "P1 = " << GetVertex(1) << G4en    100             << "P1 = " << GetVertex(1) << G4endl
100             << "P2 = " << GetVertex(2) << G4en    101             << "P2 = " << GetVertex(2) << G4endl
101             << "P3 = " << GetVertex(3) << G4en    102             << "P3 = " << GetVertex(3) << G4endl
102             << "Side1 length (P0->P1) = " << l    103             << "Side1 length (P0->P1) = " << leng1 << G4endl
103             << "Side2 length (P1->P2) = " << l    104             << "Side2 length (P1->P2) = " << leng2 << G4endl
104             << "Side3 length (P2->P3) = " << l    105             << "Side3 length (P2->P3) = " << leng3 << G4endl
105             << "Side4 length (P3->P0) = " << l    106             << "Side4 length (P3->P0) = " << leng4 << G4endl
106             << "Diagonal1 length (P0->P2) = "     107             << "Diagonal1 length (P0->P2) = " << diag1 << G4endl
107             << "Diagonal2 length (P1->P3) = "     108             << "Diagonal2 length (P1->P3) = " << diag2;
108     G4Exception("G4QuadrangularFacet::G4Quadra    109     G4Exception("G4QuadrangularFacet::G4QuadrangularFacet()",
109                 "GeomSolids1001", JustWarning,    110                 "GeomSolids1001", JustWarning, message);
110     return;                                       111     return;
111   }                                               112   }
112                                                   113 
113   // Check that vertices are not collinear        114   // Check that vertices are not collinear
114   //                                              115   //
115   G4double s1 = (e1.cross(e2)).mag()*0.5;         116   G4double s1 = (e1.cross(e2)).mag()*0.5;
116   G4double s2 = ((e2-e1).cross(e3-e2)).mag()*0    117   G4double s2 = ((e2-e1).cross(e3-e2)).mag()*0.5;
117   G4double s3 = (e2.cross(e3)).mag()*0.5;         118   G4double s3 = (e2.cross(e3)).mag()*0.5;
118   G4double s4 = (e1.cross(e3)).mag()*0.5;         119   G4double s4 = (e1.cross(e3)).mag()*0.5;
119                                                   120 
120   G4double h1 = 2.*s1 / std::max(std::max(leng    121   G4double h1 = 2.*s1 / std::max(std::max(leng1,leng2),diag1);
121   G4double h2 = 2.*s2 / std::max(std::max(leng    122   G4double h2 = 2.*s2 / std::max(std::max(leng2,leng3),diag2);
122   G4double h3 = 2.*s3 / std::max(std::max(leng    123   G4double h3 = 2.*s3 / std::max(std::max(leng3,leng4),diag1);
123   G4double h4 = 2.*s4 / std::max(std::max(leng    124   G4double h4 = 2.*s4 / std::max(std::max(leng4,leng1),diag2);
124                                                   125    
125   if (h1 <= delta || h2 <= delta || h3 <= delt    126   if (h1 <= delta || h2 <= delta || h3 <= delta || h4 <= delta )
126   {                                               127   {
127     ostringstream message;                        128     ostringstream message;
128     message << "Facet has three or more collin    129     message << "Facet has three or more collinear vertices." << G4endl
129             << "P0 = " << GetVertex(0) << G4en    130             << "P0 = " << GetVertex(0) << G4endl
130             << "P1 = " << GetVertex(1) << G4en    131             << "P1 = " << GetVertex(1) << G4endl
131             << "P2 = " << GetVertex(2) << G4en    132             << "P2 = " << GetVertex(2) << G4endl
132             << "P3 = " << GetVertex(3) << G4en    133             << "P3 = " << GetVertex(3) << G4endl
133             << "Smallest heights:" << G4endl      134             << "Smallest heights:" << G4endl
134             << "  in triangle P0-P1-P2 = " <<     135             << "  in triangle P0-P1-P2 = " << h1 << G4endl
135             << "  in triangle P1-P2-P3 = " <<     136             << "  in triangle P1-P2-P3 = " << h2 << G4endl
136             << "  in triangle P2-P3-P0 = " <<     137             << "  in triangle P2-P3-P0 = " << h3 << G4endl
137             << "  in triangle P3-P0-P1 = " <<     138             << "  in triangle P3-P0-P1 = " << h4;
138     G4Exception("G4QuadrangularFacet::G4Quadra    139     G4Exception("G4QuadrangularFacet::G4QuadrangularFacet()",
139           "GeomSolids1001", JustWarning, messa    140           "GeomSolids1001", JustWarning, message);
140     return;                                       141     return;
141   }                                               142   }
142                                                   143 
143   // Check that vertices are coplanar by compu    144   // Check that vertices are coplanar by computing minimal
144   // height of tetrahedron comprising of verti    145   // height of tetrahedron comprising of vertices
145   //                                              146   //
146   G4double smax = std::max( std::max(s1,s2), s    147   G4double smax = std::max( std::max(s1,s2), std::max(s3,s4) ); 
147   G4double hmin = 0.5 * std::fabs( e1.dot(e2.c    148   G4double hmin = 0.5 * std::fabs( e1.dot(e2.cross(e3)) ) / smax;
148   if (hmin >= epsilon)                            149   if (hmin >= epsilon)
149   {                                               150   {
150     ostringstream message;                        151     ostringstream message;
151     message << "Facet is not planar." << G4end    152     message << "Facet is not planar." << G4endl
152             << "Disrepancy = " << hmin << G4en    153             << "Disrepancy = " << hmin << G4endl
153             << "P0 = " << GetVertex(0) << G4en    154             << "P0 = " << GetVertex(0) << G4endl
154             << "P1 = " << GetVertex(1) << G4en    155             << "P1 = " << GetVertex(1) << G4endl
155             << "P2 = " << GetVertex(2) << G4en    156             << "P2 = " << GetVertex(2) << G4endl
156             << "P3 = " << GetVertex(3);           157             << "P3 = " << GetVertex(3);
157     G4Exception("G4QuadrangularFacet::G4Quadra    158     G4Exception("G4QuadrangularFacet::G4QuadrangularFacet()",
158           "GeomSolids1001", JustWarning, messa    159           "GeomSolids1001", JustWarning, message);
159     return;                                       160     return;
160   }                                               161   }
161                                                   162   
162   // Check that facet is convex by computing c    163   // Check that facet is convex by computing crosspoint 
163   // of diagonals                                 164   // of diagonals
164   //                                              165   //
165   G4ThreeVector normal = e2.cross(e3-e1);         166   G4ThreeVector normal = e2.cross(e3-e1);
166   G4double s = kInfinity, t = kInfinity, magni    167   G4double s = kInfinity, t = kInfinity, magnitude2 = normal.mag2();
167   if (magnitude2 > delta*delta) // check: magn    168   if (magnitude2 > delta*delta) // check: magnitude2 != 0.
168   {                                               169   {
169     s = normal.dot(e1.cross(e3-e1)) / magnitud    170     s = normal.dot(e1.cross(e3-e1)) / magnitude2;
170     t = normal.dot(e1.cross(e2))    / magnitud    171     t = normal.dot(e1.cross(e2))    / magnitude2;
171   }                                               172   }
172   if (s <= 0. || s >= 1. || t <= 0. || t >= 1.    173   if (s <= 0. || s >= 1. || t <= 0. || t >= 1.)
173   {                                               174   {
174     ostringstream message;                        175     ostringstream message;
175     message << "Facet is not convex." << G4end    176     message << "Facet is not convex." << G4endl
176             << "Parameters of crosspoint of di    177             << "Parameters of crosspoint of diagonals: "
177             << s << " and " << t << G4endl        178             << s << " and " << t << G4endl
178             << "should both be within (0,1) ra    179             << "should both be within (0,1) range" << G4endl
179       << "P0 = " << GetVertex(0) << G4endl        180       << "P0 = " << GetVertex(0) << G4endl
180       << "P1 = " << GetVertex(1) << G4endl        181       << "P1 = " << GetVertex(1) << G4endl
181       << "P2 = " << GetVertex(2) << G4endl        182       << "P2 = " << GetVertex(2) << G4endl
182       << "P3 = " << GetVertex(3);                 183       << "P3 = " << GetVertex(3);
183     G4Exception("G4QuadrangularFacet::G4Quadra    184     G4Exception("G4QuadrangularFacet::G4QuadrangularFacet()",
184           "GeomSolids1001", JustWarning, messa    185           "GeomSolids1001", JustWarning, message);
185     return;                                       186     return;
186   }                                               187   }
187                                                   188 
188   // Define facet                                 189   // Define facet
189   //                                              190   //
190   fFacet1 = G4TriangularFacet(GetVertex(0),Get    191   fFacet1 = G4TriangularFacet(GetVertex(0),GetVertex(1),GetVertex(2),ABSOLUTE);
191   fFacet2 = G4TriangularFacet(GetVertex(0),Get    192   fFacet2 = G4TriangularFacet(GetVertex(0),GetVertex(2),GetVertex(3),ABSOLUTE);
192                                                   193 
193   normal = normal.unit();                         194   normal = normal.unit();
194   fFacet1.SetSurfaceNormal(normal);               195   fFacet1.SetSurfaceNormal(normal);
195   fFacet2.SetSurfaceNormal(normal);               196   fFacet2.SetSurfaceNormal(normal);
196                                                   197   
197   G4ThreeVector vtmp = 0.5 * (e1 + e2);           198   G4ThreeVector vtmp = 0.5 * (e1 + e2);
198   fCircumcentre = GetVertex(0) + vtmp;            199   fCircumcentre = GetVertex(0) + vtmp;
199   G4double radiusSqr = vtmp.mag2();               200   G4double radiusSqr = vtmp.mag2();
200   fRadius = std::sqrt(radiusSqr);                 201   fRadius = std::sqrt(radiusSqr);
201   // 29.02.2016 Remark by E.Tcherniaev: comput    202   // 29.02.2016 Remark by E.Tcherniaev: computation
202   // of fCircumcenter and fRadius is wrong, ho    203   // of fCircumcenter and fRadius is wrong, however
203   // it did not create any problem till now.      204   // it did not create any problem till now.
204   // Bizarre! Need to investigate!                205   // Bizarre! Need to investigate!
205 }                                                 206 }
206                                                   207 
207 //////////////////////////////////////////////    208 ///////////////////////////////////////////////////////////////////////////////
208 //                                                209 //
209 G4QuadrangularFacet::~G4QuadrangularFacet () = << 210 G4QuadrangularFacet::~G4QuadrangularFacet ()
                                                   >> 211 {
                                                   >> 212 }
210                                                   213 
211 //////////////////////////////////////////////    214 ///////////////////////////////////////////////////////////////////////////////
212 //                                                215 //
213 G4QuadrangularFacet::G4QuadrangularFacet (cons    216 G4QuadrangularFacet::G4QuadrangularFacet (const G4QuadrangularFacet& rhs)
214   : G4VFacet(rhs)                                 217   : G4VFacet(rhs)
215 {                                                 218 {
216   fFacet1 = rhs.fFacet1;                          219   fFacet1 = rhs.fFacet1;
217   fFacet2 = rhs.fFacet2;                          220   fFacet2 = rhs.fFacet2;
218   fRadius = 0.0;                                  221   fRadius = 0.0;
219 }                                                 222 }
220                                                   223 
221 //////////////////////////////////////////////    224 ///////////////////////////////////////////////////////////////////////////////
222 //                                                225 //
223 G4QuadrangularFacet &                             226 G4QuadrangularFacet &
224 G4QuadrangularFacet::operator=(const G4Quadran    227 G4QuadrangularFacet::operator=(const G4QuadrangularFacet& rhs)
225 {                                                 228 {
226   if (this == &rhs)  return *this;                229   if (this == &rhs)  return *this;
227                                                   230 
228   fFacet1 = rhs.fFacet1;                          231   fFacet1 = rhs.fFacet1;
229   fFacet2 = rhs.fFacet2;                          232   fFacet2 = rhs.fFacet2;
230   fRadius = 0.0;                                  233   fRadius = 0.0;
231                                                   234 
232   return *this;                                   235   return *this;
233 }                                                 236 }
234                                                   237 
235 //////////////////////////////////////////////    238 ///////////////////////////////////////////////////////////////////////////////
236 //                                                239 //
237 G4VFacet* G4QuadrangularFacet::GetClone ()        240 G4VFacet* G4QuadrangularFacet::GetClone ()
238 {                                                 241 {
239   auto c = new G4QuadrangularFacet (GetVertex( << 242   G4QuadrangularFacet *c = new G4QuadrangularFacet (GetVertex(0), GetVertex(1),
240                                     GetVertex( << 243                                                     GetVertex(2), GetVertex(3),
                                                   >> 244                                                     ABSOLUTE);
241   return c;                                       245   return c;
242 }                                                 246 }
243                                                   247 
244 //////////////////////////////////////////////    248 ///////////////////////////////////////////////////////////////////////////////
245 //                                                249 //
246 G4ThreeVector G4QuadrangularFacet::Distance (c    250 G4ThreeVector G4QuadrangularFacet::Distance (const G4ThreeVector& p)
247 {                                                 251 {
248   G4ThreeVector v1 = fFacet1.Distance(p);         252   G4ThreeVector v1 = fFacet1.Distance(p);
249   G4ThreeVector v2 = fFacet2.Distance(p);         253   G4ThreeVector v2 = fFacet2.Distance(p);
250                                                   254 
251   if (v1.mag2() < v2.mag2()) return v1;           255   if (v1.mag2() < v2.mag2()) return v1;
252   else return v2;                                 256   else return v2;
253 }                                                 257 }
254                                                   258 
255 //////////////////////////////////////////////    259 ///////////////////////////////////////////////////////////////////////////////
256 //                                                260 //
257 G4double G4QuadrangularFacet::Distance (const     261 G4double G4QuadrangularFacet::Distance (const G4ThreeVector& p,
258                                                   262                                               G4double)
259 {                                                 263 {  
260   G4double dist = Distance(p).mag();              264   G4double dist = Distance(p).mag();
261   return dist;                                    265   return dist;
262 }                                                 266 }
263                                                   267 
264 //////////////////////////////////////////////    268 ///////////////////////////////////////////////////////////////////////////////
265 //                                                269 //
266 G4double G4QuadrangularFacet::Distance (const     270 G4double G4QuadrangularFacet::Distance (const G4ThreeVector& p, G4double,
267                                         const     271                                         const G4bool outgoing)
268 {                                                 272 {
269   G4double dist;                                  273   G4double dist;
270                                                   274 
271   G4ThreeVector v = Distance(p);                  275   G4ThreeVector v = Distance(p);
272   G4double dir = v.dot(GetSurfaceNormal());       276   G4double dir = v.dot(GetSurfaceNormal());
273   if ( ((dir > dirTolerance) && (!outgoing))      277   if ( ((dir > dirTolerance) && (!outgoing))
274     || ((dir < -dirTolerance) && outgoing))       278     || ((dir < -dirTolerance) && outgoing))
275     dist = kInfinity;                             279     dist = kInfinity;
276   else                                            280   else 
277     dist = v.mag();                               281     dist = v.mag();
278   return dist;                                    282   return dist;
279 }                                                 283 }
280                                                   284 
281 //////////////////////////////////////////////    285 ///////////////////////////////////////////////////////////////////////////////
282 //                                                286 //
283 G4double G4QuadrangularFacet::Extent (const G4    287 G4double G4QuadrangularFacet::Extent (const G4ThreeVector axis)
284 {                                                 288 {
285   G4double ss  = 0;                               289   G4double ss  = 0;
286                                                   290 
287   for (G4int i = 0; i <= 3; ++i)                  291   for (G4int i = 0; i <= 3; ++i)
288   {                                               292   {
289     G4double sp = GetVertex(i).dot(axis);         293     G4double sp = GetVertex(i).dot(axis);
290     if (sp > ss) ss = sp;                         294     if (sp > ss) ss = sp;
291   }                                               295   }
292   return ss;                                      296   return ss;
293 }                                                 297 }
294                                                   298 
295 //////////////////////////////////////////////    299 ///////////////////////////////////////////////////////////////////////////////
296 //                                                300 //
297 G4bool G4QuadrangularFacet::Intersect (const G    301 G4bool G4QuadrangularFacet::Intersect (const G4ThreeVector& p,
298                                        const G    302                                        const G4ThreeVector& v,
299                                              G    303                                              G4bool outgoing,
300                                              G    304                                              G4double& distance,
301                                              G    305                                              G4double& distFromSurface,
302                                              G    306                                              G4ThreeVector& normal)
303 {                                                 307 {
304   G4bool intersect =                              308   G4bool intersect =
305     fFacet1.Intersect(p,v,outgoing,distance,di    309     fFacet1.Intersect(p,v,outgoing,distance,distFromSurface,normal);
306   if (!intersect) intersect =                     310   if (!intersect) intersect =
307     fFacet2.Intersect(p,v,outgoing,distance,di    311     fFacet2.Intersect(p,v,outgoing,distance,distFromSurface,normal);
308   if (!intersect)                                 312   if (!intersect)
309   {                                               313   {
310     distance = distFromSurface = kInfinity;       314     distance = distFromSurface = kInfinity;
311     normal.set(0,0,0);                            315     normal.set(0,0,0);
312   }                                               316   }
313   return intersect;                               317   return intersect;
314 }                                                 318 }
315                                                   319 
316 //////////////////////////////////////////////    320 ///////////////////////////////////////////////////////////////////////////////
317 //                                                321 //
318 // Auxiliary method to get a uniform random po    322 // Auxiliary method to get a uniform random point on the facet
319 //                                                323 //
320 G4ThreeVector G4QuadrangularFacet::GetPointOnF    324 G4ThreeVector G4QuadrangularFacet::GetPointOnFace() const
321 {                                                 325 {
322   G4double s1 = fFacet1.GetArea();                326   G4double s1 = fFacet1.GetArea();
323   G4double s2 = fFacet2.GetArea();                327   G4double s2 = fFacet2.GetArea();
324   return ((s1+s2)*G4UniformRand() < s1) ?         328   return ((s1+s2)*G4UniformRand() < s1) ?
325     fFacet1.GetPointOnFace() : fFacet2.GetPoin    329     fFacet1.GetPointOnFace() : fFacet2.GetPointOnFace();
326 }                                                 330 }
327                                                   331 
328 //////////////////////////////////////////////    332 ///////////////////////////////////////////////////////////////////////////////
329 //                                                333 //
330 // Auxiliary method for returning the surface     334 // Auxiliary method for returning the surface area
331 //                                                335 //
332 G4double G4QuadrangularFacet::GetArea() const     336 G4double G4QuadrangularFacet::GetArea() const
333 {                                                 337 {
334   G4double area = fFacet1.GetArea() + fFacet2.    338   G4double area = fFacet1.GetArea() + fFacet2.GetArea();
335   return area;                                    339   return area;
336 }                                                 340 }
337                                                   341 
338 //////////////////////////////////////////////    342 ///////////////////////////////////////////////////////////////////////////////
339 //                                                343 //
340 G4String G4QuadrangularFacet::GetEntityType ()    344 G4String G4QuadrangularFacet::GetEntityType () const
341 {                                                 345 {
342   return "G4QuadrangularFacet";                   346   return "G4QuadrangularFacet";
343 }                                                 347 }
344                                                   348 
345 //////////////////////////////////////////////    349 ///////////////////////////////////////////////////////////////////////////////
346 //                                                350 //
347 G4ThreeVector G4QuadrangularFacet::GetSurfaceN    351 G4ThreeVector G4QuadrangularFacet::GetSurfaceNormal () const
348 {                                                 352 {
349   return fFacet1.GetSurfaceNormal();              353   return fFacet1.GetSurfaceNormal();
350 }                                                 354 }
351                                                   355