| 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 // Implementation of G4Polycone, a CSG polycon
27 //
28 // Author: David C. Williams (davidw@scipp.ucs
29 // -------------------------------------------
30
31 #include "G4Polycone.hh"
32
33 #if !defined(G4GEOM_USE_UPOLYCONE)
34
35 #include "G4PolyconeSide.hh"
36 #include "G4PolyPhiFace.hh"
37
38 #include "G4GeomTools.hh"
39 #include "G4VoxelLimits.hh"
40 #include "G4AffineTransform.hh"
41 #include "G4BoundingEnvelope.hh"
42
43 #include "G4QuickRand.hh"
44
45 #include "G4EnclosingCylinder.hh"
46 #include "G4ReduciblePolygon.hh"
47 #include "G4VPVParameterisation.hh"
48
49 namespace
50 {
51 G4Mutex surface_elementsMutex = G4MUTEX_INIT
52 }
53
54 using namespace CLHEP;
55
56 // Constructor (GEANT3 style parameters)
57 //
58 G4Polycone::G4Polycone( const G4String& name,
59 G4double phiStar
60 G4double phiTota
61 G4int numZPlanes
62 const G4double zPlane[
63 const G4double rInner[
64 const G4double rOuter[
65 : G4VCSGfaceted( name )
66 {
67 //
68 // Some historical ugliness
69 //
70 original_parameters = new G4PolyconeHistoric
71
72 original_parameters->Start_angle = phiStart;
73 original_parameters->Opening_angle = phiTota
74 original_parameters->Num_z_planes = numZPlan
75 original_parameters->Z_values = new G4double
76 original_parameters->Rmin = new G4double[num
77 original_parameters->Rmax = new G4double[num
78
79 for (G4int i=0; i<numZPlanes; ++i)
80 {
81 if(rInner[i]>rOuter[i])
82 {
83 DumpInfo();
84 std::ostringstream message;
85 message << "Cannot create a Polycone wit
86 << G4endl
87 << " rInner > rOuter for
88 << " rMin[" << i << "] =
89 << " -- rMax[" << i << "] = " <<
90 G4Exception("G4Polycone::G4Polycone()",
91 FatalErrorInArgument, messag
92 }
93 if (( i < numZPlanes-1) && ( zPlane[i] ==
94 {
95 if( (rInner[i] > rOuter[i+1])
96 ||(rInner[i+1] > rOuter[i]) )
97 {
98 DumpInfo();
99 std::ostringstream message;
100 message << "Cannot create a Polycone w
101 << G4endl
102 << " Segments are not c
103 << " rMin[" << i << "]
104 << " -- rMax[" << i+1 << "] =
105 << " rMin[" << i+1 << "
106 << " -- rMax[" << i << "] = "
107 G4Exception("G4Polycone::G4Polycone()"
108 FatalErrorInArgument, mess
109 }
110 }
111 original_parameters->Z_values[i] = zPlane[
112 original_parameters->Rmin[i] = rInner[i];
113 original_parameters->Rmax[i] = rOuter[i];
114 }
115
116 //
117 // Build RZ polygon using special PCON/PGON
118 //
119 auto rz = new G4ReduciblePolygon( rInner, rO
120
121 //
122 // Do the real work
123 //
124 Create( phiStart, phiTotal, rz );
125
126 delete rz;
127 }
128
129 // Constructor (generic parameters)
130 //
131 G4Polycone::G4Polycone( const G4String& name,
132 G4double phiStar
133 G4double phiTota
134 G4int numRZ,
135 const G4double r[],
136 const G4double z[] )
137 : G4VCSGfaceted( name )
138 {
139
140 auto rz = new G4ReduciblePolygon( r, z, numR
141
142 Create( phiStart, phiTotal, rz );
143
144 // Set original_parameters struct for consis
145 //
146
147 G4bool convertible = SetOriginalParameters(r
148
149 if(!convertible)
150 {
151 std::ostringstream message;
152 message << "Polycone " << GetName() << "ca
153 << "to Polycone with (Rmin,Rmaz,Z)
154 G4Exception("G4Polycone::G4Polycone()", "G
155 FatalException, message, "Use
156 }
157 else
158 {
159 G4cout << "INFO: Converting polycone " <<
160 << "to optimized polycone with (Rmi
161 << G4endl;
162 }
163 delete rz;
164 }
165
166 // Create
167 //
168 // Generic create routine, called by each cons
169 // conversion of arguments
170 //
171 void G4Polycone::Create( G4double phiStart,
172 G4double phiTotal,
173 G4ReduciblePolygon* r
174 {
175 //
176 // Perform checks of rz values
177 //
178 if (rz->Amin() < 0.0)
179 {
180 std::ostringstream message;
181 message << "Illegal input parameters - " <
182 << " All R values must be >
183 G4Exception("G4Polycone::Create()", "GeomS
184 FatalErrorInArgument, message)
185 }
186
187 G4double rzArea = rz->Area();
188 if (rzArea < -kCarTolerance)
189 {
190 rz->ReverseOrder();
191 }
192 else if (rzArea < kCarTolerance)
193 {
194 std::ostringstream message;
195 message << "Illegal input parameters - " <
196 << " R/Z cross section is z
197 G4Exception("G4Polycone::Create()", "GeomS
198 FatalErrorInArgument, message)
199 }
200
201 if ( (!rz->RemoveDuplicateVertices( kCarTole
202 || (!rz->RemoveRedundantVertices( kCarTole
203 {
204 std::ostringstream message;
205 message << "Illegal input parameters - " <
206 << " Too few unique R/Z val
207 G4Exception("G4Polycone::Create()", "GeomS
208 FatalErrorInArgument, message)
209 }
210
211 if (rz->CrossesItself(1/kInfinity))
212 {
213 std::ostringstream message;
214 message << "Illegal input parameters - " <
215 << " R/Z segments cross !";
216 G4Exception("G4Polycone::Create()", "GeomS
217 FatalErrorInArgument, message)
218 }
219
220 numCorner = rz->NumVertices();
221
222 startPhi = phiStart;
223 while( startPhi < 0. ) // Loop checking,
224 startPhi += twopi;
225 //
226 // Phi opening? Account for some possible ro
227 // nonsense value as representing no phi ope
228 //
229 if ( (phiTotal <= 0) || (phiTotal > twopi*(1
230 {
231 phiIsOpen = false;
232 startPhi = 0.;
233 endPhi = twopi;
234 }
235 else
236 {
237 phiIsOpen = true;
238 endPhi = startPhi + phiTotal;
239 }
240
241 //
242 // Allocate corner array.
243 //
244 corners = new G4PolyconeSideRZ[numCorner];
245
246 //
247 // Copy corners
248 //
249 G4ReduciblePolygonIterator iterRZ(rz);
250
251 G4PolyconeSideRZ *next = corners;
252 iterRZ.Begin();
253 do // Loop checking, 13.08.2015, G.Cosmo
254 {
255 next->r = iterRZ.GetA();
256 next->z = iterRZ.GetB();
257 } while( ++next, iterRZ.Next() );
258
259 //
260 // Allocate face pointer array
261 //
262 numFace = phiIsOpen ? numCorner+2 : numCorne
263 faces = new G4VCSGface*[numFace];
264
265 //
266 // Construct conical faces
267 //
268 // But! Don't construct a face if both point
269 //
270 G4PolyconeSideRZ* corner = corners,
271 * prev = corners + numCorner
272 * nextNext;
273 G4VCSGface **face = faces;
274 do // Loop checking, 13.08.2015, G.Cosmo
275 {
276 next = corner+1;
277 if (next >= corners+numCorner) next = corn
278 nextNext = next+1;
279 if (nextNext >= corners+numCorner) nextNex
280
281 if (corner->r < 1/kInfinity && next->r < 1
282
283 //
284 // We must decide here if we can dare decl
285 // as having a "valid" normal (i.e. allBeh
286 // is never possible if the face faces "in
287 //
288 G4bool allBehind;
289 if (corner->z > next->z)
290 {
291 allBehind = false;
292 }
293 else
294 {
295 //
296 // Otherwise, it is only true if the lin
297 // through the two points of the segment
298 // split the r/z cross section
299 //
300 allBehind = !rz->BisectedBy( corner->r,
301 next->r, next->z, kCarToleran
302 }
303
304 *face++ = new G4PolyconeSide( prev, corner
305 startPhi, endPhi-startPhi, phi
306 } while( prev=corner, corner=next, corner >
307
308 if (phiIsOpen)
309 {
310 //
311 // Construct phi open edges
312 //
313 *face++ = new G4PolyPhiFace( rz, startPhi,
314 *face++ = new G4PolyPhiFace( rz, endPhi,
315 }
316
317 //
318 // We might have dropped a face or two: reca
319 //
320 numFace = (G4int)(face-faces);
321
322 //
323 // Make enclosingCylinder
324 //
325 enclosingCylinder =
326 new G4EnclosingCylinder( rz, phiIsOpen, ph
327 }
328
329 // Fake default constructor - sets only member
330 // for usage restri
331 //
332 G4Polycone::G4Polycone( __void__& a )
333 : G4VCSGfaceted(a), startPhi(0.), endPhi(0.
334 {
335 }
336
337
338 //
339 // Destructor
340 //
341 G4Polycone::~G4Polycone()
342 {
343 delete [] corners;
344 delete original_parameters;
345 delete enclosingCylinder;
346 delete fElements;
347 delete fpPolyhedron;
348 corners = nullptr;
349 original_parameters = nullptr;
350 enclosingCylinder = nullptr;
351 fElements = nullptr;
352 fpPolyhedron = nullptr;
353 }
354
355 // Copy constructor
356 //
357 G4Polycone::G4Polycone( const G4Polycone& sour
358 : G4VCSGfaceted( source )
359 {
360 CopyStuff( source );
361 }
362
363 // Assignment operator
364 //
365 G4Polycone &G4Polycone::operator=( const G4Pol
366 {
367 if (this == &source) return *this;
368
369 G4VCSGfaceted::operator=( source );
370
371 delete [] corners;
372 delete original_parameters;
373
374 delete enclosingCylinder;
375
376 CopyStuff( source );
377
378 return *this;
379 }
380
381 // CopyStuff
382 //
383 void G4Polycone::CopyStuff( const G4Polycone&
384 {
385 //
386 // Simple stuff
387 //
388 startPhi = source.startPhi;
389 endPhi = source.endPhi;
390 phiIsOpen = source.phiIsOpen;
391 numCorner = source.numCorner;
392
393 //
394 // The corner array
395 //
396 corners = new G4PolyconeSideRZ[numCorner];
397
398 G4PolyconeSideRZ* corn = corners,
399 * sourceCorn = source.corner
400 do // Loop checking, 13.08.2015, G.Cosmo
401 {
402 *corn = *sourceCorn;
403 } while( ++sourceCorn, ++corn < corners+numC
404
405 //
406 // Original parameters
407 //
408 if (source.original_parameters != nullptr)
409 {
410 original_parameters =
411 new G4PolyconeHistorical( *source.origin
412 }
413
414 //
415 // Enclosing cylinder
416 //
417 enclosingCylinder = new G4EnclosingCylinder(
418
419 //
420 // Surface elements
421 //
422 delete fElements;
423 fElements = nullptr;
424
425 //
426 // Polyhedron
427 //
428 fRebuildPolyhedron = false;
429 delete fpPolyhedron;
430 fpPolyhedron = nullptr;
431 }
432
433 // Reset
434 //
435 G4bool G4Polycone::Reset()
436 {
437 //
438 // Clear old setup
439 //
440 G4VCSGfaceted::DeleteStuff();
441 delete [] corners;
442 delete enclosingCylinder;
443 delete fElements;
444 corners = nullptr;
445 fElements = nullptr;
446 enclosingCylinder = nullptr;
447
448 //
449 // Rebuild polycone
450 //
451 auto rz = new G4ReduciblePolygon( original_p
452 original_p
453 original_p
454 original_p
455 Create( original_parameters->Start_angle,
456 original_parameters->Opening_angle,
457 delete rz;
458
459 return false;
460 }
461
462 // Inside
463 //
464 // This is an override of G4VCSGfaceted::Insid
465 // to speed things up by first checking with G
466 //
467 EInside G4Polycone::Inside( const G4ThreeVecto
468 {
469 //
470 // Quick test
471 //
472 if (enclosingCylinder->MustBeOutside(p)) ret
473
474 //
475 // Long answer
476 //
477 return G4VCSGfaceted::Inside(p);
478 }
479
480 // DistanceToIn
481 //
482 // This is an override of G4VCSGfaceted::Insid
483 // to speed things up by first checking with G
484 //
485 G4double G4Polycone::DistanceToIn( const G4Thr
486 const G4Thr
487 {
488 //
489 // Quick test
490 //
491 if (enclosingCylinder->ShouldMiss(p,v))
492 return kInfinity;
493
494 //
495 // Long answer
496 //
497 return G4VCSGfaceted::DistanceToIn( p, v );
498 }
499
500 // DistanceToIn
501 //
502 G4double G4Polycone::DistanceToIn( const G4Thr
503 {
504 return G4VCSGfaceted::DistanceToIn(p);
505 }
506
507 // Get bounding box
508 //
509 void G4Polycone::BoundingLimits(G4ThreeVector&
510 G4ThreeVector&
511 {
512 G4double rmin = kInfinity, rmax = -kInfinity
513 G4double zmin = kInfinity, zmax = -kInfinity
514
515 for (G4int i=0; i<GetNumRZCorner(); ++i)
516 {
517 G4PolyconeSideRZ corner = GetCorner(i);
518 if (corner.r < rmin) rmin = corner.r;
519 if (corner.r > rmax) rmax = corner.r;
520 if (corner.z < zmin) zmin = corner.z;
521 if (corner.z > zmax) zmax = corner.z;
522 }
523
524 if (IsOpen())
525 {
526 G4TwoVector vmin,vmax;
527 G4GeomTools::DiskExtent(rmin,rmax,
528 GetSinStartPhi(),G
529 GetSinEndPhi(),Get
530 vmin,vmax);
531 pMin.set(vmin.x(),vmin.y(),zmin);
532 pMax.set(vmax.x(),vmax.y(),zmax);
533 }
534 else
535 {
536 pMin.set(-rmax,-rmax, zmin);
537 pMax.set( rmax, rmax, zmax);
538 }
539
540 // Check correctness of the bounding box
541 //
542 if (pMin.x() >= pMax.x() || pMin.y() >= pMax
543 {
544 std::ostringstream message;
545 message << "Bad bounding box (min >= max)
546 << GetName() << " !"
547 << "\npMin = " << pMin
548 << "\npMax = " << pMax;
549 G4Exception("G4Polycone::BoundingLimits()"
550 JustWarning, message);
551 DumpInfo();
552 }
553 }
554
555 // Calculate extent under transform and specif
556 //
557 G4bool G4Polycone::CalculateExtent(const EAxis
558 const G4Vox
559 const G4Aff
560 G4dou
561 {
562 G4ThreeVector bmin, bmax;
563 G4bool exist;
564
565 // Check bounding box (bbox)
566 //
567 BoundingLimits(bmin,bmax);
568 G4BoundingEnvelope bbox(bmin,bmax);
569 #ifdef G4BBOX_EXTENT
570 return bbox.CalculateExtent(pAxis,pVoxelLimi
571 #endif
572 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox
573 {
574 return exist = pMin < pMax;
575 }
576
577 // To find the extent, RZ contour of the pol
578 // in triangles. The extent is calculated as
579 // all sub-polycones formed by rotation of t
580 //
581 G4TwoVectorList contourRZ;
582 G4TwoVectorList triangles;
583 std::vector<G4int> iout;
584 G4double eminlim = pVoxelLimit.GetMinExtent(
585 G4double emaxlim = pVoxelLimit.GetMaxExtent(
586
587 // get RZ contour, ensure anticlockwise orde
588 for (G4int i=0; i<GetNumRZCorner(); ++i)
589 {
590 G4PolyconeSideRZ corner = GetCorner(i);
591 contourRZ.emplace_back(corner.r,corner.z);
592 }
593 G4GeomTools::RemoveRedundantVertices(contour
594 G4double area = G4GeomTools::PolygonArea(con
595 if (area < 0.) std::reverse(contourRZ.begin(
596
597 // triangulate RZ countour
598 if (!G4GeomTools::TriangulatePolygon(contour
599 {
600 std::ostringstream message;
601 message << "Triangulation of RZ contour ha
602 << GetName() << " !"
603 << "\nExtent has been calculated u
604 G4Exception("G4Polycone::CalculateExtent()
605 "GeomMgt1002", JustWarning, me
606 return bbox.CalculateExtent(pAxis,pVoxelLi
607 }
608
609 // set trigonometric values
610 const G4int NSTEPS = 24; // numbe
611 G4double astep = twopi/NSTEPS; // max a
612
613 G4double sphi = GetStartPhi();
614 G4double ephi = GetEndPhi();
615 G4double dphi = IsOpen() ? ephi-sphi : two
616 G4int ksteps = (dphi <= astep) ? 1 : (G4i
617 G4double ang = dphi/ksteps;
618
619 G4double sinHalf = std::sin(0.5*ang);
620 G4double cosHalf = std::cos(0.5*ang);
621 G4double sinStep = 2.*sinHalf*cosHalf;
622 G4double cosStep = 1. - 2.*sinHalf*sinHalf;
623
624 G4double sinStart = GetSinStartPhi();
625 G4double cosStart = GetCosStartPhi();
626 G4double sinEnd = GetSinEndPhi();
627 G4double cosEnd = GetCosEndPhi();
628
629 // define vectors and arrays
630 std::vector<const G4ThreeVectorList *> polyg
631 polygons.resize(ksteps+2);
632 G4ThreeVectorList pols[NSTEPS+2];
633 for (G4int k=0; k<ksteps+2; ++k) pols[k].res
634 for (G4int k=0; k<ksteps+2; ++k) polygons[k]
635 G4double r0[6],z0[6]; // contour with origin
636 G4double r1[6]; // shifted radii of ex
637
638 // main loop along triangles
639 pMin = kInfinity;
640 pMax =-kInfinity;
641 G4int ntria = (G4int)triangles.size()/3;
642 for (G4int i=0; i<ntria; ++i)
643 {
644 G4int i3 = i*3;
645 for (G4int k=0; k<3; ++k)
646 {
647 G4int e0 = i3+k, e1 = (k<2) ? e0+1 : i3;
648 G4int k2 = k*2;
649 // set contour with original edges of tr
650 r0[k2+0] = triangles[e0].x(); z0[k2+0] =
651 r0[k2+1] = triangles[e1].x(); z0[k2+1] =
652 // set shifted radii
653 r1[k2+0] = r0[k2+0];
654 r1[k2+1] = r0[k2+1];
655 if (z0[k2+1] - z0[k2+0] <= 0) continue;
656 r1[k2+0] /= cosHalf;
657 r1[k2+1] /= cosHalf;
658 }
659
660 // rotate countour, set sequence of 6-side
661 G4double sinCur = sinStart*cosHalf + cosSt
662 G4double cosCur = cosStart*cosHalf - sinSt
663 for (G4int j=0; j<6; ++j) pols[0][j].set(r
664 for (G4int k=1; k<ksteps+1; ++k)
665 {
666 for (G4int j=0; j<6; ++j) pols[k][j].set
667 G4double sinTmp = sinCur;
668 sinCur = sinCur*cosStep + cosCur*sinStep
669 cosCur = cosCur*cosStep - sinTmp*sinStep
670 }
671 for (G4int j=0; j<6; ++j) pols[ksteps+1][j
672
673 // set sub-envelope and adjust extent
674 G4double emin,emax;
675 G4BoundingEnvelope benv(polygons);
676 if (!benv.CalculateExtent(pAxis,pVoxelLimi
677 if (emin < pMin) pMin = emin;
678 if (emax > pMax) pMax = emax;
679 if (eminlim > pMin && emaxlim < pMax) retu
680 }
681 return (pMin < pMax);
682 }
683
684 // ComputeDimensions
685 //
686 void G4Polycone::ComputeDimensions( G4VP
687 const G4in
688 const G4VP
689 {
690 p->ComputeDimensions(*this,n,pRep);
691 }
692
693 // GetEntityType
694 //
695 G4GeometryType G4Polycone::GetEntityType() co
696 {
697 return {"G4Polycone"};
698 }
699
700 // Make a clone of the object
701 //
702 G4VSolid* G4Polycone::Clone() const
703 {
704 return new G4Polycone(*this);
705 }
706
707 //
708 // Stream object contents to an output stream
709 //
710 std::ostream& G4Polycone::StreamInfo( std::ost
711 {
712 G4long oldprc = os.precision(16);
713 os << "-------------------------------------
714 << " *** Dump for solid - " << GetName
715 << " =================================
716 << " Solid type: G4Polycone\n"
717 << " Parameters: \n"
718 << " starting phi angle : " << startPh
719 << " ending phi angle : " << endPhi/
720 G4int i=0;
721
722 G4int numPlanes = original_parameters->Num
723 os << " number of Z planes: " << numPla
724 << " Z values: \n";
725 for (i=0; i<numPlanes; ++i)
726 {
727 os << " Z plane " << i << "
728 << original_parameters->Z_values[i] <
729 }
730 os << " Tangent distances to
731 for (i=0; i<numPlanes; ++i)
732 {
733 os << " Z plane " << i << "
734 << original_parameters->Rmin[i] << "\
735 }
736 os << " Tangent distances to
737 for (i=0; i<numPlanes; ++i)
738 {
739 os << " Z plane " << i << "
740 << original_parameters->Rmax[i] << "\
741 }
742
743 os << " number of RZ points: " << numCorn
744 << " RZ values (corners): \n
745 for (i=0; i<numCorner; ++i)
746 {
747 os << " "
748 << corners[i].r << ", " << corners[i
749 }
750 os << "-------------------------------------
751 os.precision(oldprc);
752
753 return os;
754 }
755
756 //////////////////////////////////////////////
757 //
758 // Return volume
759
760 G4double G4Polycone::GetCubicVolume()
761 {
762 if (fCubicVolume == 0.)
763 {
764 G4double total = 0.;
765 G4int nrz = GetNumRZCorner();
766 G4PolyconeSideRZ a = GetCorner(nrz - 1);
767 for (G4int i=0; i<nrz; ++i)
768 {
769 G4PolyconeSideRZ b = GetCorner(i);
770 total += (b.r*b.r + b.r*a.r + a.r*a.r)*(
771 a = b;
772 }
773 fCubicVolume = std::abs(total)*(GetEndPhi(
774 }
775 return fCubicVolume;
776 }
777
778 //////////////////////////////////////////////
779 //
780 // Return surface area
781
782 G4double G4Polycone::GetSurfaceArea()
783 {
784 if (fSurfaceArea == 0.)
785 {
786 // phi cut area
787 G4int nrz = GetNumRZCorner();
788 G4double scut = 0.;
789 if (IsOpen())
790 {
791 G4PolyconeSideRZ a = GetCorner(nrz - 1);
792 for (G4int i=0; i<nrz; ++i)
793 {
794 G4PolyconeSideRZ b = GetCorner(i);
795 scut += a.r*b.z - a.z*b.r;
796 a = b;
797 }
798 scut = std::abs(scut);
799 }
800 // lateral surface area
801 G4double slat = 0;
802 G4PolyconeSideRZ a = GetCorner(nrz - 1);
803 for (G4int i=0; i<nrz; ++i)
804 {
805 G4PolyconeSideRZ b = GetCorner(i);
806 G4double h = std::sqrt((b.r - a.r)*(b.r
807 slat += (b.r + a.r)*h;
808 a = b;
809 }
810 slat *= (GetEndPhi() - GetStartPhi())/2.;
811 fSurfaceArea = scut + slat;
812 }
813 return fSurfaceArea;
814 }
815
816 //////////////////////////////////////////////
817 //
818 // Set vector of surface elements, auxiliary m
819 // random points on surface
820
821 void G4Polycone::SetSurfaceElements() const
822 {
823 fElements = new std::vector<G4Polycone::surf
824 G4double total = 0.;
825 G4int nrz = GetNumRZCorner();
826
827 // set lateral surface elements
828 G4double dphi = GetEndPhi() - GetStartPhi();
829 G4int ia = nrz - 1;
830 for (G4int ib=0; ib<nrz; ++ib)
831 {
832 G4PolyconeSideRZ a = GetCorner(ia);
833 G4PolyconeSideRZ b = GetCorner(ib);
834 G4Polycone::surface_element selem;
835 selem.i0 = ia;
836 selem.i1 = ib;
837 selem.i2 = -1;
838 ia = ib;
839 if (a.r == 0. && b.r == 0.) continue;
840 G4double h = std::sqrt((b.r - a.r)*(b.r -
841 total += 0.5*dphi*(b.r + a.r)*h;
842 selem.area = total;
843 fElements->push_back(selem);
844 }
845
846 // set elements for phi cuts
847 if (IsOpen())
848 {
849 G4TwoVectorList contourRZ;
850 std::vector<G4int> triangles;
851 for (G4int i=0; i<nrz; ++i)
852 {
853 G4PolyconeSideRZ corner = GetCorner(i);
854 contourRZ.emplace_back(corner.r, corner.
855 }
856 G4GeomTools::TriangulatePolygon(contourRZ,
857 auto ntria = (G4int)triangles.size();
858 for (G4int i=0; i<ntria; i+=3)
859 {
860 G4Polycone::surface_element selem;
861 selem.i0 = triangles[i];
862 selem.i1 = triangles[i+1];
863 selem.i2 = triangles[i+2];
864 G4PolyconeSideRZ a = GetCorner(selem.i0)
865 G4PolyconeSideRZ b = GetCorner(selem.i1)
866 G4PolyconeSideRZ c = GetCorner(selem.i2)
867 G4double stria =
868 std::abs(G4GeomTools::TriangleArea(a.r
869 total += stria;
870 selem.area = total;
871 fElements->push_back(selem); // start ph
872 total += stria;
873 selem.area = total;
874 selem.i0 += nrz;
875 fElements->push_back(selem); // end phi
876 }
877 }
878 }
879
880 //////////////////////////////////////////////
881 //
882 // Generate random point on surface
883
884 G4ThreeVector G4Polycone::GetPointOnSurface()
885 {
886 // Set surface elements
887 if (fElements == nullptr)
888 {
889 G4AutoLock l(&surface_elementsMutex);
890 SetSurfaceElements();
891 l.unlock();
892 }
893
894 // Select surface element
895 G4Polycone::surface_element selem;
896 selem = fElements->back();
897 G4double select = selem.area*G4QuickRand();
898 auto it = std::lower_bound(fElements->begin(
899 [](const G4Polyco
900 -> G4bool { retur
901
902 // Generate random point
903 G4double r = 0, z = 0, phi = 0;
904 G4double u = G4QuickRand();
905 G4double v = G4QuickRand();
906 G4int i0 = (*it).i0;
907 G4int i1 = (*it).i1;
908 G4int i2 = (*it).i2;
909 if (i2 < 0) // lateral surface
910 {
911 G4PolyconeSideRZ p0 = GetCorner(i0);
912 G4PolyconeSideRZ p1 = GetCorner(i1);
913 if (p1.r < p0.r)
914 {
915 p0 = GetCorner(i1);
916 p1 = GetCorner(i0);
917 }
918 if (p1.r - p0.r < kCarTolerance) // cylind
919 {
920 r = (p1.r - p0.r)*u + p0.r;
921 z = (p1.z - p0.z)*u + p0.z;
922 }
923 else // conical surface
924 {
925 r = std::sqrt(p1.r*p1.r*u + p0.r*p0.r*(1
926 z = p0.z + (p1.z - p0.z)*(r - p0.r)/(p1.
927 }
928 phi = (GetEndPhi() - GetStartPhi())*v + Ge
929 }
930 else // phi cut
931 {
932 G4int nrz = GetNumRZCorner();
933 phi = (i0 < nrz) ? GetStartPhi() : GetEndP
934 if (i0 >= nrz) { i0 -= nrz; }
935 G4PolyconeSideRZ p0 = GetCorner(i0);
936 G4PolyconeSideRZ p1 = GetCorner(i1);
937 G4PolyconeSideRZ p2 = GetCorner(i2);
938 if (u + v > 1.) { u = 1. - u; v = 1. - v;
939 r = (p1.r - p0.r)*u + (p2.r - p0.r)*v + p
940 z = (p1.z - p0.z)*u + (p2.z - p0.z)*v + p
941 }
942 return { r*std::cos(phi), r*std::sin(phi), z
943 }
944
945 //////////////////////////////////////////////
946 //
947 // CreatePolyhedron
948
949 G4Polyhedron* G4Polycone::CreatePolyhedron() c
950 {
951 std::vector<G4TwoVector> rz(numCorner);
952 for (G4int i = 0; i < numCorner; ++i)
953 rz[i].set(corners[i].r, corners[i].z);
954 return new G4PolyhedronPcon(startPhi, endPhi
955 }
956
957 // SetOriginalParameters
958 //
959 G4bool G4Polycone::SetOriginalParameters(G4Re
960 {
961 G4int numPlanes = numCorner;
962 G4bool isConvertible = true;
963 G4double Zmax=rz->Bmax();
964 rz->StartWithZMin();
965
966 // Prepare vectors for storage
967 //
968 std::vector<G4double> Z;
969 std::vector<G4double> Rmin;
970 std::vector<G4double> Rmax;
971
972 G4int countPlanes=1;
973 G4int icurr=0;
974 G4int icurl=0;
975
976 // first plane Z=Z[0]
977 //
978 Z.push_back(corners[0].z);
979 G4double Zprev=Z[0];
980 if (Zprev == corners[1].z)
981 {
982 Rmin.push_back(corners[0].r);
983 Rmax.push_back (corners[1].r);icurr=1;
984 }
985 else if (Zprev == corners[numPlanes-1].z)
986 {
987 Rmin.push_back(corners[numPlanes-1].r);
988 Rmax.push_back (corners[0].r);
989 icurl=numPlanes-1;
990 }
991 else
992 {
993 Rmin.push_back(corners[0].r);
994 Rmax.push_back (corners[0].r);
995 }
996
997 // next planes until last
998 //
999 G4int inextr=0, inextl=0;
1000 for (G4int i=0; i < numPlanes-2; ++i)
1001 {
1002 inextr=1+icurr;
1003 inextl=(icurl <= 0)? numPlanes-1 : icurl-
1004
1005 if((corners[inextr].z >= Zmax) & (corners
1006
1007 G4double Zleft = corners[inextl].z;
1008 G4double Zright = corners[inextr].z;
1009 if(Zright > Zleft) // Next plane will be
1010 {
1011 Z.push_back(Zleft);
1012 countPlanes++;
1013 G4double difZr=corners[inextr].z - corn
1014 G4double difZl=corners[inextl].z - corn
1015
1016 if(std::fabs(difZl) < kCarTolerance)
1017 {
1018 if(std::fabs(difZr) < kCarTolerance)
1019 {
1020 Rmin.push_back(corners[inextl].r);
1021 Rmax.push_back(corners[icurr].r);
1022 }
1023 else
1024 {
1025 Rmin.push_back(corners[inextl].r);
1026 Rmax.push_back(corners[icurr].r + (
1027 *(corners[ine
1028 }
1029 }
1030 else if (difZl >= kCarTolerance)
1031 {
1032 if(std::fabs(difZr) < kCarTolerance)
1033 {
1034 Rmin.push_back(corners[icurl].r);
1035 Rmax.push_back(corners[icurr].r);
1036 }
1037 else
1038 {
1039 Rmin.push_back(corners[icurl].r);
1040 Rmax.push_back(corners[icurr].r + (
1041 *(corners[ine
1042 }
1043 }
1044 else
1045 {
1046 isConvertible=false; break;
1047 }
1048 icurl=(icurl == 0)? numPlanes-1 : icurl
1049 }
1050 else if(std::fabs(Zright-Zleft)<kCarToler
1051 {
1052 Z.push_back(Zleft);
1053 ++countPlanes;
1054 ++icurr;
1055
1056 icurl=(icurl == 0)? numPlanes-1 : icurl
1057
1058 Rmin.push_back(corners[inextl].r);
1059 Rmax.push_back(corners[inextr].r);
1060 }
1061 else // Zright<Zleft
1062 {
1063 Z.push_back(Zright);
1064 ++countPlanes;
1065
1066 G4double difZr=corners[inextr].z - corn
1067 G4double difZl=corners[inextl].z - corn
1068 if(std::fabs(difZr) < kCarTolerance)
1069 {
1070 if(std::fabs(difZl) < kCarTolerance)
1071 {
1072 Rmax.push_back(corners[inextr].r);
1073 Rmin.push_back(corners[icurr].r);
1074 }
1075 else
1076 {
1077 Rmin.push_back(corners[icurl].r + (
1078 *(corners[ine
1079 Rmax.push_back(corners[inextr].r);
1080 }
1081 ++icurr;
1082 } // plate
1083 else if (difZr >= kCarTolerance)
1084 {
1085 if(std::fabs(difZl) < kCarTolerance)
1086 {
1087 Rmax.push_back(corners[inextr].r);
1088 Rmin.push_back (corners[icurr].r);
1089 }
1090 else
1091 {
1092 Rmax.push_back(corners[inextr].r);
1093 Rmin.push_back (corners[icurl].r+(Z
1094 * (corners[
1095 }
1096 ++icurr;
1097 }
1098 else
1099 {
1100 isConvertible=false; break;
1101 }
1102 }
1103 } // end for loop
1104
1105 // last plane Z=Zmax
1106 //
1107 Z.push_back(Zmax);
1108 ++countPlanes;
1109 inextr=1+icurr;
1110 inextl=(icurl <= 0)? numPlanes-1 : icurl-1;
1111
1112 Rmax.push_back(corners[inextr].r);
1113 Rmin.push_back(corners[inextl].r);
1114
1115 // Set original parameters Rmin,Rmax,Z
1116 //
1117 if(isConvertible)
1118 {
1119 original_parameters = new G4PolyconeHistor
1120 original_parameters->Z_values = new G4doub
1121 original_parameters->Rmin = new G4double[c
1122 original_parameters->Rmax = new G4double[c
1123
1124 for(G4int j=0; j < countPlanes; ++j)
1125 {
1126 original_parameters->Z_values[j] = Z[j];
1127 original_parameters->Rmax[j] = Rmax[j];
1128 original_parameters->Rmin[j] = Rmin[j];
1129 }
1130 original_parameters->Start_angle = startPh
1131 original_parameters->Opening_angle = endPh
1132 original_parameters->Num_z_planes = countP
1133
1134 }
1135 else // Set parameters(r,z) with Rmin==0 a
1136 {
1137 #ifdef G4SPECSDEBUG
1138 std::ostringstream message;
1139 message << "Polycone " << GetName() << G4
1140 << "cannot be converted to Polyco
1141 G4Exception("G4Polycone::SetOriginalParam
1142 JustWarning, message);
1143 #endif
1144 original_parameters = new G4PolyconeHisto
1145 original_parameters->Z_values = new G4dou
1146 original_parameters->Rmin = new G4double[
1147 original_parameters->Rmax = new G4double[
1148
1149 for(G4int j=0; j < numPlanes; ++j)
1150 {
1151 original_parameters->Z_values[j] = corn
1152 original_parameters->Rmax[j] = corners[
1153 original_parameters->Rmin[j] = 0.0;
1154 }
1155 original_parameters->Start_angle = startP
1156 original_parameters->Opening_angle = endP
1157 original_parameters->Num_z_planes = numPl
1158 }
1159 return isConvertible;
1160 }
1161
1162 #endif
1163