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 for G4UTorus wrapper class 27 // 28 // 19.08.15 Guilherme Lima, FNAL 29 // -------------------------------------------------------------------- 30 31 #include "G4Torus.hh" 32 #include "G4UTorus.hh" 33 34 #if ( defined(G4GEOM_USE_USOLIDS) || defined(G4GEOM_USE_PARTIAL_USOLIDS) ) 35 36 #include "G4TwoVector.hh" 37 #include "G4GeomTools.hh" 38 #include "G4AffineTransform.hh" 39 #include "G4BoundingEnvelope.hh" 40 41 #include "G4VPVParameterisation.hh" 42 43 using namespace CLHEP; 44 45 //////////////////////////////////////////////////////////////////////// 46 // 47 // Constructor - check & set half widths 48 49 50 G4UTorus::G4UTorus(const G4String& pName, 51 G4double rmin, G4double rmax, G4double rtor, 52 G4double sphi, G4double dphi) 53 : Base_t(pName, rmin, rmax, rtor, sphi, dphi) 54 { } 55 56 ////////////////////////////////////////////////////////////////////////// 57 // 58 // Fake default constructor - sets only member data and allocates memory 59 // for usage restricted to object persistency. 60 61 G4UTorus::G4UTorus( __void__& a ) 62 : Base_t(a) 63 { } 64 65 ////////////////////////////////////////////////////////////////////////// 66 // 67 // Destructor 68 69 G4UTorus::~G4UTorus() = default; 70 71 ////////////////////////////////////////////////////////////////////////// 72 // 73 // Copy constructor 74 75 G4UTorus::G4UTorus(const G4UTorus& rhs) 76 : Base_t(rhs) 77 { } 78 79 ////////////////////////////////////////////////////////////////////////// 80 // 81 // Assignment operator 82 83 G4UTorus& G4UTorus::operator = (const G4UTorus& rhs) 84 { 85 // Check assignment to self 86 // 87 if (this == &rhs) { return *this; } 88 89 // Copy base class data 90 // 91 Base_t::operator=(rhs); 92 93 return *this; 94 } 95 96 ////////////////////////////////////////////////////////////////////////// 97 // 98 // Accessors & modifiers 99 100 G4double G4UTorus::GetRmin() const 101 { 102 return rmin(); 103 } 104 105 G4double G4UTorus::GetRmax() const 106 { 107 return rmax(); 108 } 109 110 G4double G4UTorus::GetRtor() const 111 { 112 return rtor(); 113 } 114 115 G4double G4UTorus::GetSPhi() const 116 { 117 return sphi(); 118 } 119 120 G4double G4UTorus::GetDPhi() const 121 { 122 return dphi(); 123 } 124 125 G4double G4UTorus::GetSinStartPhi() const 126 { 127 return std::sin(sphi()); 128 } 129 130 G4double G4UTorus::GetCosStartPhi() const 131 { 132 return std::cos(sphi()); 133 } 134 135 G4double G4UTorus::GetSinEndPhi() const 136 { 137 return std::sin(sphi()+dphi()); 138 } 139 140 G4double G4UTorus::GetCosEndPhi() const 141 { 142 return std::cos(sphi()+dphi()); 143 } 144 145 void G4UTorus::SetRmin(G4double arg) 146 { 147 Base_t::SetRMin(arg); 148 fRebuildPolyhedron = true; 149 } 150 151 void G4UTorus::SetRmax(G4double arg) 152 { 153 Base_t::SetRMax(arg); 154 fRebuildPolyhedron = true; 155 } 156 157 void G4UTorus::SetRtor(G4double arg) 158 { 159 Base_t::SetRTor(arg); 160 fRebuildPolyhedron = true; 161 } 162 163 void G4UTorus::SetSPhi(G4double arg) 164 { 165 Base_t::SetSPhi(arg); 166 fRebuildPolyhedron = true; 167 } 168 169 void G4UTorus::SetDPhi(G4double arg) 170 { 171 Base_t::SetDPhi(arg); 172 fRebuildPolyhedron = true; 173 } 174 175 void G4UTorus::SetAllParameters(G4double arg1, G4double arg2, 176 G4double arg3, G4double arg4, G4double arg5) 177 { 178 SetRmin(arg1); 179 SetRmax(arg2); 180 SetRtor(arg3); 181 SetSPhi(arg4); 182 SetDPhi(arg5); 183 fRebuildPolyhedron = true; 184 } 185 186 //////////////////////////////////////////////////////////////////////// 187 // 188 // Dispatch to parameterisation for replication mechanism dimension 189 // computation & modification. 190 191 void G4UTorus::ComputeDimensions(G4VPVParameterisation* p, 192 const G4int n, 193 const G4VPhysicalVolume* pRep) 194 { 195 p->ComputeDimensions(*(G4Torus*)this,n,pRep); 196 } 197 198 ////////////////////////////////////////////////////////////////////////// 199 // 200 // Make a clone of the object 201 202 G4VSolid* G4UTorus::Clone() const 203 { 204 return new G4UTorus(*this); 205 } 206 207 ////////////////////////////////////////////////////////////////////////// 208 // 209 // Get bounding box 210 211 void G4UTorus::BoundingLimits(G4ThreeVector& pMin, G4ThreeVector& pMax) const 212 { 213 static G4bool checkBBox = true; 214 215 G4double rmax = GetRmax(); 216 G4double rtor = GetRtor(); 217 G4double rint = rtor - rmax; 218 G4double rext = rtor + rmax; 219 G4double dz = rmax; 220 221 // Find bounding box 222 // 223 if (GetDPhi() >= twopi) 224 { 225 pMin.set(-rext,-rext,-dz); 226 pMax.set( rext, rext, dz); 227 } 228 else 229 { 230 G4TwoVector vmin,vmax; 231 G4GeomTools::DiskExtent(rint,rext, 232 GetSinStartPhi(),GetCosStartPhi(), 233 GetSinEndPhi(),GetCosEndPhi(), 234 vmin,vmax); 235 pMin.set(vmin.x(),vmin.y(),-dz); 236 pMax.set(vmax.x(),vmax.y(), dz); 237 } 238 239 // Check correctness of the bounding box 240 // 241 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z()) 242 { 243 std::ostringstream message; 244 message << "Bad bounding box (min >= max) for solid: " 245 << GetName() << " !" 246 << "\npMin = " << pMin 247 << "\npMax = " << pMax; 248 G4Exception("G4UTorus::BoundingLimits()", "GeomMgt0001", 249 JustWarning, message); 250 StreamInfo(G4cout); 251 } 252 253 // Check consistency of bounding boxes 254 // 255 if (checkBBox) 256 { 257 U3Vector vmin, vmax; 258 Base_t::Extent(vmin,vmax); 259 if (std::abs(pMin.x()-vmin.x()) > kCarTolerance || 260 std::abs(pMin.y()-vmin.y()) > kCarTolerance || 261 std::abs(pMin.z()-vmin.z()) > kCarTolerance || 262 std::abs(pMax.x()-vmax.x()) > kCarTolerance || 263 std::abs(pMax.y()-vmax.y()) > kCarTolerance || 264 std::abs(pMax.z()-vmax.z()) > kCarTolerance) 265 { 266 std::ostringstream message; 267 message << "Inconsistency in bounding boxes for solid: " 268 << GetName() << " !" 269 << "\nBBox min: wrapper = " << pMin << " solid = " << vmin 270 << "\nBBox max: wrapper = " << pMax << " solid = " << vmax; 271 G4Exception("G4UTorus::BoundingLimits()", "GeomMgt0001", 272 JustWarning, message); 273 checkBBox = false; 274 } 275 } 276 } 277 278 ////////////////////////////////////////////////////////////////////////// 279 // 280 // Calculate extent under transform and specified limit 281 282 G4bool 283 G4UTorus::CalculateExtent(const EAxis pAxis, 284 const G4VoxelLimits& pVoxelLimit, 285 const G4AffineTransform& pTransform, 286 G4double& pMin, G4double& pMax) const 287 { 288 G4ThreeVector bmin, bmax; 289 G4bool exist; 290 291 // Get bounding box 292 BoundingLimits(bmin,bmax); 293 294 // Check bounding box 295 G4BoundingEnvelope bbox(bmin,bmax); 296 #ifdef G4BBOX_EXTENT 297 if (true) return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax); 298 #endif 299 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax)) 300 { 301 return exist = pMin < pMax; 302 } 303 304 // Get parameters of the solid 305 G4double rmin = GetRmin(); 306 G4double rmax = GetRmax(); 307 G4double rtor = GetRtor(); 308 G4double dphi = GetDPhi(); 309 G4double sinStart = GetSinStartPhi(); 310 G4double cosStart = GetCosStartPhi(); 311 G4double sinEnd = GetSinEndPhi(); 312 G4double cosEnd = GetCosEndPhi(); 313 G4double rint = rtor - rmax; 314 G4double rext = rtor + rmax; 315 316 // Find bounding envelope and calculate extent 317 // 318 static const G4int NPHI = 24; // number of steps for whole torus 319 static const G4int NDISK = 16; // number of steps for disk 320 static const G4double sinHalfDisk = std::sin(pi/NDISK); 321 static const G4double cosHalfDisk = std::cos(pi/NDISK); 322 static const G4double sinStepDisk = 2.*sinHalfDisk*cosHalfDisk; 323 static const G4double cosStepDisk = 1. - 2.*sinHalfDisk*sinHalfDisk; 324 325 G4double astep = (360/NPHI)*deg; // max angle for one slice in phi 326 G4int kphi = (dphi <= astep) ? 1 : (G4int)((dphi-deg)/astep) + 1; 327 G4double ang = dphi/kphi; 328 329 G4double sinHalf = std::sin(0.5*ang); 330 G4double cosHalf = std::cos(0.5*ang); 331 G4double sinStep = 2.*sinHalf*cosHalf; 332 G4double cosStep = 1. - 2.*sinHalf*sinHalf; 333 334 // define vectors for bounding envelope 335 G4ThreeVectorList pols[NDISK+1]; 336 for (auto & pol : pols) pol.resize(4); 337 338 std::vector<const G4ThreeVectorList *> polygons; 339 polygons.resize(NDISK+1); 340 for (G4int k=0; k<NDISK+1; ++k) polygons[k] = &pols[k]; 341 342 // set internal and external reference circles 343 G4TwoVector rzmin[NDISK]; 344 G4TwoVector rzmax[NDISK]; 345 346 if ((rtor-rmin*sinHalfDisk)/cosHalf > (rtor+rmin*sinHalfDisk)) rmin = 0; 347 rmax /= cosHalfDisk; 348 G4double sinCurDisk = sinHalfDisk; 349 G4double cosCurDisk = cosHalfDisk; 350 for (G4int k=0; k<NDISK; ++k) 351 { 352 G4double rmincur = rtor + rmin*cosCurDisk; 353 if (cosCurDisk < 0 && rmin > 0) rmincur /= cosHalf; 354 rzmin[k].set(rmincur,rmin*sinCurDisk); 355 356 G4double rmaxcur = rtor + rmax*cosCurDisk; 357 if (cosCurDisk > 0) rmaxcur /= cosHalf; 358 rzmax[k].set(rmaxcur,rmax*sinCurDisk); 359 360 G4double sinTmpDisk = sinCurDisk; 361 sinCurDisk = sinCurDisk*cosStepDisk + cosCurDisk*sinStepDisk; 362 cosCurDisk = cosCurDisk*cosStepDisk - sinTmpDisk*sinStepDisk; 363 } 364 365 // Loop along slices in Phi. The extent is calculated as cumulative 366 // extent of the slices 367 pMin = kInfinity; 368 pMax = -kInfinity; 369 G4double eminlim = pVoxelLimit.GetMinExtent(pAxis); 370 G4double emaxlim = pVoxelLimit.GetMaxExtent(pAxis); 371 G4double sinCur1 = 0, cosCur1 = 0, sinCur2 = 0, cosCur2 = 0; 372 for (G4int i=0; i<kphi+1; ++i) 373 { 374 if (i == 0) 375 { 376 sinCur1 = sinStart; 377 cosCur1 = cosStart; 378 sinCur2 = sinCur1*cosHalf + cosCur1*sinHalf; 379 cosCur2 = cosCur1*cosHalf - sinCur1*sinHalf; 380 } 381 else 382 { 383 sinCur1 = sinCur2; 384 cosCur1 = cosCur2; 385 sinCur2 = (i == kphi) ? sinEnd : sinCur1*cosStep + cosCur1*sinStep; 386 cosCur2 = (i == kphi) ? cosEnd : cosCur1*cosStep - sinCur1*sinStep; 387 } 388 for (G4int k=0; k<NDISK; ++k) 389 { 390 G4double r1 = rzmin[k].x(), r2 = rzmax[k].x(); 391 G4double z1 = rzmin[k].y(), z2 = rzmax[k].y(); 392 pols[k][0].set(r1*cosCur1,r1*sinCur1,z1); 393 pols[k][1].set(r2*cosCur1,r2*sinCur1,z2); 394 pols[k][2].set(r2*cosCur2,r2*sinCur2,z2); 395 pols[k][3].set(r1*cosCur2,r1*sinCur2,z1); 396 } 397 pols[NDISK] = pols[0]; 398 399 // get bounding box of current slice 400 G4TwoVector vmin,vmax; 401 G4GeomTools:: 402 DiskExtent(rint,rext,sinCur1,cosCur1,sinCur2,cosCur2,vmin,vmax); 403 bmin.setX(vmin.x()); bmin.setY(vmin.y()); 404 bmax.setX(vmax.x()); bmax.setY(vmax.y()); 405 406 // set bounding envelope for current slice and adjust extent 407 G4double emin,emax; 408 G4BoundingEnvelope benv(bmin,bmax,polygons); 409 if (!benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,emin,emax)) continue; 410 if (emin < pMin) pMin = emin; 411 if (emax > pMax) pMax = emax; 412 if (eminlim > pMin && emaxlim < pMax) break; // max possible extent 413 } 414 return (pMin < pMax); 415 } 416 417 ////////////////////////////////////////////////////////////////////////// 418 // 419 // Create polyhedron for visualization 420 421 G4Polyhedron* G4UTorus::CreatePolyhedron() const 422 { 423 return new G4PolyhedronTorus(GetRmin(), 424 GetRmax(), 425 GetRtor(), 426 GetSPhi(), 427 GetDPhi()); 428 } 429 430 #endif // G4GEOM_USE_USOLIDS 431