70 using CLHEP::perMillion;
86 for (
G4int k=0; k<4; k++) {
87 ostr <<
" " << facet.edge[k].v <<
"/" << facet.edge[k].f;
94 ostr <<
"Nverteces=" << ph.
nvert <<
", Nfacets=" << ph.
nface << std::endl;
96 for (i=1; i<=ph.
nvert; i++) {
97 ostr <<
"xyz(" << i <<
")="
98 << ph.
pV[i].
x() <<
' ' << ph.
pV[i].
y() <<
' ' << ph.
pV[i].
z()
101 for (i=1; i<=ph.
nface; i++) {
102 ostr <<
"face(" << i <<
")=" << ph.
pF[i] << std::endl;
114 : nvert(0), nface(0), pV(0), pF(0)
151 for (i=0; i<4; i++) {
152 if (iNode == std::abs(pF[iFace].edge[i].
v))
break;
156 <<
"HepPolyhedron::FindNeighbour: face " << iFace
157 <<
" has no node " << iNode
163 if (pF[iFace].edge[i].
v == 0) i = 2;
165 return (pF[iFace].edge[i].
v > 0) ? 0 : pF[iFace].edge[i].f;
179 G4int k = iFace, iOrder = 1,
n = 1;
182 k = FindNeighbour(k, iNode, iOrder);
183 if (k == iFace)
break;
186 normal += GetUnitNormal(k);
188 if (iOrder < 0)
break;
193 return normal.
unit();
219 const G4int nMin = 3;
222 <<
"HepPolyhedron::SetNumberOfRotationSteps: attempt to set the\n"
223 <<
"number of steps per circle < " << nMin <<
"; forced to " << nMin
258 if (
pV != 0)
delete []
pV;
259 if (
pF != 0)
delete []
pF;
260 if (Nvert > 0 && Nface > 0) {
280 enum {DUMMY, BOTTOM, LEFT, BACK, RIGHT, FRONT, TOP};
282 pF[1] =
G4Facet(1,LEFT, 4,BACK, 3,RIGHT, 2,FRONT);
283 pF[2] =
G4Facet(5,TOP, 8,BACK, 4,BOTTOM, 1,FRONT);
284 pF[3] =
G4Facet(8,TOP, 7,RIGHT, 3,BOTTOM, 4,LEFT);
285 pF[4] =
G4Facet(7,TOP, 6,FRONT, 2,BOTTOM, 3,BACK);
286 pF[5] =
G4Facet(6,TOP, 5,LEFT, 1,BOTTOM, 2,RIGHT);
287 pF[6] =
G4Facet(5,FRONT, 6,RIGHT, 7,BACK, 8,LEFT);
312 if (r1 == 0. && r2 == 0)
return;
317 G4int ii1 = ifWholeCircle ? i1 : i1+nds;
318 G4int ii2 = ifWholeCircle ? i2 : i2+nds;
319 G4int vv = ifWholeCircle ? vEdge : 1;
323 pF[kface++] =
G4Facet(i1,0, v2*i2,0, (i2+1),0);
324 }
else if (r2 == 0.) {
325 pF[kface++] =
G4Facet(i1,0, i2,0, v1*(i1+1),0);
327 pF[kface++] =
G4Facet(i1,0, v2*i2,0, (i2+1),0, v1*(i1+1),0);
331 pF[kface++] =
G4Facet(vv*i1,0, v2*i2,0, vEdge*(i2+1),0);
332 for (i2++,i=1; i<nds-1; i2++,i++) {
333 pF[kface++] =
G4Facet(vEdge*i1,0, v2*i2,0, vEdge*(i2+1),0);
335 pF[kface++] =
G4Facet(vEdge*i1,0, v2*i2,0, vv*ii2,0);
336 }
else if (r2 == 0.) {
337 pF[kface++] =
G4Facet(vv*i1,0, vEdge*i2,0, v1*(i1+1),0);
338 for (i1++,i=1; i<nds-1; i1++,i++) {
339 pF[kface++] =
G4Facet(vEdge*i1,0, vEdge*i2,0, v1*(i1+1),0);
341 pF[kface++] =
G4Facet(vEdge*i1,0, vv*i2,0, v1*ii1,0);
343 pF[kface++] =
G4Facet(vv*i1,0, v2*i2,0, vEdge*(i2+1),0,v1*(i1+1),0);
344 for (i1++,i2++,i=1; i<nds-1; i1++,i2++,i++) {
345 pF[kface++] =
G4Facet(vEdge*i1,0, v2*i2,0, vEdge*(i2+1),0,v1*(i1+1),0);
347 pF[kface++] =
G4Facet(vEdge*i1,0, v2*i2,0, vv*ii2,0, v1*ii1,0);
372 G4int k1, k2, k3, k4;
375 for (
G4int i=0; i<4; i++) {
377 k2 = (i == 3) ? ii[0] : ii[i+1];
378 if (r[k1] == 0. && r[k2] == 0.) vv[i] = -1;
382 if (ii[1] == ii[2]) {
386 pF[kface++] =
G4Facet(vv[0]*k1,0, vv[2]*k2,0, vv[3]*k3,0);
387 if (r[ii[0]] != 0.) k1 += nds;
388 if (r[ii[2]] != 0.) k2 += nds;
389 if (r[ii[3]] != 0.) k3 += nds;
390 pF[kface++] =
G4Facet(vv[2]*k3,0, vv[0]*k2,0, vv[3]*k1,0);
391 }
else if (kk[ii[0]] == kk[ii[1]]) {
395 pF[kface++] =
G4Facet(vv[1]*k1,0, vv[2]*k2,0, vv[3]*k3,0);
396 if (r[ii[0]] != 0.) k1 += nds;
397 if (r[ii[2]] != 0.) k2 += nds;
398 if (r[ii[3]] != 0.) k3 += nds;
399 pF[kface++] =
G4Facet(vv[2]*k3,0, vv[1]*k2,0, vv[3]*k1,0);
400 }
else if (kk[ii[2]] == kk[ii[3]]) {
404 pF[kface++] =
G4Facet(vv[0]*k1,0, vv[1]*k2,0, vv[3]*k3,0);
405 if (r[ii[0]] != 0.) k1 += nds;
406 if (r[ii[1]] != 0.) k2 += nds;
407 if (r[ii[2]] != 0.) k3 += nds;
408 pF[kface++] =
G4Facet(vv[1]*k3,0, vv[0]*k2,0, vv[3]*k1,0);
414 pF[kface++] =
G4Facet(vv[0]*k1,0, vv[1]*k2,0, vv[2]*k3,0, vv[3]*k4,0);
415 if (r[ii[0]] != 0.) k1 += nds;
416 if (r[ii[1]] != 0.) k2 += nds;
417 if (r[ii[2]] != 0.) k3 += nds;
418 if (r[ii[3]] != 0.) k4 += nds;
419 pF[kface++] =
G4Facet(vv[2]*k4,0, vv[1]*k3,0, vv[0]*k2,0, vv[3]*k1,0);
453 G4bool ifWholeCircle = (std::abs(dphi-wholeCircle) <
perMillion) ?
true :
false;
454 G4double delPhi = ifWholeCircle ? wholeCircle : dphi;
455 G4int nSphi = (nstep > 0) ?
457 if (nSphi == 0) nSphi = 1;
458 G4int nVphi = ifWholeCircle ? nSphi : nSphi+1;
459 G4bool ifClosed = np1 > 0 ?
false :
true;
463 G4int absNp1 = std::abs(np1);
464 G4int absNp2 = std::abs(np2);
466 G4int i1end = absNp1-1;
467 G4int i2beg = absNp1;
468 G4int i2end = absNp1+absNp2-1;
471 for(i=i1beg; i<=i2end; i++) {
476 for (i=i1beg; i<=i1end; i++) {
477 j += (r[i] == 0.) ? 1 : nVphi;
483 if (r[i2beg] != r[i1beg] || z[i2beg] != z[i1beg]) {
484 j += (r[i2beg] == 0.) ? 1 : nVphi;
488 for(i=i2beg+1; i<i2end; i++) {
489 j += (r[i] == 0.) ? 1 : nVphi;
492 if (r[i2end] != r[i1end] || z[i2end] != z[i1end]) {
493 if (absNp2 > 1) j += (r[i2end] == 0.) ? 1 : nVphi;
499 k = ifClosed ? absNp1*nSphi : (absNp1-1)*nSphi;
502 for(i=i2beg; i<i2end; i++) {
503 if (r[i] > 0. || r[i+1] > 0.) k += nSphi;
507 if (r[i2end] > 0. || r[i2beg] > 0.) k += nSphi;
512 if (ifSide1 && (r[i1beg] > 0. || r[i2beg] > 0.)) k += nSphi;
513 if (ifSide2 && (r[i1end] > 0. || r[i2end] > 0.)) k += nSphi;
516 if (!ifWholeCircle) {
517 k += ifClosed ? 2*absNp1 : 2*(absNp1-1);
527 kk =
new G4int[absNp1+absNp2];
530 for(i=i1beg; i<=i1end; i++) {
533 {
pV[k++] =
G4Point3D(0, 0, z[i]); }
else { k += nVphi; }
540 {
pV[k++] =
G4Point3D(0, 0, z[i]); }
else { k += nVphi; }
545 for(i=i2beg+1; i<i2end; i++) {
548 {
pV[k++] =
G4Point3D(0, 0, z[i]); }
else { k += nVphi; }
563 for(j=0; j<nVphi; j++) {
564 cosPhi = std::cos(phi+j*delPhi/nSphi);
565 sinPhi = std::sin(phi+j*delPhi/nSphi);
566 for(i=i1beg; i<=i2end; i++) {
568 pV[kk[i]+j] =
G4Point3D(r[i]*cosPhi,r[i]*sinPhi,z[i]);
577 v2 = ifClosed ? nodeVis : 1;
578 for(i=i1beg; i<i1end; i++) {
580 if (!ifClosed && i == i1end-1) {
583 v2 = (r[i] == r[i+1] && r[i+1] == r[i+2]) ? -1 : nodeVis;
585 RotateEdge(kk[i], kk[i+1], r[i], r[i+1], v1, v2,
586 edgeVis, ifWholeCircle, nSphi, k);
589 RotateEdge(kk[i1end], kk[i1beg], r[i1end],r[i1beg], nodeVis, nodeVis,
590 edgeVis, ifWholeCircle, nSphi, k);
596 v2 = ifClosed ? nodeVis : 1;
597 for(i=i2beg; i<i2end; i++) {
599 if (!ifClosed && i==i2end-1) {
602 v2 = (r[i] == r[i+1] && r[i+1] == r[i+2]) ? -1 : nodeVis;
604 RotateEdge(kk[i+1], kk[i], r[i+1], r[i], v2, v1,
605 edgeVis, ifWholeCircle, nSphi, k);
608 RotateEdge(kk[i2beg], kk[i2end], r[i2beg], r[i2end], nodeVis, nodeVis,
609 edgeVis, ifWholeCircle, nSphi, k);
617 RotateEdge(kk[i2beg], kk[i1beg], r[i2beg], r[i1beg], 1, 1,
618 -1, ifWholeCircle, nSphi, k);
621 RotateEdge(kk[i1end], kk[i2end], r[i1end], r[i2end], 1, 1,
622 -1, ifWholeCircle, nSphi, k);
628 if (!ifWholeCircle) {
633 for (i=i1beg; i<=i1end; i++) {
635 ii[3] = (i == i1end) ? i1beg : i+1;
636 ii[1] = (absNp2 == 1) ? i2beg : ii[0]+absNp1;
637 ii[2] = (absNp2 == 1) ? i2beg : ii[3]+absNp1;
645 for (i=i1beg; i<i1end; i++) {
648 ii[1] = (absNp2 == 1) ? i2beg : ii[0]+absNp1;
649 ii[2] = (absNp2 == 1) ? i2beg : ii[3]+absNp1;
650 vv[0] = (i == i1beg) ? 1 : -1;
652 vv[2] = (i == i1end-1) ? 1 : -1;
663 <<
"Polyhedron::RotateAroundZ: number of generated faces ("
664 << k-1 <<
") is not equal to the number of allocated faces ("
680 if (
nface <= 0)
return;
682 struct edgeListMember {
683 edgeListMember *next;
687 } *edgeList, *freeList, **headList;
692 edgeList =
new edgeListMember[2*
nface];
693 headList =
new edgeListMember*[
nvert];
696 for (i=0; i<
nvert; i++) {
700 for (i=0; i<2*
nface-1; i++) {
701 edgeList[i].next = &edgeList[i+1];
703 edgeList[2*nface-1].next = 0;
707 G4int iface, iedge, nedge, i1, i2, k1, k2;
708 edgeListMember *prev, *cur;
710 for(iface=1; iface<=
nface; iface++) {
711 nedge = (
pF[iface].edge[3].v == 0) ? 3 : 4;
712 for (iedge=0; iedge<nedge; iedge++) {
714 i2 = (iedge < nedge-1) ? iedge+1 : 0;
715 i1 = std::abs(
pF[iface].edge[i1].
v);
716 i2 = std::abs(
pF[iface].edge[i2].v);
717 k1 = (i1 < i2) ? i1 : i2;
718 k2 = (i1 > i2) ? i1 : i2;
723 headList[k1] = freeList;
724 freeList = freeList->next;
734 headList[k1] = cur->next;
735 cur->next = freeList;
737 pF[iface].edge[iedge].f = cur->iface;
738 pF[cur->iface].edge[cur->iedge].f = iface;
739 i1 = (
pF[iface].edge[iedge].v < 0) ? -1 : 1;
740 i2 = (
pF[cur->iface].edge[cur->iedge].v < 0) ? -1 : 1;
743 <<
"Polyhedron::SetReferences: different edge visibility "
744 << iface <<
"/" << iedge <<
"/"
745 <<
pF[iface].edge[iedge].v <<
" and "
746 << cur->iface <<
"/" << cur->iedge <<
"/"
747 <<
pF[cur->iface].edge[cur->iedge].v
758 prev->next = freeList;
759 freeList = freeList->next;
769 prev->next = cur->next;
770 cur->next = freeList;
772 pF[iface].edge[iedge].f = cur->iface;
773 pF[cur->iface].edge[cur->iedge].f = iface;
774 i1 = (
pF[iface].edge[iedge].v < 0) ? -1 : 1;
775 i2 = (
pF[cur->iface].edge[cur->iedge].v < 0) ? -1 : 1;
778 <<
"Polyhedron::SetReferences: different edge visibility "
779 << iface <<
"/" << iedge <<
"/"
780 <<
pF[iface].edge[iedge].v <<
" and "
781 << cur->iface <<
"/" << cur->iedge <<
"/"
782 <<
pF[cur->iface].edge[cur->iedge].v
793 for (i=0; i<
nvert; i++) {
794 if (headList[i] != 0) {
796 <<
"Polyhedron::SetReferences: List " << i <<
" is not empty"
817 if (
nface <= 0)
return;
819 for (i=1; i<=
nface; i++) {
820 nnode = (
pF[i].edge[3].v == 0) ? 3 : 4;
821 for (k=0; k<nnode; k++) {
822 v[k] = (k+1 == nnode) ?
pF[i].edge[0].v :
pF[i].edge[k+1].v;
823 if (v[k] *
pF[i].edge[k].v < 0) v[k] = -v[k];
824 f[k] =
pF[i].edge[k].f;
826 for (k=0; k<nnode; k++) {
827 pF[i].edge[nnode-1-k].v = v[k];
828 pF[i].edge[nnode-1-k].f = f[k];
870 G4int vIndex = pF[iFace].edge[iQVertex].v;
872 edgeFlag = (vIndex > 0) ? 1 : 0;
873 index = std::abs(vIndex);
875 if (iQVertex >= 3 || pF[iFace].edge[iQVertex+1].
v == 0) {
877 if (++iFace > nface) iFace = 1;
895 if (index <= 0 || index > nvert) {
897 <<
"HepPolyhedron::GetVertex: irrelevant index " <<
index
918 G4bool rep = GetNextVertexIndex(index, edgeFlag);
939 if (nface == 0)
return false;
941 G4int k = pF[iFace].edge[iNode].v;
942 if (k > 0) { edgeFlag = 1; }
else { edgeFlag = -1; k = -k; }
944 normal = FindNodeNormal(iFace,k);
945 if (iNode >= 3 || pF[iFace].edge[iNode+1].
v == 0) {
947 if (++iFace > nface) iFace = 1;
971 G4int k1, k2, kflag, kface1, kface2;
973 if (iFace == 1 && iQVertex == 0) {
974 k2 = pF[nface].edge[0].v;
975 k1 = pF[nface].edge[3].v;
976 if (k1 == 0) k1 = pF[nface].edge[2].v;
977 if (std::abs(k1) > std::abs(k2)) iOrder = -1;
981 k1 = pF[iFace].edge[iQVertex].v;
985 kface2 = pF[iFace].edge[iQVertex].f;
986 if (iQVertex >= 3 || pF[iFace].edge[iQVertex+1].
v == 0) {
988 k2 = std::abs(pF[iFace].edge[iQVertex].
v);
992 k2 = std::abs(pF[iFace].edge[iQVertex].
v);
994 }
while (iOrder*k1 > iOrder*k2);
996 i1 = k1; i2 = k2; edgeFlag = (kflag > 0) ? 1 : 0;
997 iface1 = kface1; iface2 = kface2;
1000 iFace = 1; iOrder = 1;
1019 G4int kface1, kface2;
1020 return GetNextEdgeIndeces(i1, i2, edgeFlag, kface1, kface2);
1026 G4int &edgeFlag)
const
1038 G4bool rep = GetNextEdgeIndeces(i1,i2,edgeFlag);
1059 G4bool rep = GetNextEdgeIndeces(i1,i2,edgeFlag,iface1,iface2);
1076 if (iFace < 1 || iFace > nface) {
1078 <<
"HepPolyhedron::GetFacet: irrelevant index " << iFace
1083 for (i=0; i<4; i++) {
1084 k = pF[iFace].edge[i].v;
1086 if (iFaces != 0) iFaces[i] = pF[iFace].edge[i].f;
1089 if (edgeFlags != 0) edgeFlags[i] = 1;
1092 if (edgeFlags != 0) edgeFlags[i] = -1;
1111 GetFacet(
index,
n, iNodes, edgeFlags);
1113 for (
G4int i=0; i<
n; i++) {
1114 nodes[i] = pV[iNodes[i]];
1115 if (normals != 0) normals[i] = FindNodeNormal(
index,iNodes[i]);
1135 if (edgeFlags == 0) {
1136 GetFacet(iFace,
n, nodes);
1137 }
else if (normals == 0) {
1138 GetFacet(iFace,
n, nodes, edgeFlags);
1140 GetFacet(iFace,
n, nodes, edgeFlags, normals);
1143 if (++iFace > nface) {
1161 if (iFace < 1 || iFace > nface) {
1163 <<
"HepPolyhedron::GetNormal: irrelevant index " << iFace
1168 G4int i0 = std::abs(pF[iFace].edge[0].
v);
1169 G4int i1 = std::abs(pF[iFace].edge[1].v);
1170 G4int i2 = std::abs(pF[iFace].edge[2].v);
1171 G4int i3 = std::abs(pF[iFace].edge[3].v);
1172 if (i3 == 0) i3 = i0;
1173 return (pV[i2] - pV[i0]).cross(pV[i3] - pV[i1]);
1186 if (iFace < 1 || iFace > nface) {
1188 <<
"HepPolyhedron::GetUnitNormal: irrelevant index " << iFace
1193 G4int i0 = std::abs(pF[iFace].edge[0].
v);
1194 G4int i1 = std::abs(pF[iFace].edge[1].v);
1195 G4int i2 = std::abs(pF[iFace].edge[2].v);
1196 G4int i3 = std::abs(pF[iFace].edge[3].v);
1197 if (i3 == 0) i3 = i0;
1198 return ((pV[i2] - pV[i0]).cross(pV[i3] - pV[i1])).unit();
1213 normal = GetNormal(iFace);
1214 if (++iFace > nface) {
1233 G4bool rep = GetNextNormal(normal);
1234 normal = normal.unit();
1250 G4int i0 = std::abs(
pF[iFace].edge[0].
v);
1251 G4int i1 = std::abs(
pF[iFace].edge[1].v);
1252 G4int i2 = std::abs(
pF[iFace].edge[2].v);
1253 G4int i3 = std::abs(
pF[iFace].edge[3].v);
1254 if (i3 == 0) i3 = i0;
1255 srf += ((
pV[i2] -
pV[i0]).cross(
pV[i3] -
pV[i1])).mag();
1272 G4int i0 = std::abs(
pF[iFace].edge[0].v);
1273 G4int i1 = std::abs(
pF[iFace].edge[1].v);
1274 G4int i2 = std::abs(
pF[iFace].edge[2].v);
1275 G4int i3 = std::abs(
pF[iFace].edge[3].v);
1279 pt = (
pV[i0]+
pV[i1]+
pV[i2]) * (1./3.);
1281 pt = (
pV[i0]+
pV[i1]+
pV[i2]+
pV[i3]) * 0.25;
1283 v += ((
pV[i2] -
pV[i0]).cross(
pV[i3] -
pV[i1])).dot(pt);
1323 enum {DUMMY, BOTTOM,
1324 LEFT_BOTTOM, LEFT_FRONT, LEFT_TOP, LEFT_BACK,
1325 BACK_BOTTOM, BACK_LEFT, BACK_TOP, BACK_RIGHT,
1326 RIGHT_BOTTOM, RIGHT_BACK, RIGHT_TOP, RIGHT_FRONT,
1327 FRONT_BOTTOM, FRONT_RIGHT, FRONT_TOP, FRONT_LEFT,
1330 pF[ 1]=
G4Facet(1,LEFT_BOTTOM, 4,BACK_BOTTOM, 3,RIGHT_BOTTOM, 2,FRONT_BOTTOM);
1332 pF[ 2]=
G4Facet(4,BOTTOM, -1,LEFT_FRONT, -12,LEFT_BACK, 0,0);
1333 pF[ 3]=
G4Facet(1,FRONT_LEFT, -5,LEFT_TOP, -12,LEFT_BOTTOM, 0,0);
1334 pF[ 4]=
G4Facet(5,TOP, -8,LEFT_BACK, -12,LEFT_FRONT, 0,0);
1335 pF[ 5]=
G4Facet(8,BACK_LEFT, -4,LEFT_BOTTOM, -12,LEFT_TOP, 0,0);
1337 pF[ 6]=
G4Facet(3,BOTTOM, -4,BACK_LEFT, -11,BACK_RIGHT, 0,0);
1338 pF[ 7]=
G4Facet(4,LEFT_BACK, -8,BACK_TOP, -11,BACK_BOTTOM, 0,0);
1339 pF[ 8]=
G4Facet(8,TOP, -7,BACK_RIGHT, -11,BACK_LEFT, 0,0);
1340 pF[ 9]=
G4Facet(7,RIGHT_BACK, -3,BACK_BOTTOM, -11,BACK_TOP, 0,0);
1342 pF[10]=
G4Facet(2,BOTTOM, -3,RIGHT_BACK, -10,RIGHT_FRONT, 0,0);
1343 pF[11]=
G4Facet(3,BACK_RIGHT, -7,RIGHT_TOP, -10,RIGHT_BOTTOM, 0,0);
1344 pF[12]=
G4Facet(7,TOP, -6,RIGHT_FRONT, -10,RIGHT_BACK, 0,0);
1345 pF[13]=
G4Facet(6,FRONT_RIGHT,-2,RIGHT_BOTTOM,-10,RIGHT_TOP, 0,0);
1347 pF[14]=
G4Facet(1,BOTTOM, -2,FRONT_RIGHT, -9,FRONT_LEFT, 0,0);
1348 pF[15]=
G4Facet(2,RIGHT_FRONT,-6,FRONT_TOP, -9,FRONT_BOTTOM, 0,0);
1349 pF[16]=
G4Facet(6,TOP, -5,FRONT_LEFT, -9,FRONT_RIGHT, 0,0);
1350 pF[17]=
G4Facet(5,LEFT_FRONT, -1,FRONT_BOTTOM, -9,FRONT_TOP, 0,0);
1352 pF[18]=
G4Facet(5,FRONT_TOP, 6,RIGHT_TOP, 7,BACK_TOP, 8,LEFT_TOP);
1360 const G4int faces[][4])
1376 if (
nvert == 0)
return 1;
1378 for (
G4int i=0; i<Nnodes; i++) {
1379 pV[i+1] =
G4Point3D(xyz[i][0], xyz[i][1], xyz[i][2]);
1381 for (
G4int k=0; k<Nfaces; k++) {
1382 pF[k+1] =
G4Facet(faces[k][0],0,faces[k][1],0,faces[k][2],0,faces[k][3],0);
1467 G4double DzTthetaCphi = Dz*std::tan(Theta)*std::cos(Phi);
1468 G4double DzTthetaSphi = Dz*std::tan(Theta)*std::sin(Phi);
1469 G4double Dy1Talp1 = Dy1*std::tan(Alp1);
1470 G4double Dy2Talp2 = Dy2*std::tan(Alp2);
1474 pV[1] =
G4Point3D(-DzTthetaCphi-Dy1Talp1-Dx1,-DzTthetaSphi-Dy1,-Dz);
1475 pV[2] =
G4Point3D(-DzTthetaCphi-Dy1Talp1+Dx1,-DzTthetaSphi-Dy1,-Dz);
1476 pV[3] =
G4Point3D(-DzTthetaCphi+Dy1Talp1+Dx2,-DzTthetaSphi+Dy1,-Dz);
1477 pV[4] =
G4Point3D(-DzTthetaCphi+Dy1Talp1-Dx2,-DzTthetaSphi+Dy1,-Dz);
1478 pV[5] =
G4Point3D( DzTthetaCphi-Dy2Talp2-Dx3, DzTthetaSphi-Dy2, Dz);
1479 pV[6] =
G4Point3D( DzTthetaCphi-Dy2Talp2+Dx3, DzTthetaSphi-Dy2, Dz);
1480 pV[7] =
G4Point3D( DzTthetaCphi+Dy2Talp2+Dx4, DzTthetaSphi+Dy2, Dz);
1481 pV[8] =
G4Point3D( DzTthetaCphi+Dy2Talp2-Dx4, DzTthetaSphi+Dy2, Dz);
1491 :
HepPolyhedronTrap(Dz, Theta, Phi, Dy, Dx, Dx, Alpha, Dy, Dx, Dx, Alpha) {}
1520 if (r1 < 0. || r2 <= 0.) k = 1;
1522 if (dz <= 0.) k += 2;
1528 phi2 = sPhi; phi1 = phi2 + dPhi;
1532 phi1 = sPhi; phi2 = phi1 + wholeCircle;
1536 phi1 = sPhi; phi2 = phi1 + dPhi;
1540 if (std::abs(dphi-wholeCircle) <
perMillion) dphi = wholeCircle;
1541 if (dphi > wholeCircle) k += 4;
1544 std::cerr <<
"HepPolyhedronParaboloid: error in input parameters";
1545 if ((k & 1) != 0) std::cerr <<
" (radiuses)";
1546 if ((k & 2) != 0) std::cerr <<
" (half-length)";
1547 if ((k & 4) != 0) std::cerr <<
" (angles)";
1548 std::cerr << std::endl;
1549 std::cerr <<
" r1=" << r1;
1550 std::cerr <<
" r2=" << r2;
1551 std::cerr <<
" dz=" << dz <<
" sPhi=" << sPhi <<
" dPhi=" << dPhi
1560 G4double k1 = (r2*r2 - r1*r1) / 2 / dz;
1568 for(
G4int i = 1; i < n - 1; i++)
1570 rr[i] = rr[i-1] - dl;
1571 zz[i] = (rr[i]*rr[i] - k2) / k1;
1624 if (r2 < 0. || r1 < 0. ) k = 1;
1625 if (r1 > r2 ) k = 1;
1626 if (r1 == r2) k = 1;
1628 if (halfZ <= 0.) k += 2;
1630 if (sqrtan1<0.||sqrtan2<0.) k += 4;
1634 std::cerr <<
"HepPolyhedronHype: error in input parameters";
1635 if ((k & 1) != 0) std::cerr <<
" (radiuses)";
1636 if ((k & 2) != 0) std::cerr <<
" (half-length)";
1637 if ((k & 4) != 0) std::cerr <<
" (angles)";
1638 std::cerr << std::endl;
1639 std::cerr <<
" r1=" << r1 <<
" r2=" << r2;
1640 std::cerr <<
" halfZ=" << halfZ <<
" sqrTan1=" << sqrtan1
1641 <<
" sqrTan2=" << sqrtan2
1656 rr[0] = std::sqrt(sqrtan2*halfZ*halfZ+k2);
1658 for(
G4int i = 1; i < n-1; i++)
1660 zz[i] = zz[i-1] - dz;
1661 rr[i] =std::sqrt(sqrtan2*zz[i]*zz[i]+k2);
1668 rr[
n] = std::sqrt(sqrtan1*halfZ*halfZ+k1);
1670 for(
G4int i = n+1; i < n+
n; i++)
1672 zz[i] = zz[i-1] - dz;
1673 rr[i] =std::sqrt(sqrtan1*zz[i]*zz[i]+k1);
1716 if (Rmn1 < 0. || Rmx1 < 0. || Rmn2 < 0. || Rmx2 < 0.) k = 1;
1717 if (Rmn1 > Rmx1 || Rmn2 > Rmx2) k = 1;
1718 if (Rmn1 == Rmx1 && Rmn2 == Rmx2) k = 1;
1720 if (Dz <= 0.) k += 2;
1724 phi2 = Phi1; phi1 = phi2 - Dphi;
1725 }
else if (Dphi == 0.) {
1726 phi1 = Phi1; phi2 = phi1 + wholeCircle;
1728 phi1 = Phi1; phi2 = phi1 + Dphi;
1731 if (std::abs(dphi-wholeCircle) <
perMillion) dphi = wholeCircle;
1732 if (dphi > wholeCircle) k += 4;
1735 std::cerr <<
"HepPolyhedronCone(s)/Tube(s): error in input parameters";
1736 if ((k & 1) != 0) std::cerr <<
" (radiuses)";
1737 if ((k & 2) != 0) std::cerr <<
" (half-length)";
1738 if ((k & 4) != 0) std::cerr <<
" (angles)";
1739 std::cerr << std::endl;
1740 std::cerr <<
" Rmn1=" << Rmn1 <<
" Rmx1=" << Rmx1;
1741 std::cerr <<
" Rmn2=" << Rmn2 <<
" Rmx2=" << Rmx2;
1742 std::cerr <<
" Dz=" << Dz <<
" Phi1=" << Phi1 <<
" Dphi=" << Dphi
1813 if (dphi <= 0. || dphi >
twopi) {
1815 <<
"HepPolyhedronPgon/Pcon: wrong delta phi = " << dphi
1822 <<
"HepPolyhedronPgon/Pcon: number of z-planes less than two = " << nz
1829 <<
"HepPolyhedronPgon/Pcon: error in number of phi-steps =" << npdv
1835 for (i=0; i<nz; i++) {
1836 if (rmin[i] < 0. || rmax[i] < 0. || rmin[i] > rmax[i]) {
1838 <<
"HepPolyhedronPgon: error in radiuses rmin[" << i <<
"]="
1839 << rmin[i] <<
" rmax[" << i <<
"]=" << rmax[i]
1851 if (z[0] > z[nz-1]) {
1852 for (i=0; i<nz; i++) {
1859 for (i=0; i<nz; i++) {
1861 rr[i] = rmax[nz-i-1];
1862 zz[i+nz] = z[nz-i-1];
1863 rr[i+nz] = rmin[nz-i-1];
1869 RotateAroundZ(npdv, phi, dphi, nz, nz, zz, rr, -1, (npdv == 0) ? -1 : 1);
1907 if (dphi <= 0. || dphi >
twopi) {
1909 <<
"HepPolyhedronSphere: wrong delta phi = " << dphi
1914 if (the < 0. || the >
pi) {
1916 <<
"HepPolyhedronSphere: wrong theta = " << the
1921 if (dthe <= 0. || dthe > pi) {
1923 <<
"HepPolyhedronSphere: wrong delta theta = " << dthe
1928 if (the+dthe > pi) {
1930 <<
"HepPolyhedronSphere: wrong theta + delta theta = "
1931 << the <<
" " << dthe
1936 if (rmin < 0. || rmin >= rmax) {
1938 <<
"HepPolyhedronSphere: error in radiuses"
1939 <<
" rmin=" << rmin <<
" rmax=" << rmax
1948 if (np1 <= 1) np1 = 2;
1957 for (
G4int i=0; i<np1; i++) {
1958 cosa = std::cos(the+i*a);
1959 sina = std::sin(the+i*a);
1963 zz[i+np1] = rmin*cosa;
1964 rr[i+np1] = rmin*sina;
2005 if (dphi <= 0. || dphi >
twopi) {
2007 <<
"HepPolyhedronTorus: wrong delta phi = " << dphi
2012 if (rmin < 0. || rmin >= rmax || rmax >= rtor) {
2014 <<
"HepPolyhedronTorus: error in radiuses"
2015 <<
" rmin=" << rmin <<
" rmax=" << rmax <<
" rtorus=" << rtor
2031 for (
G4int i=0; i<np1; i++) {
2032 cosa = std::cos(i*a);
2033 sina = std::sin(i*a);
2035 rr[i] = rtor+rmax*sina;
2037 zz[i+np1] = rmin*cosa;
2038 rr[i+np1] = rtor+rmin*sina;
2078 if (zCut1 >= cz || zCut2 <= -cz || zCut1 > zCut2) {
2079 std::cerr <<
"HepPolyhedronEllipsoid: wrong zCut1 = " << zCut1
2080 <<
" zCut2 = " << zCut2
2081 <<
" for given cz = " << cz << std::endl;
2085 std::cerr <<
"HepPolyhedronEllipsoid: bad z semi-axis: cz = " << cz
2100 sthe= std::acos(zCut2/cz);
2109 dthe= std::acos(zCut1/cz)-sthe;
2124 G4Exception(
"HepPolyhedronEllipsoid::HepPolyhedronEllipsoid",
2137 for (
G4int i=0; i<np1-cutflag; i++) {
2138 cosa = std::cos(sthe+i*a);
2139 sina = std::sin(sthe+i*a);
2152 std::cerr <<
"Logic error in HepPolyhedronEllipsoid, memory corrupted!"
2157 std::cerr <<
"Warning: logic error in HepPolyhedronEllipsoid."
2177 p->
setX( p->
x() * ax/cz );
2178 p->
setY( p->
y() * by/cz );
2205 if ( (ax <= 0.) || (ay <= 0.) || (h <= 0.) || (zTopCut <= 0.) ) { k = 1; }
2208 std::cerr <<
"HepPolyhedronCone: error in input parameters";
2209 std::cerr << std::endl;
2215 zTopCut = (h >= zTopCut ? zTopCut : h);
2224 rr[0] = (h-zTopCut);
2225 rr[1] = (h+zTopCut);
2241 p->
setX( p->
x() * ax );
2242 p->
setY( p->
y() * ay );
2259 #include "BooleanProcessor.src"
2273 return processor.execute(OP_UNION, *
this,
p,ierr);
2288 return processor.execute(OP_INTERSECTION, *
this,
p,ierr);
2303 return processor.execute(OP_SUBTRACTION, *
this,
p,ierr);
2311 #include "HepPolyhedronProcessor.src"
virtual ~HepPolyhedronTorus()
#define DEFAULT_NUMBER_OF_STEPS
void AllocateMemory(G4int Nvert, G4int Nface)
const G4double spatialTolerance
void GetFacet(G4int iFace, G4int &n, G4int *iNodes, G4int *edgeFlags=0, G4int *iFaces=0) const
HepPolyhedronEllipticalCone(G4double dx, G4double dy, G4double z, G4double zcut1)
HepPolyhedronPgon(G4double phi, G4double dphi, G4int npdv, G4int nz, const G4double *z, const G4double *rmin, const G4double *rmax)
G4Normal3D GetUnitNormal(G4int iFace) const
virtual ~HepPolyhedronSphere()
G4int createTwistedTrap(G4double Dz, const G4double xy1[][2], const G4double xy2[][2])
HepPolyhedronParaboloid(G4double r1, G4double r2, G4double dz, G4double Phi1, G4double Dphi)
virtual ~HepPolyhedronEllipsoid()
G4bool GetNextUnitNormal(G4Normal3D &normal) const
void SetSideFacets(G4int ii[4], G4int vv[4], G4int *kk, G4double *r, G4double dphi, G4int ns, G4int &kface)
G4bool GetNextNormal(G4Normal3D &normal) const
G4int createPolyhedron(G4int Nnodes, G4int Nfaces, const G4double xyz[][3], const G4int faces[][4])
for(Int_t i=0;i< nentries;i++)
G4bool GetNextEdgeIndeces(G4int &i1, G4int &i2, G4int &edgeFlag, G4int &iface1, G4int &iface2) const
G4bool GetNextEdge(G4Point3D &p1, G4Point3D &p2, G4int &edgeFlag) const
BasicVector3D< T > unit() const
virtual ~HepPolyhedronPara()
virtual ~HepPolyhedronCone()
virtual ~HepPolyhedronTrap()
HepGeom::Point3D< G4double > G4Point3D
HepPolyhedronTrd1(G4double Dx1, G4double Dx2, G4double Dy, G4double Dz)
HepGeom::Vector3D< G4double > G4Vector3D
HepPolyhedron subtract(const HepPolyhedron &p) const
HepPolyhedronSphere(G4double rmin, G4double rmax, G4double phi, G4double dphi, G4double the, G4double dthe)
HepPolyhedron intersect(const HepPolyhedron &p) const
G4Normal3D FindNodeNormal(G4int iFace, G4int iNode) const
const G4ThreeVector const G4double const
HepPolyhedron & Transform(const G4Transform3D &t)
HepPolyhedronTorus(G4double rmin, G4double rmax, G4double rtor, G4double phi, G4double dphi)
virtual ~HepPolyhedronTubs()
HepPolyhedronEllipsoid(G4double dx, G4double dy, G4double dz, G4double zcut1, G4double zcut2)
G4int FindNeighbour(G4int iFace, G4int iNode, G4int iOrder) const
virtual ~HepPolyhedronHype()
HepPolyhedronCone(G4double Rmn1, G4double Rmx1, G4double Rmn2, G4double Rmx2, G4double Dz)
HepPolyhedronCons(G4double Rmn1, G4double Rmx1, G4double Rmn2, G4double Rmx2, G4double Dz, G4double Phi1, G4double Dphi)
static G4ThreadLocal G4int fNumberOfRotationSteps
HepPolyhedron add(const HepPolyhedron &p) const
G4Normal3D GetNormal(G4int iFace) const
virtual ~HepPolyhedronTube()
G4bool GetNextFacet(G4int &n, G4Point3D *nodes, G4int *edgeFlags=0, G4Normal3D *normals=0) const
HepPolyhedronPara(G4double Dx, G4double Dy, G4double Dz, G4double Alpha, G4double Theta, G4double Phi)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
HepPolyhedronHype(G4double r1, G4double r2, G4double tan1, G4double tan2, G4double halfZ)
virtual ~HepPolyhedronTrd2()
HepPolyhedronBox(G4double Dx, G4double Dy, G4double Dz)
void RotateAroundZ(G4int nstep, G4double phi, G4double dphi, G4int np1, G4int np2, const G4double *z, G4double *r, G4int nodeVis, G4int edgeVis)
virtual ~HepPolyhedronEllipticalCone()
static G4int GetNumberOfRotationSteps()
HepPolyhedronPcon(G4double phi, G4double dphi, G4int nz, const G4double *z, const G4double *rmin, const G4double *rmax)
virtual ~HepPolyhedronBox()
HepPolyhedronTubs(G4double Rmin, G4double Rmax, G4double Dz, G4double Phi1, G4double Dphi)
HepPolyhedronTrap(G4double Dz, G4double Theta, G4double Phi, G4double Dy1, G4double Dx1, G4double Dx2, G4double Alp1, G4double Dy2, G4double Dx3, G4double Dx4, G4double Alp2)
std::ostream & operator<<(std::ostream &, const BasicVector3D< float > &)
HepPolyhedron & operator=(const HepPolyhedron &from)
G4bool GetNextVertexIndex(G4int &index, G4int &edgeFlag) const
HepPolyhedronTrd2(G4double Dx1, G4double Dx2, G4double Dy1, G4double Dy2, G4double Dz)
void RotateEdge(G4int k1, G4int k2, G4double r1, G4double r2, G4int v1, G4int v2, G4int vEdge, G4bool ifWholeCircle, G4int ns, G4int &kface)
virtual ~HepPolyhedronPgon()
static void SetNumberOfRotationSteps(G4int n)
virtual ~HepPolyhedronTrd1()
G4double GetVolume() const
G4bool GetNextVertex(G4Point3D &vertex, G4int &edgeFlag) const
virtual ~HepPolyhedronCons()
G4double GetSurfaceArea() const
BasicVector3D< T > cross(const BasicVector3D< T > &v) const
HepPolyhedronTube(G4double Rmin, G4double Rmax, G4double Dz)
G4Point3D GetVertex(G4int index) const
virtual ~HepPolyhedronPcon()
HepGeom::Normal3D< G4double > G4Normal3D
static void ResetNumberOfRotationSteps()
virtual ~HepPolyhedronParaboloid()