Geant4_10
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 
146  ComputeProjectionParameters();
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 
221  ComputeProjectionParameters();
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;
262  fTriangles = rhs.fTriangles; fIsConvex = rhs.fIsConvex;
263  fGeometryType = rhs.fGeometryType; fKScales = rhs.fKScales;
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 
279 void G4ExtrudedSolid::ComputeProjectionParameters()
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 
326 G4TwoVector G4ExtrudedSolid::ProjectPoint(const G4ThreeVector& point) const
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 
356 G4bool G4ExtrudedSolid::IsSameLine(G4TwoVector p,
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 
385 G4bool G4ExtrudedSolid::IsSameLineSegment(G4TwoVector p,
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 
404 G4bool G4ExtrudedSolid::IsSameSide(G4TwoVector p1, G4TwoVector p2,
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 
417 G4bool G4ExtrudedSolid::IsPointInside(G4TwoVector a, G4TwoVector b,
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
446 G4ExtrudedSolid::GetAngle(G4TwoVector po, G4TwoVector pa, G4TwoVector pb) const
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*
463 G4ExtrudedSolid::MakeDownFacet(G4int ind1, G4int ind2, G4int ind3) const
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*
497 G4ExtrudedSolid::MakeUpFacet(G4int ind1, G4int ind2, G4int ind3) const
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 
530 G4bool G4ExtrudedSolid::AddGeneralPolygonFacets()
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 
647 G4bool G4ExtrudedSolid::MakeFacets()
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 
723 G4bool G4ExtrudedSolid::IsConvex() const
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
double y() const
TTree * t1
Definition: plottest35.C:26
G4GeometryType GetEntityType() const
Float_t tmp
Definition: plot.C:37
void SetSolidClosed(const G4bool t)
double x() const
tuple a
Definition: test.py:11
CLHEP::Hep3Vector G4ThreeVector
G4VSolid * Clone() const
double x() const
G4ExtrudedSolid(const G4String &pName, std::vector< G4TwoVector > polygon, std::vector< ZSection > zsections)
const char * p
Definition: xmltok.h:285
G4double GetMaxXExtent() const
G4double G4NeutronHPJENDLHEData::G4double result
TCanvas * c1
Definition: plotHisto.C:7
G4double GetMinXExtent() const
tuple x
Definition: test.py:50
virtual ~G4ExtrudedSolid()
EInside Inside(const G4ThreeVector &p) const
G4double GetMaxZExtent() const
int G4int
Definition: G4Types.hh:78
double z() const
G4TwoVector GetVertex(G4int index) const
tuple b
Definition: test.py:12
Char_t n[5]
G4double GetMinZExtent() const
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
G4bool AddFacet(G4VFacet *aFacet)
tuple v
Definition: test.py:18
std::ostream & StreamInfo(std::ostream &os) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
double phi() const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
EInside
Definition: geomdefs.hh:58
double y() const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
TTree * t2
Definition: plottest35.C:36
#define G4endl
Definition: G4ios.hh:61
G4double kCarTolerance
Definition: G4VSolid.hh:305
double G4double
Definition: G4Types.hh:76
tuple c
Definition: test.py:13
G4double GetMaxYExtent() const
G4ExtrudedSolid & operator=(const G4ExtrudedSolid &rhs)
G4double GetMinYExtent() const
G4TessellatedSolid & operator=(const G4TessellatedSolid &right)