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 G4UCutTubs wrapper class 27 // 28 // 07.07.17 G.Cosmo, CERN/PH 29 // ------------------------------------------- 30 31 #include "G4CutTubs.hh" 32 #include "G4UCutTubs.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 G4UCutTubs::G4UCutTubs( const G4String& pName, 49 G4double pRMin, 50 G4double pDz, 51 G4double pSPhi, 52 const G4ThreeVec 53 const G4ThreeVec 54 : Base_t(pName, pRMin, pRMax, pDz, pSPhi, pD 55 pLowNorm.x(), pLowNorm.y(), pLowNor 56 pHighNorm.x(), pHighNorm.y(), pHigh 57 { 58 } 59 60 ////////////////////////////////////////////// 61 // 62 // Fake default constructor - sets only member 63 // for usage restri 64 // 65 G4UCutTubs::G4UCutTubs( __void__& a ) 66 : Base_t(a) 67 { 68 } 69 70 ////////////////////////////////////////////// 71 // 72 // Destructor 73 74 G4UCutTubs::~G4UCutTubs() = default; 75 76 ////////////////////////////////////////////// 77 // 78 // Copy constructor 79 80 G4UCutTubs::G4UCutTubs(const G4UCutTubs& rhs) 81 : Base_t(rhs) 82 { 83 } 84 85 ////////////////////////////////////////////// 86 // 87 // Assignment operator 88 89 G4UCutTubs& G4UCutTubs::operator = (const G4UC 90 { 91 // Check assignment to self 92 // 93 if (this == &rhs) { return *this; } 94 95 // Copy base class data 96 // 97 Base_t::operator=(rhs); 98 99 return *this; 100 } 101 102 ////////////////////////////////////////////// 103 // 104 // Accessors and modifiers 105 106 G4double G4UCutTubs::GetInnerRadius() const 107 { 108 return rmin(); 109 } 110 G4double G4UCutTubs::GetOuterRadius() const 111 { 112 return rmax(); 113 } 114 G4double G4UCutTubs::GetZHalfLength() const 115 { 116 return z(); 117 } 118 G4double G4UCutTubs::GetStartPhiAngle() const 119 { 120 return sphi(); 121 } 122 G4double G4UCutTubs::GetDeltaPhiAngle() const 123 { 124 return dphi(); 125 } 126 G4double G4UCutTubs::GetSinStartPhi() const 127 { 128 return std::sin(GetStartPhiAngle()); 129 } 130 G4double G4UCutTubs::GetCosStartPhi() const 131 { 132 return std::cos(GetStartPhiAngle()); 133 } 134 G4double G4UCutTubs::GetSinEndPhi() const 135 { 136 return std::sin(GetStartPhiAngle()+GetDeltaP 137 } 138 G4double G4UCutTubs::GetCosEndPhi() const 139 { 140 return std::cos(GetStartPhiAngle()+GetDeltaP 141 } 142 G4ThreeVector G4UCutTubs::GetLowNorm () const 143 { 144 U3Vector lc = BottomNormal(); 145 return {lc.x(), lc.y(), lc.z()}; 146 } 147 G4ThreeVector G4UCutTubs::GetHighNorm () const 148 { 149 U3Vector hc = TopNormal(); 150 return {hc.x(), hc.y(), hc.z()}; 151 } 152 153 void G4UCutTubs::SetInnerRadius(G4double newRM 154 { 155 SetRMin(newRMin); 156 fRebuildPolyhedron = true; 157 } 158 void G4UCutTubs::SetOuterRadius(G4double newRM 159 { 160 SetRMax(newRMax); 161 fRebuildPolyhedron = true; 162 } 163 void G4UCutTubs::SetZHalfLength(G4double newDz 164 { 165 SetDz(newDz); 166 fRebuildPolyhedron = true; 167 } 168 void G4UCutTubs::SetStartPhiAngle(G4double new 169 { 170 SetSPhi(newSPhi); 171 fRebuildPolyhedron = true; 172 } 173 void G4UCutTubs::SetDeltaPhiAngle(G4double new 174 { 175 SetDPhi(newDPhi); 176 fRebuildPolyhedron = true; 177 } 178 179 ////////////////////////////////////////////// 180 // 181 // Make a clone of the object 182 183 G4VSolid* G4UCutTubs::Clone() const 184 { 185 return new G4UCutTubs(*this); 186 } 187 188 ////////////////////////////////////////////// 189 // 190 // Get bounding box 191 192 void G4UCutTubs::BoundingLimits(G4ThreeVector& 193 { 194 static G4bool checkBBox = true; 195 196 G4double rmin = GetInnerRadius(); 197 G4double rmax = GetOuterRadius(); 198 G4double dz = GetZHalfLength(); 199 G4double dphi = GetDeltaPhiAngle(); 200 201 G4double sinSphi = GetSinStartPhi(); 202 G4double cosSphi = GetCosStartPhi(); 203 G4double sinEphi = GetSinEndPhi(); 204 G4double cosEphi = GetCosEndPhi(); 205 206 G4ThreeVector norm; 207 G4double mag, topx, topy, dists, diste; 208 G4bool iftop; 209 210 // Find Zmin 211 // 212 G4double zmin; 213 norm = GetLowNorm(); 214 mag = std::sqrt(norm.x()*norm.x() + norm.y( 215 topx = (mag == 0) ? 0 : -rmax*norm.x()/mag; 216 topy = (mag == 0) ? 0 : -rmax*norm.y()/mag; 217 dists = sinSphi*topx - cosSphi*topy; 218 diste = -sinEphi*topx + cosEphi*topy; 219 if (dphi > pi) 220 { 221 iftop = true; 222 if (dists > 0 && diste > 0)iftop = false; 223 } 224 else 225 { 226 iftop = false; 227 if (dists <= 0 && diste <= 0) iftop = true 228 } 229 if (iftop) 230 { 231 zmin = -(norm.x()*topx + norm.y()*topy)/no 232 } 233 else 234 { 235 G4double z1 = -rmin*(norm.x()*cosSphi + no 236 G4double z2 = -rmin*(norm.x()*cosEphi + no 237 G4double z3 = -rmax*(norm.x()*cosSphi + no 238 G4double z4 = -rmax*(norm.x()*cosEphi + no 239 zmin = std::min(std::min(std::min(z1,z2),z 240 } 241 242 // Find Zmax 243 // 244 G4double zmax; 245 norm = GetHighNorm(); 246 mag = std::sqrt(norm.x()*norm.x() + norm.y( 247 topx = (mag == 0) ? 0 : -rmax*norm.x()/mag; 248 topy = (mag == 0) ? 0 : -rmax*norm.y()/mag; 249 dists = sinSphi*topx - cosSphi*topy; 250 diste = -sinEphi*topx + cosEphi*topy; 251 if (dphi > pi) 252 { 253 iftop = true; 254 if (dists > 0 && diste > 0) iftop = false; 255 } 256 else 257 { 258 iftop = false; 259 if (dists <= 0 && diste <= 0) iftop = true 260 } 261 if (iftop) 262 { 263 zmax = -(norm.x()*topx + norm.y()*topy)/no 264 } 265 else 266 { 267 G4double z1 = -rmin*(norm.x()*cosSphi + no 268 G4double z2 = -rmin*(norm.x()*cosEphi + no 269 G4double z3 = -rmax*(norm.x()*cosSphi + no 270 G4double z4 = -rmax*(norm.x()*cosEphi + no 271 zmax = std::max(std::max(std::max(z1,z2),z 272 } 273 274 // Find bounding box 275 // 276 if (GetDeltaPhiAngle() < twopi) 277 { 278 G4TwoVector vmin,vmax; 279 G4GeomTools::DiskExtent(rmin,rmax, 280 GetSinStartPhi(),G 281 GetSinEndPhi(),Get 282 vmin,vmax); 283 pMin.set(vmin.x(),vmin.y(), zmin); 284 pMax.set(vmax.x(),vmax.y(), zmax); 285 } 286 else 287 { 288 pMin.set(-rmax,-rmax, zmin); 289 pMax.set( rmax, rmax, zmax); 290 } 291 292 // Check correctness of the bounding box 293 // 294 if (pMin.x() >= pMax.x() || pMin.y() >= pMax 295 { 296 std::ostringstream message; 297 message << "Bad bounding box (min >= max) 298 << GetName() << " !" 299 << "\npMin = " << pMin 300 << "\npMax = " << pMax; 301 G4Exception("G4CUutTubs::BoundingLimits()" 302 JustWarning, message); 303 StreamInfo(G4cout); 304 } 305 306 // Check consistency of bounding boxes 307 // 308 if (checkBBox) 309 { 310 U3Vector vmin, vmax; 311 Extent(vmin,vmax); 312 if (std::abs(pMin.x()-vmin.x()) > kCarTole 313 std::abs(pMin.y()-vmin.y()) > kCarTole 314 std::abs(pMin.z()-vmin.z()) > kCarTole 315 std::abs(pMax.x()-vmax.x()) > kCarTole 316 std::abs(pMax.y()-vmax.y()) > kCarTole 317 std::abs(pMax.z()-vmax.z()) > kCarTole 318 { 319 std::ostringstream message; 320 message << "Inconsistency in bounding bo 321 << GetName() << " !" 322 << "\nBBox min: wrapper = " << p 323 << "\nBBox max: wrapper = " << p 324 G4Exception("G4UCutTubs::BoundingLimits( 325 JustWarning, message); 326 checkBBox = false; 327 } 328 } 329 } 330 331 ////////////////////////////////////////////// 332 // 333 // Calculate extent under transform and specif 334 335 G4bool 336 G4UCutTubs::CalculateExtent(const EAxis pAxis, 337 const G4VoxelLimit 338 const G4AffineTran 339 G4double& pM 340 { 341 G4ThreeVector bmin, bmax; 342 G4bool exist; 343 344 // Get bounding box 345 BoundingLimits(bmin,bmax); 346 347 // Check bounding box 348 G4BoundingEnvelope bbox(bmin,bmax); 349 #ifdef G4BBOX_EXTENT 350 if (true) return bbox.CalculateExtent(pAxis, 351 #endif 352 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox 353 { 354 return exist = pMin < pMax; 355 } 356 357 // Get parameters of the solid 358 G4double rmin = GetInnerRadius(); 359 G4double rmax = GetOuterRadius(); 360 G4double dphi = GetDeltaPhiAngle(); 361 G4double zmin = bmin.z(); 362 G4double zmax = bmax.z(); 363 364 // Find bounding envelope and calculate exte 365 // 366 const G4int NSTEPS = 24; // numbe 367 G4double astep = twopi/NSTEPS; // max a 368 G4int ksteps = (dphi <= astep) ? 1 : (G4i 369 G4double ang = dphi/ksteps; 370 371 G4double sinHalf = std::sin(0.5*ang); 372 G4double cosHalf = std::cos(0.5*ang); 373 G4double sinStep = 2.*sinHalf*cosHalf; 374 G4double cosStep = 1. - 2.*sinHalf*sinHalf; 375 G4double rext = rmax/cosHalf; 376 377 // bounding envelope for full cylinder consi 378 // in other cases it is a sequence of quadri 379 if (rmin == 0 && dphi == twopi) 380 { 381 G4double sinCur = sinHalf; 382 G4double cosCur = cosHalf; 383 384 G4ThreeVectorList baseA(NSTEPS),baseB(NSTE 385 for (G4int k=0; k<NSTEPS; ++k) 386 { 387 baseA[k].set(rext*cosCur,rext*sinCur,zmi 388 baseB[k].set(rext*cosCur,rext*sinCur,zma 389 390 G4double sinTmp = sinCur; 391 sinCur = sinCur*cosStep + cosCur*sinStep 392 cosCur = cosCur*cosStep - sinTmp*sinStep 393 } 394 std::vector<const G4ThreeVectorList *> pol 395 polygons[0] = &baseA; 396 polygons[1] = &baseB; 397 G4BoundingEnvelope benv(bmin,bmax,polygons 398 exist = benv.CalculateExtent(pAxis,pVoxelL 399 } 400 else 401 { 402 G4double sinStart = GetSinStartPhi(); 403 G4double cosStart = GetCosStartPhi(); 404 G4double sinEnd = GetSinEndPhi(); 405 G4double cosEnd = GetCosEndPhi(); 406 G4double sinCur = sinStart*cosHalf + cos 407 G4double cosCur = cosStart*cosHalf - sin 408 409 // set quadrilaterals 410 G4ThreeVectorList pols[NSTEPS+2]; 411 for (G4int k=0; k<ksteps+2; ++k) pols[k].r 412 pols[0][0].set(rmin*cosStart,rmin*sinStart 413 pols[0][1].set(rmin*cosStart,rmin*sinStart 414 pols[0][2].set(rmax*cosStart,rmax*sinStart 415 pols[0][3].set(rmax*cosStart,rmax*sinStart 416 for (G4int k=1; k<ksteps+1; ++k) 417 { 418 pols[k][0].set(rmin*cosCur,rmin*sinCur,z 419 pols[k][1].set(rmin*cosCur,rmin*sinCur,z 420 pols[k][2].set(rext*cosCur,rext*sinCur,z 421 pols[k][3].set(rext*cosCur,rext*sinCur,z 422 423 G4double sinTmp = sinCur; 424 sinCur = sinCur*cosStep + cosCur*sinStep 425 cosCur = cosCur*cosStep - sinTmp*sinStep 426 } 427 pols[ksteps+1][0].set(rmin*cosEnd,rmin*sin 428 pols[ksteps+1][1].set(rmin*cosEnd,rmin*sin 429 pols[ksteps+1][2].set(rmax*cosEnd,rmax*sin 430 pols[ksteps+1][3].set(rmax*cosEnd,rmax*sin 431 432 // set envelope and calculate extent 433 std::vector<const G4ThreeVectorList *> pol 434 polygons.resize(ksteps+2); 435 for (G4int k=0; k<ksteps+2; ++k) polygons[ 436 G4BoundingEnvelope benv(bmin,bmax,polygons 437 exist = benv.CalculateExtent(pAxis,pVoxelL 438 } 439 return exist; 440 } 441 442 ////////////////////////////////////////////// 443 // 444 // Return real Z coordinate of point on Cutted 445 446 G4double G4UCutTubs::GetCutZ(const G4ThreeVect 447 { 448 G4double newz = p.z(); // p.z() should be e 449 G4ThreeVector fLowNorm = GetLowNorm(); 450 G4ThreeVector fHighNorm = GetHighNorm(); 451 452 if (p.z()<0) 453 { 454 if(fLowNorm.z()!=0.) 455 { 456 newz = -GetZHalfLength() 457 - (p.x()*fLowNorm.x()+p.y()*fLowNo 458 } 459 } 460 else 461 { 462 if(fHighNorm.z()!=0.) 463 { 464 newz = GetZHalfLength() 465 - (p.x()*fHighNorm.x()+p.y()*fHigh 466 } 467 } 468 return newz; 469 } 470 471 ////////////////////////////////////////////// 472 // 473 // Create polyhedron for visualization 474 // 475 G4Polyhedron* G4UCutTubs::CreatePolyhedron() c 476 { 477 typedef G4double G4double3[3]; 478 typedef G4int G4int4[4]; 479 480 auto ph = new G4Polyhedron; 481 G4Polyhedron *ph1 = new G4PolyhedronTubs(Get 482 Get 483 Get 484 Get 485 Get 486 G4int nn=ph1->GetNoVertices(); 487 G4int nf=ph1->GetNoFacets(); 488 auto xyz = new G4double3[nn]; // number of 489 auto faces = new G4int4[nf] ; // number of 490 G4double fDz = GetZHalfLength(); 491 492 for(G4int i=0; i<nn; ++i) 493 { 494 xyz[i][0]=ph1->GetVertex(i+1).x(); 495 xyz[i][1]=ph1->GetVertex(i+1).y(); 496 G4double tmpZ=ph1->GetVertex(i+1).z(); 497 if (tmpZ>=fDz-kCarTolerance) 498 { 499 xyz[i][2]=GetCutZ(G4ThreeVector(xyz[i][0 500 } 501 else if(tmpZ<=-fDz+kCarTolerance) 502 { 503 xyz[i][2]=GetCutZ(G4ThreeVector(xyz[i][0 504 } 505 else 506 { 507 xyz[i][2]=tmpZ; 508 } 509 } 510 G4int iNodes[4]; 511 G4int* iEdge=nullptr; 512 G4int n; 513 for(G4int i=0; i<nf; ++i) 514 { 515 ph1->GetFacet(i+1,n,iNodes,iEdge); 516 for(G4int k=0; k<n; ++k) 517 { 518 faces[i][k]=iNodes[k]; 519 } 520 for(G4int k=n; k<4; ++k) 521 { 522 faces[i][k]=0; 523 } 524 } 525 ph->createPolyhedron(nn,nf,xyz,faces); 526 527 delete [] xyz; 528 delete [] faces; 529 delete ph1; 530 531 return ph; 532 } 533 534 #endif // G4GEOM_USE_USOLIDS 535