Geant4 Cross Reference |
1 // 2 // ******************************************************************** 3 // * License and Disclaimer * 4 // * * 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. * 10 // * * 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitation of liability. * 17 // * * 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************************************** 25 // 26 // Implementation of G4Polycone, a CSG polycone 27 // 28 // Author: David C. Williams (davidw@scipp.ucsc.edu) 29 // -------------------------------------------------------------------- 30 31 #include "G4Polycone.hh" 32 33 #if !defined(G4GEOM_USE_UPOLYCONE) 34 35 #include "G4PolyconeSide.hh" 36 #include "G4PolyPhiFace.hh" 37 38 #include "G4GeomTools.hh" 39 #include "G4VoxelLimits.hh" 40 #include "G4AffineTransform.hh" 41 #include "G4BoundingEnvelope.hh" 42 43 #include "G4QuickRand.hh" 44 45 #include "G4EnclosingCylinder.hh" 46 #include "G4ReduciblePolygon.hh" 47 #include "G4VPVParameterisation.hh" 48 49 namespace 50 { 51 G4Mutex surface_elementsMutex = G4MUTEX_INITIALIZER; 52 } 53 54 using namespace CLHEP; 55 56 // Constructor (GEANT3 style parameters) 57 // 58 G4Polycone::G4Polycone( const G4String& name, 59 G4double phiStart, 60 G4double phiTotal, 61 G4int numZPlanes, 62 const G4double zPlane[], 63 const G4double rInner[], 64 const G4double rOuter[] ) 65 : G4VCSGfaceted( name ) 66 { 67 // 68 // Some historical ugliness 69 // 70 original_parameters = new G4PolyconeHistorical(); 71 72 original_parameters->Start_angle = phiStart; 73 original_parameters->Opening_angle = phiTotal; 74 original_parameters->Num_z_planes = numZPlanes; 75 original_parameters->Z_values = new G4double[numZPlanes]; 76 original_parameters->Rmin = new G4double[numZPlanes]; 77 original_parameters->Rmax = new G4double[numZPlanes]; 78 79 for (G4int i=0; i<numZPlanes; ++i) 80 { 81 if(rInner[i]>rOuter[i]) 82 { 83 DumpInfo(); 84 std::ostringstream message; 85 message << "Cannot create a Polycone with rInner > rOuter for the same Z" 86 << G4endl 87 << " rInner > rOuter for the same Z !" << G4endl 88 << " rMin[" << i << "] = " << rInner[i] 89 << " -- rMax[" << i << "] = " << rOuter[i]; 90 G4Exception("G4Polycone::G4Polycone()", "GeomSolids0002", 91 FatalErrorInArgument, message); 92 } 93 if (( i < numZPlanes-1) && ( zPlane[i] == zPlane[i+1] )) 94 { 95 if( (rInner[i] > rOuter[i+1]) 96 ||(rInner[i+1] > rOuter[i]) ) 97 { 98 DumpInfo(); 99 std::ostringstream message; 100 message << "Cannot create a Polycone with no contiguous segments." 101 << G4endl 102 << " Segments are not contiguous !" << G4endl 103 << " rMin[" << i << "] = " << rInner[i] 104 << " -- rMax[" << i+1 << "] = " << rOuter[i+1] << G4endl 105 << " rMin[" << i+1 << "] = " << rInner[i+1] 106 << " -- rMax[" << i << "] = " << rOuter[i]; 107 G4Exception("G4Polycone::G4Polycone()", "GeomSolids0002", 108 FatalErrorInArgument, message); 109 } 110 } 111 original_parameters->Z_values[i] = zPlane[i]; 112 original_parameters->Rmin[i] = rInner[i]; 113 original_parameters->Rmax[i] = rOuter[i]; 114 } 115 116 // 117 // Build RZ polygon using special PCON/PGON GEANT3 constructor 118 // 119 auto rz = new G4ReduciblePolygon( rInner, rOuter, zPlane, numZPlanes ); 120 121 // 122 // Do the real work 123 // 124 Create( phiStart, phiTotal, rz ); 125 126 delete rz; 127 } 128 129 // Constructor (generic parameters) 130 // 131 G4Polycone::G4Polycone( const G4String& name, 132 G4double phiStart, 133 G4double phiTotal, 134 G4int numRZ, 135 const G4double r[], 136 const G4double z[] ) 137 : G4VCSGfaceted( name ) 138 { 139 140 auto rz = new G4ReduciblePolygon( r, z, numRZ ); 141 142 Create( phiStart, phiTotal, rz ); 143 144 // Set original_parameters struct for consistency 145 // 146 147 G4bool convertible = SetOriginalParameters(rz); 148 149 if(!convertible) 150 { 151 std::ostringstream message; 152 message << "Polycone " << GetName() << "cannot be converted" << G4endl 153 << "to Polycone with (Rmin,Rmaz,Z) parameters!"; 154 G4Exception("G4Polycone::G4Polycone()", "GeomSolids0002", 155 FatalException, message, "Use G4GenericPolycone instead!"); 156 } 157 else 158 { 159 G4cout << "INFO: Converting polycone " << GetName() << G4endl 160 << "to optimized polycone with (Rmin,Rmaz,Z) parameters !" 161 << G4endl; 162 } 163 delete rz; 164 } 165 166 // Create 167 // 168 // Generic create routine, called by each constructor after 169 // conversion of arguments 170 // 171 void G4Polycone::Create( G4double phiStart, 172 G4double phiTotal, 173 G4ReduciblePolygon* rz ) 174 { 175 // 176 // Perform checks of rz values 177 // 178 if (rz->Amin() < 0.0) 179 { 180 std::ostringstream message; 181 message << "Illegal input parameters - " << GetName() << G4endl 182 << " All R values must be >= 0 !"; 183 G4Exception("G4Polycone::Create()", "GeomSolids0002", 184 FatalErrorInArgument, message); 185 } 186 187 G4double rzArea = rz->Area(); 188 if (rzArea < -kCarTolerance) 189 { 190 rz->ReverseOrder(); 191 } 192 else if (rzArea < kCarTolerance) 193 { 194 std::ostringstream message; 195 message << "Illegal input parameters - " << GetName() << G4endl 196 << " R/Z cross section is zero or near zero: " << rzArea; 197 G4Exception("G4Polycone::Create()", "GeomSolids0002", 198 FatalErrorInArgument, message); 199 } 200 201 if ( (!rz->RemoveDuplicateVertices( kCarTolerance )) 202 || (!rz->RemoveRedundantVertices( kCarTolerance )) ) 203 { 204 std::ostringstream message; 205 message << "Illegal input parameters - " << GetName() << G4endl 206 << " Too few unique R/Z values !"; 207 G4Exception("G4Polycone::Create()", "GeomSolids0002", 208 FatalErrorInArgument, message); 209 } 210 211 if (rz->CrossesItself(1/kInfinity)) 212 { 213 std::ostringstream message; 214 message << "Illegal input parameters - " << GetName() << G4endl 215 << " R/Z segments cross !"; 216 G4Exception("G4Polycone::Create()", "GeomSolids0002", 217 FatalErrorInArgument, message); 218 } 219 220 numCorner = rz->NumVertices(); 221 222 startPhi = phiStart; 223 while( startPhi < 0. ) // Loop checking, 13.08.2015, G.Cosmo 224 startPhi += twopi; 225 // 226 // Phi opening? Account for some possible roundoff, and interpret 227 // nonsense value as representing no phi opening 228 // 229 if ( (phiTotal <= 0) || (phiTotal > twopi*(1-DBL_EPSILON)) ) 230 { 231 phiIsOpen = false; 232 startPhi = 0.; 233 endPhi = twopi; 234 } 235 else 236 { 237 phiIsOpen = true; 238 endPhi = startPhi + phiTotal; 239 } 240 241 // 242 // Allocate corner array. 243 // 244 corners = new G4PolyconeSideRZ[numCorner]; 245 246 // 247 // Copy corners 248 // 249 G4ReduciblePolygonIterator iterRZ(rz); 250 251 G4PolyconeSideRZ *next = corners; 252 iterRZ.Begin(); 253 do // Loop checking, 13.08.2015, G.Cosmo 254 { 255 next->r = iterRZ.GetA(); 256 next->z = iterRZ.GetB(); 257 } while( ++next, iterRZ.Next() ); 258 259 // 260 // Allocate face pointer array 261 // 262 numFace = phiIsOpen ? numCorner+2 : numCorner; 263 faces = new G4VCSGface*[numFace]; 264 265 // 266 // Construct conical faces 267 // 268 // But! Don't construct a face if both points are at zero radius! 269 // 270 G4PolyconeSideRZ* corner = corners, 271 * prev = corners + numCorner-1, 272 * nextNext; 273 G4VCSGface **face = faces; 274 do // Loop checking, 13.08.2015, G.Cosmo 275 { 276 next = corner+1; 277 if (next >= corners+numCorner) next = corners; 278 nextNext = next+1; 279 if (nextNext >= corners+numCorner) nextNext = corners; 280 281 if (corner->r < 1/kInfinity && next->r < 1/kInfinity) continue; 282 283 // 284 // We must decide here if we can dare declare one of our faces 285 // as having a "valid" normal (i.e. allBehind = true). This 286 // is never possible if the face faces "inward" in r. 287 // 288 G4bool allBehind; 289 if (corner->z > next->z) 290 { 291 allBehind = false; 292 } 293 else 294 { 295 // 296 // Otherwise, it is only true if the line passing 297 // through the two points of the segment do not 298 // split the r/z cross section 299 // 300 allBehind = !rz->BisectedBy( corner->r, corner->z, 301 next->r, next->z, kCarTolerance ); 302 } 303 304 *face++ = new G4PolyconeSide( prev, corner, next, nextNext, 305 startPhi, endPhi-startPhi, phiIsOpen, allBehind ); 306 } while( prev=corner, corner=next, corner > corners ); 307 308 if (phiIsOpen) 309 { 310 // 311 // Construct phi open edges 312 // 313 *face++ = new G4PolyPhiFace( rz, startPhi, 0, endPhi ); 314 *face++ = new G4PolyPhiFace( rz, endPhi, 0, startPhi ); 315 } 316 317 // 318 // We might have dropped a face or two: recalculate numFace 319 // 320 numFace = (G4int)(face-faces); 321 322 // 323 // Make enclosingCylinder 324 // 325 enclosingCylinder = 326 new G4EnclosingCylinder( rz, phiIsOpen, phiStart, phiTotal ); 327 } 328 329 // Fake default constructor - sets only member data and allocates memory 330 // for usage restricted to object persistency. 331 // 332 G4Polycone::G4Polycone( __void__& a ) 333 : G4VCSGfaceted(a), startPhi(0.), endPhi(0.), numCorner(0) 334 { 335 } 336 337 338 // 339 // Destructor 340 // 341 G4Polycone::~G4Polycone() 342 { 343 delete [] corners; 344 delete original_parameters; 345 delete enclosingCylinder; 346 delete fElements; 347 delete fpPolyhedron; 348 corners = nullptr; 349 original_parameters = nullptr; 350 enclosingCylinder = nullptr; 351 fElements = nullptr; 352 fpPolyhedron = nullptr; 353 } 354 355 // Copy constructor 356 // 357 G4Polycone::G4Polycone( const G4Polycone& source ) 358 : G4VCSGfaceted( source ) 359 { 360 CopyStuff( source ); 361 } 362 363 // Assignment operator 364 // 365 G4Polycone &G4Polycone::operator=( const G4Polycone& source ) 366 { 367 if (this == &source) return *this; 368 369 G4VCSGfaceted::operator=( source ); 370 371 delete [] corners; 372 delete original_parameters; 373 374 delete enclosingCylinder; 375 376 CopyStuff( source ); 377 378 return *this; 379 } 380 381 // CopyStuff 382 // 383 void G4Polycone::CopyStuff( const G4Polycone& source ) 384 { 385 // 386 // Simple stuff 387 // 388 startPhi = source.startPhi; 389 endPhi = source.endPhi; 390 phiIsOpen = source.phiIsOpen; 391 numCorner = source.numCorner; 392 393 // 394 // The corner array 395 // 396 corners = new G4PolyconeSideRZ[numCorner]; 397 398 G4PolyconeSideRZ* corn = corners, 399 * sourceCorn = source.corners; 400 do // Loop checking, 13.08.2015, G.Cosmo 401 { 402 *corn = *sourceCorn; 403 } while( ++sourceCorn, ++corn < corners+numCorner ); 404 405 // 406 // Original parameters 407 // 408 if (source.original_parameters != nullptr) 409 { 410 original_parameters = 411 new G4PolyconeHistorical( *source.original_parameters ); 412 } 413 414 // 415 // Enclosing cylinder 416 // 417 enclosingCylinder = new G4EnclosingCylinder( *source.enclosingCylinder ); 418 419 // 420 // Surface elements 421 // 422 delete fElements; 423 fElements = nullptr; 424 425 // 426 // Polyhedron 427 // 428 fRebuildPolyhedron = false; 429 delete fpPolyhedron; 430 fpPolyhedron = nullptr; 431 } 432 433 // Reset 434 // 435 G4bool G4Polycone::Reset() 436 { 437 // 438 // Clear old setup 439 // 440 G4VCSGfaceted::DeleteStuff(); 441 delete [] corners; 442 delete enclosingCylinder; 443 delete fElements; 444 corners = nullptr; 445 fElements = nullptr; 446 enclosingCylinder = nullptr; 447 448 // 449 // Rebuild polycone 450 // 451 auto rz = new G4ReduciblePolygon( original_parameters->Rmin, 452 original_parameters->Rmax, 453 original_parameters->Z_values, 454 original_parameters->Num_z_planes ); 455 Create( original_parameters->Start_angle, 456 original_parameters->Opening_angle, rz ); 457 delete rz; 458 459 return false; 460 } 461 462 // Inside 463 // 464 // This is an override of G4VCSGfaceted::Inside, created in order 465 // to speed things up by first checking with G4EnclosingCylinder. 466 // 467 EInside G4Polycone::Inside( const G4ThreeVector& p ) const 468 { 469 // 470 // Quick test 471 // 472 if (enclosingCylinder->MustBeOutside(p)) return kOutside; 473 474 // 475 // Long answer 476 // 477 return G4VCSGfaceted::Inside(p); 478 } 479 480 // DistanceToIn 481 // 482 // This is an override of G4VCSGfaceted::Inside, created in order 483 // to speed things up by first checking with G4EnclosingCylinder. 484 // 485 G4double G4Polycone::DistanceToIn( const G4ThreeVector& p, 486 const G4ThreeVector& v ) const 487 { 488 // 489 // Quick test 490 // 491 if (enclosingCylinder->ShouldMiss(p,v)) 492 return kInfinity; 493 494 // 495 // Long answer 496 // 497 return G4VCSGfaceted::DistanceToIn( p, v ); 498 } 499 500 // DistanceToIn 501 // 502 G4double G4Polycone::DistanceToIn( const G4ThreeVector& p ) const 503 { 504 return G4VCSGfaceted::DistanceToIn(p); 505 } 506 507 // Get bounding box 508 // 509 void G4Polycone::BoundingLimits(G4ThreeVector& pMin, 510 G4ThreeVector& pMax) const 511 { 512 G4double rmin = kInfinity, rmax = -kInfinity; 513 G4double zmin = kInfinity, zmax = -kInfinity; 514 515 for (G4int i=0; i<GetNumRZCorner(); ++i) 516 { 517 G4PolyconeSideRZ corner = GetCorner(i); 518 if (corner.r < rmin) rmin = corner.r; 519 if (corner.r > rmax) rmax = corner.r; 520 if (corner.z < zmin) zmin = corner.z; 521 if (corner.z > zmax) zmax = corner.z; 522 } 523 524 if (IsOpen()) 525 { 526 G4TwoVector vmin,vmax; 527 G4GeomTools::DiskExtent(rmin,rmax, 528 GetSinStartPhi(),GetCosStartPhi(), 529 GetSinEndPhi(),GetCosEndPhi(), 530 vmin,vmax); 531 pMin.set(vmin.x(),vmin.y(),zmin); 532 pMax.set(vmax.x(),vmax.y(),zmax); 533 } 534 else 535 { 536 pMin.set(-rmax,-rmax, zmin); 537 pMax.set( rmax, rmax, zmax); 538 } 539 540 // Check correctness of the bounding box 541 // 542 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z()) 543 { 544 std::ostringstream message; 545 message << "Bad bounding box (min >= max) for solid: " 546 << GetName() << " !" 547 << "\npMin = " << pMin 548 << "\npMax = " << pMax; 549 G4Exception("G4Polycone::BoundingLimits()", "GeomMgt0001", 550 JustWarning, message); 551 DumpInfo(); 552 } 553 } 554 555 // Calculate extent under transform and specified limit 556 // 557 G4bool G4Polycone::CalculateExtent(const EAxis pAxis, 558 const G4VoxelLimits& pVoxelLimit, 559 const G4AffineTransform& pTransform, 560 G4double& pMin, G4double& pMax) const 561 { 562 G4ThreeVector bmin, bmax; 563 G4bool exist; 564 565 // Check bounding box (bbox) 566 // 567 BoundingLimits(bmin,bmax); 568 G4BoundingEnvelope bbox(bmin,bmax); 569 #ifdef G4BBOX_EXTENT 570 return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax); 571 #endif 572 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax)) 573 { 574 return exist = pMin < pMax; 575 } 576 577 // To find the extent, RZ contour of the polycone is subdivided 578 // in triangles. The extent is calculated as cumulative extent of 579 // all sub-polycones formed by rotation of triangles around Z 580 // 581 G4TwoVectorList contourRZ; 582 G4TwoVectorList triangles; 583 std::vector<G4int> iout; 584 G4double eminlim = pVoxelLimit.GetMinExtent(pAxis); 585 G4double emaxlim = pVoxelLimit.GetMaxExtent(pAxis); 586 587 // get RZ contour, ensure anticlockwise order of corners 588 for (G4int i=0; i<GetNumRZCorner(); ++i) 589 { 590 G4PolyconeSideRZ corner = GetCorner(i); 591 contourRZ.emplace_back(corner.r,corner.z); 592 } 593 G4GeomTools::RemoveRedundantVertices(contourRZ,iout,2*kCarTolerance); 594 G4double area = G4GeomTools::PolygonArea(contourRZ); 595 if (area < 0.) std::reverse(contourRZ.begin(),contourRZ.end()); 596 597 // triangulate RZ countour 598 if (!G4GeomTools::TriangulatePolygon(contourRZ,triangles)) 599 { 600 std::ostringstream message; 601 message << "Triangulation of RZ contour has failed for solid: " 602 << GetName() << " !" 603 << "\nExtent has been calculated using boundary box"; 604 G4Exception("G4Polycone::CalculateExtent()", 605 "GeomMgt1002", JustWarning, message); 606 return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax); 607 } 608 609 // set trigonometric values 610 const G4int NSTEPS = 24; // number of steps for whole circle 611 G4double astep = twopi/NSTEPS; // max angle for one step 612 613 G4double sphi = GetStartPhi(); 614 G4double ephi = GetEndPhi(); 615 G4double dphi = IsOpen() ? ephi-sphi : twopi; 616 G4int ksteps = (dphi <= astep) ? 1 : (G4int)((dphi-deg)/astep) + 1; 617 G4double ang = dphi/ksteps; 618 619 G4double sinHalf = std::sin(0.5*ang); 620 G4double cosHalf = std::cos(0.5*ang); 621 G4double sinStep = 2.*sinHalf*cosHalf; 622 G4double cosStep = 1. - 2.*sinHalf*sinHalf; 623 624 G4double sinStart = GetSinStartPhi(); 625 G4double cosStart = GetCosStartPhi(); 626 G4double sinEnd = GetSinEndPhi(); 627 G4double cosEnd = GetCosEndPhi(); 628 629 // define vectors and arrays 630 std::vector<const G4ThreeVectorList *> polygons; 631 polygons.resize(ksteps+2); 632 G4ThreeVectorList pols[NSTEPS+2]; 633 for (G4int k=0; k<ksteps+2; ++k) pols[k].resize(6); 634 for (G4int k=0; k<ksteps+2; ++k) polygons[k] = &pols[k]; 635 G4double r0[6],z0[6]; // contour with original edges of triangle 636 G4double r1[6]; // shifted radii of external edges of triangle 637 638 // main loop along triangles 639 pMin = kInfinity; 640 pMax =-kInfinity; 641 G4int ntria = (G4int)triangles.size()/3; 642 for (G4int i=0; i<ntria; ++i) 643 { 644 G4int i3 = i*3; 645 for (G4int k=0; k<3; ++k) 646 { 647 G4int e0 = i3+k, e1 = (k<2) ? e0+1 : i3; 648 G4int k2 = k*2; 649 // set contour with original edges of triangle 650 r0[k2+0] = triangles[e0].x(); z0[k2+0] = triangles[e0].y(); 651 r0[k2+1] = triangles[e1].x(); z0[k2+1] = triangles[e1].y(); 652 // set shifted radii 653 r1[k2+0] = r0[k2+0]; 654 r1[k2+1] = r0[k2+1]; 655 if (z0[k2+1] - z0[k2+0] <= 0) continue; 656 r1[k2+0] /= cosHalf; 657 r1[k2+1] /= cosHalf; 658 } 659 660 // rotate countour, set sequence of 6-sided polygons 661 G4double sinCur = sinStart*cosHalf + cosStart*sinHalf; 662 G4double cosCur = cosStart*cosHalf - sinStart*sinHalf; 663 for (G4int j=0; j<6; ++j) pols[0][j].set(r0[j]*cosStart,r0[j]*sinStart,z0[j]); 664 for (G4int k=1; k<ksteps+1; ++k) 665 { 666 for (G4int j=0; j<6; ++j) pols[k][j].set(r1[j]*cosCur,r1[j]*sinCur,z0[j]); 667 G4double sinTmp = sinCur; 668 sinCur = sinCur*cosStep + cosCur*sinStep; 669 cosCur = cosCur*cosStep - sinTmp*sinStep; 670 } 671 for (G4int j=0; j<6; ++j) pols[ksteps+1][j].set(r0[j]*cosEnd,r0[j]*sinEnd,z0[j]); 672 673 // set sub-envelope and adjust extent 674 G4double emin,emax; 675 G4BoundingEnvelope benv(polygons); 676 if (!benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,emin,emax)) continue; 677 if (emin < pMin) pMin = emin; 678 if (emax > pMax) pMax = emax; 679 if (eminlim > pMin && emaxlim < pMax) return true; // max possible extent 680 } 681 return (pMin < pMax); 682 } 683 684 // ComputeDimensions 685 // 686 void G4Polycone::ComputeDimensions( G4VPVParameterisation* p, 687 const G4int n, 688 const G4VPhysicalVolume* pRep ) 689 { 690 p->ComputeDimensions(*this,n,pRep); 691 } 692 693 // GetEntityType 694 // 695 G4GeometryType G4Polycone::GetEntityType() const 696 { 697 return {"G4Polycone"}; 698 } 699 700 // Make a clone of the object 701 // 702 G4VSolid* G4Polycone::Clone() const 703 { 704 return new G4Polycone(*this); 705 } 706 707 // 708 // Stream object contents to an output stream 709 // 710 std::ostream& G4Polycone::StreamInfo( std::ostream& os ) const 711 { 712 G4long oldprc = os.precision(16); 713 os << "-----------------------------------------------------------\n" 714 << " *** Dump for solid - " << GetName() << " ***\n" 715 << " ===================================================\n" 716 << " Solid type: G4Polycone\n" 717 << " Parameters: \n" 718 << " starting phi angle : " << startPhi/degree << " degrees \n" 719 << " ending phi angle : " << endPhi/degree << " degrees \n"; 720 G4int i=0; 721 722 G4int numPlanes = original_parameters->Num_z_planes; 723 os << " number of Z planes: " << numPlanes << "\n" 724 << " Z values: \n"; 725 for (i=0; i<numPlanes; ++i) 726 { 727 os << " Z plane " << i << ": " 728 << original_parameters->Z_values[i] << "\n"; 729 } 730 os << " Tangent distances to inner surface (Rmin): \n"; 731 for (i=0; i<numPlanes; ++i) 732 { 733 os << " Z plane " << i << ": " 734 << original_parameters->Rmin[i] << "\n"; 735 } 736 os << " Tangent distances to outer surface (Rmax): \n"; 737 for (i=0; i<numPlanes; ++i) 738 { 739 os << " Z plane " << i << ": " 740 << original_parameters->Rmax[i] << "\n"; 741 } 742 743 os << " number of RZ points: " << numCorner << "\n" 744 << " RZ values (corners): \n"; 745 for (i=0; i<numCorner; ++i) 746 { 747 os << " " 748 << corners[i].r << ", " << corners[i].z << "\n"; 749 } 750 os << "-----------------------------------------------------------\n"; 751 os.precision(oldprc); 752 753 return os; 754 } 755 756 ////////////////////////////////////////////////////////////////////////// 757 // 758 // Return volume 759 760 G4double G4Polycone::GetCubicVolume() 761 { 762 if (fCubicVolume == 0.) 763 { 764 G4double total = 0.; 765 G4int nrz = GetNumRZCorner(); 766 G4PolyconeSideRZ a = GetCorner(nrz - 1); 767 for (G4int i=0; i<nrz; ++i) 768 { 769 G4PolyconeSideRZ b = GetCorner(i); 770 total += (b.r*b.r + b.r*a.r + a.r*a.r)*(b.z - a.z); 771 a = b; 772 } 773 fCubicVolume = std::abs(total)*(GetEndPhi() - GetStartPhi())/6.; 774 } 775 return fCubicVolume; 776 } 777 778 ////////////////////////////////////////////////////////////////////////// 779 // 780 // Return surface area 781 782 G4double G4Polycone::GetSurfaceArea() 783 { 784 if (fSurfaceArea == 0.) 785 { 786 // phi cut area 787 G4int nrz = GetNumRZCorner(); 788 G4double scut = 0.; 789 if (IsOpen()) 790 { 791 G4PolyconeSideRZ a = GetCorner(nrz - 1); 792 for (G4int i=0; i<nrz; ++i) 793 { 794 G4PolyconeSideRZ b = GetCorner(i); 795 scut += a.r*b.z - a.z*b.r; 796 a = b; 797 } 798 scut = std::abs(scut); 799 } 800 // lateral surface area 801 G4double slat = 0; 802 G4PolyconeSideRZ a = GetCorner(nrz - 1); 803 for (G4int i=0; i<nrz; ++i) 804 { 805 G4PolyconeSideRZ b = GetCorner(i); 806 G4double h = std::sqrt((b.r - a.r)*(b.r - a.r) + (b.z - a.z)*(b.z - a.z)); 807 slat += (b.r + a.r)*h; 808 a = b; 809 } 810 slat *= (GetEndPhi() - GetStartPhi())/2.; 811 fSurfaceArea = scut + slat; 812 } 813 return fSurfaceArea; 814 } 815 816 ////////////////////////////////////////////////////////////////////////// 817 // 818 // Set vector of surface elements, auxiliary method for sampling 819 // random points on surface 820 821 void G4Polycone::SetSurfaceElements() const 822 { 823 fElements = new std::vector<G4Polycone::surface_element>; 824 G4double total = 0.; 825 G4int nrz = GetNumRZCorner(); 826 827 // set lateral surface elements 828 G4double dphi = GetEndPhi() - GetStartPhi(); 829 G4int ia = nrz - 1; 830 for (G4int ib=0; ib<nrz; ++ib) 831 { 832 G4PolyconeSideRZ a = GetCorner(ia); 833 G4PolyconeSideRZ b = GetCorner(ib); 834 G4Polycone::surface_element selem; 835 selem.i0 = ia; 836 selem.i1 = ib; 837 selem.i2 = -1; 838 ia = ib; 839 if (a.r == 0. && b.r == 0.) continue; 840 G4double h = std::sqrt((b.r - a.r)*(b.r - a.r) + (b.z - a.z)*(b.z - a.z)); 841 total += 0.5*dphi*(b.r + a.r)*h; 842 selem.area = total; 843 fElements->push_back(selem); 844 } 845 846 // set elements for phi cuts 847 if (IsOpen()) 848 { 849 G4TwoVectorList contourRZ; 850 std::vector<G4int> triangles; 851 for (G4int i=0; i<nrz; ++i) 852 { 853 G4PolyconeSideRZ corner = GetCorner(i); 854 contourRZ.emplace_back(corner.r, corner.z); 855 } 856 G4GeomTools::TriangulatePolygon(contourRZ, triangles); 857 auto ntria = (G4int)triangles.size(); 858 for (G4int i=0; i<ntria; i+=3) 859 { 860 G4Polycone::surface_element selem; 861 selem.i0 = triangles[i]; 862 selem.i1 = triangles[i+1]; 863 selem.i2 = triangles[i+2]; 864 G4PolyconeSideRZ a = GetCorner(selem.i0); 865 G4PolyconeSideRZ b = GetCorner(selem.i1); 866 G4PolyconeSideRZ c = GetCorner(selem.i2); 867 G4double stria = 868 std::abs(G4GeomTools::TriangleArea(a.r, a.z, b.r, b.z, c.r, c.z)); 869 total += stria; 870 selem.area = total; 871 fElements->push_back(selem); // start phi 872 total += stria; 873 selem.area = total; 874 selem.i0 += nrz; 875 fElements->push_back(selem); // end phi 876 } 877 } 878 } 879 880 ////////////////////////////////////////////////////////////////////////// 881 // 882 // Generate random point on surface 883 884 G4ThreeVector G4Polycone::GetPointOnSurface() const 885 { 886 // Set surface elements 887 if (fElements == nullptr) 888 { 889 G4AutoLock l(&surface_elementsMutex); 890 SetSurfaceElements(); 891 l.unlock(); 892 } 893 894 // Select surface element 895 G4Polycone::surface_element selem; 896 selem = fElements->back(); 897 G4double select = selem.area*G4QuickRand(); 898 auto it = std::lower_bound(fElements->begin(), fElements->end(), select, 899 [](const G4Polycone::surface_element& x, G4double val) 900 -> G4bool { return x.area < val; }); 901 902 // Generate random point 903 G4double r = 0, z = 0, phi = 0; 904 G4double u = G4QuickRand(); 905 G4double v = G4QuickRand(); 906 G4int i0 = (*it).i0; 907 G4int i1 = (*it).i1; 908 G4int i2 = (*it).i2; 909 if (i2 < 0) // lateral surface 910 { 911 G4PolyconeSideRZ p0 = GetCorner(i0); 912 G4PolyconeSideRZ p1 = GetCorner(i1); 913 if (p1.r < p0.r) 914 { 915 p0 = GetCorner(i1); 916 p1 = GetCorner(i0); 917 } 918 if (p1.r - p0.r < kCarTolerance) // cylindrical surface 919 { 920 r = (p1.r - p0.r)*u + p0.r; 921 z = (p1.z - p0.z)*u + p0.z; 922 } 923 else // conical surface 924 { 925 r = std::sqrt(p1.r*p1.r*u + p0.r*p0.r*(1. - u)); 926 z = p0.z + (p1.z - p0.z)*(r - p0.r)/(p1.r - p0.r); 927 } 928 phi = (GetEndPhi() - GetStartPhi())*v + GetStartPhi(); 929 } 930 else // phi cut 931 { 932 G4int nrz = GetNumRZCorner(); 933 phi = (i0 < nrz) ? GetStartPhi() : GetEndPhi(); 934 if (i0 >= nrz) { i0 -= nrz; } 935 G4PolyconeSideRZ p0 = GetCorner(i0); 936 G4PolyconeSideRZ p1 = GetCorner(i1); 937 G4PolyconeSideRZ p2 = GetCorner(i2); 938 if (u + v > 1.) { u = 1. - u; v = 1. - v; } 939 r = (p1.r - p0.r)*u + (p2.r - p0.r)*v + p0.r; 940 z = (p1.z - p0.z)*u + (p2.z - p0.z)*v + p0.z; 941 } 942 return { r*std::cos(phi), r*std::sin(phi), z }; 943 } 944 945 ////////////////////////////////////////////////////////////////////////// 946 // 947 // CreatePolyhedron 948 949 G4Polyhedron* G4Polycone::CreatePolyhedron() const 950 { 951 std::vector<G4TwoVector> rz(numCorner); 952 for (G4int i = 0; i < numCorner; ++i) 953 rz[i].set(corners[i].r, corners[i].z); 954 return new G4PolyhedronPcon(startPhi, endPhi - startPhi, rz); 955 } 956 957 // SetOriginalParameters 958 // 959 G4bool G4Polycone::SetOriginalParameters(G4ReduciblePolygon* rz) 960 { 961 G4int numPlanes = numCorner; 962 G4bool isConvertible = true; 963 G4double Zmax=rz->Bmax(); 964 rz->StartWithZMin(); 965 966 // Prepare vectors for storage 967 // 968 std::vector<G4double> Z; 969 std::vector<G4double> Rmin; 970 std::vector<G4double> Rmax; 971 972 G4int countPlanes=1; 973 G4int icurr=0; 974 G4int icurl=0; 975 976 // first plane Z=Z[0] 977 // 978 Z.push_back(corners[0].z); 979 G4double Zprev=Z[0]; 980 if (Zprev == corners[1].z) 981 { 982 Rmin.push_back(corners[0].r); 983 Rmax.push_back (corners[1].r);icurr=1; 984 } 985 else if (Zprev == corners[numPlanes-1].z) 986 { 987 Rmin.push_back(corners[numPlanes-1].r); 988 Rmax.push_back (corners[0].r); 989 icurl=numPlanes-1; 990 } 991 else 992 { 993 Rmin.push_back(corners[0].r); 994 Rmax.push_back (corners[0].r); 995 } 996 997 // next planes until last 998 // 999 G4int inextr=0, inextl=0; 1000 for (G4int i=0; i < numPlanes-2; ++i) 1001 { 1002 inextr=1+icurr; 1003 inextl=(icurl <= 0)? numPlanes-1 : icurl-1; 1004 1005 if((corners[inextr].z >= Zmax) & (corners[inextl].z >= Zmax)) { break; } 1006 1007 G4double Zleft = corners[inextl].z; 1008 G4double Zright = corners[inextr].z; 1009 if(Zright > Zleft) // Next plane will be Zleft 1010 { 1011 Z.push_back(Zleft); 1012 countPlanes++; 1013 G4double difZr=corners[inextr].z - corners[icurr].z; 1014 G4double difZl=corners[inextl].z - corners[icurl].z; 1015 1016 if(std::fabs(difZl) < kCarTolerance) 1017 { 1018 if(std::fabs(difZr) < kCarTolerance) 1019 { 1020 Rmin.push_back(corners[inextl].r); 1021 Rmax.push_back(corners[icurr].r); 1022 } 1023 else 1024 { 1025 Rmin.push_back(corners[inextl].r); 1026 Rmax.push_back(corners[icurr].r + (Zleft-corners[icurr].z)/difZr 1027 *(corners[inextr].r - corners[icurr].r)); 1028 } 1029 } 1030 else if (difZl >= kCarTolerance) 1031 { 1032 if(std::fabs(difZr) < kCarTolerance) 1033 { 1034 Rmin.push_back(corners[icurl].r); 1035 Rmax.push_back(corners[icurr].r); 1036 } 1037 else 1038 { 1039 Rmin.push_back(corners[icurl].r); 1040 Rmax.push_back(corners[icurr].r + (Zleft-corners[icurr].z)/difZr 1041 *(corners[inextr].r - corners[icurr].r)); 1042 } 1043 } 1044 else 1045 { 1046 isConvertible=false; break; 1047 } 1048 icurl=(icurl == 0)? numPlanes-1 : icurl-1; 1049 } 1050 else if(std::fabs(Zright-Zleft)<kCarTolerance) // Zright=Zleft 1051 { 1052 Z.push_back(Zleft); 1053 ++countPlanes; 1054 ++icurr; 1055 1056 icurl=(icurl == 0)? numPlanes-1 : icurl-1; 1057 1058 Rmin.push_back(corners[inextl].r); 1059 Rmax.push_back(corners[inextr].r); 1060 } 1061 else // Zright<Zleft 1062 { 1063 Z.push_back(Zright); 1064 ++countPlanes; 1065 1066 G4double difZr=corners[inextr].z - corners[icurr].z; 1067 G4double difZl=corners[inextl].z - corners[icurl].z; 1068 if(std::fabs(difZr) < kCarTolerance) 1069 { 1070 if(std::fabs(difZl) < kCarTolerance) 1071 { 1072 Rmax.push_back(corners[inextr].r); 1073 Rmin.push_back(corners[icurr].r); 1074 } 1075 else 1076 { 1077 Rmin.push_back(corners[icurl].r + (Zright-corners[icurl].z)/difZl 1078 *(corners[inextl].r - corners[icurl].r)); 1079 Rmax.push_back(corners[inextr].r); 1080 } 1081 ++icurr; 1082 } // plate 1083 else if (difZr >= kCarTolerance) 1084 { 1085 if(std::fabs(difZl) < kCarTolerance) 1086 { 1087 Rmax.push_back(corners[inextr].r); 1088 Rmin.push_back (corners[icurr].r); 1089 } 1090 else 1091 { 1092 Rmax.push_back(corners[inextr].r); 1093 Rmin.push_back (corners[icurl].r+(Zright-corners[icurl].z)/difZl 1094 * (corners[inextl].r - corners[icurl].r)); 1095 } 1096 ++icurr; 1097 } 1098 else 1099 { 1100 isConvertible=false; break; 1101 } 1102 } 1103 } // end for loop 1104 1105 // last plane Z=Zmax 1106 // 1107 Z.push_back(Zmax); 1108 ++countPlanes; 1109 inextr=1+icurr; 1110 inextl=(icurl <= 0)? numPlanes-1 : icurl-1; 1111 1112 Rmax.push_back(corners[inextr].r); 1113 Rmin.push_back(corners[inextl].r); 1114 1115 // Set original parameters Rmin,Rmax,Z 1116 // 1117 if(isConvertible) 1118 { 1119 original_parameters = new G4PolyconeHistorical; 1120 original_parameters->Z_values = new G4double[countPlanes]; 1121 original_parameters->Rmin = new G4double[countPlanes]; 1122 original_parameters->Rmax = new G4double[countPlanes]; 1123 1124 for(G4int j=0; j < countPlanes; ++j) 1125 { 1126 original_parameters->Z_values[j] = Z[j]; 1127 original_parameters->Rmax[j] = Rmax[j]; 1128 original_parameters->Rmin[j] = Rmin[j]; 1129 } 1130 original_parameters->Start_angle = startPhi; 1131 original_parameters->Opening_angle = endPhi-startPhi; 1132 original_parameters->Num_z_planes = countPlanes; 1133 1134 } 1135 else // Set parameters(r,z) with Rmin==0 as convention 1136 { 1137 #ifdef G4SPECSDEBUG 1138 std::ostringstream message; 1139 message << "Polycone " << GetName() << G4endl 1140 << "cannot be converted to Polycone with (Rmin,Rmaz,Z) parameters!"; 1141 G4Exception("G4Polycone::SetOriginalParameters()", "GeomSolids0002", 1142 JustWarning, message); 1143 #endif 1144 original_parameters = new G4PolyconeHistorical; 1145 original_parameters->Z_values = new G4double[numPlanes]; 1146 original_parameters->Rmin = new G4double[numPlanes]; 1147 original_parameters->Rmax = new G4double[numPlanes]; 1148 1149 for(G4int j=0; j < numPlanes; ++j) 1150 { 1151 original_parameters->Z_values[j] = corners[j].z; 1152 original_parameters->Rmax[j] = corners[j].r; 1153 original_parameters->Rmin[j] = 0.0; 1154 } 1155 original_parameters->Start_angle = startPhi; 1156 original_parameters->Opening_angle = endPhi-startPhi; 1157 original_parameters->Num_z_planes = numPlanes; 1158 } 1159 return isConvertible; 1160 } 1161 1162 #endif 1163