Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer << 3 // * DISCLAIMER * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers * 6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which * 7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: * 8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license * 9 // * include a list of copyright holders. << 10 // * 9 // * * 11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file << 14 // * use. * 16 // * for the full disclaimer and the limitatio << 17 // * 15 // * * 18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the * 19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. * 20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work * 21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this * 22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. * 23 // * acceptance of all terms of the Geant4 Sof << 24 // ******************************************* 21 // ******************************************************************** 25 // 22 // 26 // G4TwistedTrap implementation << 23 // $Id: G4TwistedTrap.cc,v 1.5 2004/12/10 16:22:39 gcosmo Exp $ >> 24 // GEANT4 tag $Name: geant4-07-00-patch-01 $ 27 // 25 // 28 // Author: 10/11/2004 - O.Link (Oliver.Link@ce << 26 // 29 // ------------------------------------------- 27 // -------------------------------------------------------------------- >> 28 // GEANT 4 class source file >> 29 // >> 30 // >> 31 // G4TwistedTrap.cc >> 32 // >> 33 // Author: >> 34 // >> 35 // 04-Nov-2004 - O.Link (Oliver.Link@cern.ch) >> 36 // >> 37 // -------------------------------------------------------------------- >> 38 >> 39 // #define DISTANCETOIN 30 40 31 #include "G4TwistedTrap.hh" 41 #include "G4TwistedTrap.hh" 32 #include "G4SystemOfUnits.hh" << 42 >> 43 #include "G4VoxelLimits.hh" >> 44 #include "G4AffineTransform.hh" >> 45 #include "G4SolidExtentList.hh" >> 46 #include "G4ClippablePolygon.hh" >> 47 #include "G4VPVParameterisation.hh" >> 48 #include "meshdefs.hh" >> 49 >> 50 #include "G4VGraphicsScene.hh" 33 #include "G4Polyhedron.hh" 51 #include "G4Polyhedron.hh" >> 52 #include "G4VisExtent.hh" >> 53 #include "G4NURBS.hh" >> 54 #include "G4NURBStube.hh" >> 55 #include "G4NURBScylinder.hh" >> 56 #include "G4NURBStubesector.hh" >> 57 >> 58 //===================================================================== >> 59 //* constructors ------------------------------------------------------ >> 60 >> 61 G4TwistedTrap::G4TwistedTrap(const G4String &pname, >> 62 G4double twistedangle, >> 63 G4double pDx1, >> 64 G4double pDx2, >> 65 G4double pDy, >> 66 G4double pDz) >> 67 >> 68 : G4VSolid(pname), >> 69 fLowerEndcap(0), fUpperEndcap(0), fSide0(0), >> 70 fSide90(0), fSide180(0), fSide270(0) >> 71 { >> 72 >> 73 if ( ( pDx1 > 2*kCarTolerance) >> 74 && ( pDx2 > 2*kCarTolerance) >> 75 && ( pDy > 2*kCarTolerance) >> 76 && ( pDz > 2*kCarTolerance) >> 77 && ( std::fabs(twistedangle) > 2*kAngTolerance ) >> 78 && ( std::fabs(twistedangle) < halfpi ) ) >> 79 { >> 80 >> 81 SetFields(twistedangle, pDx1, pDx2, pDy, pDz); >> 82 CreateSurfaces(); >> 83 fCubicVolume = 8 * ( pDx1 + pDx2 ) / 2 * pDy * pDz ; >> 84 >> 85 } >> 86 else >> 87 { >> 88 G4cerr << "ERROR - G4TwistedTrap()::G4TwistedTrap(): " << GetName() << G4endl >> 89 << " Dimensions too small ! - " >> 90 << pDx1 << ", " << pDx2 << ", " << pDy << ", " << pDz << G4endl >> 91 << " twistangle " << twistedangle/deg << " deg" << G4endl ; >> 92 >> 93 G4Exception("G4TwistedTrap::G4TwistedTrap()", "InvalidSetup", >> 94 FatalException, "Invalid dimensions. Too small, or twist angle too big."); >> 95 } >> 96 >> 97 } >> 98 >> 99 //===================================================================== >> 100 //* destructor -------------------------------------------------------- >> 101 >> 102 G4TwistedTrap::~G4TwistedTrap() >> 103 { >> 104 >> 105 if (fLowerEndcap) delete fLowerEndcap ; >> 106 if (fUpperEndcap) delete fUpperEndcap ; >> 107 >> 108 if (fSide0) delete fSide0 ; >> 109 if (fSide90) delete fSide90 ; >> 110 if (fSide180) delete fSide180 ; >> 111 if (fSide270) delete fSide270 ; >> 112 >> 113 >> 114 } >> 115 >> 116 //===================================================================== >> 117 //* ComputeDimensions ------------------------------------------------- >> 118 >> 119 void G4TwistedTrap::ComputeDimensions(G4VPVParameterisation* , >> 120 const G4int , >> 121 const G4VPhysicalVolume* ) >> 122 { >> 123 G4Exception("G4TwistedTubs::ComputeDimensions()", >> 124 "NotSupported", FatalException, >> 125 "G4TwistedTubs does not support Parameterisation."); >> 126 >> 127 } >> 128 >> 129 //===================================================================== >> 130 //* CalculateExtent --------------------------------------------------- >> 131 >> 132 G4bool G4TwistedTrap::CalculateExtent( const EAxis pAxis, >> 133 const G4VoxelLimits &pVoxelLimit, >> 134 const G4AffineTransform &pTransform, >> 135 G4double &pMin, >> 136 G4double &pMax ) const >> 137 { >> 138 >> 139 >> 140 G4double maxSide = ( fDx1 > fDx2 ? fDx1 : fDx2 ) ; >> 141 G4double maxRad = std::sqrt( maxSide*maxSide + fDy*fDy); >> 142 >> 143 if (!pTransform.IsRotated()) >> 144 { >> 145 // Special case handling for unrotated boxes >> 146 // Compute x/y/z mins and maxs respecting limits, with early returns >> 147 // if outside limits. Then switch() on pAxis >> 148 >> 149 G4double xoffset,xMin,xMax; >> 150 G4double yoffset,yMin,yMax; >> 151 G4double zoffset,zMin,zMax; >> 152 >> 153 xoffset = pTransform.NetTranslation().x() ; >> 154 xMin = xoffset - maxRad ; >> 155 xMax = xoffset + maxRad ; >> 156 >> 157 if (pVoxelLimit.IsXLimited()) >> 158 { >> 159 if ( xMin > pVoxelLimit.GetMaxXExtent()+kCarTolerance || >> 160 xMax < pVoxelLimit.GetMinXExtent()-kCarTolerance ) return false ; >> 161 else >> 162 { >> 163 if (xMin < pVoxelLimit.GetMinXExtent()) >> 164 { >> 165 xMin = pVoxelLimit.GetMinXExtent() ; >> 166 } >> 167 if (xMax > pVoxelLimit.GetMaxXExtent()) >> 168 { >> 169 xMax = pVoxelLimit.GetMaxXExtent() ; >> 170 } >> 171 } >> 172 } >> 173 yoffset = pTransform.NetTranslation().y() ; >> 174 yMin = yoffset - maxRad ; >> 175 yMax = yoffset + maxRad ; >> 176 >> 177 if (pVoxelLimit.IsYLimited()) >> 178 { >> 179 if ( yMin > pVoxelLimit.GetMaxYExtent()+kCarTolerance || >> 180 yMax < pVoxelLimit.GetMinYExtent()-kCarTolerance ) return false ; >> 181 else >> 182 { >> 183 if (yMin < pVoxelLimit.GetMinYExtent()) >> 184 { >> 185 yMin = pVoxelLimit.GetMinYExtent() ; >> 186 } >> 187 if (yMax > pVoxelLimit.GetMaxYExtent()) >> 188 { >> 189 yMax = pVoxelLimit.GetMaxYExtent() ; >> 190 } >> 191 } >> 192 } >> 193 zoffset = pTransform.NetTranslation().z() ; >> 194 zMin = zoffset - fDz ; >> 195 zMax = zoffset + fDz ; >> 196 >> 197 if (pVoxelLimit.IsZLimited()) >> 198 { >> 199 if ( zMin > pVoxelLimit.GetMaxZExtent()+kCarTolerance || >> 200 zMax < pVoxelLimit.GetMinZExtent()-kCarTolerance ) return false ; >> 201 else >> 202 { >> 203 if (zMin < pVoxelLimit.GetMinZExtent()) >> 204 { >> 205 zMin = pVoxelLimit.GetMinZExtent() ; >> 206 } >> 207 if (zMax > pVoxelLimit.GetMaxZExtent()) >> 208 { >> 209 zMax = pVoxelLimit.GetMaxZExtent() ; >> 210 } >> 211 } >> 212 } >> 213 switch (pAxis) >> 214 { >> 215 case kXAxis: >> 216 pMin = xMin ; >> 217 pMax = xMax ; >> 218 break ; >> 219 case kYAxis: >> 220 pMin=yMin; >> 221 pMax=yMax; >> 222 break; >> 223 case kZAxis: >> 224 pMin=zMin; >> 225 pMax=zMax; >> 226 break; >> 227 default: >> 228 break; >> 229 } >> 230 pMin -= kCarTolerance ; >> 231 pMax += kCarTolerance ; >> 232 >> 233 return true; >> 234 } >> 235 else // General rotated case - create and clip mesh to boundaries >> 236 { >> 237 G4bool existsAfterClip = false ; >> 238 G4ThreeVectorList* vertices ; >> 239 >> 240 pMin = +kInfinity ; >> 241 pMax = -kInfinity ; >> 242 >> 243 // Calculate rotated vertex coordinates >> 244 >> 245 vertices = CreateRotatedVertices(pTransform) ; >> 246 ClipCrossSection(vertices,0,pVoxelLimit,pAxis,pMin,pMax) ; >> 247 ClipCrossSection(vertices,4,pVoxelLimit,pAxis,pMin,pMax) ; >> 248 ClipBetweenSections(vertices,0,pVoxelLimit,pAxis,pMin,pMax) ; >> 249 >> 250 if (pVoxelLimit.IsLimited(pAxis) == false) >> 251 { >> 252 if ( pMin != kInfinity || pMax != -kInfinity ) >> 253 { >> 254 existsAfterClip = true ; >> 255 >> 256 // Add 2*tolerance to avoid precision troubles >> 257 >> 258 pMin -= kCarTolerance; >> 259 pMax += kCarTolerance; >> 260 } >> 261 } >> 262 else >> 263 { >> 264 G4ThreeVector clipCentre( >> 265 ( pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5, >> 266 ( pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5, >> 267 ( pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5); >> 268 >> 269 if ( pMin != kInfinity || pMax != -kInfinity ) >> 270 { >> 271 existsAfterClip = true ; >> 272 >> 273 >> 274 // Check to see if endpoints are in the solid >> 275 >> 276 clipCentre(pAxis) = pVoxelLimit.GetMinExtent(pAxis); >> 277 >> 278 if (Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside) >> 279 { >> 280 pMin = pVoxelLimit.GetMinExtent(pAxis); >> 281 } >> 282 else >> 283 { >> 284 pMin -= kCarTolerance; >> 285 } >> 286 clipCentre(pAxis) = pVoxelLimit.GetMaxExtent(pAxis); >> 287 >> 288 if (Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside) >> 289 { >> 290 pMax = pVoxelLimit.GetMaxExtent(pAxis); >> 291 } >> 292 else >> 293 { >> 294 pMax += kCarTolerance; >> 295 } >> 296 } >> 297 >> 298 // Check for case where completely enveloping clipping volume >> 299 // If point inside then we are confident that the solid completely >> 300 // envelopes the clipping volume. Hence set min/max extents according >> 301 // to clipping volume extents along the specified axis. >> 302 >> 303 else if (Inside(pTransform.Inverse().TransformPoint(clipCentre)) >> 304 != kOutside) >> 305 { >> 306 existsAfterClip = true ; >> 307 pMin = pVoxelLimit.GetMinExtent(pAxis) ; >> 308 pMax = pVoxelLimit.GetMaxExtent(pAxis) ; >> 309 } >> 310 } >> 311 delete vertices; >> 312 return existsAfterClip; >> 313 } >> 314 >> 315 >> 316 } >> 317 >> 318 G4ThreeVectorList* >> 319 G4TwistedTrap::CreateRotatedVertices(const G4AffineTransform& pTransform) const >> 320 { >> 321 G4ThreeVectorList* vertices = new G4ThreeVectorList(); >> 322 vertices->reserve(8); >> 323 >> 324 if (vertices) >> 325 { >> 326 >> 327 G4double maxSide = ( fDx1 > fDx2 ? fDx1 : fDx2 ) ; >> 328 G4double maxRad = std::sqrt( maxSide*maxSide + fDy*fDy); >> 329 >> 330 G4ThreeVector vertex0(-maxRad,-maxRad,-fDz) ; >> 331 G4ThreeVector vertex1(maxRad,-maxRad,-fDz) ; >> 332 G4ThreeVector vertex2(maxRad,maxRad,-fDz) ; >> 333 G4ThreeVector vertex3(-maxRad,maxRad,-fDz) ; >> 334 G4ThreeVector vertex4(-maxRad,-maxRad,fDz) ; >> 335 G4ThreeVector vertex5(maxRad,-maxRad,fDz) ; >> 336 G4ThreeVector vertex6(maxRad,maxRad,fDz) ; >> 337 G4ThreeVector vertex7(-maxRad,maxRad,fDz) ; >> 338 >> 339 vertices->push_back(pTransform.TransformPoint(vertex0)); >> 340 vertices->push_back(pTransform.TransformPoint(vertex1)); >> 341 vertices->push_back(pTransform.TransformPoint(vertex2)); >> 342 vertices->push_back(pTransform.TransformPoint(vertex3)); >> 343 vertices->push_back(pTransform.TransformPoint(vertex4)); >> 344 vertices->push_back(pTransform.TransformPoint(vertex5)); >> 345 vertices->push_back(pTransform.TransformPoint(vertex6)); >> 346 vertices->push_back(pTransform.TransformPoint(vertex7)); >> 347 } >> 348 else >> 349 { >> 350 DumpInfo(); >> 351 G4Exception("G4TwistedTrap::CreateRotatedVertices()", >> 352 "FatalError", FatalException, >> 353 "Error in allocation of vertices. Out of memory !"); >> 354 } >> 355 return vertices; >> 356 } 34 357 35 //============================================ 358 //===================================================================== 36 //* Constructors ----------------------------- << 359 //* Inside ------------------------------------------------------------ 37 360 38 G4TwistedTrap::G4TwistedTrap( const G4String& << 361 EInside G4TwistedTrap::Inside(const G4ThreeVector& p) const 39 G4double << 40 G4double << 41 G4double << 42 G4double << 43 G4double << 44 : G4VTwistedFaceted( pName, pPhiTwist,pDz,0. << 45 pDy, pDx1, pDx2, pDy, p << 46 { 362 { >> 363 >> 364 G4ThreeVector *tmpp; >> 365 EInside *tmpin; >> 366 if (fLastInside.p == p) { >> 367 return fLastInside.inside; >> 368 } else { >> 369 tmpp = const_cast<G4ThreeVector*>(&(fLastInside.p)); >> 370 tmpin = const_cast<EInside*>(&(fLastInside.inside)); >> 371 tmpp->set(p.x(), p.y(), p.z()); >> 372 } >> 373 >> 374 *tmpin = kOutside ; >> 375 >> 376 G4double phi = - p.z()/(2*fDz) * fPhiTwist ; // rotate the point to z=0 >> 377 G4double cphi = std::cos(phi) ; >> 378 G4double sphi = std::sin(phi) ; >> 379 G4double posx = p.x() * cphi - p.y() * sphi ; >> 380 G4double posy = p.x() * sphi + p.y() * cphi ; >> 381 G4double posz = p.z() ; >> 382 >> 383 G4double wpos = fDx2 + ( fDx1-fDx2)/2. - posy * (fDx1-fDx2)/(2*fDy) ; >> 384 >> 385 >> 386 #ifdef G4SPECSDEBUG >> 387 >> 388 G4cout << "inside called: p = " << p << G4endl ; >> 389 G4cout << "Trapezoid is 1/2*(x,y,z) : " << fDx1 << ", " << fDx2 << ", " << >> 390 fDy << ", " << fDz << G4endl ; >> 391 G4cout << "posx = " << posx << G4endl ; >> 392 G4cout << "posy = " << posy << G4endl ; >> 393 >> 394 G4cout << "wpos = " << wpos << G4endl ; >> 395 >> 396 #endif >> 397 >> 398 if ( std::fabs(posx) <= wpos - kCarTolerance*0.5 ) >> 399 { >> 400 if (std::fabs(posy) <= fDy - kCarTolerance*0.5 ) >> 401 { >> 402 if (std::fabs(posz) <= fDz - kCarTolerance*0.5 ) *tmpin = kInside ; >> 403 else if (std::fabs(posz) <= fDz + kCarTolerance*0.5 ) *tmpin = kSurface ; >> 404 } >> 405 else if (std::fabs(posy) <= fDy + kCarTolerance*0.5 ) >> 406 { >> 407 if (std::fabs(posz) <= fDz + kCarTolerance*0.5 ) *tmpin = kSurface ; >> 408 } >> 409 } >> 410 else if (std::fabs(posx) <= wpos + kCarTolerance*0.5 ) >> 411 { >> 412 if (std::fabs(posy) <= fDy + kCarTolerance*0.5 ) >> 413 { >> 414 if (std::fabs(posz) <= fDz + kCarTolerance*0.5) *tmpin = kSurface ; >> 415 } >> 416 } >> 417 >> 418 return fLastInside.inside; >> 419 47 } 420 } 48 421 49 G4TwistedTrap:: << 422 //===================================================================== 50 G4TwistedTrap(const G4String& pName, // N << 423 //* SurfaceNormal ----------------------------------------------------- 51 G4double pPhiTwist, // t << 424 52 G4double pDz, // h << 425 G4ThreeVector G4TwistedTrap::SurfaceNormal(const G4ThreeVector& p) const 53 G4double pTheta, // dire << 54 G4double pPhi, // defi << 55 G4double pDy1, // half << 56 G4double pDx1, // half << 57 G4double pDx2, // half << 58 G4double pDy2, // half << 59 G4double pDx3, // half << 60 G4double pDx4, // half << 61 G4double pAlph ) // tilt << 62 : G4VTwistedFaceted( pName, pPhiTwist, pDz, << 63 pPhi, pDy1, pDx1, pDx2, << 64 { 426 { >> 427 // >> 428 // return the normal unit vector to the Hyperbolical Surface at a point >> 429 // p on (or nearly on) the surface >> 430 // >> 431 // Which of the three or four surfaces are we closest to? >> 432 // >> 433 >> 434 >> 435 if (fLastNormal.p == p) { >> 436 >> 437 return fLastNormal.vec; >> 438 } >> 439 >> 440 G4ThreeVector *tmpp = const_cast<G4ThreeVector*>(&(fLastNormal.p)); >> 441 G4ThreeVector *tmpnormal = const_cast<G4ThreeVector*>(&(fLastNormal.vec)); >> 442 G4VSurface **tmpsurface = const_cast<G4VSurface**>(fLastNormal.surface); >> 443 tmpp->set(p.x(), p.y(), p.z()); >> 444 >> 445 G4double distance = kInfinity; >> 446 >> 447 G4VSurface *surfaces[6]; >> 448 >> 449 surfaces[0] = fSide0 ; >> 450 surfaces[1] = fSide90 ; >> 451 surfaces[2] = fSide180 ; >> 452 surfaces[3] = fSide270 ; >> 453 surfaces[4] = fLowerEndcap; >> 454 surfaces[5] = fUpperEndcap; >> 455 >> 456 G4ThreeVector xx; >> 457 G4ThreeVector bestxx; >> 458 G4int i; >> 459 G4int besti = -1; >> 460 for (i=0; i< 6; i++) { >> 461 G4double tmpdistance = surfaces[i]->DistanceTo(p, xx); >> 462 if (tmpdistance < distance) { >> 463 distance = tmpdistance; >> 464 bestxx = xx; >> 465 besti = i; >> 466 } >> 467 } >> 468 >> 469 tmpsurface[0] = surfaces[besti]; >> 470 *tmpnormal = tmpsurface[0]->GetNormal(bestxx, true); >> 471 >> 472 return fLastNormal.vec; 65 } 473 } 66 474 67 //============================================ 475 //===================================================================== 68 // Fake default constructor - sets only member << 476 //* DistanceToIn (p, v) ----------------------------------------------- 69 // for usage restri << 70 477 71 G4TwistedTrap::G4TwistedTrap( __void__& a ) << 478 G4double G4TwistedTrap::DistanceToIn (const G4ThreeVector& p, 72 : G4VTwistedFaceted(a) << 479 const G4ThreeVector& v ) const 73 { 480 { >> 481 >> 482 // DistanceToIn (p, v): >> 483 // Calculate distance to surface of shape from `outside' >> 484 // along with the v, allowing for tolerance. >> 485 // The function returns kInfinity if no intersection or >> 486 // just grazing within tolerance. >> 487 >> 488 // >> 489 // checking last value >> 490 // >> 491 >> 492 G4ThreeVector *tmpp; >> 493 G4ThreeVector *tmpv; >> 494 G4double *tmpdist; >> 495 if (fLastDistanceToInWithV.p == p && fLastDistanceToInWithV.vec == v) { >> 496 >> 497 >> 498 return fLastDistanceToIn.value; >> 499 } else { >> 500 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToInWithV.p)); >> 501 tmpv = const_cast<G4ThreeVector*>(&(fLastDistanceToInWithV.vec)); >> 502 tmpdist = const_cast<G4double*>(&(fLastDistanceToInWithV.value)); >> 503 tmpp->set(p.x(), p.y(), p.z()); >> 504 tmpv->set(v.x(), v.y(), v.z()); >> 505 } >> 506 >> 507 // >> 508 // Calculate DistanceToIn(p,v) >> 509 // >> 510 >> 511 EInside currentside = Inside(p); >> 512 >> 513 if (currentside == kInside) { >> 514 >> 515 } else if (currentside == kSurface) { >> 516 // particle is just on a boundary. >> 517 // if the particle is entering to the volume, return 0. >> 518 G4ThreeVector normal = SurfaceNormal(p); >> 519 if (normal*v < 0) { >> 520 >> 521 *tmpdist = 0; >> 522 return fLastDistanceToInWithV.value; >> 523 } >> 524 } >> 525 >> 526 // now, we can take smallest positive distance. >> 527 >> 528 // Initialize >> 529 G4double distance = kInfinity; >> 530 >> 531 // find intersections and choose nearest one. >> 532 G4VSurface *surfaces[6]; >> 533 >> 534 surfaces[0] = fSide0; >> 535 surfaces[1] = fSide90 ; >> 536 surfaces[2] = fSide180 ; >> 537 surfaces[3] = fSide270 ; >> 538 surfaces[4] = fLowerEndcap; >> 539 surfaces[5] = fUpperEndcap; >> 540 >> 541 G4ThreeVector xx; >> 542 G4ThreeVector bestxx; >> 543 G4int i; >> 544 G4int besti = -1; >> 545 for (i=0; i < 6 ; i++) { >> 546 >> 547 #ifdef G4SPECSDEBUG >> 548 G4cout << G4endl << "surface " << i << ": " << G4endl << G4endl ; >> 549 #endif >> 550 G4double tmpdistance = surfaces[i]->DistanceToIn(p, v, xx); >> 551 #ifdef G4SPECSDEBUG >> 552 G4cout << "Solid DistanceToIn : distance = " << tmpdistance << G4endl ; >> 553 G4cout << "intersection point = " << xx << G4endl ; >> 554 #endif >> 555 if (tmpdistance < distance) { >> 556 distance = tmpdistance; >> 557 bestxx = xx; >> 558 besti = i; >> 559 } >> 560 } >> 561 >> 562 #ifdef G4SPECSDEBUG >> 563 G4cout << "best distance = " << distance << G4endl ; >> 564 #endif >> 565 >> 566 *tmpdist = distance; >> 567 // timer.Stop(); >> 568 return fLastDistanceToInWithV.value; 74 } 569 } 75 570 >> 571 #ifdef DISTANCETOIN 76 //============================================ 572 //===================================================================== 77 //* Destructor ------------------------------- << 573 //* DistanceToIn (p) -------------------------------------------------- >> 574 >> 575 G4double G4TwistedTrap::DistanceToIn (const G4ThreeVector& p) const >> 576 { >> 577 // DistanceToIn(p): >> 578 // Calculate distance to surface of shape from `outside', >> 579 // allowing for tolerance >> 580 // >> 581 >> 582 // >> 583 // checking last value >> 584 // >> 585 >> 586 G4ThreeVector *tmpp; >> 587 G4double *tmpdist; >> 588 if (fLastDistanceToIn.p == p) { >> 589 return fLastDistanceToIn.value; >> 590 } else { >> 591 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToIn.p)); >> 592 tmpdist = const_cast<G4double*>(&(fLastDistanceToIn.value)); >> 593 tmpp->set(p.x(), p.y(), p.z()); >> 594 } >> 595 >> 596 >> 597 // >> 598 // Calculate DistanceToIn(p) >> 599 // >> 600 >> 601 EInside currentside = Inside(p); >> 602 >> 603 switch (currentside) { >> 604 >> 605 case (kInside) : { >> 606 } >> 607 >> 608 case (kSurface) : { >> 609 *tmpdist = 0.; >> 610 return fLastDistanceToIn.value; >> 611 } >> 612 >> 613 case (kOutside) : { >> 614 // Initialize >> 615 >> 616 G4double safex, safey, safez, safe = 0.0 ; >> 617 >> 618 G4double maxSide = ( fDx1 > fDx2 ? fDx1 : fDx2 ) ; >> 619 G4double maxRad = std::sqrt( maxSide*maxSide + fDy*fDy); >> 620 >> 621 G4cout << "maxRad = " << maxRad << G4endl ; >> 622 >> 623 safex = std::fabs(p.x()) - maxRad ; >> 624 safey = std::fabs(p.y()) - maxRad ; >> 625 safez = std::fabs(p.z()) - fDz ; >> 626 >> 627 if (safex > safe) safe = safex ; >> 628 if (safey > safe) safe = safey ; >> 629 if (safez > safe) safe = safez ; >> 630 >> 631 *tmpdist = safe; >> 632 return fLastDistanceToIn.value; >> 633 >> 634 } >> 635 >> 636 >> 637 default : { >> 638 G4Exception("G4TwistedTrap::DistanceToIn(p)", "InvalidCondition", >> 639 FatalException, "Unknown point location!"); >> 640 } >> 641 } // switch end >> 642 return kInfinity; >> 643 } >> 644 #else >> 645 >> 646 G4double G4TwistedTrap::DistanceToIn (const G4ThreeVector& p) const >> 647 { >> 648 // DistanceToIn(p): >> 649 // Calculate distance to surface of shape from `outside', >> 650 // allowing for tolerance >> 651 // >> 652 >> 653 // >> 654 // checking last value >> 655 // >> 656 >> 657 G4ThreeVector *tmpp; >> 658 G4double *tmpdist; >> 659 if (fLastDistanceToIn.p == p) { >> 660 >> 661 return fLastDistanceToIn.value; >> 662 } else { >> 663 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToIn.p)); >> 664 tmpdist = const_cast<G4double*>(&(fLastDistanceToIn.value)); >> 665 tmpp->set(p.x(), p.y(), p.z()); >> 666 } >> 667 >> 668 // >> 669 // Calculate DistanceToIn(p) >> 670 // >> 671 >> 672 EInside currentside = Inside(p); >> 673 >> 674 switch (currentside) { 78 675 79 G4TwistedTrap::~G4TwistedTrap() = default; << 676 case (kInside) : { >> 677 >> 678 } >> 679 >> 680 case (kSurface) : { >> 681 *tmpdist = 0.; >> 682 return fLastDistanceToIn.value; >> 683 } >> 684 >> 685 case (kOutside) : { >> 686 // Initialize >> 687 G4double distance = kInfinity; >> 688 >> 689 // find intersections and choose nearest one. >> 690 G4VSurface *surfaces[6]; >> 691 >> 692 surfaces[0] = fSide0; >> 693 surfaces[1] = fSide90 ; >> 694 surfaces[2] = fSide180 ; >> 695 surfaces[3] = fSide270 ; >> 696 surfaces[4] = fLowerEndcap; >> 697 surfaces[5] = fUpperEndcap; >> 698 >> 699 G4int i; >> 700 G4int besti = -1; >> 701 G4ThreeVector xx; >> 702 G4ThreeVector bestxx; >> 703 for (i=0; i< 6; i++) { >> 704 G4double tmpdistance = surfaces[i]->DistanceTo(p, xx); >> 705 if (tmpdistance < distance) { >> 706 distance = tmpdistance; >> 707 bestxx = xx; >> 708 besti = i; >> 709 } >> 710 } >> 711 >> 712 *tmpdist = distance; >> 713 return fLastDistanceToIn.value; >> 714 } >> 715 >> 716 default : { >> 717 G4Exception("G4TwistedTrap::DistanceToIn(p)", "InvalidCondition", >> 718 FatalException, "Unknown point location!"); >> 719 } >> 720 } // switch end >> 721 return kInfinity; >> 722 } >> 723 >> 724 #endif 80 725 81 //============================================ 726 //===================================================================== 82 //* Copy constructor ------------------------- << 727 //* DistanceToOut (p, v) ---------------------------------------------- 83 728 84 G4TwistedTrap::G4TwistedTrap(const G4TwistedTr << 729 G4double G4TwistedTrap::DistanceToOut( const G4ThreeVector& p, 85 : G4VTwistedFaceted(rhs) << 730 const G4ThreeVector& v, >> 731 const G4bool calcNorm, >> 732 G4bool *validNorm, G4ThreeVector *norm ) const 86 { 733 { 87 fpPolyhedron = GetPolyhedron(); << 734 // DistanceToOut (p, v): >> 735 // Calculate distance to surface of shape from `inside' >> 736 // along with the v, allowing for tolerance. >> 737 // The function returns kInfinity if no intersection or >> 738 // just grazing within tolerance. >> 739 // >> 740 >> 741 >> 742 // >> 743 // checking last value >> 744 // >> 745 >> 746 G4ThreeVector *tmpp; >> 747 G4ThreeVector *tmpv; >> 748 G4double *tmpdist; >> 749 if (fLastDistanceToOutWithV.p == p && fLastDistanceToOutWithV.vec == v ) { >> 750 return fLastDistanceToOutWithV.value; >> 751 } else { >> 752 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToOutWithV.p)); >> 753 tmpv = const_cast<G4ThreeVector*>(&(fLastDistanceToOutWithV.vec)); >> 754 tmpdist = const_cast<G4double*>(&(fLastDistanceToOutWithV.value)); >> 755 tmpp->set(p.x(), p.y(), p.z()); >> 756 tmpv->set(v.x(), v.y(), v.z()); >> 757 } >> 758 >> 759 // >> 760 // Calculate DistanceToOut(p,v) >> 761 // >> 762 >> 763 EInside currentside = Inside(p); >> 764 >> 765 if (currentside == kOutside) { >> 766 >> 767 } else if (currentside == kSurface) { >> 768 // particle is just on a boundary. >> 769 // if the particle is exiting from the volume, return 0. >> 770 G4ThreeVector normal = SurfaceNormal(p); >> 771 G4VSurface *blockedsurface = fLastNormal.surface[0]; >> 772 if (normal*v > 0) { >> 773 if (calcNorm) { >> 774 *norm = (blockedsurface->GetNormal(p, true)); >> 775 *validNorm = blockedsurface->IsValidNorm(); >> 776 } >> 777 *tmpdist = 0.; >> 778 // timer.Stop(); >> 779 return fLastDistanceToOutWithV.value; >> 780 } >> 781 } >> 782 >> 783 // now, we can take smallest positive distance. >> 784 >> 785 // Initialize >> 786 G4double distance = kInfinity; >> 787 >> 788 // find intersections and choose nearest one. >> 789 G4VSurface *surfaces[6]; >> 790 >> 791 surfaces[0] = fSide0; >> 792 surfaces[1] = fSide90 ; >> 793 surfaces[2] = fSide180 ; >> 794 surfaces[3] = fSide270 ; >> 795 surfaces[4] = fLowerEndcap; >> 796 surfaces[5] = fUpperEndcap; >> 797 >> 798 G4int i; >> 799 G4int besti = -1; >> 800 G4ThreeVector xx; >> 801 G4ThreeVector bestxx; >> 802 for (i=0; i< 6 ; i++) { >> 803 G4double tmpdistance = surfaces[i]->DistanceToOut(p, v, xx); >> 804 if (tmpdistance < distance) { >> 805 distance = tmpdistance; >> 806 bestxx = xx; >> 807 besti = i; >> 808 } >> 809 } >> 810 >> 811 if (calcNorm) { >> 812 if (besti != -1) { >> 813 *norm = (surfaces[besti]->GetNormal(p, true)); >> 814 *validNorm = surfaces[besti]->IsValidNorm(); >> 815 } >> 816 } >> 817 >> 818 *tmpdist = distance; >> 819 // timer.Stop(); >> 820 return fLastDistanceToOutWithV.value; 88 } 821 } 89 822 >> 823 #ifdef DISTANCETOIN 90 //============================================ 824 //===================================================================== 91 //* Assignment operator ---------------------- << 825 //* DistanceToOut (p) ---------------------------------------------- 92 826 93 G4TwistedTrap& G4TwistedTrap::operator = (cons << 827 G4double G4TwistedTrap::DistanceToOut( const G4ThreeVector& p ) const 94 { 828 { 95 // Check assignment to self << 829 // DistanceToOut(p): >> 830 // Calculate distance to surface of shape from `inside', >> 831 // allowing for tolerance >> 832 // >> 833 >> 834 // >> 835 // checking last value 96 // 836 // 97 if (this == &rhs) { return *this; } << 98 837 99 // Copy base class data << 838 G4ThreeVector *tmpp; >> 839 G4double *tmpdist; >> 840 if (fLastDistanceToOut.p == p) { >> 841 return fLastDistanceToOut.value; >> 842 } else { >> 843 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToOut.p)); >> 844 tmpdist = const_cast<G4double*>(&(fLastDistanceToOut.value)); >> 845 tmpp->set(p.x(), p.y(), p.z()); >> 846 } >> 847 >> 848 // >> 849 // Calculate DistanceToOut(p) >> 850 // >> 851 >> 852 EInside currentside = Inside(p); >> 853 >> 854 switch (currentside) { >> 855 case (kOutside) : { >> 856 >> 857 } >> 858 case (kSurface) : { >> 859 *tmpdist = 0.; >> 860 return fLastDistanceToOut.value; >> 861 } >> 862 >> 863 case (kInside) : { >> 864 // Initialize >> 865 >> 866 G4double rminX = ( fDx1 < fDx2 ? fDx1 : fDx2 ) ; >> 867 G4double rmin = ( rminX < fDy ? rminX : fDy ) ; >> 868 >> 869 G4double safx1,safx2,safy1,safy2,safz1,safz2,safe=0.0; >> 870 >> 871 safx1 = rmin - p.x() ; >> 872 safx2 = rmin + p.x() ; >> 873 safy1 = rmin - p.y() ; >> 874 safy2 = rmin + p.y() ; >> 875 safz1 = fDz - p.z() ; >> 876 safz2 = fDz + p.z() ; >> 877 >> 878 // shortest Dist to any boundary now MIN(safx1,safx2,safy1..) >> 879 >> 880 if (safx2 < safx1) safe = safx2 ; >> 881 else safe = safx1 ; >> 882 if (safy1 < safe) safe = safy1 ; >> 883 if (safy2 < safe) safe = safy2 ; >> 884 if (safz1 < safe) safe = safz1 ; >> 885 if (safz2 < safe) safe = safz2 ; >> 886 >> 887 if (safe < 0) safe = 0 ; >> 888 >> 889 *tmpdist = safe; >> 890 >> 891 return fLastDistanceToOut.value; >> 892 } >> 893 >> 894 default : { >> 895 G4Exception("G4TwistedTrap::DistanceToOut(p)", "InvalidCondition", >> 896 FatalException, "Unknown point location!"); >> 897 } >> 898 >> 899 } // switch end >> 900 >> 901 return 0; >> 902 >> 903 } >> 904 >> 905 #else >> 906 >> 907 G4double G4TwistedTrap::DistanceToOut( const G4ThreeVector& p ) const >> 908 { >> 909 // DistanceToOut(p): >> 910 // Calculate distance to surface of shape from `inside', >> 911 // allowing for tolerance >> 912 // >> 913 >> 914 100 // 915 // 101 G4VTwistedFaceted::operator=(rhs); << 916 // checking last value 102 fpPolyhedron = GetPolyhedron(); << 917 // >> 918 >> 919 G4ThreeVector *tmpp; >> 920 G4double *tmpdist; >> 921 if (fLastDistanceToOut.p == p) { >> 922 return fLastDistanceToOut.value; >> 923 } else { >> 924 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToOut.p)); >> 925 tmpdist = const_cast<G4double*>(&(fLastDistanceToOut.value)); >> 926 tmpp->set(p.x(), p.y(), p.z()); >> 927 } >> 928 >> 929 // >> 930 // Calculate DistanceToOut(p) >> 931 // >> 932 >> 933 EInside currentside = Inside(p); 103 934 104 return *this; << 935 switch (currentside) { >> 936 case (kOutside) : { >> 937 >> 938 } >> 939 case (kSurface) : { >> 940 *tmpdist = 0.; >> 941 return fLastDistanceToOut.value; >> 942 } >> 943 >> 944 case (kInside) : { >> 945 // Initialize >> 946 G4double distance = kInfinity; >> 947 >> 948 // find intersections and choose nearest one. >> 949 G4VSurface *surfaces[6]; >> 950 >> 951 surfaces[0] = fSide0; >> 952 surfaces[1] = fSide90 ; >> 953 surfaces[2] = fSide180 ; >> 954 surfaces[3] = fSide270 ; >> 955 surfaces[4] = fLowerEndcap; >> 956 surfaces[5] = fUpperEndcap; >> 957 >> 958 G4int i; >> 959 G4int besti = -1; >> 960 G4ThreeVector xx; >> 961 G4ThreeVector bestxx; >> 962 for (i=0; i< 6; i++) { >> 963 G4double tmpdistance = surfaces[i]->DistanceTo(p, xx); >> 964 if (tmpdistance < distance) { >> 965 distance = tmpdistance; >> 966 bestxx = xx; >> 967 besti = i; >> 968 } >> 969 } >> 970 >> 971 >> 972 *tmpdist = distance; >> 973 >> 974 return fLastDistanceToOut.value; >> 975 } >> 976 >> 977 default : { >> 978 G4Exception("G4TwistedTrap::DistanceToOut(p)", "InvalidCondition", >> 979 FatalException, "Unknown point location!"); >> 980 } >> 981 } // switch end >> 982 return 0; 105 } 983 } 106 984 >> 985 >> 986 #endif >> 987 107 //============================================ 988 //===================================================================== 108 //* StreamInfo ------------------------------- 989 //* StreamInfo -------------------------------------------------------- 109 990 110 std::ostream& G4TwistedTrap::StreamInfo(std::o 991 std::ostream& G4TwistedTrap::StreamInfo(std::ostream& os) const 111 { 992 { 112 // 993 // 113 // Stream object contents to an output strea 994 // Stream object contents to an output stream 114 // 995 // 115 os << "------------------------------------- 996 os << "-----------------------------------------------------------\n" 116 << " *** Dump for solid - " << GetName 997 << " *** Dump for solid - " << GetName() << " ***\n" 117 << " ================================= 998 << " ===================================================\n" 118 << " Solid type: G4TwistedTrap\n" 999 << " Solid type: G4TwistedTrap\n" 119 << " Parameters: \n" 1000 << " Parameters: \n" 120 << " Twist angle = " << GetPhi << 1001 << " Twist angle : " << fPhiTwist/degree << " degrees \n" 121 << G4endl << 122 << " Polar Angle Theta = " << GetPol << 123 << G4endl << 124 << " Azimuthal Angle Phi = " << GetAzi << 125 << G4endl << 126 << " pDy1 = " << GetY1HalfLength()/cm << 127 << " pDx1 = " << GetX1HalfLength()/cm << 128 << " pDx2 = " << GetX2HalfLength()/cm << 129 << " pDy2 = " << GetY2HalfLength()/cm << 130 << " pDx3 = " << GetX3HalfLength()/cm << 131 << " pDx4 = " << GetX4HalfLength()/cm << 132 << " pDz = " << GetZHalfLength()/cm < << 133 << " Tilt Angle Alpha = " << GetTil << 134 << G4endl << 135 << "------------------------------------- 1002 << "-----------------------------------------------------------\n"; 136 1003 137 return os; 1004 return os; 138 } 1005 } 139 1006 >> 1007 140 //============================================ 1008 //===================================================================== 141 //* GetEntityType ---------------------------- << 1009 //* DiscribeYourselfTo ------------------------------------------------ 142 1010 143 G4GeometryType G4TwistedTrap::GetEntityType() << 1011 void G4TwistedTrap::DescribeYourselfTo (G4VGraphicsScene& scene) const >> 1012 { >> 1013 scene.AddThis (*this); >> 1014 } >> 1015 >> 1016 //===================================================================== >> 1017 //* GetExtent --------------------------------------------------------- >> 1018 >> 1019 G4VisExtent G4TwistedTrap::GetExtent() const 144 { 1020 { 145 return {"G4TwistedTrap"}; << 1021 >> 1022 G4double maxSide = ( fDx1 > fDx2 ? fDx1 : fDx2 ) ; >> 1023 G4double maxRad = std::sqrt( maxSide*maxSide + fDy*fDy); >> 1024 >> 1025 return G4VisExtent(-maxRad, maxRad , >> 1026 -maxRad, maxRad , >> 1027 -fDz, fDz ); 146 } 1028 } 147 1029 148 //============================================ 1030 //===================================================================== 149 //* Clone ------------------------------------ << 1031 //* CreatePolyhedron -------------------------------------------------- 150 1032 151 G4VSolid* G4TwistedTrap::Clone() const << 1033 G4Polyhedron* G4TwistedTrap::CreatePolyhedron () const >> 1034 { >> 1035 // return new G4PolyhedronHype (fRMin, fRMax, fDz, fSPhi, fDPhi); >> 1036 return 0; >> 1037 } >> 1038 >> 1039 //===================================================================== >> 1040 //* CreateNUBS -------------------------------------------------------- >> 1041 >> 1042 G4NURBS* G4TwistedTrap::CreateNURBS () const >> 1043 { >> 1044 G4double maxSide = ( fDx1 > fDx2 ? fDx1 : fDx2 ) ; >> 1045 G4double maxRad = std::sqrt( maxSide*maxSide + fDy*fDy); >> 1046 >> 1047 return new G4NURBStube(maxRad, maxRad, fDz); >> 1048 // Tube for now!!! >> 1049 } >> 1050 >> 1051 //===================================================================== >> 1052 //* CreateSurfaces ---------------------------------------------------- >> 1053 >> 1054 void G4TwistedTrap::CreateSurfaces() >> 1055 { >> 1056 >> 1057 // create 6 surfaces of TwistedTub. >> 1058 >> 1059 >> 1060 fSide0 = new G4TwistedTrapSide("0deg", fPhiTwist, fDx1, fDx2, fDy, fDz, 0.*deg ) ; >> 1061 fSide180 = new G4TwistedTrapSide("180deg", fPhiTwist, fDx2, fDx1, fDy, fDz, 180.*deg) ; >> 1062 >> 1063 fSide90 = new G4TwistedBoxSide("90deg", fPhiTwist, fDy, fDx2, fDz, 90.*deg) ; >> 1064 fSide270 = new G4TwistedBoxSide("270deg", fPhiTwist, fDy, fDx1, fDz, 270.*deg) ; >> 1065 >> 1066 fUpperEndcap = new G4FlatTrapSide("UpperCap",fPhiTwist, fDx1, fDx2, fDy, fDz, 1 ) ; >> 1067 fLowerEndcap = new G4FlatTrapSide("LowerCap",fPhiTwist, fDx1, fDx2, fDy, fDz, -1 ) ; >> 1068 >> 1069 // Set neighbour surfaces >> 1070 >> 1071 fSide0->SetNeighbours( fSide270 , fLowerEndcap , fSide90 , fUpperEndcap ); >> 1072 fSide90->SetNeighbours( fSide0 , fLowerEndcap , fSide180 , fUpperEndcap ); >> 1073 fSide180->SetNeighbours(fSide90 , fLowerEndcap , fSide270 , fUpperEndcap ); >> 1074 fSide270->SetNeighbours(fSide180 , fLowerEndcap , fSide0 , fUpperEndcap ); >> 1075 fUpperEndcap->SetNeighbours( fSide180, fSide270 , fSide0 , fSide90 ); >> 1076 fLowerEndcap->SetNeighbours( fSide180, fSide270 , fSide0 , fSide90 ); >> 1077 >> 1078 } >> 1079 >> 1080 >> 1081 //===================================================================== >> 1082 //* GetEntityType ----------------------------------------------------- >> 1083 >> 1084 G4GeometryType G4TwistedTrap::GetEntityType() const 152 { 1085 { 153 return new G4TwistedTrap(*this); << 1086 return G4String("G4TwistedTrap"); 154 } 1087 } 155 1088