87 for (
G4int k=0; k<4; k++) {
88 ostr <<
" " << facet.
edge[k].
v <<
"/" << facet.
edge[k].
f;
95 ostr <<
"Nvertices=" << ph.
nvert <<
", Nfacets=" << ph.
nface << std::endl;
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()
102 for (i=1; i<=ph.
nface; i++) {
103 ostr <<
"face(" << i <<
")=" << ph.
pF[i] << std::endl;
115 : nvert(0), nface(0), pV(0), pF(0)
152 for (i=0; i<4; i++) {
153 if (iNode == std::abs(pF[iFace].edge[i].v))
break;
157 <<
"HepPolyhedron::FindNeighbour: face " << iFace
158 <<
" has no node " << iNode
164 if (pF[iFace].edge[i].v == 0) i = 2;
166 return (pF[iFace].edge[i].v > 0) ? 0 : pF[iFace].edge[i].f;
180 G4int k = iFace, iOrder = 1,
n = 1;
183 k = FindNeighbour(k, iNode, iOrder);
184 if (k == iFace)
break;
187 normal += GetUnitNormal(k);
189 if (iOrder < 0)
break;
194 return normal.
unit();
220 const G4int nMin = 3;
223 <<
"HepPolyhedron::SetNumberOfRotationSteps: attempt to set the\n"
224 <<
"number of steps per circle < " << nMin <<
"; forced to " << nMin
259 if (
pV != 0)
delete []
pV;
260 if (
pF != 0)
delete []
pF;
261 if (Nvert > 0 && Nface > 0) {
281 enum {DUMMY, BOTTOM, LEFT, BACK, RIGHT, FRONT, TOP};
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);
313 if (r1 == 0. && r2 == 0)
return;
318 G4int ii1 = ifWholeCircle ? i1 : i1+nds;
319 G4int ii2 = ifWholeCircle ? i2 : i2+nds;
320 G4int vv = ifWholeCircle ? vEdge : 1;
324 pF[kface++] =
G4Facet(i1,0, v2*i2,0, (i2+1),0);
325 }
else if (r2 == 0.) {
326 pF[kface++] =
G4Facet(i1,0, i2,0, v1*(i1+1),0);
328 pF[kface++] =
G4Facet(i1,0, v2*i2,0, (i2+1),0, v1*(i1+1),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);
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);
342 pF[kface++] =
G4Facet(vEdge*i1,0, vv*i2,0, v1*ii1,0);
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);
348 pF[kface++] =
G4Facet(vEdge*i1,0, v2*i2,0, vv*ii2,0, v1*ii1,0);
373 G4int k1, k2, k3, k4;
376 for (
G4int i=0; i<4; i++) {
378 k2 = (i == 3) ? ii[0] : ii[i+1];
379 if (r[k1] == 0. && r[k2] == 0.) vv[i] = -1;
383 if (ii[1] == ii[2]) {
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]]) {
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]]) {
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);
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);
454 G4bool ifWholeCircle = (std::abs(dphi-wholeCircle) <
perMillion) ?
true :
false;
455 G4double delPhi = ifWholeCircle ? wholeCircle : dphi;
456 G4int nSphi = (nstep > 0) ?
458 if (nSphi == 0) nSphi = 1;
459 G4int nVphi = ifWholeCircle ? nSphi : nSphi+1;
460 G4bool ifClosed = np1 > 0 ?
false :
true;
464 G4int absNp1 = std::abs(np1);
465 G4int absNp2 = std::abs(np2);
467 G4int i1end = absNp1-1;
468 G4int i2beg = absNp1;
469 G4int i2end = absNp1+absNp2-1;
472 for(i=i1beg; i<=i2end; i++) {
477 for (i=i1beg; i<=i1end; i++) {
478 j += (r[i] == 0.) ? 1 : nVphi;
484 if (r[i2beg] != r[i1beg] || z[i2beg] != z[i1beg]) {
485 j += (r[i2beg] == 0.) ? 1 : nVphi;
489 for(i=i2beg+1; i<i2end; i++) {
490 j += (r[i] == 0.) ? 1 : nVphi;
493 if (r[i2end] != r[i1end] || z[i2end] != z[i1end]) {
494 if (absNp2 > 1) j += (r[i2end] == 0.) ? 1 : nVphi;
500 k = ifClosed ? absNp1*nSphi : (absNp1-1)*nSphi;
503 for(i=i2beg; i<i2end; i++) {
504 if (r[i] > 0. || r[i+1] > 0.) k += nSphi;
508 if (r[i2end] > 0. || r[i2beg] > 0.) k += nSphi;
513 if (ifSide1 && (r[i1beg] > 0. || r[i2beg] > 0.)) k += nSphi;
514 if (ifSide2 && (r[i1end] > 0. || r[i2end] > 0.)) k += nSphi;
517 if (!ifWholeCircle) {
518 k += ifClosed ? 2*absNp1 : 2*(absNp1-1);
528 kk =
new G4int[absNp1+absNp2];
531 for(i=i1beg; i<=i1end; i++) {
534 {
pV[k++] =
G4Point3D(0, 0, z[i]); }
else { k += nVphi; }
541 {
pV[k++] =
G4Point3D(0, 0, z[i]); }
else { k += nVphi; }
546 for(i=i2beg+1; i<i2end; i++) {
549 {
pV[k++] =
G4Point3D(0, 0, z[i]); }
else { k += nVphi; }
564 for(j=0; j<nVphi; j++) {
565 cosPhi = std::cos(phi+j*delPhi/nSphi);
566 sinPhi = std::sin(phi+j*delPhi/nSphi);
567 for(i=i1beg; i<=i2end; i++) {
569 pV[kk[i]+j] =
G4Point3D(r[i]*cosPhi,r[i]*sinPhi,z[i]);
578 v2 = ifClosed ? nodeVis : 1;
579 for(i=i1beg; i<i1end; i++) {
581 if (!ifClosed && i == i1end-1) {
584 v2 = (r[i] == r[i+1] && r[i+1] == r[i+2]) ? -1 : nodeVis;
586 RotateEdge(kk[i], kk[i+1], r[i], r[i+1], v1, v2,
587 edgeVis, ifWholeCircle, nSphi, k);
590 RotateEdge(kk[i1end], kk[i1beg], r[i1end],r[i1beg], nodeVis, nodeVis,
591 edgeVis, ifWholeCircle, nSphi, k);
597 v2 = ifClosed ? nodeVis : 1;
598 for(i=i2beg; i<i2end; i++) {
600 if (!ifClosed && i==i2end-1) {
603 v2 = (r[i] == r[i+1] && r[i+1] == r[i+2]) ? -1 : nodeVis;
605 RotateEdge(kk[i+1], kk[i], r[i+1], r[i], v2, v1,
606 edgeVis, ifWholeCircle, nSphi, k);
609 RotateEdge(kk[i2beg], kk[i2end], r[i2beg], r[i2end], nodeVis, nodeVis,
610 edgeVis, ifWholeCircle, nSphi, k);
618 RotateEdge(kk[i2beg], kk[i1beg], r[i2beg], r[i1beg], 1, 1,
619 -1, ifWholeCircle, nSphi, k);
622 RotateEdge(kk[i1end], kk[i2end], r[i1end], r[i2end], 1, 1,
623 -1, ifWholeCircle, nSphi, k);
629 if (!ifWholeCircle) {
634 for (i=i1beg; i<=i1end; 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;
646 for (i=i1beg; i<i1end; i++) {
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;
653 vv[2] = (i == i1end-1) ? 1 : -1;
664 <<
"Polyhedron::RotateAroundZ: number of generated faces ("
665 << k-1 <<
") is not equal to the number of allocated faces ("
681 if (
nface <= 0)
return;
683 struct edgeListMember {
684 edgeListMember *next;
688 } *edgeList, *freeList, **headList;
693 edgeList =
new edgeListMember[2*
nface];
694 headList =
new edgeListMember*[
nvert];
697 for (i=0; i<
nvert; i++) {
701 for (i=0; i<2*
nface-1; i++) {
702 edgeList[i].next = &edgeList[i+1];
704 edgeList[2*nface-1].next = 0;
708 G4int iface, iedge, nedge, i1, i2, k1, k2;
709 edgeListMember *prev, *cur;
711 for(iface=1; iface<=
nface; iface++) {
712 nedge = (
pF[iface].
edge[3].
v == 0) ? 3 : 4;
713 for (iedge=0; iedge<nedge; 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;
719 k2 = (i1 > i2) ? i1 : i2;
724 headList[k1] = freeList;
727 <<
"Polyhedron::SetReferences: bad link "
731 freeList = freeList->next;
741 headList[k1] = cur->next;
742 cur->next = freeList;
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;
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
765 prev->next = freeList;
768 <<
"Polyhedron::SetReferences: bad link "
772 freeList = freeList->next;
782 prev->next = cur->next;
783 cur->next = freeList;
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;
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
806 for (i=0; i<
nvert; i++) {
807 if (headList[i] != 0) {
809 <<
"Polyhedron::SetReferences: List " << i <<
" is not empty"
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];
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];
883 G4int vIndex = pF[iFace].edge[iQVertex].v;
885 edgeFlag = (vIndex > 0) ? 1 : 0;
886 index = std::abs(vIndex);
888 if (iQVertex >= 3 || pF[iFace].edge[iQVertex+1].v == 0) {
890 if (++iFace > nface) iFace = 1;
908 if (index <= 0 || index > nvert) {
910 <<
"HepPolyhedron::GetVertex: irrelevant index " << index
931 G4bool rep = GetNextVertexIndex(index, edgeFlag);
952 if (nface == 0)
return false;
954 G4int k = pF[iFace].edge[iNode].v;
955 if (k > 0) { edgeFlag = 1; }
else { edgeFlag = -1; k = -k; }
957 normal = FindNodeNormal(iFace,k);
958 if (iNode >= 3 || pF[iFace].edge[iNode+1].v == 0) {
960 if (++iFace > nface) iFace = 1;
984 G4int k1, k2, kflag, kface1, kface2;
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;
994 k1 = pF[iFace].edge[iQVertex].v;
998 kface2 = pF[iFace].edge[iQVertex].f;
999 if (iQVertex >= 3 || pF[iFace].edge[iQVertex+1].v == 0) {
1001 k2 = std::abs(pF[iFace].edge[iQVertex].v);
1005 k2 = std::abs(pF[iFace].edge[iQVertex].v);
1007 }
while (iOrder*k1 > iOrder*k2);
1009 i1 = k1; i2 = k2; edgeFlag = (kflag > 0) ? 1 : 0;
1010 iface1 = kface1; iface2 = kface2;
1012 if (iFace > nface) {
1013 iFace = 1; iOrder = 1;
1032 G4int kface1, kface2;
1033 return GetNextEdgeIndices(i1, i2, edgeFlag, kface1, kface2);
1039 G4int &edgeFlag)
const
1051 G4bool rep = GetNextEdgeIndices(i1,i2,edgeFlag);
1072 G4bool rep = GetNextEdgeIndices(i1,i2,edgeFlag,iface1,iface2);
1089 if (iFace < 1 || iFace > nface) {
1091 <<
"HepPolyhedron::GetFacet: irrelevant index " << iFace
1096 for (i=0; i<4; i++) {
1097 k = pF[iFace].edge[i].v;
1099 if (iFaces != 0) iFaces[i] = pF[iFace].edge[i].f;
1102 if (edgeFlags != 0) edgeFlags[i] = 1;
1105 if (edgeFlags != 0) edgeFlags[i] = -1;
1124 GetFacet(index,
n, iNodes, edgeFlags);
1126 for (
G4int i=0; i<
n; i++) {
1127 nodes[i] = pV[iNodes[i]];
1128 if (normals != 0) normals[i] = FindNodeNormal(index,iNodes[i]);
1148 if (edgeFlags == 0) {
1149 GetFacet(iFace,
n, nodes);
1150 }
else if (normals == 0) {
1151 GetFacet(iFace,
n, nodes, edgeFlags);
1153 GetFacet(iFace,
n, nodes, edgeFlags, normals);
1156 if (++iFace > nface) {
1174 if (iFace < 1 || iFace > nface) {
1176 <<
"HepPolyhedron::GetNormal: irrelevant index " << iFace
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]);
1199 if (iFace < 1 || iFace > nface) {
1201 <<
"HepPolyhedron::GetUnitNormal: irrelevant index " << iFace
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();
1226 normal = GetNormal(iFace);
1227 if (++iFace > nface) {
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();
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);
1292 pt = (
pV[i0]+
pV[i1]+
pV[i2]) * (1./3.);
1294 pt = (
pV[i0]+
pV[i1]+
pV[i2]+
pV[i3]) * 0.25;
1296 v += ((
pV[i2] -
pV[i0]).cross(
pV[i3] -
pV[i1])).dot(pt);
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,
1343 pF[ 1]=
G4Facet(1,LEFT_BOTTOM, 4,BACK_BOTTOM, 3,RIGHT_BOTTOM, 2,FRONT_BOTTOM);
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);
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);
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);
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);
1365 pF[18]=
G4Facet(5,FRONT_TOP, 6,RIGHT_TOP, 7,BACK_TOP, 8,LEFT_TOP);
1373 const G4int faces[][4])
1389 if (
nvert == 0)
return 1;
1391 for (
G4int i=0; i<Nnodes; i++) {
1392 pV[i+1] =
G4Point3D(xyz[i][0], xyz[i][1], xyz[i][2]);
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);
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);
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);
1504 :
HepPolyhedronTrap(Dz, Theta, Phi, Dy, Dx, Dx, Alpha, Dy, Dx, Dx, Alpha) {}
1533 if (r1 < 0. || r2 <= 0.) k = 1;
1535 if (dz <= 0.) k += 2;
1541 phi2 = sPhi; phi1 = phi2 + dPhi;
1545 phi1 = sPhi; phi2 = phi1 + wholeCircle;
1549 phi1 = sPhi; phi2 = phi1 + dPhi;
1553 if (std::abs(dphi-wholeCircle) <
perMillion) dphi = wholeCircle;
1554 if (dphi > wholeCircle) k += 4;
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
1573 G4double k1 = (r2*r2 - r1*r1) / 2 / dz;
1581 for(
G4int i = 1; i < n - 1; i++)
1583 rr[i] = rr[i-1] - dl;
1584 zz[i] = (rr[i]*rr[i] - k2) / k1;
1637 if (r2 < 0. || r1 < 0. ) k = 1;
1638 if (r1 > r2 ) k = 1;
1639 if (r1 == r2) k = 1;
1641 if (halfZ <= 0.) k += 2;
1643 if (sqrtan1<0.||sqrtan2<0.) k += 4;
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
1669 rr[0] = std::sqrt(sqrtan2*halfZ*halfZ+k2);
1671 for(
G4int i = 1; i < n-1; i++)
1673 zz[i] = zz[i-1] - dz;
1674 rr[i] =std::sqrt(sqrtan2*zz[i]*zz[i]+k2);
1681 rr[
n] = std::sqrt(sqrtan1*halfZ*halfZ+k1);
1683 for(
G4int i = n+1; i < n+
n; i++)
1685 zz[i] = zz[i-1] - dz;
1686 rr[i] =std::sqrt(sqrtan1*zz[i]*zz[i]+k1);
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;
1733 if (Dz <= 0.) k += 2;
1737 phi2 = Phi1; phi1 = phi2 - Dphi;
1738 }
else if (Dphi == 0.) {
1739 phi1 = Phi1; phi2 = phi1 + wholeCircle;
1741 phi1 = Phi1; phi2 = phi1 + Dphi;
1744 if (std::abs(dphi-wholeCircle) <
perMillion) dphi = wholeCircle;
1745 if (dphi > wholeCircle) k += 4;
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
1826 if (dphi <= 0. || dphi >
twopi) {
1828 <<
"HepPolyhedronPgon/Pcon: wrong delta phi = " << dphi
1835 <<
"HepPolyhedronPgon/Pcon: number of z-planes less than two = " << nz
1842 <<
"HepPolyhedronPgon/Pcon: error in number of phi-steps =" << npdv
1848 for (i=0; i<nz; i++) {
1849 if (rmin[i] < 0. || rmax[i] < 0. || rmin[i] > rmax[i]) {
1851 <<
"HepPolyhedronPgon: error in radiuses rmin[" << i <<
"]="
1852 << rmin[i] <<
" rmax[" << i <<
"]=" << rmax[i]
1864 if (z[0] > z[nz-1]) {
1865 for (i=0; i<nz; i++) {
1872 for (i=0; i<nz; i++) {
1874 rr[i] = rmax[nz-i-1];
1875 zz[i+nz] = z[nz-i-1];
1876 rr[i+nz] = rmin[nz-i-1];
1882 RotateAroundZ(npdv, phi, dphi, nz, nz, zz, rr, -1, (npdv == 0) ? -1 : 1);
1920 if (dphi <= 0. || dphi >
twopi) {
1922 <<
"HepPolyhedronSphere: wrong delta phi = " << dphi
1927 if (the < 0. || the >
pi) {
1929 <<
"HepPolyhedronSphere: wrong theta = " << the
1934 if (dthe <= 0. || dthe > pi) {
1936 <<
"HepPolyhedronSphere: wrong delta theta = " << dthe
1941 if (the+dthe > pi) {
1943 <<
"HepPolyhedronSphere: wrong theta + delta theta = "
1944 << the <<
" " << dthe
1949 if (rmin < 0. || rmin >= rmax) {
1951 <<
"HepPolyhedronSphere: error in radiuses"
1952 <<
" rmin=" << rmin <<
" rmax=" << rmax
1961 if (np1 <= 1) np1 = 2;
1970 for (
G4int i=0; i<np1; i++) {
1971 cosa = std::cos(the+i*a);
1972 sina = std::sin(the+i*a);
1976 zz[i+np1] = rmin*cosa;
1977 rr[i+np1] = rmin*sina;
2018 if (dphi <= 0. || dphi >
twopi) {
2020 <<
"HepPolyhedronTorus: wrong delta phi = " << dphi
2025 if (rmin < 0. || rmin >= rmax || rmax >= rtor) {
2027 <<
"HepPolyhedronTorus: error in radiuses"
2028 <<
" rmin=" << rmin <<
" rmax=" << rmax <<
" rtorus=" << rtor
2045 for (
G4int i=0; i<np1; i++) {
2046 cosa = std::cos(i*a);
2047 sina = std::sin(i*a);
2049 rr[i] = rtor+rmax*sina;
2051 zz[i+np1] = rmin*cosa;
2052 rr[i+np1] = rtor+rmin*sina;
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;
2099 std::cerr <<
"HepPolyhedronEllipsoid: bad z semi-axis: cz = " << cz
2114 sthe= std::acos(zCut2/cz);
2123 dthe= std::acos(zCut1/cz)-sthe;
2138 G4Exception(
"HepPolyhedronEllipsoid::HepPolyhedronEllipsoid",
2151 for (
G4int i=0; i<np1-cutflag; i++) {
2152 cosa = std::cos(sthe+i*a);
2153 sina = std::sin(sthe+i*a);
2166 std::cerr <<
"Logic error in HepPolyhedronEllipsoid, memory corrupted!"
2171 std::cerr <<
"Warning: logic error in HepPolyhedronEllipsoid."
2191 p->
setX( p->
x() * ax/cz );
2192 p->
setY( p->
y() * by/cz );
2219 if ( (ax <= 0.) || (ay <= 0.) || (h <= 0.) || (zTopCut <= 0.) ) { k = 1; }
2222 std::cerr <<
"HepPolyhedronCone: error in input parameters";
2223 std::cerr << std::endl;
2229 zTopCut = (h >= zTopCut ? zTopCut : h);
2238 rr[0] = (h-zTopCut);
2239 rr[1] = (h+zTopCut);
2255 p->
setX( p->
x() * ax );
2256 p->
setY( p->
y() * ay );
2273 #include "BooleanProcessor.src"
2287 return processor.execute(OP_UNION, *
this,
p,ierr);
2302 return processor.execute(OP_INTERSECTION, *
this,
p,ierr);
2317 return processor.execute(OP_SUBTRACTION, *
this,
p,ierr);
2325 #include "HepPolyhedronProcessor.src"
HepPolyhedronTubs(G4double Rmin, G4double Rmax, G4double Dz, G4double Phi1, G4double Dphi)
G4Point3D GetVertex(G4int index) const
const G4double spatialTolerance
HepPolyhedronTorus(G4double rmin, G4double rmax, G4double rtor, G4double phi, G4double dphi)
G4int createPolyhedron(G4int Nnodes, G4int Nfaces, const G4double xyz[][3], const G4int faces[][4])
HepPolyhedronPcon(G4double phi, G4double dphi, G4int nz, const G4double *z, const G4double *rmin, const G4double *rmax)
virtual ~HepPolyhedronBox()
G4bool GetNextFacet(G4int &n, G4Point3D *nodes, G4int *edgeFlags=0, G4Normal3D *normals=0) const
std::vector< ExP01TrackerHit * > a
G4Normal3D FindNodeNormal(G4int iFace, G4int iNode) const
std::ostream & operator<<(std::ostream &, const BasicVector3D< float > &)
HepPolyhedron subtract(const HepPolyhedron &p) const
G4double GetSurfaceArea() const
static constexpr double nm
virtual ~HepPolyhedronCone()
HepPolyhedronCone(G4double Rmn1, G4double Rmx1, G4double Rmn2, G4double Rmx2, G4double Dz)
static void ResetNumberOfRotationSteps()
HepPolyhedronTrap(G4double Dz, G4double Theta, G4double Phi, G4double Dy1, G4double Dx1, G4double Dx2, G4double Alp1, G4double Dy2, G4double Dx3, G4double Dx4, G4double Alp2)
BasicVector3D< T > unit() const
HepPolyhedronEllipticalCone(G4double dx, G4double dy, G4double z, G4double zcut1)
void RotateAroundZ(G4int nstep, G4double phi, G4double dphi, G4int np1, G4int np2, const G4double *z, G4double *r, G4int nodeVis, G4int edgeVis)
virtual ~HepPolyhedronTrd1()
BasicVector3D< T > cross(const BasicVector3D< T > &v) const
HepPolyhedronTube(G4double Rmin, G4double Rmax, G4double Dz)
HepPolyhedronBox(G4double Dx, G4double Dy, G4double Dz)
static G4ThreadLocal G4int fNumberOfRotationSteps
static constexpr double perMillion
static double normal(HepRandomEngine *eptr)
G4int FindNeighbour(G4int iFace, G4int iNode, G4int iOrder) const
HepPolyhedronEllipsoid(G4double dx, G4double dy, G4double dz, G4double zcut1, G4double zcut2)
virtual ~HepPolyhedronTrd2()
static void SetNumberOfRotationSteps(G4int n)
virtual ~HepPolyhedronPgon()
virtual ~HepPolyhedronPcon()
G4bool GetNextVertexIndex(G4int &index, G4int &edgeFlag) const
HepPolyhedron add(const HepPolyhedron &p) const
G4bool GetNextEdge(G4Point3D &p1, G4Point3D &p2, G4int &edgeFlag) const
void GetFacet(G4int iFace, G4int &n, G4int *iNodes, G4int *edgeFlags=0, G4int *iFaces=0) const
static constexpr double deg
HepPolyhedronPgon(G4double phi, G4double dphi, G4int npdv, G4int nz, const G4double *z, const G4double *rmin, const G4double *rmax)
HepPolyhedron & Transform(const G4Transform3D &t)
HepGeom::Normal3D< G4double > G4Normal3D
static constexpr double twopi
void RotateEdge(G4int k1, G4int k2, G4double r1, G4double r2, G4int v1, G4int v2, G4int vEdge, G4bool ifWholeCircle, G4int ns, G4int &kface)
G4bool GetNextUnitNormal(G4Normal3D &normal) const
virtual ~HepPolyhedronTubs()
virtual ~HepPolyhedronTube()
G4bool GetNextVertex(G4Point3D &vertex, G4int &edgeFlag) const
void SetSideFacets(G4int ii[4], G4int vv[4], G4int *kk, G4double *r, G4double dphi, G4int ns, G4int &kface)
static constexpr double perMillion
virtual ~HepPolyhedronPara()
const G4ThreeVector const G4double const
virtual ~HepPolyhedronTorus()
virtual ~HepPolyhedronParaboloid()
HepGeom::Vector3D< G4double > G4Vector3D
static G4int GetNumberOfRotationSteps()
virtual ~HepPolyhedronEllipticalCone()
static constexpr double twopi
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
G4bool GetNextNormal(G4Normal3D &normal) const
G4Normal3D GetUnitNormal(G4int iFace) const
static constexpr double nm
G4int createTwistedTrap(G4double Dz, const G4double xy1[][2], const G4double xy2[][2])
G4Normal3D GetNormal(G4int iFace) const
virtual ~HepPolyhedronCons()
HepPolyhedronSphere(G4double rmin, G4double rmax, G4double phi, G4double dphi, G4double the, G4double dthe)
virtual ~HepPolyhedronTrap()
static constexpr double deg
HepPolyhedronParaboloid(G4double r1, G4double r2, G4double dz, G4double Phi1, G4double Dphi)
virtual ~HepPolyhedronSphere()
G4double GetVolume() const
static constexpr double pi
virtual ~HepPolyhedronEllipsoid()
void AllocateMemory(G4int Nvert, G4int Nface)
HepPolyhedronTrd1(G4double Dx1, G4double Dx2, G4double Dy, G4double Dz)
G4bool GetNextEdgeIndices(G4int &i1, G4int &i2, G4int &edgeFlag, G4int &iface1, G4int &iface2) const
HepPolyhedron intersect(const HepPolyhedron &p) const
HepGeom::Point3D< G4double > G4Point3D
virtual ~HepPolyhedronHype()
HepPolyhedronPara(G4double Dx, G4double Dy, G4double Dz, G4double Alpha, G4double Theta, G4double Phi)
static constexpr double pi
HepPolyhedron & operator=(const HepPolyhedron &from)
#define DEFAULT_NUMBER_OF_STEPS
HepPolyhedronCons(G4double Rmn1, G4double Rmx1, G4double Rmn2, G4double Rmx2, G4double Dz, G4double Phi1, G4double Dphi)
HepPolyhedronTrd2(G4double Dx1, G4double Dx2, G4double Dy1, G4double Dy2, G4double Dz)
HepPolyhedronHype(G4double r1, G4double r2, G4double tan1, G4double tan2, G4double halfZ)