00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045 #include "G4TwistTubsSide.hh"
00046
00047
00048
00049
00050 G4TwistTubsSide::G4TwistTubsSide(const G4String &name,
00051 const G4RotationMatrix &rot,
00052 const G4ThreeVector &tlate,
00053 G4int handedness,
00054 const G4double kappa,
00055 const EAxis axis0,
00056 const EAxis axis1,
00057 G4double axis0min,
00058 G4double axis1min,
00059 G4double axis0max,
00060 G4double axis1max)
00061 : G4VTwistSurface(name, rot, tlate, handedness, axis0, axis1,
00062 axis0min, axis1min, axis0max, axis1max),
00063 fKappa(kappa)
00064 {
00065 if (axis0 == kZAxis && axis1 == kXAxis) {
00066 G4Exception("G4TwistTubsSide::G4TwistTubsSide()", "GeomSolids0002",
00067 FatalErrorInArgument, "Should swap axis0 and axis1!");
00068 }
00069 fIsValidNorm = false;
00070 SetCorners();
00071 SetBoundaries();
00072 }
00073
00074 G4TwistTubsSide::G4TwistTubsSide(const G4String &name,
00075 G4double EndInnerRadius[2],
00076 G4double EndOuterRadius[2],
00077 G4double DPhi,
00078 G4double EndPhi[2],
00079 G4double EndZ[2],
00080 G4double InnerRadius,
00081 G4double OuterRadius,
00082 G4double Kappa,
00083 G4int handedness)
00084 : G4VTwistSurface(name)
00085 {
00086 fHandedness = handedness;
00087 fAxis[0] = kXAxis;
00088 fAxis[1] = kZAxis;
00089 fAxisMin[0] = InnerRadius;
00090 fAxisMax[0] = OuterRadius;
00091 fAxisMin[1] = EndZ[0];
00092 fAxisMax[1] = EndZ[1];
00093
00094 fKappa = Kappa;
00095 fRot.rotateZ( fHandedness > 0
00096 ? -0.5*DPhi
00097 : 0.5*DPhi );
00098 fTrans.set(0, 0, 0);
00099 fIsValidNorm = false;
00100
00101 SetCorners( EndInnerRadius, EndOuterRadius, EndPhi, EndZ) ;
00102 SetBoundaries();
00103 }
00104
00105
00106
00107
00108
00109 G4TwistTubsSide::G4TwistTubsSide( __void__& a )
00110 : G4VTwistSurface(a), fKappa(0.)
00111 {
00112 }
00113
00114
00115
00116
00117
00118 G4TwistTubsSide::~G4TwistTubsSide()
00119 {
00120 }
00121
00122
00123
00124
00125 G4ThreeVector G4TwistTubsSide::GetNormal(const G4ThreeVector &tmpxx,
00126 G4bool isGlobal)
00127 {
00128
00129
00130
00131
00132 G4ThreeVector xx;
00133 if (isGlobal) {
00134 xx = ComputeLocalPoint(tmpxx);
00135 if ((xx - fCurrentNormal.p).mag() < 0.5 * kCarTolerance) {
00136 return ComputeGlobalDirection(fCurrentNormal.normal);
00137 }
00138 } else {
00139 xx = tmpxx;
00140 if (xx == fCurrentNormal.p) {
00141 return fCurrentNormal.normal;
00142 }
00143 }
00144
00145 G4ThreeVector er(1, fKappa * xx.z(), 0);
00146 G4ThreeVector ez(0, fKappa * xx.x(), 1);
00147 G4ThreeVector normal = fHandedness*(er.cross(ez));
00148
00149 if (isGlobal) {
00150 fCurrentNormal.normal = ComputeGlobalDirection(normal.unit());
00151 } else {
00152 fCurrentNormal.normal = normal.unit();
00153 }
00154 return fCurrentNormal.normal;
00155 }
00156
00157
00158
00159
00160 G4int G4TwistTubsSide::DistanceToSurface(const G4ThreeVector &gp,
00161 const G4ThreeVector &gv,
00162 G4ThreeVector gxx[],
00163 G4double distance[],
00164 G4int areacode[],
00165 G4bool isvalid[],
00166 EValidate validate)
00167 {
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212 fCurStatWithV.ResetfDone(validate, &gp, &gv);
00213
00214 if (fCurStatWithV.IsDone()) {
00215 G4int i;
00216 for (i=0; i<fCurStatWithV.GetNXX(); i++) {
00217 gxx[i] = fCurStatWithV.GetXX(i);
00218 distance[i] = fCurStatWithV.GetDistance(i);
00219 areacode[i] = fCurStatWithV.GetAreacode(i);
00220 isvalid[i] = fCurStatWithV.IsValid(i);
00221 }
00222 return fCurStatWithV.GetNXX();
00223 } else {
00224
00225 G4int i;
00226 for (i=0; i<2; i++) {
00227 distance[i] = kInfinity;
00228 areacode[i] = sOutside;
00229 isvalid[i] = false;
00230 gxx[i].set(kInfinity, kInfinity, kInfinity);
00231 }
00232 }
00233
00234 G4ThreeVector p = ComputeLocalPoint(gp);
00235 G4ThreeVector v = ComputeLocalDirection(gv);
00236 G4ThreeVector xx[2];
00237
00238
00239
00240
00241
00242
00243
00244 G4double absvz = std::fabs(v.z());
00245
00246 if ((absvz < DBL_MIN) && (std::fabs(p.x() * v.y() - p.y() * v.x()) < DBL_MIN)) {
00247
00248
00249 isvalid[0] = false;
00250 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
00251 isvalid[0], 0, validate, &gp, &gv);
00252 return 0;
00253 }
00254
00255
00256
00257
00258
00259
00260 G4double a = fKappa * v.x() * v.z();
00261 G4double b = fKappa * (v.x() * p.z() + v.z() * p.x()) - v.y();
00262 G4double c = fKappa * p.x() * p.z() - p.y();
00263 G4double D = b * b - 4 * a * c;
00264 G4int vout = 0;
00265
00266 if (std::fabs(a) < DBL_MIN) {
00267 if (std::fabs(b) > DBL_MIN) {
00268
00269
00270
00271 distance[0] = - c / b;
00272 xx[0] = p + distance[0]*v;
00273 gxx[0] = ComputeGlobalPoint(xx[0]);
00274
00275 if (validate == kValidateWithTol) {
00276 areacode[0] = GetAreaCode(xx[0]);
00277 if (!IsOutside(areacode[0])) {
00278 if (distance[0] >= 0) isvalid[0] = true;
00279 }
00280 } else if (validate == kValidateWithoutTol) {
00281 areacode[0] = GetAreaCode(xx[0], false);
00282 if (IsInside(areacode[0])) {
00283 if (distance[0] >= 0) isvalid[0] = true;
00284 }
00285 } else {
00286
00287 if (xx[0].x() > 0) {
00288 areacode[0] = sInside;
00289 if (distance[0] >= 0) isvalid[0] = true;
00290 } else {
00291 distance[0] = kInfinity;
00292 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0],
00293 areacode[0], isvalid[0],
00294 0, validate, &gp, &gv);
00295 return vout;
00296 }
00297 }
00298
00299 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
00300 isvalid[0], 1, validate, &gp, &gv);
00301 vout = 1;
00302
00303 } else {
00304
00305
00306
00307
00308
00309
00310
00311 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
00312 isvalid[0], 0, validate, &gp, &gv);
00313 }
00314
00315 } else if (D > DBL_MIN) {
00316
00317
00318
00319 D = std::sqrt(D);
00320 G4double factor = 0.5/a;
00321 G4double tmpdist[2] = {kInfinity, kInfinity};
00322 G4ThreeVector tmpxx[2];
00323 G4int tmpareacode[2] = {sOutside, sOutside};
00324 G4bool tmpisvalid[2] = {false, false};
00325 G4int i;
00326
00327 for (i=0; i<2; i++) {
00328 G4double bminusD = - b - D;
00329
00330
00331
00332 G4double protection = 0;
00333 if ( b * D < 0 && std::fabs(bminusD / D) < protection ) {
00334 G4double acovbb = (a*c)/(b*b);
00335 tmpdist[i] = - c/b * ( 1 - acovbb * (1 + 2*acovbb));
00336 } else {
00337 tmpdist[i] = factor * bminusD;
00338 }
00339
00340 D = -D;
00341 tmpxx[i] = p + tmpdist[i]*v;
00342
00343 if (validate == kValidateWithTol) {
00344 tmpareacode[i] = GetAreaCode(tmpxx[i]);
00345 if (!IsOutside(tmpareacode[i])) {
00346 if (tmpdist[i] >= 0) tmpisvalid[i] = true;
00347 continue;
00348 }
00349 } else if (validate == kValidateWithoutTol) {
00350 tmpareacode[i] = GetAreaCode(tmpxx[i], false);
00351 if (IsInside(tmpareacode[i])) {
00352 if (tmpdist[i] >= 0) tmpisvalid[i] = true;
00353 continue;
00354 }
00355 } else {
00356
00357 if (tmpxx[i].x() > 0) {
00358 tmpareacode[i] = sInside;
00359 if (tmpdist[i] >= 0) tmpisvalid[i] = true;
00360 continue;
00361 } else {
00362 tmpdist[i] = kInfinity;
00363 continue;
00364 }
00365 }
00366 }
00367
00368 if (tmpdist[0] <= tmpdist[1]) {
00369 distance[0] = tmpdist[0];
00370 distance[1] = tmpdist[1];
00371 xx[0] = tmpxx[0];
00372 xx[1] = tmpxx[1];
00373 gxx[0] = ComputeGlobalPoint(tmpxx[0]);
00374 gxx[1] = ComputeGlobalPoint(tmpxx[1]);
00375 areacode[0] = tmpareacode[0];
00376 areacode[1] = tmpareacode[1];
00377 isvalid[0] = tmpisvalid[0];
00378 isvalid[1] = tmpisvalid[1];
00379 } else {
00380 distance[0] = tmpdist[1];
00381 distance[1] = tmpdist[0];
00382 xx[0] = tmpxx[1];
00383 xx[1] = tmpxx[0];
00384 gxx[0] = ComputeGlobalPoint(tmpxx[1]);
00385 gxx[1] = ComputeGlobalPoint(tmpxx[0]);
00386 areacode[0] = tmpareacode[1];
00387 areacode[1] = tmpareacode[0];
00388 isvalid[0] = tmpisvalid[1];
00389 isvalid[1] = tmpisvalid[0];
00390 }
00391
00392 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
00393 isvalid[0], 2, validate, &gp, &gv);
00394 fCurStatWithV.SetCurrentStatus(1, gxx[1], distance[1], areacode[1],
00395 isvalid[1], 2, validate, &gp, &gv);
00396
00397
00398
00399 for (G4int k=0; k<2; k++) {
00400 if (!isvalid[k]) continue;
00401
00402 G4ThreeVector xxonsurface(xx[k].x(), fKappa * std::fabs(xx[k].x())
00403 * xx[k].z() , xx[k].z());
00404 G4double deltaY = (xx[k] - xxonsurface).mag();
00405
00406 if ( deltaY > 0.5*kCarTolerance ) {
00407
00408 G4int maxcount = 10;
00409 G4int l;
00410 G4double lastdeltaY = deltaY;
00411 for (l=0; l<maxcount; l++) {
00412 G4ThreeVector surfacenormal = GetNormal(xxonsurface);
00413 distance[k] = DistanceToPlaneWithV(p, v, xxonsurface,
00414 surfacenormal, xx[k]);
00415 deltaY = (xx[k] - xxonsurface).mag();
00416 if (deltaY > lastdeltaY) {
00417
00418 }
00419 gxx[k] = ComputeGlobalPoint(xx[k]);
00420
00421 if (deltaY <= 0.5*kCarTolerance) {
00422
00423 break;
00424 }
00425 xxonsurface.set(xx[k].x(),
00426 fKappa * std::fabs(xx[k].x()) * xx[k].z(),
00427 xx[k].z());
00428 }
00429 if (l == maxcount) {
00430 std::ostringstream message;
00431 message << "Exceeded maxloop count!" << G4endl
00432 << " maxloop count " << maxcount;
00433 G4Exception("G4TwistTubsFlatSide::DistanceToSurface(p,v)",
00434 "GeomSolids0003", FatalException, message);
00435 }
00436 }
00437
00438 }
00439 vout = 2;
00440 } else {
00441
00442
00443
00444 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
00445 isvalid[0], 0, validate, &gp, &gv);
00446 }
00447
00448 return vout;
00449 }
00450
00451
00452
00453
00454 G4int G4TwistTubsSide::DistanceToSurface(const G4ThreeVector &gp,
00455 G4ThreeVector gxx[],
00456 G4double distance[],
00457 G4int areacode[])
00458 {
00459 fCurStat.ResetfDone(kDontValidate, &gp);
00460 G4int i = 0;
00461 if (fCurStat.IsDone()) {
00462 for (i=0; i<fCurStat.GetNXX(); i++) {
00463 gxx[i] = fCurStat.GetXX(i);
00464 distance[i] = fCurStat.GetDistance(i);
00465 areacode[i] = fCurStat.GetAreacode(i);
00466 }
00467 return fCurStat.GetNXX();
00468 } else {
00469
00470 for (i=0; i<2; i++) {
00471 distance[i] = kInfinity;
00472 areacode[i] = sOutside;
00473 gxx[i].set(kInfinity, kInfinity, kInfinity);
00474 }
00475 }
00476
00477 static const G4double halftol = 0.5 * kCarTolerance;
00478
00479 G4ThreeVector p = ComputeLocalPoint(gp);
00480 G4ThreeVector xx;
00481 G4int parity = (fKappa >= 0 ? 1 : -1);
00482
00483
00484
00485
00486
00487
00488
00489
00490 G4ThreeVector lastgxx[2];
00491 for (i=0; i<2; i++) {
00492 lastgxx[i] = fCurStatWithV.GetXX(i);
00493 }
00494
00495 if ((gp - lastgxx[0]).mag() < halftol
00496 || (gp - lastgxx[1]).mag() < halftol) {
00497
00498 xx = p;
00499 distance[0] = 0;
00500 gxx[0] = gp;
00501
00502 G4bool isvalid = true;
00503 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
00504 isvalid, 1, kDontValidate, &gp);
00505 return 1;
00506 }
00507
00508 if (p.getRho() == 0) {
00509
00510
00511
00512
00513 G4bool isvalid = true;
00514 if (fAxis[0] == kXAxis && fAxis[1] == kZAxis) {
00515 distance[0] = DistanceToBoundary(sAxis0 & sAxisMin, xx, p);
00516 areacode[0] = sInside;
00517 } else {
00518 distance[0] = 0;
00519 xx.set(0., 0., 0.);
00520 }
00521 gxx[0] = ComputeGlobalPoint(xx);
00522 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
00523 isvalid, 0, kDontValidate, &gp);
00524 return 1;
00525 }
00526
00527
00528
00529
00530
00531
00532
00533 G4ThreeVector A;
00534 G4ThreeVector C;
00535 G4ThreeVector B;
00536 G4ThreeVector D;
00537
00538
00539 DistanceToBoundary(sAxis0 & sAxisMin, A, p);
00540
00541 DistanceToBoundary(sAxis0 & sAxisMax, C, p);
00542
00543
00544
00545 if (A.z() > C.z()) {
00546 if (p.z() > A.z()) {
00547 A = GetBoundaryAtPZ(sAxis0 & sAxisMin, p);
00548 } else if (p.z() < C.z()) {
00549 C = GetBoundaryAtPZ(sAxis0 & sAxisMax, p);
00550 }
00551 } else {
00552 if (p.z() > C.z()) {
00553 C = GetBoundaryAtPZ(sAxis0 & sAxisMax, p);
00554 } else if (p.z() < A.z()) {
00555 A = GetBoundaryAtPZ(sAxis0 & sAxisMin, p);
00556 }
00557 }
00558
00559 G4ThreeVector d[2];
00560 G4ThreeVector x0[2];
00561 G4int btype[2];
00562
00563 for (i=0; i<2; i++) {
00564 if (i == 0) {
00565 GetBoundaryParameters((sAxis0 & sAxisMax), d[i], x0[i], btype[i]);
00566 B = x0[i] + ((A.z() - x0[i].z()) / d[i].z()) * d[i];
00567
00568 } else {
00569 GetBoundaryParameters((sAxis0 & sAxisMin), d[i], x0[i], btype[i]);
00570 D = x0[i] + ((C.z() - x0[i].z()) / d[i].z()) * d[i];
00571 }
00572 }
00573
00574
00575 G4ThreeVector pt(p.x(), p.y(), 0.);
00576 G4double rc = std::fabs(p.x());
00577 G4ThreeVector surfacevector(rc, rc * fKappa * p.z(), 0.);
00578 G4int pside = AmIOnLeftSide(pt, surfacevector);
00579 G4double test = (A.z() - C.z()) * parity * pside;
00580
00581 if (test == 0) {
00582 if (pside == 0) {
00583
00584 xx = p;
00585 distance[0] = 0;
00586 gxx[0] = gp;
00587
00588 G4bool isvalid = true;
00589 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
00590 isvalid, 1, kDontValidate, &gp);
00591 return 1;
00592 } else {
00593
00594 d[0] = C - A;
00595 distance[0] = DistanceToLine(p, A, d[0], xx);
00596 areacode[0] = sInside;
00597 gxx[0] = ComputeGlobalPoint(xx);
00598 G4bool isvalid = true;
00599 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
00600 isvalid, 1, kDontValidate, &gp);
00601 return 1;
00602 }
00603
00604 } else if (test < 0) {
00605
00606
00607
00608 G4ThreeVector tmp;
00609 tmp = A;
00610 A = D;
00611 D = tmp;
00612 tmp = C;
00613 C = B;
00614 B = tmp;
00615
00616 } else {
00617
00618 }
00619
00620
00621
00622
00623
00624 G4ThreeVector xxacb;
00625 G4ThreeVector nacb;
00626 G4ThreeVector xxcad;
00627 G4ThreeVector ncad;
00628 G4ThreeVector AB(A.x(), A.y(), 0);
00629 G4ThreeVector DC(C.x(), C.y(), 0);
00630
00631 G4double distToACB = G4VTwistSurface::DistanceToPlane(p, A, C-A, AB, xxacb, nacb) * parity;
00632 G4double distToCAD = G4VTwistSurface::DistanceToPlane(p, C, C-A, DC, xxcad, ncad) * parity;
00633
00634
00635
00636 if (std::fabs(distToACB) <= halftol || std::fabs(distToCAD) <= halftol) {
00637 xx = (std::fabs(distToACB) < std::fabs(distToCAD) ? xxacb : xxcad);
00638 areacode[0] = sInside;
00639 gxx[0] = ComputeGlobalPoint(xx);
00640 distance[0] = 0;
00641 G4bool isvalid = true;
00642 fCurStat.SetCurrentStatus(0, gxx[0], distance[0] , areacode[0],
00643 isvalid, 1, kDontValidate, &gp);
00644 return 1;
00645 }
00646
00647 if (distToACB * distToCAD > 0 && distToACB < 0) {
00648
00649
00650 G4ThreeVector normal;
00651 distance[0] = DistanceToPlane(p, A, B, C, D, parity, xx, normal);
00652 } else {
00653 if (distToACB * distToCAD > 0) {
00654
00655
00656 if (distToACB <= distToCAD) {
00657 distance[0] = distToACB;
00658 xx = xxacb;
00659 } else {
00660 distance[0] = distToCAD;
00661 xx = xxcad;
00662 }
00663 } else {
00664
00665
00666 if (distToACB > 0) {
00667 distance[0] = distToACB;
00668 xx = xxacb;
00669 } else {
00670 distance[0] = distToCAD;
00671 xx = xxcad;
00672 }
00673 }
00674
00675 }
00676 areacode[0] = sInside;
00677 gxx[0] = ComputeGlobalPoint(xx);
00678 G4bool isvalid = true;
00679 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
00680 isvalid, 1, kDontValidate, &gp);
00681 return 1;
00682 }
00683
00684
00685
00686
00687 G4double G4TwistTubsSide::DistanceToPlane(const G4ThreeVector &p,
00688 const G4ThreeVector &A,
00689 const G4ThreeVector &B,
00690 const G4ThreeVector &C,
00691 const G4ThreeVector &D,
00692 const G4int parity,
00693 G4ThreeVector &xx,
00694 G4ThreeVector &n)
00695 {
00696 static const G4double halftol = 0.5 * kCarTolerance;
00697
00698 G4ThreeVector M = 0.5*(A + B);
00699 G4ThreeVector N = 0.5*(C + D);
00700 G4ThreeVector xxanm;
00701 G4ThreeVector nanm;
00702 G4ThreeVector xxcmn;
00703 G4ThreeVector ncmn;
00704
00705 G4double distToanm = G4VTwistSurface::DistanceToPlane(p, A, (N - A), (M - A), xxanm, nanm) * parity;
00706 G4double distTocmn = G4VTwistSurface::DistanceToPlane(p, C, (M - C), (N - C), xxcmn, ncmn) * parity;
00707
00708
00709 if (distToanm * distTocmn > 0 && distToanm < 0) {
00710 G4Exception("G4TwistTubsSide::DistanceToPlane()",
00711 "GeomSolids0003", FatalException,
00712 "Point p is behind the surfaces.");
00713 }
00714
00715
00716 if (std::fabs(distToanm) <= halftol) {
00717 xx = xxanm;
00718 n = nanm * parity;
00719 return 0;
00720 } else if (std::fabs(distTocmn) <= halftol) {
00721 xx = xxcmn;
00722 n = ncmn * parity;
00723 return 0;
00724 }
00725
00726 if (distToanm <= distTocmn) {
00727 if (distToanm > 0) {
00728
00729 xx = xxanm;
00730 n = nanm * parity;
00731 return distToanm;
00732 } else {
00733
00734 return DistanceToPlane(p, A, M, N, D, parity, xx, n);
00735 }
00736 } else {
00737 if (distTocmn > 0) {
00738
00739 xx = xxcmn;
00740 n = ncmn * parity;
00741 return distTocmn;
00742 } else {
00743
00744 return DistanceToPlane(p, C, N, M, B, parity, xx, n);
00745 }
00746 }
00747 }
00748
00749
00750
00751
00752 G4int G4TwistTubsSide::GetAreaCode(const G4ThreeVector &xx,
00753 G4bool withTol)
00754 {
00755
00756
00757
00758 static const G4double ctol = 0.5 * kCarTolerance;
00759 G4int areacode = sInside;
00760
00761 if (fAxis[0] == kXAxis && fAxis[1] == kZAxis) {
00762 G4int xaxis = 0;
00763 G4int zaxis = 1;
00764
00765 if (withTol) {
00766
00767 G4bool isoutside = false;
00768
00769
00770
00771 if (xx.x() < fAxisMin[xaxis] + ctol) {
00772 areacode |= (sAxis0 & (sAxisX | sAxisMin)) | sBoundary;
00773 if (xx.x() <= fAxisMin[xaxis] - ctol) isoutside = true;
00774
00775 } else if (xx.x() > fAxisMax[xaxis] - ctol) {
00776 areacode |= (sAxis0 & (sAxisX | sAxisMax)) | sBoundary;
00777 if (xx.x() >= fAxisMax[xaxis] + ctol) isoutside = true;
00778 }
00779
00780
00781
00782 if (xx.z() < fAxisMin[zaxis] + ctol) {
00783 areacode |= (sAxis1 & (sAxisZ | sAxisMin));
00784
00785 if (areacode & sBoundary) areacode |= sCorner;
00786 else areacode |= sBoundary;
00787 if (xx.z() <= fAxisMin[zaxis] - ctol) isoutside = true;
00788
00789 } else if (xx.z() > fAxisMax[zaxis] - ctol) {
00790 areacode |= (sAxis1 & (sAxisZ | sAxisMax));
00791
00792 if (areacode & sBoundary) areacode |= sCorner;
00793 else areacode |= sBoundary;
00794 if (xx.z() >= fAxisMax[zaxis] + ctol) isoutside = true;
00795 }
00796
00797
00798
00799
00800 if (isoutside) {
00801 G4int tmpareacode = areacode & (~sInside);
00802 areacode = tmpareacode;
00803 } else if ((areacode & sBoundary) != sBoundary) {
00804 areacode |= (sAxis0 & sAxisX) | (sAxis1 & sAxisZ);
00805 }
00806
00807 } else {
00808
00809
00810
00811 if (xx.x() < fAxisMin[xaxis] ) {
00812 areacode |= (sAxis0 & (sAxisX | sAxisMin)) | sBoundary;
00813 } else if (xx.x() > fAxisMax[xaxis]) {
00814 areacode |= (sAxis0 & (sAxisX | sAxisMax)) | sBoundary;
00815 }
00816
00817
00818
00819 if (xx.z() < fAxisMin[zaxis]) {
00820 areacode |= (sAxis1 & (sAxisZ | sAxisMin));
00821 if (areacode & sBoundary) areacode |= sCorner;
00822 else areacode |= sBoundary;
00823
00824 } else if (xx.z() > fAxisMax[zaxis]) {
00825 areacode |= (sAxis1 & (sAxisZ | sAxisMax)) ;
00826 if (areacode & sBoundary) areacode |= sCorner;
00827 else areacode |= sBoundary;
00828 }
00829
00830 if ((areacode & sBoundary) != sBoundary) {
00831 areacode |= (sAxis0 & sAxisX) | (sAxis1 & sAxisZ);
00832 }
00833 }
00834 return areacode;
00835 } else {
00836 G4Exception("G4TwistTubsSide::GetAreaCode()",
00837 "GeomSolids0001", FatalException,
00838 "Feature NOT implemented !");
00839 }
00840 return areacode;
00841 }
00842
00843
00844
00845
00846 void G4TwistTubsSide::SetCorners(
00847 G4double endInnerRad[2],
00848 G4double endOuterRad[2],
00849 G4double endPhi[2],
00850 G4double endZ[2])
00851 {
00852
00853
00854 if (fAxis[0] == kXAxis && fAxis[1] == kZAxis) {
00855
00856 G4int zmin = 0 ;
00857 G4int zmax = 1 ;
00858
00859 G4double x, y, z;
00860
00861
00862 x = endInnerRad[zmin]*std::cos(endPhi[zmin]);
00863 y = endInnerRad[zmin]*std::sin(endPhi[zmin]);
00864 z = endZ[zmin];
00865 SetCorner(sC0Min1Min, x, y, z);
00866
00867
00868 x = endOuterRad[zmin]*std::cos(endPhi[zmin]);
00869 y = endOuterRad[zmin]*std::sin(endPhi[zmin]);
00870 z = endZ[zmin];
00871 SetCorner(sC0Max1Min, x, y, z);
00872
00873
00874 x = endOuterRad[zmax]*std::cos(endPhi[zmax]);
00875 y = endOuterRad[zmax]*std::sin(endPhi[zmax]);
00876 z = endZ[zmax];
00877 SetCorner(sC0Max1Max, x, y, z);
00878
00879
00880 x = endInnerRad[zmax]*std::cos(endPhi[zmax]);
00881 y = endInnerRad[zmax]*std::sin(endPhi[zmax]);
00882 z = endZ[zmax];
00883 SetCorner(sC0Min1Max, x, y, z);
00884
00885 } else {
00886 std::ostringstream message;
00887 message << "Feature NOT implemented !" << G4endl
00888 << " fAxis[0] = " << fAxis[0] << G4endl
00889 << " fAxis[1] = " << fAxis[1];
00890 G4Exception("G4TwistTubsSide::SetCorners()",
00891 "GeomSolids0001", FatalException, message);
00892 }
00893 }
00894
00895
00896
00897
00898 void G4TwistTubsSide::SetCorners()
00899 {
00900 G4Exception("G4TwistTubsSide::SetCorners()",
00901 "GeomSolids0001", FatalException,
00902 "Method NOT implemented !");
00903 }
00904
00905
00906
00907
00908 void G4TwistTubsSide::SetBoundaries()
00909 {
00910
00911
00912 G4ThreeVector direction;
00913
00914 if (fAxis[0] == kXAxis && fAxis[1] == kZAxis) {
00915
00916
00917 direction = GetCorner(sC0Min1Max) - GetCorner(sC0Min1Min);
00918 direction = direction.unit();
00919 SetBoundary(sAxis0 & (sAxisX | sAxisMin), direction,
00920 GetCorner(sC0Min1Min), sAxisZ) ;
00921
00922
00923 direction = GetCorner(sC0Max1Max) - GetCorner(sC0Max1Min);
00924 direction = direction.unit();
00925 SetBoundary(sAxis0 & (sAxisX | sAxisMax), direction,
00926 GetCorner(sC0Max1Min), sAxisZ);
00927
00928
00929 direction = GetCorner(sC0Max1Min) - GetCorner(sC0Min1Min);
00930 direction = direction.unit();
00931 SetBoundary(sAxis1 & (sAxisZ | sAxisMin), direction,
00932 GetCorner(sC0Min1Min), sAxisX);
00933
00934
00935 direction = GetCorner(sC0Max1Max) - GetCorner(sC0Min1Max);
00936 direction = direction.unit();
00937 SetBoundary(sAxis1 & (sAxisZ | sAxisMax), direction,
00938 GetCorner(sC0Min1Max), sAxisX);
00939
00940 } else {
00941 std::ostringstream message;
00942 message << "Feature NOT implemented !" << G4endl
00943 << " fAxis[0] = " << fAxis[0] << G4endl
00944 << " fAxis[1] = " << fAxis[1];
00945 G4Exception("G4TwistTubsSide::SetCorners()",
00946 "GeomSolids0001", FatalException, message);
00947 }
00948 }
00949
00950
00951
00952
00953 void G4TwistTubsSide::GetFacets( G4int k, G4int n, G4double xyz[][3],
00954 G4int faces[][4], G4int iside )
00955 {
00956
00957 G4double z ;
00958 G4double x,xmin,xmax ;
00959
00960 G4ThreeVector p ;
00961
00962 G4int nnode ;
00963 G4int nface ;
00964
00965
00966
00967 G4int i,j ;
00968
00969 for ( i = 0 ; i<n ; i++ )
00970 {
00971
00972 z = fAxisMin[1] + i*(fAxisMax[1]-fAxisMin[1])/(n-1) ;
00973
00974 for ( j = 0 ; j<k ; j++ ) {
00975
00976 nnode = GetNode(i,j,k,n,iside) ;
00977
00978 xmin = GetBoundaryMin(z) ;
00979 xmax = GetBoundaryMax(z) ;
00980
00981 if (fHandedness < 0) {
00982 x = xmin + j*(xmax-xmin)/(k-1) ;
00983 } else {
00984 x = xmax - j*(xmax-xmin)/(k-1) ;
00985 }
00986
00987 p = SurfacePoint(x,z,true) ;
00988
00989 xyz[nnode][0] = p.x() ;
00990 xyz[nnode][1] = p.y() ;
00991 xyz[nnode][2] = p.z() ;
00992
00993 if ( i<n-1 && j<k-1 ) {
00994
00995 nface = GetFace(i,j,k,n,iside) ;
00996
00997 faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,1) * ( GetNode(i ,j ,k,n,iside)+1) ;
00998 faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,1) * ( GetNode(i+1,j ,k,n,iside)+1) ;
00999 faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,1) * ( GetNode(i+1,j+1,k,n,iside)+1) ;
01000 faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,1) * ( GetNode(i ,j+1,k,n,iside)+1) ;
01001
01002 }
01003 }
01004 }
01005 }