Bug Summary

File:scratch/gluex/scan-build-work/hdgeant4/src/G4fixes/HepPolyhedron.cc
Location:line 324, column 7
Description:Called C++ object pointer is null

Annotated Source Code

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