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 10.0.p1)


  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.   << 
 28 //                                                 27 //
 29 // 31 October 2004 P R Truscott, QinetiQ Ltd,  <<  28 // $Id: G4QuadrangularFacet.cc 66819 2013-01-12 16:20:10Z gcosmo $
 30 // 12 October 2012 M Gayer, CERN               <<  29 //
 31 //                 New implementation reducing <<  30 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 32 //                 and considerable CPU speedu <<  31 //
 33 //                 implementation of G4Tessell <<  32 // CHANGE HISTORY
 34 // 29 February 2016 E Tcherniaev, CERN         <<  33 // --------------
 35 //                 Added exhaustive tests to c <<  34 //
 36 //                 quadrangular facet: colline <<  35 // 31 October 2004, P R Truscott, QinetiQ Ltd, UK - Created.
 37 //                 degenerate, concave or self <<  36 // 12 October 2012, M Gayer, CERN
 38 // ------------------------------------------- <<  37 //                  New implementation reducing memory requirements by 50%,
                                                   >>  38 //                  and considerable CPU speedup together with the new
                                                   >>  39 //                  implementation of G4TessellatedSolid.
                                                   >>  40 //
                                                   >>  41 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 39                                                    42 
 40 #include "G4QuadrangularFacet.hh"                  43 #include "G4QuadrangularFacet.hh"
 41 #include "geomdefs.hh"                             44 #include "geomdefs.hh"
 42 #include "Randomize.hh"                            45 #include "Randomize.hh"
 43                                                    46  
 44 using namespace std;                               47 using namespace std;
 45                                                    48 
 46 //////////////////////////////////////////////     49 ///////////////////////////////////////////////////////////////////////////////
 47 //                                                 50 //
 48 // Constructing two adjacent G4TriangularFacet <<  51 // !!!THIS IS A FUDGE!!!  IT'S TWO ADJACENT G4TRIANGULARFACETS
 49 // Not efficient, but practical...             <<  52 // --- NOT EFFICIENT BUT PRACTICAL.
 50 //                                                 53 //
 51 G4QuadrangularFacet::G4QuadrangularFacet (cons <<  54 G4QuadrangularFacet::G4QuadrangularFacet (const G4ThreeVector &vt0,
 52                                           cons <<  55                                           const G4ThreeVector &vt1,
 53                                           cons <<  56                                           const G4ThreeVector &vt2,
 54                                           cons <<  57                                           const G4ThreeVector &vt3,
 55                                                    58                                                 G4FacetVertexType vertexType)
 56 {                                                  59 {
 57   G4double delta   =  1.0 * kCarTolerance; //  << 
 58   G4double epsilon = 0.01 * kCarTolerance; //  << 
 59                                                << 
 60   G4ThreeVector e1, e2, e3;                        60   G4ThreeVector e1, e2, e3;
                                                   >>  61 
 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   G4double length1 = e1.mag();
 83   // Check length of sides and diagonals       <<  84   G4double length2 = (GetVertex(2)-GetVertex(1)).mag();
 84   //                                           <<  85   G4double length3 = (GetVertex(3)-GetVertex(2)).mag();
 85   G4double leng1 = e1.mag();                   <<  86   G4double length4 = e3.mag();
 86   G4double leng2 = (e2-e1).mag();              <<  87 
 87   G4double leng3 = (e3-e2).mag();              <<  88   G4ThreeVector normal1 = e1.cross(e2).unit();
 88   G4double leng4 = e3.mag();                   <<  89   G4ThreeVector normal2 = e2.cross(e3).unit(); 
 89                                                <<  90 
 90   G4double diag1 = e2.mag();                   <<  91   bool isDefined = (length1 > kCarTolerance && length2 > kCarTolerance &&
 91   G4double diag2 = (e3-e1).mag();              <<  92     length3 > kCarTolerance && length4 > kCarTolerance &&
 92                                                <<  93     normal1.dot(normal2) >= 0.9999999999);
 93   if (leng1 <= delta || leng2 <= delta || leng <<  94 
 94       diag1 <= delta || diag2 <= delta)        <<  95   if (isDefined)
 95   {                                            <<  96   {
 96     ostringstream message;                     <<  97     fFacet1 = G4TriangularFacet (GetVertex(0),GetVertex(1),
 97     message << "Sides/diagonals of facet are t <<  98                                  GetVertex(2),ABSOLUTE);
 98             << "P0 = " << GetVertex(0) << G4en <<  99     fFacet2 = G4TriangularFacet (GetVertex(0),GetVertex(2),
 99             << "P1 = " << GetVertex(1) << G4en << 100                                  GetVertex(3),ABSOLUTE);
100             << "P2 = " << GetVertex(2) << G4en << 101 
101             << "P3 = " << GetVertex(3) << G4en << 102     G4TriangularFacet facet3 (GetVertex(0),GetVertex(1),GetVertex(3),ABSOLUTE);
102             << "Side1 length (P0->P1) = " << l << 103     G4TriangularFacet facet4 (GetVertex(1),GetVertex(2),GetVertex(3),ABSOLUTE);
103             << "Side2 length (P1->P2) = " << l << 104 
104             << "Side3 length (P2->P3) = " << l << 105     G4ThreeVector normal12 = fFacet1.GetSurfaceNormal()
105             << "Side4 length (P3->P0) = " << l << 106                            + fFacet2.GetSurfaceNormal();
106             << "Diagonal1 length (P0->P2) = "  << 107     G4ThreeVector normal34 = facet3.GetSurfaceNormal()
107             << "Diagonal2 length (P1->P3) = "  << 108                            + facet4.GetSurfaceNormal();
108     G4Exception("G4QuadrangularFacet::G4Quadra << 109     G4ThreeVector normal = 0.25 * (normal12 + normal34);
109                 "GeomSolids1001", JustWarning, << 110 
110     return;                                    << 111     fFacet1.SetSurfaceNormal (normal);
111   }                                            << 112     fFacet2.SetSurfaceNormal (normal);
112                                                << 113 
113   // Check that vertices are not collinear     << 114     G4ThreeVector vtmp = 0.5 * (e1 + e2);
114   //                                           << 115     fCircumcentre = GetVertex(0) + vtmp;
115   G4double s1 = (e1.cross(e2)).mag()*0.5;      << 116     G4double radiusSqr = vtmp.mag2();
116   G4double s2 = ((e2-e1).cross(e3-e2)).mag()*0 << 117     fRadius = std::sqrt(radiusSqr);
117   G4double s3 = (e2.cross(e3)).mag()*0.5;      << 
118   G4double s4 = (e1.cross(e3)).mag()*0.5;      << 
119                                                << 
120   G4double h1 = 2.*s1 / std::max(std::max(leng << 
121   G4double h2 = 2.*s2 / std::max(std::max(leng << 
122   G4double h3 = 2.*s3 / std::max(std::max(leng << 
123   G4double h4 = 2.*s4 / std::max(std::max(leng << 
124                                                << 
125   if (h1 <= delta || h2 <= delta || h3 <= delt << 
126   {                                            << 
127     ostringstream message;                     << 
128     message << "Facet has three or more collin << 
129             << "P0 = " << GetVertex(0) << G4en << 
130             << "P1 = " << GetVertex(1) << G4en << 
131             << "P2 = " << GetVertex(2) << G4en << 
132             << "P3 = " << GetVertex(3) << G4en << 
133             << "Smallest heights:" << G4endl   << 
134             << "  in triangle P0-P1-P2 = " <<  << 
135             << "  in triangle P1-P2-P3 = " <<  << 
136             << "  in triangle P2-P3-P0 = " <<  << 
137             << "  in triangle P3-P0-P1 = " <<  << 
138     G4Exception("G4QuadrangularFacet::G4Quadra << 
139           "GeomSolids1001", JustWarning, messa << 
140     return;                                    << 
141   }                                            << 
142                                                << 
143   // Check that vertices are coplanar by compu << 
144   // height of tetrahedron comprising of verti << 
145   //                                           << 
146   G4double smax = std::max( std::max(s1,s2), s << 
147   G4double hmin = 0.5 * std::fabs( e1.dot(e2.c << 
148   if (hmin >= epsilon)                         << 
149   {                                            << 
150     ostringstream message;                     << 
151     message << "Facet is not planar." << G4end << 
152             << "Disrepancy = " << hmin << G4en << 
153             << "P0 = " << GetVertex(0) << G4en << 
154             << "P1 = " << GetVertex(1) << G4en << 
155             << "P2 = " << GetVertex(2) << G4en << 
156             << "P3 = " << GetVertex(3);        << 
157     G4Exception("G4QuadrangularFacet::G4Quadra << 
158           "GeomSolids1001", JustWarning, messa << 
159     return;                                    << 
160   }                                               118   }
161                                                << 119   else
162   // Check that facet is convex by computing c << 
163   // of diagonals                              << 
164   //                                           << 
165   G4ThreeVector normal = e2.cross(e3-e1);      << 
166   G4double s = kInfinity, t = kInfinity, magni << 
167   if (magnitude2 > delta*delta) // check: magn << 
168   {                                            << 
169     s = normal.dot(e1.cross(e3-e1)) / magnitud << 
170     t = normal.dot(e1.cross(e2))    / magnitud << 
171   }                                            << 
172   if (s <= 0. || s >= 1. || t <= 0. || t >= 1. << 
173   {                                               120   {
174     ostringstream message;                     << 
175     message << "Facet is not convex." << G4end << 
176             << "Parameters of crosspoint of di << 
177             << s << " and " << t << G4endl     << 
178             << "should both be within (0,1) ra << 
179       << "P0 = " << GetVertex(0) << G4endl     << 
180       << "P1 = " << GetVertex(1) << G4endl     << 
181       << "P2 = " << GetVertex(2) << G4endl     << 
182       << "P3 = " << GetVertex(3);              << 
183     G4Exception("G4QuadrangularFacet::G4Quadra    121     G4Exception("G4QuadrangularFacet::G4QuadrangularFacet()",
184           "GeomSolids1001", JustWarning, messa << 122                 "GeomSolids0002", JustWarning,
185     return;                                    << 123                 "Length of sides of facet are too small or sides not planar.");
                                                   >> 124     G4cout << G4endl;
                                                   >> 125     G4cout << "P0 = " << GetVertex(0) << G4endl;
                                                   >> 126     G4cout << "P1 = " << GetVertex(1) << G4endl;
                                                   >> 127     G4cout << "P2 = " << GetVertex(2) << G4endl;
                                                   >> 128     G4cout << "P3 = " << GetVertex(3) << G4endl;
                                                   >> 129     G4cout << "Side lengths = P0->P1" << length1 << G4endl;    
                                                   >> 130     G4cout << "Side lengths = P1->P2" << length2 << G4endl;    
                                                   >> 131     G4cout << "Side lengths = P2->P3" << length3 << G4endl;    
                                                   >> 132     G4cout << "Side lengths = P3->P0" << length4 << G4endl;    
                                                   >> 133     G4cout << G4endl;
                                                   >> 134     fRadius = 0.0;
186   }                                               135   }
187                                                << 
188   // Define facet                              << 
189   //                                           << 
190   fFacet1 = G4TriangularFacet(GetVertex(0),Get << 
191   fFacet2 = G4TriangularFacet(GetVertex(0),Get << 
192                                                << 
193   normal = normal.unit();                      << 
194   fFacet1.SetSurfaceNormal(normal);            << 
195   fFacet2.SetSurfaceNormal(normal);            << 
196                                                << 
197   G4ThreeVector vtmp = 0.5 * (e1 + e2);        << 
198   fCircumcentre = GetVertex(0) + vtmp;         << 
199   G4double radiusSqr = vtmp.mag2();            << 
200   fRadius = std::sqrt(radiusSqr);              << 
201   // 29.02.2016 Remark by E.Tcherniaev: comput << 
202   // of fCircumcenter and fRadius is wrong, ho << 
203   // it did not create any problem till now.   << 
204   // Bizarre! Need to investigate!             << 
205 }                                                 136 }
206                                                   137 
207 //////////////////////////////////////////////    138 ///////////////////////////////////////////////////////////////////////////////
208 //                                                139 //
209 G4QuadrangularFacet::~G4QuadrangularFacet () = << 140 G4QuadrangularFacet::~G4QuadrangularFacet ()
                                                   >> 141 {
                                                   >> 142 }
210                                                   143 
211 //////////////////////////////////////////////    144 ///////////////////////////////////////////////////////////////////////////////
212 //                                                145 //
213 G4QuadrangularFacet::G4QuadrangularFacet (cons << 146 G4QuadrangularFacet::G4QuadrangularFacet (const G4QuadrangularFacet &rhs)
214   : G4VFacet(rhs)                                 147   : G4VFacet(rhs)
215 {                                                 148 {
216   fFacet1 = rhs.fFacet1;                          149   fFacet1 = rhs.fFacet1;
217   fFacet2 = rhs.fFacet2;                          150   fFacet2 = rhs.fFacet2;
218   fRadius = 0.0;                                  151   fRadius = 0.0;
219 }                                                 152 }
220                                                   153 
221 //////////////////////////////////////////////    154 ///////////////////////////////////////////////////////////////////////////////
222 //                                                155 //
223 G4QuadrangularFacet &                             156 G4QuadrangularFacet &
224 G4QuadrangularFacet::operator=(const G4Quadran << 157 G4QuadrangularFacet::operator=(const G4QuadrangularFacet &rhs)
225 {                                                 158 {
226   if (this == &rhs)  return *this;             << 159   if (this == &rhs)
                                                   >> 160     return *this;
227                                                   161 
228   fFacet1 = rhs.fFacet1;                          162   fFacet1 = rhs.fFacet1;
229   fFacet2 = rhs.fFacet2;                          163   fFacet2 = rhs.fFacet2;
230   fRadius = 0.0;                                  164   fRadius = 0.0;
231                                                   165 
232   return *this;                                   166   return *this;
233 }                                                 167 }
234                                                   168 
235 //////////////////////////////////////////////    169 ///////////////////////////////////////////////////////////////////////////////
236 //                                                170 //
237 G4VFacet* G4QuadrangularFacet::GetClone ()     << 171 G4VFacet *G4QuadrangularFacet::GetClone ()
238 {                                                 172 {
239   auto c = new G4QuadrangularFacet (GetVertex( << 173   G4QuadrangularFacet *c = new G4QuadrangularFacet (GetVertex(0), GetVertex(1),
240                                     GetVertex( << 174                                                     GetVertex(2), GetVertex(3),
                                                   >> 175                                                     ABSOLUTE);
241   return c;                                       176   return c;
242 }                                                 177 }
243                                                   178 
244 //////////////////////////////////////////////    179 ///////////////////////////////////////////////////////////////////////////////
245 //                                                180 //
246 G4ThreeVector G4QuadrangularFacet::Distance (c << 181 G4ThreeVector G4QuadrangularFacet::Distance (const G4ThreeVector &p)
247 {                                                 182 {
248   G4ThreeVector v1 = fFacet1.Distance(p);         183   G4ThreeVector v1 = fFacet1.Distance(p);
249   G4ThreeVector v2 = fFacet2.Distance(p);         184   G4ThreeVector v2 = fFacet2.Distance(p);
250                                                   185 
251   if (v1.mag2() < v2.mag2()) return v1;           186   if (v1.mag2() < v2.mag2()) return v1;
252   else return v2;                                 187   else return v2;
253 }                                                 188 }
254                                                   189 
255 //////////////////////////////////////////////    190 ///////////////////////////////////////////////////////////////////////////////
256 //                                                191 //
257 G4double G4QuadrangularFacet::Distance (const  << 192 G4double G4QuadrangularFacet::Distance (const G4ThreeVector &p,
258                                                   193                                               G4double)
259 {                                                 194 {  
260   G4double dist = Distance(p).mag();              195   G4double dist = Distance(p).mag();
261   return dist;                                    196   return dist;
262 }                                                 197 }
263                                                   198 
264 //////////////////////////////////////////////    199 ///////////////////////////////////////////////////////////////////////////////
265 //                                                200 //
266 G4double G4QuadrangularFacet::Distance (const  << 201 G4double G4QuadrangularFacet::Distance (const G4ThreeVector &p, G4double,
267                                         const     202                                         const G4bool outgoing)
268 {                                                 203 {
269   G4double dist;                                  204   G4double dist;
270                                                   205 
271   G4ThreeVector v = Distance(p);                  206   G4ThreeVector v = Distance(p);
272   G4double dir = v.dot(GetSurfaceNormal());       207   G4double dir = v.dot(GetSurfaceNormal());
273   if ( ((dir > dirTolerance) && (!outgoing))      208   if ( ((dir > dirTolerance) && (!outgoing))
274     || ((dir < -dirTolerance) && outgoing))       209     || ((dir < -dirTolerance) && outgoing))
275     dist = kInfinity;                             210     dist = kInfinity;
276   else                                            211   else 
277     dist = v.mag();                               212     dist = v.mag();
278   return dist;                                    213   return dist;
279 }                                                 214 }
280                                                   215 
281 //////////////////////////////////////////////    216 ///////////////////////////////////////////////////////////////////////////////
282 //                                                217 //
283 G4double G4QuadrangularFacet::Extent (const G4    218 G4double G4QuadrangularFacet::Extent (const G4ThreeVector axis)
284 {                                                 219 {
285   G4double ss  = 0;                               220   G4double ss  = 0;
286                                                   221 
287   for (G4int i = 0; i <= 3; ++i)                  222   for (G4int i = 0; i <= 3; ++i)
288   {                                               223   {
289     G4double sp = GetVertex(i).dot(axis);         224     G4double sp = GetVertex(i).dot(axis);
290     if (sp > ss) ss = sp;                         225     if (sp > ss) ss = sp;
291   }                                               226   }
292   return ss;                                      227   return ss;
293 }                                                 228 }
294                                                   229 
295 //////////////////////////////////////////////    230 ///////////////////////////////////////////////////////////////////////////////
296 //                                                231 //
297 G4bool G4QuadrangularFacet::Intersect (const G << 232 G4bool G4QuadrangularFacet::Intersect (const G4ThreeVector &p,
298                                        const G << 233                                        const G4ThreeVector &v,
299                                              G    234                                              G4bool outgoing,
300                                              G << 235                                              G4double &distance,
301                                              G << 236                                              G4double &distFromSurface,
302                                              G << 237                                              G4ThreeVector &normal)
303 {                                                 238 {
304   G4bool intersect =                              239   G4bool intersect =
305     fFacet1.Intersect(p,v,outgoing,distance,di    240     fFacet1.Intersect(p,v,outgoing,distance,distFromSurface,normal);
306   if (!intersect) intersect =                     241   if (!intersect) intersect =
307     fFacet2.Intersect(p,v,outgoing,distance,di    242     fFacet2.Intersect(p,v,outgoing,distance,distFromSurface,normal);
308   if (!intersect)                                 243   if (!intersect)
309   {                                               244   {
310     distance = distFromSurface = kInfinity;       245     distance = distFromSurface = kInfinity;
311     normal.set(0,0,0);                            246     normal.set(0,0,0);
312   }                                               247   }
313   return intersect;                               248   return intersect;
314 }                                                 249 }
315                                                   250 
316 //////////////////////////////////////////////    251 ///////////////////////////////////////////////////////////////////////////////
317 //                                                252 //
318 // Auxiliary method to get a uniform random po << 253 // Auxiliary method to get a random point on surface
319 //                                                254 //
320 G4ThreeVector G4QuadrangularFacet::GetPointOnF    255 G4ThreeVector G4QuadrangularFacet::GetPointOnFace() const
321 {                                                 256 {
322   G4double s1 = fFacet1.GetArea();             << 257   G4ThreeVector pr = (G4RandFlat::shoot(0.,1.) < 0.5)
323   G4double s2 = fFacet2.GetArea();             << 258                    ? fFacet1.GetPointOnFace() : fFacet2.GetPointOnFace();
324   return ((s1+s2)*G4UniformRand() < s1) ?      << 259   return pr;
325     fFacet1.GetPointOnFace() : fFacet2.GetPoin << 
326 }                                                 260 }
327                                                   261 
328 //////////////////////////////////////////////    262 ///////////////////////////////////////////////////////////////////////////////
329 //                                                263 //
330 // Auxiliary method for returning the surface     264 // Auxiliary method for returning the surface area
331 //                                                265 //
332 G4double G4QuadrangularFacet::GetArea() const  << 266 G4double G4QuadrangularFacet::GetArea()
333 {                                                 267 {
334   G4double area = fFacet1.GetArea() + fFacet2.    268   G4double area = fFacet1.GetArea() + fFacet2.GetArea();
335   return area;                                    269   return area;
336 }                                                 270 }
337                                                   271 
338 //////////////////////////////////////////////    272 ///////////////////////////////////////////////////////////////////////////////
339 //                                                273 //
340 G4String G4QuadrangularFacet::GetEntityType ()    274 G4String G4QuadrangularFacet::GetEntityType () const
341 {                                                 275 {
342   return "G4QuadrangularFacet";                   276   return "G4QuadrangularFacet";
343 }                                                 277 }
344                                                   278 
345 //////////////////////////////////////////////    279 ///////////////////////////////////////////////////////////////////////////////
346 //                                                280 //
347 G4ThreeVector G4QuadrangularFacet::GetSurfaceN    281 G4ThreeVector G4QuadrangularFacet::GetSurfaceNormal () const
348 {                                                 282 {
349   return fFacet1.GetSurfaceNormal();              283   return fFacet1.GetSurfaceNormal();
350 }                                                 284 }
351                                                   285