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