Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4QuadrangularFacet.cc
이 파일의 문서화 페이지로 가기
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration and of QinetiQ Ltd, *
20 // * subject to DEFCON 705 IPR conditions. *
21 // * By using, copying, modifying or distributing the software (or *
22 // * any work based on the software) you agree to acknowledge its *
23 // * use in resulting scientific publications, and indicate your *
24 // * acceptance of all terms of the Geant4 Software license. *
25 // ********************************************************************
26 //
27 //
28 // $Id: G4QuadrangularFacet.cc 95945 2016-03-03 09:54:38Z gcosmo $
29 //
30 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
31 //
32 // CHANGE HISTORY
33 // --------------
34 //
35 // 31 October 2004 P R Truscott, QinetiQ Ltd, UK - Created.
36 //
37 // 12 October 2012 M Gayer, CERN
38 // New implementation reducing memory requirements by 50%,
39 // and considerable CPU speedup together with the new
40 // implementation of G4TessellatedSolid.
41 //
42 // 29 February 2016 E Tcherniaev, CERN
43 // Added exhaustive tests to catch various problems with a
44 // quadrangular facet: collinear vertices, non planar surface,
45 // degenerate, concave or self intersecting quadrilateral.
46 //
47 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
48 
49 #include "G4QuadrangularFacet.hh"
50 #include "geomdefs.hh"
51 #include "Randomize.hh"
52 
53 using namespace std;
54 
56 //
57 // !!!THIS IS A FUDGE!!! IT'S TWO ADJACENT G4TRIANGULARFACETS
58 // --- NOT EFFICIENT BUT PRACTICAL.
59 //
61  const G4ThreeVector &vt1,
62  const G4ThreeVector &vt2,
63  const G4ThreeVector &vt3,
64  G4FacetVertexType vertexType)
65 {
66  G4double delta = 1.0 * kCarTolerance; // dimension tolerance
67  G4double epsilon = 0.01 * kCarTolerance; // planarity tolerance
68 
69  fRadius = 0.0;
70 
71  G4ThreeVector e1, e2, e3;
72  SetVertex(0, vt0);
73  if (vertexType == ABSOLUTE)
74  {
75  SetVertex(1, vt1);
76  SetVertex(2, vt2);
77  SetVertex(3, vt3);
78 
79  e1 = vt1 - vt0;
80  e2 = vt2 - vt0;
81  e3 = vt3 - vt0;
82  }
83  else
84  {
85  SetVertex(1, vt0 + vt1);
86  SetVertex(2, vt0 + vt2);
87  SetVertex(3, vt0 + vt3);
88 
89  e1 = vt1;
90  e2 = vt2;
91  e3 = vt3;
92  }
93 
94  // Check length of sides and diagonals
95  //
96  G4double leng1 = e1.mag();
97  G4double leng2 = (e2-e1).mag();
98  G4double leng3 = (e3-e2).mag();
99  G4double leng4 = e3.mag();
100 
101  G4double diag1 = e2.mag();
102  G4double diag2 = (e3-e1).mag();
103 
104  if (leng1 <= delta || leng2 <= delta || leng3 <= delta || leng4 <= delta ||
105  diag1 <= delta || diag2 <= delta)
106  {
107  ostringstream message;
108  message << "Sides/diagonals of facet are too small." << G4endl
109  << "P0 = " << GetVertex(0) << G4endl
110  << "P1 = " << GetVertex(1) << G4endl
111  << "P2 = " << GetVertex(2) << G4endl
112  << "P3 = " << GetVertex(3) << G4endl
113  << "Side1 length (P0->P1) = " << leng1 << G4endl
114  << "Side2 length (P1->P2) = " << leng2 << G4endl
115  << "Side3 length (P2->P3) = " << leng3 << G4endl
116  << "Side4 length (P3->P0) = " << leng4 << G4endl
117  << "Diagonal1 length (P0->P2) = " << diag1 << G4endl
118  << "Diagonal2 length (P1->P3) = " << diag2;
119  G4Exception("G4QuadrangularFacet::G4QuadrangularFacet()",
120  "GeomSolids1001", JustWarning, message);
121  return;
122  }
123 
124  // Check that vertices are not collinear
125  //
126  G4double s1 = (e1.cross(e2)).mag()*0.5;
127  G4double s2 = ((e2-e1).cross(e3-e2)).mag()*0.5;
128  G4double s3 = (e2.cross(e3)).mag()*0.5;
129  G4double s4 = (e1.cross(e3)).mag()*0.5;
130 
131  G4double h1 = 2.*s1 / std::max(std::max(leng1,leng2),diag1);
132  G4double h2 = 2.*s2 / std::max(std::max(leng2,leng3),diag2);
133  G4double h3 = 2.*s3 / std::max(std::max(leng3,leng4),diag1);
134  G4double h4 = 2.*s4 / std::max(std::max(leng4,leng1),diag2);
135 
136  if (h1 <= delta || h2 <= delta || h3 <= delta || h4 <= delta )
137  {
138  ostringstream message;
139  message << "Facet has three or more collinear vertices." << G4endl
140  << "P0 = " << GetVertex(0) << G4endl
141  << "P1 = " << GetVertex(1) << G4endl
142  << "P2 = " << GetVertex(2) << G4endl
143  << "P3 = " << GetVertex(3) << G4endl
144  << "Height in P0-P1-P2 = " << h1 << G4endl
145  << "Height in P1-P2-P3 = " << h2 << G4endl
146  << "Height in P2-P3-P4 = " << h3 << G4endl
147  << "Height in P4-P0-P1 = " << h4;
148  G4Exception("G4QuadrangularFacet::G4QuadrangularFacet()",
149  "GeomSolids1001", JustWarning, message);
150  return;
151  }
152 
153  // Check that vertices are coplanar by computing minimal
154  // height of tetrahedron comprising of vertices
155  //
156  G4double smax = std::max( std::max(s1,s2), std::max(s3,s4) );
157  G4double hmin = 0.5 * std::fabs( e1.dot(e2.cross(e3)) ) / smax;
158  if (hmin >= epsilon)
159  {
160  ostringstream message;
161  message << "Facet is not planar." << G4endl
162  << "Disrepancy = " << hmin << G4endl
163  << "P0 = " << GetVertex(0) << G4endl
164  << "P1 = " << GetVertex(1) << G4endl
165  << "P2 = " << GetVertex(2) << G4endl
166  << "P3 = " << GetVertex(3);
167  G4Exception("G4QuadrangularFacet::G4QuadrangularFacet()",
168  "GeomSolids1001", JustWarning, message);
169  return;
170  }
171 
172  // Check that facet is convex by computing crosspoint
173  // of diagonals
174  //
175  G4ThreeVector normal = e2.cross(e3-e1);
176  G4double s = kInfinity, t = kInfinity, magnitude2 = normal.mag2();
177  if (magnitude2 > delta*delta) // check: magnitude2 != 0.
178  {
179  s = normal.dot(e1.cross(e3-e1)) / magnitude2;
180  t = normal.dot(e1.cross(e2)) / magnitude2;
181  }
182  if (s <= 0. || s >= 1. || t <= 0. || t >= 1.)
183  {
184  ostringstream message;
185  message << "Facet is not convex." << G4endl
186  << "Parameters of crosspoint of diagonals: "
187  << s << " and " << t << G4endl
188  << "should both be within (0,1) range" << G4endl
189  << "P0 = " << GetVertex(0) << G4endl
190  << "P1 = " << GetVertex(1) << G4endl
191  << "P2 = " << GetVertex(2) << G4endl
192  << "P3 = " << GetVertex(3);
193  G4Exception("G4QuadrangularFacet::G4QuadrangularFacet()",
194  "GeomSolids1001", JustWarning, message);
195  return;
196  }
197 
198  // Define facet
199  //
200  fFacet1 = G4TriangularFacet(GetVertex(0),GetVertex(1),GetVertex(2),ABSOLUTE);
201  fFacet2 = G4TriangularFacet(GetVertex(0),GetVertex(2),GetVertex(3),ABSOLUTE);
202 
203  normal = normal.unit();
204  fFacet1.SetSurfaceNormal(normal);
205  fFacet2.SetSurfaceNormal(normal);
206 
207  G4ThreeVector vtmp = 0.5 * (e1 + e2);
208  fCircumcentre = GetVertex(0) + vtmp;
209  G4double radiusSqr = vtmp.mag2();
210  fRadius = std::sqrt(radiusSqr);
211  // 29.02.2016 Remark by E.Tcherniaev: computation
212  // of fCircumcenter and fRadius is wrong, however
213  // it did not create any problem till now.
214  // Bizarre! Need to investigate!
215 }
216 
218 //
220 {
221 }
222 
224 //
226  : G4VFacet(rhs)
227 {
228  fFacet1 = rhs.fFacet1;
229  fFacet2 = rhs.fFacet2;
230  fRadius = 0.0;
231 }
232 
234 //
237 {
238  if (this == &rhs)
239  return *this;
240 
241  fFacet1 = rhs.fFacet1;
242  fFacet2 = rhs.fFacet2;
243  fRadius = 0.0;
244 
245  return *this;
246 }
247 
249 //
251 {
253  GetVertex(2), GetVertex(3),
254  ABSOLUTE);
255  return c;
256 }
257 
259 //
261 {
264 
265  if (v1.mag2() < v2.mag2()) return v1;
266  else return v2;
267 }
268 
270 //
272  G4double)
273 {
274  G4double dist = Distance(p).mag();
275  return dist;
276 }
277 
279 //
281  const G4bool outgoing)
282 {
283  G4double dist;
284 
285  G4ThreeVector v = Distance(p);
287  if ( ((dir > dirTolerance) && (!outgoing))
288  || ((dir < -dirTolerance) && outgoing))
289  dist = kInfinity;
290  else
291  dist = v.mag();
292  return dist;
293 }
294 
296 //
298 {
299  G4double ss = 0;
300 
301  for (G4int i = 0; i <= 3; ++i)
302  {
303  G4double sp = GetVertex(i).dot(axis);
304  if (sp > ss) ss = sp;
305  }
306  return ss;
307 }
308 
310 //
312  const G4ThreeVector &v,
313  G4bool outgoing,
314  G4double &distance,
315  G4double &distFromSurface,
317 {
318  G4bool intersect =
319  fFacet1.Intersect(p,v,outgoing,distance,distFromSurface,normal);
320  if (!intersect) intersect =
321  fFacet2.Intersect(p,v,outgoing,distance,distFromSurface,normal);
322  if (!intersect)
323  {
324  distance = distFromSurface = kInfinity;
325  normal.set(0,0,0);
326  }
327  return intersect;
328 }
329 
331 //
332 // Auxiliary method to get a uniform random point on the facet
333 //
335 {
336  G4double s1 = fFacet1.GetArea();
337  G4double s2 = fFacet2.GetArea();
338  return ((s1+s2)*G4UniformRand() < s1) ?
340 }
341 
343 //
344 // Auxiliary method for returning the surface area
345 //
347 {
348  G4double area = fFacet1.GetArea() + fFacet2.GetArea();
349  return area;
350 }
351 
353 //
355 {
356  return "G4QuadrangularFacet";
357 }
358 
360 //
362 {
363  return fFacet1.GetSurfaceNormal();
364 }
void set(double x, double y, double z)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
static const G4double kInfinity
Definition: geomdefs.hh:42
TH1F * h3
G4FacetVertexType
Definition: G4VFacet.hh:56
#define G4endl
Definition: G4ios.hh:61
const char * p
Definition: xmltok.h:285
void message(RunManager *runmanager)
Definition: ts_scorers.cc:72
G4ThreeVector Distance(const G4ThreeVector &p)
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
double dot(const Hep3Vector &) const
const G4double kCarTolerance
static const G4double dirTolerance
Definition: G4VFacet.hh:101
G4GeometryType GetEntityType() const
G4bool Intersect(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool outgoing, G4double &distance, G4double &distFromSurface, G4ThreeVector &normal)
G4ThreeVector GetSurfaceNormal() const
const XML_Char * s
Definition: expat.h:262
G4TriangularFacet fFacet1
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
G4double GetArea() const
G4ThreeVector GetPointOnFace() const
G4QuadrangularFacet(const G4ThreeVector &Pt0, const G4ThreeVector &vt1, const G4ThreeVector &vt2, const G4ThreeVector &vt3, G4FacetVertexType)
G4TriangularFacet fFacet2
G4QuadrangularFacet & operator=(const G4QuadrangularFacet &right)
#define G4UniformRand()
Definition: Randomize.hh:53
double mag2() const
Hep3Vector unit() const
Hep3Vector cross(const Hep3Vector &) const
G4ThreeVector GetVertex(G4int i) const
TH1F * h2
G4bool Intersect(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool outgoing, G4double &distance, G4double &distFromSurface, G4ThreeVector &normal)
double epsilon(double density, double temperature)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
double mag() const
int G4int
Definition: G4Types.hh:78
TDirectory * dir
G4ThreeVector GetPointOnFace() const
const G4int smax
G4ThreeVector Distance(const G4ThreeVector &p)
TH1F * h1
G4double GetArea() const
G4double Extent(const G4ThreeVector axis)
G4ThreeVector GetSurfaceNormal() const