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 9.4.p3)


  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,v 1.9 2010-09-23 10:27:25 gcosmo Exp $
 30 // 12 October 2012 M Gayer, CERN               <<  29 // GEANT4 tag $Name: not supported by cvs2svn $
 31 //                 New implementation reducing <<  30 //
 32 //                 and considerable CPU speedu <<  31 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 33 //                 implementation of G4Tessell <<  32 //
 34 // 29 February 2016 E Tcherniaev, CERN         <<  33 // MODULE:              G4QuadrangularFacet.cc
 35 //                 Added exhaustive tests to c <<  34 //
 36 //                 quadrangular facet: colline <<  35 // Date:                15/06/2005
 37 //                 degenerate, concave or self <<  36 // Author:              P R Truscott
 38 // ------------------------------------------- <<  37 // Organisation:        QinetiQ Ltd, UK
                                                   >>  38 // Customer:            UK Ministry of Defence : RAO CRP TD Electronic Systems
                                                   >>  39 // Contract:            C/MAT/N03517
                                                   >>  40 //
                                                   >>  41 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                                   >>  42 //
                                                   >>  43 // CHANGE HISTORY
                                                   >>  44 // --------------
                                                   >>  45 //
                                                   >>  46 // 31 October 2004, P R Truscott, QinetiQ Ltd, UK - Created.
                                                   >>  47 //
                                                   >>  48 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 39                                                    49 
 40 #include "G4QuadrangularFacet.hh"                  50 #include "G4QuadrangularFacet.hh"
 41 #include "geomdefs.hh"                         <<  51 #include "globals.hh"
 42 #include "Randomize.hh"                            52 #include "Randomize.hh"
 43                                                << 
 44 using namespace std;                           << 
 45                                                    53 
 46 //////////////////////////////////////////////     54 ///////////////////////////////////////////////////////////////////////////////
 47 //                                                 55 //
 48 // Constructing two adjacent G4TriangularFacet <<  56 // !!!THIS IS A FUDGE!!!  IT'S TWO ADJACENT G4TRIANGULARFACETS
 49 // Not efficient, but practical...             <<  57 // --- NOT EFFICIENT BUT PRACTICAL.
 50 //                                                 58 //
 51 G4QuadrangularFacet::G4QuadrangularFacet (cons <<  59 G4QuadrangularFacet::G4QuadrangularFacet (const G4ThreeVector Pt0,
 52                                           cons <<  60                  const G4ThreeVector vt1, const G4ThreeVector vt2,
 53                                           cons <<  61                  const G4ThreeVector vt3, G4FacetVertexType vertexType)
 54                                           cons <<  62   : G4VFacet(), facet1(0), facet2(0)
 55                                                << 
 56 {                                                  63 {
 57   G4double delta   =  1.0 * kCarTolerance; //  <<  64   P0        = Pt0;
 58   G4double epsilon = 0.01 * kCarTolerance; //  <<  65   nVertices = 4;
 59                                                << 
 60   G4ThreeVector e1, e2, e3;                    << 
 61   SetVertex(0, vt0);                           << 
 62   if (vertexType == ABSOLUTE)                      66   if (vertexType == ABSOLUTE)
 63   {                                                67   {
 64     SetVertex(1, vt1);                         <<  68     P.push_back(vt1);
 65     SetVertex(2, vt2);                         <<  69     P.push_back(vt2);
 66     SetVertex(3, vt3);                         <<  70     P.push_back(vt3);
 67                                                <<  71   
 68     e1 = vt1 - vt0;                            <<  72     E.push_back(vt1 - P0);
 69     e2 = vt2 - vt0;                            <<  73     E.push_back(vt2 - P0);
 70     e3 = vt3 - vt0;                            <<  74     E.push_back(vt3 - P0);
 71   }                                                75   }
 72   else                                             76   else
 73   {                                                77   {
 74     SetVertex(1, vt0 + vt1);                   <<  78     P.push_back(P0 + vt1);
 75     SetVertex(2, vt0 + vt2);                   <<  79     P.push_back(P0 + vt2);
 76     SetVertex(3, vt0 + vt3);                   <<  80     P.push_back(P0 + vt3);
 77                                                <<  81   
 78     e1 = vt1;                                  <<  82     E.push_back(vt1);
 79     e2 = vt2;                                  <<  83     E.push_back(vt2);
 80     e3 = vt3;                                  <<  84     E.push_back(vt3);
 81   }                                            << 
 82                                                << 
 83   // Check length of sides and diagonals       << 
 84   //                                           << 
 85   G4double leng1 = e1.mag();                   << 
 86   G4double leng2 = (e2-e1).mag();              << 
 87   G4double leng3 = (e3-e2).mag();              << 
 88   G4double leng4 = e3.mag();                   << 
 89                                                << 
 90   G4double diag1 = e2.mag();                   << 
 91   G4double diag2 = (e3-e1).mag();              << 
 92                                                << 
 93   if (leng1 <= delta || leng2 <= delta || leng << 
 94       diag1 <= delta || diag2 <= delta)        << 
 95   {                                            << 
 96     ostringstream message;                     << 
 97     message << "Sides/diagonals of facet are t << 
 98             << "P0 = " << GetVertex(0) << G4en << 
 99             << "P1 = " << GetVertex(1) << G4en << 
100             << "P2 = " << GetVertex(2) << G4en << 
101             << "P3 = " << GetVertex(3) << G4en << 
102             << "Side1 length (P0->P1) = " << l << 
103             << "Side2 length (P1->P2) = " << l << 
104             << "Side3 length (P2->P3) = " << l << 
105             << "Side4 length (P3->P0) = " << l << 
106             << "Diagonal1 length (P0->P2) = "  << 
107             << "Diagonal2 length (P1->P3) = "  << 
108     G4Exception("G4QuadrangularFacet::G4Quadra << 
109                 "GeomSolids1001", JustWarning, << 
110     return;                                    << 
111   }                                            << 
112                                                << 
113   // Check that vertices are not collinear     << 
114   //                                           << 
115   G4double s1 = (e1.cross(e2)).mag()*0.5;      << 
116   G4double s2 = ((e2-e1).cross(e3-e2)).mag()*0 << 
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   }                                                85   }
142                                                    86 
143   // Check that vertices are coplanar by compu <<  87   G4double length1 = E[0].mag();
144   // height of tetrahedron comprising of verti <<  88   G4double length2 = (P[1]-P[0]).mag();
145   //                                           <<  89   G4double length3 = (P[2]-P[1]).mag();
146   G4double smax = std::max( std::max(s1,s2), s <<  90   G4double length4 = E[2].mag();
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   }                                            << 
161                                                    91   
162   // Check that facet is convex by computing c <<  92   G4ThreeVector normal1 = E[0].cross(E[1]).unit();
163   // of diagonals                              <<  93   G4ThreeVector normal2 = E[1].cross(E[2]).unit(); 
164   //                                           <<  94   
165   G4ThreeVector normal = e2.cross(e3-e1);      <<  95   if (length1 <= kCarTolerance || length2 <= kCarTolerance ||
166   G4double s = kInfinity, t = kInfinity, magni <<  96       length3 <= kCarTolerance || length4 <= kCarTolerance ||
167   if (magnitude2 > delta*delta) // check: magn <<  97       normal1.dot(normal2) < 0.9999999999)
168   {                                                98   {
169     s = normal.dot(e1.cross(e3-e1)) / magnitud <<  99     G4Exception("G4QuadrangularFacet::G4QuadrangularFacet()",
170     t = normal.dot(e1.cross(e2))    / magnitud << 100                 "InvalidSetup", JustWarning,
                                                   >> 101                 "Length of sides of facet are too small or sides not planar.");
                                                   >> 102     G4cerr << G4endl;
                                                   >> 103     G4cerr << "P0 = " << P0   << G4endl;
                                                   >> 104     G4cerr << "P1 = " << P[0] << G4endl;
                                                   >> 105     G4cerr << "P2 = " << P[1] << G4endl;
                                                   >> 106     G4cerr << "P3 = " << P[2] << G4endl;
                                                   >> 107     G4cerr << "Side lengths = P0->P1" << length1 << G4endl;    
                                                   >> 108     G4cerr << "Side lengths = P1->P2" << length2 << G4endl;    
                                                   >> 109     G4cerr << "Side lengths = P2->P3" << length3 << G4endl;    
                                                   >> 110     G4cerr << "Side lengths = P3->P0" << length4 << G4endl;    
                                                   >> 111     G4cerr << G4endl;
                                                   >> 112     
                                                   >> 113     isDefined     = false;
                                                   >> 114     geometryType  = "G4QuadragularFacet";
                                                   >> 115     surfaceNormal = G4ThreeVector(0.0,0.0,0.0);
171   }                                               116   }
172   if (s <= 0. || s >= 1. || t <= 0. || t >= 1. << 117   else
173   {                                               118   {
174     ostringstream message;                     << 119     isDefined     = true;
175     message << "Facet is not convex." << G4end << 120     geometryType  = "G4QuadrangularFacet";
176             << "Parameters of crosspoint of di << 121     
177             << s << " and " << t << G4endl     << 122     facet1 = new G4TriangularFacet(P0,P[0],P[1],ABSOLUTE);
178             << "should both be within (0,1) ra << 123     facet2 = new G4TriangularFacet(P0,P[1],P[2],ABSOLUTE);
179       << "P0 = " << GetVertex(0) << G4endl     << 124     surfaceNormal = normal1;
180       << "P1 = " << GetVertex(1) << G4endl     << 125     
181       << "P2 = " << GetVertex(2) << G4endl     << 126     G4ThreeVector vtmp = 0.5 * (E[0] + E[1]);
182       << "P3 = " << GetVertex(3);              << 127     circumcentre       = P0 + vtmp;
183     G4Exception("G4QuadrangularFacet::G4Quadra << 128     radiusSqr          = vtmp.mag2();
184           "GeomSolids1001", JustWarning, messa << 129     radius             = std::sqrt(radiusSqr);
185     return;                                    << 130   
                                                   >> 131     for (size_t i=0; i<4; i++) I.push_back(0);
186   }                                               132   }
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 }                                                 133 }
206                                                   134 
207 //////////////////////////////////////////////    135 ///////////////////////////////////////////////////////////////////////////////
208 //                                                136 //
209 G4QuadrangularFacet::~G4QuadrangularFacet () = << 137 G4QuadrangularFacet::~G4QuadrangularFacet ()
                                                   >> 138 {
                                                   >> 139   delete facet1;
                                                   >> 140   delete facet2;
                                                   >> 141   
                                                   >> 142   P.clear();
                                                   >> 143   E.clear();
                                                   >> 144   I.clear();
                                                   >> 145 }
210                                                   146 
211 //////////////////////////////////////////////    147 ///////////////////////////////////////////////////////////////////////////////
212 //                                                148 //
213 G4QuadrangularFacet::G4QuadrangularFacet (cons << 149 G4QuadrangularFacet::G4QuadrangularFacet (const G4QuadrangularFacet &rhs)
214   : G4VFacet(rhs)                                 150   : G4VFacet(rhs)
215 {                                                 151 {
216   fFacet1 = rhs.fFacet1;                       << 152   facet1 = new G4TriangularFacet(*(rhs.facet1));
217   fFacet2 = rhs.fFacet2;                       << 153   facet2 = new G4TriangularFacet(*(rhs.facet2));
218   fRadius = 0.0;                               << 
219 }                                                 154 }
220                                                   155 
221 //////////////////////////////////////////////    156 ///////////////////////////////////////////////////////////////////////////////
222 //                                                157 //
223 G4QuadrangularFacet &                          << 158 const G4QuadrangularFacet &
224 G4QuadrangularFacet::operator=(const G4Quadran << 159 G4QuadrangularFacet::operator=(G4QuadrangularFacet &rhs)
225 {                                                 160 {
226   if (this == &rhs)  return *this;             << 161    // Check assignment to self
                                                   >> 162    //
                                                   >> 163    if (this == &rhs)  { return *this; }
227                                                   164 
228   fFacet1 = rhs.fFacet1;                       << 165    // Copy base class data
229   fFacet2 = rhs.fFacet2;                       << 166    //
230   fRadius = 0.0;                               << 167    G4VFacet::operator=(rhs);
231                                                   168 
232   return *this;                                << 169    // Copy data
                                                   >> 170    //
                                                   >> 171    delete facet1; facet1 = new G4TriangularFacet(*(rhs.facet1));
                                                   >> 172    delete facet2; facet2 = new G4TriangularFacet(*(rhs.facet2));
                                                   >> 173 
                                                   >> 174    return *this;
233 }                                                 175 }
234                                                   176 
235 //////////////////////////////////////////////    177 ///////////////////////////////////////////////////////////////////////////////
236 //                                                178 //
237 G4VFacet* G4QuadrangularFacet::GetClone ()     << 179 G4VFacet *G4QuadrangularFacet::GetClone ()
238 {                                                 180 {
239   auto c = new G4QuadrangularFacet (GetVertex( << 181   G4QuadrangularFacet *c =
240                                     GetVertex( << 182     new G4QuadrangularFacet (P0, P[0], P[1], P[2], ABSOLUTE);
241   return c;                                    << 183   G4VFacet *cc         = 0;
                                                   >> 184   cc                   = c;
                                                   >> 185   return cc;
242 }                                                 186 }
243                                                   187 
244 //////////////////////////////////////////////    188 ///////////////////////////////////////////////////////////////////////////////
245 //                                                189 //
246 G4ThreeVector G4QuadrangularFacet::Distance (c << 190 G4ThreeVector G4QuadrangularFacet::Distance (const G4ThreeVector &p)
247 {                                                 191 {
248   G4ThreeVector v1 = fFacet1.Distance(p);      << 192   G4ThreeVector v1 = facet1->Distance(p);
249   G4ThreeVector v2 = fFacet2.Distance(p);      << 193   G4ThreeVector v2 = facet2->Distance(p);
250                                                << 194   
251   if (v1.mag2() < v2.mag2()) return v1;           195   if (v1.mag2() < v2.mag2()) return v1;
252   else return v2;                                 196   else return v2;
253 }                                                 197 }
254                                                   198 
255 //////////////////////////////////////////////    199 ///////////////////////////////////////////////////////////////////////////////
256 //                                                200 //
257 G4double G4QuadrangularFacet::Distance (const  << 201 G4double G4QuadrangularFacet::Distance (const G4ThreeVector &p,
258                                                << 202   const G4double)
259 {                                              << 203 {
260   G4double dist = Distance(p).mag();           << 204   /*G4ThreeVector D  = P0 - p;
                                                   >> 205   G4double d       = E[0].dot(D);
                                                   >> 206   G4double e       = E[1].dot(D);
                                                   >> 207   G4double s       = b*e - c*d;
                                                   >> 208   G4double t       = b*d - a*e;*/
                                                   >> 209   G4double dist = kInfinity;
                                                   >> 210   
                                                   >> 211   /*if (s+t > 1.0 || s < 0.0 || t < 0.0)
                                                   >> 212   {
                                                   >> 213     G4ThreeVector D0 = P0 - p;
                                                   >> 214     G4ThreeVector D1 = P[0] - p;
                                                   >> 215     G4ThreeVector D2 = P[1] - p;
                                                   >> 216     
                                                   >> 217     G4double d0 = D0.mag();
                                                   >> 218     G4double d1 = D1.mag();
                                                   >> 219     G4double d2 = D2.mag();
                                                   >> 220     
                                                   >> 221     dist = min(d0, min(d1, d2));
                                                   >> 222     if (dist > minDist) return kInfinity;
                                                   >> 223   }*/
                                                   >> 224   
                                                   >> 225   dist = Distance(p).mag();
                                                   >> 226   
261   return dist;                                    227   return dist;
262 }                                                 228 }
263                                                   229 
264 //////////////////////////////////////////////    230 ///////////////////////////////////////////////////////////////////////////////
265 //                                                231 //
266 G4double G4QuadrangularFacet::Distance (const  << 232 G4double G4QuadrangularFacet::Distance (const G4ThreeVector &p,
267                                         const  << 233                                         const G4double, const G4bool outgoing)
268 {                                                 234 {
269   G4double dist;                               << 235   /*G4ThreeVector D  = P0 - p;
270                                                << 236   G4double d       = E[0].dot(D);
                                                   >> 237   G4double e       = E[1].dot(D);
                                                   >> 238   G4double s       = b*e - c*d;
                                                   >> 239   G4double t       = b*d - a*e;*/
                                                   >> 240   G4double dist = kInfinity;
                                                   >> 241   
                                                   >> 242   /*if (s+t > 1.0 || s < 0.0 || t < 0.0)
                                                   >> 243   {
                                                   >> 244     G4ThreeVector D0 = P0 - p;
                                                   >> 245     G4ThreeVector D1 = P[0] - p;
                                                   >> 246     G4ThreeVector D2 = P[1] - p;
                                                   >> 247     
                                                   >> 248     G4double d0 = D0.mag();
                                                   >> 249     G4double d1 = D1.mag();
                                                   >> 250     G4double d2 = D2.mag();
                                                   >> 251     
                                                   >> 252     dist = min(d0, min(d1, d2));
                                                   >> 253     if (dist > minDist ||
                                                   >> 254       (D0.dot(surfaceNormal) > 0.0 && !outgoing) ||
                                                   >> 255       (D0.dot(surfaceNormal) < 0.0 && outgoing)) return kInfinity;
                                                   >> 256   }*/
                                                   >> 257   
271   G4ThreeVector v = Distance(p);                  258   G4ThreeVector v = Distance(p);
272   G4double dir = v.dot(GetSurfaceNormal());    << 259   G4double dir    = v.dot(surfaceNormal);
273   if ( ((dir > dirTolerance) && (!outgoing))   << 260   if ((dir > dirTolerance && !outgoing) ||
274     || ((dir < -dirTolerance) && outgoing))    << 261       (dir <-dirTolerance && outgoing)) dist = kInfinity;
275     dist = kInfinity;                          << 262   else dist = v.mag();
276   else                                         << 263   
277     dist = v.mag();                            << 
278   return dist;                                    264   return dist;
279 }                                                 265 }
280                                                   266 
281 //////////////////////////////////////////////    267 ///////////////////////////////////////////////////////////////////////////////
282 //                                                268 //
283 G4double G4QuadrangularFacet::Extent (const G4    269 G4double G4QuadrangularFacet::Extent (const G4ThreeVector axis)
284 {                                                 270 {
285   G4double ss  = 0;                            << 271   G4double s  = P0.dot(axis);
286                                                << 272   for (G4ThreeVectorList::iterator it=P.begin(); it!=P.end(); it++)
287   for (G4int i = 0; i <= 3; ++i)               << 
288   {                                               273   {
289     G4double sp = GetVertex(i).dot(axis);      << 274     G4double sp = it->dot(axis);
290     if (sp > ss) ss = sp;                      << 275     if (sp > s) s = sp;
291   }                                               276   }
292   return ss;                                   << 277 
                                                   >> 278   return s;
293 }                                                 279 }
294                                                   280 
295 //////////////////////////////////////////////    281 ///////////////////////////////////////////////////////////////////////////////
296 //                                                282 //
297 G4bool G4QuadrangularFacet::Intersect (const G << 283 G4bool G4QuadrangularFacet::Intersect (const G4ThreeVector &p,
298                                        const G << 284   const G4ThreeVector &v, G4bool outgoing, G4double &distance,
299                                              G << 285   G4double &distFromSurface, G4ThreeVector &normal)
300                                              G << 
301                                              G << 
302                                              G << 
303 {                                                 286 {
304   G4bool intersect =                              287   G4bool intersect =
305     fFacet1.Intersect(p,v,outgoing,distance,di << 288     facet1->Intersect(p,v,outgoing,distance,distFromSurface,normal);
306   if (!intersect) intersect =                  << 
307     fFacet2.Intersect(p,v,outgoing,distance,di << 
308   if (!intersect)                                 289   if (!intersect)
309   {                                               290   {
310     distance = distFromSurface = kInfinity;    << 291     intersect = facet2->Intersect(p,v,outgoing,distance,distFromSurface,normal);
311     normal.set(0,0,0);                         << 
312   }                                               292   }
                                                   >> 293   
                                                   >> 294   if (!intersect)
                                                   >> 295   {
                                                   >> 296     distance        = kInfinity;
                                                   >> 297     distFromSurface = kInfinity;
                                                   >> 298     normal          = G4ThreeVector(0.0,0.0,0.0);
                                                   >> 299   }
                                                   >> 300   
313   return intersect;                               301   return intersect;
314 }                                                 302 }
315                                                   303 
316 ////////////////////////////////////////////// << 304 ////////////////////////////////////////////////////////////////////////
317 //                                                305 //
318 // Auxiliary method to get a uniform random po << 306 // GetPointOnFace
319 //                                                307 //
                                                   >> 308 // Auxiliary method for get a random point on surface
                                                   >> 309 
320 G4ThreeVector G4QuadrangularFacet::GetPointOnF    310 G4ThreeVector G4QuadrangularFacet::GetPointOnFace() const
321 {                                                 311 {
322   G4double s1 = fFacet1.GetArea();             << 312   G4ThreeVector pr;
323   G4double s2 = fFacet2.GetArea();             << 313 
324   return ((s1+s2)*G4UniformRand() < s1) ?      << 314   if ( G4UniformRand() < 0.5 )
325     fFacet1.GetPointOnFace() : fFacet2.GetPoin << 315   {
                                                   >> 316     pr = facet1->GetPointOnFace();
                                                   >> 317   }
                                                   >> 318   else
                                                   >> 319   {
                                                   >> 320     pr = facet2->GetPointOnFace();
                                                   >> 321   }
                                                   >> 322 
                                                   >> 323   return pr;
326 }                                                 324 }
327                                                   325 
328 ////////////////////////////////////////////// << 326 ////////////////////////////////////////////////////////////////////////
329 //                                                327 //
330 // Auxiliary method for returning the surface  << 328 // GetArea
331 //                                                329 //
332 G4double G4QuadrangularFacet::GetArea() const  << 330 // Auxiliary method for returning the surface area
333 {                                              << 
334   G4double area = fFacet1.GetArea() + fFacet2. << 
335   return area;                                 << 
336 }                                              << 
337                                                   331 
338 ////////////////////////////////////////////// << 332 G4double G4QuadrangularFacet::GetArea()
339 //                                             << 
340 G4String G4QuadrangularFacet::GetEntityType () << 
341 {                                                 333 {
342   return "G4QuadrangularFacet";                << 334   if (!area)  { area = facet1->GetArea() + facet2->GetArea(); }
343 }                                              << 
344                                                   335 
345 ////////////////////////////////////////////// << 336   return area;
346 //                                             << 
347 G4ThreeVector G4QuadrangularFacet::GetSurfaceN << 
348 {                                              << 
349   return fFacet1.GetSurfaceNormal();           << 
350 }                                                 337 }
351                                                   338