Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // Implementation for G4Trd class 26 // Implementation for G4Trd class 27 // 27 // 28 // 12.01.95 P.Kent: First version 28 // 12.01.95 P.Kent: First version 29 // 28.04.05 V.Grichine: new SurfaceNormal acco 29 // 28.04.05 V.Grichine: new SurfaceNormal according to J.Apostolakis proposal 30 // 25.05.17 E.Tcherniaev: complete revision, s 30 // 25.05.17 E.Tcherniaev: complete revision, speed-up 31 // ------------------------------------------- 31 // -------------------------------------------------------------------- 32 32 33 #include "G4Trd.hh" 33 #include "G4Trd.hh" 34 34 35 #if !defined(G4GEOM_USE_UTRD) 35 #if !defined(G4GEOM_USE_UTRD) 36 36 37 #include "G4GeomTools.hh" 37 #include "G4GeomTools.hh" 38 38 39 #include "G4VoxelLimits.hh" 39 #include "G4VoxelLimits.hh" 40 #include "G4AffineTransform.hh" 40 #include "G4AffineTransform.hh" 41 #include "G4BoundingEnvelope.hh" 41 #include "G4BoundingEnvelope.hh" 42 #include "G4QuickRand.hh" 42 #include "G4QuickRand.hh" 43 43 44 #include "G4VPVParameterisation.hh" 44 #include "G4VPVParameterisation.hh" 45 45 46 #include "G4VGraphicsScene.hh" 46 #include "G4VGraphicsScene.hh" 47 47 48 using namespace CLHEP; 48 using namespace CLHEP; 49 49 50 ////////////////////////////////////////////// 50 ////////////////////////////////////////////////////////////////////////// 51 // 51 // 52 // Constructor - set & check half widths 52 // Constructor - set & check half widths 53 53 54 G4Trd::G4Trd(const G4String& pName, 54 G4Trd::G4Trd(const G4String& pName, 55 G4double pdx1, G4double pdx 55 G4double pdx1, G4double pdx2, 56 G4double pdy1, G4double pdy 56 G4double pdy1, G4double pdy2, 57 G4double pdz) 57 G4double pdz) 58 : G4CSGSolid(pName), halfCarTolerance(0.5*kC 58 : G4CSGSolid(pName), halfCarTolerance(0.5*kCarTolerance), 59 fDx1(pdx1), fDx2(pdx2), fDy1(pdy1), fDy2(p 59 fDx1(pdx1), fDx2(pdx2), fDy1(pdy1), fDy2(pdy2), fDz(pdz) 60 { 60 { 61 CheckParameters(); 61 CheckParameters(); 62 MakePlanes(); 62 MakePlanes(); 63 } 63 } 64 64 65 ////////////////////////////////////////////// 65 ////////////////////////////////////////////////////////////////////////// 66 // 66 // 67 // Fake default constructor - sets only member 67 // Fake default constructor - sets only member data and allocates memory 68 // for usage restri 68 // for usage restricted to object persistency 69 // 69 // 70 G4Trd::G4Trd( __void__& a ) 70 G4Trd::G4Trd( __void__& a ) 71 : G4CSGSolid(a), halfCarTolerance(0.5*kCarTo 71 : G4CSGSolid(a), halfCarTolerance(0.5*kCarTolerance), 72 fDx1(1.), fDx2(1.), fDy1(1.), fDy2(1.), fD 72 fDx1(1.), fDx2(1.), fDy1(1.), fDy2(1.), fDz(1.) 73 { 73 { 74 MakePlanes(); 74 MakePlanes(); 75 } 75 } 76 76 77 ////////////////////////////////////////////// 77 ////////////////////////////////////////////////////////////////////////// 78 // 78 // 79 // Destructor 79 // Destructor 80 80 81 G4Trd::~G4Trd() = default; << 81 G4Trd::~G4Trd() >> 82 { >> 83 } 82 84 83 ////////////////////////////////////////////// 85 ////////////////////////////////////////////////////////////////////////// 84 // 86 // 85 // Copy constructor 87 // Copy constructor 86 88 87 G4Trd::G4Trd(const G4Trd& rhs) 89 G4Trd::G4Trd(const G4Trd& rhs) 88 : G4CSGSolid(rhs), halfCarTolerance(rhs.half 90 : G4CSGSolid(rhs), halfCarTolerance(rhs.halfCarTolerance), 89 fDx1(rhs.fDx1), fDx2(rhs.fDx2), 91 fDx1(rhs.fDx1), fDx2(rhs.fDx2), 90 fDy1(rhs.fDy1), fDy2(rhs.fDy2), fDz(rhs.fD 92 fDy1(rhs.fDy1), fDy2(rhs.fDy2), fDz(rhs.fDz), 91 fHx(rhs.fHx), fHy(rhs.fHy) 93 fHx(rhs.fHx), fHy(rhs.fHy) 92 { 94 { 93 for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs 95 for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs.fPlanes[i]; } 94 } 96 } 95 97 96 ////////////////////////////////////////////// 98 ////////////////////////////////////////////////////////////////////////// 97 // 99 // 98 // Assignment operator 100 // Assignment operator 99 101 100 G4Trd& G4Trd::operator = (const G4Trd& rhs) 102 G4Trd& G4Trd::operator = (const G4Trd& rhs) 101 { 103 { 102 // Check assignment to self 104 // Check assignment to self 103 // 105 // 104 if (this == &rhs) { return *this; } 106 if (this == &rhs) { return *this; } 105 107 106 // Copy base class data 108 // Copy base class data 107 // 109 // 108 G4CSGSolid::operator=(rhs); 110 G4CSGSolid::operator=(rhs); 109 111 110 // Copy data 112 // Copy data 111 // 113 // 112 halfCarTolerance = rhs.halfCarTolerance; 114 halfCarTolerance = rhs.halfCarTolerance; 113 fDx1 = rhs.fDx1; fDx2 = rhs.fDx2; 115 fDx1 = rhs.fDx1; fDx2 = rhs.fDx2; 114 fDy1 = rhs.fDy1; fDy2 = rhs.fDy2; 116 fDy1 = rhs.fDy1; fDy2 = rhs.fDy2; 115 fDz = rhs.fDz; 117 fDz = rhs.fDz; 116 fHx = rhs.fHx; fHy = rhs.fHy; 118 fHx = rhs.fHx; fHy = rhs.fHy; 117 for (G4int i=0; i<4; ++i) { fPlanes[i] = rh 119 for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs.fPlanes[i]; } 118 120 119 return *this; 121 return *this; 120 } 122 } 121 123 122 ////////////////////////////////////////////// 124 ////////////////////////////////////////////////////////////////////////// 123 // 125 // 124 // Set all parameters, as for constructor - se 126 // Set all parameters, as for constructor - set and check half-widths 125 127 126 void G4Trd::SetAllParameters(G4double pdx1, G4 128 void G4Trd::SetAllParameters(G4double pdx1, G4double pdx2, 127 G4double pdy1, G4 129 G4double pdy1, G4double pdy2, G4double pdz) 128 { 130 { 129 // Reset data of the base class 131 // Reset data of the base class 130 fCubicVolume = 0.; 132 fCubicVolume = 0.; 131 fSurfaceArea = 0.; 133 fSurfaceArea = 0.; 132 fRebuildPolyhedron = true; 134 fRebuildPolyhedron = true; 133 135 134 // Set parameters 136 // Set parameters 135 fDx1 = pdx1; fDx2 = pdx2; 137 fDx1 = pdx1; fDx2 = pdx2; 136 fDy1 = pdy1; fDy2 = pdy2; 138 fDy1 = pdy1; fDy2 = pdy2; 137 fDz = pdz; 139 fDz = pdz; 138 140 139 CheckParameters(); 141 CheckParameters(); 140 MakePlanes(); 142 MakePlanes(); 141 } 143 } 142 144 143 ////////////////////////////////////////////// 145 ////////////////////////////////////////////////////////////////////////// 144 // 146 // 145 // Check dimensions 147 // Check dimensions 146 148 147 void G4Trd::CheckParameters() 149 void G4Trd::CheckParameters() 148 { 150 { 149 G4double dmin = 2*kCarTolerance; 151 G4double dmin = 2*kCarTolerance; 150 if ((fDx1 < 0 || fDx2 < 0 || fDy1 < 0 || fDy 152 if ((fDx1 < 0 || fDx2 < 0 || fDy1 < 0 || fDy2 < 0 || fDz < dmin) || 151 (fDx1 < dmin && fDx2 < dmin) || 153 (fDx1 < dmin && fDx2 < dmin) || 152 (fDy1 < dmin && fDy2 < dmin)) 154 (fDy1 < dmin && fDy2 < dmin)) 153 { 155 { 154 std::ostringstream message; 156 std::ostringstream message; 155 message << "Invalid (too small or negative 157 message << "Invalid (too small or negative) dimensions for Solid: " 156 << GetName() 158 << GetName() 157 << "\n X - " << fDx1 << ", " << f 159 << "\n X - " << fDx1 << ", " << fDx2 158 << "\n Y - " << fDy1 << ", " << f 160 << "\n Y - " << fDy1 << ", " << fDy2 159 << "\n Z - " << fDz; 161 << "\n Z - " << fDz; 160 G4Exception("G4Trd::CheckParameters()", "G 162 G4Exception("G4Trd::CheckParameters()", "GeomSolids0002", 161 FatalException, message); 163 FatalException, message); 162 } 164 } 163 } 165 } 164 166 165 ////////////////////////////////////////////// 167 ////////////////////////////////////////////////////////////////////////// 166 // 168 // 167 // Set side planes 169 // Set side planes 168 170 169 void G4Trd::MakePlanes() 171 void G4Trd::MakePlanes() 170 { 172 { 171 G4double dx = fDx1 - fDx2; 173 G4double dx = fDx1 - fDx2; 172 G4double dy = fDy1 - fDy2; 174 G4double dy = fDy1 - fDy2; 173 G4double dz = 2*fDz; 175 G4double dz = 2*fDz; 174 fHx = std::sqrt(dy*dy + dz*dz); 176 fHx = std::sqrt(dy*dy + dz*dz); 175 fHy = std::sqrt(dx*dx + dz*dz); 177 fHy = std::sqrt(dx*dx + dz*dz); 176 178 177 // Set X planes at -Y & +Y 179 // Set X planes at -Y & +Y 178 // 180 // 179 fPlanes[0].a = 0.; 181 fPlanes[0].a = 0.; 180 fPlanes[0].b = -dz/fHx; 182 fPlanes[0].b = -dz/fHx; 181 fPlanes[0].c = dy/fHx; 183 fPlanes[0].c = dy/fHx; 182 fPlanes[0].d = fPlanes[0].b*fDy1 + fPlanes[0 184 fPlanes[0].d = fPlanes[0].b*fDy1 + fPlanes[0].c*fDz; 183 185 184 fPlanes[1].a = fPlanes[0].a; 186 fPlanes[1].a = fPlanes[0].a; 185 fPlanes[1].b = -fPlanes[0].b; 187 fPlanes[1].b = -fPlanes[0].b; 186 fPlanes[1].c = fPlanes[0].c; 188 fPlanes[1].c = fPlanes[0].c; 187 fPlanes[1].d = fPlanes[0].d; 189 fPlanes[1].d = fPlanes[0].d; 188 190 189 // Set Y planes at -X & +X 191 // Set Y planes at -X & +X 190 // 192 // 191 fPlanes[2].a = -dz/fHy; 193 fPlanes[2].a = -dz/fHy; 192 fPlanes[2].b = 0.; 194 fPlanes[2].b = 0.; 193 fPlanes[2].c = dx/fHy; 195 fPlanes[2].c = dx/fHy; 194 fPlanes[2].d = fPlanes[2].a*fDx1 + fPlanes[2 196 fPlanes[2].d = fPlanes[2].a*fDx1 + fPlanes[2].c*fDz; 195 197 196 fPlanes[3].a = -fPlanes[2].a; 198 fPlanes[3].a = -fPlanes[2].a; 197 fPlanes[3].b = fPlanes[2].b; 199 fPlanes[3].b = fPlanes[2].b; 198 fPlanes[3].c = fPlanes[2].c; 200 fPlanes[3].c = fPlanes[2].c; 199 fPlanes[3].d = fPlanes[2].d; 201 fPlanes[3].d = fPlanes[2].d; 200 } 202 } 201 203 202 ////////////////////////////////////////////// 204 ////////////////////////////////////////////////////////////////////////// 203 // 205 // 204 // Get volume 206 // Get volume 205 207 206 G4double G4Trd::GetCubicVolume() 208 G4double G4Trd::GetCubicVolume() 207 { 209 { 208 if (fCubicVolume == 0.) 210 if (fCubicVolume == 0.) 209 { 211 { 210 fCubicVolume = 2*fDz*( (fDx1+fDx2)*(fDy1+f 212 fCubicVolume = 2*fDz*( (fDx1+fDx2)*(fDy1+fDy2) + 211 (fDx2-fDx1)*(fDy2-f 213 (fDx2-fDx1)*(fDy2-fDy1)/3 ); 212 } 214 } 213 return fCubicVolume; 215 return fCubicVolume; 214 } 216 } 215 217 216 ////////////////////////////////////////////// 218 ////////////////////////////////////////////////////////////////////////// 217 // 219 // 218 // Get surface area 220 // Get surface area 219 221 220 G4double G4Trd::GetSurfaceArea() 222 G4double G4Trd::GetSurfaceArea() 221 { 223 { 222 if (fSurfaceArea == 0.) 224 if (fSurfaceArea == 0.) 223 { 225 { 224 fSurfaceArea = 226 fSurfaceArea = 225 4*(fDx1*fDy1 + fDx2*fDy2) + 2*(fDx1+fDx2 227 4*(fDx1*fDy1 + fDx2*fDy2) + 2*(fDx1+fDx2)*fHx + 2*(fDy1+fDy2)*fHy; 226 } 228 } 227 return fSurfaceArea; 229 return fSurfaceArea; 228 } 230 } 229 231 230 ////////////////////////////////////////////// 232 ////////////////////////////////////////////////////////////////////////// 231 // 233 // 232 // Dispatch to parameterisation for replicatio 234 // Dispatch to parameterisation for replication mechanism dimension 233 // computation & modification 235 // computation & modification 234 236 235 void G4Trd::ComputeDimensions( G4VPVPara 237 void G4Trd::ComputeDimensions( G4VPVParameterisation* p, 236 const G4int n, 238 const G4int n, 237 const G4VPhysic 239 const G4VPhysicalVolume* pRep ) 238 { 240 { 239 p->ComputeDimensions(*this,n,pRep); 241 p->ComputeDimensions(*this,n,pRep); 240 } 242 } 241 243 242 ////////////////////////////////////////////// 244 ////////////////////////////////////////////////////////////////////////// 243 // 245 // 244 // Get bounding box 246 // Get bounding box 245 247 246 void G4Trd::BoundingLimits(G4ThreeVector& pMin 248 void G4Trd::BoundingLimits(G4ThreeVector& pMin, G4ThreeVector& pMax) const 247 { 249 { 248 G4double dx1 = GetXHalfLength1(); 250 G4double dx1 = GetXHalfLength1(); 249 G4double dx2 = GetXHalfLength2(); 251 G4double dx2 = GetXHalfLength2(); 250 G4double dy1 = GetYHalfLength1(); 252 G4double dy1 = GetYHalfLength1(); 251 G4double dy2 = GetYHalfLength2(); 253 G4double dy2 = GetYHalfLength2(); 252 G4double dz = GetZHalfLength(); 254 G4double dz = GetZHalfLength(); 253 255 254 G4double xmax = std::max(dx1,dx2); 256 G4double xmax = std::max(dx1,dx2); 255 G4double ymax = std::max(dy1,dy2); 257 G4double ymax = std::max(dy1,dy2); 256 pMin.set(-xmax,-ymax,-dz); 258 pMin.set(-xmax,-ymax,-dz); 257 pMax.set( xmax, ymax, dz); 259 pMax.set( xmax, ymax, dz); 258 260 259 // Check correctness of the bounding box 261 // Check correctness of the bounding box 260 // 262 // 261 if (pMin.x() >= pMax.x() || pMin.y() >= pMax 263 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z()) 262 { 264 { 263 std::ostringstream message; 265 std::ostringstream message; 264 message << "Bad bounding box (min >= max) 266 message << "Bad bounding box (min >= max) for solid: " 265 << GetName() << " !" 267 << GetName() << " !" 266 << "\npMin = " << pMin 268 << "\npMin = " << pMin 267 << "\npMax = " << pMax; 269 << "\npMax = " << pMax; 268 G4Exception("G4Trd::BoundingLimits()", "Ge 270 G4Exception("G4Trd::BoundingLimits()", "GeomMgt0001", JustWarning, message); 269 DumpInfo(); 271 DumpInfo(); 270 } 272 } 271 } 273 } 272 274 273 ////////////////////////////////////////////// 275 ////////////////////////////////////////////////////////////////////////// 274 // 276 // 275 // Calculate extent under transform and specif 277 // Calculate extent under transform and specified limit 276 278 277 G4bool G4Trd::CalculateExtent( const EAxis pAx 279 G4bool G4Trd::CalculateExtent( const EAxis pAxis, 278 const G4VoxelLi 280 const G4VoxelLimits& pVoxelLimit, 279 const G4AffineT 281 const G4AffineTransform& pTransform, 280 G4double& 282 G4double& pMin, G4double& pMax ) const 281 { 283 { 282 G4ThreeVector bmin, bmax; 284 G4ThreeVector bmin, bmax; 283 G4bool exist; 285 G4bool exist; 284 286 285 // Check bounding box (bbox) 287 // Check bounding box (bbox) 286 // 288 // 287 BoundingLimits(bmin,bmax); 289 BoundingLimits(bmin,bmax); 288 G4BoundingEnvelope bbox(bmin,bmax); 290 G4BoundingEnvelope bbox(bmin,bmax); 289 #ifdef G4BBOX_EXTENT 291 #ifdef G4BBOX_EXTENT 290 return bbox.CalculateExtent(pAxis,pVoxelLimi 292 return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax); 291 #endif 293 #endif 292 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox 294 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax)) 293 { 295 { 294 return exist = pMin < pMax; << 296 return exist = (pMin < pMax) ? true : false; 295 } 297 } 296 298 297 // Set bounding envelope (benv) and calculat 299 // Set bounding envelope (benv) and calculate extent 298 // 300 // 299 G4double dx1 = GetXHalfLength1(); 301 G4double dx1 = GetXHalfLength1(); 300 G4double dx2 = GetXHalfLength2(); 302 G4double dx2 = GetXHalfLength2(); 301 G4double dy1 = GetYHalfLength1(); 303 G4double dy1 = GetYHalfLength1(); 302 G4double dy2 = GetYHalfLength2(); 304 G4double dy2 = GetYHalfLength2(); 303 G4double dz = GetZHalfLength(); 305 G4double dz = GetZHalfLength(); 304 306 305 G4ThreeVectorList baseA(4), baseB(4); 307 G4ThreeVectorList baseA(4), baseB(4); 306 baseA[0].set(-dx1,-dy1,-dz); 308 baseA[0].set(-dx1,-dy1,-dz); 307 baseA[1].set( dx1,-dy1,-dz); 309 baseA[1].set( dx1,-dy1,-dz); 308 baseA[2].set( dx1, dy1,-dz); 310 baseA[2].set( dx1, dy1,-dz); 309 baseA[3].set(-dx1, dy1,-dz); 311 baseA[3].set(-dx1, dy1,-dz); 310 baseB[0].set(-dx2,-dy2, dz); 312 baseB[0].set(-dx2,-dy2, dz); 311 baseB[1].set( dx2,-dy2, dz); 313 baseB[1].set( dx2,-dy2, dz); 312 baseB[2].set( dx2, dy2, dz); 314 baseB[2].set( dx2, dy2, dz); 313 baseB[3].set(-dx2, dy2, dz); 315 baseB[3].set(-dx2, dy2, dz); 314 316 315 std::vector<const G4ThreeVectorList *> polyg 317 std::vector<const G4ThreeVectorList *> polygons(2); 316 polygons[0] = &baseA; 318 polygons[0] = &baseA; 317 polygons[1] = &baseB; 319 polygons[1] = &baseB; 318 320 319 G4BoundingEnvelope benv(bmin,bmax,polygons); 321 G4BoundingEnvelope benv(bmin,bmax,polygons); 320 exist = benv.CalculateExtent(pAxis,pVoxelLim 322 exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax); 321 return exist; 323 return exist; 322 } 324 } 323 325 324 ////////////////////////////////////////////// 326 ////////////////////////////////////////////////////////////////////////// 325 // 327 // 326 // Return whether point inside/outside/on surf 328 // Return whether point inside/outside/on surface, using tolerance 327 329 328 EInside G4Trd::Inside( const G4ThreeVector& p 330 EInside G4Trd::Inside( const G4ThreeVector& p ) const 329 { 331 { 330 G4double dx = fPlanes[3].a*std::abs(p.x())+f 332 G4double dx = fPlanes[3].a*std::abs(p.x())+fPlanes[3].c*p.z()+fPlanes[3].d; 331 G4double dy = fPlanes[1].b*std::abs(p.y())+f 333 G4double dy = fPlanes[1].b*std::abs(p.y())+fPlanes[1].c*p.z()+fPlanes[1].d; 332 G4double dxy = std::max(dx,dy); 334 G4double dxy = std::max(dx,dy); 333 335 334 G4double dz = std::abs(p.z())-fDz; 336 G4double dz = std::abs(p.z())-fDz; 335 G4double dist = std::max(dz,dxy); 337 G4double dist = std::max(dz,dxy); 336 338 337 return (dist > halfCarTolerance) ? kOutside 339 return (dist > halfCarTolerance) ? kOutside : 338 ((dist > -halfCarTolerance) ? kSurface : k 340 ((dist > -halfCarTolerance) ? kSurface : kInside); 339 } 341 } 340 342 341 ////////////////////////////////////////////// 343 ////////////////////////////////////////////////////////////////////////// 342 // 344 // 343 // Determine side where point is, and return c 345 // Determine side where point is, and return corresponding normal 344 346 345 G4ThreeVector G4Trd::SurfaceNormal( const G4Th 347 G4ThreeVector G4Trd::SurfaceNormal( const G4ThreeVector& p ) const 346 { 348 { 347 G4int nsurf = 0; // number of surfaces where 349 G4int nsurf = 0; // number of surfaces where p is placed 348 350 349 // Check Z faces 351 // Check Z faces 350 // 352 // 351 G4double nz = 0; 353 G4double nz = 0; 352 G4double dz = std::abs(p.z()) - fDz; 354 G4double dz = std::abs(p.z()) - fDz; 353 if (std::abs(dz) <= halfCarTolerance) 355 if (std::abs(dz) <= halfCarTolerance) 354 { 356 { 355 nz = (p.z() < 0) ? -1 : 1; 357 nz = (p.z() < 0) ? -1 : 1; 356 ++nsurf; 358 ++nsurf; 357 } 359 } 358 360 359 // Check Y faces 361 // Check Y faces 360 // 362 // 361 G4double ny = 0; 363 G4double ny = 0; 362 G4double dy1 = fPlanes[0].b*p.y(); 364 G4double dy1 = fPlanes[0].b*p.y(); 363 G4double dy2 = fPlanes[0].c*p.z() + fPlanes[ 365 G4double dy2 = fPlanes[0].c*p.z() + fPlanes[0].d; 364 if (std::abs(dy2 + dy1) <= halfCarTolerance) 366 if (std::abs(dy2 + dy1) <= halfCarTolerance) 365 { 367 { 366 ny += fPlanes[0].b; 368 ny += fPlanes[0].b; 367 nz += fPlanes[0].c; 369 nz += fPlanes[0].c; 368 ++nsurf; 370 ++nsurf; 369 } 371 } 370 if (std::abs(dy2 - dy1) <= halfCarTolerance) 372 if (std::abs(dy2 - dy1) <= halfCarTolerance) 371 { 373 { 372 ny += fPlanes[1].b; 374 ny += fPlanes[1].b; 373 nz += fPlanes[1].c; 375 nz += fPlanes[1].c; 374 ++nsurf; 376 ++nsurf; 375 } 377 } 376 378 377 // Check X faces 379 // Check X faces 378 // 380 // 379 G4double nx = 0; 381 G4double nx = 0; 380 G4double dx1 = fPlanes[2].a*p.x(); 382 G4double dx1 = fPlanes[2].a*p.x(); 381 G4double dx2 = fPlanes[2].c*p.z() + fPlanes[ 383 G4double dx2 = fPlanes[2].c*p.z() + fPlanes[2].d; 382 if (std::abs(dx2 + dx1) <= halfCarTolerance) 384 if (std::abs(dx2 + dx1) <= halfCarTolerance) 383 { 385 { 384 nx += fPlanes[2].a; 386 nx += fPlanes[2].a; 385 nz += fPlanes[2].c; 387 nz += fPlanes[2].c; 386 ++nsurf; 388 ++nsurf; 387 } 389 } 388 if (std::abs(dx2 - dx1) <= halfCarTolerance) 390 if (std::abs(dx2 - dx1) <= halfCarTolerance) 389 { 391 { 390 nx += fPlanes[3].a; 392 nx += fPlanes[3].a; 391 nz += fPlanes[3].c; 393 nz += fPlanes[3].c; 392 ++nsurf; 394 ++nsurf; 393 } 395 } 394 396 395 // Return normal 397 // Return normal 396 // 398 // 397 if (nsurf == 1) return {nx,ny,nz}; << 399 if (nsurf == 1) return G4ThreeVector(nx,ny,nz); 398 else if (nsurf != 0) return G4ThreeVector(nx 400 else if (nsurf != 0) return G4ThreeVector(nx,ny,nz).unit(); // edge or corner 399 else 401 else 400 { 402 { 401 // Point is not on the surface 403 // Point is not on the surface 402 // 404 // 403 #ifdef G4CSGDEBUG 405 #ifdef G4CSGDEBUG 404 std::ostringstream message; 406 std::ostringstream message; 405 G4long oldprc = message.precision(16); << 407 G4int oldprc = message.precision(16); 406 message << "Point p is not on surface (!?) 408 message << "Point p is not on surface (!?) of solid: " 407 << GetName() << G4endl; 409 << GetName() << G4endl; 408 message << "Position:\n"; 410 message << "Position:\n"; 409 message << " p.x() = " << p.x()/mm << " 411 message << " p.x() = " << p.x()/mm << " mm\n"; 410 message << " p.y() = " << p.y()/mm << " 412 message << " p.y() = " << p.y()/mm << " mm\n"; 411 message << " p.z() = " << p.z()/mm << " 413 message << " p.z() = " << p.z()/mm << " mm"; 412 G4cout.precision(oldprc) ; 414 G4cout.precision(oldprc) ; 413 G4Exception("G4Trd::SurfaceNormal(p)", "Ge 415 G4Exception("G4Trd::SurfaceNormal(p)", "GeomSolids1002", 414 JustWarning, message ); 416 JustWarning, message ); 415 DumpInfo(); 417 DumpInfo(); 416 #endif 418 #endif 417 return ApproxSurfaceNormal(p); 419 return ApproxSurfaceNormal(p); 418 } 420 } 419 } 421 } 420 422 421 ////////////////////////////////////////////// 423 ////////////////////////////////////////////////////////////////////////// 422 // 424 // 423 // Algorithm for SurfaceNormal() following the 425 // Algorithm for SurfaceNormal() following the original specification 424 // for points not on the surface 426 // for points not on the surface 425 427 426 G4ThreeVector G4Trd::ApproxSurfaceNormal( cons 428 G4ThreeVector G4Trd::ApproxSurfaceNormal( const G4ThreeVector& p ) const 427 { 429 { 428 G4double dist = -DBL_MAX; 430 G4double dist = -DBL_MAX; 429 G4int iside = 0; 431 G4int iside = 0; 430 for (G4int i=0; i<4; ++i) 432 for (G4int i=0; i<4; ++i) 431 { 433 { 432 G4double d = fPlanes[i].a*p.x() + 434 G4double d = fPlanes[i].a*p.x() + 433 fPlanes[i].b*p.y() + 435 fPlanes[i].b*p.y() + 434 fPlanes[i].c*p.z() + fPlanes[ 436 fPlanes[i].c*p.z() + fPlanes[i].d; 435 if (d > dist) { dist = d; iside = i; } 437 if (d > dist) { dist = d; iside = i; } 436 } 438 } 437 439 438 G4double distz = std::abs(p.z()) - fDz; 440 G4double distz = std::abs(p.z()) - fDz; 439 if (dist > distz) 441 if (dist > distz) 440 return { fPlanes[iside].a, fPlanes[iside]. << 442 return G4ThreeVector(fPlanes[iside].a, fPlanes[iside].b, fPlanes[iside].c); 441 else 443 else 442 return { 0, 0, (G4double)((p.z() < 0) ? -1 << 444 return G4ThreeVector(0, 0, (p.z() < 0) ? -1 : 1); 443 } 445 } 444 446 445 ////////////////////////////////////////////// 447 ////////////////////////////////////////////////////////////////////////// 446 // 448 // 447 // Calculate distance to shape from outside 449 // Calculate distance to shape from outside 448 // - return kInfinity if no intersection 450 // - return kInfinity if no intersection 449 451 450 G4double G4Trd::DistanceToIn(const G4ThreeVect 452 G4double G4Trd::DistanceToIn(const G4ThreeVector& p, 451 const G4ThreeVect 453 const G4ThreeVector& v ) const 452 { 454 { 453 // Z intersections 455 // Z intersections 454 // 456 // 455 if ((std::abs(p.z()) - fDz) >= -halfCarToler 457 if ((std::abs(p.z()) - fDz) >= -halfCarTolerance && p.z()*v.z() >= 0) 456 return kInfinity; 458 return kInfinity; 457 G4double invz = (-v.z() == 0) ? DBL_MAX : -1 459 G4double invz = (-v.z() == 0) ? DBL_MAX : -1./v.z(); 458 G4double dz = (invz < 0) ? fDz : -fDz; 460 G4double dz = (invz < 0) ? fDz : -fDz; 459 G4double tzmin = (p.z() + dz)*invz; 461 G4double tzmin = (p.z() + dz)*invz; 460 G4double tzmax = (p.z() - dz)*invz; 462 G4double tzmax = (p.z() - dz)*invz; 461 463 462 // Y intersections 464 // Y intersections 463 // 465 // 464 G4double tmin0 = tzmin, tmax0 = tzmax; 466 G4double tmin0 = tzmin, tmax0 = tzmax; 465 G4double ya = fPlanes[0].b*v.y(), yb = fPlan 467 G4double ya = fPlanes[0].b*v.y(), yb = fPlanes[0].c*v.z(); 466 G4double yc = fPlanes[0].b*p.y(), yd = fPlan 468 G4double yc = fPlanes[0].b*p.y(), yd = fPlanes[0].c*p.z()+fPlanes[0].d; 467 G4double cos0 = yb + ya; 469 G4double cos0 = yb + ya; 468 G4double dis0 = yd + yc; 470 G4double dis0 = yd + yc; 469 if (dis0 >= -halfCarTolerance) 471 if (dis0 >= -halfCarTolerance) 470 { 472 { 471 if (cos0 >= 0) return kInfinity; 473 if (cos0 >= 0) return kInfinity; 472 G4double tmp = -dis0/cos0; 474 G4double tmp = -dis0/cos0; 473 if (tmin0 < tmp) tmin0 = tmp; 475 if (tmin0 < tmp) tmin0 = tmp; 474 } 476 } 475 else if (cos0 > 0) 477 else if (cos0 > 0) 476 { 478 { 477 G4double tmp = -dis0/cos0; 479 G4double tmp = -dis0/cos0; 478 if (tmax0 > tmp) tmax0 = tmp; 480 if (tmax0 > tmp) tmax0 = tmp; 479 } 481 } 480 482 481 G4double tmin1 = tmin0, tmax1 = tmax0; 483 G4double tmin1 = tmin0, tmax1 = tmax0; 482 G4double cos1 = yb - ya; 484 G4double cos1 = yb - ya; 483 G4double dis1 = yd - yc; 485 G4double dis1 = yd - yc; 484 if (dis1 >= -halfCarTolerance) 486 if (dis1 >= -halfCarTolerance) 485 { 487 { 486 if (cos1 >= 0) return kInfinity; 488 if (cos1 >= 0) return kInfinity; 487 G4double tmp = -dis1/cos1; 489 G4double tmp = -dis1/cos1; 488 if (tmin1 < tmp) tmin1 = tmp; 490 if (tmin1 < tmp) tmin1 = tmp; 489 } 491 } 490 else if (cos1 > 0) 492 else if (cos1 > 0) 491 { 493 { 492 G4double tmp = -dis1/cos1; 494 G4double tmp = -dis1/cos1; 493 if (tmax1 > tmp) tmax1 = tmp; 495 if (tmax1 > tmp) tmax1 = tmp; 494 } 496 } 495 497 496 // X intersections 498 // X intersections 497 // 499 // 498 G4double tmin2 = tmin1, tmax2 = tmax1; 500 G4double tmin2 = tmin1, tmax2 = tmax1; 499 G4double xa = fPlanes[2].a*v.x(), xb = fPlan 501 G4double xa = fPlanes[2].a*v.x(), xb = fPlanes[2].c*v.z(); 500 G4double xc = fPlanes[2].a*p.x(), xd = fPlan 502 G4double xc = fPlanes[2].a*p.x(), xd = fPlanes[2].c*p.z()+fPlanes[2].d; 501 G4double cos2 = xb + xa; 503 G4double cos2 = xb + xa; 502 G4double dis2 = xd + xc; 504 G4double dis2 = xd + xc; 503 if (dis2 >= -halfCarTolerance) 505 if (dis2 >= -halfCarTolerance) 504 { 506 { 505 if (cos2 >= 0) return kInfinity; 507 if (cos2 >= 0) return kInfinity; 506 G4double tmp = -dis2/cos2; 508 G4double tmp = -dis2/cos2; 507 if (tmin2 < tmp) tmin2 = tmp; 509 if (tmin2 < tmp) tmin2 = tmp; 508 } 510 } 509 else if (cos2 > 0) 511 else if (cos2 > 0) 510 { 512 { 511 G4double tmp = -dis2/cos2; 513 G4double tmp = -dis2/cos2; 512 if (tmax2 > tmp) tmax2 = tmp; 514 if (tmax2 > tmp) tmax2 = tmp; 513 } 515 } 514 516 515 G4double tmin3 = tmin2, tmax3 = tmax2; 517 G4double tmin3 = tmin2, tmax3 = tmax2; 516 G4double cos3 = xb - xa; 518 G4double cos3 = xb - xa; 517 G4double dis3 = xd - xc; 519 G4double dis3 = xd - xc; 518 if (dis3 >= -halfCarTolerance) 520 if (dis3 >= -halfCarTolerance) 519 { 521 { 520 if (cos3 >= 0) return kInfinity; 522 if (cos3 >= 0) return kInfinity; 521 G4double tmp = -dis3/cos3; 523 G4double tmp = -dis3/cos3; 522 if (tmin3 < tmp) tmin3 = tmp; 524 if (tmin3 < tmp) tmin3 = tmp; 523 } 525 } 524 else if (cos3 > 0) 526 else if (cos3 > 0) 525 { 527 { 526 G4double tmp = -dis3/cos3; 528 G4double tmp = -dis3/cos3; 527 if (tmax3 > tmp) tmax3 = tmp; 529 if (tmax3 > tmp) tmax3 = tmp; 528 } 530 } 529 531 530 // Find distance 532 // Find distance 531 // 533 // 532 G4double tmin = tmin3, tmax = tmax3; 534 G4double tmin = tmin3, tmax = tmax3; 533 if (tmax <= tmin + halfCarTolerance) return 535 if (tmax <= tmin + halfCarTolerance) return kInfinity; // touch or no hit 534 return (tmin < halfCarTolerance ) ? 0. : tmi 536 return (tmin < halfCarTolerance ) ? 0. : tmin; 535 } 537 } 536 538 537 ////////////////////////////////////////////// 539 ////////////////////////////////////////////////////////////////////////// 538 // 540 // 539 // Calculate exact shortest distance to any bo 541 // Calculate exact shortest distance to any boundary from outside 540 // This is the best fast estimation of the sho 542 // This is the best fast estimation of the shortest distance to trap 541 // - returns 0 if point is inside 543 // - returns 0 if point is inside 542 544 543 G4double G4Trd::DistanceToIn( const G4ThreeVec 545 G4double G4Trd::DistanceToIn( const G4ThreeVector& p ) const 544 { 546 { 545 G4double dx = fPlanes[3].a*std::abs(p.x())+f 547 G4double dx = fPlanes[3].a*std::abs(p.x())+fPlanes[3].c*p.z()+fPlanes[3].d; 546 G4double dy = fPlanes[1].b*std::abs(p.y())+f 548 G4double dy = fPlanes[1].b*std::abs(p.y())+fPlanes[1].c*p.z()+fPlanes[1].d; 547 G4double dxy = std::max(dx,dy); 549 G4double dxy = std::max(dx,dy); 548 550 549 G4double dz = std::abs(p.z())-fDz; 551 G4double dz = std::abs(p.z())-fDz; 550 G4double dist = std::max(dz,dxy); 552 G4double dist = std::max(dz,dxy); 551 553 552 return (dist > 0) ? dist : 0.; 554 return (dist > 0) ? dist : 0.; 553 } 555 } 554 556 555 ////////////////////////////////////////////// 557 ////////////////////////////////////////////////////////////////////////// 556 // 558 // 557 // Calculate distance to surface of shape from 559 // Calculate distance to surface of shape from inside and 558 // find normal at exit point, if required 560 // find normal at exit point, if required 559 // - when leaving the surface, return 0 561 // - when leaving the surface, return 0 560 562 561 G4double G4Trd::DistanceToOut(const G4ThreeVec 563 G4double G4Trd::DistanceToOut(const G4ThreeVector& p, const G4ThreeVector& v, 562 const G4bool cal 564 const G4bool calcNorm, 563 G4bool* va 565 G4bool* validNorm, G4ThreeVector* n) const 564 { 566 { 565 // Z intersections 567 // Z intersections 566 // 568 // 567 if ((std::abs(p.z()) - fDz) >= -halfCarToler 569 if ((std::abs(p.z()) - fDz) >= -halfCarTolerance && p.z()*v.z() > 0) 568 { 570 { 569 if (calcNorm) 571 if (calcNorm) 570 { 572 { 571 *validNorm = true; 573 *validNorm = true; 572 n->set(0, 0, (p.z() < 0) ? -1 : 1); 574 n->set(0, 0, (p.z() < 0) ? -1 : 1); 573 } 575 } 574 return 0; 576 return 0; 575 } 577 } 576 G4double vz = v.z(); 578 G4double vz = v.z(); 577 G4double tmax = (vz == 0) ? DBL_MAX : (std:: 579 G4double tmax = (vz == 0) ? DBL_MAX : (std::copysign(fDz,vz) - p.z())/vz; 578 G4int iside = (vz < 0) ? -4 : -2; // little 580 G4int iside = (vz < 0) ? -4 : -2; // little trick: (-4+3)=-1, (-2+3)=+1 579 581 580 // Y intersections 582 // Y intersections 581 // 583 // 582 G4int i = 0; 584 G4int i = 0; 583 for ( ; i<2; ++i) 585 for ( ; i<2; ++i) 584 { 586 { 585 G4double cosa = fPlanes[i].b*v.y() + fPlan 587 G4double cosa = fPlanes[i].b*v.y() + fPlanes[i].c*v.z(); 586 if (cosa > 0) 588 if (cosa > 0) 587 { 589 { 588 G4double dist = fPlanes[i].b*p.y()+fPlan 590 G4double dist = fPlanes[i].b*p.y()+fPlanes[i].c*p.z()+fPlanes[i].d; 589 if (dist >= -halfCarTolerance) 591 if (dist >= -halfCarTolerance) 590 { 592 { 591 if (calcNorm) 593 if (calcNorm) 592 { 594 { 593 *validNorm = true; 595 *validNorm = true; 594 n->set(0, fPlanes[i].b, fPlanes[i].c 596 n->set(0, fPlanes[i].b, fPlanes[i].c); 595 } 597 } 596 return 0; 598 return 0; 597 } 599 } 598 G4double tmp = -dist/cosa; 600 G4double tmp = -dist/cosa; 599 if (tmax > tmp) { tmax = tmp; iside = i; 601 if (tmax > tmp) { tmax = tmp; iside = i; } 600 } 602 } 601 } 603 } 602 604 603 // X intersections 605 // X intersections 604 // 606 // 605 for ( ; i<4; ++i) 607 for ( ; i<4; ++i) 606 { 608 { 607 G4double cosa = fPlanes[i].a*v.x()+fPlanes 609 G4double cosa = fPlanes[i].a*v.x()+fPlanes[i].c*v.z(); 608 if (cosa > 0) 610 if (cosa > 0) 609 { 611 { 610 G4double dist = fPlanes[i].a*p.x()+fPlan 612 G4double dist = fPlanes[i].a*p.x()+fPlanes[i].c*p.z()+fPlanes[i].d; 611 if (dist >= -halfCarTolerance) 613 if (dist >= -halfCarTolerance) 612 { 614 { 613 if (calcNorm) 615 if (calcNorm) 614 { 616 { 615 *validNorm = true; 617 *validNorm = true; 616 n->set(fPlanes[i].a, fPlanes[i].b, 618 n->set(fPlanes[i].a, fPlanes[i].b, fPlanes[i].c); 617 } 619 } 618 return 0; 620 return 0; 619 } 621 } 620 G4double tmp = -dist/cosa; 622 G4double tmp = -dist/cosa; 621 if (tmax > tmp) { tmax = tmp; iside = i; 623 if (tmax > tmp) { tmax = tmp; iside = i; } 622 } 624 } 623 } 625 } 624 626 625 // Set normal, if required, and return dista 627 // Set normal, if required, and return distance 626 // 628 // 627 if (calcNorm) 629 if (calcNorm) 628 { 630 { 629 *validNorm = true; 631 *validNorm = true; 630 if (iside < 0) 632 if (iside < 0) 631 n->set(0, 0, iside + 3); // (-4+3)=-1, ( 633 n->set(0, 0, iside + 3); // (-4+3)=-1, (-2+3)=+1 632 else 634 else 633 n->set(fPlanes[iside].a, fPlanes[iside]. 635 n->set(fPlanes[iside].a, fPlanes[iside].b, fPlanes[iside].c); 634 } 636 } 635 return tmax; 637 return tmax; 636 } 638 } 637 639 638 ////////////////////////////////////////////// 640 ////////////////////////////////////////////////////////////////////////// 639 // 641 // 640 // Calculate exact shortest distance to any bo 642 // Calculate exact shortest distance to any boundary from inside 641 // - returns 0 if point is outside 643 // - returns 0 if point is outside 642 644 643 G4double G4Trd::DistanceToOut( const G4ThreeVe 645 G4double G4Trd::DistanceToOut( const G4ThreeVector& p ) const 644 { 646 { 645 #ifdef G4CSGDEBUG 647 #ifdef G4CSGDEBUG 646 if( Inside(p) == kOutside ) 648 if( Inside(p) == kOutside ) 647 { 649 { 648 std::ostringstream message; 650 std::ostringstream message; 649 G4long oldprc = message.precision(16); << 651 G4int oldprc = message.precision(16); 650 message << "Point p is outside (!?) of sol 652 message << "Point p is outside (!?) of solid: " << GetName() << G4endl; 651 message << "Position:\n"; 653 message << "Position:\n"; 652 message << " p.x() = " << p.x()/mm << " 654 message << " p.x() = " << p.x()/mm << " mm\n"; 653 message << " p.y() = " << p.y()/mm << " 655 message << " p.y() = " << p.y()/mm << " mm\n"; 654 message << " p.z() = " << p.z()/mm << " 656 message << " p.z() = " << p.z()/mm << " mm"; 655 G4cout.precision(oldprc); 657 G4cout.precision(oldprc); 656 G4Exception("G4Trd::DistanceToOut(p)", "Ge 658 G4Exception("G4Trd::DistanceToOut(p)", "GeomSolids1002", 657 JustWarning, message ); 659 JustWarning, message ); 658 DumpInfo(); 660 DumpInfo(); 659 } 661 } 660 #endif 662 #endif 661 G4double dx = fPlanes[3].a*std::abs(p.x())+f 663 G4double dx = fPlanes[3].a*std::abs(p.x())+fPlanes[3].c*p.z()+fPlanes[3].d; 662 G4double dy = fPlanes[1].b*std::abs(p.y())+f 664 G4double dy = fPlanes[1].b*std::abs(p.y())+fPlanes[1].c*p.z()+fPlanes[1].d; 663 G4double dxy = std::max(dx,dy); 665 G4double dxy = std::max(dx,dy); 664 666 665 G4double dz = std::abs(p.z())-fDz; 667 G4double dz = std::abs(p.z())-fDz; 666 G4double dist = std::max(dz,dxy); 668 G4double dist = std::max(dz,dxy); 667 669 668 return (dist < 0) ? -dist : 0.; 670 return (dist < 0) ? -dist : 0.; 669 } 671 } 670 672 671 ////////////////////////////////////////////// 673 ////////////////////////////////////////////////////////////////////////// 672 // 674 // 673 // GetEntityType 675 // GetEntityType 674 676 675 G4GeometryType G4Trd::GetEntityType() const 677 G4GeometryType G4Trd::GetEntityType() const 676 { 678 { 677 return {"G4Trd"}; << 679 return G4String("G4Trd"); 678 } << 679 << 680 ////////////////////////////////////////////// << 681 // << 682 // IsFaceted << 683 << 684 G4bool G4Trd::IsFaceted() const << 685 { << 686 return true; << 687 } 680 } 688 681 689 ////////////////////////////////////////////// 682 ////////////////////////////////////////////////////////////////////////// 690 // 683 // 691 // Make a clone of the object 684 // Make a clone of the object 692 // 685 // 693 G4VSolid* G4Trd::Clone() const 686 G4VSolid* G4Trd::Clone() const 694 { 687 { 695 return new G4Trd(*this); 688 return new G4Trd(*this); 696 } 689 } 697 690 698 ////////////////////////////////////////////// 691 ////////////////////////////////////////////////////////////////////////// 699 // 692 // 700 // Stream object contents to an output stream 693 // Stream object contents to an output stream 701 694 702 std::ostream& G4Trd::StreamInfo( std::ostream& 695 std::ostream& G4Trd::StreamInfo( std::ostream& os ) const 703 { 696 { 704 G4long oldprc = os.precision(16); << 697 G4int oldprc = os.precision(16); 705 os << "------------------------------------- 698 os << "-----------------------------------------------------------\n" 706 << " *** Dump for solid - " << GetName 699 << " *** Dump for solid - " << GetName() << " ***\n" 707 << " ================================= 700 << " ===================================================\n" 708 << " Solid type: G4Trd\n" 701 << " Solid type: G4Trd\n" 709 << " Parameters: \n" 702 << " Parameters: \n" 710 << " half length X, surface -dZ: " << 703 << " half length X, surface -dZ: " << fDx1/mm << " mm \n" 711 << " half length X, surface +dZ: " << 704 << " half length X, surface +dZ: " << fDx2/mm << " mm \n" 712 << " half length Y, surface -dZ: " << 705 << " half length Y, surface -dZ: " << fDy1/mm << " mm \n" 713 << " half length Y, surface +dZ: " << 706 << " half length Y, surface +dZ: " << fDy2/mm << " mm \n" 714 << " half length Z : " << 707 << " half length Z : " << fDz/mm << " mm \n" 715 << "------------------------------------- 708 << "-----------------------------------------------------------\n"; 716 os.precision(oldprc); 709 os.precision(oldprc); 717 710 718 return os; 711 return os; 719 } 712 } 720 713 721 ////////////////////////////////////////////// 714 ////////////////////////////////////////////////////////////////////////// 722 // 715 // 723 // Return a point randomly and uniformly selec 716 // Return a point randomly and uniformly selected on the solid surface 724 717 725 G4ThreeVector G4Trd::GetPointOnSurface() const 718 G4ThreeVector G4Trd::GetPointOnSurface() const 726 { 719 { 727 // Set areas 720 // Set areas 728 // 721 // 729 G4double sxz = (fDx1 + fDx2)*fHx; 722 G4double sxz = (fDx1 + fDx2)*fHx; 730 G4double syz = (fDy1 + fDy2)*fHy; 723 G4double syz = (fDy1 + fDy2)*fHy; 731 G4double ssurf[6] = { 4.*fDx1*fDy1, sxz, sxz 724 G4double ssurf[6] = { 4.*fDx1*fDy1, sxz, sxz, syz, syz, 4.*fDx2*fDy2 }; 732 ssurf[1] += ssurf[0]; 725 ssurf[1] += ssurf[0]; 733 ssurf[2] += ssurf[1]; 726 ssurf[2] += ssurf[1]; 734 ssurf[3] += ssurf[2]; 727 ssurf[3] += ssurf[2]; 735 ssurf[4] += ssurf[3]; 728 ssurf[4] += ssurf[3]; 736 ssurf[5] += ssurf[4]; 729 ssurf[5] += ssurf[4]; 737 730 738 // Select face 731 // Select face 739 // 732 // 740 G4double select = ssurf[5]*G4QuickRand(); 733 G4double select = ssurf[5]*G4QuickRand(); 741 G4int k = 5; 734 G4int k = 5; 742 k -= (G4int)(select <= ssurf[4]); << 735 k -= (select <= ssurf[4]); 743 k -= (G4int)(select <= ssurf[3]); << 736 k -= (select <= ssurf[3]); 744 k -= (G4int)(select <= ssurf[2]); << 737 k -= (select <= ssurf[2]); 745 k -= (G4int)(select <= ssurf[1]); << 738 k -= (select <= ssurf[1]); 746 k -= (G4int)(select <= ssurf[0]); << 739 k -= (select <= ssurf[0]); 747 740 748 // Generate point on selected surface 741 // Generate point on selected surface 749 // 742 // 750 G4double u = G4QuickRand(); 743 G4double u = G4QuickRand(); 751 G4double v = G4QuickRand(); 744 G4double v = G4QuickRand(); 752 switch(k) 745 switch(k) 753 { 746 { 754 case 0: // base at -Z 747 case 0: // base at -Z 755 { 748 { 756 return { (2.*u - 1.)*fDx1, (2.*v - 1.)*f << 749 return G4ThreeVector((2.*u - 1.)*fDx1, (2.*v - 1.)*fDy1, -fDz); 757 } 750 } 758 case 1: // X face at -Y 751 case 1: // X face at -Y 759 { 752 { 760 if (u + v > 1.) { u = 1. - u; v = 1. - v 753 if (u + v > 1.) { u = 1. - u; v = 1. - v; } 761 G4ThreeVector p0(-fDx1,-fDy1,-fDz); 754 G4ThreeVector p0(-fDx1,-fDy1,-fDz); 762 G4ThreeVector p1( fDx2,-fDy2, fDz); 755 G4ThreeVector p1( fDx2,-fDy2, fDz); 763 return (select <= ssurf[0] + fDx1*fHx) ? 756 return (select <= ssurf[0] + fDx1*fHx) ? 764 (1. - u - v)*p0 + u*p1 + v*G4ThreeVect 757 (1. - u - v)*p0 + u*p1 + v*G4ThreeVector( fDx1,-fDy1,-fDz) : 765 (1. - u - v)*p0 + u*p1 + v*G4ThreeVect 758 (1. - u - v)*p0 + u*p1 + v*G4ThreeVector(-fDx2,-fDy2, fDz); 766 } 759 } 767 case 2: // X face at +Y 760 case 2: // X face at +Y 768 { 761 { 769 if (u + v > 1.) { u = 1. - u; v = 1. - v 762 if (u + v > 1.) { u = 1. - u; v = 1. - v; } 770 G4ThreeVector p0( fDx1, fDy1,-fDz); 763 G4ThreeVector p0( fDx1, fDy1,-fDz); 771 G4ThreeVector p1(-fDx2, fDy2, fDz); 764 G4ThreeVector p1(-fDx2, fDy2, fDz); 772 return (select <= ssurf[1] + fDx1*fHx) ? 765 return (select <= ssurf[1] + fDx1*fHx) ? 773 (1. - u - v)*p0 + u*p1 + v*G4ThreeVect 766 (1. - u - v)*p0 + u*p1 + v*G4ThreeVector(-fDx1, fDy1,-fDz) : 774 (1. - u - v)*p0 + u*p1 + v*G4ThreeVect 767 (1. - u - v)*p0 + u*p1 + v*G4ThreeVector( fDx2, fDy2, fDz); 775 } 768 } 776 case 3: // Y face at -X 769 case 3: // Y face at -X 777 { 770 { 778 if (u + v > 1.) { u = 1. - u; v = 1. - v 771 if (u + v > 1.) { u = 1. - u; v = 1. - v; } 779 G4ThreeVector p0(-fDx1, fDy1,-fDz); 772 G4ThreeVector p0(-fDx1, fDy1,-fDz); 780 G4ThreeVector p1(-fDx2,-fDy2, fDz); 773 G4ThreeVector p1(-fDx2,-fDy2, fDz); 781 return (select <= ssurf[2] + fDy1*fHy) ? 774 return (select <= ssurf[2] + fDy1*fHy) ? 782 (1. - u - v)*p0 + u*p1 + v*G4ThreeVect 775 (1. - u - v)*p0 + u*p1 + v*G4ThreeVector(-fDx1,-fDy1,-fDz) : 783 (1. - u - v)*p0 + u*p1 + v*G4ThreeVect 776 (1. - u - v)*p0 + u*p1 + v*G4ThreeVector(-fDx2, fDy2, fDz); 784 } 777 } 785 case 4: // Y face at +X 778 case 4: // Y face at +X 786 { 779 { 787 if (u + v > 1.) { u = 1. - u; v = 1. - v 780 if (u + v > 1.) { u = 1. - u; v = 1. - v; } 788 G4ThreeVector p0( fDx1,-fDy1,-fDz); 781 G4ThreeVector p0( fDx1,-fDy1,-fDz); 789 G4ThreeVector p1( fDx2, fDy2, fDz); 782 G4ThreeVector p1( fDx2, fDy2, fDz); 790 return (select <= ssurf[3] + fDy1*fHy) ? 783 return (select <= ssurf[3] + fDy1*fHy) ? 791 (1. - u - v)*p0 + u*p1 + v*G4ThreeVect 784 (1. - u - v)*p0 + u*p1 + v*G4ThreeVector( fDx1, fDy1,-fDz) : 792 (1. - u - v)*p0 + u*p1 + v*G4ThreeVect 785 (1. - u - v)*p0 + u*p1 + v*G4ThreeVector( fDx2,-fDy2, fDz); 793 } 786 } 794 case 5: // base at +Z 787 case 5: // base at +Z 795 { 788 { 796 return { (2.*u - 1.)*fDx2, (2.*v - 1.)*f << 789 return G4ThreeVector((2.*u - 1.)*fDx2, (2.*v - 1.)*fDy2, fDz); 797 } 790 } 798 } 791 } 799 return {0., 0., 0.}; << 792 return G4ThreeVector(0., 0., 0.); 800 } 793 } 801 794 802 ////////////////////////////////////////////// 795 ////////////////////////////////////////////////////////////////////////// 803 // 796 // 804 // Methods for visualisation 797 // Methods for visualisation 805 798 806 void G4Trd::DescribeYourselfTo ( G4VGraphicsSc 799 void G4Trd::DescribeYourselfTo ( G4VGraphicsScene& scene ) const 807 { 800 { 808 scene.AddSolid (*this); 801 scene.AddSolid (*this); 809 } 802 } 810 803 811 G4Polyhedron* G4Trd::CreatePolyhedron () const 804 G4Polyhedron* G4Trd::CreatePolyhedron () const 812 { 805 { 813 return new G4PolyhedronTrd2 (fDx1, fDx2, fDy 806 return new G4PolyhedronTrd2 (fDx1, fDx2, fDy1, fDy2, fDz); 814 } 807 } 815 808 816 #endif 809 #endif 817 810