Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // 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(G 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, G4doub 52 G4double sphi, G4doub 53 : Base_t(pName, rmin, rmax, rtor, sphi, dphi 54 { } 55 56 ////////////////////////////////////////////// 57 // 58 // Fake default constructor - sets only member 59 // for usage restri 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 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, 176 G4double arg3, 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 replicatio 189 // computation & modification. 190 191 void G4UTorus::ComputeDimensions(G4VPVParamete 192 const G4int n 193 const G4VPhys 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& p 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(),G 233 GetSinEndPhi(),Get 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 242 { 243 std::ostringstream message; 244 message << "Bad bounding box (min >= max) 245 << GetName() << " !" 246 << "\npMin = " << pMin 247 << "\npMax = " << pMax; 248 G4Exception("G4UTorus::BoundingLimits()", 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()) > kCarTole 260 std::abs(pMin.y()-vmin.y()) > kCarTole 261 std::abs(pMin.z()-vmin.z()) > kCarTole 262 std::abs(pMax.x()-vmax.x()) > kCarTole 263 std::abs(pMax.y()-vmax.y()) > kCarTole 264 std::abs(pMax.z()-vmax.z()) > kCarTole 265 { 266 std::ostringstream message; 267 message << "Inconsistency in bounding bo 268 << GetName() << " !" 269 << "\nBBox min: wrapper = " << p 270 << "\nBBox max: wrapper = " << p 271 G4Exception("G4UTorus::BoundingLimits()" 272 JustWarning, message); 273 checkBBox = false; 274 } 275 } 276 } 277 278 ////////////////////////////////////////////// 279 // 280 // Calculate extent under transform and specif 281 282 G4bool 283 G4UTorus::CalculateExtent(const EAxis pAxis, 284 const G4VoxelLimits& 285 const G4AffineTransf 286 G4double& pMin 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, 298 #endif 299 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox 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 exte 317 // 318 static const G4int NPHI = 24; // number of 319 static const G4int NDISK = 16; // number of 320 static const G4double sinHalfDisk = std::sin 321 static const G4double cosHalfDisk = std::cos 322 static const G4double sinStepDisk = 2.*sinHa 323 static const G4double cosStepDisk = 1. - 2.* 324 325 G4double astep = (360/NPHI)*deg; // max angl 326 G4int kphi = (dphi <= astep) ? 1 : (G4in 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 *> polyg 339 polygons.resize(NDISK+1); 340 for (G4int k=0; k<NDISK+1; ++k) polygons[k] 341 342 // set internal and external reference circl 343 G4TwoVector rzmin[NDISK]; 344 G4TwoVector rzmax[NDISK]; 345 346 if ((rtor-rmin*sinHalfDisk)/cosHalf > (rtor+ 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 /= 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 + cosC 362 cosCurDisk = cosCurDisk*cosStepDisk - sinT 363 } 364 365 // Loop along slices in Phi. The extent is c 366 // extent of the slices 367 pMin = kInfinity; 368 pMax = -kInfinity; 369 G4double eminlim = pVoxelLimit.GetMinExtent( 370 G4double emaxlim = pVoxelLimit.GetMaxExtent( 371 G4double sinCur1 = 0, cosCur1 = 0, sinCur2 = 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*sinH 379 cosCur2 = cosCur1*cosHalf - sinCur1*sinH 380 } 381 else 382 { 383 sinCur1 = sinCur2; 384 cosCur1 = cosCur2; 385 sinCur2 = (i == kphi) ? sinEnd : sinCur1 386 cosCur2 = (i == kphi) ? cosEnd : cosCur1 387 } 388 for (G4int k=0; k<NDISK; ++k) 389 { 390 G4double r1 = rzmin[k].x(), r2 = rzmax[k 391 G4double z1 = rzmin[k].y(), z2 = rzmax[k 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,sin 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 407 G4double emin,emax; 408 G4BoundingEnvelope benv(bmin,bmax,polygons 409 if (!benv.CalculateExtent(pAxis,pVoxelLimi 410 if (emin < pMin) pMin = emin; 411 if (emax > pMax) pMax = emax; 412 if (eminlim > pMin && emaxlim < pMax) brea 413 } 414 return (pMin < pMax); 415 } 416 417 ////////////////////////////////////////////// 418 // 419 // Create polyhedron for visualization 420 421 G4Polyhedron* G4UTorus::CreatePolyhedron() con 422 { 423 return new G4PolyhedronTorus(GetRmin(), 424 GetRmax(), 425 GetRtor(), 426 GetSPhi(), 427 GetDPhi()); 428 } 429 430 #endif // G4GEOM_USE_USOLIDS 431