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 G4Para class 26 // Implementation for G4Para class 27 // 27 // 28 // 21.03.95 P.Kent: Modified for `tolerant' ge 28 // 21.03.95 P.Kent: Modified for `tolerant' geom 29 // 31.10.96 V.Grichine: Modifications accordin 29 // 31.10.96 V.Grichine: Modifications according G4Box/Tubs before to commit 30 // 28.04.05 V.Grichine: new SurfaceNormal acco 30 // 28.04.05 V.Grichine: new SurfaceNormal according to J. Apostolakis proposal 31 // 29.05.17 E.Tcherniaev: complete revision, s 31 // 29.05.17 E.Tcherniaev: complete revision, speed-up 32 ////////////////////////////////////////////// 32 //////////////////////////////////////////////////////////////////////////// 33 33 34 #include "G4Para.hh" 34 #include "G4Para.hh" 35 35 36 #if !defined(G4GEOM_USE_UPARA) 36 #if !defined(G4GEOM_USE_UPARA) 37 37 38 #include "G4VoxelLimits.hh" 38 #include "G4VoxelLimits.hh" 39 #include "G4AffineTransform.hh" 39 #include "G4AffineTransform.hh" 40 #include "G4BoundingEnvelope.hh" 40 #include "G4BoundingEnvelope.hh" 41 #include "Randomize.hh" 41 #include "Randomize.hh" 42 42 43 #include "G4VPVParameterisation.hh" 43 #include "G4VPVParameterisation.hh" 44 44 45 #include "G4VGraphicsScene.hh" 45 #include "G4VGraphicsScene.hh" 46 46 47 using namespace CLHEP; 47 using namespace CLHEP; 48 48 49 ////////////////////////////////////////////// 49 ////////////////////////////////////////////////////////////////////////// 50 // 50 // 51 // Constructor - set & check half widths 51 // Constructor - set & check half widths 52 52 53 G4Para::G4Para(const G4String& pName, 53 G4Para::G4Para(const G4String& pName, 54 G4double pDx, G4double pD 54 G4double pDx, G4double pDy, G4double pDz, 55 G4double pAlpha, G4double 55 G4double pAlpha, G4double pTheta, G4double pPhi) 56 : G4CSGSolid(pName), halfCarTolerance(0.5*kC 56 : G4CSGSolid(pName), halfCarTolerance(0.5*kCarTolerance) 57 { 57 { 58 SetAllParameters(pDx, pDy, pDz, pAlpha, pThe 58 SetAllParameters(pDx, pDy, pDz, pAlpha, pTheta, pPhi); 59 fRebuildPolyhedron = false; // default valu 59 fRebuildPolyhedron = false; // default value for G4CSGSolid 60 } 60 } 61 61 62 ////////////////////////////////////////////// 62 ////////////////////////////////////////////////////////////////////////// 63 // 63 // 64 // Constructor - design of trapezoid based on 64 // Constructor - design of trapezoid based on 8 vertices 65 65 66 G4Para::G4Para( const G4String& pName, 66 G4Para::G4Para( const G4String& pName, 67 const G4ThreeVector pt[8] ) 67 const G4ThreeVector pt[8] ) 68 : G4CSGSolid(pName), halfCarTolerance(0.5*kC 68 : G4CSGSolid(pName), halfCarTolerance(0.5*kCarTolerance) 69 { 69 { 70 // Find dimensions and trigonometric values 70 // Find dimensions and trigonometric values 71 // 71 // 72 fDx = (pt[3].x() - pt[2].x())*0.5; 72 fDx = (pt[3].x() - pt[2].x())*0.5; 73 fDy = (pt[2].y() - pt[1].y())*0.5; 73 fDy = (pt[2].y() - pt[1].y())*0.5; 74 fDz = pt[7].z(); 74 fDz = pt[7].z(); 75 CheckParameters(); // check dimensions 75 CheckParameters(); // check dimensions 76 76 77 fTalpha = (pt[2].x() + pt[3].x() - pt[1].x() 77 fTalpha = (pt[2].x() + pt[3].x() - pt[1].x() - pt[0].x())*0.25/fDy; 78 fTthetaCphi = (pt[4].x() + fDy*fTalpha + fDx 78 fTthetaCphi = (pt[4].x() + fDy*fTalpha + fDx)/fDz; 79 fTthetaSphi = (pt[4].y() + fDy)/fDz; 79 fTthetaSphi = (pt[4].y() + fDy)/fDz; 80 MakePlanes(); 80 MakePlanes(); 81 81 82 // Recompute vertices 82 // Recompute vertices 83 // 83 // 84 G4ThreeVector v[8]; 84 G4ThreeVector v[8]; 85 G4double DyTalpha = fDy*fTalpha; 85 G4double DyTalpha = fDy*fTalpha; 86 G4double DzTthetaSphi = fDz*fTthetaSphi; 86 G4double DzTthetaSphi = fDz*fTthetaSphi; 87 G4double DzTthetaCphi = fDz*fTthetaCphi; 87 G4double DzTthetaCphi = fDz*fTthetaCphi; 88 v[0].set(-DzTthetaCphi-DyTalpha-fDx, -DzTthe 88 v[0].set(-DzTthetaCphi-DyTalpha-fDx, -DzTthetaSphi-fDy, -fDz); 89 v[1].set(-DzTthetaCphi-DyTalpha+fDx, -DzTthe 89 v[1].set(-DzTthetaCphi-DyTalpha+fDx, -DzTthetaSphi-fDy, -fDz); 90 v[2].set(-DzTthetaCphi+DyTalpha-fDx, -DzTthe 90 v[2].set(-DzTthetaCphi+DyTalpha-fDx, -DzTthetaSphi+fDy, -fDz); 91 v[3].set(-DzTthetaCphi+DyTalpha+fDx, -DzTthe 91 v[3].set(-DzTthetaCphi+DyTalpha+fDx, -DzTthetaSphi+fDy, -fDz); 92 v[4].set( DzTthetaCphi-DyTalpha-fDx, DzTthe 92 v[4].set( DzTthetaCphi-DyTalpha-fDx, DzTthetaSphi-fDy, fDz); 93 v[5].set( DzTthetaCphi-DyTalpha+fDx, DzTthe 93 v[5].set( DzTthetaCphi-DyTalpha+fDx, DzTthetaSphi-fDy, fDz); 94 v[6].set( DzTthetaCphi+DyTalpha-fDx, DzTthe 94 v[6].set( DzTthetaCphi+DyTalpha-fDx, DzTthetaSphi+fDy, fDz); 95 v[7].set( DzTthetaCphi+DyTalpha+fDx, DzTthe 95 v[7].set( DzTthetaCphi+DyTalpha+fDx, DzTthetaSphi+fDy, fDz); 96 96 97 // Compare with original vertices 97 // Compare with original vertices 98 // 98 // 99 for (G4int i=0; i<8; ++i) 99 for (G4int i=0; i<8; ++i) 100 { 100 { 101 G4double delx = std::abs(pt[i].x() - v[i]. 101 G4double delx = std::abs(pt[i].x() - v[i].x()); 102 G4double dely = std::abs(pt[i].y() - v[i]. 102 G4double dely = std::abs(pt[i].y() - v[i].y()); 103 G4double delz = std::abs(pt[i].z() - v[i]. 103 G4double delz = std::abs(pt[i].z() - v[i].z()); 104 G4double discrepancy = std::max(std::max(d 104 G4double discrepancy = std::max(std::max(delx,dely),delz); 105 if (discrepancy > 0.1*kCarTolerance) 105 if (discrepancy > 0.1*kCarTolerance) 106 { 106 { 107 std::ostringstream message; 107 std::ostringstream message; 108 G4long oldprc = message.precision(16); << 108 G4int oldprc = message.precision(16); 109 message << "Invalid vertice coordinates 109 message << "Invalid vertice coordinates for Solid: " << GetName() 110 << "\nVertix #" << i << ", discr 110 << "\nVertix #" << i << ", discrepancy = " << discrepancy 111 << "\n original : " << pt[i] 111 << "\n original : " << pt[i] 112 << "\n recomputed : " << v[i]; 112 << "\n recomputed : " << v[i]; 113 G4cout.precision(oldprc); 113 G4cout.precision(oldprc); 114 G4Exception("G4Para::G4Para()", "GeomSol 114 G4Exception("G4Para::G4Para()", "GeomSolids0002", 115 FatalException, message); 115 FatalException, message); 116 116 117 } 117 } 118 } 118 } 119 } 119 } 120 120 121 ////////////////////////////////////////////// 121 ////////////////////////////////////////////////////////////////////////// 122 // 122 // 123 // Fake default constructor - sets only member 123 // Fake default constructor - sets only member data and allocates memory 124 // for usage restri 124 // for usage restricted to object persistency 125 125 126 G4Para::G4Para( __void__& a ) 126 G4Para::G4Para( __void__& a ) 127 : G4CSGSolid(a), halfCarTolerance(0.5*kCarTo 127 : G4CSGSolid(a), halfCarTolerance(0.5*kCarTolerance) 128 { 128 { 129 SetAllParameters(1., 1., 1., 0., 0., 0.); 129 SetAllParameters(1., 1., 1., 0., 0., 0.); 130 fRebuildPolyhedron = false; // default value 130 fRebuildPolyhedron = false; // default value for G4CSGSolid 131 } 131 } 132 132 133 ////////////////////////////////////////////// 133 ////////////////////////////////////////////////////////////////////////// 134 // 134 // 135 // Destructor 135 // Destructor 136 136 137 G4Para::~G4Para() = default; << 137 G4Para::~G4Para() >> 138 { >> 139 } 138 140 139 ////////////////////////////////////////////// 141 ////////////////////////////////////////////////////////////////////////// 140 // 142 // 141 // Copy constructor 143 // Copy constructor 142 144 143 G4Para::G4Para(const G4Para& rhs) 145 G4Para::G4Para(const G4Para& rhs) 144 : G4CSGSolid(rhs), halfCarTolerance(rhs.half 146 : G4CSGSolid(rhs), halfCarTolerance(rhs.halfCarTolerance), 145 fDx(rhs.fDx), fDy(rhs.fDy), fDz(rhs.fDz), 147 fDx(rhs.fDx), fDy(rhs.fDy), fDz(rhs.fDz), fTalpha(rhs.fTalpha), 146 fTthetaCphi(rhs.fTthetaCphi),fTthetaSphi(r 148 fTthetaCphi(rhs.fTthetaCphi),fTthetaSphi(rhs.fTthetaSphi) 147 { 149 { 148 for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs 150 for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs.fPlanes[i]; } 149 } 151 } 150 152 151 ////////////////////////////////////////////// 153 ////////////////////////////////////////////////////////////////////////// 152 // 154 // 153 // Assignment operator 155 // Assignment operator 154 156 155 G4Para& G4Para::operator = (const G4Para& rhs) 157 G4Para& G4Para::operator = (const G4Para& rhs) 156 { 158 { 157 // Check assignment to self 159 // Check assignment to self 158 // 160 // 159 if (this == &rhs) { return *this; } 161 if (this == &rhs) { return *this; } 160 162 161 // Copy base class data 163 // Copy base class data 162 // 164 // 163 G4CSGSolid::operator=(rhs); 165 G4CSGSolid::operator=(rhs); 164 166 165 // Copy data 167 // Copy data 166 // 168 // 167 halfCarTolerance = rhs.halfCarTolerance; 169 halfCarTolerance = rhs.halfCarTolerance; 168 fDx = rhs.fDx; 170 fDx = rhs.fDx; 169 fDy = rhs.fDy; 171 fDy = rhs.fDy; 170 fDz = rhs.fDz; 172 fDz = rhs.fDz; 171 fTalpha = rhs.fTalpha; 173 fTalpha = rhs.fTalpha; 172 fTthetaCphi = rhs.fTthetaCphi; 174 fTthetaCphi = rhs.fTthetaCphi; 173 fTthetaSphi = rhs.fTthetaSphi; 175 fTthetaSphi = rhs.fTthetaSphi; 174 for (G4int i=0; i<4; ++i) { fPlanes[i] = rh 176 for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs.fPlanes[i]; } 175 177 176 return *this; 178 return *this; 177 } 179 } 178 180 179 ////////////////////////////////////////////// 181 ////////////////////////////////////////////////////////////////////////// 180 // 182 // 181 // Set all parameters, as for constructor - se 183 // Set all parameters, as for constructor - set and check half-widths 182 184 183 void G4Para::SetAllParameters(G4double pDx, G4 185 void G4Para::SetAllParameters(G4double pDx, G4double pDy, G4double pDz, 184 G4double pAlpha, 186 G4double pAlpha, G4double pTheta, G4double pPhi) 185 { 187 { 186 // Reset data of the base class 188 // Reset data of the base class 187 fCubicVolume = 0; 189 fCubicVolume = 0; 188 fSurfaceArea = 0; 190 fSurfaceArea = 0; 189 fRebuildPolyhedron = true; 191 fRebuildPolyhedron = true; 190 192 191 // Set parameters 193 // Set parameters 192 fDx = pDx; 194 fDx = pDx; 193 fDy = pDy; 195 fDy = pDy; 194 fDz = pDz; 196 fDz = pDz; 195 fTalpha = std::tan(pAlpha); 197 fTalpha = std::tan(pAlpha); 196 fTthetaCphi = std::tan(pTheta)*std::cos(pPhi 198 fTthetaCphi = std::tan(pTheta)*std::cos(pPhi); 197 fTthetaSphi = std::tan(pTheta)*std::sin(pPhi 199 fTthetaSphi = std::tan(pTheta)*std::sin(pPhi); 198 200 199 CheckParameters(); 201 CheckParameters(); 200 MakePlanes(); 202 MakePlanes(); 201 } 203 } 202 204 203 ////////////////////////////////////////////// 205 ////////////////////////////////////////////////////////////////////////// 204 // 206 // 205 // Check dimensions 207 // Check dimensions 206 208 207 void G4Para::CheckParameters() 209 void G4Para::CheckParameters() 208 { 210 { 209 if (fDx < 2*kCarTolerance || 211 if (fDx < 2*kCarTolerance || 210 fDy < 2*kCarTolerance || 212 fDy < 2*kCarTolerance || 211 fDz < 2*kCarTolerance) 213 fDz < 2*kCarTolerance) 212 { 214 { 213 std::ostringstream message; 215 std::ostringstream message; 214 message << "Invalid (too small or negative 216 message << "Invalid (too small or negative) dimensions for Solid: " 215 << GetName() 217 << GetName() 216 << "\n X - " << fDx 218 << "\n X - " << fDx 217 << "\n Y - " << fDy 219 << "\n Y - " << fDy 218 << "\n Z - " << fDz; 220 << "\n Z - " << fDz; 219 G4Exception("G4Para::CheckParameters()", " 221 G4Exception("G4Para::CheckParameters()", "GeomSolids0002", 220 FatalException, message); 222 FatalException, message); 221 } 223 } 222 } 224 } 223 225 224 ////////////////////////////////////////////// 226 ////////////////////////////////////////////////////////////////////////// 225 // 227 // 226 // Set side planes 228 // Set side planes 227 229 228 void G4Para::MakePlanes() 230 void G4Para::MakePlanes() 229 { 231 { 230 G4ThreeVector vx(1, 0, 0); 232 G4ThreeVector vx(1, 0, 0); 231 G4ThreeVector vy(fTalpha, 1, 0); 233 G4ThreeVector vy(fTalpha, 1, 0); 232 G4ThreeVector vz(fTthetaCphi, fTthetaSphi, 1 234 G4ThreeVector vz(fTthetaCphi, fTthetaSphi, 1); 233 235 234 // Set -Y & +Y planes 236 // Set -Y & +Y planes 235 // 237 // 236 G4ThreeVector ynorm = (vx.cross(vz)).unit(); 238 G4ThreeVector ynorm = (vx.cross(vz)).unit(); 237 239 238 fPlanes[0].a = 0.; 240 fPlanes[0].a = 0.; 239 fPlanes[0].b = ynorm.y(); 241 fPlanes[0].b = ynorm.y(); 240 fPlanes[0].c = ynorm.z(); 242 fPlanes[0].c = ynorm.z(); 241 fPlanes[0].d = fPlanes[0].b*fDy; // point (0 243 fPlanes[0].d = fPlanes[0].b*fDy; // point (0,fDy,0) is on plane 242 244 243 fPlanes[1].a = 0.; 245 fPlanes[1].a = 0.; 244 fPlanes[1].b = -fPlanes[0].b; 246 fPlanes[1].b = -fPlanes[0].b; 245 fPlanes[1].c = -fPlanes[0].c; 247 fPlanes[1].c = -fPlanes[0].c; 246 fPlanes[1].d = fPlanes[0].d; 248 fPlanes[1].d = fPlanes[0].d; 247 249 248 // Set -X & +X planes 250 // Set -X & +X planes 249 // 251 // 250 G4ThreeVector xnorm = (vz.cross(vy)).unit(); 252 G4ThreeVector xnorm = (vz.cross(vy)).unit(); 251 253 252 fPlanes[2].a = xnorm.x(); 254 fPlanes[2].a = xnorm.x(); 253 fPlanes[2].b = xnorm.y(); 255 fPlanes[2].b = xnorm.y(); 254 fPlanes[2].c = xnorm.z(); 256 fPlanes[2].c = xnorm.z(); 255 fPlanes[2].d = fPlanes[2].a*fDx; // point (f 257 fPlanes[2].d = fPlanes[2].a*fDx; // point (fDx,0,0) is on plane 256 258 257 fPlanes[3].a = -fPlanes[2].a; 259 fPlanes[3].a = -fPlanes[2].a; 258 fPlanes[3].b = -fPlanes[2].b; 260 fPlanes[3].b = -fPlanes[2].b; 259 fPlanes[3].c = -fPlanes[2].c; 261 fPlanes[3].c = -fPlanes[2].c; 260 fPlanes[3].d = fPlanes[2].d; 262 fPlanes[3].d = fPlanes[2].d; 261 } 263 } 262 264 263 ////////////////////////////////////////////// 265 ////////////////////////////////////////////////////////////////////////// 264 // 266 // 265 // Get volume 267 // Get volume 266 268 267 G4double G4Para::GetCubicVolume() 269 G4double G4Para::GetCubicVolume() 268 { 270 { 269 // It is like G4Box, since para transformati 271 // It is like G4Box, since para transformations keep the volume to be const 270 if (fCubicVolume == 0) 272 if (fCubicVolume == 0) 271 { 273 { 272 fCubicVolume = 8*fDx*fDy*fDz; 274 fCubicVolume = 8*fDx*fDy*fDz; 273 } 275 } 274 return fCubicVolume; 276 return fCubicVolume; 275 } 277 } 276 278 277 ////////////////////////////////////////////// 279 ////////////////////////////////////////////////////////////////////////// 278 // 280 // 279 // Get surface area 281 // Get surface area 280 282 281 G4double G4Para::GetSurfaceArea() 283 G4double G4Para::GetSurfaceArea() 282 { 284 { 283 if(fSurfaceArea == 0) 285 if(fSurfaceArea == 0) 284 { 286 { 285 G4ThreeVector vx(fDx, 0, 0); 287 G4ThreeVector vx(fDx, 0, 0); 286 G4ThreeVector vy(fDy*fTalpha, fDy, 0); 288 G4ThreeVector vy(fDy*fTalpha, fDy, 0); 287 G4ThreeVector vz(fDz*fTthetaCphi, fDz*fTth 289 G4ThreeVector vz(fDz*fTthetaCphi, fDz*fTthetaSphi, fDz); 288 290 289 G4double sxy = fDx*fDy; // (vx.cross(vy)). 291 G4double sxy = fDx*fDy; // (vx.cross(vy)).mag(); 290 G4double sxz = (vx.cross(vz)).mag(); 292 G4double sxz = (vx.cross(vz)).mag(); 291 G4double syz = (vy.cross(vz)).mag(); 293 G4double syz = (vy.cross(vz)).mag(); 292 294 293 fSurfaceArea = 8*(sxy+sxz+syz); 295 fSurfaceArea = 8*(sxy+sxz+syz); 294 } 296 } 295 return fSurfaceArea; 297 return fSurfaceArea; 296 } 298 } 297 299 298 ////////////////////////////////////////////// 300 ////////////////////////////////////////////////////////////////////////// 299 // 301 // 300 // Dispatch to parameterisation for replicatio 302 // Dispatch to parameterisation for replication mechanism dimension 301 // computation & modification 303 // computation & modification 302 304 303 void G4Para::ComputeDimensions( G4VPVPara 305 void G4Para::ComputeDimensions( G4VPVParameterisation* p, 304 const G4int n, 306 const G4int n, 305 const G4VPhysi 307 const G4VPhysicalVolume* pRep ) 306 { 308 { 307 p->ComputeDimensions(*this,n,pRep); 309 p->ComputeDimensions(*this,n,pRep); 308 } 310 } 309 311 310 ////////////////////////////////////////////// 312 ////////////////////////////////////////////////////////////////////////// 311 // 313 // 312 // Get bounding box 314 // Get bounding box 313 315 314 void G4Para::BoundingLimits(G4ThreeVector& pMi 316 void G4Para::BoundingLimits(G4ThreeVector& pMin, G4ThreeVector& pMax) const 315 { 317 { 316 G4double dz = GetZHalfLength(); 318 G4double dz = GetZHalfLength(); 317 G4double dx = GetXHalfLength(); 319 G4double dx = GetXHalfLength(); 318 G4double dy = GetYHalfLength(); 320 G4double dy = GetYHalfLength(); 319 321 320 G4double x0 = dz*fTthetaCphi; 322 G4double x0 = dz*fTthetaCphi; 321 G4double x1 = dy*GetTanAlpha(); 323 G4double x1 = dy*GetTanAlpha(); 322 G4double xmin = 324 G4double xmin = 323 std::min( 325 std::min( 324 std::min( 326 std::min( 325 std::min(-x0-x1-dx,-x0+x1-dx),x0-x1-dx),x0 327 std::min(-x0-x1-dx,-x0+x1-dx),x0-x1-dx),x0+x1-dx); 326 G4double xmax = 328 G4double xmax = 327 std::max( 329 std::max( 328 std::max( 330 std::max( 329 std::max(-x0-x1+dx,-x0+x1+dx),x0-x1+dx),x0 331 std::max(-x0-x1+dx,-x0+x1+dx),x0-x1+dx),x0+x1+dx); 330 332 331 G4double y0 = dz*fTthetaSphi; 333 G4double y0 = dz*fTthetaSphi; 332 G4double ymin = std::min(-y0-dy,y0-dy); 334 G4double ymin = std::min(-y0-dy,y0-dy); 333 G4double ymax = std::max(-y0+dy,y0+dy); 335 G4double ymax = std::max(-y0+dy,y0+dy); 334 336 335 pMin.set(xmin,ymin,-dz); 337 pMin.set(xmin,ymin,-dz); 336 pMax.set(xmax,ymax, dz); 338 pMax.set(xmax,ymax, dz); 337 339 338 // Check correctness of the bounding box 340 // Check correctness of the bounding box 339 // 341 // 340 if (pMin.x() >= pMax.x() || pMin.y() >= pMax 342 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z()) 341 { 343 { 342 std::ostringstream message; 344 std::ostringstream message; 343 message << "Bad bounding box (min >= max) 345 message << "Bad bounding box (min >= max) for solid: " 344 << GetName() << " !" 346 << GetName() << " !" 345 << "\npMin = " << pMin 347 << "\npMin = " << pMin 346 << "\npMax = " << pMax; 348 << "\npMax = " << pMax; 347 G4Exception("G4Para::BoundingLimits()", "G 349 G4Exception("G4Para::BoundingLimits()", "GeomMgt0001", 348 JustWarning, message); 350 JustWarning, message); 349 DumpInfo(); 351 DumpInfo(); 350 } 352 } 351 } 353 } 352 354 353 ////////////////////////////////////////////// 355 ////////////////////////////////////////////////////////////////////////// 354 // 356 // 355 // Calculate extent under transform and specif 357 // Calculate extent under transform and specified limit 356 358 357 G4bool G4Para::CalculateExtent( const EAxis pA 359 G4bool G4Para::CalculateExtent( const EAxis pAxis, 358 const G4VoxelL 360 const G4VoxelLimits& pVoxelLimit, 359 const G4Affine 361 const G4AffineTransform& pTransform, 360 G4double& 362 G4double& pMin, G4double& pMax ) const 361 { 363 { 362 G4ThreeVector bmin, bmax; 364 G4ThreeVector bmin, bmax; 363 G4bool exist; 365 G4bool exist; 364 366 365 // Check bounding box (bbox) 367 // Check bounding box (bbox) 366 // 368 // 367 BoundingLimits(bmin,bmax); 369 BoundingLimits(bmin,bmax); 368 G4BoundingEnvelope bbox(bmin,bmax); 370 G4BoundingEnvelope bbox(bmin,bmax); 369 #ifdef G4BBOX_EXTENT 371 #ifdef G4BBOX_EXTENT 370 return bbox.CalculateExtent(pAxis,pVoxelLimi 372 return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax); 371 #endif 373 #endif 372 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox 374 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax)) 373 { 375 { 374 return exist = pMin < pMax; << 376 return exist = (pMin < pMax) ? true : false; 375 } 377 } 376 378 377 // Set bounding envelope (benv) and calculat 379 // Set bounding envelope (benv) and calculate extent 378 // 380 // 379 G4double dz = GetZHalfLength(); 381 G4double dz = GetZHalfLength(); 380 G4double dx = GetXHalfLength(); 382 G4double dx = GetXHalfLength(); 381 G4double dy = GetYHalfLength(); 383 G4double dy = GetYHalfLength(); 382 384 383 G4double x0 = dz*fTthetaCphi; 385 G4double x0 = dz*fTthetaCphi; 384 G4double x1 = dy*GetTanAlpha(); 386 G4double x1 = dy*GetTanAlpha(); 385 G4double y0 = dz*fTthetaSphi; 387 G4double y0 = dz*fTthetaSphi; 386 388 387 G4ThreeVectorList baseA(4), baseB(4); 389 G4ThreeVectorList baseA(4), baseB(4); 388 baseA[0].set(-x0-x1-dx,-y0-dy,-dz); 390 baseA[0].set(-x0-x1-dx,-y0-dy,-dz); 389 baseA[1].set(-x0-x1+dx,-y0-dy,-dz); 391 baseA[1].set(-x0-x1+dx,-y0-dy,-dz); 390 baseA[2].set(-x0+x1+dx,-y0+dy,-dz); 392 baseA[2].set(-x0+x1+dx,-y0+dy,-dz); 391 baseA[3].set(-x0+x1-dx,-y0+dy,-dz); 393 baseA[3].set(-x0+x1-dx,-y0+dy,-dz); 392 394 393 baseB[0].set(+x0-x1-dx, y0-dy, dz); 395 baseB[0].set(+x0-x1-dx, y0-dy, dz); 394 baseB[1].set(+x0-x1+dx, y0-dy, dz); 396 baseB[1].set(+x0-x1+dx, y0-dy, dz); 395 baseB[2].set(+x0+x1+dx, y0+dy, dz); 397 baseB[2].set(+x0+x1+dx, y0+dy, dz); 396 baseB[3].set(+x0+x1-dx, y0+dy, dz); 398 baseB[3].set(+x0+x1-dx, y0+dy, dz); 397 399 398 std::vector<const G4ThreeVectorList *> polyg 400 std::vector<const G4ThreeVectorList *> polygons(2); 399 polygons[0] = &baseA; 401 polygons[0] = &baseA; 400 polygons[1] = &baseB; 402 polygons[1] = &baseB; 401 403 402 G4BoundingEnvelope benv(bmin,bmax,polygons); 404 G4BoundingEnvelope benv(bmin,bmax,polygons); 403 exist = benv.CalculateExtent(pAxis,pVoxelLim 405 exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax); 404 return exist; 406 return exist; 405 } 407 } 406 408 407 ////////////////////////////////////////////// 409 ////////////////////////////////////////////////////////////////////////// 408 // 410 // 409 // Determine where is point p, inside/on_surfa 411 // Determine where is point p, inside/on_surface/outside 410 // 412 // 411 413 412 EInside G4Para::Inside( const G4ThreeVector& p 414 EInside G4Para::Inside( const G4ThreeVector& p ) const 413 { 415 { 414 G4double xx = fPlanes[2].a*p.x()+fPlanes[2]. 416 G4double xx = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()+fPlanes[2].c*p.z(); 415 G4double dx = std::abs(xx) + fPlanes[2].d; 417 G4double dx = std::abs(xx) + fPlanes[2].d; 416 418 417 G4double yy = fPlanes[0].b*p.y()+fPlanes[0]. 419 G4double yy = fPlanes[0].b*p.y()+fPlanes[0].c*p.z(); 418 G4double dy = std::abs(yy) + fPlanes[0].d; 420 G4double dy = std::abs(yy) + fPlanes[0].d; 419 G4double dxy = std::max(dx,dy); 421 G4double dxy = std::max(dx,dy); 420 422 421 G4double dz = std::abs(p.z())-fDz; 423 G4double dz = std::abs(p.z())-fDz; 422 G4double dist = std::max(dxy,dz); 424 G4double dist = std::max(dxy,dz); 423 425 424 if (dist > halfCarTolerance) return kOutside 426 if (dist > halfCarTolerance) return kOutside; 425 return (dist > -halfCarTolerance) ? kSurface 427 return (dist > -halfCarTolerance) ? kSurface : kInside; 426 } 428 } 427 429 428 ////////////////////////////////////////////// 430 ////////////////////////////////////////////////////////////////////////// 429 // 431 // 430 // Determine side where point is, and return c 432 // Determine side where point is, and return corresponding normal 431 433 432 G4ThreeVector G4Para::SurfaceNormal( const G4T 434 G4ThreeVector G4Para::SurfaceNormal( const G4ThreeVector& p ) const 433 { 435 { 434 G4int nsurf = 0; // number of surfaces where 436 G4int nsurf = 0; // number of surfaces where p is placed 435 437 436 // Check Z faces 438 // Check Z faces 437 // 439 // 438 G4double nz = 0; 440 G4double nz = 0; 439 G4double dz = std::abs(p.z()) - fDz; 441 G4double dz = std::abs(p.z()) - fDz; 440 if (std::abs(dz) <= halfCarTolerance) 442 if (std::abs(dz) <= halfCarTolerance) 441 { 443 { 442 nz = (p.z() < 0) ? -1 : 1; 444 nz = (p.z() < 0) ? -1 : 1; 443 ++nsurf; 445 ++nsurf; 444 } 446 } 445 447 446 // Check Y faces 448 // Check Y faces 447 // 449 // 448 G4double ny = 0; 450 G4double ny = 0; 449 G4double yy = fPlanes[0].b*p.y()+fPlanes[0]. 451 G4double yy = fPlanes[0].b*p.y()+fPlanes[0].c*p.z(); 450 if (std::abs(fPlanes[0].d + yy) <= halfCarTo 452 if (std::abs(fPlanes[0].d + yy) <= halfCarTolerance) 451 { 453 { 452 ny = fPlanes[0].b; 454 ny = fPlanes[0].b; 453 nz += fPlanes[0].c; 455 nz += fPlanes[0].c; 454 ++nsurf; 456 ++nsurf; 455 } 457 } 456 else if (std::abs(fPlanes[1].d - yy) <= half 458 else if (std::abs(fPlanes[1].d - yy) <= halfCarTolerance) 457 { 459 { 458 ny = fPlanes[1].b; 460 ny = fPlanes[1].b; 459 nz += fPlanes[1].c; 461 nz += fPlanes[1].c; 460 ++nsurf; 462 ++nsurf; 461 } 463 } 462 464 463 // Check X faces 465 // Check X faces 464 // 466 // 465 G4double nx = 0; 467 G4double nx = 0; 466 G4double xx = fPlanes[2].a*p.x()+fPlanes[2]. 468 G4double xx = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()+fPlanes[2].c*p.z(); 467 if (std::abs(fPlanes[2].d + xx) <= halfCarTo 469 if (std::abs(fPlanes[2].d + xx) <= halfCarTolerance) 468 { 470 { 469 nx = fPlanes[2].a; 471 nx = fPlanes[2].a; 470 ny += fPlanes[2].b; 472 ny += fPlanes[2].b; 471 nz += fPlanes[2].c; 473 nz += fPlanes[2].c; 472 ++nsurf; 474 ++nsurf; 473 } 475 } 474 else if (std::abs(fPlanes[3].d - xx) <= half 476 else if (std::abs(fPlanes[3].d - xx) <= halfCarTolerance) 475 { 477 { 476 nx = fPlanes[3].a; 478 nx = fPlanes[3].a; 477 ny += fPlanes[3].b; 479 ny += fPlanes[3].b; 478 nz += fPlanes[3].c; 480 nz += fPlanes[3].c; 479 ++nsurf; 481 ++nsurf; 480 } 482 } 481 483 482 // Return normal 484 // Return normal 483 // 485 // 484 if (nsurf == 1) return {nx,ny,nz}; << 486 if (nsurf == 1) return G4ThreeVector(nx,ny,nz); 485 else if (nsurf != 0) return G4ThreeVector(nx 487 else if (nsurf != 0) return G4ThreeVector(nx,ny,nz).unit(); // edge or corner 486 else 488 else 487 { 489 { 488 // Point is not on the surface 490 // Point is not on the surface 489 // 491 // 490 #ifdef G4CSGDEBUG 492 #ifdef G4CSGDEBUG 491 std::ostringstream message; 493 std::ostringstream message; 492 G4int oldprc = message.precision(16); 494 G4int oldprc = message.precision(16); 493 message << "Point p is not on surface (!?) 495 message << "Point p is not on surface (!?) of solid: " 494 << GetName() << G4endl; 496 << GetName() << G4endl; 495 message << "Position:\n"; 497 message << "Position:\n"; 496 message << " p.x() = " << p.x()/mm << " 498 message << " p.x() = " << p.x()/mm << " mm\n"; 497 message << " p.y() = " << p.y()/mm << " 499 message << " p.y() = " << p.y()/mm << " mm\n"; 498 message << " p.z() = " << p.z()/mm << " 500 message << " p.z() = " << p.z()/mm << " mm"; 499 G4cout.precision(oldprc) ; 501 G4cout.precision(oldprc) ; 500 G4Exception("G4Para::SurfaceNormal(p)", "G 502 G4Exception("G4Para::SurfaceNormal(p)", "GeomSolids1002", 501 JustWarning, message ); 503 JustWarning, message ); 502 DumpInfo(); 504 DumpInfo(); 503 #endif 505 #endif 504 return ApproxSurfaceNormal(p); 506 return ApproxSurfaceNormal(p); 505 } 507 } 506 } 508 } 507 509 508 ////////////////////////////////////////////// 510 ////////////////////////////////////////////////////////////////////////// 509 // 511 // 510 // Algorithm for SurfaceNormal() following the 512 // Algorithm for SurfaceNormal() following the original specification 511 // for points not on the surface 513 // for points not on the surface 512 514 513 G4ThreeVector G4Para::ApproxSurfaceNormal( con 515 G4ThreeVector G4Para::ApproxSurfaceNormal( const G4ThreeVector& p ) const 514 { 516 { 515 G4double dist = -DBL_MAX; 517 G4double dist = -DBL_MAX; 516 G4int iside = 0; 518 G4int iside = 0; 517 for (G4int i=0; i<4; ++i) 519 for (G4int i=0; i<4; ++i) 518 { 520 { 519 G4double d = fPlanes[i].a*p.x() + 521 G4double d = fPlanes[i].a*p.x() + 520 fPlanes[i].b*p.y() + 522 fPlanes[i].b*p.y() + 521 fPlanes[i].c*p.z() + fPlanes[ 523 fPlanes[i].c*p.z() + fPlanes[i].d; 522 if (d > dist) { dist = d; iside = i; } 524 if (d > dist) { dist = d; iside = i; } 523 } 525 } 524 526 525 G4double distz = std::abs(p.z()) - fDz; 527 G4double distz = std::abs(p.z()) - fDz; 526 if (dist > distz) 528 if (dist > distz) 527 return { fPlanes[iside].a, fPlanes[iside]. << 529 return G4ThreeVector(fPlanes[iside].a, fPlanes[iside].b, fPlanes[iside].c); 528 else 530 else 529 return { 0, 0, (G4double)((p.z() < 0) ? -1 << 531 return G4ThreeVector(0, 0, (p.z() < 0) ? -1 : 1); 530 } 532 } 531 533 532 ////////////////////////////////////////////// 534 ////////////////////////////////////////////////////////////////////////// 533 // 535 // 534 // Calculate distance to shape from outside 536 // Calculate distance to shape from outside 535 // - return kInfinity if no intersection 537 // - return kInfinity if no intersection 536 538 537 G4double G4Para::DistanceToIn(const G4ThreeVec 539 G4double G4Para::DistanceToIn(const G4ThreeVector& p, 538 const G4ThreeVec 540 const G4ThreeVector& v ) const 539 { 541 { 540 // Z intersections 542 // Z intersections 541 // 543 // 542 if ((std::abs(p.z()) - fDz) >= -halfCarToler 544 if ((std::abs(p.z()) - fDz) >= -halfCarTolerance && p.z()*v.z() >= 0) 543 return kInfinity; 545 return kInfinity; 544 G4double invz = (-v.z() == 0) ? DBL_MAX : -1 546 G4double invz = (-v.z() == 0) ? DBL_MAX : -1./v.z(); 545 G4double dz = (invz < 0) ? fDz : -fDz; 547 G4double dz = (invz < 0) ? fDz : -fDz; 546 G4double tzmin = (p.z() + dz)*invz; 548 G4double tzmin = (p.z() + dz)*invz; 547 G4double tzmax = (p.z() - dz)*invz; 549 G4double tzmax = (p.z() - dz)*invz; 548 550 549 // Y intersections 551 // Y intersections 550 // 552 // 551 G4double tmin0 = tzmin, tmax0 = tzmax; 553 G4double tmin0 = tzmin, tmax0 = tzmax; 552 G4double cos0 = fPlanes[0].b*v.y() + fPlanes 554 G4double cos0 = fPlanes[0].b*v.y() + fPlanes[0].c*v.z(); 553 G4double disy = fPlanes[0].b*p.y() + fPlanes 555 G4double disy = fPlanes[0].b*p.y() + fPlanes[0].c*p.z(); 554 G4double dis0 = fPlanes[0].d + disy; 556 G4double dis0 = fPlanes[0].d + disy; 555 if (dis0 >= -halfCarTolerance) 557 if (dis0 >= -halfCarTolerance) 556 { 558 { 557 if (cos0 >= 0) return kInfinity; 559 if (cos0 >= 0) return kInfinity; 558 G4double tmp = -dis0/cos0; 560 G4double tmp = -dis0/cos0; 559 if (tmin0 < tmp) tmin0 = tmp; 561 if (tmin0 < tmp) tmin0 = tmp; 560 } 562 } 561 else if (cos0 > 0) 563 else if (cos0 > 0) 562 { 564 { 563 G4double tmp = -dis0/cos0; 565 G4double tmp = -dis0/cos0; 564 if (tmax0 > tmp) tmax0 = tmp; 566 if (tmax0 > tmp) tmax0 = tmp; 565 } 567 } 566 568 567 G4double tmin1 = tmin0, tmax1 = tmax0; 569 G4double tmin1 = tmin0, tmax1 = tmax0; 568 G4double cos1 = -cos0; 570 G4double cos1 = -cos0; 569 G4double dis1 = fPlanes[1].d - disy; 571 G4double dis1 = fPlanes[1].d - disy; 570 if (dis1 >= -halfCarTolerance) 572 if (dis1 >= -halfCarTolerance) 571 { 573 { 572 if (cos1 >= 0) return kInfinity; 574 if (cos1 >= 0) return kInfinity; 573 G4double tmp = -dis1/cos1; 575 G4double tmp = -dis1/cos1; 574 if (tmin1 < tmp) tmin1 = tmp; 576 if (tmin1 < tmp) tmin1 = tmp; 575 } 577 } 576 else if (cos1 > 0) 578 else if (cos1 > 0) 577 { 579 { 578 G4double tmp = -dis1/cos1; 580 G4double tmp = -dis1/cos1; 579 if (tmax1 > tmp) tmax1 = tmp; 581 if (tmax1 > tmp) tmax1 = tmp; 580 } 582 } 581 583 582 // X intersections 584 // X intersections 583 // 585 // 584 G4double tmin2 = tmin1, tmax2 = tmax1; 586 G4double tmin2 = tmin1, tmax2 = tmax1; 585 G4double cos2 = fPlanes[2].a*v.x() + fPlanes 587 G4double cos2 = fPlanes[2].a*v.x() + fPlanes[2].b*v.y() + fPlanes[2].c*v.z(); 586 G4double disx = fPlanes[2].a*p.x() + fPlanes 588 G4double disx = fPlanes[2].a*p.x() + fPlanes[2].b*p.y() + fPlanes[2].c*p.z(); 587 G4double dis2 = fPlanes[2].d + disx; 589 G4double dis2 = fPlanes[2].d + disx; 588 if (dis2 >= -halfCarTolerance) 590 if (dis2 >= -halfCarTolerance) 589 { 591 { 590 if (cos2 >= 0) return kInfinity; 592 if (cos2 >= 0) return kInfinity; 591 G4double tmp = -dis2/cos2; 593 G4double tmp = -dis2/cos2; 592 if (tmin2 < tmp) tmin2 = tmp; 594 if (tmin2 < tmp) tmin2 = tmp; 593 } 595 } 594 else if (cos2 > 0) 596 else if (cos2 > 0) 595 { 597 { 596 G4double tmp = -dis2/cos2; 598 G4double tmp = -dis2/cos2; 597 if (tmax2 > tmp) tmax2 = tmp; 599 if (tmax2 > tmp) tmax2 = tmp; 598 } 600 } 599 601 600 G4double tmin3 = tmin2, tmax3 = tmax2; 602 G4double tmin3 = tmin2, tmax3 = tmax2; 601 G4double cos3 = -cos2; 603 G4double cos3 = -cos2; 602 G4double dis3 = fPlanes[3].d - disx; 604 G4double dis3 = fPlanes[3].d - disx; 603 if (dis3 >= -halfCarTolerance) 605 if (dis3 >= -halfCarTolerance) 604 { 606 { 605 if (cos3 >= 0) return kInfinity; 607 if (cos3 >= 0) return kInfinity; 606 G4double tmp = -dis3/cos3; 608 G4double tmp = -dis3/cos3; 607 if (tmin3 < tmp) tmin3 = tmp; 609 if (tmin3 < tmp) tmin3 = tmp; 608 } 610 } 609 else if (cos3 > 0) 611 else if (cos3 > 0) 610 { 612 { 611 G4double tmp = -dis3/cos3; 613 G4double tmp = -dis3/cos3; 612 if (tmax3 > tmp) tmax3 = tmp; 614 if (tmax3 > tmp) tmax3 = tmp; 613 } 615 } 614 616 615 // Find distance 617 // Find distance 616 // 618 // 617 G4double tmin = tmin3, tmax = tmax3; 619 G4double tmin = tmin3, tmax = tmax3; 618 if (tmax <= tmin + halfCarTolerance) return 620 if (tmax <= tmin + halfCarTolerance) return kInfinity; // touch or no hit 619 return (tmin < halfCarTolerance ) ? 0. : tmi 621 return (tmin < halfCarTolerance ) ? 0. : tmin; 620 } 622 } 621 623 622 ////////////////////////////////////////////// 624 ////////////////////////////////////////////////////////////////////////// 623 // 625 // 624 // Calculate exact shortest distance to any bo 626 // Calculate exact shortest distance to any boundary from outside 625 // - returns 0 is point inside 627 // - returns 0 is point inside 626 628 627 G4double G4Para::DistanceToIn( const G4ThreeVe 629 G4double G4Para::DistanceToIn( const G4ThreeVector& p ) const 628 { 630 { 629 G4double xx = fPlanes[2].a*p.x()+fPlanes[2]. 631 G4double xx = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()+fPlanes[2].c*p.z(); 630 G4double dx = std::abs(xx) + fPlanes[2].d; 632 G4double dx = std::abs(xx) + fPlanes[2].d; 631 633 632 G4double yy = fPlanes[0].b*p.y()+fPlanes[0]. 634 G4double yy = fPlanes[0].b*p.y()+fPlanes[0].c*p.z(); 633 G4double dy = std::abs(yy) + fPlanes[0].d; 635 G4double dy = std::abs(yy) + fPlanes[0].d; 634 G4double dxy = std::max(dx,dy); 636 G4double dxy = std::max(dx,dy); 635 637 636 G4double dz = std::abs(p.z())-fDz; 638 G4double dz = std::abs(p.z())-fDz; 637 G4double dist = std::max(dxy,dz); 639 G4double dist = std::max(dxy,dz); 638 640 639 return (dist > 0) ? dist : 0.; 641 return (dist > 0) ? dist : 0.; 640 } 642 } 641 643 642 ////////////////////////////////////////////// 644 ////////////////////////////////////////////////////////////////////////// 643 // 645 // 644 // Calculate distance to surface of shape from 646 // Calculate distance to surface of shape from inside and, if required, 645 // find normal at exit point 647 // find normal at exit point 646 // - when leaving the surface, return 0 648 // - when leaving the surface, return 0 647 649 648 G4double G4Para::DistanceToOut(const G4ThreeVe 650 G4double G4Para::DistanceToOut(const G4ThreeVector& p, const G4ThreeVector& v, 649 const G4bool ca 651 const G4bool calcNorm, 650 G4bool* v 652 G4bool* validNorm, G4ThreeVector* n) const 651 { 653 { 652 // Z intersections 654 // Z intersections 653 // 655 // 654 if ((std::abs(p.z()) - fDz) >= -halfCarToler 656 if ((std::abs(p.z()) - fDz) >= -halfCarTolerance && p.z()*v.z() > 0) 655 { 657 { 656 if (calcNorm) 658 if (calcNorm) 657 { 659 { 658 *validNorm = true; 660 *validNorm = true; 659 n->set(0, 0, (p.z() < 0) ? -1 : 1); 661 n->set(0, 0, (p.z() < 0) ? -1 : 1); 660 } 662 } 661 return 0.; 663 return 0.; 662 } 664 } 663 G4double vz = v.z(); 665 G4double vz = v.z(); 664 G4double tmax = (vz == 0) ? DBL_MAX : (std:: 666 G4double tmax = (vz == 0) ? DBL_MAX : (std::copysign(fDz,vz) - p.z())/vz; 665 G4int iside = (vz < 0) ? -4 : -2; // little 667 G4int iside = (vz < 0) ? -4 : -2; // little trick: (-4+3)=-1, (-2+3)=+1 666 668 667 // Y intersections 669 // Y intersections 668 // 670 // 669 G4double cos0 = fPlanes[0].b*v.y() + fPlanes 671 G4double cos0 = fPlanes[0].b*v.y() + fPlanes[0].c*v.z(); 670 if (cos0 > 0) 672 if (cos0 > 0) 671 { 673 { 672 G4double dis0 = fPlanes[0].b*p.y() + fPlan 674 G4double dis0 = fPlanes[0].b*p.y() + fPlanes[0].c*p.z() + fPlanes[0].d; 673 if (dis0 >= -halfCarTolerance) 675 if (dis0 >= -halfCarTolerance) 674 { 676 { 675 if (calcNorm) 677 if (calcNorm) 676 { 678 { 677 *validNorm = true; 679 *validNorm = true; 678 n->set(0, fPlanes[0].b, fPlanes[0].c); 680 n->set(0, fPlanes[0].b, fPlanes[0].c); 679 } 681 } 680 return 0.; 682 return 0.; 681 } 683 } 682 G4double tmp = -dis0/cos0; 684 G4double tmp = -dis0/cos0; 683 if (tmax > tmp) { tmax = tmp; iside = 0; } 685 if (tmax > tmp) { tmax = tmp; iside = 0; } 684 } 686 } 685 687 686 G4double cos1 = -cos0; 688 G4double cos1 = -cos0; 687 if (cos1 > 0) 689 if (cos1 > 0) 688 { 690 { 689 G4double dis1 = fPlanes[1].b*p.y() + fPlan 691 G4double dis1 = fPlanes[1].b*p.y() + fPlanes[1].c*p.z() + fPlanes[1].d; 690 if (dis1 >= -halfCarTolerance) 692 if (dis1 >= -halfCarTolerance) 691 { 693 { 692 if (calcNorm) 694 if (calcNorm) 693 { 695 { 694 *validNorm = true; 696 *validNorm = true; 695 n->set(0, fPlanes[1].b, fPlanes[1].c); 697 n->set(0, fPlanes[1].b, fPlanes[1].c); 696 } 698 } 697 return 0.; 699 return 0.; 698 } 700 } 699 G4double tmp = -dis1/cos1; 701 G4double tmp = -dis1/cos1; 700 if (tmax > tmp) { tmax = tmp; iside = 1; } 702 if (tmax > tmp) { tmax = tmp; iside = 1; } 701 } 703 } 702 704 703 // X intersections 705 // X intersections 704 // 706 // 705 G4double cos2 = fPlanes[2].a*v.x() + fPlanes 707 G4double cos2 = fPlanes[2].a*v.x() + fPlanes[2].b*v.y() + fPlanes[2].c*v.z(); 706 if (cos2 > 0) 708 if (cos2 > 0) 707 { 709 { 708 G4double dis2 = fPlanes[2].a*p.x()+fPlanes 710 G4double dis2 = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()+fPlanes[2].c*p.z()+fPlanes[2].d; 709 if (dis2 >= -halfCarTolerance) 711 if (dis2 >= -halfCarTolerance) 710 { 712 { 711 if (calcNorm) 713 if (calcNorm) 712 { 714 { 713 *validNorm = true; 715 *validNorm = true; 714 n->set(fPlanes[2].a, fPlanes[2].b, fP 716 n->set(fPlanes[2].a, fPlanes[2].b, fPlanes[2].c); 715 } 717 } 716 return 0.; 718 return 0.; 717 } 719 } 718 G4double tmp = -dis2/cos2; 720 G4double tmp = -dis2/cos2; 719 if (tmax > tmp) { tmax = tmp; iside = 2; } 721 if (tmax > tmp) { tmax = tmp; iside = 2; } 720 } 722 } 721 723 722 G4double cos3 = -cos2; 724 G4double cos3 = -cos2; 723 if (cos3 > 0) 725 if (cos3 > 0) 724 { 726 { 725 G4double dis3 = fPlanes[3].a*p.x()+fPlanes 727 G4double dis3 = fPlanes[3].a*p.x()+fPlanes[3].b*p.y()+fPlanes[3].c*p.z()+fPlanes[3].d; 726 if (dis3 >= -halfCarTolerance) 728 if (dis3 >= -halfCarTolerance) 727 { 729 { 728 if (calcNorm) 730 if (calcNorm) 729 { 731 { 730 *validNorm = true; 732 *validNorm = true; 731 n->set(fPlanes[3].a, fPlanes[3].b, fP 733 n->set(fPlanes[3].a, fPlanes[3].b, fPlanes[3].c); 732 } 734 } 733 return 0.; 735 return 0.; 734 } 736 } 735 G4double tmp = -dis3/cos3; 737 G4double tmp = -dis3/cos3; 736 if (tmax > tmp) { tmax = tmp; iside = 3; } 738 if (tmax > tmp) { tmax = tmp; iside = 3; } 737 } 739 } 738 740 739 // Set normal, if required, and return dista 741 // Set normal, if required, and return distance 740 // 742 // 741 if (calcNorm) 743 if (calcNorm) 742 { 744 { 743 *validNorm = true; 745 *validNorm = true; 744 if (iside < 0) 746 if (iside < 0) 745 n->set(0, 0, iside + 3); // (-4+3)=-1, ( 747 n->set(0, 0, iside + 3); // (-4+3)=-1, (-2+3)=+1 746 else 748 else 747 n->set(fPlanes[iside].a, fPlanes[iside]. 749 n->set(fPlanes[iside].a, fPlanes[iside].b, fPlanes[iside].c); 748 } 750 } 749 return tmax; 751 return tmax; 750 } 752 } 751 753 752 ////////////////////////////////////////////// 754 ////////////////////////////////////////////////////////////////////////// 753 // 755 // 754 // Calculate exact shortest distance to any bo 756 // Calculate exact shortest distance to any boundary from inside 755 // - returns 0 is point outside 757 // - returns 0 is point outside 756 758 757 G4double G4Para::DistanceToOut( const G4ThreeV 759 G4double G4Para::DistanceToOut( const G4ThreeVector& p ) const 758 { 760 { 759 #ifdef G4CSGDEBUG 761 #ifdef G4CSGDEBUG 760 if( Inside(p) == kOutside ) 762 if( Inside(p) == kOutside ) 761 { 763 { 762 std::ostringstream message; 764 std::ostringstream message; 763 G4int oldprc = message.precision(16); 765 G4int oldprc = message.precision(16); 764 message << "Point p is outside (!?) of sol 766 message << "Point p is outside (!?) of solid: " << GetName() << G4endl; 765 message << "Position:\n"; 767 message << "Position:\n"; 766 message << " p.x() = " << p.x()/mm << " 768 message << " p.x() = " << p.x()/mm << " mm\n"; 767 message << " p.y() = " << p.y()/mm << " 769 message << " p.y() = " << p.y()/mm << " mm\n"; 768 message << " p.z() = " << p.z()/mm << " 770 message << " p.z() = " << p.z()/mm << " mm"; 769 G4cout.precision(oldprc) ; 771 G4cout.precision(oldprc) ; 770 G4Exception("G4Para::DistanceToOut(p)", "G 772 G4Exception("G4Para::DistanceToOut(p)", "GeomSolids1002", 771 JustWarning, message ); 773 JustWarning, message ); 772 DumpInfo(); 774 DumpInfo(); 773 } 775 } 774 #endif 776 #endif 775 G4double xx = fPlanes[2].a*p.x()+fPlanes[2]. 777 G4double xx = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()+fPlanes[2].c*p.z(); 776 G4double dx = std::abs(xx) + fPlanes[2].d; 778 G4double dx = std::abs(xx) + fPlanes[2].d; 777 779 778 G4double yy = fPlanes[0].b*p.y()+fPlanes[0]. 780 G4double yy = fPlanes[0].b*p.y()+fPlanes[0].c*p.z(); 779 G4double dy = std::abs(yy) + fPlanes[0].d; 781 G4double dy = std::abs(yy) + fPlanes[0].d; 780 G4double dxy = std::max(dx,dy); 782 G4double dxy = std::max(dx,dy); 781 783 782 G4double dz = std::abs(p.z())-fDz; 784 G4double dz = std::abs(p.z())-fDz; 783 G4double dist = std::max(dxy,dz); 785 G4double dist = std::max(dxy,dz); 784 786 785 return (dist < 0) ? -dist : 0.; 787 return (dist < 0) ? -dist : 0.; 786 } 788 } 787 789 788 ////////////////////////////////////////////// 790 ////////////////////////////////////////////////////////////////////////// 789 // 791 // 790 // GetEntityType 792 // GetEntityType 791 793 792 G4GeometryType G4Para::GetEntityType() const 794 G4GeometryType G4Para::GetEntityType() const 793 { 795 { 794 return {"G4Para"}; << 796 return G4String("G4Para"); 795 } << 796 << 797 ////////////////////////////////////////////// << 798 // << 799 // IsFaceted << 800 << 801 G4bool G4Para::IsFaceted() const << 802 { << 803 return true; << 804 } 797 } 805 798 806 ////////////////////////////////////////////// 799 ////////////////////////////////////////////////////////////////////////// 807 // 800 // 808 // Make a clone of the object 801 // Make a clone of the object 809 // 802 // 810 G4VSolid* G4Para::Clone() const 803 G4VSolid* G4Para::Clone() const 811 { 804 { 812 return new G4Para(*this); 805 return new G4Para(*this); 813 } 806 } 814 807 815 ////////////////////////////////////////////// 808 ////////////////////////////////////////////////////////////////////////// 816 // 809 // 817 // Stream object contents to an output stream 810 // Stream object contents to an output stream 818 811 819 std::ostream& G4Para::StreamInfo( std::ostream 812 std::ostream& G4Para::StreamInfo( std::ostream& os ) const 820 { 813 { 821 G4double alpha = std::atan(fTalpha); 814 G4double alpha = std::atan(fTalpha); 822 G4double theta = std::atan(std::sqrt(fTtheta 815 G4double theta = std::atan(std::sqrt(fTthetaCphi*fTthetaCphi + 823 fTtheta 816 fTthetaSphi*fTthetaSphi)); 824 G4double phi = std::atan2(fTthetaSphi,fTth 817 G4double phi = std::atan2(fTthetaSphi,fTthetaCphi); 825 818 826 G4long oldprc = os.precision(16); << 819 G4int oldprc = os.precision(16); 827 os << "------------------------------------- 820 os << "-----------------------------------------------------------\n" 828 << " *** Dump for solid - " << GetName 821 << " *** Dump for solid - " << GetName() << " ***\n" 829 << " ================================= 822 << " ===================================================\n" 830 << " Solid type: G4Para\n" 823 << " Solid type: G4Para\n" 831 << " Parameters:\n" 824 << " Parameters:\n" 832 << " half length X: " << fDx/mm << " m 825 << " half length X: " << fDx/mm << " mm\n" 833 << " half length Y: " << fDy/mm << " m 826 << " half length Y: " << fDy/mm << " mm\n" 834 << " half length Z: " << fDz/mm << " m 827 << " half length Z: " << fDz/mm << " mm\n" 835 << " alpha: " << alpha/degree << "degr 828 << " alpha: " << alpha/degree << "degrees\n" 836 << " theta: " << theta/degree << "degr 829 << " theta: " << theta/degree << "degrees\n" 837 << " phi: " << phi/degree << "degrees\ 830 << " phi: " << phi/degree << "degrees\n" 838 << "------------------------------------- 831 << "-----------------------------------------------------------\n"; 839 os.precision(oldprc); 832 os.precision(oldprc); 840 833 841 return os; 834 return os; 842 } 835 } 843 836 844 ////////////////////////////////////////////// 837 ////////////////////////////////////////////////////////////////////////// 845 // 838 // 846 // Return a point randomly and uniformly selec 839 // Return a point randomly and uniformly selected on the solid surface 847 840 848 G4ThreeVector G4Para::GetPointOnSurface() cons 841 G4ThreeVector G4Para::GetPointOnSurface() const 849 { 842 { 850 G4double DyTalpha = fDy*fTalpha; 843 G4double DyTalpha = fDy*fTalpha; 851 G4double DzTthetaSphi = fDz*fTthetaSphi; 844 G4double DzTthetaSphi = fDz*fTthetaSphi; 852 G4double DzTthetaCphi = fDz*fTthetaCphi; 845 G4double DzTthetaCphi = fDz*fTthetaCphi; 853 846 854 // Set vertices 847 // Set vertices 855 // 848 // 856 G4ThreeVector pt[8]; 849 G4ThreeVector pt[8]; 857 pt[0].set(-DzTthetaCphi-DyTalpha-fDx, -DzTth 850 pt[0].set(-DzTthetaCphi-DyTalpha-fDx, -DzTthetaSphi-fDy, -fDz); 858 pt[1].set(-DzTthetaCphi-DyTalpha+fDx, -DzTth 851 pt[1].set(-DzTthetaCphi-DyTalpha+fDx, -DzTthetaSphi-fDy, -fDz); 859 pt[2].set(-DzTthetaCphi+DyTalpha-fDx, -DzTth 852 pt[2].set(-DzTthetaCphi+DyTalpha-fDx, -DzTthetaSphi+fDy, -fDz); 860 pt[3].set(-DzTthetaCphi+DyTalpha+fDx, -DzTth 853 pt[3].set(-DzTthetaCphi+DyTalpha+fDx, -DzTthetaSphi+fDy, -fDz); 861 pt[4].set( DzTthetaCphi-DyTalpha-fDx, DzTth 854 pt[4].set( DzTthetaCphi-DyTalpha-fDx, DzTthetaSphi-fDy, fDz); 862 pt[5].set( DzTthetaCphi-DyTalpha+fDx, DzTth 855 pt[5].set( DzTthetaCphi-DyTalpha+fDx, DzTthetaSphi-fDy, fDz); 863 pt[6].set( DzTthetaCphi+DyTalpha-fDx, DzTth 856 pt[6].set( DzTthetaCphi+DyTalpha-fDx, DzTthetaSphi+fDy, fDz); 864 pt[7].set( DzTthetaCphi+DyTalpha+fDx, DzTth 857 pt[7].set( DzTthetaCphi+DyTalpha+fDx, DzTthetaSphi+fDy, fDz); 865 858 866 // Set areas (-Z, -Y, +Y, -X, +X, +Z) 859 // Set areas (-Z, -Y, +Y, -X, +X, +Z) 867 // 860 // 868 G4ThreeVector vx(fDx, 0, 0); 861 G4ThreeVector vx(fDx, 0, 0); 869 G4ThreeVector vy(DyTalpha, fDy, 0); 862 G4ThreeVector vy(DyTalpha, fDy, 0); 870 G4ThreeVector vz(DzTthetaCphi, DzTthetaSphi, 863 G4ThreeVector vz(DzTthetaCphi, DzTthetaSphi, fDz); 871 864 872 G4double sxy = fDx*fDy; // (vx.cross(vy)).ma 865 G4double sxy = fDx*fDy; // (vx.cross(vy)).mag(); 873 G4double sxz = (vx.cross(vz)).mag(); 866 G4double sxz = (vx.cross(vz)).mag(); 874 G4double syz = (vy.cross(vz)).mag(); 867 G4double syz = (vy.cross(vz)).mag(); 875 868 876 G4double sface[6] = { sxy, syz, syz, sxz, sx 869 G4double sface[6] = { sxy, syz, syz, sxz, sxz, sxy }; 877 for (G4int i=1; i<6; ++i) { sface[i] += sfac 870 for (G4int i=1; i<6; ++i) { sface[i] += sface[i-1]; } 878 871 879 // Select face 872 // Select face 880 // 873 // 881 G4double select = sface[5]*G4UniformRand(); 874 G4double select = sface[5]*G4UniformRand(); 882 G4int k = 5; 875 G4int k = 5; 883 if (select <= sface[4]) k = 4; 876 if (select <= sface[4]) k = 4; 884 if (select <= sface[3]) k = 3; 877 if (select <= sface[3]) k = 3; 885 if (select <= sface[2]) k = 2; 878 if (select <= sface[2]) k = 2; 886 if (select <= sface[1]) k = 1; 879 if (select <= sface[1]) k = 1; 887 if (select <= sface[0]) k = 0; 880 if (select <= sface[0]) k = 0; 888 881 889 // Generate point 882 // Generate point 890 // 883 // 891 G4int ip[6][3] = {{0,1,2}, {0,4,1}, {2,3,6}, 884 G4int ip[6][3] = {{0,1,2}, {0,4,1}, {2,3,6}, {0,2,4}, {1,5,3}, {4,6,5}}; 892 G4double u = G4UniformRand(); 885 G4double u = G4UniformRand(); 893 G4double v = G4UniformRand(); 886 G4double v = G4UniformRand(); 894 return (1.-u-v)*pt[ip[k][0]] + u*pt[ip[k][1] 887 return (1.-u-v)*pt[ip[k][0]] + u*pt[ip[k][1]] + v*pt[ip[k][2]]; 895 } 888 } 896 889 897 ////////////////////////////////////////////// 890 ////////////////////////////////////////////////////////////////////////// 898 // 891 // 899 // Methods for visualisation 892 // Methods for visualisation 900 893 901 void G4Para::DescribeYourselfTo ( G4VGraphicsS 894 void G4Para::DescribeYourselfTo ( G4VGraphicsScene& scene ) const 902 { 895 { 903 scene.AddSolid (*this); 896 scene.AddSolid (*this); 904 } 897 } 905 898 906 G4Polyhedron* G4Para::CreatePolyhedron () cons 899 G4Polyhedron* G4Para::CreatePolyhedron () const 907 { 900 { 908 G4double phi = std::atan2(fTthetaSphi, fTthe 901 G4double phi = std::atan2(fTthetaSphi, fTthetaCphi); 909 G4double alpha = std::atan(fTalpha); 902 G4double alpha = std::atan(fTalpha); 910 G4double theta = std::atan(std::sqrt(fTtheta 903 G4double theta = std::atan(std::sqrt(fTthetaCphi*fTthetaCphi + 911 fTtheta 904 fTthetaSphi*fTthetaSphi)); 912 905 913 return new G4PolyhedronPara(fDx, fDy, fDz, a 906 return new G4PolyhedronPara(fDx, fDy, fDz, alpha, theta, phi); 914 } 907 } 915 #endif 908 #endif 916 909