Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // G4TwistTubsSide implementation 27 // 28 // 01-Aug-2002 - Kotoyo Hoshina (hoshina@hepbu 29 // 13-Nov-2003 - O.Link (Oliver.Link@cern.ch), 30 // from original version in Jupi 31 // ------------------------------------------- 32 33 #include "G4TwistTubsSide.hh" 34 35 //============================================ 36 //* constructors ----------------------------- 37 38 G4TwistTubsSide::G4TwistTubsSide(const G4Strin 39 const G4Rotat 40 const G4Three 41 G4int 42 const G4doubl 43 const EAxis 44 const EAxis 45 G4doubl 46 G4doubl 47 G4doubl 48 G4doubl 49 : G4VTwistSurface(name, rot, tlate, handedn 50 axis0min, axis1min, axis0 51 fKappa(kappa) 52 { 53 if (axis0 == kZAxis && axis1 == kXAxis) 54 { 55 G4Exception("G4TwistTubsSide::G4TwistTub 56 FatalErrorInArgument, "Shoul 57 } 58 fIsValidNorm = false; 59 SetCorners(); 60 SetBoundaries(); 61 } 62 63 G4TwistTubsSide::G4TwistTubsSide(const G4Strin 64 G4doubl 65 G4doubl 66 G4doubl 67 G4doubl 68 G4doubl 69 G4doubl 70 G4doubl 71 G4doubl 72 G4int 73 : G4VTwistSurface(name) 74 { 75 fHandedness = handedness; // +z = +ve, -z 76 fAxis[0] = kXAxis; // in local coordinat 77 fAxis[1] = kZAxis; 78 fAxisMin[0] = InnerRadius; // Inner-hype r 79 fAxisMax[0] = OuterRadius; // Outer-hype r 80 fAxisMin[1] = EndZ[0]; 81 fAxisMax[1] = EndZ[1]; 82 83 fKappa = Kappa; 84 fRot.rotateZ( fHandedness > 0 85 ? -0.5*DPhi 86 : 0.5*DPhi ); 87 fTrans.set(0, 0, 0); 88 fIsValidNorm = false; 89 90 SetCorners( EndInnerRadius, EndOuterRadius, 91 SetBoundaries(); 92 } 93 94 //============================================ 95 //* Fake default constructor ----------------- 96 97 G4TwistTubsSide::G4TwistTubsSide( __void__& a 98 : G4VTwistSurface(a) 99 { 100 } 101 102 103 //============================================ 104 //* destructor ------------------------------- 105 106 G4TwistTubsSide::~G4TwistTubsSide() = default; 107 108 //============================================ 109 //* GetNormal -------------------------------- 110 111 G4ThreeVector G4TwistTubsSide::GetNormal(const 112 113 { 114 // GetNormal returns a normal vector at a s 115 // to surface) point at tmpxx. 116 // If isGlobal=true, it returns the normal 117 // 118 G4ThreeVector xx; 119 if (isGlobal) 120 { 121 xx = ComputeLocalPoint(tmpxx); 122 if ((xx - fCurrentNormal.p).mag() < 0.5 123 { 124 return ComputeGlobalDirection(fCurren 125 } 126 } 127 else 128 { 129 xx = tmpxx; 130 if (xx == fCurrentNormal.p) 131 { 132 return fCurrentNormal.normal; 133 } 134 } 135 136 G4ThreeVector er(1, fKappa * xx.z(), 0); 137 G4ThreeVector ez(0, fKappa * xx.x(), 1); 138 G4ThreeVector normal = fHandedness*(er.cros 139 140 if (isGlobal) 141 { 142 fCurrentNormal.normal = ComputeGlobalDir 143 } 144 else 145 { 146 fCurrentNormal.normal = normal.unit(); 147 } 148 return fCurrentNormal.normal; 149 } 150 151 //============================================ 152 //* DistanceToSurface ------------------------ 153 154 G4int G4TwistTubsSide::DistanceToSurface(const 155 const 156 157 158 159 160 161 { 162 // Coordinate system: 163 // 164 // The coordinate system is so chosen th 165 // the twisted surface with the z=0 plan 166 // x-axis. 167 // Rotation matrix from this coordinate 168 // to global system is saved in fRot fie 169 // So the (global) particle position and 170 // p and v, should be rotated fRot.inver 171 // to local vectors. 172 // 173 // Equation of a twisted surface: 174 // 175 // x(rho(z=0), z) = rho(z=0) 176 // y(rho(z=0), z) = rho(z=0)*K*z 177 // z(rho(z=0), z) = z 178 // with 179 // K = std::tan(fPhiTwist/2)/fZHalfLe 180 // 181 // Equation of a line: 182 // 183 // gxx = p + t*v 184 // with 185 // p = fRot.inverse()*gp 186 // v = fRot.inverse()*gv 187 // 188 // Solution for intersection: 189 // 190 // Required time for crossing is given b 191 // following quadratic equation: 192 // 193 // a*t^2 + b*t + c = 0 194 // 195 // where 196 // 197 // a = K*v_x*v_z 198 // b = K*(v_x*p_z + v_z*p_x) - v_y 199 // c = K*p_x*p_z - p_y 200 // 201 // Out of the possible two solutions you 202 // the one that gives a positive rho(z=0 203 // 204 // 205 206 fCurStatWithV.ResetfDone(validate, &gp, &gv 207 208 if (fCurStatWithV.IsDone()) 209 { 210 for (G4int i=0; i<fCurStatWithV.GetNXX() 211 { 212 gxx[i] = fCurStatWithV.GetXX(i); 213 distance[i] = fCurStatWithV.GetDistan 214 areacode[i] = fCurStatWithV.GetAreaco 215 isvalid[i] = fCurStatWithV.IsValid(i 216 } 217 return fCurStatWithV.GetNXX(); 218 } 219 else // initialize 220 { 221 for (auto i=0; i<2; ++i) 222 { 223 distance[i] = kInfinity; 224 areacode[i] = sOutside; 225 isvalid[i] = false; 226 gxx[i].set(kInfinity, kInfinity, kInf 227 } 228 } 229 230 G4ThreeVector p = ComputeLocalPoint(gp); 231 G4ThreeVector v = ComputeLocalDirection(gv) 232 G4ThreeVector xx[2]; 233 234 // 235 // special case! 236 // p is origin or 237 // 238 239 G4double absvz = std::fabs(v.z()); 240 241 if ((absvz<DBL_MIN) && (std::fabs(p.x() * v 242 { 243 // no intersection 244 245 isvalid[0] = false; 246 fCurStat.SetCurrentStatus(0, gxx[0], dis 247 isvalid[0], 0, 248 return 0; 249 } 250 251 // 252 // special case end 253 // 254 255 G4double a = fKappa * v.x() * v.z(); 256 G4double b = fKappa * (v.x() * p.z() + v.z( 257 G4double c = fKappa * p.x() * p.z() - p.y() 258 G4double D = b * b - 4 * a * c; 259 G4int vout = 0; 260 261 if (std::fabs(a) < DBL_MIN) 262 { 263 if (std::fabs(b) > DBL_MIN) 264 { 265 // single solution 266 267 distance[0] = - c / b; 268 xx[0] = p + distance[0]*v; 269 gxx[0] = ComputeGlobalPoint(xx[0 270 271 if (validate == kValidateWithTol) 272 { 273 areacode[0] = GetAreaCode(xx[0]); 274 if (!IsOutside(areacode[0])) 275 { 276 if (distance[0] >= 0) isvalid[0 277 } 278 } 279 else if (validate == kValidateWithout 280 { 281 areacode[0] = GetAreaCode(xx[0], f 282 if (IsInside(areacode[0])) 283 { 284 if (distance[0] >= 0) isvalid[0 285 } 286 } 287 else // kDontValidate 288 { 289 // we must omit x(rho,z) = rho(z=0 290 if (xx[0].x() > 0) 291 { 292 areacode[0] = sInside; 293 if (distance[0] >= 0) isvalid[0 294 } 295 else 296 { 297 distance[0] = kInfinity; 298 fCurStatWithV.SetCurrentStatus( 299 300 301 return vout; 302 } 303 } 304 305 fCurStatWithV.SetCurrentStatus(0, gxx 306 isvali 307 vout = 1; 308 } 309 else 310 { 311 // if a=b=0 , v.y=0 and (v.x=0 && p.x 312 // if v.x=0 && p.x=0, no intersect 313 // (in that case, v is paralell to 314 // if v.z=0 && p.z=0, no intersect 315 // (in that case, v is paralell to 316 // return distance = infinity. 317 318 fCurStatWithV.SetCurrentStatus(0, gxx 319 isvali 320 } 321 } 322 else if (D > DBL_MIN) 323 { 324 // double solutions 325 326 D = std::sqrt(D); 327 G4double factor = 0.5/a; 328 G4double tmpdist[2] = {kInfinity, k 329 G4ThreeVector tmpxx[2]; 330 G4int tmpareacode[2] = {sOutside 331 G4bool tmpisvalid[2] = {false, f 332 333 for (auto i=0; i<2; ++i) 334 { 335 G4double bminusD = - b - D; 336 337 // protection against round off error 338 //G4double protection = 1.0e-6; 339 G4double protection = 0; 340 if ( b * D < 0 && std::fabs(bminusD / 341 { 342 G4double acovbb = (a*c)/(b*b); 343 tmpdist[i] = - c/b * ( 1 - acovbb 344 } 345 else 346 { 347 tmpdist[i] = factor * bminusD; 348 } 349 350 D = -D; 351 tmpxx[i] = p + tmpdist[i]*v; 352 353 if (validate == kValidateWithTol) 354 { 355 tmpareacode[i] = GetAreaCode(tmpxx 356 if (!IsOutside(tmpareacode[i])) 357 { 358 if (tmpdist[i] >= 0) tmpisvalid 359 continue; 360 } 361 } 362 else if (validate == kValidateWithout 363 { 364 tmpareacode[i] = GetAreaCode(tmpxx 365 if (IsInside(tmpareacode[i])) 366 { 367 if (tmpdist[i] >= 0) tmpisvalid 368 continue; 369 } 370 } 371 else // kDontValidate 372 { 373 // we must choose x(rho,z) = rho(z 374 if (tmpxx[i].x() > 0) 375 { 376 tmpareacode[i] = sInside; 377 if (tmpdist[i] >= 0) tmpisvalid 378 continue; 379 } else { 380 tmpdist[i] = kInfinity; 381 continue; 382 } 383 } 384 } 385 386 if (tmpdist[0] <= tmpdist[1]) 387 { 388 distance[0] = tmpdist[0]; 389 distance[1] = tmpdist[1]; 390 xx[0] = tmpxx[0]; 391 xx[1] = tmpxx[1]; 392 gxx[0] = ComputeGlobalPoint(tmpx 393 gxx[1] = ComputeGlobalPoint(tmpx 394 areacode[0] = tmpareacode[0]; 395 areacode[1] = tmpareacode[1]; 396 isvalid[0] = tmpisvalid[0]; 397 isvalid[1] = tmpisvalid[1]; 398 } 399 else 400 { 401 distance[0] = tmpdist[1]; 402 distance[1] = tmpdist[0]; 403 xx[0] = tmpxx[1]; 404 xx[1] = tmpxx[0]; 405 gxx[0] = ComputeGlobalPoint(tmpx 406 gxx[1] = ComputeGlobalPoint(tmpx 407 areacode[0] = tmpareacode[1]; 408 areacode[1] = tmpareacode[0]; 409 isvalid[0] = tmpisvalid[1]; 410 isvalid[1] = tmpisvalid[0]; 411 } 412 413 fCurStatWithV.SetCurrentStatus(0, gxx[0] 414 isvalid[0 415 fCurStatWithV.SetCurrentStatus(1, gxx[1] 416 isvalid[1 417 418 // protection against roundoff error 419 420 for (G4int k=0; k<2; ++k) 421 { 422 if (!isvalid[k]) continue; 423 424 G4ThreeVector xxonsurface(xx[k].x(), 425 426 G4double deltaY = (xx[k] - xxo 427 428 if ( deltaY > 0.5*kCarTolerance ) 429 { 430 G4int maxcount = 10; 431 G4int l; 432 G4double lastdeltaY = deltaY; 433 for (l=0; l<maxcount; ++l) 434 { 435 G4ThreeVector surfacenormal = Get 436 distance[k] = DistanceToPlaneWith 437 438 deltaY = (xx[k] - xxonsurfac 439 if (deltaY > lastdeltaY) { } // 440 gxx[k] = ComputeGlobalPoint( 441 442 if (deltaY <= 0.5*kCarTolerance) 443 xxonsurface.set(xx[k].x(), 444 fKappa * std::fab 445 xx[k].z()); 446 } 447 if (l == maxcount) 448 { 449 std::ostringstream message; 450 message << "Exceeded maxloop coun 451 << " maxloop count 452 G4Exception("G4TwistTubsFlatSide: 453 "GeomSolids0003", Fa 454 } 455 } 456 } 457 vout = 2; 458 } 459 else 460 { 461 // if D<0, no solution 462 // if D=0, just grazing the surfaces, re 463 464 fCurStatWithV.SetCurrentStatus(0, gxx[0] 465 isvalid[0 466 } 467 468 return vout; 469 } 470 471 //============================================ 472 //* DistanceToSurface ------------------------ 473 474 G4int G4TwistTubsSide::DistanceToSurface(const 475 476 477 478 { 479 fCurStat.ResetfDone(kDontValidate, &gp); 480 if (fCurStat.IsDone()) 481 { 482 for (G4int i=0; i<fCurStat.GetNXX(); ++i 483 { 484 gxx[i] = fCurStat.GetXX(i); 485 distance[i] = fCurStat.GetDistance(i) 486 areacode[i] = fCurStat.GetAreacode(i) 487 } 488 return fCurStat.GetNXX(); 489 } 490 else // initialize 491 { 492 for (auto i=0; i<2; ++i) 493 { 494 distance[i] = kInfinity; 495 areacode[i] = sOutside; 496 gxx[i].set(kInfinity, kInfinity, kInf 497 } 498 } 499 500 const G4double halftol = 0.5 * kCarToleranc 501 502 G4ThreeVector p = ComputeLocalPoint( 503 G4ThreeVector xx; 504 G4int parity = (fKappa >= 0 ? 1 : 505 506 // 507 // special case! 508 // If p is on surface, or 509 // p is on z-axis, 510 // return here immediatery. 511 // 512 513 G4ThreeVector lastgxx[2]; 514 for (auto i=0; i<2; ++i) 515 { 516 lastgxx[i] = fCurStatWithV.GetXX(i); 517 } 518 519 if ((gp - lastgxx[0]).mag() < halftol 520 || (gp - lastgxx[1]).mag() < halftol) 521 { 522 // last winner, or last poststep point i 523 xx = p; 524 distance[0] = 0; 525 gxx[0] = gp; 526 527 G4bool isvalid = true; 528 fCurStat.SetCurrentStatus(0, gxx[0], dis 529 isvalid, 1, kDont 530 return 1; 531 } 532 533 if (p.getRho() == 0) 534 { 535 // p is on z-axis. Namely, p is on twist 536 // We must return here, however, returni 537 // boundary is better than return 0-dist 538 // 539 G4bool isvalid = true; 540 if (fAxis[0] == kXAxis && fAxis[1] == kZ 541 { 542 distance[0] = DistanceToBoundary(sAxi 543 areacode[0] = sInside; 544 } 545 else 546 { 547 distance[0] = 0; 548 xx.set(0., 0., 0.); 549 } 550 gxx[0] = ComputeGlobalPoint(xx); 551 fCurStat.SetCurrentStatus(0, gxx[0], dis 552 isvalid, 0, kD 553 return 1; 554 } 555 556 // 557 // special case end 558 // 559 560 // set corner points of quadrangle try area 561 562 G4ThreeVector A; // foot of normal from p 563 G4ThreeVector C; // foot of normal from p 564 G4ThreeVector B; // point on boundary 565 G4ThreeVector D; // point on boundary 566 567 // G4double distToA; // distance from 568 DistanceToBoundary(sAxis0 & sAxisMin, A, p) 569 // G4double distToC; // distance from 570 DistanceToBoundary(sAxis0 & sAxisMax, C, p) 571 572 // is p.z between a.z and c.z? 573 // p.z must be bracketed a.z and c.z. 574 if (A.z() > C.z()) 575 { 576 if (p.z() > A.z()) 577 { 578 A = GetBoundaryAtPZ(sAxis0 & sAxisMin 579 } 580 else if (p.z() < C.z()) 581 { 582 C = GetBoundaryAtPZ(sAxis0 & sAxisMax 583 } 584 } 585 else 586 { 587 if (p.z() > C.z()) 588 { 589 C = GetBoundaryAtPZ(sAxis0 & sAxisMax 590 } 591 else if (p.z() < A.z()) 592 { 593 A = GetBoundaryAtPZ(sAxis0 & sAxisMin 594 } 595 } 596 597 G4ThreeVector d[2]; // direction vecto 598 G4ThreeVector x0[2]; // foot of normal 599 G4int btype[2]; // boundary type 600 601 for (auto i=0; i<2; ++i) 602 { 603 if (i == 0) 604 { 605 GetBoundaryParameters((sAxis0 & sAxis 606 B = x0[i] + ((A.z() - x0[i].z()) / d[ 607 // x0 + t*d , d is direction unit vec 608 } 609 else 610 { 611 GetBoundaryParameters((sAxis0 & sAxis 612 D = x0[i] + ((C.z() - x0[i].z()) / d[ 613 } 614 } 615 616 // In order to set correct diagonal, swap A 617 G4ThreeVector pt(p.x(), p.y(), 0.); 618 G4double rc = std::fabs(p.x()); 619 G4ThreeVector surfacevector(rc, rc * fKappa 620 G4int pside = AmIOnLeftSide(pt, sur 621 G4double test = (A.z() - C.z()) * par 622 623 if (test == 0) 624 { 625 if (pside == 0) 626 { 627 // p is on surface. 628 xx = p; 629 distance[0] = 0; 630 gxx[0] = gp; 631 632 G4bool isvalid = true; 633 fCurStat.SetCurrentStatus(0, gxx[0], 634 isvalid, 1, 635 return 1; 636 } 637 else 638 { 639 // A.z = C.z(). return distance to li 640 d[0] = C - A; 641 distance[0] = DistanceToLine(p, A, d[ 642 areacode[0] = sInside; 643 gxx[0] = ComputeGlobalPoint(xx); 644 G4bool isvalid = true; 645 fCurStat.SetCurrentStatus(0, gxx[0], 646 isvalid, 1, 647 return 1; 648 } 649 } 650 else if (test < 0) // wrong diagonal. vect 651 { // swap A and D, C and 652 G4ThreeVector tmp; 653 tmp = A; 654 A = D; 655 D = tmp; 656 tmp = C; 657 C = B; 658 B = tmp; 659 660 } 661 else // correct diagonal. nothing to do. 662 { 663 } 664 665 // Now, we chose correct diagonal. 666 // First try. divide quadrangle into double 667 // calculate distance to both surfaces. 668 669 G4ThreeVector xxacb; // foot of normal fr 670 G4ThreeVector nacb; // normal of plane A 671 G4ThreeVector xxcad; // foot of normal fr 672 G4ThreeVector ncad; // normal of plane C 673 G4ThreeVector AB(A.x(), A.y(), 0); 674 G4ThreeVector DC(C.x(), C.y(), 0); 675 676 G4double distToACB = G4VTwistSurface::Dista 677 678 G4double distToCAD = G4VTwistSurface::Dista 679 680 // if calculated distance = 0, return 681 682 if (std::fabs(distToACB) <= halftol || std: 683 { 684 xx = (std::fabs(distToACB) < std::fabs(d 685 areacode[0] = sInside; 686 gxx[0] = ComputeGlobalPoint(xx); 687 distance[0] = 0; 688 G4bool isvalid = true; 689 fCurStat.SetCurrentStatus(0, gxx[0], dis 690 isvalid, 1, kD 691 return 1; 692 } 693 694 if (distToACB * distToCAD > 0 && distToACB 695 { 696 // both distToACB and distToCAD are nega 697 // divide quadrangle into double triangl 698 G4ThreeVector normal; 699 distance[0] = DistanceToPlane(p, A, B, C 700 } 701 else 702 { 703 if (distToACB * distToCAD > 0) 704 { 705 // both distToACB and distToCAD are p 706 // Take smaller one. 707 if (distToACB <= distToCAD) 708 { 709 distance[0] = distToACB; 710 xx = xxacb; 711 } 712 else 713 { 714 distance[0] = distToCAD; 715 xx = xxcad; 716 } 717 } 718 else 719 { 720 // distToACB * distToCAD is negative. 721 // take positive one 722 if (distToACB > 0) 723 { 724 distance[0] = distToACB; 725 xx = xxacb; 726 } 727 else 728 { 729 distance[0] = distToCAD; 730 xx = xxcad; 731 } 732 } 733 } 734 areacode[0] = sInside; 735 gxx[0] = ComputeGlobalPoint(xx); 736 G4bool isvalid = true; 737 fCurStat.SetCurrentStatus(0, gxx[0], distan 738 isvalid, 1, kDont 739 return 1; 740 } 741 742 //============================================ 743 //* DistanceToPlane -------------------------- 744 745 G4double G4TwistTubsSide::DistanceToPlane(cons 746 cons 747 cons 748 cons 749 cons 750 cons 751 752 753 { 754 const G4double halftol = 0.5 * kCarToleranc 755 756 G4ThreeVector M = 0.5*(A + B); 757 G4ThreeVector N = 0.5*(C + D); 758 G4ThreeVector xxanm; // foot of normal fro 759 G4ThreeVector nanm; // normal of plane AN 760 G4ThreeVector xxcmn; // foot of normal fro 761 G4ThreeVector ncmn; // normal of plane CM 762 763 G4double distToanm = G4VTwistSurface::Dista 764 765 G4double distTocmn = G4VTwistSurface::Dista 766 767 #ifdef G4SPECSDEBUG 768 // if p is behind of both surfaces, abort. 769 if (distToanm * distTocmn > 0 && distToanm 770 { 771 G4Exception("G4TwistTubsSide::DistanceToP 772 "GeomSolids0003", FatalExcept 773 "Point p is behind the surfac 774 } 775 #endif 776 // if p is on surface, return 0. 777 if (std::fabs(distToanm) <= halftol) 778 { 779 xx = xxanm; 780 n = nanm * parity; 781 return 0; 782 } 783 else if (std::fabs(distTocmn) <= halftol) 784 { 785 xx = xxcmn; 786 n = ncmn * parity; 787 return 0; 788 } 789 790 if (distToanm <= distTocmn) 791 { 792 if (distToanm > 0) 793 { 794 // both distanses are positive. take 795 xx = xxanm; 796 n = nanm * parity; 797 return distToanm; 798 } 799 else 800 { 801 // take -ve distance and call the fun 802 return DistanceToPlane(p, A, M, N, D, 803 } 804 } 805 else 806 { 807 if (distTocmn > 0) 808 { 809 // both distanses are positive. take 810 xx = xxcmn; 811 n = ncmn * parity; 812 return distTocmn; 813 } 814 else 815 { 816 // take -ve distance and call the fun 817 return DistanceToPlane(p, C, N, M, B, 818 } 819 } 820 } 821 822 //============================================ 823 //* GetAreaCode ------------------------------ 824 825 G4int G4TwistTubsSide::GetAreaCode(const G4Thr 826 G4boo 827 { 828 // We must use the function in local coordi 829 // See the description of DistanceToSurface 830 831 const G4double ctol = 0.5 * kCarTolerance; 832 G4int areacode = sInside; 833 834 if (fAxis[0] == kXAxis && fAxis[1] == kZAxi 835 { 836 G4int xaxis = 0; 837 G4int zaxis = 1; 838 839 if (withTol) 840 { 841 G4bool isoutside = false; 842 843 // test boundary of xaxis 844 845 if (xx.x() < fAxisMin[xaxis] + ctol) 846 { 847 areacode |= (sAxis0 & (sAxisX | sA 848 if (xx.x() <= fAxisMin[xaxis] - ct 849 850 } 851 else if (xx.x() > fAxisMax[xaxis] - c 852 { 853 areacode |= (sAxis0 & (sAxisX | sA 854 if (xx.x() >= fAxisMax[xaxis] + ct 855 } 856 857 // test boundary of z-axis 858 859 if (xx.z() < fAxisMin[zaxis] + ctol) 860 { 861 areacode |= (sAxis1 & (sAxisZ | sA 862 863 if ((areacode & sBoundary) != 0) 864 else areaco 865 if (xx.z() <= fAxisMin[zaxis] - ct 866 867 } 868 else if (xx.z() > fAxisMax[zaxis] - c 869 { 870 areacode |= (sAxis1 & (sAxisZ | sA 871 872 if ((areacode & sBoundary) != 0) 873 else areaco 874 if (xx.z() >= fAxisMax[zaxis] + ct 875 } 876 877 // if isoutside = true, clear inside 878 // if not on boundary, add axis infor 879 880 if (isoutside) 881 { 882 G4int tmpareacode = areacode & (~s 883 areacode = tmpareacode; 884 } 885 else if ((areacode & sBoundary) != sB 886 { 887 areacode |= (sAxis0 & sAxisX) | (s 888 } 889 } 890 else 891 { 892 // boundary of x-axis 893 894 if (xx.x() < fAxisMin[xaxis] ) 895 { 896 areacode |= (sAxis0 & (sAxisX | sA 897 } 898 else if (xx.x() > fAxisMax[xaxis]) 899 { 900 areacode |= (sAxis0 & (sAxisX | sA 901 } 902 903 // boundary of z-axis 904 905 if (xx.z() < fAxisMin[zaxis]) 906 { 907 areacode |= (sAxis1 & (sAxisZ | sA 908 if ((areacode & sBoundary) != 0) 909 else areaco 910 911 } 912 else if (xx.z() > fAxisMax[zaxis]) 913 { 914 areacode |= (sAxis1 & (sAxisZ | sA 915 if ((areacode & sBoundary) != 0) 916 else areaco 917 } 918 919 if ((areacode & sBoundary) != sBounda 920 { 921 areacode |= (sAxis0 & sAxisX) | (s 922 } 923 } 924 return areacode; 925 } 926 else 927 { 928 G4Exception("G4TwistTubsSide::GetAreaCod 929 "GeomSolids0001", FatalExcep 930 "Feature NOT implemented !") 931 } 932 return areacode; 933 } 934 935 //============================================ 936 //* SetCorners( arglist ) -------------------- 937 938 void G4TwistTubsSide::SetCorners( G4double end 939 G4double end 940 G4double end 941 G4double end 942 { 943 // Set Corner points in local coodinate. 944 945 if (fAxis[0] == kXAxis && fAxis[1] == kZAxi 946 { 947 G4int zmin = 0 ; // at -ve z 948 G4int zmax = 1 ; // at +ve z 949 950 G4double x, y, z; 951 952 // corner of Axis0min and Axis1min 953 x = endInnerRad[zmin]*std::cos(endPhi[zm 954 y = endInnerRad[zmin]*std::sin(endPhi[zm 955 z = endZ[zmin]; 956 SetCorner(sC0Min1Min, x, y, z); 957 958 // corner of Axis0max and Axis1min 959 x = endOuterRad[zmin]*std::cos(endPhi[zm 960 y = endOuterRad[zmin]*std::sin(endPhi[zm 961 z = endZ[zmin]; 962 SetCorner(sC0Max1Min, x, y, z); 963 964 // corner of Axis0max and Axis1max 965 x = endOuterRad[zmax]*std::cos(endPhi[zm 966 y = endOuterRad[zmax]*std::sin(endPhi[zm 967 z = endZ[zmax]; 968 SetCorner(sC0Max1Max, x, y, z); 969 970 // corner of Axis0min and Axis1max 971 x = endInnerRad[zmax]*std::cos(endPhi[zm 972 y = endInnerRad[zmax]*std::sin(endPhi[zm 973 z = endZ[zmax]; 974 SetCorner(sC0Min1Max, x, y, z); 975 976 } 977 else 978 { 979 std::ostringstream message; 980 message << "Feature NOT implemented !" < 981 << " fAxis[0] = " << fAxi 982 << " fAxis[1] = " << fAxi 983 G4Exception("G4TwistTubsSide::SetCorners 984 "GeomSolids0001", FatalExcep 985 } 986 } 987 988 //============================================ 989 //* SetCorners() ----------------------------- 990 991 void G4TwistTubsSide::SetCorners() 992 { 993 G4Exception("G4TwistTubsSide::SetCorners()" 994 "GeomSolids0001", FatalExceptio 995 "Method NOT implemented !"); 996 } 997 998 //============================================ 999 //* SetBoundaries() -------------------------- 1000 1001 void G4TwistTubsSide::SetBoundaries() 1002 { 1003 // Set direction-unit vector of boundary-l 1004 // 1005 G4ThreeVector direction; 1006 1007 if (fAxis[0] == kXAxis && fAxis[1] == kZAx 1008 { 1009 // sAxis0 & sAxisMin 1010 direction = GetCorner(sC0Min1Max) - Get 1011 direction = direction.unit(); 1012 SetBoundary(sAxis0 & (sAxisX | sAxisMin 1013 GetCorner(sC0Min1Min), sAxi 1014 1015 // sAxis0 & sAxisMax 1016 direction = GetCorner(sC0Max1Max) - Get 1017 direction = direction.unit(); 1018 SetBoundary(sAxis0 & (sAxisX | sAxisMax 1019 GetCorner(sC0Max1Min), sAxi 1020 1021 // sAxis1 & sAxisMin 1022 direction = GetCorner(sC0Max1Min) - Get 1023 direction = direction.unit(); 1024 SetBoundary(sAxis1 & (sAxisZ | sAxisMin 1025 GetCorner(sC0Min1Min), sAxi 1026 1027 // sAxis1 & sAxisMax 1028 direction = GetCorner(sC0Max1Max) - Get 1029 direction = direction.unit(); 1030 SetBoundary(sAxis1 & (sAxisZ | sAxisMax 1031 GetCorner(sC0Min1Max), sAxi 1032 1033 } 1034 else 1035 { 1036 std::ostringstream message; 1037 message << "Feature NOT implemented !" 1038 << " fAxis[0] = " << fAx 1039 << " fAxis[1] = " << fAx 1040 G4Exception("G4TwistTubsSide::SetCorner 1041 "GeomSolids0001", FatalExce 1042 } 1043 } 1044 1045 //=========================================== 1046 //* GetFacets() ----------------------------- 1047 1048 void G4TwistTubsSide::GetFacets( G4int k, G4i 1049 G4int faces[ 1050 { 1051 G4double z ; // the two parameters for 1052 G4double x,xmin,xmax ; 1053 1054 G4ThreeVector p ; // a point on the surfac 1055 1056 G4int nnode ; 1057 G4int nface ; 1058 1059 // calculate the (n-1)*(k-1) vertices 1060 1061 for ( G4int i = 0 ; i<n ; ++i ) 1062 { 1063 z = fAxisMin[1] + i*(fAxisMax[1]-fAxisMin 1064 1065 for ( G4int j = 0 ; j<k ; ++j ) 1066 { 1067 nnode = GetNode(i,j,k,n,iside) ; 1068 1069 xmin = GetBoundaryMin(z) ; 1070 xmax = GetBoundaryMax(z) ; 1071 1072 if (fHandedness < 0) 1073 { 1074 x = xmin + j*(xmax-xmin)/(k-1) ; 1075 } 1076 else 1077 { 1078 x = xmax - j*(xmax-xmin)/(k-1) ; 1079 } 1080 1081 p = SurfacePoint(x,z,true) ; // surfac 1082 1083 xyz[nnode][0] = p.x() ; 1084 xyz[nnode][1] = p.y() ; 1085 xyz[nnode][2] = p.z() ; 1086 1087 if ( i<n-1 && j<k-1 ) // clock wise f 1088 { 1089 nface = GetFace(i,j,k,n,iside) ; 1090 1091 faces[nface][0] = GetEdgeVisibility(i,j,k,n 1092 * ( GetNode(i ,j ,k 1093 faces[nface][1] = GetEdgeVisibility(i,j,k,n 1094 * ( GetNode(i+1,j ,k 1095 faces[nface][2] = GetEdgeVisibility(i,j,k,n 1096 * ( GetNode(i+1,j+1,k 1097 faces[nface][3] = GetEdgeVisibility(i,j,k,n 1098 * ( GetNode(i ,j+1,k 1099 } 1100 } 1101 } 1102 } 1103