Geant4-11
HepPolyhedron.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// G4 Polyhedron library
27//
28// History:
29// 23.07.96 E.Chernyaev <Evgueni.Tcherniaev@cern.ch> - initial version
30//
31// 30.09.96 E.Chernyaev
32// - added GetNextVertexIndex, GetVertex by Yasuhide Sawada
33// - added GetNextUnitNormal, GetNextEdgeIndices, GetNextEdge
34//
35// 15.12.96 E.Chernyaev
36// - added GetNumberOfRotationSteps, RotateEdge, RotateAroundZ, SetReferences
37// - rewritten G4PolyhedronCons;
38// - added G4PolyhedronPara, ...Trap, ...Pgon, ...Pcon, ...Sphere, ...Torus
39//
40// 01.06.97 E.Chernyaev
41// - modified RotateAroundZ, added SetSideFacets
42//
43// 19.03.00 E.Chernyaev
44// - implemented boolean operations (add, subtract, intersect) on polyhedra;
45//
46// 25.05.01 E.Chernyaev
47// - added GetSurfaceArea() and GetVolume();
48//
49// 05.11.02 E.Chernyaev
50// - added createTwistedTrap() and createPolyhedron();
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, HepPolyhedronHyberbolicMirror
60//
61// 12.05.21 E.Chernyaev
62// - added TriangulatePolygon(), RotateContourAroundZ()
63// - added HepPolyhedronPgon, HepPolyhedronPcon given by rz-contour
64//
65
66#include "HepPolyhedron.h"
67#include "G4PhysicalConstants.hh"
68#include "G4Vector3D.hh"
69
70#include <cstdlib> // Required on some compilers for std::abs(int) ...
71#include <cmath>
72#include <algorithm>
73
74using CLHEP::perMillion;
75using CLHEP::deg;
76using CLHEP::pi;
77using CLHEP::twopi;
78using CLHEP::nm;
79const G4double spatialTolerance = 0.01*nm;
80
81/***********************************************************************
82 * *
83 * Name: HepPolyhedron operator << Date: 09.05.96 *
84 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
85 * *
86 * Function: Print contents of G4 polyhedron *
87 * *
88 ***********************************************************************/
89std::ostream & operator<<(std::ostream & ostr, const G4Facet & facet) {
90 for (G4int k=0; k<4; k++) {
91 ostr << " " << facet.edge[k].v << "/" << facet.edge[k].f;
92 }
93 return ostr;
94}
95
96std::ostream & operator<<(std::ostream & ostr, const HepPolyhedron & ph) {
97 ostr << std::endl;
98 ostr << "Nvertices=" << ph.nvert << ", Nfacets=" << ph.nface << std::endl;
99 G4int i;
100 for (i=1; i<=ph.nvert; i++) {
101 ostr << "xyz(" << i << ")="
102 << ph.pV[i].x() << ' ' << ph.pV[i].y() << ' ' << ph.pV[i].z()
103 << std::endl;
104 }
105 for (i=1; i<=ph.nface; i++) {
106 ostr << "face(" << i << ")=" << ph.pF[i] << std::endl;
107 }
108 return ostr;
109}
110
111HepPolyhedron::HepPolyhedron(const HepPolyhedron &from)
112/***********************************************************************
113 * *
114 * Name: HepPolyhedron copy constructor Date: 23.07.96 *
115 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
116 * *
117 ***********************************************************************/
118: nvert(0), nface(0), pV(0), pF(0)
119{
120 AllocateMemory(from.nvert, from.nface);
121 for (G4int i=1; i<=nvert; i++) pV[i] = from.pV[i];
122 for (G4int k=1; k<=nface; k++) pF[k] = from.pF[k];
123}
124
125HepPolyhedron::HepPolyhedron(HepPolyhedron&& from)
126/***********************************************************************
127 * *
128 * Name: HepPolyhedron move constructor Date: 04.11.2019 *
129 * Author: E.Tcherniaev (E.Chernyaev) Revised: *
130 * *
131 ***********************************************************************/
132: nvert(0), nface(0), pV(nullptr), pF(nullptr)
133{
134 nvert = from.nvert;
135 nface = from.nface;
136 pV = from.pV;
137 pF = from.pF;
138
139 // Release the data from the source object
140 from.nvert = 0;
141 from.nface = 0;
142 from.pV = nullptr;
143 from.pF = nullptr;
144}
145
146HepPolyhedron & HepPolyhedron::operator=(const HepPolyhedron &from)
147/***********************************************************************
148 * *
149 * Name: HepPolyhedron operator = Date: 23.07.96 *
150 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
151 * *
152 * Function: Copy contents of one polyhedron to another *
153 * *
154 ***********************************************************************/
155{
156 if (this != &from) {
157 AllocateMemory(from.nvert, from.nface);
158 for (G4int i=1; i<=nvert; i++) pV[i] = from.pV[i];
159 for (G4int k=1; k<=nface; k++) pF[k] = from.pF[k];
160 }
161 return *this;
162}
163
164HepPolyhedron & HepPolyhedron::operator=(HepPolyhedron&& from)
165/***********************************************************************
166 * *
167 * Name: HepPolyhedron move operator = Date: 04.11.2019 *
168 * Author: E.Tcherniaev (E.Chernyaev) Revised: *
169 * *
170 * Function: Move contents of one polyhedron to another *
171 * *
172 ***********************************************************************/
173{
174 if (this != &from) {
175 delete [] pV;
176 delete [] pF;
177 nvert = from.nvert;
178 nface = from.nface;
179 pV = from.pV;
180 pF = from.pF;
181
182 // Release the data from the source object
183 from.nvert = 0;
184 from.nface = 0;
185 from.pV = nullptr;
186 from.pF = nullptr;
187 }
188 return *this;
189}
190
191G4int
192HepPolyhedron::FindNeighbour(G4int iFace, G4int iNode, G4int iOrder) const
193/***********************************************************************
194 * *
195 * Name: HepPolyhedron::FindNeighbour Date: 22.11.99 *
196 * Author: E.Chernyaev Revised: *
197 * *
198 * Function: Find neighbouring face *
199 * *
200 ***********************************************************************/
201{
202 G4int i;
203 for (i=0; i<4; i++) {
204 if (iNode == std::abs(pF[iFace].edge[i].v)) break;
205 }
206 if (i == 4) {
207 std::cerr
208 << "HepPolyhedron::FindNeighbour: face " << iFace
209 << " has no node " << iNode
210 << std::endl;
211 return 0;
212 }
213 if (iOrder < 0) {
214 if ( --i < 0) i = 3;
215 if (pF[iFace].edge[i].v == 0) i = 2;
216 }
217 return (pF[iFace].edge[i].v > 0) ? 0 : pF[iFace].edge[i].f;
218}
219
220G4Normal3D HepPolyhedron::FindNodeNormal(G4int iFace, G4int iNode) const
221/***********************************************************************
222 * *
223 * Name: HepPolyhedron::FindNodeNormal Date: 22.11.99 *
224 * Author: E.Chernyaev Revised: *
225 * *
226 * Function: Find normal at given node *
227 * *
228 ***********************************************************************/
229{
230 G4Normal3D normal = GetUnitNormal(iFace);
231 G4int k = iFace, iOrder = 1, n = 1;
232
233 for(;;) {
234 k = FindNeighbour(k, iNode, iOrder);
235 if (k == iFace) break;
236 if (k > 0) {
237 n++;
238 normal += GetUnitNormal(k);
239 }else{
240 if (iOrder < 0) break;
241 k = iFace;
242 iOrder = -iOrder;
243 }
244 }
245 return normal.unit();
246}
247
248G4int HepPolyhedron::GetNumberOfRotationSteps()
249/***********************************************************************
250 * *
251 * Name: HepPolyhedron::GetNumberOfRotationSteps Date: 24.06.97 *
252 * Author: J.Allison (Manchester University) Revised: *
253 * *
254 * Function: Get number of steps for whole circle *
255 * *
256 ***********************************************************************/
257{
258 return fNumberOfRotationSteps;
259}
260
261void HepPolyhedron::SetNumberOfRotationSteps(G4int n)
262/***********************************************************************
263 * *
264 * Name: HepPolyhedron::SetNumberOfRotationSteps Date: 24.06.97 *
265 * Author: J.Allison (Manchester University) Revised: *
266 * *
267 * Function: Set number of steps for whole circle *
268 * *
269 ***********************************************************************/
270{
271 const G4int nMin = 3;
272 if (n < nMin) {
273 std::cerr
274 << "HepPolyhedron::SetNumberOfRotationSteps: attempt to set the\n"
275 << "number of steps per circle < " << nMin << "; forced to " << nMin
276 << std::endl;
277 fNumberOfRotationSteps = nMin;
278 }else{
279 fNumberOfRotationSteps = n;
280 }
281}
282
283void HepPolyhedron::ResetNumberOfRotationSteps()
284/***********************************************************************
285 * *
286 * Name: HepPolyhedron::GetNumberOfRotationSteps Date: 24.06.97 *
287 * Author: J.Allison (Manchester University) Revised: *
288 * *
289 * Function: Reset number of steps for whole circle to default value *
290 * *
291 ***********************************************************************/
292{
293 fNumberOfRotationSteps = DEFAULT_NUMBER_OF_STEPS;
294}
295
296void HepPolyhedron::AllocateMemory(G4int Nvert, G4int Nface)
297/***********************************************************************
298 * *
299 * Name: HepPolyhedron::AllocateMemory Date: 19.06.96 *
300 * Author: E.Chernyaev (IHEP/Protvino) Revised: 05.11.02 *
301 * *
302 * Function: Allocate memory for GEANT4 polyhedron *
303 * *
304 * Input: Nvert - number of nodes *
305 * Nface - number of faces *
306 * *
307 ***********************************************************************/
308{
309 if (nvert == Nvert && nface == Nface) return;
310 if (pV != 0) delete [] pV;
311 if (pF != 0) delete [] pF;
312 if (Nvert > 0 && Nface > 0) {
313 nvert = Nvert;
314 nface = Nface;
315 pV = new G4Point3D[nvert+1];
316 pF = new G4Facet[nface+1];
317 }else{
318 nvert = 0; nface = 0; pV = 0; pF = 0;
319 }
320}
321
322void HepPolyhedron::CreatePrism()
323/***********************************************************************
324 * *
325 * Name: HepPolyhedron::CreatePrism Date: 15.07.96 *
326 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
327 * *
328 * Function: Set facets for a prism *
329 * *
330 ***********************************************************************/
331{
332 enum {DUMMY, BOTTOM, LEFT, BACK, RIGHT, FRONT, TOP};
333
334 pF[1] = G4Facet(1,LEFT, 4,BACK, 3,RIGHT, 2,FRONT);
335 pF[2] = G4Facet(5,TOP, 8,BACK, 4,BOTTOM, 1,FRONT);
336 pF[3] = G4Facet(8,TOP, 7,RIGHT, 3,BOTTOM, 4,LEFT);
337 pF[4] = G4Facet(7,TOP, 6,FRONT, 2,BOTTOM, 3,BACK);
338 pF[5] = G4Facet(6,TOP, 5,LEFT, 1,BOTTOM, 2,RIGHT);
339 pF[6] = G4Facet(5,FRONT, 6,RIGHT, 7,BACK, 8,LEFT);
340}
341
342void HepPolyhedron::RotateEdge(G4int k1, G4int k2, G4double r1, G4double r2,
343 G4int v1, G4int v2, G4int vEdge,
344 G4bool ifWholeCircle, G4int nds, G4int &kface)
345/***********************************************************************
346 * *
347 * Name: HepPolyhedron::RotateEdge Date: 05.12.96 *
348 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
349 * *
350 * Function: Create set of facets by rotation of an edge around Z-axis *
351 * *
352 * Input: k1, k2 - end vertices of the edge *
353 * r1, r2 - radiuses of the end vertices *
354 * v1, v2 - visibility of edges produced by rotation of the end *
355 * vertices *
356 * vEdge - visibility of the edge *
357 * ifWholeCircle - is true in case of whole circle rotation *
358 * nds - number of discrete steps *
359 * r[] - r-coordinates *
360 * kface - current free cell in the pF array *
361 * *
362 ***********************************************************************/
363{
364 if (r1 == 0. && r2 == 0.) return;
365
366 G4int i;
367 G4int i1 = k1;
368 G4int i2 = k2;
369 G4int ii1 = ifWholeCircle ? i1 : i1+nds;
370 G4int ii2 = ifWholeCircle ? i2 : i2+nds;
371 G4int vv = ifWholeCircle ? vEdge : 1;
372
373 if (nds == 1) {
374 if (r1 == 0.) {
375 pF[kface++] = G4Facet(i1,0, v2*i2,0, (i2+1),0);
376 }else if (r2 == 0.) {
377 pF[kface++] = G4Facet(i1,0, i2,0, v1*(i1+1),0);
378 }else{
379 pF[kface++] = G4Facet(i1,0, v2*i2,0, (i2+1),0, v1*(i1+1),0);
380 }
381 }else{
382 if (r1 == 0.) {
383 pF[kface++] = G4Facet(vv*i1,0, v2*i2,0, vEdge*(i2+1),0);
384 for (i2++,i=1; i<nds-1; i2++,i++) {
385 pF[kface++] = G4Facet(vEdge*i1,0, v2*i2,0, vEdge*(i2+1),0);
386 }
387 pF[kface++] = G4Facet(vEdge*i1,0, v2*i2,0, vv*ii2,0);
388 }else if (r2 == 0.) {
389 pF[kface++] = G4Facet(vv*i1,0, vEdge*i2,0, v1*(i1+1),0);
390 for (i1++,i=1; i<nds-1; i1++,i++) {
391 pF[kface++] = G4Facet(vEdge*i1,0, vEdge*i2,0, v1*(i1+1),0);
392 }
393 pF[kface++] = G4Facet(vEdge*i1,0, vv*i2,0, v1*ii1,0);
394 }else{
395 pF[kface++] = G4Facet(vv*i1,0, v2*i2,0, vEdge*(i2+1),0,v1*(i1+1),0);
396 for (i1++,i2++,i=1; i<nds-1; i1++,i2++,i++) {
397 pF[kface++] = G4Facet(vEdge*i1,0, v2*i2,0, vEdge*(i2+1),0,v1*(i1+1),0);
398 }
399 pF[kface++] = G4Facet(vEdge*i1,0, v2*i2,0, vv*ii2,0, v1*ii1,0);
400 }
401 }
402}
403
404void HepPolyhedron::SetSideFacets(G4int ii[4], G4int vv[4],
405 G4int *kk, G4double *r,
406 G4double dphi, G4int nds, G4int &kface)
407/***********************************************************************
408 * *
409 * Name: HepPolyhedron::SetSideFacets Date: 20.05.97 *
410 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
411 * *
412 * Function: Set side facets for the case of incomplete rotation *
413 * *
414 * Input: ii[4] - indices of original vertices *
415 * vv[4] - visibility of edges *
416 * kk[] - indices of nodes *
417 * r[] - radiuses *
418 * dphi - delta phi *
419 * nds - number of discrete steps *
420 * kface - current free cell in the pF array *
421 * *
422 ***********************************************************************/
423{
424 G4int k1, k2, k3, k4;
425
426 if (std::abs(dphi-pi) < perMillion) { // half a circle
427 for (G4int i=0; i<4; i++) {
428 k1 = ii[i];
429 k2 = ii[(i+1)%4];
430 if (r[k1] == 0. && r[k2] == 0.) vv[i] = -1;
431 }
432 }
433
434 if (ii[1] == ii[2]) {
435 k1 = kk[ii[0]];
436 k2 = kk[ii[2]];
437 k3 = kk[ii[3]];
438 pF[kface++] = G4Facet(vv[0]*k1,0, vv[2]*k2,0, vv[3]*k3,0);
439 if (r[ii[0]] != 0.) k1 += nds;
440 if (r[ii[2]] != 0.) k2 += nds;
441 if (r[ii[3]] != 0.) k3 += nds;
442 pF[kface++] = G4Facet(vv[2]*k3,0, vv[0]*k2,0, vv[3]*k1,0);
443 }else if (kk[ii[0]] == kk[ii[1]]) {
444 k1 = kk[ii[0]];
445 k2 = kk[ii[2]];
446 k3 = kk[ii[3]];
447 pF[kface++] = G4Facet(vv[1]*k1,0, vv[2]*k2,0, vv[3]*k3,0);
448 if (r[ii[0]] != 0.) k1 += nds;
449 if (r[ii[2]] != 0.) k2 += nds;
450 if (r[ii[3]] != 0.) k3 += nds;
451 pF[kface++] = G4Facet(vv[2]*k3,0, vv[1]*k2,0, vv[3]*k1,0);
452 }else if (kk[ii[2]] == kk[ii[3]]) {
453 k1 = kk[ii[0]];
454 k2 = kk[ii[1]];
455 k3 = kk[ii[2]];
456 pF[kface++] = G4Facet(vv[0]*k1,0, vv[1]*k2,0, vv[3]*k3,0);
457 if (r[ii[0]] != 0.) k1 += nds;
458 if (r[ii[1]] != 0.) k2 += nds;
459 if (r[ii[2]] != 0.) k3 += nds;
460 pF[kface++] = G4Facet(vv[1]*k3,0, vv[0]*k2,0, vv[3]*k1,0);
461 }else{
462 k1 = kk[ii[0]];
463 k2 = kk[ii[1]];
464 k3 = kk[ii[2]];
465 k4 = kk[ii[3]];
466 pF[kface++] = G4Facet(vv[0]*k1,0, vv[1]*k2,0, vv[2]*k3,0, vv[3]*k4,0);
467 if (r[ii[0]] != 0.) k1 += nds;
468 if (r[ii[1]] != 0.) k2 += nds;
469 if (r[ii[2]] != 0.) k3 += nds;
470 if (r[ii[3]] != 0.) k4 += nds;
471 pF[kface++] = G4Facet(vv[2]*k4,0, vv[1]*k3,0, vv[0]*k2,0, vv[3]*k1,0);
472 }
473}
474
475void HepPolyhedron::RotateAroundZ(G4int nstep, G4double phi, G4double dphi,
476 G4int np1, G4int np2,
477 const G4double *z, G4double *r,
478 G4int nodeVis, G4int edgeVis)
479/***********************************************************************
480 * *
481 * Name: HepPolyhedron::RotateAroundZ Date: 27.11.96 *
482 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
483 * *
484 * Function: Create HepPolyhedron for a solid produced by rotation of *
485 * two polylines around Z-axis *
486 * *
487 * Input: nstep - number of discrete steps, if 0 then default *
488 * phi - starting phi angle *
489 * dphi - delta phi *
490 * np1 - number of points in external polyline *
491 * (must be negative in case of closed polyline) *
492 * np2 - number of points in internal polyline (may be 1) *
493 * z[] - z-coordinates (+z >>> -z for both polylines) *
494 * r[] - r-coordinates *
495 * nodeVis - how to Draw edges joing consecutive positions of *
496 * node during rotation *
497 * edgeVis - how to Draw edges *
498 * *
499 ***********************************************************************/
500{
501 static const G4double wholeCircle = twopi;
502
503 // S E T R O T A T I O N P A R A M E T E R S
504
505 G4bool ifWholeCircle = (std::abs(dphi-wholeCircle) < perMillion) ? true : false;
506 G4double delPhi = ifWholeCircle ? wholeCircle : dphi;
507 G4int nSphi = nstep;
508 if (nSphi <= 0) nSphi = GetNumberOfRotationSteps()*delPhi/wholeCircle + 0.5;
509 if (nSphi == 0) nSphi = 1;
510 G4int nVphi = ifWholeCircle ? nSphi : nSphi + 1;
511 G4bool ifClosed = np1 > 0 ? false : true; // true if external polyline is closed
512
513 // C O U N T V E R T I C E S
514
515 G4int absNp1 = std::abs(np1);
516 G4int absNp2 = std::abs(np2);
517 G4int i1beg = 0;
518 G4int i1end = absNp1-1;
519 G4int i2beg = absNp1;
520 G4int i2end = absNp1+absNp2-1;
521 G4int i, j, k;
522
523 for(i=i1beg; i<=i2end; i++) {
524 if (std::abs(r[i]) < spatialTolerance) r[i] = 0.;
525 }
526
527 // external polyline - check position of nodes relative to Z
528 //
529 G4int Nverts = 0;
530 for (i=i1beg; i<=i1end; i++) {
531 Nverts += (r[i] == 0.) ? 1 : nVphi;
532 }
533
534 // internal polyline
535 //
536 G4bool ifSide1 = false; // whether to create bottom faces
537 G4bool ifSide2 = false; // whether to create top faces
538
539 if (r[i2beg] != r[i1beg] || z[i2beg] != z[i1beg]) { // first node
540 Nverts += (r[i2beg] == 0.) ? 1 : nVphi;
541 ifSide1 = true;
542 }
543
544 for(i=i2beg+1; i<i2end; i++) { // intermediate nodes
545 Nverts += (r[i] == 0.) ? 1 : nVphi;
546 }
547
548 if (r[i2end] != r[i1end] || z[i2end] != z[i1end]) { // last node
549 if (absNp2 > 1) Nverts += (r[i2end] == 0.) ? 1 : nVphi;
550 ifSide2 = true;
551 }
552
553 // C O U N T F A C E S
554
555 // external lateral faces
556 //
557 G4int Nfaces = ifClosed ? absNp1*nSphi : (absNp1-1)*nSphi;
558
559 // internal lateral faces
560 //
561 if (absNp2 > 1) {
562 for(i=i2beg; i<i2end; i++) {
563 if (r[i] > 0. || r[i+1] > 0.) Nfaces += nSphi;
564 }
565
566 if (ifClosed) {
567 if (r[i2end] > 0. || r[i2beg] > 0.) Nfaces += nSphi;
568 }
569 }
570
571 // bottom and top faces
572 //
573 if (!ifClosed) {
574 if (ifSide1 && (r[i1beg] > 0. || r[i2beg] > 0.)) Nfaces += nSphi;
575 if (ifSide2 && (r[i1end] > 0. || r[i2end] > 0.)) Nfaces += nSphi;
576 }
577
578 // phi_wedge faces
579 //
580 if (!ifWholeCircle) {
581 Nfaces += ifClosed ? 2*absNp1 : 2*(absNp1-1);
582 }
583
584 // A L L O C A T E M E M O R Y
585
586 AllocateMemory(Nverts, Nfaces);
587 if (pV == nullptr || pF == nullptr) return;
588
589 // G E N E R A T E V E R T I C E S
590
591 G4int *kk; // array of start indices along polylines
592 kk = new G4int[absNp1+absNp2];
593
594 // external polyline
595 //
596 k = 1; // free position in array of vertices pV
597 for(i=i1beg; i<=i1end; i++) {
598 kk[i] = k;
599 if (r[i] == 0.)
600 { pV[k++] = G4Point3D(0, 0, z[i]); } else { k += nVphi; }
601 }
602
603 // first point of internal polyline
604 //
605 i = i2beg;
606 if (ifSide1) {
607 kk[i] = k;
608 if (r[i] == 0.)
609 { pV[k++] = G4Point3D(0, 0, z[i]); } else { k += nVphi; }
610 }else{
611 kk[i] = kk[i1beg];
612 }
613
614 // intermediate points of internal polyline
615 //
616 for(i=i2beg+1; i<i2end; i++) {
617 kk[i] = k;
618 if (r[i] == 0.)
619 { pV[k++] = G4Point3D(0, 0, z[i]); } else { k += nVphi; }
620 }
621
622 // last point of internal polyline
623 //
624 if (absNp2 > 1) {
625 i = i2end;
626 if (ifSide2) {
627 kk[i] = k;
628 if (r[i] == 0.) pV[k] = G4Point3D(0, 0, z[i]);
629 }else{
630 kk[i] = kk[i1end];
631 }
632 }
633
634 // set vertices
635 //
636 G4double cosPhi, sinPhi;
637
638 for(j=0; j<nVphi; j++) {
639 cosPhi = std::cos(phi+j*delPhi/nSphi);
640 sinPhi = std::sin(phi+j*delPhi/nSphi);
641 for(i=i1beg; i<=i2end; i++) {
642 if (r[i] != 0.)
643 pV[kk[i]+j] = G4Point3D(r[i]*cosPhi,r[i]*sinPhi,z[i]);
644 }
645 }
646
647 // G E N E R A T E F A C E S
648
649 // external faces
650 //
651 G4int v1,v2;
652
653 k = 1; // free position in array of faces pF
654 v2 = ifClosed ? nodeVis : 1;
655 for(i=i1beg; i<i1end; i++) {
656 v1 = v2;
657 if (!ifClosed && i == i1end-1) {
658 v2 = 1;
659 }else{
660 v2 = (r[i] == r[i+1] && r[i+1] == r[i+2]) ? -1 : nodeVis;
661 }
662 RotateEdge(kk[i], kk[i+1], r[i], r[i+1], v1, v2,
663 edgeVis, ifWholeCircle, nSphi, k);
664 }
665 if (ifClosed) {
666 RotateEdge(kk[i1end], kk[i1beg], r[i1end],r[i1beg], nodeVis, nodeVis,
667 edgeVis, ifWholeCircle, nSphi, k);
668 }
669
670 // internal faces
671 //
672 if (absNp2 > 1) {
673 v2 = ifClosed ? nodeVis : 1;
674 for(i=i2beg; i<i2end; i++) {
675 v1 = v2;
676 if (!ifClosed && i==i2end-1) {
677 v2 = 1;
678 }else{
679 v2 = (r[i] == r[i+1] && r[i+1] == r[i+2]) ? -1 : nodeVis;
680 }
681 RotateEdge(kk[i+1], kk[i], r[i+1], r[i], v2, v1,
682 edgeVis, ifWholeCircle, nSphi, k);
683 }
684 if (ifClosed) {
685 RotateEdge(kk[i2beg], kk[i2end], r[i2beg], r[i2end], nodeVis, nodeVis,
686 edgeVis, ifWholeCircle, nSphi, k);
687 }
688 }
689
690 // bottom and top faces
691 //
692 if (!ifClosed) {
693 if (ifSide1) {
694 RotateEdge(kk[i2beg], kk[i1beg], r[i2beg], r[i1beg], 1, 1,
695 -1, ifWholeCircle, nSphi, k);
696 }
697 if (ifSide2) {
698 RotateEdge(kk[i1end], kk[i2end], r[i1end], r[i2end], 1, 1,
699 -1, ifWholeCircle, nSphi, k);
700 }
701 }
702
703 // phi_wedge faces in case of incomplete circle
704 //
705 if (!ifWholeCircle) {
706
707 G4int ii[4], vv[4];
708
709 if (ifClosed) {
710 for (i=i1beg; i<=i1end; i++) {
711 ii[0] = i;
712 ii[3] = (i == i1end) ? i1beg : i+1;
713 ii[1] = (absNp2 == 1) ? i2beg : ii[0]+absNp1;
714 ii[2] = (absNp2 == 1) ? i2beg : ii[3]+absNp1;
715 vv[0] = -1;
716 vv[1] = 1;
717 vv[2] = -1;
718 vv[3] = 1;
719 SetSideFacets(ii, vv, kk, r, delPhi, nSphi, k);
720 }
721 }else{
722 for (i=i1beg; i<i1end; i++) {
723 ii[0] = i;
724 ii[3] = i+1;
725 ii[1] = (absNp2 == 1) ? i2beg : ii[0]+absNp1;
726 ii[2] = (absNp2 == 1) ? i2beg : ii[3]+absNp1;
727 vv[0] = (i == i1beg) ? 1 : -1;
728 vv[1] = 1;
729 vv[2] = (i == i1end-1) ? 1 : -1;
730 vv[3] = 1;
731 SetSideFacets(ii, vv, kk, r, delPhi, nSphi, k);
732 }
733 }
734 }
735
736 delete [] kk; // free memory
737
738 // final check
739 //
740 if (k-1 != nface) {
741 std::cerr
742 << "HepPolyhedron::RotateAroundZ: number of generated faces ("
743 << k-1 << ") is not equal to the number of allocated faces ("
744 << nface << ")"
745 << std::endl;
746 }
747}
748
749void
750HepPolyhedron::RotateContourAroundZ(G4int nstep,
751 G4double phi,
752 G4double dphi,
753 const std::vector<G4TwoVector> &rz,
754 G4int nodeVis,
755 G4int edgeVis)
756/***********************************************************************
757 * *
758 * Name: HepPolyhedron::RotateContourAroundZ Date: 12.05.21 *
759 * Author: E.Tcherniaev (E.Chernyaev) Revised: *
760 * *
761 * Function: Create HepPolyhedron for a solid produced by rotation of *
762 * a closed polyline (rz-contour) around Z-axis *
763 * *
764 * Input: nstep - number of discrete steps, if 0 then default *
765 * phi - starting phi angle *
766 * dphi - delta phi *
767 * rz - rz-contour *
768 * nodeVis - how to Draw edges joing consecutive positions of *
769 * node during rotation *
770 * edgeVis - how to Draw edges *
771 * *
772 ***********************************************************************/
773{
774 // S E T R O T A T I O N P A R A M E T E R S
775
776 G4bool ifWholeCircle = (std::abs(dphi - twopi) < perMillion) ? true : false;
777 G4double delPhi = (ifWholeCircle) ? twopi : dphi;
778 G4int nSphi = nstep;
779 if (nSphi <= 0) nSphi = GetNumberOfRotationSteps()*delPhi/twopi + 0.5;
780 if (nSphi == 0) nSphi = 1;
781 G4int nVphi = (ifWholeCircle) ? nSphi : nSphi + 1;
782
783 // C A L C U L A T E A R E A
784
785 G4int Nrz = rz.size();
786 G4double area = 0;
787 for (G4int i = 0; i < Nrz; ++i)
788 {
789 G4int k = (i == 0) ? Nrz - 1 : i - 1;
790 area += rz[k].x()*rz[i].y() - rz[i].x()*rz[k].y();
791 }
792
793 // P R E P A R E P O L Y L I N E
794
795 G4double *r = new G4double[Nrz];
796 G4double *z = new G4double[Nrz];
797 for (G4int i = 0; i < Nrz; ++i)
798 {
799 r[i] = rz[i].x();
800 z[i] = rz[i].y();
801 if (std::abs(r[i]) < spatialTolerance) r[i] = 0.;
802 }
803
804 // C O U N T V E R T I C E S A N D F A C E S
805
806 G4int Nverts = 0;
807 for(G4int i = 0; i < Nrz; ++i) Nverts += (r[i] == 0.) ? 1 : nVphi;
808
809 G4int Nedges = Nrz;
810 for (G4int i = 0; i < Nrz; ++i)
811 {
812 G4int k = (i == 0) ? Nrz - 1 : i - 1;
813 Nedges -= (r[k] == 0 && r[i] == 0);
814 }
815
816 G4int Nfaces = Nedges*nSphi; // lateral faces
817 if (!ifWholeCircle) Nfaces += 2*(Nrz - 2); // phi_wedge faces
818
819 // A L L O C A T E M E M O R Y
820
821 AllocateMemory(Nverts, Nfaces);
822 if (pV == nullptr || pF == nullptr)
823 {
824 delete [] r;
825 delete [] z;
826 return;
827 }
828
829 // S E T V E R T I C E S
830
831 G4int *kk = new G4int[Nrz]; // start indices along contour
832 G4int kfree = 1; // current free position in array of vertices pV
833
834 // set start indices, set vertices for nodes with r == 0
835 for(G4int i = 0; i < Nrz; ++i)
836 {
837 kk[i] = kfree;
838 if (r[i] == 0.) pV[kfree++] = G4Point3D(0, 0, z[i]);
839 if (r[i] != 0.) kfree += nVphi;
840 }
841
842 // set vertices by rotating r
843 for(G4int j = 0; j < nVphi; ++j)
844 {
845 G4double cosPhi = std::cos(phi + j*delPhi/nSphi);
846 G4double sinPhi = std::sin(phi + j*delPhi/nSphi);
847 for(G4int i = 0; i < Nrz; ++i)
848 {
849 if (r[i] != 0.)
850 pV[kk[i] + j] = G4Point3D(r[i]*cosPhi, r[i]*sinPhi, z[i]);
851 }
852 }
853
854 // S E T F A C E S
855
856 kfree = 1; // current free position in array of faces pF
857 for(G4int i = 0; i < Nrz; ++i)
858 {
859 G4int i1 = (i < Nrz - 1) ? i + 1 : 0; // inverse order if area > 0
860 G4int i2 = i;
861 if (area < 0.) std::swap(i1, i2);
862 RotateEdge(kk[i1], kk[i2], r[i1], r[i2], nodeVis, nodeVis,
863 edgeVis, ifWholeCircle, nSphi, kfree);
864 }
865
866 // S E T P H I _ W E D G E F A C E S
867
868 if (!ifWholeCircle)
869 {
870 std::vector<G4int> triangles;
871 TriangulatePolygon(rz, triangles);
872
873 G4int ii[4], vv[4];
874 G4int ntria = triangles.size()/3;
875 for (G4int i = 0; i < ntria; ++i)
876 {
877 G4int i1 = triangles[0 + i*3];
878 G4int i2 = triangles[1 + i*3];
879 G4int i3 = triangles[2 + i*3];
880 if (area < 0.) std::swap(i1, i3);
881 G4int v1 = (std::abs(i2-i1) == 1 || std::abs(i2-i1) == Nrz-1) ? 1 : -1;
882 G4int v2 = (std::abs(i3-i2) == 1 || std::abs(i3-i2) == Nrz-1) ? 1 : -1;
883 G4int v3 = (std::abs(i1-i3) == 1 || std::abs(i1-i3) == Nrz-1) ? 1 : -1;
884 ii[0] = i1; ii[1] = i2; ii[2] = i2; ii[3] = i3;
885 vv[0] = v1; vv[1] = -1; vv[2] = v2; vv[3] = v3;
886 SetSideFacets(ii, vv, kk, r, delPhi, nSphi, kfree);
887 }
888 }
889
890 // free memory
891 delete [] r;
892 delete [] z;
893 delete [] kk;
894
895 // final check
896 if (kfree - 1 != nface)
897 {
898 std::cerr
899 << "HepPolyhedron::RotateContourAroundZ: number of generated faces ("
900 << kfree-1 << ") is not equal to the number of allocated faces ("
901 << nface << ")"
902 << std::endl;
903 }
904}
905
906G4bool
907HepPolyhedron::TriangulatePolygon(const std::vector<G4TwoVector> &polygon,
908 std::vector<G4int> &result)
909/***********************************************************************
910 * *
911 * Name: HepPolyhedron::TriangulatePolygon Date: 12.05.21 *
912 * Author: E.Tcherniaev (E.Chernyaev) Revised: *
913 * *
914 * Function: Simple implementation of "ear clipping" algorithm for *
915 * triangulation of a simple contour/polygon, it places *
916 * the result in a std::vector as triplets of vertex indices *
917 * *
918 * If triangulation is sucsessfull then the function *
919 * returns true, otherwise false *
920 * *
921 * Remark: It's a copy of G4GeomTools::TriangulatePolygon() *
922 * *
923 ***********************************************************************/
924{
925 result.resize(0);
926 G4int n = polygon.size();
927 if (n < 3) return false;
928
929 // calculate area
930 //
931 G4double area = 0.;
932 for(G4int i = 0; i < n; ++i)
933 {
934 G4int k = (i == 0) ? n - 1 : i - 1;
935 area += polygon[k].x()*polygon[i].y() - polygon[i].x()*polygon[k].y();
936 }
937
938 // allocate and initialize list of Vertices
939 // we want a counter-clockwise polygon in V
940 //
941 G4int* V = new G4int[n];
942 if (area > 0.)
943 for (G4int i = 0; i < n; ++i) V[i] = i;
944 else
945 for (G4int i = 0; i < n; ++i) V[i] = (n - 1) - i;
946
947 // Triangulation: remove nv-2 Vertices, creating 1 triangle every time
948 //
949 G4int nv = n;
950 G4int count = 2*nv; // error detection counter
951 for(G4int b = nv - 1; nv > 2; )
952 {
953 // ERROR: if we loop, it is probably a non-simple polygon
954 if ((count--) <= 0)
955 {
956 delete [] V;
957 if (area < 0.) std::reverse(result.begin(),result.end());
958 return false;
959 }
960
961 // three consecutive vertices in current polygon, <a,b,c>
962 G4int a = (b < nv) ? b : 0; // previous
963 b = (a+1 < nv) ? a+1 : 0; // current
964 G4int c = (b+1 < nv) ? b+1 : 0; // next
965
966 if (CheckSnip(polygon, a,b,c, nv,V))
967 {
968 // output Triangle
969 result.push_back(V[a]);
970 result.push_back(V[b]);
971 result.push_back(V[c]);
972
973 // remove vertex b from remaining polygon
974 nv--;
975 for(G4int i = b; i < nv; ++i) V[i] = V[i+1];
976
977 count = 2*nv; // resest error detection counter
978 }
979 }
980 delete [] V;
981 if (area < 0.) std::reverse(result.begin(),result.end());
982 return true;
983}
984
985G4bool HepPolyhedron::CheckSnip(const std::vector<G4TwoVector> &contour,
986 G4int a, G4int b, G4int c,
987 G4int n, const G4int* V)
988/***********************************************************************
989 * *
990 * Name: HepPolyhedron::CheckSnip Date: 12.05.21 *
991 * Author: E.Tcherniaev (E.Chernyaev) Revised: *
992 * *
993 * Function: Check for a valid snip, *
994 * it is a helper functionfor TriangulatePolygon() *
995 * *
996 ***********************************************************************/
997{
998 static const G4double kCarTolerance = 1.e-9;
999
1000 // check orientation of Triangle
1001 G4double Ax = contour[V[a]].x(), Ay = contour[V[a]].y();
1002 G4double Bx = contour[V[b]].x(), By = contour[V[b]].y();
1003 G4double Cx = contour[V[c]].x(), Cy = contour[V[c]].y();
1004 if ((Bx-Ax)*(Cy-Ay) - (By-Ay)*(Cx-Ax) < kCarTolerance) return false;
1005
1006 // check that there is no point inside Triangle
1007 G4double xmin = std::min(std::min(Ax,Bx),Cx);
1008 G4double xmax = std::max(std::max(Ax,Bx),Cx);
1009 G4double ymin = std::min(std::min(Ay,By),Cy);
1010 G4double ymax = std::max(std::max(Ay,By),Cy);
1011
1012 for (G4int i=0; i<n; ++i)
1013 {
1014 if((i == a) || (i == b) || (i == c)) continue;
1015 G4double Px = contour[V[i]].x();
1016 if (Px < xmin || Px > xmax) continue;
1017 G4double Py = contour[V[i]].y();
1018 if (Py < ymin || Py > ymax) continue;
1019 // if (PointInTriangle(Ax,Ay,Bx,By,Cx,Cy,Px,Py)) return false;
1020 if ((Bx-Ax)*(Cy-Ay) - (By-Ay)*(Cx-Ax) > 0.)
1021 {
1022 if ((Ax-Cx)*(Py-Cy) - (Ay-Cy)*(Px-Cx) < 0.) continue;
1023 if ((Bx-Ax)*(Py-Ay) - (By-Ay)*(Px-Ax) < 0.) continue;
1024 if ((Cx-Bx)*(Py-By) - (Cy-By)*(Px-Bx) < 0.) continue;
1025 }
1026 else
1027 {
1028 if ((Ax-Cx)*(Py-Cy) - (Ay-Cy)*(Px-Cx) > 0.) continue;
1029 if ((Bx-Ax)*(Py-Ay) - (By-Ay)*(Px-Ax) > 0.) continue;
1030 if ((Cx-Bx)*(Py-By) - (Cy-By)*(Px-Bx) > 0.) continue;
1031 }
1032 return false;
1033 }
1034 return true;
1035}
1036
1037void HepPolyhedron::SetReferences()
1038/***********************************************************************
1039 * *
1040 * Name: HepPolyhedron::SetReferences Date: 04.12.96 *
1041 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
1042 * *
1043 * Function: For each edge set reference to neighbouring facet *
1044 * *
1045 ***********************************************************************/
1046{
1047 if (nface <= 0) return;
1048
1049 struct edgeListMember {
1050 edgeListMember *next;
1051 G4int v2;
1052 G4int iface;
1053 G4int iedge;
1054 } *edgeList, *freeList, **headList;
1055
1056
1057 // A L L O C A T E A N D I N I T I A T E L I S T S
1058
1059 edgeList = new edgeListMember[2*nface];
1060 headList = new edgeListMember*[nvert];
1061
1062 G4int i;
1063 for (i=0; i<nvert; i++) {
1064 headList[i] = 0;
1065 }
1066 freeList = edgeList;
1067 for (i=0; i<2*nface-1; i++) {
1068 edgeList[i].next = &edgeList[i+1];
1069 }
1070 edgeList[2*nface-1].next = 0;
1071
1072 // L O O P A L O N G E D G E S
1073
1074 G4int iface, iedge, nedge, i1, i2, k1, k2;
1075 edgeListMember *prev, *cur;
1076
1077 for(iface=1; iface<=nface; iface++) {
1078 nedge = (pF[iface].edge[3].v == 0) ? 3 : 4;
1079 for (iedge=0; iedge<nedge; iedge++) {
1080 i1 = iedge;
1081 i2 = (iedge < nedge-1) ? iedge+1 : 0;
1082 i1 = std::abs(pF[iface].edge[i1].v);
1083 i2 = std::abs(pF[iface].edge[i2].v);
1084 k1 = (i1 < i2) ? i1 : i2; // k1 = ::min(i1,i2);
1085 k2 = (i1 > i2) ? i1 : i2; // k2 = ::max(i1,i2);
1086
1087 // check head of the List corresponding to k1
1088 cur = headList[k1];
1089 if (cur == 0) {
1090 headList[k1] = freeList;
1091 if (!freeList) {
1092 std::cerr
1093 << "Polyhedron::SetReferences: bad link "
1094 << std::endl;
1095 break;
1096 }
1097 freeList = freeList->next;
1098 cur = headList[k1];
1099 cur->next = 0;
1100 cur->v2 = k2;
1101 cur->iface = iface;
1102 cur->iedge = iedge;
1103 continue;
1104 }
1105
1106 if (cur->v2 == k2) {
1107 headList[k1] = cur->next;
1108 cur->next = freeList;
1109 freeList = cur;
1110 pF[iface].edge[iedge].f = cur->iface;
1111 pF[cur->iface].edge[cur->iedge].f = iface;
1112 i1 = (pF[iface].edge[iedge].v < 0) ? -1 : 1;
1113 i2 = (pF[cur->iface].edge[cur->iedge].v < 0) ? -1 : 1;
1114 if (i1 != i2) {
1115 std::cerr
1116 << "Polyhedron::SetReferences: different edge visibility "
1117 << iface << "/" << iedge << "/"
1118 << pF[iface].edge[iedge].v << " and "
1119 << cur->iface << "/" << cur->iedge << "/"
1120 << pF[cur->iface].edge[cur->iedge].v
1121 << std::endl;
1122 }
1123 continue;
1124 }
1125
1126 // check List itself
1127 for (;;) {
1128 prev = cur;
1129 cur = prev->next;
1130 if (cur == 0) {
1131 prev->next = freeList;
1132 if (!freeList) {
1133 std::cerr
1134 << "Polyhedron::SetReferences: bad link "
1135 << std::endl;
1136 break;
1137 }
1138 freeList = freeList->next;
1139 cur = prev->next;
1140 cur->next = 0;
1141 cur->v2 = k2;
1142 cur->iface = iface;
1143 cur->iedge = iedge;
1144 break;
1145 }
1146
1147 if (cur->v2 == k2) {
1148 prev->next = cur->next;
1149 cur->next = freeList;
1150 freeList = cur;
1151 pF[iface].edge[iedge].f = cur->iface;
1152 pF[cur->iface].edge[cur->iedge].f = iface;
1153 i1 = (pF[iface].edge[iedge].v < 0) ? -1 : 1;
1154 i2 = (pF[cur->iface].edge[cur->iedge].v < 0) ? -1 : 1;
1155 if (i1 != i2) {
1156 std::cerr
1157 << "Polyhedron::SetReferences: different edge visibility "
1158 << iface << "/" << iedge << "/"
1159 << pF[iface].edge[iedge].v << " and "
1160 << cur->iface << "/" << cur->iedge << "/"
1161 << pF[cur->iface].edge[cur->iedge].v
1162 << std::endl;
1163 }
1164 break;
1165 }
1166 }
1167 }
1168 }
1169
1170 // C H E C K T H A T A L L L I S T S A R E E M P T Y
1171
1172 for (i=0; i<nvert; i++) {
1173 if (headList[i] != 0) {
1174 std::cerr
1175 << "Polyhedron::SetReferences: List " << i << " is not empty"
1176 << std::endl;
1177 }
1178 }
1179
1180 // F R E E M E M O R Y
1181
1182 delete [] edgeList;
1183 delete [] headList;
1184}
1185
1186void HepPolyhedron::InvertFacets()
1187/***********************************************************************
1188 * *
1189 * Name: HepPolyhedron::InvertFacets Date: 01.12.99 *
1190 * Author: E.Chernyaev Revised: *
1191 * *
1192 * Function: Invert the order of the nodes in the facets *
1193 * *
1194 ***********************************************************************/
1195{
1196 if (nface <= 0) return;
1197 G4int i, k, nnode, v[4],f[4];
1198 for (i=1; i<=nface; i++) {
1199 nnode = (pF[i].edge[3].v == 0) ? 3 : 4;
1200 for (k=0; k<nnode; k++) {
1201 v[k] = (k+1 == nnode) ? pF[i].edge[0].v : pF[i].edge[k+1].v;
1202 if (v[k] * pF[i].edge[k].v < 0) v[k] = -v[k];
1203 f[k] = pF[i].edge[k].f;
1204 }
1205 for (k=0; k<nnode; k++) {
1206 pF[i].edge[nnode-1-k].v = v[k];
1207 pF[i].edge[nnode-1-k].f = f[k];
1208 }
1209 }
1210}
1211
1212HepPolyhedron & HepPolyhedron::Transform(const G4Transform3D &t)
1213/***********************************************************************
1214 * *
1215 * Name: HepPolyhedron::Transform Date: 01.12.99 *
1216 * Author: E.Chernyaev Revised: *
1217 * *
1218 * Function: Make transformation of the polyhedron *
1219 * *
1220 ***********************************************************************/
1221{
1222 if (nvert > 0) {
1223 for (G4int i=1; i<=nvert; i++) { pV[i] = t * pV[i]; }
1224
1225 // C H E C K D E T E R M I N A N T A N D
1226 // I N V E R T F A C E T S I F I T I S N E G A T I V E
1227
1228 G4Vector3D d = t * G4Vector3D(0,0,0);
1229 G4Vector3D x = t * G4Vector3D(1,0,0) - d;
1230 G4Vector3D y = t * G4Vector3D(0,1,0) - d;
1231 G4Vector3D z = t * G4Vector3D(0,0,1) - d;
1232 if ((x.cross(y))*z < 0) InvertFacets();
1233 }
1234 return *this;
1235}
1236
1237G4bool HepPolyhedron::GetNextVertexIndex(G4int &index, G4int &edgeFlag) const
1238/***********************************************************************
1239 * *
1240 * Name: HepPolyhedron::GetNextVertexIndex Date: 03.09.96 *
1241 * Author: Yasuhide Sawada Revised: *
1242 * *
1243 * Function: *
1244 * *
1245 ***********************************************************************/
1246{
1247 static G4ThreadLocal G4int iFace = 1;
1248 static G4ThreadLocal G4int iQVertex = 0;
1249 G4int vIndex = pF[iFace].edge[iQVertex].v;
1250
1251 edgeFlag = (vIndex > 0) ? 1 : 0;
1252 index = std::abs(vIndex);
1253
1254 if (iQVertex >= 3 || pF[iFace].edge[iQVertex+1].v == 0) {
1255 iQVertex = 0;
1256 if (++iFace > nface) iFace = 1;
1257 return false; // Last Edge
1258 }else{
1259 ++iQVertex;
1260 return true; // not Last Edge
1261 }
1262}
1263
1264G4Point3D HepPolyhedron::GetVertex(G4int index) const
1265/***********************************************************************
1266 * *
1267 * Name: HepPolyhedron::GetVertex Date: 03.09.96 *
1268 * Author: Yasuhide Sawada Revised: 17.11.99 *
1269 * *
1270 * Function: Get vertex of the index. *
1271 * *
1272 ***********************************************************************/
1273{
1274 if (index <= 0 || index > nvert) {
1275 std::cerr
1276 << "HepPolyhedron::GetVertex: irrelevant index " << index
1277 << std::endl;
1278 return G4Point3D();
1279 }
1280 return pV[index];
1281}
1282
1283G4bool
1284HepPolyhedron::GetNextVertex(G4Point3D &vertex, G4int &edgeFlag) const
1285/***********************************************************************
1286 * *
1287 * Name: HepPolyhedron::GetNextVertex Date: 22.07.96 *
1288 * Author: John Allison Revised: *
1289 * *
1290 * Function: Get vertices of the quadrilaterals in order for each *
1291 * face in face order. Returns false when finished each *
1292 * face. *
1293 * *
1294 ***********************************************************************/
1295{
1296 G4int index;
1297 G4bool rep = GetNextVertexIndex(index, edgeFlag);
1298 vertex = pV[index];
1299 return rep;
1300}
1301
1302G4bool HepPolyhedron::GetNextVertex(G4Point3D &vertex, G4int &edgeFlag,
1303 G4Normal3D &normal) const
1304/***********************************************************************
1305 * *
1306 * Name: HepPolyhedron::GetNextVertex Date: 26.11.99 *
1307 * Author: E.Chernyaev Revised: *
1308 * *
1309 * Function: Get vertices with normals of the quadrilaterals in order *
1310 * for each face in face order. *
1311 * Returns false when finished each face. *
1312 * *
1313 ***********************************************************************/
1314{
1315 static G4ThreadLocal G4int iFace = 1;
1316 static G4ThreadLocal G4int iNode = 0;
1317
1318 if (nface == 0) return false; // empty polyhedron
1319
1320 G4int k = pF[iFace].edge[iNode].v;
1321 if (k > 0) { edgeFlag = 1; } else { edgeFlag = -1; k = -k; }
1322 vertex = pV[k];
1323 normal = FindNodeNormal(iFace,k);
1324 if (iNode >= 3 || pF[iFace].edge[iNode+1].v == 0) {
1325 iNode = 0;
1326 if (++iFace > nface) iFace = 1;
1327 return false; // last node
1328 }else{
1329 ++iNode;
1330 return true; // not last node
1331 }
1332}
1333
1334G4bool HepPolyhedron::GetNextEdgeIndices(G4int &i1, G4int &i2, G4int &edgeFlag,
1335 G4int &iface1, G4int &iface2) const
1336/***********************************************************************
1337 * *
1338 * Name: HepPolyhedron::GetNextEdgeIndices Date: 30.09.96 *
1339 * Author: E.Chernyaev Revised: 17.11.99 *
1340 * *
1341 * Function: Get indices of the next edge together with indices of *
1342 * of the faces which share the edge. *
1343 * Returns false when the last edge. *
1344 * *
1345 ***********************************************************************/
1346{
1347 static G4ThreadLocal G4int iFace = 1;
1348 static G4ThreadLocal G4int iQVertex = 0;
1349 static G4ThreadLocal G4int iOrder = 1;
1350 G4int k1, k2, kflag, kface1, kface2;
1351
1352 if (iFace == 1 && iQVertex == 0) {
1353 k2 = pF[nface].edge[0].v;
1354 k1 = pF[nface].edge[3].v;
1355 if (k1 == 0) k1 = pF[nface].edge[2].v;
1356 if (std::abs(k1) > std::abs(k2)) iOrder = -1;
1357 }
1358
1359 do {
1360 k1 = pF[iFace].edge[iQVertex].v;
1361 kflag = k1;
1362 k1 = std::abs(k1);
1363 kface1 = iFace;
1364 kface2 = pF[iFace].edge[iQVertex].f;
1365 if (iQVertex >= 3 || pF[iFace].edge[iQVertex+1].v == 0) {
1366 iQVertex = 0;
1367 k2 = std::abs(pF[iFace].edge[iQVertex].v);
1368 iFace++;
1369 }else{
1370 iQVertex++;
1371 k2 = std::abs(pF[iFace].edge[iQVertex].v);
1372 }
1373 } while (iOrder*k1 > iOrder*k2);
1374
1375 i1 = k1; i2 = k2; edgeFlag = (kflag > 0) ? 1 : 0;
1376 iface1 = kface1; iface2 = kface2;
1377
1378 if (iFace > nface) {
1379 iFace = 1; iOrder = 1;
1380 return false;
1381 }else{
1382 return true;
1383 }
1384}
1385
1386G4bool
1387HepPolyhedron::GetNextEdgeIndices(G4int &i1, G4int &i2, G4int &edgeFlag) const
1388/***********************************************************************
1389 * *
1390 * Name: HepPolyhedron::GetNextEdgeIndices Date: 17.11.99 *
1391 * Author: E.Chernyaev Revised: *
1392 * *
1393 * Function: Get indices of the next edge. *
1394 * Returns false when the last edge. *
1395 * *
1396 ***********************************************************************/
1397{
1398 G4int kface1, kface2;
1399 return GetNextEdgeIndices(i1, i2, edgeFlag, kface1, kface2);
1400}
1401
1402G4bool
1403HepPolyhedron::GetNextEdge(G4Point3D &p1,
1404 G4Point3D &p2,
1405 G4int &edgeFlag) const
1406/***********************************************************************
1407 * *
1408 * Name: HepPolyhedron::GetNextEdge Date: 30.09.96 *
1409 * Author: E.Chernyaev Revised: *
1410 * *
1411 * Function: Get next edge. *
1412 * Returns false when the last edge. *
1413 * *
1414 ***********************************************************************/
1415{
1416 G4int i1,i2;
1417 G4bool rep = GetNextEdgeIndices(i1,i2,edgeFlag);
1418 p1 = pV[i1];
1419 p2 = pV[i2];
1420 return rep;
1421}
1422
1423G4bool
1424HepPolyhedron::GetNextEdge(G4Point3D &p1, G4Point3D &p2,
1425 G4int &edgeFlag, G4int &iface1, G4int &iface2) const
1426/***********************************************************************
1427 * *
1428 * Name: HepPolyhedron::GetNextEdge Date: 17.11.99 *
1429 * Author: E.Chernyaev Revised: *
1430 * *
1431 * Function: Get next edge with indices of the faces which share *
1432 * the edge. *
1433 * Returns false when the last edge. *
1434 * *
1435 ***********************************************************************/
1436{
1437 G4int i1,i2;
1438 G4bool rep = GetNextEdgeIndices(i1,i2,edgeFlag,iface1,iface2);
1439 p1 = pV[i1];
1440 p2 = pV[i2];
1441 return rep;
1442}
1443
1444void HepPolyhedron::GetFacet(G4int iFace, G4int &n, G4int *iNodes,
1445 G4int *edgeFlags, G4int *iFaces) const
1446/***********************************************************************
1447 * *
1448 * Name: HepPolyhedron::GetFacet Date: 15.12.99 *
1449 * Author: E.Chernyaev Revised: *
1450 * *
1451 * Function: Get face by index *
1452 * *
1453 ***********************************************************************/
1454{
1455 if (iFace < 1 || iFace > nface) {
1456 std::cerr
1457 << "HepPolyhedron::GetFacet: irrelevant index " << iFace
1458 << std::endl;
1459 n = 0;
1460 }else{
1461 G4int i, k;
1462 for (i=0; i<4; i++) {
1463 k = pF[iFace].edge[i].v;
1464 if (k == 0) break;
1465 if (iFaces != 0) iFaces[i] = pF[iFace].edge[i].f;
1466 if (k > 0) {
1467 iNodes[i] = k;
1468 if (edgeFlags != 0) edgeFlags[i] = 1;
1469 }else{
1470 iNodes[i] = -k;
1471 if (edgeFlags != 0) edgeFlags[i] = -1;
1472 }
1473 }
1474 n = i;
1475 }
1476}
1477
1478void HepPolyhedron::GetFacet(G4int index, G4int &n, G4Point3D *nodes,
1479 G4int *edgeFlags, G4Normal3D *normals) const
1480/***********************************************************************
1481 * *
1482 * Name: HepPolyhedron::GetFacet Date: 17.11.99 *
1483 * Author: E.Chernyaev Revised: *
1484 * *
1485 * Function: Get face by index *
1486 * *
1487 ***********************************************************************/
1488{
1489 G4int iNodes[4];
1490 GetFacet(index, n, iNodes, edgeFlags);
1491 if (n != 0) {
1492 for (G4int i=0; i<n; i++) {
1493 nodes[i] = pV[iNodes[i]];
1494 if (normals != 0) normals[i] = FindNodeNormal(index,iNodes[i]);
1495 }
1496 }
1497}
1498
1499G4bool
1500HepPolyhedron::GetNextFacet(G4int &n, G4Point3D *nodes,
1501 G4int *edgeFlags, G4Normal3D *normals) const
1502/***********************************************************************
1503 * *
1504 * Name: HepPolyhedron::GetNextFacet Date: 19.11.99 *
1505 * Author: E.Chernyaev Revised: *
1506 * *
1507 * Function: Get next face with normals of unit length at the nodes. *
1508 * Returns false when finished all faces. *
1509 * *
1510 ***********************************************************************/
1511{
1512 static G4ThreadLocal G4int iFace = 1;
1513
1514 if (edgeFlags == 0) {
1515 GetFacet(iFace, n, nodes);
1516 }else if (normals == 0) {
1517 GetFacet(iFace, n, nodes, edgeFlags);
1518 }else{
1519 GetFacet(iFace, n, nodes, edgeFlags, normals);
1520 }
1521
1522 if (++iFace > nface) {
1523 iFace = 1;
1524 return false;
1525 }else{
1526 return true;
1527 }
1528}
1529
1530G4Normal3D HepPolyhedron::GetNormal(G4int iFace) const
1531/***********************************************************************
1532 * *
1533 * Name: HepPolyhedron::GetNormal Date: 19.11.99 *
1534 * Author: E.Chernyaev Revised: *
1535 * *
1536 * Function: Get normal of the face given by index *
1537 * *
1538 ***********************************************************************/
1539{
1540 if (iFace < 1 || iFace > nface) {
1541 std::cerr
1542 << "HepPolyhedron::GetNormal: irrelevant index " << iFace
1543 << std::endl;
1544 return G4Normal3D();
1545 }
1546
1547 G4int i0 = std::abs(pF[iFace].edge[0].v);
1548 G4int i1 = std::abs(pF[iFace].edge[1].v);
1549 G4int i2 = std::abs(pF[iFace].edge[2].v);
1550 G4int i3 = std::abs(pF[iFace].edge[3].v);
1551 if (i3 == 0) i3 = i0;
1552 return (pV[i2] - pV[i0]).cross(pV[i3] - pV[i1]);
1553}
1554
1555G4Normal3D HepPolyhedron::GetUnitNormal(G4int iFace) const
1556/***********************************************************************
1557 * *
1558 * Name: HepPolyhedron::GetNormal Date: 19.11.99 *
1559 * Author: E.Chernyaev Revised: *
1560 * *
1561 * Function: Get unit normal of the face given by index *
1562 * *
1563 ***********************************************************************/
1564{
1565 if (iFace < 1 || iFace > nface) {
1566 std::cerr
1567 << "HepPolyhedron::GetUnitNormal: irrelevant index " << iFace
1568 << std::endl;
1569 return G4Normal3D();
1570 }
1571
1572 G4int i0 = std::abs(pF[iFace].edge[0].v);
1573 G4int i1 = std::abs(pF[iFace].edge[1].v);
1574 G4int i2 = std::abs(pF[iFace].edge[2].v);
1575 G4int i3 = std::abs(pF[iFace].edge[3].v);
1576 if (i3 == 0) i3 = i0;
1577 return ((pV[i2] - pV[i0]).cross(pV[i3] - pV[i1])).unit();
1578}
1579
1580G4bool HepPolyhedron::GetNextNormal(G4Normal3D &normal) const
1581/***********************************************************************
1582 * *
1583 * Name: HepPolyhedron::GetNextNormal Date: 22.07.96 *
1584 * Author: John Allison Revised: 19.11.99 *
1585 * *
1586 * Function: Get normals of each face in face order. Returns false *
1587 * when finished all faces. *
1588 * *
1589 ***********************************************************************/
1590{
1591 static G4ThreadLocal G4int iFace = 1;
1592 normal = GetNormal(iFace);
1593 if (++iFace > nface) {
1594 iFace = 1;
1595 return false;
1596 }else{
1597 return true;
1598 }
1599}
1600
1601G4bool HepPolyhedron::GetNextUnitNormal(G4Normal3D &normal) const
1602/***********************************************************************
1603 * *
1604 * Name: HepPolyhedron::GetNextUnitNormal Date: 16.09.96 *
1605 * Author: E.Chernyaev Revised: *
1606 * *
1607 * Function: Get normals of unit length of each face in face order. *
1608 * Returns false when finished all faces. *
1609 * *
1610 ***********************************************************************/
1611{
1612 G4bool rep = GetNextNormal(normal);
1613 normal = normal.unit();
1614 return rep;
1615}
1616
1617G4double HepPolyhedron::GetSurfaceArea() const
1618/***********************************************************************
1619 * *
1620 * Name: HepPolyhedron::GetSurfaceArea Date: 25.05.01 *
1621 * Author: E.Chernyaev Revised: *
1622 * *
1623 * Function: Returns area of the surface of the polyhedron. *
1624 * *
1625 ***********************************************************************/
1626{
1627 G4double srf = 0.;
1628 for (G4int iFace=1; iFace<=nface; iFace++) {
1629 G4int i0 = std::abs(pF[iFace].edge[0].v);
1630 G4int i1 = std::abs(pF[iFace].edge[1].v);
1631 G4int i2 = std::abs(pF[iFace].edge[2].v);
1632 G4int i3 = std::abs(pF[iFace].edge[3].v);
1633 if (i3 == 0) i3 = i0;
1634 srf += ((pV[i2] - pV[i0]).cross(pV[i3] - pV[i1])).mag();
1635 }
1636 return srf/2.;
1637}
1638
1639G4double HepPolyhedron::GetVolume() const
1640/***********************************************************************
1641 * *
1642 * Name: HepPolyhedron::GetVolume Date: 25.05.01 *
1643 * Author: E.Chernyaev Revised: *
1644 * *
1645 * Function: Returns volume of the polyhedron. *
1646 * *
1647 ***********************************************************************/
1648{
1649 G4double v = 0.;
1650 for (G4int iFace=1; iFace<=nface; iFace++) {
1651 G4int i0 = std::abs(pF[iFace].edge[0].v);
1652 G4int i1 = std::abs(pF[iFace].edge[1].v);
1653 G4int i2 = std::abs(pF[iFace].edge[2].v);
1654 G4int i3 = std::abs(pF[iFace].edge[3].v);
1655 G4Point3D pt;
1656 if (i3 == 0) {
1657 i3 = i0;
1658 pt = (pV[i0]+pV[i1]+pV[i2]) * (1./3.);
1659 }else{
1660 pt = (pV[i0]+pV[i1]+pV[i2]+pV[i3]) * 0.25;
1661 }
1662 v += ((pV[i2] - pV[i0]).cross(pV[i3] - pV[i1])).dot(pt);
1663 }
1664 return v/6.;
1665}
1666
1667G4int
1668HepPolyhedron::createTwistedTrap(G4double Dz,
1669 const G4double xy1[][2],
1670 const G4double xy2[][2])
1671/***********************************************************************
1672 * *
1673 * Name: createTwistedTrap Date: 05.11.02 *
1674 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
1675 * *
1676 * Function: Creates polyhedron for twisted trapezoid *
1677 * *
1678 * Input: Dz - half-length along Z 8----7 *
1679 * xy1[2,4] - quadrilateral at Z=-Dz 5----6 ! *
1680 * xy2[2,4] - quadrilateral at Z=+Dz ! 4-!--3 *
1681 * 1----2 *
1682 * *
1683 ***********************************************************************/
1684{
1685 AllocateMemory(12,18);
1686
1687 pV[ 1] = G4Point3D(xy1[0][0],xy1[0][1],-Dz);
1688 pV[ 2] = G4Point3D(xy1[1][0],xy1[1][1],-Dz);
1689 pV[ 3] = G4Point3D(xy1[2][0],xy1[2][1],-Dz);
1690 pV[ 4] = G4Point3D(xy1[3][0],xy1[3][1],-Dz);
1691
1692 pV[ 5] = G4Point3D(xy2[0][0],xy2[0][1], Dz);
1693 pV[ 6] = G4Point3D(xy2[1][0],xy2[1][1], Dz);
1694 pV[ 7] = G4Point3D(xy2[2][0],xy2[2][1], Dz);
1695 pV[ 8] = G4Point3D(xy2[3][0],xy2[3][1], Dz);
1696
1697 pV[ 9] = (pV[1]+pV[2]+pV[5]+pV[6])/4.;
1698 pV[10] = (pV[2]+pV[3]+pV[6]+pV[7])/4.;
1699 pV[11] = (pV[3]+pV[4]+pV[7]+pV[8])/4.;
1700 pV[12] = (pV[4]+pV[1]+pV[8]+pV[5])/4.;
1701
1702 enum {DUMMY, BOTTOM,
1703 LEFT_BOTTOM, LEFT_FRONT, LEFT_TOP, LEFT_BACK,
1704 BACK_BOTTOM, BACK_LEFT, BACK_TOP, BACK_RIGHT,
1705 RIGHT_BOTTOM, RIGHT_BACK, RIGHT_TOP, RIGHT_FRONT,
1706 FRONT_BOTTOM, FRONT_RIGHT, FRONT_TOP, FRONT_LEFT,
1707 TOP};
1708
1709 pF[ 1]=G4Facet(1,LEFT_BOTTOM, 4,BACK_BOTTOM, 3,RIGHT_BOTTOM, 2,FRONT_BOTTOM);
1710
1711 pF[ 2]=G4Facet(4,BOTTOM, -1,LEFT_FRONT, -12,LEFT_BACK, 0,0);
1712 pF[ 3]=G4Facet(1,FRONT_LEFT, -5,LEFT_TOP, -12,LEFT_BOTTOM, 0,0);
1713 pF[ 4]=G4Facet(5,TOP, -8,LEFT_BACK, -12,LEFT_FRONT, 0,0);
1714 pF[ 5]=G4Facet(8,BACK_LEFT, -4,LEFT_BOTTOM, -12,LEFT_TOP, 0,0);
1715
1716 pF[ 6]=G4Facet(3,BOTTOM, -4,BACK_LEFT, -11,BACK_RIGHT, 0,0);
1717 pF[ 7]=G4Facet(4,LEFT_BACK, -8,BACK_TOP, -11,BACK_BOTTOM, 0,0);
1718 pF[ 8]=G4Facet(8,TOP, -7,BACK_RIGHT, -11,BACK_LEFT, 0,0);
1719 pF[ 9]=G4Facet(7,RIGHT_BACK, -3,BACK_BOTTOM, -11,BACK_TOP, 0,0);
1720
1721 pF[10]=G4Facet(2,BOTTOM, -3,RIGHT_BACK, -10,RIGHT_FRONT, 0,0);
1722 pF[11]=G4Facet(3,BACK_RIGHT, -7,RIGHT_TOP, -10,RIGHT_BOTTOM, 0,0);
1723 pF[12]=G4Facet(7,TOP, -6,RIGHT_FRONT, -10,RIGHT_BACK, 0,0);
1724 pF[13]=G4Facet(6,FRONT_RIGHT,-2,RIGHT_BOTTOM,-10,RIGHT_TOP, 0,0);
1725
1726 pF[14]=G4Facet(1,BOTTOM, -2,FRONT_RIGHT, -9,FRONT_LEFT, 0,0);
1727 pF[15]=G4Facet(2,RIGHT_FRONT,-6,FRONT_TOP, -9,FRONT_BOTTOM, 0,0);
1728 pF[16]=G4Facet(6,TOP, -5,FRONT_LEFT, -9,FRONT_RIGHT, 0,0);
1729 pF[17]=G4Facet(5,LEFT_FRONT, -1,FRONT_BOTTOM, -9,FRONT_TOP, 0,0);
1730
1731 pF[18]=G4Facet(5,FRONT_TOP, 6,RIGHT_TOP, 7,BACK_TOP, 8,LEFT_TOP);
1732
1733 return 0;
1734}
1735
1736G4int
1737HepPolyhedron::createPolyhedron(G4int Nnodes, G4int Nfaces,
1738 const G4double xyz[][3],
1739 const G4int faces[][4])
1740/***********************************************************************
1741 * *
1742 * Name: createPolyhedron Date: 05.11.02 *
1743 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
1744 * *
1745 * Function: Creates user defined polyhedron *
1746 * *
1747 * Input: Nnodes - number of nodes *
1748 * Nfaces - number of faces *
1749 * nodes[][3] - node coordinates *
1750 * faces[][4] - faces *
1751 * *
1752 ***********************************************************************/
1753{
1754 AllocateMemory(Nnodes, Nfaces);
1755 if (nvert == 0) return 1;
1756
1757 for (G4int i=0; i<Nnodes; i++) {
1758 pV[i+1] = G4Point3D(xyz[i][0], xyz[i][1], xyz[i][2]);
1759 }
1760 for (G4int k=0; k<Nfaces; k++) {
1761 pF[k+1] = G4Facet(faces[k][0],0,faces[k][1],0,faces[k][2],0,faces[k][3],0);
1762 }
1763 SetReferences();
1764 return 0;
1765}
1766
1767HepPolyhedronTrd2::HepPolyhedronTrd2(G4double Dx1, G4double Dx2,
1768 G4double Dy1, G4double Dy2,
1769 G4double Dz)
1770/***********************************************************************
1771 * *
1772 * Name: HepPolyhedronTrd2 Date: 22.07.96 *
1773 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
1774 * *
1775 * Function: Create GEANT4 TRD2-trapezoid *
1776 * *
1777 * Input: Dx1 - half-length along X at -Dz 8----7 *
1778 * Dx2 - half-length along X ay +Dz 5----6 ! *
1779 * Dy1 - half-length along Y ay -Dz ! 4-!--3 *
1780 * Dy2 - half-length along Y ay +Dz 1----2 *
1781 * Dz - half-length along Z *
1782 * *
1783 ***********************************************************************/
1784{
1785 AllocateMemory(8,6);
1786
1787 pV[1] = G4Point3D(-Dx1,-Dy1,-Dz);
1788 pV[2] = G4Point3D( Dx1,-Dy1,-Dz);
1789 pV[3] = G4Point3D( Dx1, Dy1,-Dz);
1790 pV[4] = G4Point3D(-Dx1, Dy1,-Dz);
1791 pV[5] = G4Point3D(-Dx2,-Dy2, Dz);
1792 pV[6] = G4Point3D( Dx2,-Dy2, Dz);
1793 pV[7] = G4Point3D( Dx2, Dy2, Dz);
1794 pV[8] = G4Point3D(-Dx2, Dy2, Dz);
1795
1796 CreatePrism();
1797}
1798
1799HepPolyhedronTrd2::~HepPolyhedronTrd2() {}
1800
1801HepPolyhedronTrd1::HepPolyhedronTrd1(G4double Dx1, G4double Dx2,
1802 G4double Dy, G4double Dz)
1803 : HepPolyhedronTrd2(Dx1, Dx2, Dy, Dy, Dz) {}
1804
1805HepPolyhedronTrd1::~HepPolyhedronTrd1() {}
1806
1807HepPolyhedronBox::HepPolyhedronBox(G4double Dx, G4double Dy, G4double Dz)
1808 : HepPolyhedronTrd2(Dx, Dx, Dy, Dy, Dz) {}
1809
1810HepPolyhedronBox::~HepPolyhedronBox() {}
1811
1812HepPolyhedronTrap::HepPolyhedronTrap(G4double Dz,
1813 G4double Theta,
1814 G4double Phi,
1815 G4double Dy1,
1816 G4double Dx1,
1817 G4double Dx2,
1818 G4double Alp1,
1819 G4double Dy2,
1820 G4double Dx3,
1821 G4double Dx4,
1822 G4double Alp2)
1823/***********************************************************************
1824 * *
1825 * Name: HepPolyhedronTrap Date: 20.11.96 *
1826 * Author: E.Chernyaev Revised: *
1827 * *
1828 * Function: Create GEANT4 TRAP-trapezoid *
1829 * *
1830 * Input: DZ - half-length in Z *
1831 * Theta,Phi - polar angles of the line joining centres of the *
1832 * faces at Z=-Dz and Z=+Dz *
1833 * Dy1 - half-length in Y of the face at Z=-Dz *
1834 * Dx1 - half-length in X of low edge of the face at Z=-Dz *
1835 * Dx2 - half-length in X of top edge of the face at Z=-Dz *
1836 * Alp1 - angle between Y-axis and the median joining top and *
1837 * low edges of the face at Z=-Dz *
1838 * Dy2 - half-length in Y of the face at Z=+Dz *
1839 * Dx3 - half-length in X of low edge of the face at Z=+Dz *
1840 * Dx4 - half-length in X of top edge of the face at Z=+Dz *
1841 * Alp2 - angle between Y-axis and the median joining top and *
1842 * low edges of the face at Z=+Dz *
1843 * *
1844 ***********************************************************************/
1845{
1846 G4double DzTthetaCphi = Dz*std::tan(Theta)*std::cos(Phi);
1847 G4double DzTthetaSphi = Dz*std::tan(Theta)*std::sin(Phi);
1848 G4double Dy1Talp1 = Dy1*std::tan(Alp1);
1849 G4double Dy2Talp2 = Dy2*std::tan(Alp2);
1850
1851 AllocateMemory(8,6);
1852
1853 pV[1] = G4Point3D(-DzTthetaCphi-Dy1Talp1-Dx1,-DzTthetaSphi-Dy1,-Dz);
1854 pV[2] = G4Point3D(-DzTthetaCphi-Dy1Talp1+Dx1,-DzTthetaSphi-Dy1,-Dz);
1855 pV[3] = G4Point3D(-DzTthetaCphi+Dy1Talp1+Dx2,-DzTthetaSphi+Dy1,-Dz);
1856 pV[4] = G4Point3D(-DzTthetaCphi+Dy1Talp1-Dx2,-DzTthetaSphi+Dy1,-Dz);
1857 pV[5] = G4Point3D( DzTthetaCphi-Dy2Talp2-Dx3, DzTthetaSphi-Dy2, Dz);
1858 pV[6] = G4Point3D( DzTthetaCphi-Dy2Talp2+Dx3, DzTthetaSphi-Dy2, Dz);
1859 pV[7] = G4Point3D( DzTthetaCphi+Dy2Talp2+Dx4, DzTthetaSphi+Dy2, Dz);
1860 pV[8] = G4Point3D( DzTthetaCphi+Dy2Talp2-Dx4, DzTthetaSphi+Dy2, Dz);
1861
1862 CreatePrism();
1863}
1864
1865HepPolyhedronTrap::~HepPolyhedronTrap() {}
1866
1867HepPolyhedronPara::HepPolyhedronPara(G4double Dx, G4double Dy, G4double Dz,
1868 G4double Alpha, G4double Theta,
1869 G4double Phi)
1870 : HepPolyhedronTrap(Dz, Theta, Phi, Dy, Dx, Dx, Alpha, Dy, Dx, Dx, Alpha) {}
1871
1872HepPolyhedronPara::~HepPolyhedronPara() {}
1873
1874HepPolyhedronParaboloid::HepPolyhedronParaboloid(G4double r1,
1875 G4double r2,
1876 G4double dz,
1877 G4double sPhi,
1878 G4double dPhi)
1879/***********************************************************************
1880 * *
1881 * Name: HepPolyhedronParaboloid Date: 28.06.07 *
1882 * Author: L.Lindroos, T.Nikitina (CERN), July 2007 Revised: 28.06.07 *
1883 * *
1884 * Function: Constructor for paraboloid *
1885 * *
1886 * Input: r1 - inside and outside radiuses at -Dz *
1887 * r2 - inside and outside radiuses at +Dz *
1888 * dz - half length in Z *
1889 * sPhi - starting angle of the segment *
1890 * dPhi - segment range *
1891 * *
1892 ***********************************************************************/
1893{
1894 static const G4double wholeCircle=twopi;
1895
1896 // C H E C K I N P U T P A R A M E T E R S
1897
1898 G4int k = 0;
1899 if (r1 < 0. || r2 <= 0.) k = 1;
1900
1901 if (dz <= 0.) k += 2;
1902
1903 G4double phi1, phi2, dphi;
1904
1905 if(dPhi < 0.)
1906 {
1907 phi2 = sPhi; phi1 = phi2 + dPhi;
1908 }
1909 else if(dPhi == 0.)
1910 {
1911 phi1 = sPhi; phi2 = phi1 + wholeCircle;
1912 }
1913 else
1914 {
1915 phi1 = sPhi; phi2 = phi1 + dPhi;
1916 }
1917 dphi = phi2 - phi1;
1918
1919 if (std::abs(dphi-wholeCircle) < perMillion) dphi = wholeCircle;
1920 if (dphi > wholeCircle) k += 4;
1921
1922 if (k != 0) {
1923 std::cerr << "HepPolyhedronParaboloid: error in input parameters";
1924 if ((k & 1) != 0) std::cerr << " (radiuses)";
1925 if ((k & 2) != 0) std::cerr << " (half-length)";
1926 if ((k & 4) != 0) std::cerr << " (angles)";
1927 std::cerr << std::endl;
1928 std::cerr << " r1=" << r1;
1929 std::cerr << " r2=" << r2;
1930 std::cerr << " dz=" << dz << " sPhi=" << sPhi << " dPhi=" << dPhi
1931 << std::endl;
1932 return;
1933 }
1934
1935 // P R E P A R E T W O P O L Y L I N E S
1936
1937 G4int n = GetNumberOfRotationSteps();
1938 G4double dl = (r2 - r1) / n;
1939 G4double k1 = (r2*r2 - r1*r1) / 2 / dz;
1940 G4double k2 = (r2*r2 + r1*r1) / 2;
1941
1942 G4double *zz = new G4double[n + 2], *rr = new G4double[n + 2];
1943
1944 zz[0] = dz;
1945 rr[0] = r2;
1946
1947 for(G4int i = 1; i < n - 1; i++)
1948 {
1949 rr[i] = rr[i-1] - dl;
1950 zz[i] = (rr[i]*rr[i] - k2) / k1;
1951 if(rr[i] < 0)
1952 {
1953 rr[i] = 0;
1954 zz[i] = 0;
1955 }
1956 }
1957
1958 zz[n-1] = -dz;
1959 rr[n-1] = r1;
1960
1961 zz[n] = dz;
1962 rr[n] = 0;
1963
1964 zz[n+1] = -dz;
1965 rr[n+1] = 0;
1966
1967 // R O T A T E P O L Y L I N E S
1968
1969 RotateAroundZ(0, phi1, dphi, n, 2, zz, rr, -1, -1);
1970 SetReferences();
1971
1972 delete [] zz;
1973 delete [] rr;
1974}
1975
1976HepPolyhedronParaboloid::~HepPolyhedronParaboloid() {}
1977
1978HepPolyhedronHype::HepPolyhedronHype(G4double r1,
1979 G4double r2,
1980 G4double sqrtan1,
1981 G4double sqrtan2,
1982 G4double halfZ)
1983/***********************************************************************
1984 * *
1985 * Name: HepPolyhedronHype Date: 14.04.08 *
1986 * Author: Tatiana Nikitina (CERN) Revised: 14.04.08 *
1987 * Evgueni Tcherniaev 01.12.20 *
1988 * *
1989 * Function: Constructor for Hype *
1990 * *
1991 * Input: r1 - inside radius at z=0 *
1992 * r2 - outside radiuses at z=0 *
1993 * sqrtan1 - sqr of tan of Inner Stereo Angle *
1994 * sqrtan2 - sqr of tan of Outer Stereo Angle *
1995 * halfZ - half length in Z *
1996 * *
1997 ***********************************************************************/
1998{
1999 static const G4double wholeCircle = twopi;
2000
2001 // C H E C K I N P U T P A R A M E T E R S
2002
2003 G4int k = 0;
2004 if (r1 < 0. || r2 < 0. || r1 >= r2) k = 1;
2005 if (halfZ <= 0.) k += 2;
2006 if (sqrtan1 < 0.|| sqrtan2 < 0.) k += 4;
2007
2008 if (k != 0)
2009 {
2010 std::cerr << "HepPolyhedronHype: error in input parameters";
2011 if ((k & 1) != 0) std::cerr << " (radiuses)";
2012 if ((k & 2) != 0) std::cerr << " (half-length)";
2013 if ((k & 4) != 0) std::cerr << " (angles)";
2014 std::cerr << std::endl;
2015 std::cerr << " r1=" << r1 << " r2=" << r2;
2016 std::cerr << " halfZ=" << halfZ << " sqrTan1=" << sqrtan1
2017 << " sqrTan2=" << sqrtan2
2018 << std::endl;
2019 return;
2020 }
2021
2022 // P R E P A R E T W O P O L Y L I N E S
2023
2024 G4int ns = std::max(3, GetNumberOfRotationSteps()/4);
2025 G4int nz1 = (sqrtan1 == 0.) ? 2 : ns + 1;
2026 G4int nz2 = (sqrtan2 == 0.) ? 2 : ns + 1;
2027 G4double* zz = new G4double[nz1 + nz2];
2028 G4double* rr = new G4double[nz1 + nz2];
2029
2030 // external polyline
2031 G4double dz2 = 2.*halfZ/(nz2 - 1);
2032 for(G4int i = 0; i < nz2; ++i)
2033 {
2034 zz[i] = halfZ - dz2*i;
2035 rr[i] = std::sqrt(sqrtan2*zz[i]*zz[i] + r2*r2);
2036 }
2037
2038 // internal polyline
2039 G4double dz1 = 2.*halfZ/(nz1 - 1);
2040 for(G4int i = 0; i < nz1; ++i)
2041 {
2042 G4int j = nz2 + i;
2043 zz[j] = halfZ - dz1*i;
2044 rr[j] = std::sqrt(sqrtan1*zz[j]*zz[j] + r1*r1);
2045 }
2046
2047 // R O T A T E P O L Y L I N E S
2048
2049 RotateAroundZ(0, 0., wholeCircle, nz2, nz1, zz, rr, -1, -1);
2050 SetReferences();
2051
2052 delete [] zz;
2053 delete [] rr;
2054}
2055
2056HepPolyhedronHype::~HepPolyhedronHype() {}
2057
2058HepPolyhedronCons::HepPolyhedronCons(G4double Rmn1,
2059 G4double Rmx1,
2060 G4double Rmn2,
2061 G4double Rmx2,
2062 G4double Dz,
2063 G4double Phi1,
2064 G4double Dphi)
2065/***********************************************************************
2066 * *
2067 * Name: HepPolyhedronCons::HepPolyhedronCons Date: 15.12.96 *
2068 * Author: E.Chernyaev (IHEP/Protvino) Revised: 15.12.96 *
2069 * *
2070 * Function: Constructor for CONS, TUBS, CONE, TUBE *
2071 * *
2072 * Input: Rmn1, Rmx1 - inside and outside radiuses at -Dz *
2073 * Rmn2, Rmx2 - inside and outside radiuses at +Dz *
2074 * Dz - half length in Z *
2075 * Phi1 - starting angle of the segment *
2076 * Dphi - segment range *
2077 * *
2078 ***********************************************************************/
2079{
2080 static const G4double wholeCircle=twopi;
2081
2082 // C H E C K I N P U T P A R A M E T E R S
2083
2084 G4int k = 0;
2085 if (Rmn1 < 0. || Rmx1 < 0. || Rmn2 < 0. || Rmx2 < 0.) k = 1;
2086 if (Rmn1 > Rmx1 || Rmn2 > Rmx2) k = 1;
2087 if (Rmn1 == Rmx1 && Rmn2 == Rmx2) k = 1;
2088
2089 if (Dz <= 0.) k += 2;
2090
2091 G4double phi1, phi2, dphi;
2092 if (Dphi < 0.) {
2093 phi2 = Phi1; phi1 = phi2 - Dphi;
2094 }else if (Dphi == 0.) {
2095 phi1 = Phi1; phi2 = phi1 + wholeCircle;
2096 }else{
2097 phi1 = Phi1; phi2 = phi1 + Dphi;
2098 }
2099 dphi = phi2 - phi1;
2100 if (std::abs(dphi-wholeCircle) < perMillion) dphi = wholeCircle;
2101 if (dphi > wholeCircle) k += 4;
2102
2103 if (k != 0) {
2104 std::cerr << "HepPolyhedronCone(s)/Tube(s): error in input parameters";
2105 if ((k & 1) != 0) std::cerr << " (radiuses)";
2106 if ((k & 2) != 0) std::cerr << " (half-length)";
2107 if ((k & 4) != 0) std::cerr << " (angles)";
2108 std::cerr << std::endl;
2109 std::cerr << " Rmn1=" << Rmn1 << " Rmx1=" << Rmx1;
2110 std::cerr << " Rmn2=" << Rmn2 << " Rmx2=" << Rmx2;
2111 std::cerr << " Dz=" << Dz << " Phi1=" << Phi1 << " Dphi=" << Dphi
2112 << std::endl;
2113 return;
2114 }
2115
2116 // P R E P A R E T W O P O L Y L I N E S
2117
2118 G4double zz[4], rr[4];
2119 zz[0] = Dz;
2120 zz[1] = -Dz;
2121 zz[2] = Dz;
2122 zz[3] = -Dz;
2123 rr[0] = Rmx2;
2124 rr[1] = Rmx1;
2125 rr[2] = Rmn2;
2126 rr[3] = Rmn1;
2127
2128 // R O T A T E P O L Y L I N E S
2129
2130 RotateAroundZ(0, phi1, dphi, 2, 2, zz, rr, -1, -1);
2131 SetReferences();
2132}
2133
2134HepPolyhedronCons::~HepPolyhedronCons() {}
2135
2136HepPolyhedronCone::HepPolyhedronCone(G4double Rmn1, G4double Rmx1,
2137 G4double Rmn2, G4double Rmx2,
2138 G4double Dz) :
2139 HepPolyhedronCons(Rmn1, Rmx1, Rmn2, Rmx2, Dz, 0*deg, 360*deg) {}
2140
2141HepPolyhedronCone::~HepPolyhedronCone() {}
2142
2143HepPolyhedronTubs::HepPolyhedronTubs(G4double Rmin, G4double Rmax,
2144 G4double Dz,
2145 G4double Phi1, G4double Dphi)
2146 : HepPolyhedronCons(Rmin, Rmax, Rmin, Rmax, Dz, Phi1, Dphi) {}
2147
2148HepPolyhedronTubs::~HepPolyhedronTubs() {}
2149
2150HepPolyhedronTube::HepPolyhedronTube (G4double Rmin, G4double Rmax,
2151 G4double Dz)
2152 : HepPolyhedronCons(Rmin, Rmax, Rmin, Rmax, Dz, 0*deg, 360*deg) {}
2153
2154HepPolyhedronTube::~HepPolyhedronTube () {}
2155
2156HepPolyhedronPgon::HepPolyhedronPgon(G4double phi,
2157 G4double dphi,
2158 G4int npdv,
2159 G4int nz,
2160 const G4double *z,
2161 const G4double *rmin,
2162 const G4double *rmax)
2163/***********************************************************************
2164 * *
2165 * Name: HepPolyhedronPgon Date: 09.12.96 *
2166 * Author: E.Chernyaev Revised: *
2167 * *
2168 * Function: Constructor of polyhedron for PGON, PCON *
2169 * *
2170 * Input: phi - initial phi *
2171 * dphi - delta phi *
2172 * npdv - number of steps along phi *
2173 * nz - number of z-planes (at least two) *
2174 * z[] - z coordinates of the slices *
2175 * rmin[] - smaller r at the slices *
2176 * rmax[] - bigger r at the slices *
2177 * *
2178 ***********************************************************************/
2179{
2180 // C H E C K I N P U T P A R A M E T E R S
2181
2182 if (dphi <= 0. || dphi > twopi) {
2183 std::cerr
2184 << "HepPolyhedronPgon/Pcon: wrong delta phi = " << dphi
2185 << std::endl;
2186 return;
2187 }
2188
2189 if (nz < 2) {
2190 std::cerr
2191 << "HepPolyhedronPgon/Pcon: number of z-planes less than two = " << nz
2192 << std::endl;
2193 return;
2194 }
2195
2196 if (npdv < 0) {
2197 std::cerr
2198 << "HepPolyhedronPgon/Pcon: error in number of phi-steps =" << npdv
2199 << std::endl;
2200 return;
2201 }
2202
2203 G4int i;
2204 for (i=0; i<nz; i++) {
2205 if (rmin[i] < 0. || rmax[i] < 0. || rmin[i] > rmax[i]) {
2206 std::cerr
2207 << "HepPolyhedronPgon: error in radiuses rmin[" << i << "]="
2208 << rmin[i] << " rmax[" << i << "]=" << rmax[i]
2209 << std::endl;
2210 return;
2211 }
2212 }
2213
2214 // P R E P A R E T W O P O L Y L I N E S
2215
2216 G4double *zz, *rr;
2217 zz = new G4double[2*nz];
2218 rr = new G4double[2*nz];
2219
2220 if (z[0] > z[nz-1]) {
2221 for (i=0; i<nz; i++) {
2222 zz[i] = z[i];
2223 rr[i] = rmax[i];
2224 zz[i+nz] = z[i];
2225 rr[i+nz] = rmin[i];
2226 }
2227 }else{
2228 for (i=0; i<nz; i++) {
2229 zz[i] = z[nz-i-1];
2230 rr[i] = rmax[nz-i-1];
2231 zz[i+nz] = z[nz-i-1];
2232 rr[i+nz] = rmin[nz-i-1];
2233 }
2234 }
2235
2236 // R O T A T E P O L Y L I N E S
2237
2238 G4int nodeVis = 1;
2239 G4int edgeVis = (npdv == 0) ? -1 : 1;
2240 RotateAroundZ(npdv, phi, dphi, nz, nz, zz, rr, nodeVis, edgeVis);
2241 SetReferences();
2242
2243 delete [] zz;
2244 delete [] rr;
2245}
2246
2247HepPolyhedronPgon::HepPolyhedronPgon(G4double phi,
2248 G4double dphi,
2249 G4int npdv,
2250 const std::vector<G4TwoVector> &rz)
2251/***********************************************************************
2252 * *
2253 * Name: HepPolyhedronPgon Date: 12.05.21 *
2254 * Author: E.Tcherniaev (E.Chernyaev) Revised: *
2255 * *
2256 * Function: Constructor of polyhedron for PGON, PCON *
2257 * *
2258 * Input: phi - initial phi *
2259 * dphi - delta phi *
2260 * npdv - number of steps along phi *
2261 * rz - rz-contour *
2262 * *
2263 ***********************************************************************/
2264{
2265 // C H E C K I N P U T P A R A M E T E R S
2266
2267 if (dphi <= 0. || dphi > twopi) {
2268 std::cerr
2269 << "HepPolyhedronPgon/Pcon: wrong delta phi = " << dphi
2270 << std::endl;
2271 return;
2272 }
2273
2274 if (npdv < 0) {
2275 std::cerr
2276 << "HepPolyhedronPgon/Pcon: error in number of phi-steps = " << npdv
2277 << std::endl;
2278 return;
2279 }
2280
2281 G4int nrz = rz.size();
2282 if (nrz < 3) {
2283 std::cerr
2284 << "HepPolyhedronPgon/Pcon: invalid number of nodes in rz-contour = " << nrz
2285 << std::endl;
2286 return;
2287 }
2288
2289 // R O T A T E P O L Y L I N E
2290
2291 G4int nodeVis = 1;
2292 G4int edgeVis = (npdv == 0) ? -1 : 1;
2293 RotateContourAroundZ(npdv, phi, dphi, rz, nodeVis, edgeVis);
2294 SetReferences();
2295}
2296
2297HepPolyhedronPgon::~HepPolyhedronPgon() {}
2298
2299HepPolyhedronPcon::HepPolyhedronPcon(G4double phi, G4double dphi, G4int nz,
2300 const G4double *z,
2301 const G4double *rmin,
2302 const G4double *rmax)
2303 : HepPolyhedronPgon(phi, dphi, 0, nz, z, rmin, rmax) {}
2304
2305HepPolyhedronPcon::HepPolyhedronPcon(G4double phi, G4double dphi,
2306 const std::vector<G4TwoVector> &rz)
2307 : HepPolyhedronPgon(phi, dphi, 0, rz) {}
2308
2309HepPolyhedronPcon::~HepPolyhedronPcon() {}
2310
2311HepPolyhedronSphere::HepPolyhedronSphere(G4double rmin, G4double rmax,
2312 G4double phi, G4double dphi,
2313 G4double the, G4double dthe)
2314/***********************************************************************
2315 * *
2316 * Name: HepPolyhedronSphere Date: 11.12.96 *
2317 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
2318 * *
2319 * Function: Constructor of polyhedron for SPHERE *
2320 * *
2321 * Input: rmin - internal radius *
2322 * rmax - external radius *
2323 * phi - initial phi *
2324 * dphi - delta phi *
2325 * the - initial theta *
2326 * dthe - delta theta *
2327 * *
2328 ***********************************************************************/
2329{
2330 // C H E C K I N P U T P A R A M E T E R S
2331
2332 if (dphi <= 0. || dphi > twopi) {
2333 std::cerr
2334 << "HepPolyhedronSphere: wrong delta phi = " << dphi
2335 << std::endl;
2336 return;
2337 }
2338
2339 if (the < 0. || the > pi) {
2340 std::cerr
2341 << "HepPolyhedronSphere: wrong theta = " << the
2342 << std::endl;
2343 return;
2344 }
2345
2346 if (dthe <= 0. || dthe > pi) {
2347 std::cerr
2348 << "HepPolyhedronSphere: wrong delta theta = " << dthe
2349 << std::endl;
2350 return;
2351 }
2352
2353 if (the+dthe > pi) {
2354 std::cerr
2355 << "HepPolyhedronSphere: wrong theta + delta theta = "
2356 << the << " " << dthe
2357 << std::endl;
2358 return;
2359 }
2360
2361 if (rmin < 0. || rmin >= rmax) {
2362 std::cerr
2363 << "HepPolyhedronSphere: error in radiuses"
2364 << " rmin=" << rmin << " rmax=" << rmax
2365 << std::endl;
2366 return;
2367 }
2368
2369 // P R E P A R E T W O P O L Y L I N E S
2370
2371 G4int nds = (GetNumberOfRotationSteps() + 1) / 2;
2372 G4int np1 = G4int(dthe*nds/pi+.5) + 1;
2373 if (np1 <= 1) np1 = 2;
2374 G4int np2 = rmin < spatialTolerance ? 1 : np1;
2375
2376 G4double *zz, *rr;
2377 zz = new G4double[np1+np2];
2378 rr = new G4double[np1+np2];
2379
2380 G4double a = dthe/(np1-1);
2381 G4double cosa, sina;
2382 for (G4int i=0; i<np1; i++) {
2383 cosa = std::cos(the+i*a);
2384 sina = std::sin(the+i*a);
2385 zz[i] = rmax*cosa;
2386 rr[i] = rmax*sina;
2387 if (np2 > 1) {
2388 zz[i+np1] = rmin*cosa;
2389 rr[i+np1] = rmin*sina;
2390 }
2391 }
2392 if (np2 == 1) {
2393 zz[np1] = 0.;
2394 rr[np1] = 0.;
2395 }
2396
2397 // R O T A T E P O L Y L I N E S
2398
2399 RotateAroundZ(0, phi, dphi, np1, np2, zz, rr, -1, -1);
2400 SetReferences();
2401
2402 delete [] zz;
2403 delete [] rr;
2404}
2405
2406HepPolyhedronSphere::~HepPolyhedronSphere() {}
2407
2408HepPolyhedronTorus::HepPolyhedronTorus(G4double rmin,
2409 G4double rmax,
2410 G4double rtor,
2411 G4double phi,
2412 G4double dphi)
2413/***********************************************************************
2414 * *
2415 * Name: HepPolyhedronTorus Date: 11.12.96 *
2416 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
2417 * *
2418 * Function: Constructor of polyhedron for TORUS *
2419 * *
2420 * Input: rmin - internal radius *
2421 * rmax - external radius *
2422 * rtor - radius of torus *
2423 * phi - initial phi *
2424 * dphi - delta phi *
2425 * *
2426 ***********************************************************************/
2427{
2428 // C H E C K I N P U T P A R A M E T E R S
2429
2430 if (dphi <= 0. || dphi > twopi) {
2431 std::cerr
2432 << "HepPolyhedronTorus: wrong delta phi = " << dphi
2433 << std::endl;
2434 return;
2435 }
2436
2437 if (rmin < 0. || rmin >= rmax || rmax >= rtor) {
2438 std::cerr
2439 << "HepPolyhedronTorus: error in radiuses"
2440 << " rmin=" << rmin << " rmax=" << rmax << " rtorus=" << rtor
2441 << std::endl;
2442 return;
2443 }
2444
2445 // P R E P A R E T W O P O L Y L I N E S
2446
2447 G4int np1 = GetNumberOfRotationSteps();
2448 G4int np2 = rmin < spatialTolerance ? 1 : np1;
2449
2450 G4double *zz, *rr;
2451 zz = new G4double[np1+np2];
2452 rr = new G4double[np1+np2];
2453
2454 G4double a = twopi/np1;
2455 G4double cosa, sina;
2456 for (G4int i=0; i<np1; i++) {
2457 cosa = std::cos(i*a);
2458 sina = std::sin(i*a);
2459 zz[i] = rmax*cosa;
2460 rr[i] = rtor+rmax*sina;
2461 if (np2 > 1) {
2462 zz[i+np1] = rmin*cosa;
2463 rr[i+np1] = rtor+rmin*sina;
2464 }
2465 }
2466 if (np2 == 1) {
2467 zz[np1] = 0.;
2468 rr[np1] = rtor;
2469 np2 = -1;
2470 }
2471
2472 // R O T A T E P O L Y L I N E S
2473
2474 RotateAroundZ(0, phi, dphi, -np1, -np2, zz, rr, -1,-1);
2475 SetReferences();
2476
2477 delete [] zz;
2478 delete [] rr;
2479}
2480
2481HepPolyhedronTorus::~HepPolyhedronTorus() {}
2482
2483HepPolyhedronTet::HepPolyhedronTet(const G4double p0[3],
2484 const G4double p1[3],
2485 const G4double p2[3],
2486 const G4double p3[3])
2487/***********************************************************************
2488 * *
2489 * Name: HepPolyhedronTet Date: 21.02.2020 *
2490 * Author: E.Tcherniaev (E.Chernyaev) Revised: *
2491 * *
2492 * Function: Constructor of polyhedron for TETrahedron *
2493 * *
2494 * Input: p0,p1,p2,p3 - vertices *
2495 * *
2496 ***********************************************************************/
2497{
2498 AllocateMemory(4,4);
2499
2500 pV[1].set(p0[0], p0[1], p0[2]);
2501 pV[2].set(p1[0], p1[1], p1[2]);
2502 pV[3].set(p2[0], p2[1], p2[2]);
2503 pV[4].set(p3[0], p3[1], p3[2]);
2504
2505 G4Vector3D v1(pV[2] - pV[1]);
2506 G4Vector3D v2(pV[3] - pV[1]);
2507 G4Vector3D v3(pV[4] - pV[1]);
2508
2509 if (v1.cross(v2).dot(v3) < 0.)
2510 {
2511 pV[3].set(p3[0], p3[1], p3[2]);
2512 pV[4].set(p2[0], p2[1], p2[2]);
2513 }
2514
2515 pF[1] = G4Facet(1,2, 3,4, 2,3);
2516 pF[2] = G4Facet(1,3, 4,4, 3,1);
2517 pF[3] = G4Facet(1,1, 2,4, 4,2);
2518 pF[4] = G4Facet(2,1, 3,2, 4,3);
2519}
2520
2521HepPolyhedronTet::~HepPolyhedronTet() {}
2522
2523HepPolyhedronEllipsoid::HepPolyhedronEllipsoid(G4double ax, G4double by,
2524 G4double cz, G4double zCut1,
2525 G4double zCut2)
2526/***********************************************************************
2527 * *
2528 * Name: HepPolyhedronEllipsoid Date: 25.02.05 *
2529 * Author: G.Guerrieri Revised: *
2530 * *
2531 * Function: Constructor of polyhedron for ELLIPSOID *
2532 * *
2533 * Input: ax - semiaxis x *
2534 * by - semiaxis y *
2535 * cz - semiaxis z *
2536 * zCut1 - lower cut plane level (solid lies above this plane) *
2537 * zCut2 - upper cut plane level (solid lies below this plane) *
2538 * *
2539 ***********************************************************************/
2540{
2541 // C H E C K I N P U T P A R A M E T E R S
2542
2543 if (zCut1 >= cz || zCut2 <= -cz || zCut1 > zCut2) {
2544 std::cerr << "HepPolyhedronEllipsoid: wrong zCut1 = " << zCut1
2545 << " zCut2 = " << zCut2
2546 << " for given cz = " << cz << std::endl;
2547 return;
2548 }
2549 if (cz <= 0.0) {
2550 std::cerr << "HepPolyhedronEllipsoid: bad z semi-axis: cz = " << cz
2551 << std::endl;
2552 return;
2553 }
2554
2555 G4double dthe;
2556 G4double sthe;
2557 G4int cutflag;
2558 cutflag= 0;
2559 if (zCut2 >= cz)
2560 {
2561 sthe= 0.0;
2562 }
2563 else
2564 {
2565 sthe= std::acos(zCut2/cz);
2566 cutflag++;
2567 }
2568 if (zCut1 <= -cz)
2569 {
2570 dthe= pi - sthe;
2571 }
2572 else
2573 {
2574 dthe= std::acos(zCut1/cz)-sthe;
2575 cutflag++;
2576 }
2577
2578 // P R E P A R E T W O P O L Y L I N E S
2579 // generate sphere of radius cz first, then rescale x and y later
2580
2581 G4int nds = (GetNumberOfRotationSteps() + 1) / 2;
2582 G4int np1 = G4int(dthe*nds/pi) + 2 + cutflag;
2583
2584 G4double *zz, *rr;
2585 zz = new G4double[np1+1];
2586 rr = new G4double[np1+1];
2587 if (!zz || !rr)
2588 {
2589 G4Exception("HepPolyhedronEllipsoid::HepPolyhedronEllipsoid",
2590 "greps1002", FatalException, "Out of memory");
2591 }
2592
2593 G4double a = dthe/(np1-cutflag-1);
2594 G4double cosa, sina;
2595 G4int j=0;
2596 if (sthe > 0.0)
2597 {
2598 zz[j]= zCut2;
2599 rr[j]= 0.;
2600 j++;
2601 }
2602 for (G4int i=0; i<np1-cutflag; i++) {
2603 cosa = std::cos(sthe+i*a);
2604 sina = std::sin(sthe+i*a);
2605 zz[j] = cz*cosa;
2606 rr[j] = cz*sina;
2607 j++;
2608 }
2609 if (j < np1)
2610 {
2611 zz[j]= zCut1;
2612 rr[j]= 0.;
2613 j++;
2614 }
2615 if (j > np1)
2616 {
2617 std::cerr << "Logic error in HepPolyhedronEllipsoid, memory corrupted!"
2618 << std::endl;
2619 }
2620 if (j < np1)
2621 {
2622 std::cerr << "Warning: logic error in HepPolyhedronEllipsoid."
2623 << std::endl;
2624 np1= j;
2625 }
2626 zz[j] = 0.;
2627 rr[j] = 0.;
2628
2629
2630 // R O T A T E P O L Y L I N E S
2631
2632 RotateAroundZ(0, 0.0, twopi, np1, 1, zz, rr, -1, 1);
2633 SetReferences();
2634
2635 delete [] zz;
2636 delete [] rr;
2637
2638 // rescale x and y vertex coordinates
2639 {
2640 G4Point3D * p= pV;
2641 for (G4int i=0; i<nvert; i++, p++) {
2642 p->setX( p->x() * ax/cz );
2643 p->setY( p->y() * by/cz );
2644 }
2645 }
2646}
2647
2648HepPolyhedronEllipsoid::~HepPolyhedronEllipsoid() {}
2649
2650HepPolyhedronEllipticalCone::HepPolyhedronEllipticalCone(G4double ax,
2651 G4double ay,
2652 G4double h,
2653 G4double zTopCut)
2654/***********************************************************************
2655 * *
2656 * Name: HepPolyhedronEllipticalCone Date: 8.9.2005 *
2657 * Author: D.Anninos Revised: 9.9.2005 *
2658 * *
2659 * Function: Constructor for EllipticalCone *
2660 * *
2661 * Input: ax, ay - X & Y semi axes at z = 0 *
2662 * h - height of full cone *
2663 * zTopCut - Top Cut in Z Axis *
2664 * *
2665 ***********************************************************************/
2666{
2667 // C H E C K I N P U T P A R A M E T E R S
2668
2669 G4int k = 0;
2670 if ( (ax <= 0.) || (ay <= 0.) || (h <= 0.) || (zTopCut <= 0.) ) { k = 1; }
2671
2672 if (k != 0) {
2673 std::cerr << "HepPolyhedronCone: error in input parameters";
2674 std::cerr << std::endl;
2675 return;
2676 }
2677
2678 // P R E P A R E T W O P O L Y L I N E S
2679
2680 zTopCut = (h >= zTopCut ? zTopCut : h);
2681
2682 G4double *zz, *rr;
2683 zz = new G4double[4];
2684 rr = new G4double[4];
2685 zz[0] = zTopCut;
2686 zz[1] = -zTopCut;
2687 zz[2] = zTopCut;
2688 zz[3] = -zTopCut;
2689 rr[0] = (h-zTopCut);
2690 rr[1] = (h+zTopCut);
2691 rr[2] = 0.;
2692 rr[3] = 0.;
2693
2694 // R O T A T E P O L Y L I N E S
2695
2696 RotateAroundZ(0, 0., twopi, 2, 2, zz, rr, -1, -1);
2697 SetReferences();
2698
2699 delete [] zz;
2700 delete [] rr;
2701
2702 // rescale x and y vertex coordinates
2703 {
2704 G4Point3D * p= pV;
2705 for (G4int i=0; i<nvert; i++, p++) {
2706 p->setX( p->x() * ax );
2707 p->setY( p->y() * ay );
2708 }
2709 }
2710}
2711
2712HepPolyhedronEllipticalCone::~HepPolyhedronEllipticalCone() {}
2713
2714HepPolyhedronHyperbolicMirror::HepPolyhedronHyperbolicMirror(G4double a,
2715 G4double h,
2716 G4double r)
2717/***********************************************************************
2718 * *
2719 * Name: HepPolyhedronHyperbolicMirror Date: 22.02.2020 *
2720 * Author: E.Tcherniaev (E.Chernyaev) Revised: *
2721 * *
2722 * Function: Create polyhedron for Hyperbolic mirror *
2723 * *
2724 * Input: a - half-separation *
2725 * h - height *
2726 * r - radius *
2727 * *
2728 ***********************************************************************/
2729{
2730 G4double H = std::abs(h);
2731 G4double R = std::abs(r);
2732 G4double A = std::abs(a);
2733 G4double B = A*R/std::sqrt(2*A*H + H*H);
2734
2735 // P R E P A R E T W O P O L Y L I N E S
2736
2737 G4int np1 = (A == 0.) ? 2 : std::max(3, GetNumberOfRotationSteps()/4) + 1;
2738 G4int np2 = 2;
2739 G4double maxAng = (A == 0.) ? 0. : std::acosh(1. + H/A);
2740 G4double delAng = maxAng/(np1 - 1);
2741
2742 G4double *zz = new G4double[np1 + np2];
2743 G4double *rr = new G4double[np1 + np2];
2744
2745 // 1st polyline
2746 zz[0] = H;
2747 rr[0] = R;
2748 for (G4int iz = 1; iz < np1 - 1; ++iz)
2749 {
2750 G4double ang = maxAng - iz*delAng;
2751 zz[iz] = A*std::cosh(ang) - A;
2752 rr[iz] = B*std::sinh(ang);
2753 }
2754 zz[np1 - 1] = 0.;
2755 rr[np1 - 1] = 0.;
2756
2757 // 2nd polyline
2758 zz[np1] = H;
2759 rr[np1] = 0.;
2760 zz[np1 + 1] = 0.;
2761 rr[np1 + 1] = 0.;
2762
2763 // R O T A T E P O L Y L I N E S
2764
2765 G4double phi = 0.;
2766 G4double dphi = CLHEP::twopi;
2767 RotateAroundZ(0, phi, dphi, np1, np2, zz, rr, -1, -1);
2768 SetReferences();
2769
2770 delete [] zz;
2771 delete [] rr;
2772}
2773
2774HepPolyhedronHyperbolicMirror::~HepPolyhedronHyperbolicMirror() {}
2775
2776G4ThreadLocal G4int HepPolyhedron::fNumberOfRotationSteps = DEFAULT_NUMBER_OF_STEPS;
2777/***********************************************************************
2778 * *
2779 * Name: HepPolyhedron::fNumberOfRotationSteps Date: 24.06.97 *
2780 * Author: J.Allison (Manchester University) Revised: *
2781 * *
2782 * Function: Number of steps for whole circle *
2783 * *
2784 ***********************************************************************/
2785
2786#include "BooleanProcessor.src"
2787
2788HepPolyhedron HepPolyhedron::add(const HepPolyhedron & p) const
2789/***********************************************************************
2790 * *
2791 * Name: HepPolyhedron::add Date: 19.03.00 *
2792 * Author: E.Chernyaev Revised: *
2793 * *
2794 * Function: Boolean "union" of two polyhedra *
2795 * *
2796 ***********************************************************************/
2797{
2798 G4int ierr;
2799 BooleanProcessor processor;
2800 return processor.execute(OP_UNION, *this, p,ierr);
2801}
2802
2803HepPolyhedron HepPolyhedron::intersect(const HepPolyhedron & p) const
2804/***********************************************************************
2805 * *
2806 * Name: HepPolyhedron::intersect Date: 19.03.00 *
2807 * Author: E.Chernyaev Revised: *
2808 * *
2809 * Function: Boolean "intersection" of two polyhedra *
2810 * *
2811 ***********************************************************************/
2812{
2813 G4int ierr;
2814 BooleanProcessor processor;
2815 return processor.execute(OP_INTERSECTION, *this, p,ierr);
2816}
2817
2818HepPolyhedron HepPolyhedron::subtract(const HepPolyhedron & p) const
2819/***********************************************************************
2820 * *
2821 * Name: HepPolyhedron::add Date: 19.03.00 *
2822 * Author: E.Chernyaev Revised: *
2823 * *
2824 * Function: Boolean "subtraction" of "p" from "this" *
2825 * *
2826 ***********************************************************************/
2827{
2828 G4int ierr;
2829 BooleanProcessor processor;
2830 return processor.execute(OP_SUBTRACTION, *this, p,ierr);
2831}
2832
2833//NOTE : include the code of HepPolyhedronProcessor here
2834// since there is no BooleanProcessor.h
2835
2836#undef INTERSECTION
2837
2838#include "HepPolyhedronProcessor.src"
kCarTolerance
const G4double kCarTolerance
Definition: G4BoundingEnvelope.cc:37
B
G4double B(G4double temperature)
Definition: G4DNAElectronHoleRecombination.cc:68
FatalException
@ FatalException
Definition: G4ExceptionSeverity.hh:61
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:35
G4Normal3D
HepGeom::Normal3D< G4double > G4Normal3D
Definition: G4Normal3D.hh:34
G4Point3D
HepGeom::Point3D< G4double > G4Point3D
Definition: G4Point3D.hh:34
perMillion
static constexpr double perMillion
Definition: G4SIunits.hh:327
twopi
static constexpr double twopi
Definition: G4SIunits.hh:56
nm
static constexpr double nm
Definition: G4SIunits.hh:92
pi
static constexpr double pi
Definition: G4SIunits.hh:55
deg
static constexpr double deg
Definition: G4SIunits.hh:132
G4double
double G4double
Definition: G4Types.hh:83
G4bool
bool G4bool
Definition: G4Types.hh:86
G4int
int G4int
Definition: G4Types.hh:85
G4Vector3D
HepGeom::Vector3D< G4double > G4Vector3D
Definition: G4Vector3D.hh:34
A
const G4double A[17]
Definition: G4WaterStopping.cc:53
spatialTolerance
const G4double spatialTolerance
Definition: HepPolyhedron.cc:79
DEFAULT_NUMBER_OF_STEPS
#define DEFAULT_NUMBER_OF_STEPS
Definition: HepPolyhedron.h:195
G4Facet::edge
G4Edge edge[4]
Definition: HepPolyhedron.h:204
HepGeom::BasicVector3D::cross
BasicVector3D< T > cross(const BasicVector3D< T > &v) const
Definition: BasicVector3D.h:281
HepGeom::BasicVector3D::operator<<
std::ostream & operator<<(std::ostream &, const BasicVector3D< float > &)
Definition: BasicVector3D.cc:106
HepGeom::BasicVector3D::set
void set(T x1, T y1, T z1)
Definition: BasicVector3D.h:167
HepGeom::BasicVector3D::dot
T dot(const BasicVector3D< T > &v) const
Definition: BasicVector3D.h:275
HepPolyhedronBox::~HepPolyhedronBox
virtual ~HepPolyhedronBox()
Definition: HepPolyhedron.cc:1810
HepPolyhedronBox::HepPolyhedronBox
HepPolyhedronBox(G4double Dx, G4double Dy, G4double Dz)
Definition: HepPolyhedron.cc:1807
HepPolyhedronCone::~HepPolyhedronCone
virtual ~HepPolyhedronCone()
Definition: HepPolyhedron.cc:2141
HepPolyhedronCone::HepPolyhedronCone
HepPolyhedronCone(G4double Rmn1, G4double Rmx1, G4double Rmn2, G4double Rmx2, G4double Dz)
Definition: HepPolyhedron.cc:2136
HepPolyhedronCons::~HepPolyhedronCons
virtual ~HepPolyhedronCons()
Definition: HepPolyhedron.cc:2134
HepPolyhedronCons::HepPolyhedronCons
HepPolyhedronCons(G4double Rmn1, G4double Rmx1, G4double Rmn2, G4double Rmx2, G4double Dz, G4double Phi1, G4double Dphi)
Definition: HepPolyhedron.cc:2058
HepPolyhedronEllipsoid::HepPolyhedronEllipsoid
HepPolyhedronEllipsoid(G4double dx, G4double dy, G4double dz, G4double zcut1, G4double zcut2)
Definition: HepPolyhedron.cc:2523
HepPolyhedronEllipticalCone::HepPolyhedronEllipticalCone
HepPolyhedronEllipticalCone(G4double dx, G4double dy, G4double z, G4double zcut1)
Definition: HepPolyhedron.cc:2650
HepPolyhedronHype::HepPolyhedronHype
HepPolyhedronHype(G4double r1, G4double r2, G4double tan1, G4double tan2, G4double halfZ)
Definition: HepPolyhedron.cc:1978
HepPolyhedronHype::~HepPolyhedronHype
virtual ~HepPolyhedronHype()
Definition: HepPolyhedron.cc:2056
HepPolyhedronHyperbolicMirror::HepPolyhedronHyperbolicMirror
HepPolyhedronHyperbolicMirror(G4double a, G4double h, G4double r)
Definition: HepPolyhedron.cc:2714
HepPolyhedronPara::~HepPolyhedronPara
virtual ~HepPolyhedronPara()
Definition: HepPolyhedron.cc:1872
HepPolyhedronPara::HepPolyhedronPara
HepPolyhedronPara(G4double Dx, G4double Dy, G4double Dz, G4double Alpha, G4double Theta, G4double Phi)
Definition: HepPolyhedron.cc:1867
HepPolyhedronParaboloid::HepPolyhedronParaboloid
HepPolyhedronParaboloid(G4double r1, G4double r2, G4double dz, G4double Phi1, G4double Dphi)
Definition: HepPolyhedron.cc:1874
HepPolyhedronPcon::~HepPolyhedronPcon
virtual ~HepPolyhedronPcon()
Definition: HepPolyhedron.cc:2309
HepPolyhedronPcon::HepPolyhedronPcon
HepPolyhedronPcon(G4double phi, G4double dphi, G4int nz, const G4double *z, const G4double *rmin, const G4double *rmax)
Definition: HepPolyhedron.cc:2299
HepPolyhedronPgon::~HepPolyhedronPgon
virtual ~HepPolyhedronPgon()
Definition: HepPolyhedron.cc:2297
HepPolyhedronPgon::HepPolyhedronPgon
HepPolyhedronPgon(G4double phi, G4double dphi, G4int npdv, G4int nz, const G4double *z, const G4double *rmin, const G4double *rmax)
Definition: HepPolyhedron.cc:2156
HepPolyhedronSphere::~HepPolyhedronSphere
virtual ~HepPolyhedronSphere()
Definition: HepPolyhedron.cc:2406
HepPolyhedronSphere::HepPolyhedronSphere
HepPolyhedronSphere(G4double rmin, G4double rmax, G4double phi, G4double dphi, G4double the, G4double dthe)
Definition: HepPolyhedron.cc:2311
HepPolyhedronTet::~HepPolyhedronTet
virtual ~HepPolyhedronTet()
Definition: HepPolyhedron.cc:2521
HepPolyhedronTet::HepPolyhedronTet
HepPolyhedronTet(const G4double p0[3], const G4double p1[3], const G4double p2[3], const G4double p3[3])
Definition: HepPolyhedron.cc:2483
HepPolyhedronTorus::HepPolyhedronTorus
HepPolyhedronTorus(G4double rmin, G4double rmax, G4double rtor, G4double phi, G4double dphi)
Definition: HepPolyhedron.cc:2408
HepPolyhedronTorus::~HepPolyhedronTorus
virtual ~HepPolyhedronTorus()
Definition: HepPolyhedron.cc:2481
HepPolyhedronTrap::HepPolyhedronTrap
HepPolyhedronTrap(G4double Dz, G4double Theta, G4double Phi, G4double Dy1, G4double Dx1, G4double Dx2, G4double Alp1, G4double Dy2, G4double Dx3, G4double Dx4, G4double Alp2)
Definition: HepPolyhedron.cc:1812
HepPolyhedronTrap::~HepPolyhedronTrap
virtual ~HepPolyhedronTrap()
Definition: HepPolyhedron.cc:1865
HepPolyhedronTrd1::HepPolyhedronTrd1
HepPolyhedronTrd1(G4double Dx1, G4double Dx2, G4double Dy, G4double Dz)
Definition: HepPolyhedron.cc:1801
HepPolyhedronTrd1::~HepPolyhedronTrd1
virtual ~HepPolyhedronTrd1()
Definition: HepPolyhedron.cc:1805
HepPolyhedronTrd2::~HepPolyhedronTrd2
virtual ~HepPolyhedronTrd2()
Definition: HepPolyhedron.cc:1799
HepPolyhedronTrd2::HepPolyhedronTrd2
HepPolyhedronTrd2(G4double Dx1, G4double Dx2, G4double Dy1, G4double Dy2, G4double Dz)
Definition: HepPolyhedron.cc:1767
HepPolyhedronTube::~HepPolyhedronTube
virtual ~HepPolyhedronTube()
Definition: HepPolyhedron.cc:2154
HepPolyhedronTube::HepPolyhedronTube
HepPolyhedronTube(G4double Rmin, G4double Rmax, G4double Dz)
Definition: HepPolyhedron.cc:2150
HepPolyhedronTubs::HepPolyhedronTubs
HepPolyhedronTubs(G4double Rmin, G4double Rmax, G4double Dz, G4double Phi1, G4double Dphi)
Definition: HepPolyhedron.cc:2143
HepPolyhedronTubs::~HepPolyhedronTubs
virtual ~HepPolyhedronTubs()
Definition: HepPolyhedron.cc:2148
HepPolyhedron::SetNumberOfRotationSteps
static void SetNumberOfRotationSteps(G4int n)
Definition: HepPolyhedron.cc:261
HepPolyhedron::RotateAroundZ
void RotateAroundZ(G4int nstep, G4double phi, G4double dphi, G4int np1, G4int np2, const G4double *z, G4double *r, G4int nodeVis, G4int edgeVis)
Definition: HepPolyhedron.cc:475
HepPolyhedron::GetNextEdgeIndices
G4bool GetNextEdgeIndices(G4int &i1, G4int &i2, G4int &edgeFlag, G4int &iface1, G4int &iface2) const
Definition: HepPolyhedron.cc:1334
HepPolyhedron::GetUnitNormal
G4Normal3D GetUnitNormal(G4int iFace) const
Definition: HepPolyhedron.cc:1555
HepPolyhedron::GetFacet
void GetFacet(G4int iFace, G4int &n, G4int *iNodes, G4int *edgeFlags=0, G4int *iFaces=0) const
Definition: HepPolyhedron.cc:1444
HepPolyhedron::createTwistedTrap
G4int createTwistedTrap(G4double Dz, const G4double xy1[][2], const G4double xy2[][2])
Definition: HepPolyhedron.cc:1668
HepPolyhedron::pV
G4Point3D * pV
Definition: HepPolyhedron.h:219
HepPolyhedron::GetNextNormal
G4bool GetNextNormal(G4Normal3D &normal) const
Definition: HepPolyhedron.cc:1580
HepPolyhedron::operator=
HepPolyhedron & operator=(const HepPolyhedron &from)
Definition: HepPolyhedron.cc:146
HepPolyhedron::FindNodeNormal
G4Normal3D FindNodeNormal(G4int iFace, G4int iNode) const
Definition: HepPolyhedron.cc:220
HepPolyhedron::GetNumberOfRotationSteps
static G4int GetNumberOfRotationSteps()
Definition: HepPolyhedron.cc:248
HepPolyhedron::GetNormal
G4Normal3D GetNormal(G4int iFace) const
Definition: HepPolyhedron.cc:1530
HepPolyhedron::FindNeighbour
G4int FindNeighbour(G4int iFace, G4int iNode, G4int iOrder) const
Definition: HepPolyhedron.cc:192
HepPolyhedron::TriangulatePolygon
G4bool TriangulatePolygon(const std::vector< G4TwoVector > &polygon, std::vector< G4int > &result)
Definition: HepPolyhedron.cc:907
HepPolyhedron::RotateContourAroundZ
void RotateContourAroundZ(G4int nstep, G4double phi, G4double dphi, const std::vector< G4TwoVector > &rz, G4int nodeVis, G4int edgeVis)
Definition: HepPolyhedron.cc:750
HepPolyhedron::GetVertex
G4Point3D GetVertex(G4int index) const
Definition: HepPolyhedron.cc:1264
HepPolyhedron::GetSurfaceArea
G4double GetSurfaceArea() const
Definition: HepPolyhedron.cc:1617
HepPolyhedron::subtract
HepPolyhedron subtract(const HepPolyhedron &p) const
Definition: HepPolyhedron.cc:2818
HepPolyhedron::GetNextVertex
G4bool GetNextVertex(G4Point3D &vertex, G4int &edgeFlag) const
Definition: HepPolyhedron.cc:1284
HepPolyhedron::GetVolume
G4double GetVolume() const
Definition: HepPolyhedron.cc:1639
HepPolyhedron::intersect
HepPolyhedron intersect(const HepPolyhedron &p) const
Definition: HepPolyhedron.cc:2803
HepPolyhedron::Transform
HepPolyhedron & Transform(const G4Transform3D &t)
Definition: HepPolyhedron.cc:1212
HepPolyhedron::createPolyhedron
G4int createPolyhedron(G4int Nnodes, G4int Nfaces, const G4double xyz[][3], const G4int faces[][4])
Definition: HepPolyhedron.cc:1737
HepPolyhedron::AllocateMemory
void AllocateMemory(G4int Nvert, G4int Nface)
Definition: HepPolyhedron.cc:296
HepPolyhedron::GetNextUnitNormal
G4bool GetNextUnitNormal(G4Normal3D &normal) const
Definition: HepPolyhedron.cc:1601
HepPolyhedron::ResetNumberOfRotationSteps
static void ResetNumberOfRotationSteps()
Definition: HepPolyhedron.cc:283
HepPolyhedron::CheckSnip
G4bool CheckSnip(const std::vector< G4TwoVector > &contour, G4int a, G4int b, G4int c, G4int n, const G4int *V)
Definition: HepPolyhedron.cc:985
HepPolyhedron::RotateEdge
void RotateEdge(G4int k1, G4int k2, G4double r1, G4double r2, G4int v1, G4int v2, G4int vEdge, G4bool ifWholeCircle, G4int ns, G4int &kface)
Definition: HepPolyhedron.cc:342
HepPolyhedron::GetNextFacet
G4bool GetNextFacet(G4int &n, G4Point3D *nodes, G4int *edgeFlags=0, G4Normal3D *normals=0) const
Definition: HepPolyhedron.cc:1500
HepPolyhedron::fNumberOfRotationSteps
static G4ThreadLocal G4int fNumberOfRotationSteps
Definition: HepPolyhedron.h:217
HepPolyhedron::GetNextVertexIndex
G4bool GetNextVertexIndex(G4int &index, G4int &edgeFlag) const
Definition: HepPolyhedron.cc:1237
HepPolyhedron::pF
G4Facet * pF
Definition: HepPolyhedron.h:220
HepPolyhedron::add
HepPolyhedron add(const HepPolyhedron &p) const
Definition: HepPolyhedron.cc:2788
HepPolyhedron::GetNextEdge
G4bool GetNextEdge(G4Point3D &p1, G4Point3D &p2, G4int &edgeFlag) const
Definition: HepPolyhedron.cc:1403
HepPolyhedron::SetSideFacets
void SetSideFacets(G4int ii[4], G4int vv[4], G4int *kk, G4double *r, G4double dphi, G4int ns, G4int &kface)
Definition: HepPolyhedron.cc:404
CLHEP::deg
static constexpr double deg
Definition: SystemOfUnits.h:133
CLHEP::twopi
static constexpr double twopi
Definition: SystemOfUnits.h:56
CLHEP::perMillion
static constexpr double perMillion
Definition: SystemOfUnits.h:321
CLHEP::normal
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:79
CLHEP::pi
static constexpr double pi
Definition: SystemOfUnits.h:55
CLHEP::nm
static constexpr double nm
Definition: SystemOfUnits.h:93
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:112
G4INCL::Math::min
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
Definition: G4INCLGlobals.hh:117
geant4_check_module_cycles.nodes
nodes
Definition: geant4_check_module_cycles.py:60
G4ThreadLocal
#define G4ThreadLocal
Definition: tls.hh:77
ns
#define ns
Definition: xmlparse.cc:614
processor
#define processor
Definition: xmlparse.cc:617
Generated by doxygen 1.9.3