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.0.p4)


  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             << "Height in P0-P1-P2 = " << h1 << G4endl
134             << "  in triangle P0-P1-P2 = " <<  << 135             << "Height in P1-P2-P3 = " << h2 << G4endl
135             << "  in triangle P1-P2-P3 = " <<  << 136             << "Height in P2-P3-P4 = " << h3 << G4endl
136             << "  in triangle P2-P3-P0 = " <<  << 137             << "Height in P4-P0-P1 = " << h4;
137             << "  in triangle P3-P0-P1 = " <<  << 
138     G4Exception("G4QuadrangularFacet::G4Quadra    138     G4Exception("G4QuadrangularFacet::G4QuadrangularFacet()",
139           "GeomSolids1001", JustWarning, messa    139           "GeomSolids1001", JustWarning, message);
140     return;                                       140     return;
141   }                                               141   }
142                                                   142 
143   // Check that vertices are coplanar by compu    143   // Check that vertices are coplanar by computing minimal
144   // height of tetrahedron comprising of verti    144   // height of tetrahedron comprising of vertices
145   //                                              145   //
146   G4double smax = std::max( std::max(s1,s2), s    146   G4double smax = std::max( std::max(s1,s2), std::max(s3,s4) ); 
147   G4double hmin = 0.5 * std::fabs( e1.dot(e2.c    147   G4double hmin = 0.5 * std::fabs( e1.dot(e2.cross(e3)) ) / smax;
148   if (hmin >= epsilon)                            148   if (hmin >= epsilon)
149   {                                               149   {
150     ostringstream message;                        150     ostringstream message;
151     message << "Facet is not planar." << G4end    151     message << "Facet is not planar." << G4endl
152             << "Disrepancy = " << hmin << G4en    152             << "Disrepancy = " << hmin << G4endl
153             << "P0 = " << GetVertex(0) << G4en    153             << "P0 = " << GetVertex(0) << G4endl
154             << "P1 = " << GetVertex(1) << G4en    154             << "P1 = " << GetVertex(1) << G4endl
155             << "P2 = " << GetVertex(2) << G4en    155             << "P2 = " << GetVertex(2) << G4endl
156             << "P3 = " << GetVertex(3);           156             << "P3 = " << GetVertex(3);
157     G4Exception("G4QuadrangularFacet::G4Quadra    157     G4Exception("G4QuadrangularFacet::G4QuadrangularFacet()",
158           "GeomSolids1001", JustWarning, messa    158           "GeomSolids1001", JustWarning, message);
159     return;                                       159     return;
160   }                                               160   }
161                                                   161   
162   // Check that facet is convex by computing c    162   // Check that facet is convex by computing crosspoint 
163   // of diagonals                                 163   // of diagonals
164   //                                              164   //
165   G4ThreeVector normal = e2.cross(e3-e1);         165   G4ThreeVector normal = e2.cross(e3-e1);
166   G4double s = kInfinity, t = kInfinity, magni    166   G4double s = kInfinity, t = kInfinity, magnitude2 = normal.mag2();
167   if (magnitude2 > delta*delta) // check: magn    167   if (magnitude2 > delta*delta) // check: magnitude2 != 0.
168   {                                               168   {
169     s = normal.dot(e1.cross(e3-e1)) / magnitud    169     s = normal.dot(e1.cross(e3-e1)) / magnitude2;
170     t = normal.dot(e1.cross(e2))    / magnitud    170     t = normal.dot(e1.cross(e2))    / magnitude2;
171   }                                               171   }
172   if (s <= 0. || s >= 1. || t <= 0. || t >= 1.    172   if (s <= 0. || s >= 1. || t <= 0. || t >= 1.)
173   {                                               173   {
174     ostringstream message;                        174     ostringstream message;
175     message << "Facet is not convex." << G4end    175     message << "Facet is not convex." << G4endl
176             << "Parameters of crosspoint of di    176             << "Parameters of crosspoint of diagonals: "
177             << s << " and " << t << G4endl        177             << s << " and " << t << G4endl
178             << "should both be within (0,1) ra    178             << "should both be within (0,1) range" << G4endl
179       << "P0 = " << GetVertex(0) << G4endl        179       << "P0 = " << GetVertex(0) << G4endl
180       << "P1 = " << GetVertex(1) << G4endl        180       << "P1 = " << GetVertex(1) << G4endl
181       << "P2 = " << GetVertex(2) << G4endl        181       << "P2 = " << GetVertex(2) << G4endl
182       << "P3 = " << GetVertex(3);                 182       << "P3 = " << GetVertex(3);
183     G4Exception("G4QuadrangularFacet::G4Quadra    183     G4Exception("G4QuadrangularFacet::G4QuadrangularFacet()",
184           "GeomSolids1001", JustWarning, messa    184           "GeomSolids1001", JustWarning, message);
185     return;                                       185     return;
186   }                                               186   }
187                                                   187 
188   // Define facet                                 188   // Define facet
189   //                                              189   //
190   fFacet1 = G4TriangularFacet(GetVertex(0),Get    190   fFacet1 = G4TriangularFacet(GetVertex(0),GetVertex(1),GetVertex(2),ABSOLUTE);
191   fFacet2 = G4TriangularFacet(GetVertex(0),Get    191   fFacet2 = G4TriangularFacet(GetVertex(0),GetVertex(2),GetVertex(3),ABSOLUTE);
192                                                   192 
193   normal = normal.unit();                         193   normal = normal.unit();
194   fFacet1.SetSurfaceNormal(normal);               194   fFacet1.SetSurfaceNormal(normal);
195   fFacet2.SetSurfaceNormal(normal);               195   fFacet2.SetSurfaceNormal(normal);
196                                                   196   
197   G4ThreeVector vtmp = 0.5 * (e1 + e2);           197   G4ThreeVector vtmp = 0.5 * (e1 + e2);
198   fCircumcentre = GetVertex(0) + vtmp;            198   fCircumcentre = GetVertex(0) + vtmp;
199   G4double radiusSqr = vtmp.mag2();               199   G4double radiusSqr = vtmp.mag2();
200   fRadius = std::sqrt(radiusSqr);                 200   fRadius = std::sqrt(radiusSqr);
201   // 29.02.2016 Remark by E.Tcherniaev: comput    201   // 29.02.2016 Remark by E.Tcherniaev: computation
202   // of fCircumcenter and fRadius is wrong, ho    202   // of fCircumcenter and fRadius is wrong, however
203   // it did not create any problem till now.      203   // it did not create any problem till now.
204   // Bizarre! Need to investigate!                204   // Bizarre! Need to investigate!
205 }                                                 205 }
206                                                   206 
207 //////////////////////////////////////////////    207 ///////////////////////////////////////////////////////////////////////////////
208 //                                                208 //
209 G4QuadrangularFacet::~G4QuadrangularFacet () = << 209 G4QuadrangularFacet::~G4QuadrangularFacet ()
                                                   >> 210 {
                                                   >> 211 }
210                                                   212 
211 //////////////////////////////////////////////    213 ///////////////////////////////////////////////////////////////////////////////
212 //                                                214 //
213 G4QuadrangularFacet::G4QuadrangularFacet (cons    215 G4QuadrangularFacet::G4QuadrangularFacet (const G4QuadrangularFacet& rhs)
214   : G4VFacet(rhs)                                 216   : G4VFacet(rhs)
215 {                                                 217 {
216   fFacet1 = rhs.fFacet1;                          218   fFacet1 = rhs.fFacet1;
217   fFacet2 = rhs.fFacet2;                          219   fFacet2 = rhs.fFacet2;
218   fRadius = 0.0;                                  220   fRadius = 0.0;
219 }                                                 221 }
220                                                   222 
221 //////////////////////////////////////////////    223 ///////////////////////////////////////////////////////////////////////////////
222 //                                                224 //
223 G4QuadrangularFacet &                             225 G4QuadrangularFacet &
224 G4QuadrangularFacet::operator=(const G4Quadran    226 G4QuadrangularFacet::operator=(const G4QuadrangularFacet& rhs)
225 {                                                 227 {
226   if (this == &rhs)  return *this;                228   if (this == &rhs)  return *this;
227                                                   229 
228   fFacet1 = rhs.fFacet1;                          230   fFacet1 = rhs.fFacet1;
229   fFacet2 = rhs.fFacet2;                          231   fFacet2 = rhs.fFacet2;
230   fRadius = 0.0;                                  232   fRadius = 0.0;
231                                                   233 
232   return *this;                                   234   return *this;
233 }                                                 235 }
234                                                   236 
235 //////////////////////////////////////////////    237 ///////////////////////////////////////////////////////////////////////////////
236 //                                                238 //
237 G4VFacet* G4QuadrangularFacet::GetClone ()        239 G4VFacet* G4QuadrangularFacet::GetClone ()
238 {                                                 240 {
239   auto c = new G4QuadrangularFacet (GetVertex( << 241   G4QuadrangularFacet *c = new G4QuadrangularFacet (GetVertex(0), GetVertex(1),
240                                     GetVertex( << 242                                                     GetVertex(2), GetVertex(3),
                                                   >> 243                                                     ABSOLUTE);
241   return c;                                       244   return c;
242 }                                                 245 }
243                                                   246 
244 //////////////////////////////////////////////    247 ///////////////////////////////////////////////////////////////////////////////
245 //                                                248 //
246 G4ThreeVector G4QuadrangularFacet::Distance (c    249 G4ThreeVector G4QuadrangularFacet::Distance (const G4ThreeVector& p)
247 {                                                 250 {
248   G4ThreeVector v1 = fFacet1.Distance(p);         251   G4ThreeVector v1 = fFacet1.Distance(p);
249   G4ThreeVector v2 = fFacet2.Distance(p);         252   G4ThreeVector v2 = fFacet2.Distance(p);
250                                                   253 
251   if (v1.mag2() < v2.mag2()) return v1;           254   if (v1.mag2() < v2.mag2()) return v1;
252   else return v2;                                 255   else return v2;
253 }                                                 256 }
254                                                   257 
255 //////////////////////////////////////////////    258 ///////////////////////////////////////////////////////////////////////////////
256 //                                                259 //
257 G4double G4QuadrangularFacet::Distance (const     260 G4double G4QuadrangularFacet::Distance (const G4ThreeVector& p,
258                                                   261                                               G4double)
259 {                                                 262 {  
260   G4double dist = Distance(p).mag();              263   G4double dist = Distance(p).mag();
261   return dist;                                    264   return dist;
262 }                                                 265 }
263                                                   266 
264 //////////////////////////////////////////////    267 ///////////////////////////////////////////////////////////////////////////////
265 //                                                268 //
266 G4double G4QuadrangularFacet::Distance (const     269 G4double G4QuadrangularFacet::Distance (const G4ThreeVector& p, G4double,
267                                         const     270                                         const G4bool outgoing)
268 {                                                 271 {
269   G4double dist;                                  272   G4double dist;
270                                                   273 
271   G4ThreeVector v = Distance(p);                  274   G4ThreeVector v = Distance(p);
272   G4double dir = v.dot(GetSurfaceNormal());       275   G4double dir = v.dot(GetSurfaceNormal());
273   if ( ((dir > dirTolerance) && (!outgoing))      276   if ( ((dir > dirTolerance) && (!outgoing))
274     || ((dir < -dirTolerance) && outgoing))       277     || ((dir < -dirTolerance) && outgoing))
275     dist = kInfinity;                             278     dist = kInfinity;
276   else                                            279   else 
277     dist = v.mag();                               280     dist = v.mag();
278   return dist;                                    281   return dist;
279 }                                                 282 }
280                                                   283 
281 //////////////////////////////////////////////    284 ///////////////////////////////////////////////////////////////////////////////
282 //                                                285 //
283 G4double G4QuadrangularFacet::Extent (const G4    286 G4double G4QuadrangularFacet::Extent (const G4ThreeVector axis)
284 {                                                 287 {
285   G4double ss  = 0;                               288   G4double ss  = 0;
286                                                   289 
287   for (G4int i = 0; i <= 3; ++i)                  290   for (G4int i = 0; i <= 3; ++i)
288   {                                               291   {
289     G4double sp = GetVertex(i).dot(axis);         292     G4double sp = GetVertex(i).dot(axis);
290     if (sp > ss) ss = sp;                         293     if (sp > ss) ss = sp;
291   }                                               294   }
292   return ss;                                      295   return ss;
293 }                                                 296 }
294                                                   297 
295 //////////////////////////////////////////////    298 ///////////////////////////////////////////////////////////////////////////////
296 //                                                299 //
297 G4bool G4QuadrangularFacet::Intersect (const G    300 G4bool G4QuadrangularFacet::Intersect (const G4ThreeVector& p,
298                                        const G    301                                        const G4ThreeVector& v,
299                                              G    302                                              G4bool outgoing,
300                                              G    303                                              G4double& distance,
301                                              G    304                                              G4double& distFromSurface,
302                                              G    305                                              G4ThreeVector& normal)
303 {                                                 306 {
304   G4bool intersect =                              307   G4bool intersect =
305     fFacet1.Intersect(p,v,outgoing,distance,di    308     fFacet1.Intersect(p,v,outgoing,distance,distFromSurface,normal);
306   if (!intersect) intersect =                     309   if (!intersect) intersect =
307     fFacet2.Intersect(p,v,outgoing,distance,di    310     fFacet2.Intersect(p,v,outgoing,distance,distFromSurface,normal);
308   if (!intersect)                                 311   if (!intersect)
309   {                                               312   {
310     distance = distFromSurface = kInfinity;       313     distance = distFromSurface = kInfinity;
311     normal.set(0,0,0);                            314     normal.set(0,0,0);
312   }                                               315   }
313   return intersect;                               316   return intersect;
314 }                                                 317 }
315                                                   318 
316 //////////////////////////////////////////////    319 ///////////////////////////////////////////////////////////////////////////////
317 //                                                320 //
318 // Auxiliary method to get a uniform random po    321 // Auxiliary method to get a uniform random point on the facet
319 //                                                322 //
320 G4ThreeVector G4QuadrangularFacet::GetPointOnF    323 G4ThreeVector G4QuadrangularFacet::GetPointOnFace() const
321 {                                                 324 {
322   G4double s1 = fFacet1.GetArea();                325   G4double s1 = fFacet1.GetArea();
323   G4double s2 = fFacet2.GetArea();                326   G4double s2 = fFacet2.GetArea();
324   return ((s1+s2)*G4UniformRand() < s1) ?         327   return ((s1+s2)*G4UniformRand() < s1) ?
325     fFacet1.GetPointOnFace() : fFacet2.GetPoin    328     fFacet1.GetPointOnFace() : fFacet2.GetPointOnFace();
326 }                                                 329 }
327                                                   330 
328 //////////////////////////////////////////////    331 ///////////////////////////////////////////////////////////////////////////////
329 //                                                332 //
330 // Auxiliary method for returning the surface     333 // Auxiliary method for returning the surface area
331 //                                                334 //
332 G4double G4QuadrangularFacet::GetArea() const     335 G4double G4QuadrangularFacet::GetArea() const
333 {                                                 336 {
334   G4double area = fFacet1.GetArea() + fFacet2.    337   G4double area = fFacet1.GetArea() + fFacet2.GetArea();
335   return area;                                    338   return area;
336 }                                                 339 }
337                                                   340 
338 //////////////////////////////////////////////    341 ///////////////////////////////////////////////////////////////////////////////
339 //                                                342 //
340 G4String G4QuadrangularFacet::GetEntityType ()    343 G4String G4QuadrangularFacet::GetEntityType () const
341 {                                                 344 {
342   return "G4QuadrangularFacet";                   345   return "G4QuadrangularFacet";
343 }                                                 346 }
344                                                   347 
345 //////////////////////////////////////////////    348 ///////////////////////////////////////////////////////////////////////////////
346 //                                                349 //
347 G4ThreeVector G4QuadrangularFacet::GetSurfaceN    350 G4ThreeVector G4QuadrangularFacet::GetSurfaceNormal () const
348 {                                                 351 {
349   return fFacet1.GetSurfaceNormal();              352   return fFacet1.GetSurfaceNormal();
350 }                                                 353 }
351                                                   354