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