| 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-cand-05 $
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