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