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 // G4GenericTrap implementation 27 // 28 // Authors: 29 // Tatiana Nikitina, CERN; Ivana Hrivnacova, 30 // Adapted from Root Arb8 implementation by 31 // 32 // 27.05.2024 - Evgueni Tcherniaev, complete 33 // ------------------------------------------- 34 35 #include "G4GenericTrap.hh" 36 37 #if !defined(G4GEOM_USE_UGENERICTRAP) 38 39 #include <iomanip> 40 41 #include "G4PhysicalConstants.hh" 42 #include "G4SystemOfUnits.hh" 43 #include "G4VoxelLimits.hh" 44 #include "G4AffineTransform.hh" 45 #include "G4BoundingEnvelope.hh" 46 #include "G4QuickRand.hh" 47 48 #include "G4GeomTools.hh" 49 50 #include "G4VGraphicsScene.hh" 51 #include "G4Polyhedron.hh" 52 #include "G4VisExtent.hh" 53 54 #include "G4AutoLock.hh" 55 56 namespace 57 { 58 G4Mutex polyhedronMutex = G4MUTEX_INITIALIZ 59 G4Mutex lateralareaMutex = G4MUTEX_INITIALIZ 60 } 61 62 ////////////////////////////////////////////// 63 // 64 // Constructor 65 // 66 G4GenericTrap::G4GenericTrap(const G4String& n 67 const std::vector 68 : G4VSolid(name) 69 { 70 halfTolerance = 0.5*kCarTolerance; 71 CheckParameters(halfZ, vertices); 72 ComputeLateralSurfaces(); 73 ComputeBoundingBox(); 74 ComputeScratchLength(); 75 } 76 77 ////////////////////////////////////////////// 78 // 79 // Fake default constructor - sets only member 80 // for usage restri 81 G4GenericTrap::G4GenericTrap(__void__& a) 82 : G4VSolid(a) 83 { 84 } 85 86 ////////////////////////////////////////////// 87 // 88 // Destructor 89 90 G4GenericTrap::~G4GenericTrap() 91 { 92 } 93 94 ////////////////////////////////////////////// 95 // 96 // Copy constructor 97 // 98 G4GenericTrap::G4GenericTrap(const G4GenericTr 99 : G4VSolid(rhs), 100 halfTolerance(rhs.halfTolerance), fScratch 101 fDz(rhs.fDz), fVertices(rhs.fVertices), fI 102 fMinBBox(rhs.fMinBBox), fMaxBBox(rhs.fMaxB 103 fVisSubdivisions(rhs.fVisSubdivisions), 104 fSurfaceArea(rhs.fSurfaceArea), fCubicVolu 105 { 106 for (auto i = 0; i < 5; ++i) { fTwist[i] = r 107 for (auto i = 0; i < 4; ++i) { fDelta[i] = r 108 for (auto i = 0; i < 8; ++i) { fPlane[i] = r 109 for (auto i = 0; i < 4; ++i) { fSurf[i] = rh 110 for (auto i = 0; i < 4; ++i) { fArea[i] = rh 111 } 112 113 ////////////////////////////////////////////// 114 // 115 // Assignment 116 // 117 G4GenericTrap& G4GenericTrap::operator=(const 118 { 119 // Check assignment to self 120 if (this == &rhs) { return *this; } 121 122 // Copy base class data 123 G4VSolid::operator=(rhs); 124 125 // Copy data 126 halfTolerance = rhs.halfTolerance; fScratch 127 fDz = rhs.fDz; fVertices = rhs.fVertices; fI 128 fMinBBox = rhs.fMinBBox; fMaxBBox = rhs.fMax 129 fVisSubdivisions = rhs.fVisSubdivisions; 130 fSurfaceArea = rhs.fSurfaceArea; fCubicVolum 131 132 for (auto i = 0; i < 8; ++i) { fVertices[i] 133 for (auto i = 0; i < 5; ++i) { fTwist[i] = r 134 for (auto i = 0; i < 4; ++i) { fDelta[i] = r 135 for (auto i = 0; i < 8; ++i) { fPlane[i] = r 136 for (auto i = 0; i < 4; ++i) { fSurf[i] = rh 137 for (auto i = 0; i < 4; ++i) { fArea[i] = rh 138 139 fRebuildPolyhedron = false; 140 delete fpPolyhedron; fpPolyhedron = nullptr; 141 142 return *this; 143 } 144 145 ////////////////////////////////////////////// 146 // 147 // Returns position of the point (inside/outsi 148 // 149 EInside G4GenericTrap::Inside(const G4ThreeVec 150 { 151 G4double px = p.x(), py = p.y(), pz = p.z(); 152 153 // intersect edges by z = pz plane 154 G4TwoVector pp[4]; 155 G4double z = (pz + fDz); 156 for (auto i = 0; i < 4; ++i) { pp[i] = fVert 157 158 // estimate distance to the solid 159 G4double dist = std::abs(pz) - fDz; 160 for (auto i = 0; i < 4; ++i) 161 { 162 if (fTwist[i] == 0.) 163 { 164 G4double dd = fSurf[i].D*px + fSurf[i].E 165 dist = std::max(dd, dist); 166 } 167 else 168 { 169 auto j = (i + 1)%4; 170 G4TwoVector a = pp[i]; 171 G4TwoVector b = pp[j]; 172 G4double dx = b.x() - a.x(); 173 G4double dy = b.y() - a.y(); 174 G4double leng = std::sqrt(dx*dx + dy*dy) 175 G4double dd = (dx*(py - a.y()) - dy*(px 176 dist = std::max(dd, dist); 177 } 178 } 179 return (dist > halfTolerance) ? kOutside : 180 ((dist > -halfTolerance) ? kSurface : kIns 181 } 182 183 ////////////////////////////////////////////// 184 // 185 // Return unit normal to surface at p 186 // 187 G4ThreeVector G4GenericTrap::SurfaceNormal( co 188 { 189 G4double halfToleranceSquared = halfToleranc 190 G4double px = p.x(), py = p.y(), pz = p.z(); 191 G4double nx = 0, ny = 0, nz = 0; 192 193 // intersect edges by z = pz plane 194 G4TwoVector pp[4]; 195 G4double tz = (pz + fDz); 196 for (auto i = 0; i < 4; ++i) { pp[i] = fVert 197 198 // check bottom and top faces 199 G4double dz = std::abs(pz) - fDz; 200 nz = std::copysign(G4double(std::abs(dz) <= 201 202 // check lateral faces 203 for (auto i = 0; i < 4; ++i) 204 { 205 if (fTwist[i] == 0.) 206 { 207 G4double dd = fSurf[i].D*px + fSurf[i].E 208 if (std::abs(dd) <= halfTolerance) 209 { 210 nx += fSurf[i].D; 211 ny += fSurf[i].E; 212 nz += fSurf[i].F; 213 } 214 } 215 else 216 { 217 auto j = (i + 1)%4; 218 G4TwoVector a = pp[i]; 219 G4TwoVector b = pp[j]; 220 G4double dx = b.x() - a.x(); 221 G4double dy = b.y() - a.y(); 222 G4double ll = dx*dx + dy*dy; 223 G4double dd = dx*(py - a.y()) - dy*(px - 224 if (dd*dd <= halfToleranceSquared*ll) 225 { 226 G4double x = fSurf[i].A*pz + fSurf[i]. 227 G4double y = fSurf[i].B*pz + fSurf[i]. 228 G4double z = fSurf[i].A*px + fSurf[i]. 229 G4double inv = 1./std::sqrt(x*x + y*y 230 nx += x*inv; 231 ny += y*inv; 232 nz += z*inv; 233 } 234 } 235 } 236 237 // return normal 238 G4double mag2 = nx*nx + ny*ny + nz*nz; 239 if (mag2 == 0.) return ApproxSurfaceNormal(p 240 G4double mag = std::sqrt(mag2); 241 G4double inv = (mag == 1.) ? 1. : 1./mag; 242 return { nx*inv, ny*inv, nz*inv }; 243 } 244 245 ////////////////////////////////////////////// 246 // 247 // Find surface nearest to the point and retur 248 // Normally this method should not be called 249 // 250 G4ThreeVector G4GenericTrap::ApproxSurfaceNorm 251 { 252 #ifdef G4SPECSDEBUG 253 std::ostringstream message; 254 message.precision(16); 255 message << "Point p is not on surface of sol 256 << "Position: " << p << " is " 257 << ((Inside(p) == kInside) ? "inside 258 StreamInfo(message); 259 G4Exception("G4GenericTrap::ApproxSurfaceNor 260 JustWarning, message ); 261 #endif 262 G4double px = p.x(), py = p.y(), pz = p.z(); 263 G4double dist = -DBL_MAX; 264 auto iside = 0; 265 266 // intersect edges by z = pz plane 267 G4TwoVector pp[4]; 268 G4double tz = (pz + fDz); 269 for (auto i = 0; i < 4; ++i) { pp[i] = fVert 270 271 // check lateral faces 272 for (auto i = 0; i < 4; ++i) 273 { 274 if (fTwist[i] == 0.) 275 { 276 G4double d = fSurf[i].D*px + fSurf[i].E* 277 if (d > dist) { dist = d; iside = i; } 278 } 279 else 280 { 281 auto j = (i + 1)%4; 282 G4TwoVector a = pp[i]; 283 G4TwoVector b = pp[j]; 284 G4double dx = b.x() - a.x(); 285 G4double dy = b.y() - a.y(); 286 G4double leng = std::sqrt(dx*dx + dy*dy) 287 G4double d = (dx*(py - a.y()) - dy*(px - 288 if (d > dist) { dist = d; iside = i; } 289 } 290 } 291 // check bottom and top faces 292 G4double distz = std::abs(pz) - fDz; 293 if (distz >= dist) return { 0., 0., std::cop 294 295 G4double x = fSurf[iside].A*pz + fSurf[iside 296 G4double y = fSurf[iside].B*pz + fSurf[iside 297 G4double z = fSurf[iside].A*px + fSurf[iside 298 G4double inv = 1./std::sqrt(x*x + y*y + z*z) 299 return { x*inv, y*inv, z*inv }; 300 } 301 302 ////////////////////////////////////////////// 303 // 304 // Compute distance to the surface from outsid 305 // return kInfinity if no intersection 306 // 307 G4double G4GenericTrap::DistanceToIn(const G4T 308 const G4T 309 { 310 G4double px = p.x(), py = p.y(), pz = p.z(); 311 G4double vx = v.x(), vy = v.y(), vz = v.z(); 312 313 // Find intersections with the bounding poly 314 // 315 if (std::abs(pz) - fDz >= -halfTolerance && 316 G4double invz = (vz == 0) ? kInfinity : -1./ 317 G4double dz = std::copysign(fDz,invz); 318 G4double xin = (pz - dz)*invz; 319 G4double xout = (pz + dz)*invz; 320 321 // Check plane faces 322 for (auto k = 0; k < 4; ++k) 323 { 324 if (fTwist[k] != 0) continue; // skip twis 325 const G4GenericTrapPlane& surf = fPlane[2* 326 G4double cosa = surf.A*vx + surf.B*vy + su 327 G4double dist = surf.A*px + surf.B*py + su 328 if (dist >= -halfTolerance) 329 { 330 if (cosa >= 0.) { return kInfinity; } // 331 G4double tmp = -dist/cosa; 332 xin = std::max(tmp, xin); 333 } 334 else 335 { 336 if (cosa > 0) { xout = std::min(-dist/co 337 if (cosa < 0) { xin = std::max(-dist/cos 338 } 339 } 340 341 // Check planes around twisted faces, each t 342 G4double tin = xin; 343 G4double tout = xout; 344 for (auto i = 0; i < 4; ++i) 345 { 346 if (fTwist[i] == 0) continue; // skip plan 347 348 // check intersection with 1st bounding pl 349 const G4GenericTrapPlane& surf1 = fPlane[2 350 G4double cosa = surf1.A*vx + surf1.B*vy + 351 G4double dist = surf1.A*px + surf1.B*py + 352 if (dist >= -halfTolerance) 353 { 354 if (cosa >= 0.) { return kInfinity; } // 355 G4double tmp = -dist/cosa; 356 tin = std::max(tmp, tin); 357 } 358 else 359 { 360 if (cosa > 0) { tout = std::min(-dist/co 361 if (cosa < 0) { tin = std::max(-dist/cos 362 } 363 364 // check intersection with 2nd bounding pl 365 const G4GenericTrapPlane& surf2 = fPlane[2 366 cosa = surf2.A*vx + surf2.B*vy + surf2.C*v 367 dist = surf2.A*px + surf2.B*py + surf2.C*p 368 if (dist >= -halfTolerance) 369 { 370 if (cosa >= 0.) { return kInfinity; } // 371 G4double tmp = -dist/cosa; 372 tin = std::max(tmp, tin); 373 } 374 else 375 { 376 if (cosa > 0) { tout = std::min(-dist/co 377 if (cosa < 0) { tin = std::max(-dist/cos 378 } 379 } 380 if (tout - tin <= halfTolerance) { return kI 381 382 // if point is outside of the bounding box 383 // then move it to the surface of the boundi 384 G4double t0 = 0., x0 = px, y0 = py, z0 = pz; 385 if (x0 < fMinBBox.x() - halfTolerance || 386 y0 < fMinBBox.y() - halfTolerance || 387 z0 < fMinBBox.z() - halfTolerance || 388 x0 > fMaxBBox.x() + halfTolerance || 389 y0 > fMaxBBox.y() + halfTolerance || 390 z0 > fMaxBBox.z() + halfTolerance) 391 { 392 t0 = tin; 393 x0 += vx*t0; 394 y0 += vy*t0; 395 z0 += vz*t0; 396 } 397 398 // Check intersections with twisted faces 399 // 400 G4double ttin[2] = { DBL_MAX, DBL_MAX }; 401 G4double ttout[2] = { tout, tout }; 402 403 if (tin - xin < halfTolerance) ttin[0] = xin 404 if (xout - tout < halfTolerance) ttout[0] = 405 G4double tminimal = tin - halfTolerance; 406 G4double tmaximal = tout + halfTolerance; 407 408 constexpr G4double EPSILON = 1000.*DBL_EPSIL 409 for (auto i = 0; i < 4; ++i) 410 { 411 if (fTwist[i] == 0) continue; // skip plan 412 413 // twisted face, solve quadratic equation 414 G4double ABC = fSurf[i].A*vx + fSurf[i].B 415 G4double ABCF = fSurf[i].A*x0 + fSurf[i].B 416 G4double A = ABC*vz; 417 G4double B = 0.5*(fSurf[i].D*vx + fSurf[i] 418 G4double C = fSurf[i].D*x0 + fSurf[i].E*y0 419 if (std::abs(A) <= EPSILON) 420 { 421 // case 1: track is parallel to the surf 422 if (std::abs(B) <= EPSILON) 423 { 424 // check position of the track relativ 425 // for the reason of precision it's be 426 auto j = (i + 1)%4; 427 G4double z = (z0 + fDz); 428 G4TwoVector a = fVertices[i] + fDelta[ 429 G4TwoVector b = fVertices[j] + fDelta[ 430 G4double dx = b.x() - a.x(); 431 G4double dy = b.y() - a.y(); 432 G4double leng = std::sqrt(dx*dx + dy*d 433 G4double dist = dx*(y0 - a.y()) - dy*( 434 if (dist >= -halfTolerance*leng) { ret 435 continue; 436 } 437 438 // case 2: single root 439 G4double tmp = t0 - 0.5*C/B; 440 // compute normal at the intersection po 441 G4double x = px + vx*tmp; 442 G4double y = py + vy*tmp; 443 G4double z = pz + vz*tmp; 444 const G4GenericTrapSurface& surf = fSurf 445 G4double nx = surf.A*z + surf.D; 446 G4double ny = surf.B*z + surf.E; 447 G4double nz = surf.A*x + surf.B*y + 2.*s 448 449 if (nx*vx + ny*vy + nz*vz >= 0.) // poin 450 { 451 auto k = (i == 3) ? 0 : i + 1; 452 G4double tz = (pz + fDz); 453 G4TwoVector a = fVertices[i] + fDelta[ 454 G4TwoVector b = fVertices[k] + fDelta[ 455 G4double dx = b.x() - a.x(); 456 G4double dy = b.y() - a.y(); 457 G4double leng = std::sqrt(dx*dx + dy*d 458 G4double dist = dx*(py - a.y()) - dy*( 459 if (dist >= -halfTolerance*leng) { ret 460 461 if (tmp < tminimal || tmp > tmaximal) 462 if (std::abs(tmp - ttout[0]) < halfTol 463 if (tmp < ttout[0]) 464 { 465 ttout[1] = ttout[0]; 466 ttout[0] = tmp; 467 } 468 else { ttout[1] = std::min(tmp, ttout[ 469 } 470 else // point is flying to inside 471 { 472 if (tmp < tminimal || tmp > tmaximal) 473 if (std::abs(tmp - ttin[0]) < halfTole 474 if (tmp < ttin[0]) 475 { 476 ttin[1] = ttin[0]; 477 ttin[0] = tmp; 478 } 479 else { ttin[1] = std::min(tmp, ttin[1] 480 } 481 continue; 482 } 483 484 // case 3: scratch or no intersection 485 G4double D = B*B - A*C; 486 if (D < 0.25*fScratch*fScratch*A*A) 487 { 488 if (A > 0) return kInfinity; 489 continue; 490 } 491 492 // case 4: two intersection points 493 G4double tmp = -B - std::copysign(std::sqr 494 G4double t1 = tmp/A + t0; 495 G4double t2 = C/tmp + t0; 496 G4double tsurfin = std::min(t1, t2); 497 G4double tsurfout = std::max(t1, t2); 498 if (A < 0) std::swap(tsurfin, tsurfout); 499 500 if (tsurfin >= tminimal && tsurfin <= tmax 501 { 502 if (std::abs(tsurfin - ttin[0]) >= halfT 503 { 504 if (tsurfin < ttin[0]) 505 { 506 ttin[1] = ttin[0]; 507 ttin[0] = tsurfin; 508 } 509 else { ttin[1] = std::min(tsurfin, tti 510 } 511 } 512 if (tsurfout >= tminimal && tsurfout <= tm 513 { 514 if (std::abs(tsurfout - ttout[0]) >= hal 515 { 516 if (tsurfout < ttout[0]) 517 { 518 ttout[1] = ttout[0]; 519 ttout[0] = tsurfout; 520 } 521 else { ttout[1] = std::min(tsurfout, t 522 } 523 } 524 } 525 526 // Compute distance to In 527 // 528 if (ttin[0] == DBL_MAX) { return kInfinity; 529 530 // single entry point 531 if (ttin[1] == DBL_MAX) 532 { 533 G4double distin = ttin[0]; 534 G4double distout = (distin >= ttout[0] - h 535 G4double dist = (distout <= halfTolerance 536 return (dist < halfTolerance) ? 0. : dist; 537 } 538 539 // two entry points 540 if (ttin[1] < ttout[0]) 541 { 542 G4double distin = ttin[1], distout = ttout 543 G4double dist = (distout <= halfTolerance 544 return (dist < halfTolerance) ? 0. : dist; 545 } 546 547 // check 1st pair of in-out points 548 G4double distin1 = ttin[0], distout1 = ttout 549 G4double dist1 = (distout1 <= halfTolerance 550 if (dist1 != kInfinity) { return (dist1 < ha 551 552 // check 2nd pair of in-out points 553 G4double distin2 = ttin[1], distout2 = ttout 554 G4double dist2 = (distout2 <= halfTolerance 555 return (dist2 < halfTolerance) ? 0. : dist2; 556 } 557 558 ////////////////////////////////////////////// 559 // 560 // Estimate safety distance to the surface fro 561 // 562 G4double G4GenericTrap::DistanceToIn(const G4T 563 { 564 G4double px = p.x(), py = p.y(), pz = p.z(); 565 566 // intersect edges by z = pz plane 567 G4TwoVector pp[4]; 568 G4double z = (pz + fDz); 569 for (auto i = 0; i < 4; ++i) { pp[i] = fVert 570 571 // estimate distance to the solid 572 G4double dist = std::abs(pz) - fDz; 573 for (auto i = 0; i < 4; ++i) 574 { 575 if (fTwist[i] == 0.) 576 { 577 G4double dd = fSurf[i].D*px + fSurf[i].E 578 dist = std::max(dd, dist); 579 } 580 else 581 { 582 // comptute distance to the wedge (two p 583 auto j = (i + 1)%4; 584 G4double cx = pp[j].x() - pp[i].x(); 585 G4double cy = pp[j].y() - pp[i].y(); 586 G4double d = (pp[i].x() - px)*cy + (py - 587 G4ThreeVector na(-cy, cx, fDelta[i].x()* 588 G4ThreeVector nb(-cy, cx, fDelta[j].x()* 589 G4double amag2 = na.mag2(); 590 G4double bmag2 = nb.mag2(); 591 G4double distab = (amag2 > bmag2) ? d/st 592 dist = std::max(distab, dist); 593 } 594 } 595 return (dist > 0.) ? dist : 0.; // return sa 596 } 597 598 ////////////////////////////////////////////// 599 // 600 // Calculate distance to surface from inside 601 // 602 G4double G4GenericTrap::DistanceToOut(const G4 603 const G4 604 const G4 605 G4 606 G4 607 { 608 G4double px = p.x(), py = p.y(), pz = p.z(); 609 G4double vx = v.x(), vy = v.y(), vz = v.z(); 610 611 // Check intersections with plane faces 612 // 613 if ((std::abs(pz) - fDz) >= -halfTolerance & 614 { 615 if (calcNorm) 616 { 617 *validNorm = true; 618 n->set(0., 0., std::copysign(1., pz)); 619 } 620 return 0.; 621 } 622 G4double tout = (vz == 0) ? DBL_MAX : (std:: 623 G4int iface = (vz < 0) ? -4 : -2; // little 624 625 for (auto i = 0; i < 4; ++i) 626 { 627 if (fTwist[i] != 0) continue; // skip twis 628 const G4GenericTrapPlane& surf = fPlane[2* 629 G4double cosa = surf.A*vx + surf.B*vy + su 630 if (cosa > 0) 631 { 632 G4double dist = surf.A*px + surf.B*py + 633 if (dist >= -halfTolerance) 634 { 635 if (calcNorm) 636 { 637 *validNorm = true; 638 n->set(surf.A, surf.B, surf.C); 639 } 640 return 0.; 641 } 642 G4double tmp = -dist/cosa; 643 if (tout > tmp) { tout = tmp; iface = i; 644 } 645 } 646 647 // Check intersections with twisted faces 648 // 649 constexpr G4double EPSILON = 1000.*DBL_EPSIL 650 for (auto i = 0; i < 4; ++i) 651 { 652 if (fTwist[i] == 0) continue; // skip plan 653 654 // twisted face, solve quadratic equation 655 const G4GenericTrapSurface& surf = fSurf[i 656 G4double ABC = surf.A*vx + surf.B*vy + su 657 G4double ABCF = surf.A*px + surf.B*py + su 658 G4double A = ABC*vz; 659 G4double B = 0.5*(surf.D*vx + surf.E*vy + 660 G4double C = surf.D*px + surf.E*py + ABCF* 661 662 if (std::abs(A) <= EPSILON) 663 { 664 // case 1: track is parallel to the surf 665 if (std::abs(B) <= EPSILON) { continue; 666 667 // case 2: single root 668 G4double tmp = -0.5*C/B; 669 // compute normal at intersection point 670 G4double x = px + vx*tmp; 671 G4double y = py + vy*tmp; 672 G4double z = pz + vz*tmp; 673 G4double nx = surf.A*z + surf.D; 674 G4double ny = surf.B*z + surf.E; 675 G4double nz = surf.A*x + surf.B*y + 2.*s 676 677 if (nx*vx + ny*vy + nz*vz > 0.) // point 678 { 679 auto k = (i + 1)%4; 680 G4double tz = (pz + fDz); 681 G4TwoVector a = fVertices[i] + fDelta[ 682 G4TwoVector b = fVertices[k] + fDelta[ 683 G4double dx = b.x() - a.x(); 684 G4double dy = b.y() - a.y(); 685 G4double leng = std::sqrt(dx*dx + dy*d 686 G4double dist = dx*(py - a.y()) - dy*( 687 if (dist >= -halfTolerance*leng) // po 688 { 689 if (calcNorm) 690 { 691 *validNorm = false; 692 G4double normx = surf.A*pz + surf. 693 G4double normy = surf.B*pz + surf. 694 G4double normz = surf.A*px + surf. 695 G4double inv = 1./std::sqrt(normx* 696 n->set(normx*inv, normy*inv, normz 697 } 698 return 0.; 699 } 700 if (tout > tmp) { tout = tmp; iface = 701 } 702 continue; 703 } 704 705 // case 3: scratch or no intersection 706 G4double D = B*B - A*C; 707 if (D < 0.25*fScratch*fScratch*A*A) 708 { 709 // check position of the point 710 auto j = (i + 1)%4; 711 G4double tz = pz + fDz; 712 G4TwoVector a = fVertices[i] + fDelta[i] 713 G4TwoVector b = fVertices[j] + fDelta[j] 714 G4double dx = b.x() - a.x(); 715 G4double dy = b.y() - a.y(); 716 G4double leng = std::sqrt(dx*dx + dy*dy) 717 G4double dist = dx*(py - a.y()) - dy*(px 718 if (dist <= -halfTolerance*leng) { cont 719 if (A > 0 || dist > halfTolerance*leng) 720 { 721 if (calcNorm) 722 { 723 *validNorm = false; 724 G4double nx = surf.A*pz + surf.D; 725 G4double ny = surf.B*pz + surf.E; 726 G4double nz = surf.A*px + surf.B*py 727 G4double inv = 1./std::sqrt(nx*nx + 728 n->set(nx*inv, ny*inv, nz*inv); 729 } 730 return 0.; 731 } 732 continue; 733 } 734 735 // case 4: two intersection points 736 G4double tmp = -B - std::copysign(std::sqr 737 G4double t1 = tmp / A; 738 G4double t2 = C / tmp; 739 G4double tmin = std::min(t1, t2); 740 G4double tmax = std::max(t1, t2); 741 742 if (A < 0) // concave profile: tmin(out) - 743 { 744 if (std::abs(tmax) < std::abs(tmin)) con 745 if (tmin <= halfTolerance) // point is o 746 { 747 G4double t = 0.5*(tmin + tmax); 748 G4double x = px + vx*t; 749 G4double y = py + vy*t; 750 G4double z = pz + vz*t; 751 752 auto j = (i + 1)%4; 753 G4double tz = z + fDz; 754 G4TwoVector a = fVertices[i] + fDelta[ 755 G4TwoVector b = fVertices[j] + fDelta[ 756 G4double dx = b.x() - a.x(); 757 G4double dy = b.y() - a.y(); 758 G4double leng = std::sqrt(dx*dx + dy*d 759 G4double dist = dx*(y - a.y()) - dy*(x 760 if (dist <= -halfTolerance*leng) cont 761 if (calcNorm) 762 { 763 *validNorm = false; 764 G4double nx = surf.A*pz + surf.D; 765 G4double ny = surf.B*pz + surf.E; 766 G4double nz = surf.A*px + surf.B*py 767 G4double inv = 1./std::sqrt(nx*nx + 768 n->set(nx*inv, ny*inv, nz*inv); 769 } 770 return 0.; 771 } 772 if (tout > tmin) { tout = tmin; iface = 773 } 774 else // convex profile: tmin(in) -> tmax(o 775 { 776 if (tmax < halfTolerance) // point is on 777 { 778 if (calcNorm) 779 { 780 *validNorm = false; 781 G4double nx = surf.A*pz + surf.D; 782 G4double ny = surf.B*pz + surf.E; 783 G4double nz = surf.A*px + surf.B*py 784 G4double inv = 1./std::sqrt(nx*nx + 785 n->set(nx*inv, ny*inv, nz*inv); 786 } 787 return 0.; 788 } 789 if (tout > tmax) { tout = tmax; iface = 790 } 791 } 792 793 // Compute normal, if required, and return d 794 // 795 if (tout < halfTolerance) tout = 0.; 796 if (calcNorm) 797 { 798 if (iface < 0) 799 { 800 *validNorm = true; 801 n->set(0, 0, iface + 3); // little trick 802 } 803 else 804 { 805 const G4GenericTrapSurface& surf = fSurf 806 if (fTwist[iface] == 0) 807 { 808 *validNorm = true; 809 n->set(surf.D, surf.E, surf.F); 810 } 811 else 812 { 813 *validNorm = false; 814 G4double x = px + vx*tout; 815 G4double y = py + vy*tout; 816 G4double z = pz + vz*tout; 817 G4double nx = surf.A*z + surf.D; 818 G4double ny = surf.B*z + surf.E; 819 G4double nz = surf.A*x + surf.B*y + 2. 820 G4double inv = 1./std::sqrt(nx*nx + ny 821 n->set(nx*inv, ny*inv, nz*inv); 822 } 823 } 824 } 825 return tout; 826 } 827 828 ////////////////////////////////////////////// 829 // 830 // Estimate safety distance to the surface fro 831 // 832 G4double G4GenericTrap::DistanceToOut(const G4 833 { 834 G4double px = p.x(), py = p.y(), pz = p.z(); 835 836 // intersect edges by z = pz plane 837 G4TwoVector pp[4]; 838 G4double z = (pz + fDz); 839 for (auto i = 0; i < 4; ++i) { pp[i] = fVert 840 841 // estimate distance to the solid 842 G4double dist = std::abs(pz) - fDz; 843 for (auto i = 0; i < 4; ++i) 844 { 845 if (fTwist[i] == 0.) 846 { 847 G4double dd = fSurf[i].D*px + fSurf[i].E 848 dist = std::max(dd, dist); 849 } 850 else 851 { 852 // comptute distance to the wedge (two p 853 auto j = (i + 1)%4; 854 G4double cx = pp[j].x() - pp[i].x(); 855 G4double cy = pp[j].y() - pp[i].y(); 856 G4double d = (pp[i].x() - px)*cy + (py - 857 G4ThreeVector na(-cy, cx, fDelta[i].x()* 858 G4ThreeVector nb(-cy, cx, fDelta[j].x()* 859 G4double amag2 = na.mag2(); 860 G4double bmag2 = nb.mag2(); 861 G4double distab = (amag2 > bmag2) ? d/st 862 dist = std::max(distab, dist); 863 } 864 } 865 return (dist < 0.) ? -dist : 0.; // return s 866 } 867 868 ////////////////////////////////////////////// 869 // 870 // Return bounding box 871 // 872 void G4GenericTrap::BoundingLimits(G4ThreeVect 873 G4ThreeVect 874 { 875 pMin = fMinBBox; 876 pMax = fMaxBBox; 877 } 878 879 ////////////////////////////////////////////// 880 // 881 // Calculate extent under transform and specif 882 // 883 G4bool 884 G4GenericTrap::CalculateExtent(const EAxis pAx 885 const G4VoxelLi 886 const G4AffineT 887 G4double& 888 { 889 G4ThreeVector bmin, bmax; 890 G4bool exist; 891 892 // Check bounding box (bbox) 893 // 894 BoundingLimits(bmin,bmax); 895 G4BoundingEnvelope bbox(bmin,bmax); 896 #ifdef G4BBOX_EXTENT 897 return bbox.CalculateExtent(pAxis,pVoxelLimi 898 #endif 899 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox 900 { 901 return exist = pMin < pMax; 902 } 903 904 // Set bounding envelope (benv) and calculat 905 // 906 // To build the bounding envelope with plane 907 // the trapezoid is subdivided in two triang 908 // duplication of vertices in the bases in a 909 // a convex polyhedron (some faces of the en 910 // 911 G4double dz = GetZHalfLength(); 912 G4ThreeVectorList baseA(8), baseB(8); 913 for (G4int i = 0; i < 4; ++i) 914 { 915 G4TwoVector va = GetVertex(i); 916 G4TwoVector vb = GetVertex(i+4); 917 baseA[2*i].set(va.x(), va.y(),-dz); 918 baseB[2*i].set(vb.x(), vb.y(), dz); 919 } 920 for (G4int i = 0; i < 4; ++i) 921 { 922 G4int k1 = 2*i, k2 = (2*i + 2)%8; 923 G4double ax = (baseA[k2].x() - baseA[k1].x 924 G4double ay = (baseA[k2].y() - baseA[k1].y 925 G4double bx = (baseB[k2].x() - baseB[k1].x 926 G4double by = (baseB[k2].y() - baseB[k1].y 927 G4double znorm = ax*by - ay*bx; 928 baseA[k1+1] = (znorm < 0.0) ? baseA[k2] : 929 baseB[k1+1] = (znorm < 0.0) ? baseB[k1] : 930 } 931 932 std::vector<const G4ThreeVectorList *> polyg 933 polygons[0] = &baseA; 934 polygons[1] = &baseB; 935 936 G4BoundingEnvelope benv(bmin, bmax, polygons 937 exist = benv.CalculateExtent(pAxis,pVoxelLim 938 return exist; 939 } 940 941 ////////////////////////////////////////////// 942 // 943 // Return type of the solid 944 // 945 G4GeometryType G4GenericTrap::GetEntityType() 946 { 947 return { "G4GenericTrap" }; 948 } 949 950 ////////////////////////////////////////////// 951 // 952 // Return true if not twisted, false otherwise 953 // 954 G4bool G4GenericTrap::IsFaceted() const 955 { 956 return (!fIsTwisted); 957 } 958 959 ////////////////////////////////////////////// 960 // 961 // Make a clone of the solid 962 // 963 G4VSolid* G4GenericTrap::Clone() const 964 { 965 return new G4GenericTrap(*this); 966 } 967 968 ////////////////////////////////////////////// 969 // 970 // Write parameters of the solid to the specif 971 // 972 std::ostream& G4GenericTrap::StreamInfo(std::o 973 { 974 G4long oldprc = os.precision(16); 975 os << "------------------------------------- 976 << " *** Dump for solid - " << GetName 977 << " ================================= 978 << "Solid geometry type: " << GetEntityTy 979 << " half length Z: " << fDz/mm << "\n" 980 << " list of vertices:\n"; 981 for (G4int i = 0; i < 8; ++i) 982 { 983 os << " #" << i << " " << fVertices[i] 984 } 985 os << "------------------------------------- 986 os.precision(oldprc); 987 return os; 988 } 989 990 ////////////////////////////////////////////// 991 // 992 // Pick up a random point on the surface 993 // 994 G4ThreeVector G4GenericTrap::GetPointOnSurface 995 { 996 if (fArea[0] + fArea[1] + fArea[2] + fArea[3 997 { 998 G4AutoLock l(&lateralareaMutex); 999 fArea[0] = GetLateralFaceArea(0); 1000 fArea[1] = GetLateralFaceArea(1); 1001 fArea[2] = GetLateralFaceArea(2); 1002 fArea[3] = GetLateralFaceArea(3); 1003 l.unlock(); 1004 } 1005 1006 constexpr G4int iface[6][4] = 1007 { {0,1,2,3}, {0,4,5,1}, {1,5,6,2}, {2,6,7 1008 1009 G4bool isTwisted[6] = {false}; 1010 for (auto i = 0; i < 4; ++i) { isTwisted[i 1011 1012 G4double ssurf[6]; 1013 G4TwoVector A = fVertices[3] - fVertices[1] 1014 G4TwoVector B = fVertices[2] - fVertices[0] 1015 G4TwoVector C = fVertices[7] - fVertices[5] 1016 G4TwoVector D = fVertices[6] - fVertices[4] 1017 ssurf[0] = (A.x()*B.y() - A.y()*B.x())*0.5; 1018 ssurf[1] = ssurf[0] + fArea[0]; 1019 ssurf[2] = ssurf[1] + fArea[1]; 1020 ssurf[3] = ssurf[2] + fArea[2]; 1021 ssurf[4] = ssurf[3] + fArea[3]; 1022 ssurf[5] = ssurf[4] + (C.x()*D.y() - C.y()* 1023 1024 G4double select = ssurf[5]*G4QuickRand(); 1025 G4int k; 1026 for (k = 0; k < 5; ++k) { if (select <= ssu 1027 1028 G4int i0 = iface[k][0]; 1029 G4int i1 = iface[k][1]; 1030 G4int i2 = iface[k][2]; 1031 G4int i3 = iface[k][3]; 1032 G4ThreeVector pp[4]; 1033 pp[0].set(fVertices[i0].x(), fVertices[i0]. 1034 pp[1].set(fVertices[i1].x(), fVertices[i1]. 1035 pp[2].set(fVertices[i2].x(), fVertices[i2]. 1036 pp[3].set(fVertices[i3].x(), fVertices[i3]. 1037 1038 G4ThreeVector point; 1039 if (isTwisted[k]) 1040 { // twisted face, rejection sampling 1041 G4double maxlength = std::max((pp[2] - pp 1042 point = (pp[0] + pp[1] + pp[2] + pp[3])*0 1043 for (auto i = 0; i < 10000; ++i) 1044 { 1045 G4double u = G4QuickRand(); 1046 G4ThreeVector p0 = pp[0] + (pp[1] - pp[ 1047 G4ThreeVector p1 = pp[3] + (pp[2] - pp[ 1048 G4double v = G4QuickRand()*(maxlength/( 1049 if (v >= 1.) continue; 1050 point = p0 + (p1 - p0)*v; 1051 break; 1052 } 1053 } 1054 else 1055 { // plane face 1056 G4double u = G4QuickRand(); 1057 G4double v = G4QuickRand(); 1058 if (u + v > 1.) { u = 1. - u; v = 1. - v; 1059 G4double ss = (((pp[2] - pp[0]).cross(pp[ 1060 G4int j = (select > ssurf[k] - ss) ? 3 : 1061 point = pp[0] + (pp[2] - pp[0])*u + (pp[j 1062 } 1063 return point; 1064 } 1065 1066 ///////////////////////////////////////////// 1067 // 1068 // Return volume of the solid 1069 // 1070 G4double G4GenericTrap::GetCubicVolume() 1071 { 1072 if (fCubicVolume == 0.0) 1073 { 1074 // diagonals 1075 G4TwoVector A = fVertices[3] - fVertices[ 1076 G4TwoVector B = fVertices[2] - fVertices[ 1077 G4TwoVector C = fVertices[7] - fVertices[ 1078 G4TwoVector D = fVertices[6] - fVertices[ 1079 1080 // kross products 1081 G4double AB = A.x()*B.y() - A.y()*B.x(); 1082 G4double CD = C.x()*D.y() - C.y()*D.x(); 1083 G4double AD = A.x()*D.y() - A.y()*D.x(); 1084 G4double CB = C.x()*B.y() - C.y()*B.x(); 1085 1086 fCubicVolume = ((AB + CD)/3. + (AD + CB)/ 1087 } 1088 return fCubicVolume; 1089 } 1090 1091 ///////////////////////////////////////////// 1092 // 1093 // Compute lateral face area 1094 // 1095 G4double G4GenericTrap::GetLateralFaceArea(G4 1096 { 1097 G4int i1 = iface, i2 = (i1 + 1)%4, i3 = i1 1098 1099 // plane face 1100 if (fTwist[iface] == 0.0) 1101 { 1102 G4ThreeVector p1(fVertices[i1].x(), fVert 1103 G4ThreeVector p2(fVertices[i2].x(), fVert 1104 G4ThreeVector p3(fVertices[i3].x(), fVert 1105 G4ThreeVector p4(fVertices[i4].x(), fVert 1106 return ((p4 - p1).cross(p3 - p2)).mag()*0 1107 } 1108 1109 // twisted face 1110 constexpr G4int NSTEP = 250; 1111 constexpr G4double dt = 1./NSTEP; 1112 1113 G4double x21 = fVertices[i2].x() - fVertice 1114 G4double y21 = fVertices[i2].y() - fVertice 1115 G4double x31 = fVertices[i3].x() - fVertice 1116 G4double y31 = fVertices[i3].y() - fVertice 1117 G4double x42 = fVertices[i4].x() - fVertice 1118 G4double y42 = fVertices[i4].y() - fVertice 1119 G4double x43 = fVertices[i4].x() - fVertice 1120 G4double y43 = fVertices[i4].y() - fVertice 1121 1122 G4double A = x21*y43 - y21*x43; 1123 G4double B0 = x21*y31 - y21*x31; 1124 G4double B1 = x42*y31 - y42*x31; 1125 G4double HH = 4*fDz*fDz; 1126 G4double invAA = 1./(A*A); 1127 G4double sqrtAA = 2.*std::abs(A); 1128 G4double invSqrtAA = 1./sqrtAA; 1129 1130 G4double area = 0.; 1131 for (G4int i = 0; i < NSTEP; ++i) 1132 { 1133 G4double t = (i + 0.5)*dt; 1134 G4double I = y21 + (y43 - y21)*t; 1135 G4double J = x21 + (x43 - x21)*t; 1136 G4double IIJJ = HH*(I*I + J*J); 1137 G4double B = B1*t + B0; 1138 1139 G4double aa = A*A; 1140 G4double bb = 2.*A*B; 1141 G4double cc = IIJJ + B*B; 1142 1143 G4double R1 = std::sqrt(aa + bb + cc); 1144 G4double R0 = std::sqrt(cc); 1145 G4double log1 = std::log(std::abs(sqrtAA* 1146 G4double log0 = std::log(std::abs(sqrtAA* 1147 1148 area += 0.5*R1 + 0.25*bb*invAA*(R1 - R0) 1149 } 1150 return area*dt; 1151 } 1152 1153 ///////////////////////////////////////////// 1154 // 1155 // Return surface area of the solid 1156 // 1157 G4double G4GenericTrap::GetSurfaceArea() 1158 { 1159 if (fSurfaceArea == 0.0) 1160 { 1161 G4TwoVector A = fVertices[3] - fVertices[ 1162 G4TwoVector B = fVertices[2] - fVertices[ 1163 G4TwoVector C = fVertices[7] - fVertices[ 1164 G4TwoVector D = fVertices[6] - fVertices[ 1165 G4double S_bot = (A.x()*B.y() - A.y()*B.x 1166 G4double S_top = (C.x()*D.y() - C.y()*D.x 1167 fArea[0] = GetLateralFaceArea(0); 1168 fArea[1] = GetLateralFaceArea(1); 1169 fArea[2] = GetLateralFaceArea(2); 1170 fArea[3] = GetLateralFaceArea(3); 1171 fSurfaceArea = S_bot + S_top + fArea[0] + 1172 } 1173 return fSurfaceArea; 1174 } 1175 1176 ///////////////////////////////////////////// 1177 // 1178 // Check parameters of the solid 1179 // 1180 void 1181 G4GenericTrap::CheckParameters(G4double halfZ 1182 const std::vec 1183 { 1184 // Check number of vertices 1185 // 1186 if (vertices.size() != 8) 1187 { 1188 std::ostringstream message; 1189 message << "Number of vertices is " << ve 1190 << " instead of 8 for Solid: " << 1191 G4Exception("G4GenericTrap::CheckParamete 1192 FatalException, message); 1193 } 1194 1195 // Check dZ 1196 // 1197 if ((fDz = halfZ) < 2.*kCarTolerance) 1198 { 1199 std::ostringstream message; 1200 message << "Z-dimension is too small or n 1201 << " for Solid: " << GetName() << 1202 G4Exception("G4GenericTrap::CheckParamete 1203 FatalException, message); 1204 } 1205 1206 // Check order of vertices 1207 // 1208 for (auto i = 0; i < 8; ++i) { fVertices.at 1209 1210 // Bottom polygon area 1211 G4TwoVector diag1 = fVertices[2] - fVertice 1212 G4TwoVector diag2 = fVertices[3] - fVertice 1213 G4double ldiagbot = std::max(diag1.mag(), d 1214 G4double zbot = diag1.x()*diag2.y() - diag1 1215 if (std::abs(zbot) < ldiagbot*kCarTolerance 1216 1217 // Top polygon area 1218 G4TwoVector diag3 = fVertices[6] - fVertice 1219 G4TwoVector diag4 = fVertices[7] - fVertice 1220 G4double ldiagtop = std::max(diag3.mag(), d 1221 G4double ztop = diag3.x()*diag4.y() - diag3 1222 if (std::abs(ztop) < ldiagtop*kCarTolerance 1223 1224 if (zbot*ztop < 0.) 1225 { 1226 std::ostringstream message; 1227 message << "Vertices of the bottom and to 1228 StreamInfo(message); 1229 G4Exception("G4GenericTrap::CheckParamete 1230 FatalException, message); 1231 } 1232 if ((zbot > 0.) || (ztop > 0.)) 1233 { 1234 std::swap(fVertices[1], fVertices[3]); 1235 std::swap(fVertices[5], fVertices[7]); 1236 std::ostringstream message; 1237 message << "Vertices re-ordered in Solid: 1238 << "Vertices of the bottom and to 1239 G4Exception("G4GenericTrap::CheckParamete 1240 JustWarning, message); 1241 } 1242 1243 // Check for degeneracy 1244 // 1245 G4int ndegenerated = 0; 1246 for (auto i = 0; i < 4; ++i) 1247 { 1248 auto j = (i + 1)%4; 1249 G4double lbot = (fVertices[j] - fVertices 1250 G4double ltop = (fVertices[j + 4] - fVert 1251 ndegenerated += (std::max(lbot, ltop) < k 1252 } 1253 if (ndegenerated > 1 || 1254 GetCubicVolume() < fDz*std::max(ldiagbo 1255 { 1256 std::ostringstream message; 1257 message << "Degenerated solid !\n"; 1258 StreamInfo(message); 1259 G4Exception("G4GenericTrap::CheckParamete 1260 FatalException, message); 1261 } 1262 1263 // Check that the polygons are convex 1264 // 1265 G4bool isConvex = true; 1266 for (auto i = 0; i < 4; ++i) 1267 { 1268 auto j = (i + 1)%4; 1269 auto k = (j + 1)%4; 1270 G4TwoVector edge1 = fVertices[j] - fVerti 1271 G4TwoVector edge2 = fVertices[k] - fVerti 1272 isConvex = ((edge1.x()*edge2.y() - edge1. 1273 if (!isConvex) break; 1274 G4TwoVector edge3 = fVertices[j + 4] - fV 1275 G4TwoVector edge4 = fVertices[k + 4] - fV 1276 isConvex = ((edge3.x()*edge4.y() - edge3. 1277 if (!isConvex) break; 1278 } 1279 if (!isConvex) 1280 { 1281 std::ostringstream message; 1282 message << "The bottom and top faces must 1283 StreamInfo(message); 1284 G4Exception("G4GenericTrap::CheckParamete 1285 FatalException, message); 1286 } 1287 } 1288 1289 ///////////////////////////////////////////// 1290 // 1291 // Compute surface equations and twist angles 1292 // 1293 void G4GenericTrap::ComputeLateralSurfaces() 1294 { 1295 for (auto i = 0; i < 4; ++i) 1296 { 1297 auto j = (i + 1)%4; 1298 G4ThreeVector p1(fVertices[j].x(), fVerti 1299 G4ThreeVector p2(fVertices[i].x(), fVerti 1300 G4ThreeVector p3(fVertices[j + 4].x(), fV 1301 G4ThreeVector p4(fVertices[i + 4].x(), fV 1302 G4ThreeVector ebot = p2 - p1; 1303 G4ThreeVector etop = p4 - p3; 1304 G4double lbot = ebot.mag(); 1305 G4double ltop = etop.mag(); 1306 G4double zcross = ebot.x()*etop.y() - ebo 1307 G4double eps = kCarTolerance*std::max(lbo 1308 if (std::min(lbot, ltop) < kCarTolerance 1309 { // plane surface: Dx + Ey + Fz + G = 0 1310 G4ThreeVector normal; 1311 if (std::max(lbot, ltop) < kCarToleranc 1312 { 1313 auto k = (j + 1)%4; 1314 auto l = (k + 1)%4; 1315 G4TwoVector vl = fVertices[l] + fVert 1316 G4TwoVector vi = fVertices[i] + fVert 1317 G4TwoVector vj = fVertices[j] + fVert 1318 G4TwoVector vk = fVertices[k] + fVert 1319 G4TwoVector vij = (vi - vl).unit() + 1320 G4ThreeVector epar = (p4 + p3 - p2 - 1321 G4ThreeVector eort = epar.cross(G4Thr 1322 normal = (eort.cross(epar)).unit(); 1323 } 1324 else 1325 { 1326 normal = ((p4 - p1).cross(p3 - p2)).u 1327 } 1328 fSurf[i].D = fPlane[2*i].A = fPlane[2*i 1329 fSurf[i].E = fPlane[2*i].B = fPlane[2*i 1330 fSurf[i].F = fPlane[2*i].C = fPlane[2*i 1331 fSurf[i].G = fPlane[2*i].D = fPlane[2*i 1332 } 1333 else 1334 { // hyperbolic paraboloid: Axz + Byz + C 1335 fIsTwisted = true; 1336 G4double angle = std::acos(ebot.dot(eto 1337 if (angle > CLHEP::halfpi) 1338 { 1339 std::ostringstream message; 1340 message << "Twist on " << angle/CLHEP 1341 << " degrees, should not be m 1342 StreamInfo(message); 1343 G4Exception("G4GenericTrap::ComputeLa 1344 FatalException, message); 1345 } 1346 fTwist[i] = std::copysign(angle, zcross 1347 // set equation of twisted surface (hyp 1348 fSurf[i].A = 2.*fDz*(p4.y() - p3.y() - 1349 fSurf[i].B =-2.*fDz*(p4.x() - p3.x() - 1350 fSurf[i].C = ((p4.x() - p2.x())*(p3.y() 1351 fSurf[i].D = 2.*fDz*fDz*(p4.y() - p3.y( 1352 fSurf[i].E =-2.*fDz*fDz*(p4.x() - p3.x( 1353 fSurf[i].F = 2.*fDz*(p4.x()*p3.y() - p3 1354 fSurf[i].G = fDz*fDz*((p4.x() + p2.x()) 1355 G4double magnitude = G4ThreeVector(fSu 1356 if (magnitude < kCarTolerance) continue 1357 fSurf[i].A /= magnitude; 1358 fSurf[i].B /= magnitude; 1359 fSurf[i].C /= magnitude; 1360 fSurf[i].D /= magnitude; 1361 fSurf[i].E /= magnitude; 1362 fSurf[i].F /= magnitude; 1363 fSurf[i].G /= magnitude; 1364 // set planes of bounding polyhedron 1365 G4ThreeVector normal1, normal2; 1366 G4ThreeVector c1, c2; 1367 if (fTwist[i] < 0.) 1368 { 1369 normal1 = ((p2 - p1).cross(p4 - p1)). 1370 normal2 = ((p3 - p4).cross(p1 - p4)). 1371 c1 = p1; 1372 c2 = p4; 1373 } 1374 else 1375 { 1376 normal1 = ((p3 - p2).cross(p1 - p2)). 1377 normal2 = ((p2 - p3).cross(p4 - p3)). 1378 c1 = p2; 1379 c2 = p3; 1380 } 1381 fPlane[2*i].A = normal1.x(); 1382 fPlane[2*i].B = normal1.y(); 1383 fPlane[2*i].C = normal1.z(); 1384 fPlane[2*i].D = -normal1.dot(c1); 1385 fPlane[2*i + 1].A = normal2.x(); 1386 fPlane[2*i + 1].B = normal2.y(); 1387 fPlane[2*i + 1].C = normal2.z(); 1388 fPlane[2*i + 1].D = -normal2.dot(c2); 1389 } 1390 fDelta[i] = (fVertices[i + 4] - fVertices 1391 } 1392 } 1393 1394 ///////////////////////////////////////////// 1395 // 1396 // Set bounding box 1397 // 1398 void G4GenericTrap::ComputeBoundingBox() 1399 { 1400 G4double minX, maxX, minY, maxY; 1401 minX = maxX = fVertices[0].x(); 1402 minY = maxY = fVertices[0].y(); 1403 for (auto i = 1; i < 8; ++i) 1404 { 1405 minX = std::min(minX, fVertices[i].x()); 1406 maxX = std::max(maxX, fVertices[i].x()); 1407 minY = std::min(minY, fVertices[i].y()); 1408 maxY = std::max(maxY, fVertices[i].y()); 1409 } 1410 fMinBBox = G4ThreeVector(minX, minY,-fDz); 1411 fMaxBBox = G4ThreeVector(maxX, maxY, fDz); 1412 1413 // Check correctness of the bounding box 1414 // 1415 if (minX >= maxX || minY >= maxY || -fDz >= 1416 { 1417 std::ostringstream message; 1418 message << "Bad bounding box (min >= max) 1419 << GetName() << " !" 1420 << "\npMin = " << fMinBBox 1421 << "\npMax = " << fMaxBBox; 1422 G4Exception("G4GenericTrap::ComputeBoundi 1423 JustWarning, message); 1424 DumpInfo(); 1425 } 1426 } 1427 1428 ///////////////////////////////////////////// 1429 // 1430 // Set max length of a scratch 1431 // 1432 void G4GenericTrap::ComputeScratchLength() 1433 { 1434 G4double scratch = kInfinity; 1435 for (auto i = 0; i < 4; ++i) 1436 { 1437 if (fTwist[i] == 0.) continue; // skip pl 1438 auto k = (i + 1)%4; 1439 G4ThreeVector p1(fVertices[i].x(), fVerti 1440 G4ThreeVector p2(fVertices[k].x(), fVerti 1441 G4ThreeVector p3(fVertices[i + 4].x(), fV 1442 G4ThreeVector p4(fVertices[k + 4].x(), fV 1443 G4ThreeVector p0 = (p1 + p2 + p3 + p4)*0. 1444 G4ThreeVector norm = SurfaceNormal(p0); 1445 G4ThreeVector pp[2]; // points inside and 1446 pp[0] = p0 - norm * halfTolerance; 1447 pp[1] = p0 + norm * halfTolerance; 1448 G4ThreeVector vv[2]; // unit vectors alon 1449 vv[0] = (p4 - p1).unit(); 1450 vv[1] = (p3 - p2).unit(); 1451 // find intersection points and compute t 1452 for (auto ip = 0; ip < 2; ++ip) 1453 { 1454 G4double px = pp[ip].x(); 1455 G4double py = pp[ip].y(); 1456 G4double pz = pp[ip].z(); 1457 for (auto iv = 0; iv < 2; ++iv) 1458 { 1459 G4double vx = vv[iv].x(); 1460 G4double vy = vv[iv].y(); 1461 G4double vz = vv[iv].z(); 1462 // solve quadratic equation 1463 G4double ABC = fSurf[i].A*vx + fSurf 1464 G4double ABCF = fSurf[i].A*px + fSurf 1465 G4double A = ABC*vz; 1466 G4double B = 0.5*(fSurf[i].D*vx + fSu 1467 G4double C = fSurf[i].D*px + fSurf[i] 1468 G4double D = B*B - A*C; 1469 if (D < 0) continue; 1470 G4double leng = 2.*std::sqrt(D)/std:: 1471 scratch = std::min(leng, scratch); 1472 } 1473 } 1474 } 1475 fScratch = std::max(kCarTolerance, scratch) 1476 } 1477 1478 ///////////////////////////////////////////// 1479 // 1480 // GetPolyhedron 1481 // 1482 G4Polyhedron* G4GenericTrap::GetPolyhedron () 1483 { 1484 if ( (fpPolyhedron == nullptr) 1485 || fRebuildPolyhedron 1486 || (fpPolyhedron->GetNumberOfRotationStep 1487 fpPolyhedron->GetNumberOfRotationStep 1488 { 1489 G4AutoLock l(&polyhedronMutex); 1490 delete fpPolyhedron; 1491 fpPolyhedron = CreatePolyhedron(); 1492 fRebuildPolyhedron = false; 1493 l.unlock(); 1494 } 1495 return fpPolyhedron; 1496 } 1497 1498 ///////////////////////////////////////////// 1499 // 1500 // Method for visualisation 1501 // 1502 void G4GenericTrap::DescribeYourselfTo(G4VGra 1503 { 1504 scene.AddSolid(*this); 1505 } 1506 1507 ///////////////////////////////////////////// 1508 // 1509 // Return VisExtent 1510 // 1511 G4VisExtent G4GenericTrap::GetExtent() const 1512 { 1513 return { fMinBBox.x(), fMaxBBox.x(), 1514 fMinBBox.y(), fMaxBBox.y(), 1515 fMinBBox.z(), fMaxBBox.z() }; 1516 } 1517 1518 // ------------------------------------------ 1519 1520 G4Polyhedron* G4GenericTrap::CreatePolyhedron 1521 { 1522 // Approximation of Twisted Side 1523 // Construct extra Points, if Twisted Side 1524 // 1525 G4Polyhedron* polyhedron; 1526 G4int nVertices, nFacets; 1527 1528 G4int subdivisions = 0; 1529 if (fIsTwisted) 1530 { 1531 if (GetVisSubdivisions() != 0) 1532 { 1533 subdivisions = GetVisSubdivisions(); 1534 } 1535 else 1536 { 1537 // Estimation of Number of Subdivisions 1538 // 1539 G4double maxTwist = 0.; 1540 for(G4int i = 0; i < 4; ++i) 1541 { 1542 if (GetTwistAngle(i) > maxTwist) { ma 1543 } 1544 1545 // Computes bounding vectors for the sh 1546 // 1547 G4double Dx, Dy; 1548 Dx = 0.5*(fMaxBBox.x() - fMinBBox.y()); 1549 Dy = 0.5*(fMaxBBox.y() - fMinBBox.y()); 1550 if (Dy > Dx) { Dx = Dy; } 1551 1552 subdivisions = 8*G4int(maxTwist/(Dx*Dx* 1553 if (subdivisions < 4) { subdivisions = 1554 if (subdivisions > 30) { subdivisions = 1555 } 1556 } 1557 G4int sub4 = 4*subdivisions; 1558 nVertices = 8 + subdivisions*4; 1559 nFacets = 6 + subdivisions*4; 1560 G4double cf = 1./(subdivisions + 1); 1561 polyhedron = new G4Polyhedron(nVertices, nF 1562 1563 // Set vertices 1564 // 1565 G4int icur = 0; 1566 for (G4int i = 0; i < 4; ++i) 1567 { 1568 G4ThreeVector v(fVertices[i].x(),fVertice 1569 polyhedron->SetVertex(++icur, v); 1570 } 1571 for (G4int i = 0; i < subdivisions; ++i) 1572 { 1573 for (G4int j = 0; j < 4; ++j) 1574 { 1575 G4TwoVector u = fVertices[j]+cf*(i+1)*( 1576 G4ThreeVector v(u.x(),u.y(),-fDz+cf*2*f 1577 polyhedron->SetVertex(++icur, v); 1578 } 1579 } 1580 for (G4int i = 4; i < 8; ++i) 1581 { 1582 G4ThreeVector v(fVertices[i].x(),fVertice 1583 polyhedron->SetVertex(++icur, v); 1584 } 1585 1586 // Set facets 1587 // 1588 icur = 0; 1589 polyhedron->SetFacet(++icur, 1, 4, 3, 2); / 1590 for (G4int i = 0; i < subdivisions + 1; ++i 1591 { 1592 G4int is = i*4; 1593 polyhedron->SetFacet(++icur, 5+is, 8+is, 1594 polyhedron->SetFacet(++icur, 8+is, 7+is, 1595 polyhedron->SetFacet(++icur, 7+is, 6+is, 1596 polyhedron->SetFacet(++icur, 6+is, 5+is, 1597 } 1598 polyhedron->SetFacet(++icur, 5+sub4, 6+sub4 1599 1600 polyhedron->SetReferences(); 1601 polyhedron->InvertFacets(); 1602 1603 return polyhedron; 1604 } 1605 1606 ///////////////////////////////////////////// 1607 // 1608 // Print out a warning if A has an unexpected 1609 // 1610 void G4GenericTrap::WarningSignA(const G4Stri 1611 const G4Thre 1612 { 1613 std::ostringstream message; 1614 message.precision(16); 1615 message << icase << " in " << GetName() << 1616 << " p" << p << "\n" 1617 << " v" << v << "\n" 1618 << " A = " << A << " (is " 1619 << ((A < 0) ? "negative, instead of 1620 StreamInfo(message); 1621 const G4String function = "G4GenericTrap::D 1622 G4Exception(function, "GeomSolids1002", Jus 1623 } 1624 1625 ///////////////////////////////////////////// 1626 // 1627 // Print out a warning if B has an unexpected 1628 // 1629 void G4GenericTrap::WarningSignB(const G4Stri 1630 G4double f, 1631 const G4Thre 1632 { 1633 std::ostringstream message; 1634 message.precision(16); 1635 message << icase << " in " << GetName() << 1636 << " p" << p << "\n" 1637 << " v" << v << "\n" 1638 << " f = " << f << " B = " << B < 1639 << ((B < 0) ? "negative, instead of 1640 StreamInfo(message); 1641 const G4String function = "G4GenericTrap::D 1642 G4Exception(function, "GeomSolids1002", Jus 1643 } 1644 1645 ///////////////////////////////////////////// 1646 // 1647 // Print out a warning in DistanceToIn(p,v) 1648 // 1649 void G4GenericTrap::WarningDistanceToIn(G4int 1650 G4dou 1651 const 1652 { 1653 G4String check = ""; 1654 if (ttin[1] != DBL_MAX) 1655 { 1656 G4double tcheck = 0.5*(ttout[0] + ttin[1] 1657 if (Inside(p + v*tcheck) != kOutside) che 1658 } 1659 1660 auto position = Inside(p); 1661 std::ostringstream message; 1662 message.precision(16); 1663 message << k << "_Unexpected sequence of in 1664 << " position = " << ((position = 1665 << " p" << p << "\n" 1666 << " v" << v << "\n" 1667 << " range : [" << tmin << ", 1668 << " ttin[2] : " 1669 << ((ttin[0] == DBL_MAX) ? kInfinit 1670 << ((ttin[1] == DBL_MAX) ? kInfinit 1671 << " ttout[2] : " 1672 << ((ttout[0] == DBL_MAX) ? kInfini 1673 << ((ttout[1] == DBL_MAX) ? kInfini 1674 StreamInfo(message); 1675 G4Exception("G4GenericTrap::DistanceToIn(p, 1676 JustWarning, message ); 1677 } 1678 1679 ///////////////////////////////////////////// 1680 // 1681 // Print out a warning in DistanceToOut(p,v) 1682 // 1683 void G4GenericTrap::WarningDistanceToOut(cons 1684 cons 1685 G4do 1686 { 1687 auto position = Inside(p); 1688 std::ostringstream message; 1689 message.precision(16); 1690 message << "Unexpected final tout = " << to 1691 << " position = " << ((position = 1692 << " p" << p << "\n" 1693 << " v" << v << "\n"; 1694 StreamInfo(message); 1695 G4Exception("G4GenericTrap::DistanceToOut(p 1696 JustWarning, message ); 1697 } 1698 1699 #endif 1700