Geant4  10.00.p03
G4ExtrudedSolid.cc
Go to the documentation of this file.
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. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id: G4ExtrudedSolid.cc 69584 2013-05-08 14:12:37Z gcosmo $
28 //
29 //
30 // --------------------------------------------------------------------
31 // GEANT 4 class source file
32 //
33 // G4ExtrudedSolid.cc
34 //
35 // Author: Ivana Hrivnacova, IPN Orsay
36 // --------------------------------------------------------------------
37 
38 #include <set>
39 #include <algorithm>
40 #include <cmath>
41 #include <iomanip>
42 
43 #include "G4ExtrudedSolid.hh"
44 #include "G4PhysicalConstants.hh"
45 #include "G4SystemOfUnits.hh"
46 #include "G4VFacet.hh"
47 #include "G4TriangularFacet.hh"
48 #include "G4QuadrangularFacet.hh"
49 
50 //_____________________________________________________________________________
51 
53  std::vector<G4TwoVector> polygon,
54  std::vector<ZSection> zsections)
55  : G4TessellatedSolid(pName),
56  fNv(polygon.size()),
57  fNz(zsections.size()),
58  fPolygon(),
59  fZSections(),
60  fTriangles(),
61  fIsConvex(false),
62  fGeometryType("G4ExtrudedSolid")
63 
64 {
65  // General constructor
66 
67  // First check input parameters
68 
69  if ( fNv < 3 )
70  {
71  std::ostringstream message;
72  message << "Number of polygon vertices < 3 - " << pName;
73  G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids0002",
74  FatalErrorInArgument, message);
75  }
76 
77  if ( fNz < 2 )
78  {
79  std::ostringstream message;
80  message << "Number of z-sides < 2 - " << pName;
81  G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids0002",
82  FatalErrorInArgument, message);
83  }
84 
85  for ( G4int i=0; i<fNz-1; ++i )
86  {
87  if ( zsections[i].fZ > zsections[i+1].fZ )
88  {
89  std::ostringstream message;
90  message << "Z-sections have to be ordered by z value (z0 < z1 < z2...) - "
91  << pName;
92  G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids0002",
93  FatalErrorInArgument, message);
94  }
95  if ( std::fabs( zsections[i+1].fZ - zsections[i].fZ ) < kCarTolerance * 0.5 )
96  {
97  std::ostringstream message;
98  message << "Z-sections with the same z position are not supported - "
99  << pName;
100  G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids0001",
101  FatalException, message);
102  }
103  }
104 
105  // Check if polygon vertices are defined clockwise
106  // (the area is positive if polygon vertices are defined anti-clockwise)
107  //
108  G4double area = 0.;
109  for ( G4int i=0; i<fNv; ++i ) {
110  G4int j = i+1;
111  if ( j == fNv ) j = 0;
112  area += 0.5 * ( polygon[i].x()*polygon[j].y() - polygon[j].x()*polygon[i].y());
113  }
114 
115  // Copy polygon
116  //
117  if ( area < 0. ) {
118  // Polygon vertices are defined clockwise, we just copy the polygon
119  for ( G4int i=0; i<fNv; ++i ) { fPolygon.push_back(polygon[i]); }
120  }
121  else {
122  // Polygon vertices are defined anti-clockwise, we revert them
123  //G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids1001",
124  // JustWarning,
125  // "Polygon vertices defined anti-clockwise, reverting polygon");
126  for ( G4int i=0; i<fNv; ++i ) { fPolygon.push_back(polygon[fNv-i-1]); }
127  }
128 
129 
130  // Copy z-sections
131  //
132  for ( G4int i=0; i<fNz; ++i ) { fZSections.push_back(zsections[i]); }
133 
134 
135  G4bool result = MakeFacets();
136  if (!result)
137  {
138  std::ostringstream message;
139  message << "Making facets failed - " << pName;
140  G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids0003",
141  FatalException, message);
142  }
143  fIsConvex = IsConvex();
144 
145 
147 }
148 
149 //_____________________________________________________________________________
150 
152  std::vector<G4TwoVector> polygon,
153  G4double dz,
154  G4TwoVector off1, G4double scale1,
155  G4TwoVector off2, G4double scale2 )
156  : G4TessellatedSolid(pName),
157  fNv(polygon.size()),
158  fNz(2),
159  fPolygon(),
160  fZSections(),
161  fTriangles(),
162  fIsConvex(false),
163  fGeometryType("G4ExtrudedSolid")
164 
165 {
166  // Special constructor for solid with 2 z-sections
167 
168  // First check input parameters
169  //
170  if ( fNv < 3 )
171  {
172  std::ostringstream message;
173  message << "Number of polygon vertices < 3 - " << pName;
174  G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids0002",
175  FatalErrorInArgument, message);
176  }
177 
178  // Check if polygon vertices are defined clockwise
179  // (the area is positive if polygon vertices are defined anti-clockwise)
180 
181  G4double area = 0.;
182  for ( G4int i=0; i<fNv; ++i )
183  {
184  G4int j = i+1;
185  if ( j == fNv ) { j = 0; }
186  area += 0.5 * ( polygon[i].x()*polygon[j].y()
187  - polygon[j].x()*polygon[i].y());
188  }
189 
190  // Copy polygon
191  //
192  if ( area < 0. )
193  {
194  // Polygon vertices are defined clockwise, we just copy the polygon
195  for ( G4int i=0; i<fNv; ++i ) { fPolygon.push_back(polygon[i]); }
196  }
197  else
198  {
199  // Polygon vertices are defined anti-clockwise, we revert them
200  //G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids1001",
201  // JustWarning,
202  // "Polygon vertices defined anti-clockwise, reverting polygon");
203  for ( G4int i=0; i<fNv; ++i ) { fPolygon.push_back(polygon[fNv-i-1]); }
204  }
205 
206  // Copy z-sections
207  //
208  fZSections.push_back(ZSection(-dz, off1, scale1));
209  fZSections.push_back(ZSection( dz, off2, scale2));
210 
211  G4bool result = MakeFacets();
212  if (!result)
213  {
214  std::ostringstream message;
215  message << "Making facets failed - " << pName;
216  G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids0003",
217  FatalException, message);
218  }
219  fIsConvex = IsConvex();
220 
222 }
223 
224 //_____________________________________________________________________________
225 
227  : G4TessellatedSolid(a), fNv(0), fNz(0), fPolygon(), fZSections(),
228  fTriangles(), fIsConvex(false), fGeometryType("G4ExtrudedSolid")
229 {
230  // Fake default constructor - sets only member data and allocates memory
231  // for usage restricted to object persistency.
232 }
233 
234 //_____________________________________________________________________________
235 
237  : G4TessellatedSolid(rhs), fNv(rhs.fNv), fNz(rhs.fNz),
238  fPolygon(rhs.fPolygon), fZSections(rhs.fZSections),
239  fTriangles(rhs.fTriangles), fIsConvex(rhs.fIsConvex),
240  fGeometryType(rhs.fGeometryType), fKScales(rhs.fKScales),
241  fScale0s(rhs.fScale0s), fKOffsets(rhs.fKOffsets), fOffset0s(rhs.fOffset0s)
242 {
243 }
244 
245 
246 //_____________________________________________________________________________
247 
249 {
250  // Check assignment to self
251  //
252  if (this == &rhs) { return *this; }
253 
254  // Copy base class data
255  //
257 
258  // Copy data
259  //
260  fNv = rhs.fNv; fNz = rhs.fNz;
261  fPolygon = rhs.fPolygon; fZSections = rhs.fZSections;
264  fScale0s = rhs.fScale0s; fKOffsets = rhs.fKOffsets;
265  fOffset0s = rhs.fOffset0s;
266 
267  return *this;
268 }
269 
270 //_____________________________________________________________________________
271 
273 {
274  // Destructor
275 }
276 
277 //_____________________________________________________________________________
278 
280 {
281  // Compute parameters for point projections p(z)
282  // to the polygon scale & offset:
283  // scale(z) = k*z + scale0
284  // offset(z) = l*z + offset0
285  // p(z) = scale(z)*p0 + offset(z)
286  // p0 = (p(z) - offset(z))/scale(z);
287  //
288 
289  for ( G4int iz=0; iz<fNz-1; ++iz)
290  {
291  G4double z1 = fZSections[iz].fZ;
292  G4double z2 = fZSections[iz+1].fZ;
293  G4double scale1 = fZSections[iz].fScale;
294  G4double scale2 = fZSections[iz+1].fScale;
295  G4TwoVector off1 = fZSections[iz].fOffset;
296  G4TwoVector off2 = fZSections[iz+1].fOffset;
297 
298  G4double kscale = (scale2 - scale1)/(z2 - z1);
299  G4double scale0 = scale2 - kscale*(z2 - z1)/2.0;
300  G4TwoVector koff = (off2 - off1)/(z2 - z1);
301  G4TwoVector off0 = off2 - koff*(z2 - z1)/2.0;
302 
303  fKScales.push_back(kscale);
304  fScale0s.push_back(scale0);
305  fKOffsets.push_back(koff);
306  fOffset0s.push_back(off0);
307  }
308 }
309 
310 
311 //_____________________________________________________________________________
312 
314 {
315  // Shift and scale vertices
316 
317  return G4ThreeVector( fPolygon[ind].x() * fZSections[iz].fScale
318  + fZSections[iz].fOffset.x(),
319  fPolygon[ind].y() * fZSections[iz].fScale
320  + fZSections[iz].fOffset.y(), fZSections[iz].fZ);
321 }
322 
323 //_____________________________________________________________________________
324 
325 
327 {
328  // Project point in the polygon scale
329  // scale(z) = k*z + scale0
330  // offset(z) = l*z + offset0
331  // p(z) = scale(z)*p0 + offset(z)
332  // p0 = (p(z) - offset(z))/scale(z);
333 
334  // Select projection (z-segment of the solid) according to p.z()
335  //
336  G4int iz = 0;
337  while ( point.z() > fZSections[iz+1].fZ && iz < fNz-2 ) { ++iz; }
338 
339  G4double z0 = ( fZSections[iz+1].fZ + fZSections[iz].fZ )/2.0;
340  G4TwoVector p2(point.x(), point.y());
341  G4double pscale = fKScales[iz]*(point.z()-z0) + fScale0s[iz];
342  G4TwoVector poffset = fKOffsets[iz]*(point.z()-z0) + fOffset0s[iz];
343 
344  // G4cout << point << " projected to "
345  // << iz << "-th z-segment polygon as "
346  // << (p2 - poffset)/pscale << G4endl;
347 
348  // pscale is always >0 as it is an interpolation between two
349  // positive scale values
350  //
351  return (p2 - poffset)/pscale;
352 }
353 
354 //_____________________________________________________________________________
355 
357  G4TwoVector l1, G4TwoVector l2) const
358 {
359  // Return true if p is on the line through l1, l2
360 
361  if ( l1.x() == l2.x() )
362  {
363  return std::fabs(p.x() - l1.x()) < kCarTolerance * 0.5;
364  }
365  G4double slope= ((l2.y() - l1.y())/(l2.x() - l1.x()));
366  G4double predy= l1.y() + slope *(p.x() - l1.x());
367  G4double dy= p.y() - predy;
368 
369  // Calculate perpendicular distance
370  //
371  // G4double perpD= std::fabs(dy) / std::sqrt( 1 + slope * slope );
372  // G4bool simpleComp= (perpD<0.5*kCarTolerance);
373 
374  // Check perpendicular distance vs tolerance 'directly'
375  //
376  const G4double tol= 0.5 * kCarTolerance ;
377  G4bool squareComp= (dy*dy < (1+slope*slope) * tol * tol);
378 
379  // return simpleComp;
380  return squareComp;
381 }
382 
383 //_____________________________________________________________________________
384 
386  G4TwoVector l1, G4TwoVector l2) const
387 {
388  // Return true if p is on the line through l1, l2 and lies between
389  // l1 and l2
390 
391  if ( p.x() < std::min(l1.x(), l2.x()) - kCarTolerance * 0.5 ||
392  p.x() > std::max(l1.x(), l2.x()) + kCarTolerance * 0.5 ||
393  p.y() < std::min(l1.y(), l2.y()) - kCarTolerance * 0.5 ||
394  p.y() > std::max(l1.y(), l2.y()) + kCarTolerance * 0.5 )
395  {
396  return false;
397  }
398 
399  return IsSameLine(p, l1, l2);
400 }
401 
402 //_____________________________________________________________________________
403 
405  G4TwoVector l1, G4TwoVector l2) const
406 {
407  // Return true if p1 and p2 are on the same side of the line through l1, l2
408 
409  return ( (p1.x() - l1.x()) * (l2.y() - l1.y())
410  - (l2.x() - l1.x()) * (p1.y() - l1.y()) )
411  * ( (p2.x() - l1.x()) * (l2.y() - l1.y())
412  - (l2.x() - l1.x()) * (p2.y() - l1.y()) ) > 0;
413 }
414 
415 //_____________________________________________________________________________
416 
418  G4TwoVector c, G4TwoVector p) const
419 {
420  // Return true if p is inside of triangle abc or on its edges,
421  // else returns false
422 
423  // Check extent first
424  //
425  if ( ( p.x() < a.x() && p.x() < b.x() && p.x() < c.x() ) ||
426  ( p.x() > a.x() && p.x() > b.x() && p.x() > c.x() ) ||
427  ( p.y() < a.y() && p.y() < b.y() && p.y() < c.y() ) ||
428  ( p.y() > a.y() && p.y() > b.y() && p.y() > c.y() ) ) return false;
429 
430  G4bool inside
431  = IsSameSide(p, a, b, c)
432  && IsSameSide(p, b, a, c)
433  && IsSameSide(p, c, a, b);
434 
435  G4bool onEdge
436  = IsSameLineSegment(p, a, b)
437  || IsSameLineSegment(p, b, c)
438  || IsSameLineSegment(p, c, a);
439 
440  return inside || onEdge;
441 }
442 
443 //_____________________________________________________________________________
444 
445 G4double
447 {
448  // Return the angle of the vertex in po
449 
450  G4TwoVector t1 = pa - po;
451  G4TwoVector t2 = pb - po;
452 
453  G4double result = (std::atan2(t1.y(), t1.x()) - std::atan2(t2.y(), t2.x()));
454 
455  if ( result < 0 ) result += 2*pi;
456 
457  return result;
458 }
459 
460 //_____________________________________________________________________________
461 
462 G4VFacet*
464 {
465  // Create a triangular facet from the polygon points given by indices
466  // forming the down side ( the normal goes in -z)
467 
468  std::vector<G4ThreeVector> vertices;
469  vertices.push_back(GetVertex(0, ind1));
470  vertices.push_back(GetVertex(0, ind2));
471  vertices.push_back(GetVertex(0, ind3));
472 
473  // first vertex most left
474  //
475  G4ThreeVector cross
476  = (vertices[1]-vertices[0]).cross(vertices[2]-vertices[1]);
477 
478  if ( cross.z() > 0.0 )
479  {
480  // vertices ardered clock wise has to be reordered
481 
482  // G4cout << "G4ExtrudedSolid::MakeDownFacet: reordering vertices "
483  // << ind1 << ", " << ind2 << ", " << ind3 << G4endl;
484 
485  G4ThreeVector tmp = vertices[1];
486  vertices[1] = vertices[2];
487  vertices[2] = tmp;
488  }
489 
490  return new G4TriangularFacet(vertices[0], vertices[1],
491  vertices[2], ABSOLUTE);
492 }
493 
494 //_____________________________________________________________________________
495 
496 G4VFacet*
498 {
499  // Creates a triangular facet from the polygon points given by indices
500  // forming the upper side ( z>0 )
501 
502  std::vector<G4ThreeVector> vertices;
503  vertices.push_back(GetVertex(fNz-1, ind1));
504  vertices.push_back(GetVertex(fNz-1, ind2));
505  vertices.push_back(GetVertex(fNz-1, ind3));
506 
507  // first vertex most left
508  //
509  G4ThreeVector cross
510  = (vertices[1]-vertices[0]).cross(vertices[2]-vertices[1]);
511 
512  if ( cross.z() < 0.0 )
513  {
514  // vertices ordered clock wise has to be reordered
515 
516  // G4cout << "G4ExtrudedSolid::MakeUpFacet: reordering vertices "
517  // << ind1 << ", " << ind2 << ", " << ind3 << G4endl;
518 
519  G4ThreeVector tmp = vertices[1];
520  vertices[1] = vertices[2];
521  vertices[2] = tmp;
522  }
523 
524  return new G4TriangularFacet(vertices[0], vertices[1],
525  vertices[2], ABSOLUTE);
526 }
527 
528 //_____________________________________________________________________________
529 
531 {
532  // Decompose polygonal sides in triangular facets
533 
534  typedef std::pair < G4TwoVector, G4int > Vertex;
535 
536  // Fill one more vector
537  //
538  std::vector< Vertex > verticesToBeDone;
539  for ( G4int i=0; i<fNv; ++i )
540  {
541  verticesToBeDone.push_back(Vertex(fPolygon[i], i));
542  }
543  std::vector< Vertex > ears;
544 
545  std::vector< Vertex >::iterator c1 = verticesToBeDone.begin();
546  std::vector< Vertex >::iterator c2 = c1+1;
547  std::vector< Vertex >::iterator c3 = c1+2;
548  while ( verticesToBeDone.size()>2 )
549  {
550 
551  // G4cout << "Looking at triangle : "
552  // << c1->second << " " << c2->second
553  // << " " << c3->second << G4endl;
554 
555  // skip concave vertices
556  //
557  G4double angle = GetAngle(c2->first, c3->first, c1->first);
558  //G4cout << "angle " << angle << G4endl;
559 
560  G4int counter = 0;
561  while ( angle > pi )
562  {
563  // G4cout << "Skipping concave vertex " << c2->second << G4endl;
564 
565  // try next three consecutive vertices
566  //
567  c1 = c2;
568  c2 = c3;
569  ++c3;
570  if ( c3 == verticesToBeDone.end() ) { c3 = verticesToBeDone.begin(); }
571 
572  // G4cout << "Looking at triangle : "
573  // << c1->second << " " << c2->second
574  // << " " << c3->second << G4endl;
575 
576  angle = GetAngle(c2->first, c3->first, c1->first);
577  //G4cout << "angle " << angle << G4endl;
578 
579  counter++;
580 
581  if ( counter > fNv) {
582  G4Exception("G4ExtrudedSolid::AddGeneralPolygonFacets",
583  "GeomSolids0003", FatalException,
584  "Triangularisation has failed.");
585  break;
586  }
587  }
588 
589  G4bool good = true;
590  std::vector< Vertex >::iterator it;
591  for ( it=verticesToBeDone.begin(); it != verticesToBeDone.end(); ++it )
592  {
593  // skip vertices of tested triangle
594  //
595  if ( it == c1 || it == c2 || it == c3 ) { continue; }
596 
597  if ( IsPointInside(c1->first, c2->first, c3->first, it->first) )
598  {
599  // G4cout << "Point " << it->second << " is inside" << G4endl;
600  good = false;
601 
602  // try next three consecutive vertices
603  //
604  c1 = c2;
605  c2 = c3;
606  ++c3;
607  if ( c3 == verticesToBeDone.end() ) { c3 = verticesToBeDone.begin(); }
608  break;
609  }
610  // else
611  // { G4cout << "Point " << it->second << " is outside" << G4endl; }
612  }
613  if ( good )
614  {
615  // all points are outside triangle, we can make a facet
616 
617  // G4cout << "Found triangle : "
618  // << c1->second << " " << c2->second
619  // << " " << c3->second << G4endl;
620 
621  G4bool result;
622  result = AddFacet( MakeDownFacet(c1->second, c2->second, c3->second) );
623  if ( ! result ) { return false; }
624 
625  result = AddFacet( MakeUpFacet(c1->second, c2->second, c3->second) );
626  if ( ! result ) { return false; }
627 
628  std::vector<G4int> triangle(3);
629  triangle[0] = c1->second;
630  triangle[1] = c2->second;
631  triangle[2] = c3->second;
632  fTriangles.push_back(triangle);
633 
634  // remove the ear point from verticesToBeDone
635  //
636  verticesToBeDone.erase(c2);
637  c1 = verticesToBeDone.begin();
638  c2 = c1+1;
639  c3 = c1+2;
640  }
641  }
642  return true;
643 }
644 
645 //_____________________________________________________________________________
646 
648 {
649  // Define facets
650 
651  G4bool good;
652 
653  // Decomposition of polygonal sides in the facets
654  //
655  if ( fNv == 3 )
656  {
657  good = AddFacet( new G4TriangularFacet( GetVertex(0, 0), GetVertex(0, 1),
658  GetVertex(0, 2), ABSOLUTE) );
659  if ( ! good ) { return false; }
660 
661  good = AddFacet( new G4TriangularFacet( GetVertex(fNz-1, 2), GetVertex(fNz-1, 1),
662  GetVertex(fNz-1, 0), ABSOLUTE) );
663  if ( ! good ) { return false; }
664 
665  std::vector<G4int> triangle(3);
666  triangle[0] = 0;
667  triangle[1] = 1;
668  triangle[2] = 2;
669  fTriangles.push_back(triangle);
670  }
671 
672  else if ( fNv == 4 )
673  {
674  good = AddFacet( new G4QuadrangularFacet( GetVertex(0, 0),GetVertex(0, 1),
675  GetVertex(0, 2),GetVertex(0, 3),
676  ABSOLUTE) );
677  if ( ! good ) { return false; }
678 
679  good = AddFacet( new G4QuadrangularFacet( GetVertex(fNz-1, 3), GetVertex(fNz-1, 2),
680  GetVertex(fNz-1, 1), GetVertex(fNz-1, 0),
681  ABSOLUTE) );
682  if ( ! good ) { return false; }
683 
684  std::vector<G4int> triangle1(3);
685  triangle1[0] = 0;
686  triangle1[1] = 1;
687  triangle1[2] = 2;
688  fTriangles.push_back(triangle1);
689 
690  std::vector<G4int> triangle2(3);
691  triangle2[0] = 0;
692  triangle2[1] = 2;
693  triangle2[2] = 3;
694  fTriangles.push_back(triangle2);
695  }
696  else
697  {
698  good = AddGeneralPolygonFacets();
699  if ( ! good ) { return false; }
700  }
701 
702  // The quadrangular sides
703  //
704  for ( G4int iz = 0; iz < fNz-1; ++iz )
705  {
706  for ( G4int i = 0; i < fNv; ++i )
707  {
708  G4int j = (i+1) % fNv;
709  good = AddFacet( new G4QuadrangularFacet
710  ( GetVertex(iz, j), GetVertex(iz, i),
711  GetVertex(iz+1, i), GetVertex(iz+1, j), ABSOLUTE) );
712  if ( ! good ) { return false; }
713  }
714  }
715 
716  SetSolidClosed(true);
717 
718  return good;
719 }
720 
721 //_____________________________________________________________________________
722 
724 {
725  // Get polygon convexity (polygon is convex if all vertex angles are < pi )
726 
727  for ( G4int i=0; i< fNv; ++i )
728  {
729  G4int j = ( i + 1 ) % fNv;
730  G4int k = ( i + 2 ) % fNv;
731  G4TwoVector v1 = fPolygon[i]-fPolygon[j];
732  G4TwoVector v2 = fPolygon[k]-fPolygon[j];
733  G4double dphi = v2.phi() - v1.phi();
734  if ( dphi < 0. ) { dphi += 2.*pi; }
735 
736  if ( dphi >= pi ) { return false; }
737  }
738 
739  return true;
740 }
741 
742 //_____________________________________________________________________________
743 
745 {
746  // Return entity type
747 
748  return fGeometryType;
749 }
750 
751 //_____________________________________________________________________________
752 
754 {
755  return new G4ExtrudedSolid(*this);
756 }
757 
758 //_____________________________________________________________________________
759 
761 {
762  // Override the base class function as it fails in case of concave polygon.
763  // Project the point in the original polygon scale and check if it is inside
764  // for each triangle.
765 
766  // Check first if outside extent
767  //
768  if ( p.x() < GetMinXExtent() - kCarTolerance * 0.5 ||
769  p.x() > GetMaxXExtent() + kCarTolerance * 0.5 ||
770  p.y() < GetMinYExtent() - kCarTolerance * 0.5 ||
771  p.y() > GetMaxYExtent() + kCarTolerance * 0.5 ||
772  p.z() < GetMinZExtent() - kCarTolerance * 0.5 ||
773  p.z() > GetMaxZExtent() + kCarTolerance * 0.5 )
774  {
775  // G4cout << "G4ExtrudedSolid::Outside extent: " << p << G4endl;
776  return kOutside;
777  }
778 
779  // Project point p(z) to the polygon scale p0
780  //
781  G4TwoVector pscaled = ProjectPoint(p);
782 
783  // Check if on surface of polygon
784  //
785  for ( G4int i=0; i<fNv; ++i )
786  {
787  G4int j = (i+1) % fNv;
788  if ( IsSameLineSegment(pscaled, fPolygon[i], fPolygon[j]) )
789  {
790  // G4cout << "G4ExtrudedSolid::Inside return Surface (on polygon) "
791  // << G4endl;
792 
793  return kSurface;
794  }
795  }
796 
797  // Now check if inside triangles
798  //
799  std::vector< std::vector<G4int> >::const_iterator it = fTriangles.begin();
800  G4bool inside = false;
801  do
802  {
803  if ( IsPointInside(fPolygon[(*it)[0]], fPolygon[(*it)[1]],
804  fPolygon[(*it)[2]], pscaled) ) { inside = true; }
805  ++it;
806  } while ( (inside == false) && (it != fTriangles.end()) );
807 
808  if ( inside )
809  {
810  // Check if on surface of z sides
811  //
812  if ( std::fabs( p.z() - fZSections[0].fZ ) < kCarTolerance * 0.5 ||
813  std::fabs( p.z() - fZSections[fNz-1].fZ ) < kCarTolerance * 0.5 )
814  {
815  // G4cout << "G4ExtrudedSolid::Inside return Surface (on z side)"
816  // << G4endl;
817 
818  return kSurface;
819  }
820 
821  // G4cout << "G4ExtrudedSolid::Inside return Inside" << G4endl;
822 
823  return kInside;
824  }
825 
826  // G4cout << "G4ExtrudedSolid::Inside return Outside " << G4endl;
827 
828  return kOutside;
829 }
830 
831 //_____________________________________________________________________________
832 
834  const G4ThreeVector &v,
835  const G4bool calcNorm,
836  G4bool *validNorm,
837  G4ThreeVector *n) const
838 {
839  // Override the base class function to redefine validNorm
840  // (the solid can be concave)
841 
842  G4double distOut =
843  G4TessellatedSolid::DistanceToOut(p, v, calcNorm, validNorm, n);
844  if (validNorm) { *validNorm = fIsConvex; }
845 
846  return distOut;
847 }
848 
849 
850 //_____________________________________________________________________________
851 
853 {
854  // Override the overloaded base class function
855 
857 }
858 
859 //_____________________________________________________________________________
860 
861 std::ostream& G4ExtrudedSolid::StreamInfo(std::ostream &os) const
862 {
863  G4int oldprc = os.precision(16);
864  os << "-----------------------------------------------------------\n"
865  << " *** Dump for solid - " << GetName() << " ***\n"
866  << " ===================================================\n"
867  << " Solid geometry type: " << fGeometryType << G4endl;
868 
869  if ( fIsConvex)
870  { os << " Convex polygon; list of vertices:" << G4endl; }
871  else
872  { os << " Concave polygon; list of vertices:" << G4endl; }
873 
874  for ( G4int i=0; i<fNv; ++i )
875  {
876  os << std::setw(5) << "#" << i
877  << " vx = " << fPolygon[i].x()/mm << " mm"
878  << " vy = " << fPolygon[i].y()/mm << " mm" << G4endl;
879  }
880 
881  os << " Sections:" << G4endl;
882  for ( G4int iz=0; iz<fNz; ++iz )
883  {
884  os << " z = " << fZSections[iz].fZ/mm << " mm "
885  << " x0= " << fZSections[iz].fOffset.x()/mm << " mm "
886  << " y0= " << fZSections[iz].fOffset.y()/mm << " mm "
887  << " scale= " << fZSections[iz].fScale << G4endl;
888  }
889 
890 /*
891  // Triangles (for debugging)
892  os << G4endl;
893  os << " Triangles:" << G4endl;
894  os << " Triangle # vertex1 vertex2 vertex3" << G4endl;
895 
896  G4int counter = 0;
897  std::vector< std::vector<G4int> >::const_iterator it;
898  for ( it = fTriangles.begin(); it != fTriangles.end(); it++ ) {
899  std::vector<G4int> triangle = *it;
900  os << std::setw(10) << counter++
901  << std::setw(10) << triangle[0] << std::setw(10) << triangle[1] << std::setw(10) << triangle[2]
902  << G4endl;
903  }
904 */
905  os.precision(oldprc);
906 
907  return os;
908 }
G4String GetName() const
static c2_factory< G4double > c2
G4bool IsPointInside(G4TwoVector a, G4TwoVector b, G4TwoVector c, G4TwoVector p) const
G4GeometryType GetEntityType() const
void SetSolidClosed(const G4bool t)
G4bool IsConvex() const
CLHEP::Hep3Vector G4ThreeVector
G4VSolid * Clone() const
std::vector< G4TwoVector > fOffset0s
std::vector< G4double > fScale0s
G4ExtrudedSolid(const G4String &pName, std::vector< G4TwoVector > polygon, std::vector< ZSection > zsections)
const G4double pi
G4double GetMaxXExtent() const
std::vector< ZSection > fZSections
G4double GetMinXExtent() const
G4double a
Definition: TRTMaterials.hh:39
virtual ~G4ExtrudedSolid()
EInside Inside(const G4ThreeVector &p) const
G4double GetMaxZExtent() const
int G4int
Definition: G4Types.hh:78
G4TwoVector ProjectPoint(const G4ThreeVector &point) const
G4bool AddGeneralPolygonFacets()
G4TwoVector GetVertex(G4int index) const
G4bool IsSameLine(G4TwoVector p, G4TwoVector l1, G4TwoVector l2) const
G4double GetAngle(G4TwoVector p0, G4TwoVector pa, G4TwoVector pb) const
G4double GetMinZExtent() const
void ComputeProjectionParameters()
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const
virtual G4double DistanceToOut(const G4ThreeVector &p) const
bool G4bool
Definition: G4Types.hh:79
G4double iz
Definition: TRTMaterials.hh:39
static const G4double c3
G4bool AddFacet(G4VFacet *aFacet)
std::vector< G4double > fKScales
const G4int n
static const G4double c1
std::ostream & StreamInfo(std::ostream &os) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
T max(const T t1, const T t2)
brief Return the largest of the two arguments
EInside
Definition: geomdefs.hh:58
G4bool IsSameLineSegment(G4TwoVector p, G4TwoVector l1, G4TwoVector l2) const
std::vector< G4TwoVector > fKOffsets
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4bool IsSameSide(G4TwoVector p1, G4TwoVector p2, G4TwoVector l1, G4TwoVector l2) const
CLHEP::Hep2Vector G4TwoVector
Definition: G4TwoVector.hh:42
#define G4endl
Definition: G4ios.hh:61
G4double kCarTolerance
Definition: G4VSolid.hh:305
double G4double
Definition: G4Types.hh:76
G4VFacet * MakeUpFacet(G4int ind1, G4int ind2, G4int ind3) const
G4GeometryType fGeometryType
static const double mm
Definition: G4SIunits.hh:102
G4double GetMaxYExtent() const
G4ExtrudedSolid & operator=(const G4ExtrudedSolid &rhs)
G4VFacet * MakeDownFacet(G4int ind1, G4int ind2, G4int ind3) const
G4double GetMinYExtent() const
std::vector< std::vector< G4int > > fTriangles
G4TessellatedSolid & operator=(const G4TessellatedSolid &right)
std::vector< G4TwoVector > fPolygon