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 G4UTubs wrapper class 27 // 28 // 30.10.13 G.Cosmo, CERN/PH 29 // ------------------------------------------- 30 31 #include "G4Tubs.hh" 32 #include "G4UTubs.hh" 33 34 #if ( defined(G4GEOM_USE_USOLIDS) || defined(G 35 36 #include "G4GeomTools.hh" 37 #include "G4AffineTransform.hh" 38 #include "G4VPVParameterisation.hh" 39 #include "G4BoundingEnvelope.hh" 40 41 using namespace CLHEP; 42 43 ////////////////////////////////////////////// 44 // 45 // Constructor - check parameters, convert ang 46 // - note if pdphi>2PI then reset 47 48 G4UTubs::G4UTubs( const G4String& pName, 49 G4double pRMin, G4doub 50 G4double pDz, 51 G4double pSPhi, G4doub 52 : Base_t(pName, pRMin, pRMax, pDz, pSPhi, pD 53 { 54 } 55 56 ////////////////////////////////////////////// 57 // 58 // Fake default constructor - sets only member 59 // for usage restri 60 // 61 G4UTubs::G4UTubs( __void__& a ) 62 : Base_t(a) 63 { 64 } 65 66 ////////////////////////////////////////////// 67 // 68 // Destructor 69 70 G4UTubs::~G4UTubs() = default; 71 72 ////////////////////////////////////////////// 73 // 74 // Copy constructor 75 76 G4UTubs::G4UTubs(const G4UTubs& rhs) 77 : Base_t(rhs) 78 { 79 } 80 81 ////////////////////////////////////////////// 82 // 83 // Assignment operator 84 85 G4UTubs& G4UTubs::operator = (const G4UTubs& r 86 { 87 // Check assignment to self 88 // 89 if (this == &rhs) { return *this; } 90 91 // Copy base class data 92 // 93 Base_t::operator=(rhs); 94 95 return *this; 96 } 97 98 ////////////////////////////////////////////// 99 // 100 // Accessors and modifiers 101 102 G4double G4UTubs::GetInnerRadius() const 103 { 104 return rmin(); 105 } 106 G4double G4UTubs::GetOuterRadius() const 107 { 108 return rmax(); 109 } 110 G4double G4UTubs::GetZHalfLength() const 111 { 112 return z(); 113 } 114 G4double G4UTubs::GetStartPhiAngle() const 115 { 116 return sphi(); 117 } 118 G4double G4UTubs::GetDeltaPhiAngle() const 119 { 120 return dphi(); 121 } 122 G4double G4UTubs::GetSinStartPhi() const 123 { 124 return std::sin(GetStartPhiAngle()); 125 } 126 G4double G4UTubs::GetCosStartPhi() const 127 { 128 return std::cos(GetStartPhiAngle()); 129 } 130 G4double G4UTubs::GetSinEndPhi() const 131 { 132 return std::sin(GetStartPhiAngle()+GetDeltaP 133 } 134 G4double G4UTubs::GetCosEndPhi() const 135 { 136 return std::cos(GetStartPhiAngle()+GetDeltaP 137 } 138 139 void G4UTubs::SetInnerRadius(G4double newRMin) 140 { 141 SetRMin(newRMin); 142 fRebuildPolyhedron = true; 143 } 144 void G4UTubs::SetOuterRadius(G4double newRMax) 145 { 146 SetRMax(newRMax); 147 fRebuildPolyhedron = true; 148 } 149 void G4UTubs::SetZHalfLength(G4double newDz) 150 { 151 SetDz(newDz); 152 fRebuildPolyhedron = true; 153 } 154 void G4UTubs::SetStartPhiAngle(G4double newSPh 155 { 156 SetSPhi(newSPhi); 157 fRebuildPolyhedron = true; 158 } 159 void G4UTubs::SetDeltaPhiAngle(G4double newDPh 160 { 161 SetDPhi(newDPhi); 162 fRebuildPolyhedron = true; 163 } 164 165 ////////////////////////////////////////////// 166 // 167 // Dispatch to parameterisation for replicatio 168 // computation & modification. 169 170 void G4UTubs::ComputeDimensions( G4VPVPar 171 const G4int n, 172 const G4VPhysi 173 { 174 p->ComputeDimensions(*(G4Tubs*)this,n,pRep) 175 } 176 177 ////////////////////////////////////////////// 178 // 179 // Make a clone of the object 180 181 G4VSolid* G4UTubs::Clone() const 182 { 183 return new G4UTubs(*this); 184 } 185 186 ////////////////////////////////////////////// 187 // 188 // Get bounding box 189 190 void G4UTubs::BoundingLimits(G4ThreeVector& pM 191 { 192 static G4bool checkBBox = true; 193 194 G4double rmin = GetInnerRadius(); 195 G4double rmax = GetOuterRadius(); 196 G4double dz = GetZHalfLength(); 197 198 // Find bounding box 199 // 200 if (GetDeltaPhiAngle() < twopi) 201 { 202 G4TwoVector vmin,vmax; 203 G4GeomTools::DiskExtent(rmin,rmax, 204 GetSinStartPhi(),G 205 GetSinEndPhi(),Get 206 vmin,vmax); 207 pMin.set(vmin.x(),vmin.y(),-dz); 208 pMax.set(vmax.x(),vmax.y(), dz); 209 } 210 else 211 { 212 pMin.set(-rmax,-rmax,-dz); 213 pMax.set( rmax, rmax, dz); 214 } 215 216 // Check correctness of the bounding box 217 // 218 if (pMin.x() >= pMax.x() || pMin.y() >= pMax 219 { 220 std::ostringstream message; 221 message << "Bad bounding box (min >= max) 222 << GetName() << " !" 223 << "\npMin = " << pMin 224 << "\npMax = " << pMax; 225 G4Exception("G4UTubs::BoundingLimits()", " 226 JustWarning, message); 227 StreamInfo(G4cout); 228 } 229 230 // Check consistency of bounding boxes 231 // 232 if (checkBBox) 233 { 234 U3Vector vmin, vmax; 235 Extent(vmin,vmax); 236 if (std::abs(pMin.x()-vmin.x()) > kCarTole 237 std::abs(pMin.y()-vmin.y()) > kCarTole 238 std::abs(pMin.z()-vmin.z()) > kCarTole 239 std::abs(pMax.x()-vmax.x()) > kCarTole 240 std::abs(pMax.y()-vmax.y()) > kCarTole 241 std::abs(pMax.z()-vmax.z()) > kCarTole 242 { 243 std::ostringstream message; 244 message << "Inconsistency in bounding bo 245 << GetName() << " !" 246 << "\nBBox min: wrapper = " << p 247 << "\nBBox max: wrapper = " << p 248 G4Exception("G4UTubs::BoundingLimits()", 249 JustWarning, message); 250 checkBBox = false; 251 } 252 } 253 } 254 255 ////////////////////////////////////////////// 256 // 257 // Calculate extent under transform and specif 258 259 G4bool 260 G4UTubs::CalculateExtent(const EAxis pAxis, 261 const G4VoxelLimits& 262 const G4AffineTransfo 263 G4double& pMin, 264 { 265 G4ThreeVector bmin, bmax; 266 G4bool exist; 267 268 // Get bounding box 269 BoundingLimits(bmin,bmax); 270 271 // Check bounding box 272 G4BoundingEnvelope bbox(bmin,bmax); 273 #ifdef G4BBOX_EXTENT 274 if (true) return bbox.CalculateExtent(pAxis, 275 #endif 276 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox 277 { 278 return exist = pMin < pMax; 279 } 280 281 // Get parameters of the solid 282 G4double rmin = GetInnerRadius(); 283 G4double rmax = GetOuterRadius(); 284 G4double dz = GetZHalfLength(); 285 G4double dphi = GetDeltaPhiAngle(); 286 287 // Find bounding envelope and calculate exte 288 // 289 const G4int NSTEPS = 24; // numbe 290 G4double astep = twopi/NSTEPS; // max a 291 G4int ksteps = (dphi <= astep) ? 1 : (G4i 292 G4double ang = dphi/ksteps; 293 294 G4double sinHalf = std::sin(0.5*ang); 295 G4double cosHalf = std::cos(0.5*ang); 296 G4double sinStep = 2.*sinHalf*cosHalf; 297 G4double cosStep = 1. - 2.*sinHalf*sinHalf; 298 G4double rext = rmax/cosHalf; 299 300 // bounding envelope for full cylinder consi 301 // in other cases it is a sequence of quadri 302 if (rmin == 0 && dphi == twopi) 303 { 304 G4double sinCur = sinHalf; 305 G4double cosCur = cosHalf; 306 307 G4ThreeVectorList baseA(NSTEPS),baseB(NSTE 308 for (G4int k=0; k<NSTEPS; ++k) 309 { 310 baseA[k].set(rext*cosCur,rext*sinCur,-dz 311 baseB[k].set(rext*cosCur,rext*sinCur, dz 312 313 G4double sinTmp = sinCur; 314 sinCur = sinCur*cosStep + cosCur*sinStep 315 cosCur = cosCur*cosStep - sinTmp*sinStep 316 } 317 std::vector<const G4ThreeVectorList *> pol 318 polygons[0] = &baseA; 319 polygons[1] = &baseB; 320 G4BoundingEnvelope benv(bmin,bmax,polygons 321 exist = benv.CalculateExtent(pAxis,pVoxelL 322 } 323 else 324 { 325 G4double sinStart = GetSinStartPhi(); 326 G4double cosStart = GetCosStartPhi(); 327 G4double sinEnd = GetSinEndPhi(); 328 G4double cosEnd = GetCosEndPhi(); 329 G4double sinCur = sinStart*cosHalf + cos 330 G4double cosCur = cosStart*cosHalf - sin 331 332 // set quadrilaterals 333 G4ThreeVectorList pols[NSTEPS+2]; 334 for (G4int k=0; k<ksteps+2; ++k) pols[k].r 335 pols[0][0].set(rmin*cosStart,rmin*sinStart 336 pols[0][1].set(rmin*cosStart,rmin*sinStart 337 pols[0][2].set(rmax*cosStart,rmax*sinStart 338 pols[0][3].set(rmax*cosStart,rmax*sinStart 339 for (G4int k=1; k<ksteps+1; ++k) 340 { 341 pols[k][0].set(rmin*cosCur,rmin*sinCur, 342 pols[k][1].set(rmin*cosCur,rmin*sinCur,- 343 pols[k][2].set(rext*cosCur,rext*sinCur,- 344 pols[k][3].set(rext*cosCur,rext*sinCur, 345 346 G4double sinTmp = sinCur; 347 sinCur = sinCur*cosStep + cosCur*sinStep 348 cosCur = cosCur*cosStep - sinTmp*sinStep 349 } 350 pols[ksteps+1][0].set(rmin*cosEnd,rmin*sin 351 pols[ksteps+1][1].set(rmin*cosEnd,rmin*sin 352 pols[ksteps+1][2].set(rmax*cosEnd,rmax*sin 353 pols[ksteps+1][3].set(rmax*cosEnd,rmax*sin 354 355 // set envelope and calculate extent 356 std::vector<const G4ThreeVectorList *> pol 357 polygons.resize(ksteps+2); 358 for (G4int k=0; k<ksteps+2; ++k) polygons[ 359 G4BoundingEnvelope benv(bmin,bmax,polygons 360 exist = benv.CalculateExtent(pAxis,pVoxelL 361 } 362 return exist; 363 } 364 365 ////////////////////////////////////////////// 366 // 367 // Create polyhedron for visualization 368 // 369 G4Polyhedron* G4UTubs::CreatePolyhedron() cons 370 { 371 return new G4PolyhedronTubs(GetInnerRadius() 372 GetOuterRadius() 373 GetZHalfLength() 374 GetStartPhiAngle 375 GetDeltaPhiAngle 376 } 377 378 #endif // G4GEOM_USE_USOLIDS 379