| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // G4ExtrudedSolid implementation
27 //
28 // Author: Ivana Hrivnacova, IPN Orsay
29 //
30 // CHANGE HISTORY
31 // --------------
32 //
33 // 31.10.2017 E.Tcherniaev: added implementati
34 // right prism
35 // 08.09.2017 E.Tcherniaev: added implementati
36 // right prism
37 // 21.10.2016 E.Tcherniaev: reimplemented Calc
38 // used G4GeomTools::PolygonArea()
39 // replaced IsConvex() with G4GeomT
40 // 02.03.2016 E.Tcherniaev: added CheckPolygon
41 // collinear and coincident points
42 // -------------------------------------------
43
44 #include "G4ExtrudedSolid.hh"
45
46 #if !defined(G4GEOM_USE_UEXTRUDEDSOLID)
47
48 #include <set>
49 #include <algorithm>
50 #include <cmath>
51 #include <iomanip>
52
53 #include "G4GeomTools.hh"
54 #include "G4VoxelLimits.hh"
55 #include "G4AffineTransform.hh"
56 #include "G4BoundingEnvelope.hh"
57
58 #include "G4GeometryTolerance.hh"
59 #include "G4PhysicalConstants.hh"
60 #include "G4SystemOfUnits.hh"
61 #include "G4TriangularFacet.hh"
62 #include "G4QuadrangularFacet.hh"
63
64 //____________________________________________
65
66 G4ExtrudedSolid::G4ExtrudedSolid( const G4Stri
67 const std::v
68 const std::v
69 : G4TessellatedSolid(pName),
70 fNv(polygon.size()),
71 fNz(zsections.size()),
72 fIsConvex(false),
73 fGeometryType("G4ExtrudedSolid"),
74 fSolidType(0)
75 {
76 // General constructor
77
78 // First check input parameters
79
80 if (fNv < 3)
81 {
82 std::ostringstream message;
83 message << "Number of vertices in polygon
84 G4Exception("G4ExtrudedSolid::G4ExtrudedSo
85 FatalErrorInArgument, message)
86 }
87
88 if (fNz < 2)
89 {
90 std::ostringstream message;
91 message << "Number of z-sides < 2 - " << p
92 G4Exception("G4ExtrudedSolid::G4ExtrudedSo
93 FatalErrorInArgument, message)
94 }
95
96 for ( std::size_t i=0; i<fNz-1; ++i )
97 {
98 if ( zsections[i].fZ > zsections[i+1].fZ )
99 {
100 std::ostringstream message;
101 message << "Z-sections have to be ordere
102 << pName;
103 G4Exception("G4ExtrudedSolid::G4Extruded
104 FatalErrorInArgument, messag
105 }
106 if ( std::fabs( zsections[i+1].fZ - zsecti
107 {
108 std::ostringstream message;
109 message << "Z-sections with the same z p
110 << pName;
111 G4Exception("G4ExtrudedSolid::G4Extruded
112 FatalException, message);
113 }
114 }
115
116 // Copy polygon
117 //
118 fPolygon = polygon;
119
120 // Remove collinear and coincident vertices,
121 //
122 std::vector<G4int> removedVertices;
123 G4GeomTools::RemoveRedundantVertices(fPolygo
124 2*kCarT
125 if (!removedVertices.empty())
126 {
127 std::size_t nremoved = removedVertices.siz
128 std::ostringstream message;
129 message << "The following "<< nremoved
130 << " vertices have been removed fr
131 << "\nas collinear or coincident w
132 << removedVertices[0];
133 for (std::size_t i=1; i<nremoved; ++i) mes
134 G4Exception("G4ExtrudedSolid::G4ExtrudedSo
135 JustWarning, message);
136 }
137
138 fNv = fPolygon.size();
139 if (fNv < 3)
140 {
141 std::ostringstream message;
142 message << "Number of vertices in polygon
143 G4Exception("G4ExtrudedSolid::G4ExtrudedSo
144 FatalErrorInArgument, message)
145 }
146
147 // Check if polygon vertices are defined clo
148 // (the area is positive if polygon vertices
149 //
150 if (G4GeomTools::PolygonArea(fPolygon) > 0.)
151 {
152 // Polygon vertices are defined anti-clock
153 // G4Exception("G4ExtrudedSolid::G4Extrude
154 // JustWarning,
155 // "Polygon vertices defined an
156 std::reverse(fPolygon.begin(),fPolygon.end
157 }
158
159 // Copy z-sections
160 //
161 fZSections = zsections;
162
163 G4bool result = MakeFacets();
164 if (!result)
165 {
166 std::ostringstream message;
167 message << "Making facets failed - " << pN
168 G4Exception("G4ExtrudedSolid::G4ExtrudedSo
169 FatalException, message);
170 }
171 fIsConvex = G4GeomTools::IsConvex(fPolygon);
172
173 ComputeProjectionParameters();
174
175 // Check if the solid is a right prism, if s
176 //
177 if ((fNz == 2)
178 && (fZSections[0].fScale == 1) && (fZSec
179 && (fZSections[0].fOffset == G4TwoVector
180 && (fZSections[1].fOffset == G4TwoVector
181 {
182 fSolidType = (fIsConvex) ? 1 : 2; // 1 - c
183 ComputeLateralPlanes();
184 }
185 }
186
187 //____________________________________________
188
189 G4ExtrudedSolid::G4ExtrudedSolid( const G4Stri
190 const std::v
191 G4doub
192 const G4TwoV
193 const G4TwoV
194 : G4TessellatedSolid(pName),
195 fNv(polygon.size()),
196 fNz(2),
197 fGeometryType("G4ExtrudedSolid")
198 {
199 // Special constructor for solid with 2 z-se
200
201 // First check input parameters
202 //
203 if (fNv < 3)
204 {
205 std::ostringstream message;
206 message << "Number of vertices in polygon
207 G4Exception("G4ExtrudedSolid::G4ExtrudedSo
208 FatalErrorInArgument, message)
209 }
210
211 // Copy polygon
212 //
213 fPolygon = polygon;
214
215 // Remove collinear and coincident vertices,
216 //
217 std::vector<G4int> removedVertices;
218 G4GeomTools::RemoveRedundantVertices(fPolygo
219 2*kCarT
220 if (!removedVertices.empty())
221 {
222 std::size_t nremoved = removedVertices.siz
223 std::ostringstream message;
224 message << "The following "<< nremoved
225 << " vertices have been removed fr
226 << "\nas collinear or coincident w
227 << removedVertices[0];
228 for (std::size_t i=1; i<nremoved; ++i) mes
229 G4Exception("G4ExtrudedSolid::G4ExtrudedSo
230 JustWarning, message);
231 }
232
233 fNv = fPolygon.size();
234 if (fNv < 3)
235 {
236 std::ostringstream message;
237 message << "Number of vertices in polygon
238 G4Exception("G4ExtrudedSolid::G4ExtrudedSo
239 FatalErrorInArgument, message)
240 }
241
242 // Check if polygon vertices are defined clo
243 // (the area is positive if polygon vertices
244 //
245 if (G4GeomTools::PolygonArea(fPolygon) > 0.)
246 {
247 // Polygon vertices are defined anti-clock
248 // G4Exception("G4ExtrudedSolid::G4Extrude
249 // JustWarning,
250 // "Polygon vertices defined an
251 std::reverse(fPolygon.begin(),fPolygon.end
252 }
253
254 // Copy z-sections
255 //
256 fZSections.emplace_back(-dz, off1, scale1);
257 fZSections.emplace_back( dz, off2, scale2);
258
259 G4bool result = MakeFacets();
260 if (!result)
261 {
262 std::ostringstream message;
263 message << "Making facets failed - " << pN
264 G4Exception("G4ExtrudedSolid::G4ExtrudedSo
265 FatalException, message);
266 }
267 fIsConvex = G4GeomTools::IsConvex(fPolygon);
268
269 ComputeProjectionParameters();
270
271 // Check if the solid is a right prism, if s
272 //
273 if ((scale1 == 1) && (scale2 == 1)
274 && (off1 == G4TwoVector(0,0)) && (off2 =
275 {
276 fSolidType = (fIsConvex) ? 1 : 2; // 1 - c
277 ComputeLateralPlanes();
278 }
279 }
280
281 //____________________________________________
282
283 G4ExtrudedSolid::G4ExtrudedSolid( __void__& a
284 : G4TessellatedSolid(a), fGeometryType("G4Ex
285 {
286 // Fake default constructor - sets only memb
287 // for usage rest
288 }
289
290 //____________________________________________
291
292 G4ExtrudedSolid::G4ExtrudedSolid(const G4Extru
293
294 //____________________________________________
295
296 G4ExtrudedSolid& G4ExtrudedSolid::operator = (
297 {
298 // Check assignment to self
299 //
300 if (this == &rhs) { return *this; }
301
302 // Copy base class data
303 //
304 G4TessellatedSolid::operator=(rhs);
305
306 // Copy data
307 //
308 fNv = rhs.fNv; fNz = rhs.fNz;
309 fPolygon = rhs.fPolygon; fZSections = rhs.f
310 fTriangles = rhs.fTriangles; fIsConvex = rh
311 fGeometryType = rhs.fGeometryType;
312 fSolidType = rhs.fSolidType; fPlanes = rhs.
313 fLines = rhs.fLines; fLengths = rhs.fLength
314 fKScales = rhs.fKScales; fScale0s = rhs.fSc
315 fKOffsets = rhs.fKOffsets; fOffset0s = rhs.
316
317 return *this;
318 }
319
320 //____________________________________________
321
322 G4ExtrudedSolid::~G4ExtrudedSolid()
323 {
324 // Destructor
325 }
326
327 //____________________________________________
328
329 void G4ExtrudedSolid::ComputeProjectionParamet
330 {
331 // Compute parameters for point projections
332 // to the polygon scale & offset:
333 // scale(z) = k*z + scale0
334 // offset(z) = l*z + offset0
335 // p(z) = scale(z)*p0 + offset(z)
336 // p0 = (p(z) - offset(z))/scale(z);
337 //
338
339 for (std::size_t iz=0; iz<fNz-1; ++iz)
340 {
341 G4double z1 = fZSections[iz].fZ;
342 G4double z2 = fZSections[iz+1].fZ;
343 G4double scale1 = fZSections[iz].fScale;
344 G4double scale2 = fZSections[iz+1].fScale
345 G4TwoVector off1 = fZSections[iz].fOffset;
346 G4TwoVector off2 = fZSections[iz+1].fOffse
347
348 G4double kscale = (scale2 - scale1)/(z2 -
349 G4double scale0 = scale2 - kscale*(z2 - z
350 G4TwoVector koff = (off2 - off1)/(z2 - z1)
351 G4TwoVector off0 = off2 - koff*(z2 - z1)/
352
353 fKScales.push_back(kscale);
354 fScale0s.push_back(scale0);
355 fKOffsets.push_back(koff);
356 fOffset0s.push_back(off0);
357 }
358 }
359
360 //____________________________________________
361
362 void G4ExtrudedSolid::ComputeLateralPlanes()
363 {
364 // Compute lateral planes: a*x + b*y + c*z +
365 //
366 std::size_t Nv = fPolygon.size();
367 fPlanes.resize(Nv);
368 for (std::size_t i=0, k=Nv-1; i<Nv; k=i++)
369 {
370 G4TwoVector norm = (fPolygon[i] - fPolygon
371 fPlanes[i].a = -norm.y();
372 fPlanes[i].b = norm.x();
373 fPlanes[i].c = 0;
374 fPlanes[i].d = norm.y()*fPolygon[i].x() -
375 }
376
377 // Compute edge equations: x = k*y + m
378 // and edge lengths
379 //
380 fLines.resize(Nv);
381 fLengths.resize(Nv);
382 for (std::size_t i=0, k=Nv-1; i<Nv; k=i++)
383 {
384 if (fPolygon[k].y() == fPolygon[i].y())
385 {
386 fLines[i].k = 0;
387 fLines[i].m = fPolygon[i].x();
388 }
389 else
390 {
391 G4double ctg = (fPolygon[k].x()-fPolygon
392 fLines[i].k = ctg;
393 fLines[i].m = fPolygon[i].x() - ctg*fPol
394 }
395 fLengths[i] = (fPolygon[i] - fPolygon[k]
396 }
397 }
398
399 //____________________________________________
400
401 G4ThreeVector G4ExtrudedSolid::GetVertex(G4int
402 {
403 // Shift and scale vertices
404
405 return { fPolygon[ind].x() * fZSections[iz].
406 + fZSections[iz].fOffset.x(),
407 fPolygon[ind].y() * fZSections[iz].
408 + fZSections[iz].fOffset.y(),
409 fZSections[iz].fZ };
410 }
411
412 //____________________________________________
413
414 G4TwoVector G4ExtrudedSolid::ProjectPoint(cons
415 {
416 // Project point in the polygon scale
417 // scale(z) = k*z + scale0
418 // offset(z) = l*z + offset0
419 // p(z) = scale(z)*p0 + offset(z)
420 // p0 = (p(z) - offset(z))/scale(z);
421
422 // Select projection (z-segment of the solid
423 //
424 std::size_t iz = 0;
425 while ( point.z() > fZSections[iz+1].fZ && i
426 // Loop checking, 13.08.2015, G.Cosmo
427
428 G4double z0 = ( fZSections[iz+1].fZ + fZSect
429 G4TwoVector p2(point.x(), point.y());
430 G4double pscale = fKScales[iz]*(point.z()-z
431 G4TwoVector poffset = fKOffsets[iz]*(point.z
432
433 // G4cout << point << " projected to "
434 // << iz << "-th z-segment polygon as
435 // << (p2 - poffset)/pscale << G4endl
436
437 // pscale is always >0 as it is an interpola
438 // positive scale values
439 //
440 return (p2 - poffset)/pscale;
441 }
442
443 //____________________________________________
444
445 G4bool G4ExtrudedSolid::IsSameLine(const G4Two
446 const G4Two
447 const G4Two
448 {
449 // Return true if p is on the line through l
450
451 if ( l1.x() == l2.x() )
452 {
453 return std::fabs(p.x() - l1.x()) < kCarTol
454 }
455 G4double slope= ((l2.y() - l1.y())/(l2.x()
456 G4double predy= l1.y() + slope *(p.x() - l
457 G4double dy= p.y() - predy;
458
459 // Calculate perpendicular distance
460 //
461 // G4double perpD= std::fabs(dy) / std::sqr
462 // G4bool simpleComp= (perpD<kCarToleranc
463
464 // Check perpendicular distance vs toleranc
465 //
466 G4bool squareComp = (dy*dy < (1+slope*slope
467 * kCarToleranceHalf * kCa
468
469 // return simpleComp;
470 return squareComp;
471 }
472
473 //____________________________________________
474
475 G4bool G4ExtrudedSolid::IsSameLineSegment(cons
476 cons
477 cons
478 {
479 // Return true if p is on the line through l
480 // l1 and l2
481
482 if ( p.x() < std::min(l1.x(), l2.x()) - kCar
483 p.x() > std::max(l1.x(), l2.x()) + kCar
484 p.y() < std::min(l1.y(), l2.y()) - kCar
485 p.y() > std::max(l1.y(), l2.y()) + kCar
486 {
487 return false;
488 }
489
490 return IsSameLine(p, l1, l2);
491 }
492
493 //____________________________________________
494
495 G4bool G4ExtrudedSolid::IsSameSide(const G4Two
496 const G4Two
497 const G4Two
498 const G4Two
499 {
500 // Return true if p1 and p2 are on the same
501
502 return ( (p1.x() - l1.x()) * (l2.y() - l1.
503 - (l2.x() - l1.x()) * (p1.y() - l1.y(
504 * ( (p2.x() - l1.x()) * (l2.y() - l1.
505 - (l2.x() - l1.x()) * (p2.y() - l1.y(
506 }
507
508 //____________________________________________
509
510 G4bool G4ExtrudedSolid::IsPointInside(const G4
511 const G4
512 const G4
513 const G4
514 {
515 // Return true if p is inside of triangle ab
516 // else returns false
517
518 // Check extent first
519 //
520 if ( ( p.x() < a.x() && p.x() < b.x() && p.x
521 ( p.x() > a.x() && p.x() > b.x() && p.x
522 ( p.y() < a.y() && p.y() < b.y() && p.y
523 ( p.y() > a.y() && p.y() > b.y() && p.y
524
525 G4bool inside
526 = IsSameSide(p, a, b, c)
527 && IsSameSide(p, b, a, c)
528 && IsSameSide(p, c, a, b);
529
530 G4bool onEdge
531 = IsSameLineSegment(p, a, b)
532 || IsSameLineSegment(p, b, c)
533 || IsSameLineSegment(p, c, a);
534
535 return inside || onEdge;
536 }
537
538 //____________________________________________
539
540 G4double
541 G4ExtrudedSolid::GetAngle(const G4TwoVector& p
542 const G4TwoVector& p
543 const G4TwoVector& p
544 {
545 // Return the angle of the vertex in po
546
547 G4TwoVector t1 = pa - po;
548 G4TwoVector t2 = pb - po;
549
550 G4double result = (std::atan2(t1.y(), t1.x()
551
552 if ( result < 0 ) result += 2*pi;
553
554 return result;
555 }
556
557 //____________________________________________
558
559 G4VFacet*
560 G4ExtrudedSolid::MakeDownFacet(G4int ind1, G4i
561 {
562 // Create a triangular facet from the polygo
563 // forming the down side ( the normal goes i
564
565 std::vector<G4ThreeVector> vertices;
566 vertices.push_back(GetVertex(0, ind1));
567 vertices.push_back(GetVertex(0, ind2));
568 vertices.push_back(GetVertex(0, ind3));
569
570 // first vertex most left
571 //
572 G4ThreeVector cross
573 = (vertices[1]-vertices[0]).cross(vertices
574
575 if ( cross.z() > 0.0 )
576 {
577 // vertices ordered clock wise has to be r
578
579 // G4cout << "G4ExtrudedSolid::MakeDownFac
580 // << ind1 << ", " << ind2 << ", "
581
582 G4ThreeVector tmp = vertices[1];
583 vertices[1] = vertices[2];
584 vertices[2] = tmp;
585 }
586
587 return new G4TriangularFacet(vertices[0], ve
588 vertices[2], AB
589 }
590
591 //____________________________________________
592
593 G4VFacet*
594 G4ExtrudedSolid::MakeUpFacet(G4int ind1, G4int
595 {
596 // Creates a triangular facet from the polyg
597 // forming the upper side ( z>0 )
598
599 std::vector<G4ThreeVector> vertices;
600 vertices.push_back(GetVertex((G4int)fNz-1, i
601 vertices.push_back(GetVertex((G4int)fNz-1, i
602 vertices.push_back(GetVertex((G4int)fNz-1, i
603
604 // first vertex most left
605 //
606 G4ThreeVector cross
607 = (vertices[1]-vertices[0]).cross(vertices
608
609 if ( cross.z() < 0.0 )
610 {
611 // vertices ordered clock wise has to be r
612
613 // G4cout << "G4ExtrudedSolid::MakeUpFacet
614 // << ind1 << ", " << ind2 << ", "
615
616 G4ThreeVector tmp = vertices[1];
617 vertices[1] = vertices[2];
618 vertices[2] = tmp;
619 }
620
621 return new G4TriangularFacet(vertices[0], ve
622 vertices[2], AB
623 }
624
625 //____________________________________________
626
627 G4bool G4ExtrudedSolid::AddGeneralPolygonFacet
628 {
629 // Decompose polygonal sides in triangular f
630
631 typedef std::pair < G4TwoVector, G4int > Ver
632
633 static const G4double kAngTolerance =
634 G4GeometryTolerance::GetInstance()->GetAng
635
636 // Fill one more vector
637 //
638 std::vector< Vertex > verticesToBeDone;
639 for ( G4int i=0; i<(G4int)fNv; ++i )
640 {
641 verticesToBeDone.emplace_back(fPolygon[i],
642 }
643 std::vector< Vertex > ears;
644
645 auto c1 = verticesToBeDone.begin();
646 auto c2 = c1+1;
647 auto c3 = c1+2;
648 while ( verticesToBeDone.size()>2 ) // Lo
649 {
650
651 // G4cout << "Looking at triangle : "
652 // << c1->second << " " << c2->se
653 // << " " << c3->second << G4endl;
654 //G4cout << "Looking at triangle : "
655 // << c1->first << " " << c2->firs
656 // << " " << c3->first << G4endl;
657
658 // skip concave vertices
659 //
660 G4double angle = GetAngle(c2->first, c3->f
661
662 //G4cout << "angle " << angle << G4endl;
663
664 std::size_t counter = 0;
665 while ( angle >= (pi-kAngTolerance) ) //
666 {
667 // G4cout << "Skipping concave vertex "
668
669 // try next three consecutive vertices
670 //
671 c1 = c2;
672 c2 = c3;
673 ++c3;
674 if ( c3 == verticesToBeDone.end() ) { c3
675
676 //G4cout << "Looking at triangle : "
677 // << c1->first << " " << c2->firs
678 // << " " << c3->first << G4endl
679
680 angle = GetAngle(c2->first, c3->first, c
681 //G4cout << "angle " << angle << G4endl
682
683 ++counter;
684
685 if ( counter > fNv )
686 {
687 G4Exception("G4ExtrudedSolid::AddGener
688 "GeomSolids0003", FatalExc
689 "Triangularisation has fai
690 break;
691 }
692 }
693
694 G4bool good = true;
695 for ( auto it=verticesToBeDone.cbegin(); i
696 {
697 // skip vertices of tested triangle
698 //
699 if ( it == c1 || it == c2 || it == c3 )
700
701 if ( IsPointInside(c1->first, c2->first,
702 {
703 // G4cout << "Point " << it->second <<
704 good = false;
705
706 // try next three consecutive vertices
707 //
708 c1 = c2;
709 c2 = c3;
710 ++c3;
711 if ( c3 == verticesToBeDone.end() ) {
712 break;
713 }
714 // else
715 // { G4cout << "Point " << it->second
716 }
717 if ( good )
718 {
719 // all points are outside triangle, we c
720
721 // G4cout << "Found triangle : "
722 // << c1->second << " " << c2->s
723 // << " " << c3->second << G4end
724
725 G4bool result;
726 result = AddFacet( MakeDownFacet(c1->sec
727 if ( ! result ) { return false; }
728
729 result = AddFacet( MakeUpFacet(c1->secon
730 if ( ! result ) { return false; }
731
732 std::vector<G4int> triangle(3);
733 triangle[0] = c1->second;
734 triangle[1] = c2->second;
735 triangle[2] = c3->second;
736 fTriangles.push_back(std::move(triangle)
737
738 // remove the ear point from verticesToB
739 //
740 verticesToBeDone.erase(c2);
741 c1 = verticesToBeDone.begin();
742 c2 = c1+1;
743 c3 = c1+2;
744 }
745 }
746 return true;
747 }
748
749 //____________________________________________
750
751 G4bool G4ExtrudedSolid::MakeFacets()
752 {
753 // Define facets
754
755 G4bool good;
756
757 // Decomposition of polygonal sides in the f
758 //
759 if ( fNv == 3 )
760 {
761 good = AddFacet( new G4TriangularFacet( Ge
762 Ge
763 if ( ! good ) { return false; }
764
765 good = AddFacet( new G4TriangularFacet( Ge
766 Ge
767 Ge
768 AB
769 if ( ! good ) { return false; }
770
771 std::vector<G4int> triangle(3);
772 triangle[0] = 0;
773 triangle[1] = 1;
774 triangle[2] = 2;
775 fTriangles.push_back(std::move(triangle));
776 }
777
778 else if ( fNv == 4 )
779 {
780 good = AddFacet( new G4QuadrangularFacet(
781
782
783 if ( ! good ) { return false; }
784
785 good = AddFacet( new G4QuadrangularFacet(
786
787
788
789
790 if ( ! good ) { return false; }
791
792 std::vector<G4int> triangle1(3);
793 triangle1[0] = 0;
794 triangle1[1] = 1;
795 triangle1[2] = 2;
796 fTriangles.push_back(std::move(triangle1))
797
798 std::vector<G4int> triangle2(3);
799 triangle2[0] = 0;
800 triangle2[1] = 2;
801 triangle2[2] = 3;
802 fTriangles.push_back(std::move(triangle2))
803 }
804 else
805 {
806 good = AddGeneralPolygonFacets();
807 if ( ! good ) { return false; }
808 }
809
810 // The quadrangular sides
811 //
812 for ( G4int iz = 0; iz < (G4int)fNz-1; ++iz
813 {
814 for ( G4int i = 0; i < (G4int)fNv; ++i )
815 {
816 G4int j = (i+1) % fNv;
817 good = AddFacet( new G4QuadrangularFacet
818 ( GetVertex(iz, j), Ge
819 GetVertex(iz+1, i),
820 if ( ! good ) { return false; }
821 }
822 }
823
824 SetSolidClosed(true);
825
826 return good;
827 }
828
829 //____________________________________________
830
831 G4GeometryType G4ExtrudedSolid::GetEntityType
832 {
833 // Return entity type
834
835 return fGeometryType;
836 }
837
838 //____________________________________________
839
840 G4bool G4ExtrudedSolid::IsFaceted () const
841 {
842 return true;
843 }
844
845 //____________________________________________
846
847 G4VSolid* G4ExtrudedSolid::Clone() const
848 {
849 return new G4ExtrudedSolid(*this);
850 }
851
852 //____________________________________________
853
854 EInside G4ExtrudedSolid::Inside(const G4ThreeV
855 {
856 switch (fSolidType)
857 {
858 case 1: // convex right prism
859 {
860 G4double dist = std::max(fZSections[0].f
861 if (dist > kCarToleranceHalf) { return
862
863 std::size_t np = fPlanes.size();
864 for (std::size_t i=0; i<np; ++i)
865 {
866 G4double dd = fPlanes[i].a*p.x() + fPl
867 if (dd > dist) { dist = dd; }
868 }
869 if (dist > kCarToleranceHalf) { return
870 return (dist > -kCarToleranceHalf) ? kSu
871 }
872 case 2: // non-convex right prism
873 {
874 G4double distz = std::max(fZSections[0].
875 if (distz > kCarToleranceHalf) { return
876
877 G4bool in = PointInPolygon(p);
878 if (distz > -kCarToleranceHalf && in) {
879
880 G4double dd = DistanceToPolygonSqr(p) -
881 if (in)
882 {
883 return (dd >= 0) ? kInside : kSurface;
884 }
885 else
886 {
887 return (dd > 0) ? kOutside : kSurface;
888 }
889 }
890 }
891
892 // Override the base class function as it f
893 // Project the point in the original polygon
894 // for each triangle.
895
896 // Check first if outside extent
897 //
898 if ( p.x() < GetMinXExtent() - kCarTolerance
899 p.x() > GetMaxXExtent() + kCarTolerance
900 p.y() < GetMinYExtent() - kCarTolerance
901 p.y() > GetMaxYExtent() + kCarTolerance
902 p.z() < GetMinZExtent() - kCarTolerance
903 p.z() > GetMaxZExtent() + kCarTolerance
904 {
905 // G4cout << "G4ExtrudedSolid::Outside ext
906 return kOutside;
907 }
908
909 // Project point p(z) to the polygon scale p
910 //
911 G4TwoVector pscaled = ProjectPoint(p);
912
913 // Check if on surface of polygon
914 //
915 for ( G4int i=0; i<(G4int)fNv; ++i )
916 {
917 G4int j = (i+1) % fNv;
918 if ( IsSameLineSegment(pscaled, fPolygon[i
919 {
920 // G4cout << "G4ExtrudedSolid::Inside re
921 // << G4endl;
922
923 return kSurface;
924 }
925 }
926
927 // Now check if inside triangles
928 //
929 auto it = fTriangles.cbegin();
930 G4bool inside = false;
931 do // Loop checking, 13.08.2015, G.Cosmo
932 {
933 if ( IsPointInside(fPolygon[(*it)[0]], fPo
934 fPolygon[(*it)[2]], psc
935 ++it;
936 } while ( (!inside) && (it != fTriangles.cen
937
938 if ( inside )
939 {
940 // Check if on surface of z sides
941 //
942 if ( std::fabs( p.z() - fZSections[0].fZ )
943 std::fabs( p.z() - fZSections[fNz-1].
944 {
945 // G4cout << "G4ExtrudedSolid::Inside re
946 // << G4endl;
947
948 return kSurface;
949 }
950
951 // G4cout << "G4ExtrudedSolid::Inside retu
952
953 return kInside;
954 }
955
956 // G4cout << "G4ExtrudedSolid::Inside return
957
958 return kOutside;
959 }
960
961 //____________________________________________
962
963 G4ThreeVector G4ExtrudedSolid::SurfaceNormal(c
964 {
965 G4int nsurf = 0;
966 G4double nx = 0., ny = 0., nz = 0.;
967 switch (fSolidType)
968 {
969 case 1: // convex right prism
970 {
971 if (std::abs(p.z() - fZSections[0].fZ) <
972 {
973 nz = -1; ++nsurf;
974 }
975 if (std::abs(p.z() - fZSections[1].fZ) <
976 {
977 nz = 1; ++nsurf;
978 }
979 for (std::size_t i=0; i<fNv; ++i)
980 {
981 G4double dd = fPlanes[i].a*p.x() + fPl
982 if (std::abs(dd) > kCarToleranceHalf)
983 nx += fPlanes[i].a;
984 ny += fPlanes[i].b;
985 ++nsurf;
986 }
987 break;
988 }
989 case 2: // non-convex right prism
990 {
991 if (std::abs(p.z() - fZSections[0].fZ) <
992 {
993 nz = -1; ++nsurf;
994 }
995 if (std::abs(p.z() - fZSections[1].fZ) <
996 {
997 nz = 1; ++nsurf;
998 }
999
1000 G4double sqrCarToleranceHalf = kCarTole
1001 for (std::size_t i=0, k=fNv-1; i<fNv; k
1002 {
1003 G4double ix = p.x() - fPolygon[i].x()
1004 G4double iy = p.y() - fPolygon[i].y()
1005 G4double u = fPlanes[i].a*iy - fPlan
1006 if (u < 0)
1007 {
1008 if (ix*ix + iy*iy > sqrCarTolerance
1009 }
1010 else if (u > fLengths[i])
1011 {
1012 G4double kx = p.x() - fPolygon[k].x
1013 G4double ky = p.y() - fPolygon[k].y
1014 if (kx*kx + ky*ky > sqrCarTolerance
1015 }
1016 else
1017 {
1018 G4double dd = fPlanes[i].a*p.x() +
1019 if (dd*dd > sqrCarToleranceHalf) co
1020 }
1021 nx += fPlanes[i].a;
1022 ny += fPlanes[i].b;
1023 ++nsurf;
1024 }
1025 break;
1026 }
1027 default:
1028 {
1029 return G4TessellatedSolid::SurfaceNorma
1030 }
1031 }
1032
1033 // Return normal (right prism)
1034 //
1035 if (nsurf == 1)
1036 {
1037 return { nx,ny,nz };
1038 }
1039 else if (nsurf != 0) // edge or corner
1040 {
1041 return G4ThreeVector(nx,ny,nz).unit();
1042 }
1043 else
1044 {
1045 // Point is not on the surface, compute a
1046 //
1047 #ifdef G4CSGDEBUG
1048 std::ostringstream message;
1049 G4long oldprc = message.precision(16);
1050 message << "Point p is not on surface (!?
1051 << GetName() << G4endl;
1052 message << "Position:\n";
1053 message << " p.x() = " << p.x()/mm << "
1054 message << " p.y() = " << p.y()/mm << "
1055 message << " p.z() = " << p.z()/mm << "
1056 G4cout.precision(oldprc) ;
1057 G4Exception("G4TesselatedSolid::SurfaceNo
1058 JustWarning, message );
1059 DumpInfo();
1060 #endif
1061 return ApproxSurfaceNormal(p);
1062 }
1063 }
1064
1065 //___________________________________________
1066
1067 G4ThreeVector G4ExtrudedSolid::ApproxSurfaceN
1068 {
1069 // This method is valid only for right pris
1070 // normally should not be called
1071
1072 if (fSolidType == 1 || fSolidType == 2)
1073 {
1074 // Find distances to z-planes
1075 //
1076 G4double dz0 = fZSections[0].fZ - p.z();
1077 G4double dz1 = p.z() - fZSections[1].fZ;
1078 G4double ddz0 = dz0*dz0;
1079 G4double ddz1 = dz1*dz1;
1080
1081 // Find nearest lateral side and distance
1082 //
1083 std::size_t iside = 0;
1084 G4double dd = DBL_MAX;
1085 for (std::size_t i=0, k=fNv-1; i<fNv; k=i
1086 {
1087 G4double ix = p.x() - fPolygon[i].x();
1088 G4double iy = p.y() - fPolygon[i].y();
1089 G4double u = fPlanes[i].a*iy - fPlanes
1090 if (u < 0)
1091 {
1092 G4double tmp = ix*ix + iy*iy;
1093 if (tmp < dd) { dd = tmp; iside = i;
1094 }
1095 else if (u > fLengths[i])
1096 {
1097 G4double kx = p.x() - fPolygon[k].x()
1098 G4double ky = p.y() - fPolygon[k].y()
1099 G4double tmp = kx*kx + ky*ky;
1100 if (tmp < dd) { dd = tmp; iside = i;
1101 }
1102 else
1103 {
1104 G4double tmp = fPlanes[i].a*p.x() + f
1105 tmp *= tmp;
1106 if (tmp < dd) { dd = tmp; iside = i;
1107 }
1108 }
1109
1110 // Find region
1111 //
1112 // 3 | 1 | 3
1113 // ----+-------+----
1114 // 2 | 0 | 2
1115 // ----+-------+----
1116 // 3 | 1 | 3
1117 //
1118 G4int iregion = 0;
1119 if (std::max(dz0,dz1) > 0) iregion = 1;
1120
1121 G4bool in = PointInPolygon(p);
1122 if (!in) iregion += 2;
1123
1124 // Return normal
1125 //
1126 switch (iregion)
1127 {
1128 case 0:
1129 {
1130 if (ddz0 <= ddz1 && ddz0 <= dd) retur
1131 if (ddz1 <= ddz0 && ddz1 <= dd) retur
1132 return { fPlanes[iside].a, fPlanes[is
1133 }
1134 case 1:
1135 {
1136 return { 0, 0, (G4double)((dz0 > dz1)
1137 }
1138 case 2:
1139 {
1140 return { fPlanes[iside].a, fPlanes[is
1141 }
1142 case 3:
1143 {
1144 G4double dzmax = std::max(dz0,dz1);
1145 if (dzmax*dzmax > dd) return { 0, 0,
1146 return { fPlanes[iside].a, fPlanes[is
1147 }
1148 }
1149 }
1150 return {0,0,0};
1151 }
1152
1153 //___________________________________________
1154
1155 G4double G4ExtrudedSolid::DistanceToIn(const
1156 const
1157 {
1158 G4double z0 = fZSections[0].fZ;
1159 G4double z1 = fZSections[fNz-1].fZ;
1160 if ((p.z() <= z0 + kCarToleranceHalf) && v.
1161 if ((p.z() >= z1 - kCarToleranceHalf) && v.
1162
1163 switch (fSolidType)
1164 {
1165 case 1: // convex right prism
1166 {
1167 // Intersection with Z planes
1168 //
1169 G4double dz = (z1 - z0)*0.5;
1170 G4double pz = p.z() - dz - z0;
1171
1172 G4double invz = (v.z() == 0) ? DBL_MAX
1173 G4double ddz = (invz < 0) ? dz : -dz;
1174 G4double tzmin = (pz + ddz)*invz;
1175 G4double tzmax = (pz - ddz)*invz;
1176
1177 // Intersection with lateral planes
1178 //
1179 std::size_t np = fPlanes.size();
1180 G4double txmin = tzmin, txmax = tzmax;
1181 for (std::size_t i=0; i<np; ++i)
1182 {
1183 G4double cosa = fPlanes[i].a*v.x()+fP
1184 G4double dist = fPlanes[i].a*p.x()+fP
1185 if (dist >= -kCarToleranceHalf)
1186 {
1187 if (cosa >= 0) { return kInfinity;
1188 G4double tmp = -dist/cosa;
1189 if (txmin < tmp) { txmin = tmp; }
1190 }
1191 else if (cosa > 0)
1192 {
1193 G4double tmp = -dist/cosa;
1194 if (txmax > tmp) { txmax = tmp; }
1195 }
1196 }
1197
1198 // Find distance
1199 //
1200 G4double tmin = txmin, tmax = txmax;
1201 if (tmax <= tmin + kCarToleranceHalf)
1202 {
1203 return kInfinity;
1204 }
1205 return (tmin < kCarToleranceHalf) ? 0.
1206 }
1207 case 2: // non-convex right prism
1208 {
1209 }
1210 }
1211 return G4TessellatedSolid::DistanceToIn(p,v
1212 }
1213
1214 //___________________________________________
1215
1216 G4double G4ExtrudedSolid::DistanceToIn (const
1217 {
1218 switch (fSolidType)
1219 {
1220 case 1: // convex right prism
1221 {
1222 G4double dist = std::max(fZSections[0].
1223 std::size_t np = fPlanes.size();
1224 for (std::size_t i=0; i<np; ++i)
1225 {
1226 G4double dd = fPlanes[i].a*p.x() + fP
1227 if (dd > dist) dist = dd;
1228 }
1229 return (dist > 0) ? dist : 0.;
1230 }
1231 case 2: // non-convex right prism
1232 {
1233 G4bool in = PointInPolygon(p);
1234 if (in)
1235 {
1236 G4double distz= std::max(fZSections[0
1237 return (distz > 0) ? distz : 0;
1238 }
1239 else
1240 {
1241 G4double distz= std::max(fZSections[0
1242 G4double dd = DistanceToPolygonSqr(p)
1243 if (distz > 0) dd += distz*distz;
1244 return std::sqrt(dd);
1245 }
1246 }
1247 }
1248
1249 // General case: use tessellated solid
1250 return G4TessellatedSolid::DistanceToIn(p);
1251 }
1252
1253 //___________________________________________
1254
1255 G4double G4ExtrudedSolid::DistanceToOut (cons
1256 cons
1257 cons
1258
1259
1260 {
1261 G4bool getnorm = calcNorm;
1262 if (getnorm) *validNorm = true;
1263
1264 G4double z0 = fZSections[0].fZ;
1265 G4double z1 = fZSections[fNz-1].fZ;
1266 if ((p.z() <= z0 + kCarToleranceHalf) && v.
1267 {
1268 if (getnorm) n->set(0,0,-1);
1269 return 0;
1270 }
1271 if ((p.z() >= z1 - kCarToleranceHalf) && v.
1272 {
1273 if (getnorm) n->set(0,0,1);
1274 return 0;
1275 }
1276
1277 switch (fSolidType)
1278 {
1279 case 1: // convex right prism
1280 {
1281 // Intersection with Z planes
1282 //
1283 G4double dz = (z1 - z0)*0.5;
1284 G4double pz = p.z() - 0.5 * (z0 + z1);
1285
1286 G4double vz = v.z();
1287 G4double tmax = (vz == 0) ? DBL_MAX : (
1288 G4int iside = (vz < 0) ? -4 : -2; // li
1289
1290 // Intersection with lateral planes
1291 //
1292 std::size_t np = fPlanes.size();
1293 for (std::size_t i=0; i<np; ++i)
1294 {
1295 G4double cosa = fPlanes[i].a*v.x()+fP
1296 if (cosa > 0)
1297 {
1298 G4double dist = fPlanes[i].a*p.x()+
1299 if (dist >= -kCarToleranceHalf)
1300 {
1301 if (getnorm) n->set(fPlanes[i].a,
1302 return 0;
1303 }
1304 G4double tmp = -dist/cosa;
1305 if (tmax > tmp) { tmax = tmp; iside
1306 }
1307 }
1308
1309 // Set normal, if required, and return
1310 //
1311 if (getnorm)
1312 {
1313 if (iside < 0)
1314 { n->set(0, 0, iside + 3); } // (-4
1315 else
1316 { n->set(fPlanes[iside].a, fPlanes[
1317 }
1318 return tmax;
1319 }
1320 case 2: // non-convex right prism
1321 {
1322 }
1323 }
1324
1325 // Override the base class function to rede
1326 // (the solid can be concave)
1327
1328 G4double distOut =
1329 G4TessellatedSolid::DistanceToOut(p, v, c
1330 if (validNorm != nullptr) { *validNorm = fI
1331
1332 return distOut;
1333 }
1334
1335 //___________________________________________
1336
1337 G4double G4ExtrudedSolid::DistanceToOut(const
1338 {
1339 switch (fSolidType)
1340 {
1341 case 1: // convex right prism
1342 {
1343 G4double dist = std::max(fZSections[0].
1344 std::size_t np = fPlanes.size();
1345 for (std::size_t i=0; i<np; ++i)
1346 {
1347 G4double dd = fPlanes[i].a*p.x() + fP
1348 if (dd > dist) dist = dd;
1349 }
1350 return (dist < 0) ? -dist : 0.;
1351 }
1352 case 2: // non-convex right prism
1353 {
1354 G4double distz = std::max(fZSections[0]
1355 G4bool in = PointInPolygon(p);
1356 if (distz >= 0 || (!in)) return 0; // p
1357 return std::min(-distz,std::sqrt(Distan
1358 }
1359 }
1360
1361 // General case: use tessellated solid
1362 return G4TessellatedSolid::DistanceToOut(p)
1363 }
1364
1365 //___________________________________________
1366 // Get bounding box
1367
1368 void G4ExtrudedSolid::BoundingLimits(G4ThreeV
1369 G4ThreeV
1370 {
1371 G4double xmin0 = kInfinity, xmax0 = -kInfin
1372 G4double ymin0 = kInfinity, ymax0 = -kInfin
1373
1374 for (G4int i=0; i<GetNofVertices(); ++i)
1375 {
1376 G4double x = fPolygon[i].x();
1377 if (x < xmin0) xmin0 = x;
1378 if (x > xmax0) xmax0 = x;
1379 G4double y = fPolygon[i].y();
1380 if (y < ymin0) ymin0 = y;
1381 if (y > ymax0) ymax0 = y;
1382 }
1383
1384 G4double xmin = kInfinity, xmax = -kInfinit
1385 G4double ymin = kInfinity, ymax = -kInfinit
1386
1387 G4int nsect = GetNofZSections();
1388 for (G4int i=0; i<nsect; ++i)
1389 {
1390 ZSection zsect = GetZSection(i);
1391 G4double dx = zsect.fOffset.x();
1392 G4double dy = zsect.fOffset.y();
1393 G4double scale = zsect.fScale;
1394 xmin = std::min(xmin,xmin0*scale+dx);
1395 xmax = std::max(xmax,xmax0*scale+dx);
1396 ymin = std::min(ymin,ymin0*scale+dy);
1397 ymax = std::max(ymax,ymax0*scale+dy);
1398 }
1399
1400 G4double zmin = GetZSection(0).fZ;
1401 G4double zmax = GetZSection(nsect-1).fZ;
1402
1403 pMin.set(xmin,ymin,zmin);
1404 pMax.set(xmax,ymax,zmax);
1405
1406 // Check correctness of the bounding box
1407 //
1408 if (pMin.x() >= pMax.x() || pMin.y() >= pMa
1409 {
1410 std::ostringstream message;
1411 message << "Bad bounding box (min >= max)
1412 << GetName() << " !"
1413 << "\npMin = " << pMin
1414 << "\npMax = " << pMax;
1415 G4Exception("G4ExtrudedSolid::BoundingLim
1416 "GeomMgt0001", JustWarning, m
1417 DumpInfo();
1418 }
1419 }
1420
1421 //___________________________________________
1422 // Calculate extent under transform and speci
1423
1424 G4bool
1425 G4ExtrudedSolid::CalculateExtent(const EAxis
1426 const G4Voxe
1427 const G4Affi
1428 G4doub
1429 {
1430 G4ThreeVector bmin, bmax;
1431 G4bool exist;
1432
1433 // Check bounding box (bbox)
1434 //
1435 BoundingLimits(bmin,bmax);
1436 G4BoundingEnvelope bbox(bmin,bmax);
1437 #ifdef G4BBOX_EXTENT
1438 return bbox.CalculateExtent(pAxis,pVoxelLim
1439 #endif
1440 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVo
1441 {
1442 return exist = pMin < pMax;
1443 }
1444
1445 // To find the extent, the base polygon is
1446 // The extent is calculated as cumulative e
1447 // formed by extrusion of the triangles
1448 //
1449 G4TwoVectorList triangles;
1450 G4double eminlim = pVoxelLimit.GetMinExtent
1451 G4double emaxlim = pVoxelLimit.GetMaxExtent
1452
1453 // triangulate the base polygon
1454 if (!G4GeomTools::TriangulatePolygon(fPolyg
1455 {
1456 std::ostringstream message;
1457 message << "Triangulation of the base pol
1458 << GetName() << " !"
1459 << "\nExtent has been calculated
1460 G4Exception("G4ExtrudedSolid::CalculateEx
1461 "GeomMgt1002",JustWarning,mes
1462 return bbox.CalculateExtent(pAxis,pVoxelL
1463 }
1464
1465 // allocate vector lists
1466 G4int nsect = GetNofZSections();
1467 std::vector<const G4ThreeVectorList *> poly
1468 polygons.resize(nsect);
1469 for (G4int k=0; k<nsect; ++k) { polygons[k]
1470
1471 // main loop along triangles
1472 pMin = kInfinity;
1473 pMax = -kInfinity;
1474 G4int ntria = (G4int)triangles.size()/3;
1475 for (G4int i=0; i<ntria; ++i)
1476 {
1477 G4int i3 = i*3;
1478 for (G4int k=0; k<nsect; ++k) // extrude
1479 {
1480 ZSection zsect = GetZSection(k);
1481 G4double z = zsect.fZ;
1482 G4double dx = zsect.fOffset.x();
1483 G4double dy = zsect.fOffset.y();
1484 G4double scale = zsect.fScale;
1485
1486 auto ptr = const_cast<G4ThreeVectorList
1487 auto iter = ptr->begin();
1488 G4double x0 = triangles[i3+0].x()*scale
1489 G4double y0 = triangles[i3+0].y()*scale
1490 iter->set(x0,y0,z);
1491 iter++;
1492 G4double x1 = triangles[i3+1].x()*scale
1493 G4double y1 = triangles[i3+1].y()*scale
1494 iter->set(x1,y1,z);
1495 iter++;
1496 G4double x2 = triangles[i3+2].x()*scale
1497 G4double y2 = triangles[i3+2].y()*scale
1498 iter->set(x2,y2,z);
1499 }
1500
1501 // set sub-envelope and adjust extent
1502 G4double emin,emax;
1503 G4BoundingEnvelope benv(polygons);
1504 if (!benv.CalculateExtent(pAxis,pVoxelLim
1505 if (emin < pMin) pMin = emin;
1506 if (emax > pMax) pMax = emax;
1507 if (eminlim > pMin && emaxlim < pMax) bre
1508 }
1509 // free memory
1510 for (G4int k=0; k<nsect; ++k) { delete poly
1511 return (pMin < pMax);
1512 }
1513
1514 //___________________________________________
1515
1516 std::ostream& G4ExtrudedSolid::StreamInfo(std
1517 {
1518 G4long oldprc = os.precision(16);
1519 os << "------------------------------------
1520 << " *** Dump for solid - " << GetNam
1521 << " ================================
1522 << " Solid geometry type: " << fGeometry
1523
1524 if ( fIsConvex)
1525 { os << " Convex polygon; list of vertice
1526 else
1527 { os << " Concave polygon; list of vertic
1528
1529 for ( std::size_t i=0; i<fNv; ++i )
1530 {
1531 os << std::setw(5) << "#" << i
1532 << " vx = " << fPolygon[i].x()/mm <<
1533 << " vy = " << fPolygon[i].y()/mm <<
1534 }
1535
1536 os << " Sections:" << G4endl;
1537 for ( std::size_t iz=0; iz<fNz; ++iz )
1538 {
1539 os << " z = " << fZSections[iz].fZ/mm
1540 << " x0= " << fZSections[iz].fOffs
1541 << " y0= " << fZSections[iz].fOffs
1542 << " scale= " << fZSections[iz].fScal
1543 }
1544
1545 /*
1546 // Triangles (for debugging)
1547 os << G4endl;
1548 os << " Triangles:" << G4endl;
1549 os << " Triangle # vertex1 vertex2 ve
1550
1551 G4int counter = 0;
1552 std::vector< std::vector<G4int> >::const_it
1553 for ( it = fTriangles.begin(); it != fTrian
1554 std::vector<G4int> triangle = *it;
1555 os << std::setw(10) << counter++
1556 << std::setw(10) << triangle[0] << st
1557 << std::setw(10) << triangle[2]
1558 << G4endl;
1559 }
1560 */
1561 os.precision(oldprc);
1562
1563 return os;
1564 }
1565
1566 #endif
1567