| 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 // G4 Polyhedron library
27 //
28 // History:
29 // 23.07.96 E.Chernyaev <Evgueni.Tcherniaev@ce
30 //
31 // 30.09.96 E.Chernyaev
32 // - added GetNextVertexIndex, GetVertex by Ya
33 // - added GetNextUnitNormal, GetNextEdgeIndic
34 //
35 // 15.12.96 E.Chernyaev
36 // - added GetNumberOfRotationSteps, RotateEdg
37 // - rewritten G4PolyhedronCons;
38 // - added G4PolyhedronPara, ...Trap, ...Pgon,
39 //
40 // 01.06.97 E.Chernyaev
41 // - modified RotateAroundZ, added SetSideFace
42 //
43 // 19.03.00 E.Chernyaev
44 // - implemented boolean operations (add, subt
45 //
46 // 25.05.01 E.Chernyaev
47 // - added GetSurfaceArea() and GetVolume()
48 //
49 // 05.11.02 E.Chernyaev
50 // - added createTwistedTrap() and createPolyh
51 //
52 // 20.06.05 G.Cosmo
53 // - added HepPolyhedronEllipsoid
54 //
55 // 18.07.07 T.Nikitina
56 // - added HepPolyhedronParaboloid
57 //
58 // 22.02.20 E.Chernyaev
59 // - added HepPolyhedronTet, HepPolyhedronHybe
60 //
61 // 12.05.21 E.Chernyaev
62 // - added TriangulatePolygon(), RotateContour
63 // - added HepPolyhedronPgon, HepPolyhedronPco
64 //
65 // 26.03.22 E.Chernyaev
66 // - added SetVertex(), SetFacet()
67 // - added HepPolyhedronTetMesh
68 //
69 // 04.04.22 E.Chernyaev
70 // - added JoinCoplanarFacets()
71 //
72 // 07.04.22 E.Chernyaev
73 // - added HepPolyhedronBoxMesh
74
75 #include "HepPolyhedron.h"
76 #include "G4PhysicalConstants.hh"
77 #include "G4Vector3D.hh"
78
79 #include <cstdlib> // Required on some compil
80 #include <cmath>
81 #include <algorithm>
82
83 using CLHEP::perMillion;
84 using CLHEP::deg;
85 using CLHEP::pi;
86 using CLHEP::twopi;
87 using CLHEP::nm;
88 const G4double spatialTolerance = 0.01*nm;
89
90 /*********************************************
91 *
92 * Name: HepPolyhedron operator <<
93 * Author: E.Chernyaev (IHEP/Protvino)
94 *
95 * Function: Print contents of G4 polyhedron
96 *
97 *********************************************
98 std::ostream & operator<<(std::ostream & ostr,
99 for (const auto& edge : facet.edge) {
100 ostr << " " << edge.v << "/" << edge.f;
101 }
102 return ostr;
103 }
104
105 std::ostream & operator<<(std::ostream & ostr,
106 ostr << std::endl;
107 ostr << "Nvertices=" << ph.nvert << ", Nface
108 G4int i;
109 for (i=1; i<=ph.nvert; i++) {
110 ostr << "xyz(" << i << ")="
111 << ph.pV[i].x() << ' ' << ph.pV[i].y
112 << std::endl;
113 }
114 for (i=1; i<=ph.nface; i++) {
115 ostr << "face(" << i << ")=" << ph.pF[i] <
116 }
117 return ostr;
118 }
119
120 HepPolyhedron::HepPolyhedron(G4int Nvert, G4in
121 /*********************************************
122 *
123 * Name: HepPolyhedron constructor with
124 * allocation of memory
125 * Author: E.Tcherniaev (E.Chernyaev)
126 *
127 *********************************************
128 : nvert(0), nface(0), pV(nullptr), pF(nullptr)
129 {
130 AllocateMemory(Nvert, Nface);
131 }
132
133 HepPolyhedron::HepPolyhedron(const HepPolyhedr
134 /*********************************************
135 *
136 * Name: HepPolyhedron copy constructor
137 * Author: E.Chernyaev (IHEP/Protvino)
138 *
139 *********************************************
140 : nvert(0), nface(0), pV(nullptr), pF(nullptr)
141 {
142 AllocateMemory(from.nvert, from.nface);
143 for (G4int i=1; i<=nvert; i++) pV[i] = from.
144 for (G4int k=1; k<=nface; k++) pF[k] = from.
145 }
146
147 HepPolyhedron::HepPolyhedron(HepPolyhedron&& f
148 /*********************************************
149 *
150 * Name: HepPolyhedron move constructor
151 * Author: E.Tcherniaev (E.Chernyaev)
152 *
153 *********************************************
154 : nvert(0), nface(0), pV(nullptr), pF(nullptr)
155 {
156 nvert = from.nvert;
157 nface = from.nface;
158 pV = from.pV;
159 pF = from.pF;
160
161 // Release the data from the source object
162 from.nvert = 0;
163 from.nface = 0;
164 from.pV = nullptr;
165 from.pF = nullptr;
166 }
167
168 HepPolyhedron & HepPolyhedron::operator=(const
169 /*********************************************
170 *
171 * Name: HepPolyhedron operator =
172 * Author: E.Chernyaev (IHEP/Protvino)
173 *
174 * Function: Copy contents of one polyhedron t
175 *
176 *********************************************
177 {
178 if (this != &from) {
179 AllocateMemory(from.nvert, from.nface);
180 for (G4int i=1; i<=nvert; i++) pV[i] = fro
181 for (G4int k=1; k<=nface; k++) pF[k] = fro
182 }
183 return *this;
184 }
185
186 HepPolyhedron & HepPolyhedron::operator=(HepPo
187 /*********************************************
188 *
189 * Name: HepPolyhedron move operator =
190 * Author: E.Tcherniaev (E.Chernyaev)
191 *
192 * Function: Move contents of one polyhedron t
193 *
194 *********************************************
195 {
196 if (this != &from) {
197 delete [] pV;
198 delete [] pF;
199 nvert = from.nvert;
200 nface = from.nface;
201 pV = from.pV;
202 pF = from.pF;
203
204 // Release the data from the source object
205 from.nvert = 0;
206 from.nface = 0;
207 from.pV = nullptr;
208 from.pF = nullptr;
209 }
210 return *this;
211 }
212
213 G4int
214 HepPolyhedron::FindNeighbour(G4int iFace, G4in
215 /*********************************************
216 *
217 * Name: HepPolyhedron::FindNeighbour
218 * Author: E.Chernyaev
219 *
220 * Function: Find neighbouring face
221 *
222 *********************************************
223 {
224 G4int i;
225 for (i=0; i<4; i++) {
226 if (iNode == std::abs(pF[iFace].edge[i].v)
227 }
228 if (i == 4) {
229 std::cerr
230 << "HepPolyhedron::FindNeighbour: face "
231 << " has no node " << iNode
232 << std::endl;
233 return 0;
234 }
235 if (iOrder < 0) {
236 if ( --i < 0) i = 3;
237 if (pF[iFace].edge[i].v == 0) i = 2;
238 }
239 return (pF[iFace].edge[i].v > 0) ? 0 : pF[iF
240 }
241
242 G4Normal3D HepPolyhedron::FindNodeNormal(G4int
243 /*********************************************
244 *
245 * Name: HepPolyhedron::FindNodeNormal
246 * Author: E.Chernyaev
247 *
248 * Function: Find normal at given node
249 *
250 *********************************************
251 {
252 G4Normal3D normal = GetUnitNormal(iFace);
253 G4int k = iFace, iOrder = 1;
254
255 for(;;) {
256 k = FindNeighbour(k, iNode, iOrder);
257 if (k == iFace) break;
258 if (k > 0) {
259 normal += GetUnitNormal(k);
260 }else{
261 if (iOrder < 0) break;
262 k = iFace;
263 iOrder = -iOrder;
264 }
265 }
266 return normal.unit();
267 }
268
269 G4int HepPolyhedron::GetNumberOfRotationSteps(
270 /*********************************************
271 *
272 * Name: HepPolyhedron::GetNumberOfRotationSte
273 * Author: J.Allison (Manchester University)
274 *
275 * Function: Get number of steps for whole cir
276 *
277 *********************************************
278 {
279 return fNumberOfRotationSteps;
280 }
281
282 void HepPolyhedron::SetVertex(G4int index, con
283 /*********************************************
284 *
285 * Name: HepPolyhedron::SetVertex
286 * Author: E.Tcherniaev (E.Chernyaev)
287 *
288 * Function: Set vertex
289 *
290 *********************************************
291 {
292 if (index < 1 || index > nvert)
293 {
294 std::cerr
295 << "HepPolyhedron::SetVertex: vertex ind
296 << " is out of range\n"
297 << " N. of vertices = " << nvert << "\
298 << " N. of facets = " << nface << std:
299 return;
300 }
301 pV[index] = v;
302 }
303
304 void
305 HepPolyhedron::SetFacet(G4int index, G4int iv1
306 /*********************************************
307 *
308 * Name: HepPolyhedron::SetFacet
309 * Author: E.Tcherniaev (E.Chernyaev)
310 *
311 * Function: Set facet
312 *
313 *********************************************
314 {
315 if (index < 1 || index > nface)
316 {
317 std::cerr
318 << "HepPolyhedron::SetFacet: facet index
319 << " is out of range\n"
320 << " N. of vertices = " << nvert << "\
321 << " N. of facets = " << nface << std:
322 return;
323 }
324 if (iv1 < 1 || iv1 > nvert ||
325 iv2 < 1 || iv2 > nvert ||
326 iv3 < 1 || iv3 > nvert ||
327 iv4 < 0 || iv4 > nvert)
328 {
329 std::cerr
330 << "HepPolyhedron::SetFacet: incorrectly
331 << " (" << iv1 << ", " << iv2 << ", " <<
332 << " N. of vertices = " << nvert << "\
333 << " N. of facets = " << nface << std:
334 return;
335 }
336 pF[index] = G4Facet(iv1, 0, iv2, 0, iv3, 0,
337 }
338
339 void HepPolyhedron::SetNumberOfRotationSteps(G
340 /*********************************************
341 *
342 * Name: HepPolyhedron::SetNumberOfRotationSte
343 * Author: J.Allison (Manchester University)
344 *
345 * Function: Set number of steps for whole cir
346 *
347 *********************************************
348 {
349 const G4int nMin = 3;
350 if (n < nMin) {
351 std::cerr
352 << "HepPolyhedron::SetNumberOfRotationSt
353 << "number of steps per circle < " << nM
354 << std::endl;
355 fNumberOfRotationSteps = nMin;
356 }else{
357 fNumberOfRotationSteps = n;
358 }
359 }
360
361 void HepPolyhedron::ResetNumberOfRotationSteps
362 /*********************************************
363 *
364 * Name: HepPolyhedron::GetNumberOfRotationSte
365 * Author: J.Allison (Manchester University)
366 *
367 * Function: Reset number of steps for whole c
368 *
369 *********************************************
370 {
371 fNumberOfRotationSteps = DEFAULT_NUMBER_OF_S
372 }
373
374 void HepPolyhedron::AllocateMemory(G4int Nvert
375 /*********************************************
376 *
377 * Name: HepPolyhedron::AllocateMemory
378 * Author: E.Chernyaev (IHEP/Protvino)
379 *
380 * Function: Allocate memory for GEANT4 polyhe
381 *
382 * Input: Nvert - number of nodes
383 * Nface - number of faces
384 *
385 *********************************************
386 {
387 if (nvert == Nvert && nface == Nface) return
388 delete [] pV;
389 delete [] pF;
390 if (Nvert > 0 && Nface > 0) {
391 nvert = Nvert;
392 nface = Nface;
393 pV = new G4Point3D[nvert+1];
394 pF = new G4Facet[nface+1];
395 }else{
396 nvert = 0; nface = 0; pV = nullptr; pF = n
397 }
398 }
399
400 void HepPolyhedron::CreatePrism()
401 /*********************************************
402 *
403 * Name: HepPolyhedron::CreatePrism
404 * Author: E.Chernyaev (IHEP/Protvino)
405 *
406 * Function: Set facets for a prism
407 *
408 *********************************************
409 {
410 enum {DUMMY, BOTTOM, LEFT, BACK, RIGHT, FRON
411
412 pF[1] = G4Facet(1,LEFT, 4,BACK, 3,RIGHT,
413 pF[2] = G4Facet(5,TOP, 8,BACK, 4,BOTTOM,
414 pF[3] = G4Facet(8,TOP, 7,RIGHT, 3,BOTTOM,
415 pF[4] = G4Facet(7,TOP, 6,FRONT, 2,BOTTOM,
416 pF[5] = G4Facet(6,TOP, 5,LEFT, 1,BOTTOM,
417 pF[6] = G4Facet(5,FRONT, 6,RIGHT, 7,BACK,
418 }
419
420 void HepPolyhedron::RotateEdge(G4int k1, G4int
421 G4int v1, G4int
422 G4bool ifWholeCi
423 /*********************************************
424 *
425 * Name: HepPolyhedron::RotateEdge
426 * Author: E.Chernyaev (IHEP/Protvino)
427 *
428 * Function: Create set of facets by rotation
429 *
430 * Input: k1, k2 - end vertices of the edge
431 * r1, r2 - radiuses of the end vertice
432 * v1, v2 - visibility of edges produce
433 * vertices
434 * vEdge - visibility of the edge
435 * ifWholeCircle - is true in case of w
436 * nds - number of discrete steps
437 * r[] - r-coordinates
438 * kface - current free cell in the pF
439 *
440 *********************************************
441 {
442 if (r1 == 0. && r2 == 0.) return;
443
444 G4int i;
445 G4int i1 = k1;
446 G4int i2 = k2;
447 G4int ii1 = ifWholeCircle ? i1 : i1+nds;
448 G4int ii2 = ifWholeCircle ? i2 : i2+nds;
449 G4int vv = ifWholeCircle ? vEdge : 1;
450
451 if (nds == 1) {
452 if (r1 == 0.) {
453 pF[kface++] = G4Facet(i1,0, v2*i2,0
454 }else if (r2 == 0.) {
455 pF[kface++] = G4Facet(i1,0, i2,0,
456 }else{
457 pF[kface++] = G4Facet(i1,0, v2*i2,0
458 }
459 }else{
460 if (r1 == 0.) {
461 pF[kface++] = G4Facet(vv*i1,0, v2*i
462 for (i2++,i=1; i<nds-1; i2++,i++) {
463 pF[kface++] = G4Facet(vEdge*i1,0, v2*i
464 }
465 pF[kface++] = G4Facet(vEdge*i1,0, v2*i
466 }else if (r2 == 0.) {
467 pF[kface++] = G4Facet(vv*i1,0, vEdg
468 for (i1++,i=1; i<nds-1; i1++,i++) {
469 pF[kface++] = G4Facet(vEdge*i1,0, vEdg
470 }
471 pF[kface++] = G4Facet(vEdge*i1,0, vv*i
472 }else{
473 pF[kface++] = G4Facet(vv*i1,0, v2*i
474 for (i1++,i2++,i=1; i<nds-1; i1++,i2++,i
475 pF[kface++] = G4Facet(vEdge*i1,0, v2*i
476 }
477 pF[kface++] = G4Facet(vEdge*i1,0, v2*i
478 }
479 }
480 }
481
482 void HepPolyhedron::SetSideFacets(G4int ii[4],
483 G4int *kk, G4
484 G4double dphi
485 /*********************************************
486 *
487 * Name: HepPolyhedron::SetSideFacets
488 * Author: E.Chernyaev (IHEP/Protvino)
489 *
490 * Function: Set side facets for the case of i
491 *
492 * Input: ii[4] - indices of original vertices
493 * vv[4] - visibility of edges
494 * kk[] - indices of nodes
495 * r[] - radiuses
496 * dphi - delta phi
497 * nds - number of discrete steps
498 * kface - current free cell in the pF
499 *
500 *********************************************
501 {
502 G4int k1, k2, k3, k4;
503
504 if (std::abs(dphi-pi) < perMillion) { // hal
505 for (G4int i=0; i<4; i++) {
506 k1 = ii[i];
507 k2 = ii[(i+1)%4];
508 if (r[k1] == 0. && r[k2] == 0.) vv[i] =
509 }
510 }
511
512 if (ii[1] == ii[2]) {
513 k1 = kk[ii[0]];
514 k2 = kk[ii[2]];
515 k3 = kk[ii[3]];
516 pF[kface++] = G4Facet(vv[0]*k1,0, vv[2]*k2
517 if (r[ii[0]] != 0.) k1 += nds;
518 if (r[ii[2]] != 0.) k2 += nds;
519 if (r[ii[3]] != 0.) k3 += nds;
520 pF[kface++] = G4Facet(vv[2]*k3,0, vv[0]*k2
521 }else if (kk[ii[0]] == kk[ii[1]]) {
522 k1 = kk[ii[0]];
523 k2 = kk[ii[2]];
524 k3 = kk[ii[3]];
525 pF[kface++] = G4Facet(vv[1]*k1,0, vv[2]*k2
526 if (r[ii[0]] != 0.) k1 += nds;
527 if (r[ii[2]] != 0.) k2 += nds;
528 if (r[ii[3]] != 0.) k3 += nds;
529 pF[kface++] = G4Facet(vv[2]*k3,0, vv[1]*k2
530 }else if (kk[ii[2]] == kk[ii[3]]) {
531 k1 = kk[ii[0]];
532 k2 = kk[ii[1]];
533 k3 = kk[ii[2]];
534 pF[kface++] = G4Facet(vv[0]*k1,0, vv[1]*k2
535 if (r[ii[0]] != 0.) k1 += nds;
536 if (r[ii[1]] != 0.) k2 += nds;
537 if (r[ii[2]] != 0.) k3 += nds;
538 pF[kface++] = G4Facet(vv[1]*k3,0, vv[0]*k2
539 }else{
540 k1 = kk[ii[0]];
541 k2 = kk[ii[1]];
542 k3 = kk[ii[2]];
543 k4 = kk[ii[3]];
544 pF[kface++] = G4Facet(vv[0]*k1,0, vv[1]*k2
545 if (r[ii[0]] != 0.) k1 += nds;
546 if (r[ii[1]] != 0.) k2 += nds;
547 if (r[ii[2]] != 0.) k3 += nds;
548 if (r[ii[3]] != 0.) k4 += nds;
549 pF[kface++] = G4Facet(vv[2]*k4,0, vv[1]*k3
550 }
551 }
552
553 void HepPolyhedron::RotateAroundZ(G4int nstep,
554 G4int np1, G4
555 const G4doubl
556 G4int nodeVis
557 /*********************************************
558 *
559 * Name: HepPolyhedron::RotateAroundZ
560 * Author: E.Chernyaev (IHEP/Protvino)
561 *
562 * Function: Create HepPolyhedron for a solid
563 * two polylines around Z-axis
564 *
565 * Input: nstep - number of discrete steps, if
566 * phi - starting phi angle
567 * dphi - delta phi
568 * np1 - number of points in external
569 * (must be negative in case of
570 * np2 - number of points in internal
571 * z[] - z-coordinates (+z >>> -z for
572 * r[] - r-coordinates
573 * nodeVis - how to Draw edges joing co
574 * node during rotation
575 * edgeVis - how to Draw edges
576 *
577 *********************************************
578 {
579 static const G4double wholeCircle = twopi;
580
581 // S E T R O T A T I O N P A R A M E T
582
583 G4bool ifWholeCircle = std::abs(dphi-wholeCi
584 G4double delPhi = ifWholeCircle ? wholeCircl
585 G4int nSphi = nstep;
586 if (nSphi <= 0) nSphi = GetNumberOfRotationS
587 if (nSphi == 0) nSphi = 1;
588 G4int nVphi = ifWholeCircle ? nSphi : nSphi
589 G4bool ifClosed = np1 <= 0; // true if exter
590
591 // C O U N T V E R T I C E S
592
593 G4int absNp1 = std::abs(np1);
594 G4int absNp2 = std::abs(np2);
595 G4int i1beg = 0;
596 G4int i1end = absNp1-1;
597 G4int i2beg = absNp1;
598 G4int i2end = absNp1+absNp2-1;
599 G4int i, j, k;
600
601 for(i=i1beg; i<=i2end; i++) {
602 if (std::abs(r[i]) < spatialTolerance) r[i
603 }
604
605 // external polyline - check position of nod
606 //
607 G4int Nverts = 0;
608 for (i=i1beg; i<=i1end; i++) {
609 Nverts += (r[i] == 0.) ? 1 : nVphi;
610 }
611
612 // internal polyline
613 //
614 G4bool ifSide1 = false; // whether to create
615 G4bool ifSide2 = false; // whether to create
616
617 if (r[i2beg] != r[i1beg] || z[i2beg] != z[i1
618 Nverts += (r[i2beg] == 0.) ? 1 : nVphi;
619 ifSide1 = true;
620 }
621
622 for(i=i2beg+1; i<i2end; i++) { // intermedia
623 Nverts += (r[i] == 0.) ? 1 : nVphi;
624 }
625
626 if (r[i2end] != r[i1end] || z[i2end] != z[i1
627 if (absNp2 > 1) Nverts += (r[i2end] == 0.)
628 ifSide2 = true;
629 }
630
631 // C O U N T F A C E S
632
633 // external lateral faces
634 //
635 G4int Nfaces = ifClosed ? absNp1*nSphi : (ab
636
637 // internal lateral faces
638 //
639 if (absNp2 > 1) {
640 for(i=i2beg; i<i2end; i++) {
641 if (r[i] > 0. || r[i+1] > 0.) Nfaces +=
642 }
643
644 if (ifClosed) {
645 if (r[i2end] > 0. || r[i2beg] > 0.) Nfac
646 }
647 }
648
649 // bottom and top faces
650 //
651 if (!ifClosed) {
652 if (ifSide1 && (r[i1beg] > 0. || r[i2beg]
653 if (ifSide2 && (r[i1end] > 0. || r[i2end]
654 }
655
656 // phi_wedge faces
657 //
658 if (!ifWholeCircle) {
659 Nfaces += ifClosed ? 2*absNp1 : 2*(absNp1-
660 }
661
662 // A L L O C A T E M E M O R Y
663
664 AllocateMemory(Nverts, Nfaces);
665 if (pV == nullptr || pF == nullptr) return;
666
667 // G E N E R A T E V E R T I C E S
668
669 G4int *kk; // array of start indices along p
670 kk = new G4int[absNp1+absNp2];
671
672 // external polyline
673 //
674 k = 1; // free position in array of vertices
675 for(i=i1beg; i<=i1end; i++) {
676 kk[i] = k;
677 if (r[i] == 0.)
678 { pV[k++] = G4Point3D(0, 0, z[i]); } else
679 }
680
681 // first point of internal polyline
682 //
683 i = i2beg;
684 if (ifSide1) {
685 kk[i] = k;
686 if (r[i] == 0.)
687 { pV[k++] = G4Point3D(0, 0, z[i]); } else
688 }else{
689 kk[i] = kk[i1beg];
690 }
691
692 // intermediate points of internal polyline
693 //
694 for(i=i2beg+1; i<i2end; i++) {
695 kk[i] = k;
696 if (r[i] == 0.)
697 { pV[k++] = G4Point3D(0, 0, z[i]); } else
698 }
699
700 // last point of internal polyline
701 //
702 if (absNp2 > 1) {
703 i = i2end;
704 if (ifSide2) {
705 kk[i] = k;
706 if (r[i] == 0.) pV[k] = G4Point3D(0, 0,
707 }else{
708 kk[i] = kk[i1end];
709 }
710 }
711
712 // set vertices
713 //
714 G4double cosPhi, sinPhi;
715
716 for(j=0; j<nVphi; j++) {
717 cosPhi = std::cos(phi+j*delPhi/nSphi);
718 sinPhi = std::sin(phi+j*delPhi/nSphi);
719 for(i=i1beg; i<=i2end; i++) {
720 if (r[i] != 0.)
721 pV[kk[i]+j] = G4Point3D(r[i]*cosPhi,r[
722 }
723 }
724
725 // G E N E R A T E F A C E S
726
727 // external faces
728 //
729 G4int v1,v2;
730
731 k = 1; // free position in array of faces pF
732 v2 = ifClosed ? nodeVis : 1;
733 for(i=i1beg; i<i1end; i++) {
734 v1 = v2;
735 if (!ifClosed && i == i1end-1) {
736 v2 = 1;
737 }else{
738 v2 = (r[i] == r[i+1] && r[i+1] == r[i+2]
739 }
740 RotateEdge(kk[i], kk[i+1], r[i], r[i+1], v
741 edgeVis, ifWholeCircle, nSphi,
742 }
743 if (ifClosed) {
744 RotateEdge(kk[i1end], kk[i1beg], r[i1end],
745 edgeVis, ifWholeCircle, nSphi,
746 }
747
748 // internal faces
749 //
750 if (absNp2 > 1) {
751 v2 = ifClosed ? nodeVis : 1;
752 for(i=i2beg; i<i2end; i++) {
753 v1 = v2;
754 if (!ifClosed && i==i2end-1) {
755 v2 = 1;
756 }else{
757 v2 = (r[i] == r[i+1] && r[i+1] == r[i+
758 }
759 RotateEdge(kk[i+1], kk[i], r[i+1], r[i],
760 edgeVis, ifWholeCircle, nSphi
761 }
762 if (ifClosed) {
763 RotateEdge(kk[i2beg], kk[i2end], r[i2beg
764 edgeVis, ifWholeCircle, nSphi
765 }
766 }
767
768 // bottom and top faces
769 //
770 if (!ifClosed) {
771 if (ifSide1) {
772 RotateEdge(kk[i2beg], kk[i1beg], r[i2beg
773 -1, ifWholeCircle, nSphi, k);
774 }
775 if (ifSide2) {
776 RotateEdge(kk[i1end], kk[i2end], r[i1end
777 -1, ifWholeCircle, nSphi, k);
778 }
779 }
780
781 // phi_wedge faces in case of incomplete cir
782 //
783 if (!ifWholeCircle) {
784
785 G4int ii[4], vv[4];
786
787 if (ifClosed) {
788 for (i=i1beg; i<=i1end; i++) {
789 ii[0] = i;
790 ii[3] = (i == i1end) ? i1beg : i+1;
791 ii[1] = (absNp2 == 1) ? i2beg : ii[0]+
792 ii[2] = (absNp2 == 1) ? i2beg : ii[3]+
793 vv[0] = -1;
794 vv[1] = 1;
795 vv[2] = -1;
796 vv[3] = 1;
797 SetSideFacets(ii, vv, kk, r, delPhi, n
798 }
799 }else{
800 for (i=i1beg; i<i1end; i++) {
801 ii[0] = i;
802 ii[3] = i+1;
803 ii[1] = (absNp2 == 1) ? i2beg : ii[0]+
804 ii[2] = (absNp2 == 1) ? i2beg : ii[3]+
805 vv[0] = (i == i1beg) ? 1 : -1;
806 vv[1] = 1;
807 vv[2] = (i == i1end-1) ? 1 : -1;
808 vv[3] = 1;
809 SetSideFacets(ii, vv, kk, r, delPhi, n
810 }
811 }
812 }
813
814 delete [] kk; // free memory
815
816 // final check
817 //
818 if (k-1 != nface) {
819 std::cerr
820 << "HepPolyhedron::RotateAroundZ: number
821 << k-1 << ") is not equal to the number
822 << nface << ")"
823 << std::endl;
824 }
825 }
826
827 void
828 HepPolyhedron::RotateContourAroundZ(G4int nste
829 G4double p
830 G4double d
831 const std:
832 G4int node
833 G4int edge
834 /*********************************************
835 *
836 * Name: HepPolyhedron::RotateContourAroundZ
837 * Author: E.Tcherniaev (E.Chernyaev)
838 *
839 * Function: Create HepPolyhedron for a solid
840 * a closed polyline (rz-contour) ar
841 *
842 * Input: nstep - number of discrete steps, if
843 * phi - starting phi angle
844 * dphi - delta phi
845 * rz - rz-contour
846 * nodeVis - how to Draw edges joing co
847 * node during rotation
848 * edgeVis - how to Draw edges
849 *
850 *********************************************
851 {
852 // S E T R O T A T I O N P A R A M E T
853
854 G4bool ifWholeCircle = std::abs(dphi - twopi
855 G4double delPhi = (ifWholeCircle) ? twopi :
856 G4int nSphi = nstep;
857 if (nSphi <= 0) nSphi = GetNumberOfRotationS
858 if (nSphi == 0) nSphi = 1;
859 G4int nVphi = (ifWholeCircle) ? nSphi : nSph
860
861 // C A L C U L A T E A R E A
862
863 G4int Nrz = (G4int)rz.size();
864 G4double area = 0;
865 for (G4int i = 0; i < Nrz; ++i)
866 {
867 G4int k = (i == 0) ? Nrz - 1 : i - 1;
868 area += rz[k].x()*rz[i].y() - rz[i].x()*rz
869 }
870
871 // P R E P A R E P O L Y L I N E
872
873 auto r = new G4double[Nrz];
874 auto z = new G4double[Nrz];
875 for (G4int i = 0; i < Nrz; ++i)
876 {
877 r[i] = rz[i].x();
878 z[i] = rz[i].y();
879 if (std::abs(r[i]) < spatialTolerance) r[i
880 }
881
882 // C O U N T V E R T I C E S A N D F
883
884 G4int Nverts = 0;
885 for(G4int i = 0; i < Nrz; ++i) Nverts += (r[
886
887 G4int Nedges = Nrz;
888 for (G4int i = 0; i < Nrz; ++i)
889 {
890 G4int k = (i == 0) ? Nrz - 1 : i - 1;
891 Nedges -= static_cast<int>(r[k] == 0 && r[
892 }
893
894 G4int Nfaces = Nedges*nSphi; /
895 if (!ifWholeCircle) Nfaces += 2*(Nrz - 2); /
896
897 // A L L O C A T E M E M O R Y
898
899 AllocateMemory(Nverts, Nfaces);
900 if (pV == nullptr || pF == nullptr)
901 {
902 delete [] r;
903 delete [] z;
904 return;
905 }
906
907 // S E T V E R T I C E S
908
909 auto kk = new G4int[Nrz]; // start indices a
910 G4int kfree = 1; // current free position in
911
912 // set start indices, set vertices for nodes
913 for(G4int i = 0; i < Nrz; ++i)
914 {
915 kk[i] = kfree;
916 if (r[i] == 0.) pV[kfree++] = G4Point3D(0,
917 if (r[i] != 0.) kfree += nVphi;
918 }
919
920 // set vertices by rotating r
921 for(G4int j = 0; j < nVphi; ++j)
922 {
923 G4double cosPhi = std::cos(phi + j*delPhi/
924 G4double sinPhi = std::sin(phi + j*delPhi/
925 for(G4int i = 0; i < Nrz; ++i)
926 {
927 if (r[i] != 0.)
928 pV[kk[i] + j] = G4Point3D(r[i]*cosPhi,
929 }
930 }
931
932 // S E T F A C E S
933
934 kfree = 1; // current free position in array
935 for(G4int i = 0; i < Nrz; ++i)
936 {
937 G4int i1 = (i < Nrz - 1) ? i + 1 : 0; // i
938 G4int i2 = i;
939 if (area < 0.) std::swap(i1, i2);
940 RotateEdge(kk[i1], kk[i2], r[i1], r[i2], n
941 edgeVis, ifWholeCircle, nSphi,
942 }
943
944 // S E T P H I _ W E D G E F A C E S
945
946 if (!ifWholeCircle)
947 {
948 std::vector<G4int> triangles;
949 TriangulatePolygon(rz, triangles);
950
951 G4int ii[4], vv[4];
952 G4int ntria = G4int(triangles.size()/3);
953 for (G4int i = 0; i < ntria; ++i)
954 {
955 G4int i1 = triangles[0 + i*3];
956 G4int i2 = triangles[1 + i*3];
957 G4int i3 = triangles[2 + i*3];
958 if (area < 0.) std::swap(i1, i3);
959 G4int v1 = (std::abs(i2-i1) == 1 || std:
960 G4int v2 = (std::abs(i3-i2) == 1 || std:
961 G4int v3 = (std::abs(i1-i3) == 1 || std:
962 ii[0] = i1; ii[1] = i2; ii[2] = i2; ii[3
963 vv[0] = v1; vv[1] = -1; vv[2] = v2; vv[3
964 SetSideFacets(ii, vv, kk, r, delPhi, nSp
965 }
966 }
967
968 // free memory
969 delete [] r;
970 delete [] z;
971 delete [] kk;
972
973 // final check
974 if (kfree - 1 != nface)
975 {
976 std::cerr
977 << "HepPolyhedron::RotateContourAroundZ:
978 << kfree-1 << ") is not equal to the num
979 << nface << ")"
980 << std::endl;
981 }
982 }
983
984 G4bool
985 HepPolyhedron::TriangulatePolygon(const std::v
986 std::vector<
987 /*********************************************
988 *
989 * Name: HepPolyhedron::TriangulatePolygon
990 * Author: E.Tcherniaev (E.Chernyaev)
991 *
992 * Function: Simple implementation of "ear cli
993 * triangulation of a simple contour
994 * the result in a std::vector as tr
995 *
996 * If triangulation is sucsessfull t
997 * returns true, otherwise false
998 *
999 * Remark: It's a copy of G4GeomTools::Trian
1000 *
1001 ********************************************
1002 {
1003 result.resize(0);
1004 G4int n = (G4int)polygon.size();
1005 if (n < 3) return false;
1006
1007 // calculate area
1008 //
1009 G4double area = 0.;
1010 for(G4int i = 0; i < n; ++i)
1011 {
1012 G4int k = (i == 0) ? n - 1 : i - 1;
1013 area += polygon[k].x()*polygon[i].y() - p
1014 }
1015
1016 // allocate and initialize list of Vertices
1017 // we want a counter-clockwise polygon in V
1018 //
1019 auto V = new G4int[n];
1020 if (area > 0.)
1021 for (G4int i = 0; i < n; ++i) V[i] = i;
1022 else
1023 for (G4int i = 0; i < n; ++i) V[i] = (n -
1024
1025 // Triangulation: remove nv-2 Vertices, cr
1026 //
1027 G4int nv = n;
1028 G4int count = 2*nv; // error detection coun
1029 for(G4int b = nv - 1; nv > 2; )
1030 {
1031 // ERROR: if we loop, it is probably a no
1032 if ((count--) <= 0)
1033 {
1034 delete [] V;
1035 if (area < 0.) std::reverse(result.begi
1036 return false;
1037 }
1038
1039 // three consecutive vertices in current
1040 G4int a = (b < nv) ? b : 0; // previo
1041 b = (a+1 < nv) ? a+1 : 0; // curren
1042 G4int c = (b+1 < nv) ? b+1 : 0; // next
1043
1044 if (CheckSnip(polygon, a,b,c, nv,V))
1045 {
1046 // output Triangle
1047 result.push_back(V[a]);
1048 result.push_back(V[b]);
1049 result.push_back(V[c]);
1050
1051 // remove vertex b from remaining polyg
1052 nv--;
1053 for(G4int i = b; i < nv; ++i) V[i] = V[
1054
1055 count = 2*nv; // resest error detection
1056 }
1057 }
1058 delete [] V;
1059 if (area < 0.) std::reverse(result.begin(),
1060 return true;
1061 }
1062
1063 G4bool HepPolyhedron::CheckSnip(const std::ve
1064 G4int a, G4in
1065 G4int n, cons
1066 /********************************************
1067 *
1068 * Name: HepPolyhedron::CheckSnip
1069 * Author: E.Tcherniaev (E.Chernyaev)
1070 *
1071 * Function: Check for a valid snip,
1072 * it is a helper functionfor Trian
1073 *
1074 ********************************************
1075 {
1076 static const G4double kCarTolerance = 1.e-9
1077
1078 // check orientation of Triangle
1079 G4double Ax = contour[V[a]].x(), Ay = conto
1080 G4double Bx = contour[V[b]].x(), By = conto
1081 G4double Cx = contour[V[c]].x(), Cy = conto
1082 if ((Bx-Ax)*(Cy-Ay) - (By-Ay)*(Cx-Ax) < kCa
1083
1084 // check that there is no point inside Tria
1085 G4double xmin = std::min(std::min(Ax,Bx),Cx
1086 G4double xmax = std::max(std::max(Ax,Bx),Cx
1087 G4double ymin = std::min(std::min(Ay,By),Cy
1088 G4double ymax = std::max(std::max(Ay,By),Cy
1089
1090 for (G4int i=0; i<n; ++i)
1091 {
1092 if((i == a) || (i == b) || (i == c)) cont
1093 G4double Px = contour[V[i]].x();
1094 if (Px < xmin || Px > xmax) continue;
1095 G4double Py = contour[V[i]].y();
1096 if (Py < ymin || Py > ymax) continue;
1097 // if (PointInTriangle(Ax,Ay,Bx,By,Cx,Cy,
1098 if ((Bx-Ax)*(Cy-Ay) - (By-Ay)*(Cx-Ax) > 0
1099 {
1100 if ((Ax-Cx)*(Py-Cy) - (Ay-Cy)*(Px-Cx) <
1101 if ((Bx-Ax)*(Py-Ay) - (By-Ay)*(Px-Ax) <
1102 if ((Cx-Bx)*(Py-By) - (Cy-By)*(Px-Bx) <
1103 }
1104 else
1105 {
1106 if ((Ax-Cx)*(Py-Cy) - (Ay-Cy)*(Px-Cx) >
1107 if ((Bx-Ax)*(Py-Ay) - (By-Ay)*(Px-Ax) >
1108 if ((Cx-Bx)*(Py-By) - (Cy-By)*(Px-Bx) >
1109 }
1110 return false;
1111 }
1112 return true;
1113 }
1114
1115 void HepPolyhedron::SetReferences()
1116 /********************************************
1117 *
1118 * Name: HepPolyhedron::SetReferences
1119 * Author: E.Chernyaev (IHEP/Protvino)
1120 *
1121 * Function: For each edge set reference to n
1122 *
1123 ********************************************
1124 {
1125 if (nface <= 0) return;
1126
1127 struct edgeListMember {
1128 edgeListMember *next;
1129 G4int v2;
1130 G4int iface;
1131 G4int iedge;
1132 } *edgeList, *freeList, **headList;
1133
1134
1135 // A L L O C A T E A N D I N I T I A
1136
1137 edgeList = new edgeListMember[2*nface];
1138 headList = new edgeListMember*[nvert];
1139
1140 G4int i;
1141 for (i=0; i<nvert; i++) {
1142 headList[i] = nullptr;
1143 }
1144 freeList = edgeList;
1145 for (i=0; i<2*nface-1; i++) {
1146 edgeList[i].next = &edgeList[i+1];
1147 }
1148 edgeList[2*nface-1].next = nullptr;
1149
1150 // L O O P A L O N G E D G E S
1151
1152 G4int iface, iedge, nedge, i1, i2, k1, k2;
1153 edgeListMember *prev, *cur;
1154
1155 for(iface=1; iface<=nface; iface++) {
1156 nedge = (pF[iface].edge[3].v == 0) ? 3 :
1157 for (iedge=0; iedge<nedge; iedge++) {
1158 i1 = iedge;
1159 i2 = (iedge < nedge-1) ? iedge+1 : 0;
1160 i1 = std::abs(pF[iface].edge[i1].v);
1161 i2 = std::abs(pF[iface].edge[i2].v);
1162 k1 = (i1 < i2) ? i1 : i2; // k
1163 k2 = (i1 > i2) ? i1 : i2; // k
1164
1165 // check head of the List corresponding
1166 cur = headList[k1];
1167 if (cur == nullptr) {
1168 headList[k1] = freeList;
1169 if (freeList == nullptr) {
1170 std::cerr
1171 << "Polyhedron::SetReferences: bad
1172 << std::endl;
1173 break;
1174 }
1175 freeList = freeList->next;
1176 cur = headList[k1];
1177 cur->next = nullptr;
1178 cur->v2 = k2;
1179 cur->iface = iface;
1180 cur->iedge = iedge;
1181 continue;
1182 }
1183
1184 if (cur->v2 == k2) {
1185 headList[k1] = cur->next;
1186 cur->next = freeList;
1187 freeList = cur;
1188 pF[iface].edge[iedge].f = cur->iface;
1189 pF[cur->iface].edge[cur->iedge].f = i
1190 i1 = (pF[iface].edge[iedge].v < 0) ?
1191 i2 = (pF[cur->iface].edge[cur->iedge]
1192 if (i1 != i2) {
1193 std::cerr
1194 << "Polyhedron::SetReferences: di
1195 << iface << "/" << iedge << "/"
1196 << pF[iface].edge[iedge].v << " a
1197 << cur->iface << "/" << cur->iedg
1198 << pF[cur->iface].edge[cur->iedge
1199 << std::endl;
1200 }
1201 continue;
1202 }
1203
1204 // check List itself
1205 for (;;) {
1206 prev = cur;
1207 cur = prev->next;
1208 if (cur == nullptr) {
1209 prev->next = freeList;
1210 if (freeList == nullptr) {
1211 std::cerr
1212 << "Polyhedron::SetReferences: ba
1213 << std::endl;
1214 break;
1215 }
1216 freeList = freeList->next;
1217 cur = prev->next;
1218 cur->next = nullptr;
1219 cur->v2 = k2;
1220 cur->iface = iface;
1221 cur->iedge = iedge;
1222 break;
1223 }
1224
1225 if (cur->v2 == k2) {
1226 prev->next = cur->next;
1227 cur->next = freeList;
1228 freeList = cur;
1229 pF[iface].edge[iedge].f = cur->ifac
1230 pF[cur->iface].edge[cur->iedge].f =
1231 i1 = (pF[iface].edge[iedge].v < 0)
1232 i2 = (pF[cur->iface].edge[cur->iedg
1233 if (i1 != i2) {
1234 std::cerr
1235 << "Polyhedron::SetReferences
1236 << iface << "/" << iedge << "
1237 << pF[iface].edge[iedge].v <<
1238 << cur->iface << "/" << cur->
1239 << pF[cur->iface].edge[cur->i
1240 << std::endl;
1241 }
1242 break;
1243 }
1244 }
1245 }
1246 }
1247
1248 // C H E C K T H A T A L L L I S T S
1249
1250 for (i=0; i<nvert; i++) {
1251 if (headList[i] != nullptr) {
1252 std::cerr
1253 << "Polyhedron::SetReferences: List "
1254 << std::endl;
1255 }
1256 }
1257
1258 // F R E E M E M O R Y
1259
1260 delete [] edgeList;
1261 delete [] headList;
1262 }
1263
1264 void HepPolyhedron::JoinCoplanarFacets(G4doub
1265 /********************************************
1266 *
1267 * Name: HepPolyhedron::JoinCoplanarFacets
1268 * Author: E.Tcherniaev (E.Chernyaev)
1269 *
1270 * Function: Join couples of triangular facet
1271 * where it is possible
1272 *
1273 ********************************************
1274 {
1275 G4int njoin = 0;
1276 for (G4int icur = 1; icur <= nface; ++icur)
1277 {
1278 // skip if already joined or quadrangle
1279 if (pF[icur].edge[0].v == 0) continue;
1280 if (pF[icur].edge[3].v != 0) continue;
1281 // skip if all references point to alread
1282 if (pF[icur].edge[0].f < icur &&
1283 pF[icur].edge[1].f < icur &&
1284 pF[icur].edge[2].f < icur) continue;
1285 // compute plane equation
1286 G4Normal3D norm = GetUnitNormal(icur);
1287 G4double dd = norm.dot(pV[pF[icur].edge[0
1288 G4int vcur0 = std::abs(pF[icur].edge[0].v
1289 G4int vcur1 = std::abs(pF[icur].edge[1].v
1290 G4int vcur2 = std::abs(pF[icur].edge[2].v
1291 // select neighbouring facet
1292 G4int kcheck = 0, icheck = 0, vcheck = 0;
1293 G4double dist = DBL_MAX;
1294 for (G4int k = 0; k < 3; ++k)
1295 {
1296 G4int itmp = pF[icur].edge[k].f;
1297 // skip if already checked, joined or q
1298 if (itmp < icur) continue;
1299 if (pF[itmp].edge[0].v == 0 ||
1300 pF[itmp].edge[3].v != 0) continue;
1301 // get candidate vertex
1302 G4int vtmp = 0;
1303 for (G4int j = 0; j < 3; ++j)
1304 {
1305 vtmp = std::abs(pF[itmp].edge[j].v);
1306 if (vtmp != vcur0 && vtmp != vcur1 && vtmp
1307 }
1308 // check distance to the plane
1309 G4double dtmp = std::abs(norm.dot(pV[vt
1310 if (dtmp > tolerance || dtmp >= dist) c
1311 dist = dtmp;
1312 kcheck = k;
1313 icheck = itmp;
1314 vcheck = vtmp;
1315 }
1316 if (icheck == 0) continue; // no facet se
1317 // join facets
1318 njoin++;
1319 pF[icheck].edge[0].v = 0; // mark facet a
1320 if (kcheck == 0)
1321 {
1322 pF[icur].edge[3].v = pF[icur].edge[2].v
1323 pF[icur].edge[2].v = pF[icur].edge[1].v
1324 pF[icur].edge[1].v = vcheck;
1325 }
1326 else if (kcheck == 1)
1327 {
1328 pF[icur].edge[3].v = pF[icur].edge[2].v
1329 pF[icur].edge[2].v = vcheck;
1330 }
1331 else
1332 {
1333 pF[icur].edge[3].v = vcheck;
1334 }
1335 }
1336 if (njoin == 0) return; // no joined facets
1337
1338 // restructure facets
1339 G4int nnew = 0;
1340 for (G4int icur = 1; icur <= nface; ++icur)
1341 {
1342 if (pF[icur].edge[0].v == 0) continue;
1343 nnew++;
1344 pF[nnew].edge[0].v = pF[icur].edge[0].v;
1345 pF[nnew].edge[1].v = pF[icur].edge[1].v;
1346 pF[nnew].edge[2].v = pF[icur].edge[2].v;
1347 pF[nnew].edge[3].v = pF[icur].edge[3].v;
1348 }
1349 nface = nnew;
1350 SetReferences();
1351 }
1352
1353 void HepPolyhedron::InvertFacets()
1354 /********************************************
1355 *
1356 * Name: HepPolyhedron::InvertFacets
1357 * Author: E.Chernyaev
1358 *
1359 * Function: Invert the order of the nodes in
1360 *
1361 ********************************************
1362 {
1363 if (nface <= 0) return;
1364 G4int i, k, nnode, v[4],f[4];
1365 for (i=1; i<=nface; i++) {
1366 nnode = (pF[i].edge[3].v == 0) ? 3 : 4;
1367 for (k=0; k<nnode; k++) {
1368 v[k] = (k+1 == nnode) ? pF[i].edge[0].v
1369 if (v[k] * pF[i].edge[k].v < 0) v[k] =
1370 f[k] = pF[i].edge[k].f;
1371 }
1372 for (k=0; k<nnode; k++) {
1373 pF[i].edge[nnode-1-k].v = v[k];
1374 pF[i].edge[nnode-1-k].f = f[k];
1375 }
1376 }
1377 }
1378
1379 HepPolyhedron & HepPolyhedron::Transform(cons
1380 /********************************************
1381 *
1382 * Name: HepPolyhedron::Transform
1383 * Author: E.Chernyaev
1384 *
1385 * Function: Make transformation of the polyh
1386 *
1387 ********************************************
1388 {
1389 if (nvert > 0) {
1390 for (G4int i=1; i<=nvert; i++) { pV[i] =
1391
1392 // C H E C K D E T E R M I N A N T A
1393 // I N V E R T F A C E T S I F I T
1394
1395 G4Vector3D d = t * G4Vector3D(0,0,0);
1396 G4Vector3D x = t * G4Vector3D(1,0,0) - d;
1397 G4Vector3D y = t * G4Vector3D(0,1,0) - d;
1398 G4Vector3D z = t * G4Vector3D(0,0,1) - d;
1399 if ((x.cross(y))*z < 0) InvertFacets();
1400 }
1401 return *this;
1402 }
1403
1404 G4bool HepPolyhedron::GetNextVertexIndex(G4in
1405 /********************************************
1406 *
1407 * Name: HepPolyhedron::GetNextVertexIndex
1408 * Author: Yasuhide Sawada
1409 *
1410 * Function:
1411 *
1412 ********************************************
1413 {
1414 static G4ThreadLocal G4int iFace = 1;
1415 static G4ThreadLocal G4int iQVertex = 0;
1416 G4int vIndex = pF[iFace].edge[iQVertex].v;
1417
1418 edgeFlag = (vIndex > 0) ? 1 : 0;
1419 index = std::abs(vIndex);
1420
1421 if (iQVertex >= 3 || pF[iFace].edge[iQVerte
1422 iQVertex = 0;
1423 if (++iFace > nface) iFace = 1;
1424 return false; // Last Edge
1425 }
1426
1427 ++iQVertex;
1428 return true; // not Last Edge
1429 }
1430
1431 G4Point3D HepPolyhedron::GetVertex(G4int inde
1432 /********************************************
1433 *
1434 * Name: HepPolyhedron::GetVertex
1435 * Author: Yasuhide Sawada
1436 *
1437 * Function: Get vertex of the index.
1438 *
1439 ********************************************
1440 {
1441 if (index <= 0 || index > nvert) {
1442 std::cerr
1443 << "HepPolyhedron::GetVertex: irrelevan
1444 << std::endl;
1445 return G4Point3D();
1446 }
1447 return pV[index];
1448 }
1449
1450 G4bool
1451 HepPolyhedron::GetNextVertex(G4Point3D &verte
1452 /********************************************
1453 *
1454 * Name: HepPolyhedron::GetNextVertex
1455 * Author: John Allison
1456 *
1457 * Function: Get vertices of the quadrilatera
1458 * face in face order. Returns fal
1459 * face.
1460 *
1461 ********************************************
1462 {
1463 G4int index;
1464 G4bool rep = GetNextVertexIndex(index, edge
1465 vertex = pV[index];
1466 return rep;
1467 }
1468
1469 G4bool HepPolyhedron::GetNextVertex(G4Point3D
1470 G4Normal3D
1471 /********************************************
1472 *
1473 * Name: HepPolyhedron::GetNextVertex
1474 * Author: E.Chernyaev
1475 *
1476 * Function: Get vertices with normals of the
1477 * for each face in face order.
1478 * Returns false when finished each
1479 *
1480 ********************************************
1481 {
1482 static G4ThreadLocal G4int iFace = 1;
1483 static G4ThreadLocal G4int iNode = 0;
1484
1485 if (nface == 0) return false; // empty pol
1486
1487 G4int k = pF[iFace].edge[iNode].v;
1488 if (k > 0) { edgeFlag = 1; } else { edgeFla
1489 vertex = pV[k];
1490 normal = FindNodeNormal(iFace,k);
1491 if (iNode >= 3 || pF[iFace].edge[iNode+1].v
1492 iNode = 0;
1493 if (++iFace > nface) iFace = 1;
1494 return false; // last node
1495 }
1496 ++iNode;
1497 return true; // not last no
1498 }
1499
1500 G4bool HepPolyhedron::GetNextEdgeIndices(G4in
1501 G4int
1502 /********************************************
1503 *
1504 * Name: HepPolyhedron::GetNextEdgeIndices
1505 * Author: E.Chernyaev
1506 *
1507 * Function: Get indices of the next edge tog
1508 * of the faces which share the edg
1509 * Returns false when the last edge
1510 *
1511 ********************************************
1512 {
1513 static G4ThreadLocal G4int iFace = 1;
1514 static G4ThreadLocal G4int iQVertex = 0;
1515 static G4ThreadLocal G4int iOrder = 1;
1516 G4int k1, k2, kflag, kface1, kface2;
1517
1518 if (iFace == 1 && iQVertex == 0) {
1519 k2 = pF[nface].edge[0].v;
1520 k1 = pF[nface].edge[3].v;
1521 if (k1 == 0) k1 = pF[nface].edge[2].v;
1522 if (std::abs(k1) > std::abs(k2)) iOrder =
1523 }
1524
1525 do {
1526 k1 = pF[iFace].edge[iQVertex].v;
1527 kflag = k1;
1528 k1 = std::abs(k1);
1529 kface1 = iFace;
1530 kface2 = pF[iFace].edge[iQVertex].f;
1531 if (iQVertex >= 3 || pF[iFace].edge[iQVer
1532 iQVertex = 0;
1533 k2 = std::abs(pF[iFace].edge[iQVertex].
1534 iFace++;
1535 }else{
1536 iQVertex++;
1537 k2 = std::abs(pF[iFace].edge[iQVertex].
1538 }
1539 } while (iOrder*k1 > iOrder*k2);
1540
1541 i1 = k1; i2 = k2; edgeFlag = (kflag > 0) ?
1542 iface1 = kface1; iface2 = kface2;
1543
1544 if (iFace > nface) {
1545 iFace = 1; iOrder = 1;
1546 return false;
1547 }
1548
1549 return true;
1550 }
1551
1552 G4bool
1553 HepPolyhedron::GetNextEdgeIndices(G4int &i1,
1554 /********************************************
1555 *
1556 * Name: HepPolyhedron::GetNextEdgeIndices
1557 * Author: E.Chernyaev
1558 *
1559 * Function: Get indices of the next edge.
1560 * Returns false when the last edge
1561 *
1562 ********************************************
1563 {
1564 G4int kface1, kface2;
1565 return GetNextEdgeIndices(i1, i2, edgeFlag,
1566 }
1567
1568 G4bool
1569 HepPolyhedron::GetNextEdge(G4Point3D &p1,
1570 G4Point3D &p2,
1571 G4int &edgeFlag) c
1572 /********************************************
1573 *
1574 * Name: HepPolyhedron::GetNextEdge
1575 * Author: E.Chernyaev
1576 *
1577 * Function: Get next edge.
1578 * Returns false when the last edge
1579 *
1580 ********************************************
1581 {
1582 G4int i1,i2;
1583 G4bool rep = GetNextEdgeIndices(i1,i2,edgeF
1584 p1 = pV[i1];
1585 p2 = pV[i2];
1586 return rep;
1587 }
1588
1589 G4bool
1590 HepPolyhedron::GetNextEdge(G4Point3D &p1, G4P
1591 G4int &edgeFlag, G4
1592 /********************************************
1593 *
1594 * Name: HepPolyhedron::GetNextEdge
1595 * Author: E.Chernyaev
1596 *
1597 * Function: Get next edge with indices of th
1598 * the edge.
1599 * Returns false when the last edge
1600 *
1601 ********************************************
1602 {
1603 G4int i1,i2;
1604 G4bool rep = GetNextEdgeIndices(i1,i2,edgeF
1605 p1 = pV[i1];
1606 p2 = pV[i2];
1607 return rep;
1608 }
1609
1610 void HepPolyhedron::GetFacet(G4int iFace, G4i
1611 G4int *edgeFlags,
1612 /********************************************
1613 *
1614 * Name: HepPolyhedron::GetFacet
1615 * Author: E.Chernyaev
1616 *
1617 * Function: Get face by index
1618 *
1619 ********************************************
1620 {
1621 if (iFace < 1 || iFace > nface) {
1622 std::cerr
1623 << "HepPolyhedron::GetFacet: irrelevant
1624 << std::endl;
1625 n = 0;
1626 }else{
1627 G4int i, k;
1628 for (i=0; i<4; i++) {
1629 k = pF[iFace].edge[i].v;
1630 if (k == 0) break;
1631 if (iFaces != nullptr) iFaces[i] = pF[i
1632 if (k > 0) {
1633 iNodes[i] = k;
1634 if (edgeFlags != nullptr) edgeFlags[i
1635 }else{
1636 iNodes[i] = -k;
1637 if (edgeFlags != nullptr) edgeFlags[i
1638 }
1639 }
1640 n = i;
1641 }
1642 }
1643
1644 void HepPolyhedron::GetFacet(G4int index, G4i
1645 G4int *edgeFlags
1646 /********************************************
1647 *
1648 * Name: HepPolyhedron::GetFacet
1649 * Author: E.Chernyaev
1650 *
1651 * Function: Get face by index
1652 *
1653 ********************************************
1654 {
1655 G4int iNodes[4];
1656 GetFacet(index, n, iNodes, edgeFlags);
1657 if (n != 0) {
1658 for (G4int i=0; i<n; i++) {
1659 nodes[i] = pV[iNodes[i]];
1660 if (normals != nullptr) normals[i] = Fi
1661 }
1662 }
1663 }
1664
1665 G4bool
1666 HepPolyhedron::GetNextFacet(G4int &n, G4Point
1667 G4int *edgeFlags,
1668 /********************************************
1669 *
1670 * Name: HepPolyhedron::GetNextFacet
1671 * Author: E.Chernyaev
1672 *
1673 * Function: Get next face with normals of un
1674 * Returns false when finished all
1675 *
1676 ********************************************
1677 {
1678 static G4ThreadLocal G4int iFace = 1;
1679
1680 if (edgeFlags == nullptr) {
1681 GetFacet(iFace, n, nodes);
1682 }else if (normals == nullptr) {
1683 GetFacet(iFace, n, nodes, edgeFlags);
1684 }else{
1685 GetFacet(iFace, n, nodes, edgeFlags, norm
1686 }
1687
1688 if (++iFace > nface) {
1689 iFace = 1;
1690 return false;
1691 }
1692
1693 return true;
1694 }
1695
1696 G4Normal3D HepPolyhedron::GetNormal(G4int iFa
1697 /********************************************
1698 *
1699 * Name: HepPolyhedron::GetNormal
1700 * Author: E.Chernyaev
1701 *
1702 * Function: Get normal of the face given by
1703 *
1704 ********************************************
1705 {
1706 if (iFace < 1 || iFace > nface) {
1707 std::cerr
1708 << "HepPolyhedron::GetNormal: irrelevan
1709 << std::endl;
1710 return G4Normal3D();
1711 }
1712
1713 G4int i0 = std::abs(pF[iFace].edge[0].v);
1714 G4int i1 = std::abs(pF[iFace].edge[1].v);
1715 G4int i2 = std::abs(pF[iFace].edge[2].v);
1716 G4int i3 = std::abs(pF[iFace].edge[3].v);
1717 if (i3 == 0) i3 = i0;
1718 return (pV[i2] - pV[i0]).cross(pV[i3] - pV[
1719 }
1720
1721 G4Normal3D HepPolyhedron::GetUnitNormal(G4int
1722 /********************************************
1723 *
1724 * Name: HepPolyhedron::GetNormal
1725 * Author: E.Chernyaev
1726 *
1727 * Function: Get unit normal of the face give
1728 *
1729 ********************************************
1730 {
1731 if (iFace < 1 || iFace > nface) {
1732 std::cerr
1733 << "HepPolyhedron::GetUnitNormal: irrel
1734 << std::endl;
1735 return G4Normal3D();
1736 }
1737
1738 G4int i0 = std::abs(pF[iFace].edge[0].v);
1739 G4int i1 = std::abs(pF[iFace].edge[1].v);
1740 G4int i2 = std::abs(pF[iFace].edge[2].v);
1741 G4int i3 = std::abs(pF[iFace].edge[3].v);
1742 if (i3 == 0) i3 = i0;
1743 return ((pV[i2] - pV[i0]).cross(pV[i3] - pV
1744 }
1745
1746 G4bool HepPolyhedron::GetNextNormal(G4Normal3
1747 /********************************************
1748 *
1749 * Name: HepPolyhedron::GetNextNormal
1750 * Author: John Allison
1751 *
1752 * Function: Get normals of each face in face
1753 * when finished all faces.
1754 *
1755 ********************************************
1756 {
1757 static G4ThreadLocal G4int iFace = 1;
1758 normal = GetNormal(iFace);
1759 if (++iFace > nface) {
1760 iFace = 1;
1761 return false;
1762 }
1763 return true;
1764 }
1765
1766 G4bool HepPolyhedron::GetNextUnitNormal(G4Nor
1767 /********************************************
1768 *
1769 * Name: HepPolyhedron::GetNextUnitNormal
1770 * Author: E.Chernyaev
1771 *
1772 * Function: Get normals of unit length of ea
1773 * Returns false when finished all
1774 *
1775 ********************************************
1776 {
1777 G4bool rep = GetNextNormal(normal);
1778 normal = normal.unit();
1779 return rep;
1780 }
1781
1782 G4double HepPolyhedron::GetSurfaceArea() cons
1783 /********************************************
1784 *
1785 * Name: HepPolyhedron::GetSurfaceArea
1786 * Author: E.Chernyaev
1787 *
1788 * Function: Returns area of the surface of t
1789 *
1790 ********************************************
1791 {
1792 G4double srf = 0.;
1793 for (G4int iFace=1; iFace<=nface; iFace++)
1794 G4int i0 = std::abs(pF[iFace].edge[0].v);
1795 G4int i1 = std::abs(pF[iFace].edge[1].v);
1796 G4int i2 = std::abs(pF[iFace].edge[2].v);
1797 G4int i3 = std::abs(pF[iFace].edge[3].v);
1798 if (i3 == 0) i3 = i0;
1799 srf += ((pV[i2] - pV[i0]).cross(pV[i3] -
1800 }
1801 return srf/2.;
1802 }
1803
1804 G4double HepPolyhedron::GetVolume() const
1805 /********************************************
1806 *
1807 * Name: HepPolyhedron::GetVolume
1808 * Author: E.Chernyaev
1809 *
1810 * Function: Returns volume of the polyhedron
1811 *
1812 ********************************************
1813 {
1814 G4double v = 0.;
1815 for (G4int iFace=1; iFace<=nface; iFace++)
1816 G4int i0 = std::abs(pF[iFace].edge[0].v);
1817 G4int i1 = std::abs(pF[iFace].edge[1].v);
1818 G4int i2 = std::abs(pF[iFace].edge[2].v);
1819 G4int i3 = std::abs(pF[iFace].edge[3].v);
1820 G4Point3D pt;
1821 if (i3 == 0) {
1822 i3 = i0;
1823 pt = (pV[i0]+pV[i1]+pV[i2]) * (1./3.);
1824 }else{
1825 pt = (pV[i0]+pV[i1]+pV[i2]+pV[i3]) * 0.
1826 }
1827 v += ((pV[i2] - pV[i0]).cross(pV[i3] - pV
1828 }
1829 return v/6.;
1830 }
1831
1832 G4int
1833 HepPolyhedron::createTwistedTrap(G4double Dz,
1834 const G4doub
1835 const G4doub
1836 /********************************************
1837 *
1838 * Name: createTwistedTrap
1839 * Author: E.Chernyaev (IHEP/Protvino)
1840 *
1841 * Function: Creates polyhedron for twisted t
1842 *
1843 * Input: Dz - half-length along Z
1844 * xy1[2,4] - quadrilateral at Z=-Dz
1845 * xy2[2,4] - quadrilateral at Z=+Dz
1846 *
1847 *
1848 ********************************************
1849 {
1850 AllocateMemory(12,18);
1851
1852 pV[ 1] = G4Point3D(xy1[0][0],xy1[0][1],-Dz)
1853 pV[ 2] = G4Point3D(xy1[1][0],xy1[1][1],-Dz)
1854 pV[ 3] = G4Point3D(xy1[2][0],xy1[2][1],-Dz)
1855 pV[ 4] = G4Point3D(xy1[3][0],xy1[3][1],-Dz)
1856
1857 pV[ 5] = G4Point3D(xy2[0][0],xy2[0][1], Dz)
1858 pV[ 6] = G4Point3D(xy2[1][0],xy2[1][1], Dz)
1859 pV[ 7] = G4Point3D(xy2[2][0],xy2[2][1], Dz)
1860 pV[ 8] = G4Point3D(xy2[3][0],xy2[3][1], Dz)
1861
1862 pV[ 9] = (pV[1]+pV[2]+pV[5]+pV[6])/4.;
1863 pV[10] = (pV[2]+pV[3]+pV[6]+pV[7])/4.;
1864 pV[11] = (pV[3]+pV[4]+pV[7]+pV[8])/4.;
1865 pV[12] = (pV[4]+pV[1]+pV[8]+pV[5])/4.;
1866
1867 enum {DUMMY, BOTTOM,
1868 LEFT_BOTTOM, LEFT_FRONT, LEFT_TOP,
1869 BACK_BOTTOM, BACK_LEFT, BACK_TOP,
1870 RIGHT_BOTTOM, RIGHT_BACK, RIGHT_TOP
1871 FRONT_BOTTOM, FRONT_RIGHT, FRONT_TOP
1872 TOP};
1873
1874 pF[ 1]=G4Facet(1,LEFT_BOTTOM, 4,BACK_BOTTOM
1875
1876 pF[ 2]=G4Facet(4,BOTTOM, -1,LEFT_FRONT,
1877 pF[ 3]=G4Facet(1,FRONT_LEFT, -5,LEFT_TOP,
1878 pF[ 4]=G4Facet(5,TOP, -8,LEFT_BACK,
1879 pF[ 5]=G4Facet(8,BACK_LEFT, -4,LEFT_BOTTOM
1880
1881 pF[ 6]=G4Facet(3,BOTTOM, -4,BACK_LEFT,
1882 pF[ 7]=G4Facet(4,LEFT_BACK, -8,BACK_TOP,
1883 pF[ 8]=G4Facet(8,TOP, -7,BACK_RIGHT,
1884 pF[ 9]=G4Facet(7,RIGHT_BACK, -3,BACK_BOTTOM
1885
1886 pF[10]=G4Facet(2,BOTTOM, -3,RIGHT_BACK,
1887 pF[11]=G4Facet(3,BACK_RIGHT, -7,RIGHT_TOP,
1888 pF[12]=G4Facet(7,TOP, -6,RIGHT_FRONT
1889 pF[13]=G4Facet(6,FRONT_RIGHT,-2,RIGHT_BOTTO
1890
1891 pF[14]=G4Facet(1,BOTTOM, -2,FRONT_RIGHT
1892 pF[15]=G4Facet(2,RIGHT_FRONT,-6,FRONT_TOP,
1893 pF[16]=G4Facet(6,TOP, -5,FRONT_LEFT,
1894 pF[17]=G4Facet(5,LEFT_FRONT, -1,FRONT_BOTTO
1895
1896 pF[18]=G4Facet(5,FRONT_TOP, 6,RIGHT_TOP, 7,
1897
1898 return 0;
1899 }
1900
1901 G4int
1902 HepPolyhedron::createPolyhedron(G4int Nnodes,
1903 const G4doubl
1904 const G4int
1905 /********************************************
1906 *
1907 * Name: createPolyhedron
1908 * Author: E.Chernyaev (IHEP/Protvino)
1909 *
1910 * Function: Creates user defined polyhedron
1911 *
1912 * Input: Nnodes - number of nodes
1913 * Nfaces - number of faces
1914 * nodes[][3] - node coordinates
1915 * faces[][4] - faces
1916 *
1917 ********************************************
1918 {
1919 AllocateMemory(Nnodes, Nfaces);
1920 if (nvert == 0) return 1;
1921
1922 for (G4int i=0; i<Nnodes; i++) {
1923 pV[i+1] = G4Point3D(xyz[i][0], xyz[i][1],
1924 }
1925 for (G4int k=0; k<Nfaces; k++) {
1926 pF[k+1] = G4Facet(faces[k][0],0,faces[k][
1927 }
1928 SetReferences();
1929 return 0;
1930 }
1931
1932 G4Point3D HepPolyhedron::vertexUnweightedMean
1933 /******************************************
1934 *
1935 * Name: vertexUnweightedMean
1936 * Author: S. Boogert (Manchester)
1937 *
1938 * Function: Calculate the unweighted mean
1939 * in the polyhedron. Not to be confused wi
1940 * centre of mass
1941 ******************************************
1942
1943 auto centre = G4Point3D();
1944 for(int i=1;i<=nvert;i++) {
1945 centre += pV[i];
1946 }
1947 centre = centre/nvert;
1948 return centre;
1949 }
1950
1951 HepPolyhedronTrd2::HepPolyhedronTrd2(G4double
1952 G4double
1953 G4double
1954 /********************************************
1955 *
1956 * Name: HepPolyhedronTrd2
1957 * Author: E.Chernyaev (IHEP/Protvino)
1958 *
1959 * Function: Create GEANT4 TRD2-trapezoid
1960 *
1961 * Input: Dx1 - half-length along X at -Dz
1962 * Dx2 - half-length along X ay +Dz
1963 * Dy1 - half-length along Y ay -Dz
1964 * Dy2 - half-length along Y ay +Dz
1965 * Dz - half-length along Z
1966 *
1967 ********************************************
1968 {
1969 AllocateMemory(8,6);
1970
1971 pV[1] = G4Point3D(-Dx1,-Dy1,-Dz);
1972 pV[2] = G4Point3D( Dx1,-Dy1,-Dz);
1973 pV[3] = G4Point3D( Dx1, Dy1,-Dz);
1974 pV[4] = G4Point3D(-Dx1, Dy1,-Dz);
1975 pV[5] = G4Point3D(-Dx2,-Dy2, Dz);
1976 pV[6] = G4Point3D( Dx2,-Dy2, Dz);
1977 pV[7] = G4Point3D( Dx2, Dy2, Dz);
1978 pV[8] = G4Point3D(-Dx2, Dy2, Dz);
1979
1980 CreatePrism();
1981 }
1982
1983 HepPolyhedronTrd2::~HepPolyhedronTrd2() = def
1984
1985 HepPolyhedronTrd1::HepPolyhedronTrd1(G4double
1986 G4double
1987 : HepPolyhedronTrd2(Dx1, Dx2, Dy, Dy, Dz) {
1988
1989 HepPolyhedronTrd1::~HepPolyhedronTrd1() = def
1990
1991 HepPolyhedronBox::HepPolyhedronBox(G4double D
1992 : HepPolyhedronTrd2(Dx, Dx, Dy, Dy, Dz) {}
1993
1994 HepPolyhedronBox::~HepPolyhedronBox() = defau
1995
1996 HepPolyhedronTrap::HepPolyhedronTrap(G4double
1997 G4double
1998 G4double
1999 G4double
2000 G4double
2001 G4double
2002 G4double
2003 G4double
2004 G4double
2005 G4double
2006 G4double
2007 /********************************************
2008 *
2009 * Name: HepPolyhedronTrap
2010 * Author: E.Chernyaev
2011 *
2012 * Function: Create GEANT4 TRAP-trapezoid
2013 *
2014 * Input: DZ - half-length in Z
2015 * Theta,Phi - polar angles of the lin
2016 * faces at Z=-Dz and Z=+D
2017 * Dy1 - half-length in Y of the face
2018 * Dx1 - half-length in X of low edge
2019 * Dx2 - half-length in X of top edge
2020 * Alp1 - angle between Y-axis and the
2021 * low edges of the face at Z=-
2022 * Dy2 - half-length in Y of the face
2023 * Dx3 - half-length in X of low edge
2024 * Dx4 - half-length in X of top edge
2025 * Alp2 - angle between Y-axis and the
2026 * low edges of the face at Z=+
2027 *
2028 ********************************************
2029 {
2030 G4double DzTthetaCphi = Dz*std::tan(Theta)*
2031 G4double DzTthetaSphi = Dz*std::tan(Theta)*
2032 G4double Dy1Talp1 = Dy1*std::tan(Alp1);
2033 G4double Dy2Talp2 = Dy2*std::tan(Alp2);
2034
2035 AllocateMemory(8,6);
2036
2037 pV[1] = G4Point3D(-DzTthetaCphi-Dy1Talp1-Dx
2038 pV[2] = G4Point3D(-DzTthetaCphi-Dy1Talp1+Dx
2039 pV[3] = G4Point3D(-DzTthetaCphi+Dy1Talp1+Dx
2040 pV[4] = G4Point3D(-DzTthetaCphi+Dy1Talp1-Dx
2041 pV[5] = G4Point3D( DzTthetaCphi-Dy2Talp2-Dx
2042 pV[6] = G4Point3D( DzTthetaCphi-Dy2Talp2+Dx
2043 pV[7] = G4Point3D( DzTthetaCphi+Dy2Talp2+Dx
2044 pV[8] = G4Point3D( DzTthetaCphi+Dy2Talp2-Dx
2045
2046 CreatePrism();
2047 }
2048
2049 HepPolyhedronTrap::~HepPolyhedronTrap() = def
2050
2051 HepPolyhedronPara::HepPolyhedronPara(G4double
2052 G4double
2053 G4double
2054 : HepPolyhedronTrap(Dz, Theta, Phi, Dy, Dx,
2055
2056 HepPolyhedronPara::~HepPolyhedronPara() = def
2057
2058 HepPolyhedronParaboloid::HepPolyhedronParabol
2059
2060
2061
2062
2063 /********************************************
2064 *
2065 * Name: HepPolyhedronParaboloid
2066 * Author: L.Lindroos, T.Nikitina (CERN), Jul
2067 *
2068 * Function: Constructor for paraboloid
2069 *
2070 * Input: r1 - inside and outside radiuses
2071 * r2 - inside and outside radiuses
2072 * dz - half length in Z
2073 * sPhi - starting angle of the segme
2074 * dPhi - segment range
2075 *
2076 ********************************************
2077 {
2078 static const G4double wholeCircle=twopi;
2079
2080 // C H E C K I N P U T P A R A M E T
2081
2082 G4int k = 0;
2083 if (r1 < 0. || r2 <= 0.) k = 1;
2084
2085 if (dz <= 0.) k += 2;
2086
2087 G4double phi1, phi2, dphi;
2088
2089 if(dPhi < 0.)
2090 {
2091 phi2 = sPhi; phi1 = phi2 + dPhi;
2092 }
2093 else if(dPhi == 0.)
2094 {
2095 phi1 = sPhi; phi2 = phi1 + wholeCircle;
2096 }
2097 else
2098 {
2099 phi1 = sPhi; phi2 = phi1 + dPhi;
2100 }
2101 dphi = phi2 - phi1;
2102
2103 if (std::abs(dphi-wholeCircle) < perMillion
2104 if (dphi > wholeCircle) k += 4;
2105
2106 if (k != 0) {
2107 std::cerr << "HepPolyhedronParaboloid: er
2108 if ((k & 1) != 0) std::cerr << " (radiuse
2109 if ((k & 2) != 0) std::cerr << " (half-le
2110 if ((k & 4) != 0) std::cerr << " (angles)
2111 std::cerr << std::endl;
2112 std::cerr << " r1=" << r1;
2113 std::cerr << " r2=" << r2;
2114 std::cerr << " dz=" << dz << " sPhi=" <<
2115 << std::endl;
2116 return;
2117 }
2118
2119 // P R E P A R E T W O P O L Y L I N
2120
2121 G4int n = GetNumberOfRotationSteps();
2122 G4double dl = (r2 - r1) / n;
2123 G4double k1 = (r2*r2 - r1*r1) / 2 / dz;
2124 G4double k2 = (r2*r2 + r1*r1) / 2;
2125
2126 auto zz = new G4double[n + 2], rr = new G4d
2127
2128 zz[0] = dz;
2129 rr[0] = r2;
2130
2131 for(G4int i = 1; i < n - 1; i++)
2132 {
2133 rr[i] = rr[i-1] - dl;
2134 zz[i] = (rr[i]*rr[i] - k2) / k1;
2135 if(rr[i] < 0)
2136 {
2137 rr[i] = 0;
2138 zz[i] = 0;
2139 }
2140 }
2141
2142 zz[n-1] = -dz;
2143 rr[n-1] = r1;
2144
2145 zz[n] = dz;
2146 rr[n] = 0;
2147
2148 zz[n+1] = -dz;
2149 rr[n+1] = 0;
2150
2151 // R O T A T E P O L Y L I N E S
2152
2153 RotateAroundZ(0, phi1, dphi, n, 2, zz, rr,
2154 SetReferences();
2155
2156 delete [] zz;
2157 delete [] rr;
2158 }
2159
2160 HepPolyhedronParaboloid::~HepPolyhedronParabo
2161
2162 HepPolyhedronHype::HepPolyhedronHype(G4double
2163 G4double
2164 G4double
2165 G4double
2166 G4double
2167 /********************************************
2168 *
2169 * Name: HepPolyhedronHype
2170 * Author: Tatiana Nikitina (CERN)
2171 * Evgueni Tcherniaev
2172 *
2173 * Function: Constructor for Hype
2174 *
2175 * Input: r1 - inside radius at z=0
2176 * r2 - outside radiuses at z=0
2177 * sqrtan1 - sqr of tan of Inner Ster
2178 * sqrtan2 - sqr of tan of Outer Ster
2179 * halfZ - half length in Z
2180 *
2181 ********************************************
2182 {
2183 static const G4double wholeCircle = twopi;
2184
2185 // C H E C K I N P U T P A R A M E T
2186
2187 G4int k = 0;
2188 if (r1 < 0. || r2 < 0. || r1 >= r2) k = 1;
2189 if (halfZ <= 0.) k += 2;
2190 if (sqrtan1 < 0.|| sqrtan2 < 0.) k += 4;
2191
2192 if (k != 0)
2193 {
2194 std::cerr << "HepPolyhedronHype: error in
2195 if ((k & 1) != 0) std::cerr << " (radiuse
2196 if ((k & 2) != 0) std::cerr << " (half-le
2197 if ((k & 4) != 0) std::cerr << " (angles)
2198 std::cerr << std::endl;
2199 std::cerr << " r1=" << r1 << " r2=" << r2
2200 std::cerr << " halfZ=" << halfZ << " sqrT
2201 << " sqrTan2=" << sqrtan2
2202 << std::endl;
2203 return;
2204 }
2205
2206 // P R E P A R E T W O P O L Y L I N
2207
2208 G4int ns = std::max(3, GetNumberOfRotationS
2209 G4int nz1 = (sqrtan1 == 0.) ? 2 : ns + 1;
2210 G4int nz2 = (sqrtan2 == 0.) ? 2 : ns + 1;
2211 auto zz = new G4double[nz1 + nz2];
2212 auto rr = new G4double[nz1 + nz2];
2213
2214 // external polyline
2215 G4double dz2 = 2.*halfZ/(nz2 - 1);
2216 for(G4int i = 0; i < nz2; ++i)
2217 {
2218 zz[i] = halfZ - dz2*i;
2219 rr[i] = std::sqrt(sqrtan2*zz[i]*zz[i] + r
2220 }
2221
2222 // internal polyline
2223 G4double dz1 = 2.*halfZ/(nz1 - 1);
2224 for(G4int i = 0; i < nz1; ++i)
2225 {
2226 G4int j = nz2 + i;
2227 zz[j] = halfZ - dz1*i;
2228 rr[j] = std::sqrt(sqrtan1*zz[j]*zz[j] + r
2229 }
2230
2231 // R O T A T E P O L Y L I N E S
2232
2233 RotateAroundZ(0, 0., wholeCircle, nz2, nz1,
2234 SetReferences();
2235
2236 delete [] zz;
2237 delete [] rr;
2238 }
2239
2240 HepPolyhedronHype::~HepPolyhedronHype() = def
2241
2242 HepPolyhedronCons::HepPolyhedronCons(G4double
2243 G4double
2244 G4double
2245 G4double
2246 G4double
2247 G4double
2248 G4double
2249 /********************************************
2250 *
2251 * Name: HepPolyhedronCons::HepPolyhedronCons
2252 * Author: E.Chernyaev (IHEP/Protvino)
2253 *
2254 * Function: Constructor for CONS, TUBS, CONE
2255 *
2256 * Input: Rmn1, Rmx1 - inside and outside rad
2257 * Rmn2, Rmx2 - inside and outside rad
2258 * Dz - half length in Z
2259 * Phi1 - starting angle of the
2260 * Dphi - segment range
2261 *
2262 ********************************************
2263 {
2264 static const G4double wholeCircle=twopi;
2265
2266 // C H E C K I N P U T P A R A M E T
2267
2268 G4int k = 0;
2269 if (Rmn1 < 0. || Rmx1 < 0. || Rmn2 < 0. ||
2270 if (Rmn1 > Rmx1 || Rmn2 > Rmx2)
2271 if (Rmn1 == Rmx1 && Rmn2 == Rmx2)
2272
2273 if (Dz <= 0.) k += 2;
2274
2275 G4double phi1, phi2, dphi;
2276 if (Dphi < 0.) {
2277 phi2 = Phi1; phi1 = phi2 - Dphi;
2278 }else if (Dphi == 0.) {
2279 phi1 = Phi1; phi2 = phi1 + wholeCircle;
2280 }else{
2281 phi1 = Phi1; phi2 = phi1 + Dphi;
2282 }
2283 dphi = phi2 - phi1;
2284 if (std::abs(dphi-wholeCircle) < perMillion
2285 if (dphi > wholeCircle) k += 4;
2286
2287 if (k != 0) {
2288 std::cerr << "HepPolyhedronCone(s)/Tube(s
2289 if ((k & 1) != 0) std::cerr << " (radiuse
2290 if ((k & 2) != 0) std::cerr << " (half-le
2291 if ((k & 4) != 0) std::cerr << " (angles)
2292 std::cerr << std::endl;
2293 std::cerr << " Rmn1=" << Rmn1 << " Rmx1="
2294 std::cerr << " Rmn2=" << Rmn2 << " Rmx2="
2295 std::cerr << " Dz=" << Dz << " Phi1=" <<
2296 << std::endl;
2297 return;
2298 }
2299
2300 // P R E P A R E T W O P O L Y L I N
2301
2302 G4double zz[4], rr[4];
2303 zz[0] = Dz;
2304 zz[1] = -Dz;
2305 zz[2] = Dz;
2306 zz[3] = -Dz;
2307 rr[0] = Rmx2;
2308 rr[1] = Rmx1;
2309 rr[2] = Rmn2;
2310 rr[3] = Rmn1;
2311
2312 // R O T A T E P O L Y L I N E S
2313
2314 RotateAroundZ(0, phi1, dphi, 2, 2, zz, rr,
2315 SetReferences();
2316 }
2317
2318 HepPolyhedronCons::~HepPolyhedronCons() = def
2319
2320 HepPolyhedronCone::HepPolyhedronCone(G4double
2321 G4double
2322 G4double
2323 HepPolyhedronCons(Rmn1, Rmx1, Rmn2, Rmx2, D
2324
2325 HepPolyhedronCone::~HepPolyhedronCone() = def
2326
2327 HepPolyhedronTubs::HepPolyhedronTubs(G4double
2328 G4double
2329 G4double
2330 : HepPolyhedronCons(Rmin, Rmax, Rmin, Rma
2331
2332 HepPolyhedronTubs::~HepPolyhedronTubs() = def
2333
2334 HepPolyhedronTube::HepPolyhedronTube (G4doubl
2335 G4doubl
2336 : HepPolyhedronCons(Rmin, Rmax, Rmin, Rmax,
2337
2338 HepPolyhedronTube::~HepPolyhedronTube () = de
2339
2340 HepPolyhedronPgon::HepPolyhedronPgon(G4double
2341 G4double
2342 G4int np
2343 G4int nz
2344 const G4
2345 const G4
2346 const G4
2347 /********************************************
2348 *
2349 * Name: HepPolyhedronPgon
2350 * Author: E.Chernyaev
2351 *
2352 * Function: Constructor of polyhedron for PG
2353 *
2354 * Input: phi - initial phi
2355 * dphi - delta phi
2356 * npdv - number of steps along phi
2357 * nz - number of z-planes (at least
2358 * z[] - z coordinates of the slices
2359 * rmin[] - smaller r at the slices
2360 * rmax[] - bigger r at the slices
2361 *
2362 ********************************************
2363 {
2364 // C H E C K I N P U T P A R A M E T
2365
2366 if (dphi <= 0. || dphi > twopi) {
2367 std::cerr
2368 << "HepPolyhedronPgon/Pcon: wrong delta
2369 << std::endl;
2370 return;
2371 }
2372
2373 if (nz < 2) {
2374 std::cerr
2375 << "HepPolyhedronPgon/Pcon: number of z
2376 << std::endl;
2377 return;
2378 }
2379
2380 if (npdv < 0) {
2381 std::cerr
2382 << "HepPolyhedronPgon/Pcon: error in nu
2383 << std::endl;
2384 return;
2385 }
2386
2387 G4int i;
2388 for (i=0; i<nz; i++) {
2389 if (rmin[i] < 0. || rmax[i] < 0. || rmin[
2390 std::cerr
2391 << "HepPolyhedronPgon: error in radiu
2392 << rmin[i] << " rmax[" << i << "]=" <
2393 << std::endl;
2394 return;
2395 }
2396 }
2397
2398 // P R E P A R E T W O P O L Y L I N
2399
2400 G4double *zz, *rr;
2401 zz = new G4double[2*nz];
2402 rr = new G4double[2*nz];
2403
2404 if (z[0] > z[nz-1]) {
2405 for (i=0; i<nz; i++) {
2406 zz[i] = z[i];
2407 rr[i] = rmax[i];
2408 zz[i+nz] = z[i];
2409 rr[i+nz] = rmin[i];
2410 }
2411 }else{
2412 for (i=0; i<nz; i++) {
2413 zz[i] = z[nz-i-1];
2414 rr[i] = rmax[nz-i-1];
2415 zz[i+nz] = z[nz-i-1];
2416 rr[i+nz] = rmin[nz-i-1];
2417 }
2418 }
2419
2420 // R O T A T E P O L Y L I N E S
2421
2422 G4int nodeVis = 1;
2423 G4int edgeVis = (npdv == 0) ? -1 : 1;
2424 RotateAroundZ(npdv, phi, dphi, nz, nz, zz,
2425 SetReferences();
2426
2427 delete [] zz;
2428 delete [] rr;
2429 }
2430
2431 HepPolyhedronPgon::HepPolyhedronPgon(G4double
2432 G4double
2433 G4int np
2434 const st
2435 /********************************************
2436 *
2437 * Name: HepPolyhedronPgon
2438 * Author: E.Tcherniaev (E.Chernyaev)
2439 *
2440 * Function: Constructor of polyhedron for PG
2441 *
2442 * Input: phi - initial phi
2443 * dphi - delta phi
2444 * npdv - number of steps along phi
2445 * rz - rz-contour
2446 *
2447 ********************************************
2448 {
2449 // C H E C K I N P U T P A R A M E T
2450
2451 if (dphi <= 0. || dphi > twopi) {
2452 std::cerr
2453 << "HepPolyhedronPgon/Pcon: wrong delta
2454 << std::endl;
2455 return;
2456 }
2457
2458 if (npdv < 0) {
2459 std::cerr
2460 << "HepPolyhedronPgon/Pcon: error in nu
2461 << std::endl;
2462 return;
2463 }
2464
2465 G4int nrz = (G4int)rz.size();
2466 if (nrz < 3) {
2467 std::cerr
2468 << "HepPolyhedronPgon/Pcon: invalid num
2469 << std::endl;
2470 return;
2471 }
2472
2473 // R O T A T E P O L Y L I N E
2474
2475 G4int nodeVis = 1;
2476 G4int edgeVis = (npdv == 0) ? -1 : 1;
2477 RotateContourAroundZ(npdv, phi, dphi, rz, n
2478 SetReferences();
2479 }
2480
2481 HepPolyhedronPgon::~HepPolyhedronPgon() = def
2482
2483 HepPolyhedronPcon::HepPolyhedronPcon(G4double
2484 const G4
2485 const G4
2486 const G4
2487 : HepPolyhedronPgon(phi, dphi, 0, nz, z, rm
2488
2489 HepPolyhedronPcon::HepPolyhedronPcon(G4double
2490 const st
2491 : HepPolyhedronPgon(phi, dphi, 0, rz) {}
2492
2493 HepPolyhedronPcon::~HepPolyhedronPcon() = def
2494
2495 HepPolyhedronSphere::HepPolyhedronSphere(G4do
2496 G4do
2497 G4do
2498 /********************************************
2499 *
2500 * Name: HepPolyhedronSphere
2501 * Author: E.Chernyaev (IHEP/Protvino)
2502 *
2503 * Function: Constructor of polyhedron for SP
2504 *
2505 * Input: rmin - internal radius
2506 * rmax - external radius
2507 * phi - initial phi
2508 * dphi - delta phi
2509 * the - initial theta
2510 * dthe - delta theta
2511 *
2512 ********************************************
2513 {
2514 // C H E C K I N P U T P A R A M E T
2515
2516 if (dphi <= 0. || dphi > twopi) {
2517 std::cerr
2518 << "HepPolyhedronSphere: wrong delta ph
2519 << std::endl;
2520 return;
2521 }
2522
2523 if (the < 0. || the > pi) {
2524 std::cerr
2525 << "HepPolyhedronSphere: wrong theta =
2526 << std::endl;
2527 return;
2528 }
2529
2530 if (dthe <= 0. || dthe > pi) {
2531 std::cerr
2532 << "HepPolyhedronSphere: wrong delta th
2533 << std::endl;
2534 return;
2535 }
2536
2537 if (the+dthe > pi) {
2538 std::cerr
2539 << "HepPolyhedronSphere: wrong theta +
2540 << the << " " << dthe
2541 << std::endl;
2542 return;
2543 }
2544
2545 if (rmin < 0. || rmin >= rmax) {
2546 std::cerr
2547 << "HepPolyhedronSphere: error in radiu
2548 << " rmin=" << rmin << " rmax=" << rmax
2549 << std::endl;
2550 return;
2551 }
2552
2553 // P R E P A R E T W O P O L Y L I N
2554
2555 G4int nds = (GetNumberOfRotationSteps() + 1
2556 G4int np1 = G4int(dthe*nds/pi+.5) + 1;
2557 if (np1 <= 1) np1 = 2;
2558 G4int np2 = rmin < spatialTolerance ? 1 : n
2559
2560 G4double *zz, *rr;
2561 zz = new G4double[np1+np2];
2562 rr = new G4double[np1+np2];
2563
2564 G4double a = dthe/(np1-1);
2565 G4double cosa, sina;
2566 for (G4int i=0; i<np1; i++) {
2567 cosa = std::cos(the+i*a);
2568 sina = std::sin(the+i*a);
2569 zz[i] = rmax*cosa;
2570 rr[i] = rmax*sina;
2571 if (np2 > 1) {
2572 zz[i+np1] = rmin*cosa;
2573 rr[i+np1] = rmin*sina;
2574 }
2575 }
2576 if (np2 == 1) {
2577 zz[np1] = 0.;
2578 rr[np1] = 0.;
2579 }
2580
2581 // R O T A T E P O L Y L I N E S
2582
2583 RotateAroundZ(0, phi, dphi, np1, np2, zz, r
2584 SetReferences();
2585
2586 delete [] zz;
2587 delete [] rr;
2588 }
2589
2590 HepPolyhedronSphere::~HepPolyhedronSphere() =
2591
2592 HepPolyhedronTorus::HepPolyhedronTorus(G4doub
2593 G4doub
2594 G4doub
2595 G4doub
2596 G4doub
2597 /********************************************
2598 *
2599 * Name: HepPolyhedronTorus
2600 * Author: E.Chernyaev (IHEP/Protvino)
2601 *
2602 * Function: Constructor of polyhedron for TO
2603 *
2604 * Input: rmin - internal radius
2605 * rmax - external radius
2606 * rtor - radius of torus
2607 * phi - initial phi
2608 * dphi - delta phi
2609 *
2610 ********************************************
2611 {
2612 // C H E C K I N P U T P A R A M E T
2613
2614 if (dphi <= 0. || dphi > twopi) {
2615 std::cerr
2616 << "HepPolyhedronTorus: wrong delta phi
2617 << std::endl;
2618 return;
2619 }
2620
2621 if (rmin < 0. || rmin >= rmax || rmax >= rt
2622 std::cerr
2623 << "HepPolyhedronTorus: error in radius
2624 << " rmin=" << rmin << " rmax=" << rmax
2625 << std::endl;
2626 return;
2627 }
2628
2629 // P R E P A R E T W O P O L Y L I N
2630
2631 G4int np1 = GetNumberOfRotationSteps();
2632 G4int np2 = rmin < spatialTolerance ? 1 : n
2633
2634 G4double *zz, *rr;
2635 zz = new G4double[np1+np2];
2636 rr = new G4double[np1+np2];
2637
2638 G4double a = twopi/np1;
2639 G4double cosa, sina;
2640 for (G4int i=0; i<np1; i++) {
2641 cosa = std::cos(i*a);
2642 sina = std::sin(i*a);
2643 zz[i] = rmax*cosa;
2644 rr[i] = rtor+rmax*sina;
2645 if (np2 > 1) {
2646 zz[i+np1] = rmin*cosa;
2647 rr[i+np1] = rtor+rmin*sina;
2648 }
2649 }
2650 if (np2 == 1) {
2651 zz[np1] = 0.;
2652 rr[np1] = rtor;
2653 np2 = -1;
2654 }
2655
2656 // R O T A T E P O L Y L I N E S
2657
2658 RotateAroundZ(0, phi, dphi, -np1, -np2, zz,
2659 SetReferences();
2660
2661 delete [] zz;
2662 delete [] rr;
2663 }
2664
2665 HepPolyhedronTorus::~HepPolyhedronTorus() = d
2666
2667 HepPolyhedronTet::HepPolyhedronTet(const G4do
2668 const G4do
2669 const G4do
2670 const G4do
2671 /********************************************
2672 *
2673 * Name: HepPolyhedronTet
2674 * Author: E.Tcherniaev (E.Chernyaev)
2675 *
2676 * Function: Constructor of polyhedron for TE
2677 *
2678 * Input: p0,p1,p2,p3 - vertices
2679 *
2680 ********************************************
2681 {
2682 AllocateMemory(4,4);
2683
2684 pV[1].set(p0[0], p0[1], p0[2]);
2685 pV[2].set(p1[0], p1[1], p1[2]);
2686 pV[3].set(p2[0], p2[1], p2[2]);
2687 pV[4].set(p3[0], p3[1], p3[2]);
2688
2689 G4Vector3D v1(pV[2] - pV[1]);
2690 G4Vector3D v2(pV[3] - pV[1]);
2691 G4Vector3D v3(pV[4] - pV[1]);
2692
2693 if (v1.cross(v2).dot(v3) < 0.)
2694 {
2695 pV[3].set(p3[0], p3[1], p3[2]);
2696 pV[4].set(p2[0], p2[1], p2[2]);
2697 }
2698
2699 pF[1] = G4Facet(1,2, 3,4, 2,3);
2700 pF[2] = G4Facet(1,3, 4,4, 3,1);
2701 pF[3] = G4Facet(1,1, 2,4, 4,2);
2702 pF[4] = G4Facet(2,1, 3,2, 4,3);
2703 }
2704
2705 HepPolyhedronTet::~HepPolyhedronTet() = defau
2706
2707 HepPolyhedronEllipsoid::HepPolyhedronEllipsoi
2708
2709
2710 /********************************************
2711 *
2712 * Name: HepPolyhedronEllipsoid
2713 * Author: G.Guerrieri
2714 * Evgueni Tcherniaev
2715 *
2716 * Function: Constructor of polyhedron for EL
2717 *
2718 * Input: ax - semiaxis x
2719 * by - semiaxis y
2720 * cz - semiaxis z
2721 * zCut1 - lower cut plane level (soli
2722 * zCut2 - upper cut plane level (soli
2723 *
2724 ********************************************
2725 {
2726 // C H E C K I N P U T P A R A M E T
2727
2728 if (zCut1 >= cz || zCut2 <= -cz || zCut1 >
2729 std::cerr << "HepPolyhedronEllipsoid: wro
2730 << " zCut2 = " << zCut2
2731 << " for given cz = " << cz << std
2732 return;
2733 }
2734 if (cz <= 0.0) {
2735 std::cerr << "HepPolyhedronEllipsoid: bad
2736 << std::endl;
2737 return;
2738 }
2739
2740 // P R E P A R E T W O P O L Y L I N
2741 // generate sphere of radius cz first, th
2742
2743 G4double sthe = std::acos(zCut2/cz);
2744 G4double dthe = std::acos(zCut1/cz) - sthe;
2745 G4int nds = (GetNumberOfRotationSteps() + 1
2746 G4int np1 = G4int(dthe*nds/pi + 0.5) + 1;
2747 if (np1 <= 1) np1 = 2;
2748 G4int np2 = 2;
2749
2750 G4double *zz, *rr;
2751 zz = new G4double[np1 + np2];
2752 rr = new G4double[np1 + np2];
2753 if ((zz == nullptr) || (rr == nullptr))
2754 {
2755 G4Exception("HepPolyhedronEllipsoid::HepP
2756 "greps1002", FatalException,
2757 }
2758
2759 G4double a = dthe/(np1 - 1);
2760 G4double cosa, sina;
2761 for (G4int i = 0; i < np1; ++i)
2762 {
2763 cosa = std::cos(sthe + i*a);
2764 sina = std::sin(sthe + i*a);
2765 zz[i] = cz*cosa;
2766 rr[i] = cz*sina;
2767 }
2768 zz[np1 + 0] = zCut2;
2769 rr[np1 + 0] = 0.;
2770 zz[np1 + 1] = zCut1;
2771 rr[np1 + 1] = 0.;
2772
2773 // R O T A T E P O L Y L I N E S
2774
2775 RotateAroundZ(0, 0., twopi, np1, np2, zz, r
2776 SetReferences();
2777
2778 delete [] zz;
2779 delete [] rr;
2780
2781 // rescale x and y vertex coordinates
2782 G4double kx = ax/cz;
2783 G4double ky = by/cz;
2784 G4Point3D* p = pV;
2785 for (G4int i = 0; i < nvert; ++i, ++p)
2786 {
2787 p->setX(p->x()*kx);
2788 p->setY(p->y()*ky);
2789 }
2790 }
2791
2792 HepPolyhedronEllipsoid::~HepPolyhedronEllipso
2793
2794 HepPolyhedronEllipticalCone::HepPolyhedronEll
2795
2796
2797
2798 /********************************************
2799 *
2800 * Name: HepPolyhedronEllipticalCone
2801 * Author: D.Anninos
2802 *
2803 * Function: Constructor for EllipticalCone
2804 *
2805 * Input: ax, ay - X & Y semi axes at z =
2806 * h - height of full cone
2807 * zTopCut - Top Cut in Z Axis
2808 *
2809 ********************************************
2810 {
2811 // C H E C K I N P U T P A R A M E T
2812
2813 G4int k = 0;
2814 if ( (ax <= 0.) || (ay <= 0.) || (h <= 0.)
2815
2816 if (k != 0) {
2817 std::cerr << "HepPolyhedronCone: error in
2818 std::cerr << std::endl;
2819 return;
2820 }
2821
2822 // P R E P A R E T W O P O L Y L I N
2823
2824 zTopCut = (h >= zTopCut ? zTopCut : h);
2825
2826 G4double *zz, *rr;
2827 zz = new G4double[4];
2828 rr = new G4double[4];
2829 zz[0] = zTopCut;
2830 zz[1] = -zTopCut;
2831 zz[2] = zTopCut;
2832 zz[3] = -zTopCut;
2833 rr[0] = (h-zTopCut);
2834 rr[1] = (h+zTopCut);
2835 rr[2] = 0.;
2836 rr[3] = 0.;
2837
2838 // R O T A T E P O L Y L I N E S
2839
2840 RotateAroundZ(0, 0., twopi, 2, 2, zz, rr, -
2841 SetReferences();
2842
2843 delete [] zz;
2844 delete [] rr;
2845
2846 // rescale x and y vertex coordinates
2847 {
2848 G4Point3D * p= pV;
2849 for (G4int i=0; i<nvert; i++, p++) {
2850 p->setX( p->x() * ax );
2851 p->setY( p->y() * ay );
2852 }
2853 }
2854 }
2855
2856 HepPolyhedronEllipticalCone::~HepPolyhedronEl
2857
2858 HepPolyhedronHyperbolicMirror::HepPolyhedronH
2859
2860
2861 /********************************************
2862 *
2863 * Name: HepPolyhedronHyperbolicMirror
2864 * Author: E.Tcherniaev (E.Chernyaev)
2865 *
2866 * Function: Create polyhedron for Hyperbolic
2867 *
2868 * Input: a - half-separation
2869 * h - height
2870 * r - radius
2871 *
2872 ********************************************
2873 {
2874 G4double H = std::abs(h);
2875 G4double R = std::abs(r);
2876 G4double A = std::abs(a);
2877 G4double B = A*R/std::sqrt(2*A*H + H*H);
2878
2879 // P R E P A R E T W O P O L Y L I N
2880
2881 G4int np1 = (A == 0.) ? 2 : std::max(3, Get
2882 G4int np2 = 2;
2883 G4double maxAng = (A == 0.) ? 0. : std::aco
2884 G4double delAng = maxAng/(np1 - 1);
2885
2886 auto zz = new G4double[np1 + np2];
2887 auto rr = new G4double[np1 + np2];
2888
2889 // 1st polyline
2890 zz[0] = H;
2891 rr[0] = R;
2892 for (G4int iz = 1; iz < np1 - 1; ++iz)
2893 {
2894 G4double ang = maxAng - iz*delAng;
2895 zz[iz] = A*std::cosh(ang) - A;
2896 rr[iz] = B*std::sinh(ang);
2897 }
2898 zz[np1 - 1] = 0.;
2899 rr[np1 - 1] = 0.;
2900
2901 // 2nd polyline
2902 zz[np1] = H;
2903 rr[np1] = 0.;
2904 zz[np1 + 1] = 0.;
2905 rr[np1 + 1] = 0.;
2906
2907 // R O T A T E P O L Y L I N E S
2908
2909 G4double phi = 0.;
2910 G4double dphi = CLHEP::twopi;
2911 RotateAroundZ(0, phi, dphi, np1, np2, zz, r
2912 SetReferences();
2913
2914 delete [] zz;
2915 delete [] rr;
2916 }
2917
2918 HepPolyhedronHyperbolicMirror::~HepPolyhedron
2919
2920 HepPolyhedronTetMesh::
2921 HepPolyhedronTetMesh(const std::vector<G4Thre
2922 /********************************************
2923 *
2924 * Name: HepPolyhedronTetMesh
2925 * Author: E.Tcherniaev (E.Chernyaev)
2926 *
2927 * Function: Create polyhedron for tetrahedro
2928 *
2929 * Input: tetrahedra - array of tetrahedron v
2930 * per tetrahedron
2931 *
2932 ********************************************
2933 {
2934 // Check size of input vector
2935 G4int nnodes = (G4int)tetrahedra.size();
2936 if (nnodes == 0)
2937 {
2938 std::cerr
2939 << "HepPolyhedronTetMesh: Empty tetrahe
2940 return;
2941 }
2942 G4int ntet = nnodes/4;
2943 if (nnodes != ntet*4)
2944 {
2945 std::cerr << "HepPolyhedronTetMesh: Numbe
2946 << " in tetrahedron mesh is NOT
2947 << std::endl;
2948 return;
2949 }
2950
2951 // Find coincident vertices using hash tabl
2952 // This could be done using std::unordered_
2953 // below runs faster.
2954 std::vector<G4int> iheads(nnodes, -1);
2955 std::vector<std::pair<G4int,G4int>> ipairs(
2956 for (G4int i = 0; i < nnodes; ++i)
2957 {
2958 // Generate hash key
2959 G4ThreeVector point = tetrahedra[i];
2960 auto key = std::hash<G4double>()(point.x(
2961 key ^= std::hash<G4double>()(point.y());
2962 key ^= std::hash<G4double>()(point.z());
2963 key %= nnodes;
2964 // Check head of the list
2965 if (iheads[key] < 0)
2966 {
2967 iheads[key] = i;
2968 ipairs[i].first = i;
2969 continue;
2970 }
2971 // Loop along the list
2972 for (G4int icur = iheads[key], iprev = 0;
2973 {
2974 G4int icheck = ipairs[icur].first;
2975 if (tetrahedra[icheck] == point)
2976 {
2977 ipairs[i].first = icheck; // coincide
2978 break;
2979 }
2980 iprev = icur;
2981 icur = ipairs[icur].second;
2982 // Append vertex to the list
2983 if (icur < 0)
2984 {
2985 ipairs[i].first = i;
2986 ipairs[iprev].second = i;
2987 break;
2988 }
2989 }
2990 }
2991
2992 // Create vector of original facets
2993 struct facet
2994 {
2995 G4int i1, i2, i3;
2996 facet() : i1(0), i2(0), i3(0) {};
2997 facet(G4int k1, G4int k2, G4int k3) : i1(
2998 };
2999 G4int nfacets = nnodes;
3000 std::vector<facet> ifacets(nfacets);
3001 for (G4int i = 0; i < nfacets; i += 4)
3002 {
3003 G4int i0 = ipairs[i + 0].first;
3004 G4int i1 = ipairs[i + 1].first;
3005 G4int i2 = ipairs[i + 2].first;
3006 G4int i3 = ipairs[i + 3].first;
3007 if (i0 > i1) std::swap(i0, i1);
3008 if (i0 > i2) std::swap(i0, i2);
3009 if (i0 > i3) std::swap(i0, i3);
3010 if (i1 > i2) std::swap(i1, i2);
3011 if (i1 > i3) std::swap(i1, i3);
3012 G4ThreeVector e1 = tetrahedra[i1] - tetra
3013 G4ThreeVector e2 = tetrahedra[i2] - tetra
3014 G4ThreeVector e3 = tetrahedra[i3] - tetra
3015 G4double volume = (e1.cross(e2)).dot(e3);
3016 if (volume > 0.) std::swap(i2, i3);
3017 ifacets[i + 0] = facet(i0, i1, i2);
3018 ifacets[i + 1] = facet(i0, i2, i3);
3019 ifacets[i + 2] = facet(i0, i3, i1);
3020 ifacets[i + 3] = facet(i1, i3, i2);
3021 }
3022
3023 // Find shared facets
3024 std::fill(iheads.begin(), iheads.end(), -1)
3025 std::fill(ipairs.begin(), ipairs.end(), std
3026 for (G4int i = 0; i < nfacets; ++i)
3027 {
3028 // Check head of the list
3029 G4int key = ifacets[i].i1;
3030 if (iheads[key] < 0)
3031 {
3032 iheads[key] = i;
3033 ipairs[i].first = i;
3034 continue;
3035 }
3036 // Loop along the list
3037 G4int i2 = ifacets[i].i2, i3 = ifacets[i]
3038 for (G4int icur = iheads[key], iprev = -1
3039 {
3040 G4int icheck = ipairs[icur].first;
3041 if (ifacets[icheck].i2 == i3 && ifacets
3042 {
3043 if (iprev < 0)
3044 {
3045 iheads[key] = ipairs[icur].second;
3046 }
3047 else
3048 {
3049 ipairs[iprev].second = ipairs[icur]
3050 }
3051 ipairs[icur].first = -1; // shared fa
3052 ipairs[icur].second = -1;
3053 break;
3054 }
3055 iprev = icur;
3056 icur = ipairs[icur].second;
3057 // Append facet to the list
3058 if (icur < 0)
3059 {
3060 ipairs[i].first = i;
3061 ipairs[iprev].second = i;
3062 break;
3063 }
3064 }
3065 }
3066
3067 // Count vertices and facets skipping share
3068 std::fill(iheads.begin(), iheads.end(), -1)
3069 G4int nver = 0, nfac = 0;
3070 for (G4int i = 0; i < nfacets; ++i)
3071 {
3072 if (ipairs[i].first < 0) continue;
3073 G4int i1 = ifacets[i].i1;
3074 G4int i2 = ifacets[i].i2;
3075 G4int i3 = ifacets[i].i3;
3076 if (iheads[i1] < 0) iheads[i1] = nver++;
3077 if (iheads[i2] < 0) iheads[i2] = nver++;
3078 if (iheads[i3] < 0) iheads[i3] = nver++;
3079 nfac++;
3080 }
3081
3082 // Construct polyhedron
3083 AllocateMemory(nver, nfac);
3084 for (G4int i = 0; i < nnodes; ++i)
3085 {
3086 G4int k = iheads[i];
3087 if (k >= 0) SetVertex(k + 1, tetrahedra[i
3088 }
3089 for (G4int i = 0, k = 0; i < nfacets; ++i)
3090 {
3091 if (ipairs[i].first < 0) continue;
3092 G4int i1 = iheads[ifacets[i].i1] + 1;
3093 G4int i2 = iheads[ifacets[i].i2] + 1;
3094 G4int i3 = iheads[ifacets[i].i3] + 1;
3095 SetFacet(++k, i1, i2, i3);
3096 }
3097 SetReferences();
3098 }
3099
3100 HepPolyhedronTetMesh::~HepPolyhedronTetMesh()
3101
3102 HepPolyhedronBoxMesh::
3103 HepPolyhedronBoxMesh(G4double sizeX, G4double
3104 const std::vector<G4Thre
3105 /********************************************
3106 *
3107 * Name: HepPolyhedronBoxMesh
3108 * Author: E.Tcherniaev (E.Chernyaev)
3109 *
3110 * Function: Create polyhedron for box mesh
3111 *
3112 * Input: sizeX, sizeY, sizeZ - dimensions of
3113 * positions - vector of cell centres
3114 *
3115 ********************************************
3116 {
3117 G4int nbox = (G4int)positions.size();
3118 if (nbox == 0)
3119 {
3120 std::cerr << "HepPolyhedronBoxMesh: Empty
3121 return;
3122 }
3123 // compute inverse dimensions
3124 G4double invx = 1./sizeX, invy = 1./sizeY,
3125 // find mesh bounding box
3126 G4ThreeVector pmin = positions[0], pmax = p
3127 for (const auto& p: positions)
3128 {
3129 if (pmin.x() > p.x()) pmin.setX(p.x());
3130 if (pmin.y() > p.y()) pmin.setY(p.y());
3131 if (pmin.z() > p.z()) pmin.setZ(p.z());
3132 if (pmax.x() < p.x()) pmax.setX(p.x());
3133 if (pmax.y() < p.y()) pmax.setY(p.y());
3134 if (pmax.z() < p.z()) pmax.setZ(p.z());
3135 }
3136 // find number of voxels
3137 G4int nx = (pmax.x() - pmin.x())*invx + 1.5
3138 G4int ny = (pmax.y() - pmin.y())*invy + 1.5
3139 G4int nz = (pmax.z() - pmin.z())*invz + 1.5
3140 // create structures for voxels and node in
3141 std::vector<char> voxels(nx*ny*nz, 0);
3142 std::vector<G4int> indices((nx+1)*(ny+1)*(n
3143 // mark voxels listed in positions
3144 G4int kx = ny*nz, ky = nz;
3145 for (const auto& p: positions)
3146 {
3147 G4int ix = (p.x() - pmin.x())*invx + 0.5;
3148 G4int iy = (p.y() - pmin.y())*invy + 0.5;
3149 G4int iz = (p.z() - pmin.z())*invz + 0.5;
3150 G4int i = ix*kx + iy*ky + iz;
3151 voxels[i] = 1;
3152 }
3153 // count number of vertices and facets
3154 // set indices
3155 G4int kvx = (ny + 1)*(nz + 1), kvy = nz + 1
3156 G4int nver = 0, nfac = 0;
3157 for (const auto& p: positions)
3158 {
3159 G4int ix = (p.x() - pmin.x())*invx + 0.5;
3160 G4int iy = (p.y() - pmin.y())*invy + 0.5;
3161 G4int iz = (p.z() - pmin.z())*invz + 0.5;
3162 //
3163 // 011 111
3164 // +---–---+
3165 // | 001 | 101
3166 // | +---–---+
3167 // | | | |
3168 // +---|---+ |
3169 // 010 | 110 |
3170 // +-------+
3171 // 000 100
3172 //
3173 G4int vcheck = 0;
3174 // check (ix - 1) side
3175 vcheck = (ix == 0) ? 0 : voxels[(ix-1)*kx
3176 if (vcheck == 0)
3177 {
3178 nfac++;
3179 G4int i1 = (ix+0)*kvx + (iy+0)*kvy + (i
3180 G4int i2 = (ix+0)*kvx + (iy+0)*kvy + (i
3181 G4int i3 = (ix+0)*kvx + (iy+1)*kvy + (i
3182 G4int i4 = (ix+0)*kvx + (iy+1)*kvy + (i
3183 if (indices[i1] == 0) indices[i1] = ++n
3184 if (indices[i2] == 0) indices[i2] = ++n
3185 if (indices[i3] == 0) indices[i3] = ++n
3186 if (indices[i4] == 0) indices[i4] = ++n
3187 }
3188 // check (ix + 1) side
3189 vcheck = (ix == nx - 1) ? 0 : voxels[(ix+
3190 if (vcheck == 0)
3191 {
3192 nfac++;
3193 G4int i1 = (ix+1)*kvx + (iy+1)*kvy + (i
3194 G4int i2 = (ix+1)*kvx + (iy+1)*kvy + (i
3195 G4int i3 = (ix+1)*kvx + (iy+0)*kvy + (i
3196 G4int i4 = (ix+1)*kvx + (iy+0)*kvy + (i
3197 if (indices[i1] == 0) indices[i1] = ++n
3198 if (indices[i2] == 0) indices[i2] = ++n
3199 if (indices[i3] == 0) indices[i3] = ++n
3200 if (indices[i4] == 0) indices[i4] = ++n
3201 }
3202 // check (iy - 1) side
3203 vcheck = (iy == 0) ? 0 : voxels[ix*kx + (
3204 if (vcheck == 0)
3205 {
3206 nfac++;
3207 G4int i1 = (ix+0)*kvx + (iy+0)*kvy + (i
3208 G4int i2 = (ix+1)*kvx + (iy+0)*kvy + (i
3209 G4int i3 = (ix+1)*kvx + (iy+0)*kvy + (i
3210 G4int i4 = (ix+0)*kvx + (iy+0)*kvy + (i
3211 if (indices[i1] == 0) indices[i1] = ++n
3212 if (indices[i2] == 0) indices[i2] = ++n
3213 if (indices[i3] == 0) indices[i3] = ++n
3214 if (indices[i4] == 0) indices[i4] = ++n
3215 }
3216 // check (iy + 1) side
3217 vcheck = (iy == ny - 1) ? 0 : voxels[ix*k
3218 if (vcheck == 0)
3219 {
3220 nfac++;
3221 G4int i1 = (ix+0)*kvx + (iy+1)*kvy + (i
3222 G4int i2 = (ix+0)*kvx + (iy+1)*kvy + (i
3223 G4int i3 = (ix+1)*kvx + (iy+1)*kvy + (i
3224 G4int i4 = (ix+1)*kvx + (iy+1)*kvy + (i
3225 if (indices[i1] == 0) indices[i1] = ++n
3226 if (indices[i2] == 0) indices[i2] = ++n
3227 if (indices[i3] == 0) indices[i3] = ++n
3228 if (indices[i4] == 0) indices[i4] = ++n
3229 }
3230 // check (iz - 1) side
3231 vcheck = (iz == 0) ? 0 : voxels[ix*kx + i
3232 if (vcheck == 0)
3233 {
3234 nfac++;
3235 G4int i1 = (ix+0)*kvx + (iy+0)*kvy + (i
3236 G4int i2 = (ix+0)*kvx + (iy+1)*kvy + (i
3237 G4int i3 = (ix+1)*kvx + (iy+1)*kvy + (i
3238 G4int i4 = (ix+1)*kvx + (iy+0)*kvy + (i
3239 if (indices[i1] == 0) indices[i1] = ++n
3240 if (indices[i2] == 0) indices[i2] = ++n
3241 if (indices[i3] == 0) indices[i3] = ++n
3242 if (indices[i4] == 0) indices[i4] = ++n
3243 }
3244 // check (iz + 1) side
3245 vcheck = (iz == nz - 1) ? 0 : voxels[ix*k
3246 if (vcheck == 0)
3247 {
3248 nfac++;
3249 G4int i1 = (ix+0)*kvx + (iy+0)*kvy + (i
3250 G4int i2 = (ix+1)*kvx + (iy+0)*kvy + (i
3251 G4int i3 = (ix+1)*kvx + (iy+1)*kvy + (i
3252 G4int i4 = (ix+0)*kvx + (iy+1)*kvy + (i
3253 if (indices[i1] == 0) indices[i1] = ++n
3254 if (indices[i2] == 0) indices[i2] = ++n
3255 if (indices[i3] == 0) indices[i3] = ++n
3256 if (indices[i4] == 0) indices[i4] = ++n
3257 }
3258 }
3259 // Construct polyhedron
3260 AllocateMemory(nver, nfac);
3261 G4ThreeVector p0(pmin.x() - 0.5*sizeX, pmin
3262 for (G4int ix = 0; ix <= nx; ++ix)
3263 {
3264 for (G4int iy = 0; iy <= ny; ++iy)
3265 {
3266 for (G4int iz = 0; iz <= nz; ++iz)
3267 {
3268 G4int i = ix*kvx + iy*kvy + iz;
3269 if (indices[i] == 0) continue;
3270 SetVertex(indices[i], p0 + G4ThreeVector(ix
3271 }
3272 }
3273 }
3274 nfac = 0;
3275 for (const auto& p: positions)
3276 {
3277 G4int ix = (p.x() - pmin.x())*invx + 0.5;
3278 G4int iy = (p.y() - pmin.y())*invy + 0.5;
3279 G4int iz = (p.z() - pmin.z())*invz + 0.5;
3280 G4int vcheck = 0;
3281 // check (ix - 1) side
3282 vcheck = (ix == 0) ? 0 : voxels[(ix-1)*kx
3283 if (vcheck == 0)
3284 {
3285 G4int i1 = (ix+0)*kvx + (iy+0)*kvy + (i
3286 G4int i2 = (ix+0)*kvx + (iy+0)*kvy + (i
3287 G4int i3 = (ix+0)*kvx + (iy+1)*kvy + (i
3288 G4int i4 = (ix+0)*kvx + (iy+1)*kvy + (i
3289 SetFacet(++nfac, indices[i1], indices[i
3290 }
3291 // check (ix + 1) side
3292 vcheck = (ix == nx - 1) ? 0 : voxels[(ix+
3293 if (vcheck == 0)
3294 {
3295 G4int i1 = (ix+1)*kvx + (iy+1)*kvy + (i
3296 G4int i2 = (ix+1)*kvx + (iy+1)*kvy + (i
3297 G4int i3 = (ix+1)*kvx + (iy+0)*kvy + (i
3298 G4int i4 = (ix+1)*kvx + (iy+0)*kvy + (i
3299 SetFacet(++nfac, indices[i1], indices[i
3300
3301 }
3302 // check (iy - 1) side
3303 vcheck = (iy == 0) ? 0 : voxels[ix*kx + (
3304 if (vcheck == 0)
3305 {
3306 G4int i1 = (ix+0)*kvx + (iy+0)*kvy + (i
3307 G4int i2 = (ix+1)*kvx + (iy+0)*kvy + (i
3308 G4int i3 = (ix+1)*kvx + (iy+0)*kvy + (i
3309 G4int i4 = (ix+0)*kvx + (iy+0)*kvy + (i
3310 SetFacet(++nfac, indices[i1], indices[i
3311 }
3312 // check (iy + 1) side
3313 vcheck = (iy == ny - 1) ? 0 : voxels[ix*k
3314 if (vcheck == 0)
3315 {
3316 G4int i1 = (ix+0)*kvx + (iy+1)*kvy + (i
3317 G4int i2 = (ix+0)*kvx + (iy+1)*kvy + (i
3318 G4int i3 = (ix+1)*kvx + (iy+1)*kvy + (i
3319 G4int i4 = (ix+1)*kvx + (iy+1)*kvy + (i
3320 SetFacet(++nfac, indices[i1], indices[i
3321 }
3322 // check (iz - 1) side
3323 vcheck = (iz == 0) ? 0 : voxels[ix*kx + i
3324 if (vcheck == 0)
3325 {
3326 G4int i1 = (ix+0)*kvx + (iy+0)*kvy + (i
3327 G4int i2 = (ix+0)*kvx + (iy+1)*kvy + (i
3328 G4int i3 = (ix+1)*kvx + (iy+1)*kvy + (i
3329 G4int i4 = (ix+1)*kvx + (iy+0)*kvy + (i
3330 SetFacet(++nfac, indices[i1], indices[i
3331 }
3332 // check (iz + 1) side
3333 vcheck = (iz == nz - 1) ? 0 : voxels[ix*k
3334 if (vcheck == 0)
3335 {
3336 G4int i1 = (ix+0)*kvx + (iy+0)*kvy + (i
3337 G4int i2 = (ix+1)*kvx + (iy+0)*kvy + (i
3338 G4int i3 = (ix+1)*kvx + (iy+1)*kvy + (i
3339 G4int i4 = (ix+0)*kvx + (iy+1)*kvy + (i
3340 SetFacet(++nfac, indices[i1], indices[i
3341 }
3342 }
3343 SetReferences();
3344 }
3345
3346 HepPolyhedronBoxMesh::~HepPolyhedronBoxMesh()
3347
3348 G4ThreadLocal
3349 G4int HepPolyhedron::fNumberOfRotationSteps =
3350 /********************************************
3351 *
3352 * Name: HepPolyhedron::fNumberOfRotationStep
3353 * Author: J.Allison (Manchester University)
3354 *
3355 * Function: Number of steps for whole circle
3356 *
3357 ********************************************
3358
3359 #include "BooleanProcessor.src"
3360
3361 HepPolyhedron HepPolyhedron::add(const HepPol
3362 /********************************************
3363 *
3364 * Name: HepPolyhedron::add
3365 * Author: E.Chernyaev
3366 *
3367 * Function: Boolean "union" of two polyhedra
3368 *
3369 ********************************************
3370 {
3371 G4int ierr;
3372 BooleanProcessor processor;
3373 return processor.execute(OP_UNION, *this, p
3374 }
3375
3376 HepPolyhedron HepPolyhedron::intersect(const
3377 /********************************************
3378 *
3379 * Name: HepPolyhedron::intersect
3380 * Author: E.Chernyaev
3381 *
3382 * Function: Boolean "intersection" of two po
3383 *
3384 ********************************************
3385 {
3386 G4int ierr;
3387 BooleanProcessor processor;
3388 return processor.execute(OP_INTERSECTION, *
3389 }
3390
3391 HepPolyhedron HepPolyhedron::subtract(const H
3392 /********************************************
3393 *
3394 * Name: HepPolyhedron::add
3395 * Author: E.Chernyaev
3396 *
3397 * Function: Boolean "subtraction" of "p" fro
3398 *
3399 ********************************************
3400 {
3401 G4int ierr;
3402 BooleanProcessor processor;
3403 return processor.execute(OP_SUBTRACTION, *t
3404 }
3405
3406 //NOTE : include the code of HepPolyhedronPro
3407 // since there is no BooleanProcessor.h
3408
3409 #undef INTERSECTION
3410
3411 #include "HepPolyhedronProcessor.src"
3412