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 // >> 27 // $Id$ >> 28 // >> 29 // 26 // Implementation for G4Trd class 30 // Implementation for G4Trd class 27 // 31 // 28 // 12.01.95 P.Kent: First version << 32 // History: 29 // 28.04.05 V.Grichine: new SurfaceNormal acco << 33 // 30 // 25.05.17 E.Tcherniaev: complete revision, s << 34 // 28.04.05 V.Grichine: new SurfaceNormal according to J. Apostolakis proposal 31 // ------------------------------------------- << 35 // 26.04.05, V.Grichine, new SurfaceNoramal is default >> 36 // 07.12.04, V.Grichine, SurfaceNoramal with edges/vertices. >> 37 // 07.05.00, V.Grichine, in d = DistanceToIn(p,v), if d<0.5*kCarTolerance, d=0 >> 38 // ~1996, V.Grichine, 1st implementation based on old code of P.Kent >> 39 // >> 40 ////////////////////////////////////////////////////////////////////////////// 32 41 33 #include "G4Trd.hh" 42 #include "G4Trd.hh" 34 43 35 #if !defined(G4GEOM_USE_UTRD) << 44 #include "G4VPVParameterisation.hh" 36 << 37 #include "G4GeomTools.hh" << 38 << 39 #include "G4VoxelLimits.hh" 45 #include "G4VoxelLimits.hh" 40 #include "G4AffineTransform.hh" 46 #include "G4AffineTransform.hh" 41 #include "G4BoundingEnvelope.hh" << 47 #include "Randomize.hh" 42 #include "G4QuickRand.hh" << 43 << 44 #include "G4VPVParameterisation.hh" << 45 48 46 #include "G4VGraphicsScene.hh" 49 #include "G4VGraphicsScene.hh" >> 50 #include "G4Polyhedron.hh" >> 51 #include "G4NURBS.hh" >> 52 #include "G4NURBSbox.hh" 47 53 48 using namespace CLHEP; 54 using namespace CLHEP; 49 55 50 ////////////////////////////////////////////// << 56 ///////////////////////////////////////////////////////////////////////// 51 // 57 // 52 // Constructor - set & check half widths << 58 // Constructor - check & set half widths 53 59 54 G4Trd::G4Trd(const G4String& pName, << 60 G4Trd::G4Trd( const G4String& pName, 55 G4double pdx1, G4double pdx << 61 G4double pdx1, G4double pdx2, 56 G4double pdy1, G4double pdy << 62 G4double pdy1, G4double pdy2, 57 G4double pdz) << 63 G4double pdz ) 58 : G4CSGSolid(pName), halfCarTolerance(0.5*kC << 64 : G4CSGSolid(pName) 59 fDx1(pdx1), fDx2(pdx2), fDy1(pdy1), fDy2(p << 60 { 65 { 61 CheckParameters(); << 66 CheckAndSetAllParameters (pdx1, pdx2, pdy1, pdy2, pdz); 62 MakePlanes(); << 63 } 67 } 64 68 65 ////////////////////////////////////////////// << 69 ///////////////////////////////////////////////////////////////////////// >> 70 // >> 71 // Set and check (coplanarity) of trd parameters >> 72 >> 73 void G4Trd::CheckAndSetAllParameters ( G4double pdx1, G4double pdx2, >> 74 G4double pdy1, G4double pdy2, >> 75 G4double pdz ) >> 76 { >> 77 if ( pdx1>0&&pdx2>0&&pdy1>0&&pdy2>0&&pdz>0 ) >> 78 { >> 79 fDx1=pdx1; fDx2=pdx2; >> 80 fDy1=pdy1; fDy2=pdy2; >> 81 fDz=pdz; >> 82 } >> 83 else >> 84 { >> 85 if ( pdx1>=0 && pdx2>=0 && pdy1>=0 && pdy2>=0 && pdz>=0 ) >> 86 { >> 87 // G4double Minimum_length= (1+per_thousand) * kCarTolerance/2.; >> 88 // FIX-ME : temporary solution for ZERO or very-small parameters >> 89 // >> 90 G4double Minimum_length= kCarTolerance/2.; >> 91 fDx1=std::max(pdx1,Minimum_length); >> 92 fDx2=std::max(pdx2,Minimum_length); >> 93 fDy1=std::max(pdy1,Minimum_length); >> 94 fDy2=std::max(pdy2,Minimum_length); >> 95 fDz=std::max(pdz,Minimum_length); >> 96 } >> 97 else >> 98 { >> 99 std::ostringstream message; >> 100 message << "Invalid negative dimensions for Solid: " << GetName() >> 101 << G4endl >> 102 << " X - " << pdx1 << ", " << pdx2 << G4endl >> 103 << " Y - " << pdy1 << ", " << pdy2 << G4endl >> 104 << " Z - " << pdz; >> 105 G4Exception("G4Trd::CheckAndSetAllParameters()", >> 106 "GeomSolids0002", FatalException, message); >> 107 } >> 108 } >> 109 fCubicVolume= 0.; >> 110 fSurfaceArea= 0.; >> 111 fpPolyhedron = 0; >> 112 } >> 113 >> 114 /////////////////////////////////////////////////////////////////////// 66 // 115 // 67 // Fake default constructor - sets only member 116 // Fake default constructor - sets only member data and allocates memory 68 // for usage restri << 117 // for usage restricted to object persistency. 69 // 118 // 70 G4Trd::G4Trd( __void__& a ) 119 G4Trd::G4Trd( __void__& a ) 71 : G4CSGSolid(a), halfCarTolerance(0.5*kCarTo << 120 : G4CSGSolid(a), fDx1(0.), fDx2(0.), fDy1(0.), fDy2(0.), fDz(0.) 72 fDx1(1.), fDx2(1.), fDy1(1.), fDy2(1.), fD << 73 { 121 { 74 MakePlanes(); << 75 } 122 } 76 123 77 ////////////////////////////////////////////// 124 ////////////////////////////////////////////////////////////////////////// 78 // 125 // 79 // Destructor 126 // Destructor 80 127 81 G4Trd::~G4Trd() = default; << 128 G4Trd::~G4Trd() >> 129 { >> 130 } 82 131 83 ////////////////////////////////////////////// 132 ////////////////////////////////////////////////////////////////////////// 84 // 133 // 85 // Copy constructor 134 // Copy constructor 86 135 87 G4Trd::G4Trd(const G4Trd& rhs) 136 G4Trd::G4Trd(const G4Trd& rhs) 88 : G4CSGSolid(rhs), halfCarTolerance(rhs.half << 137 : G4CSGSolid(rhs), fDx1(rhs.fDx1), fDx2(rhs.fDx2), 89 fDx1(rhs.fDx1), fDx2(rhs.fDx2), << 138 fDy1(rhs.fDy1), fDy2(rhs.fDy2), fDz(rhs.fDz) 90 fDy1(rhs.fDy1), fDy2(rhs.fDy2), fDz(rhs.fD << 91 fHx(rhs.fHx), fHy(rhs.fHy) << 92 { 139 { 93 for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs << 94 } 140 } 95 141 96 ////////////////////////////////////////////// 142 ////////////////////////////////////////////////////////////////////////// 97 // 143 // 98 // Assignment operator 144 // Assignment operator 99 145 100 G4Trd& G4Trd::operator = (const G4Trd& rhs) << 146 G4Trd& G4Trd::operator = (const G4Trd& rhs) 101 { 147 { 102 // Check assignment to self 148 // Check assignment to self 103 // 149 // 104 if (this == &rhs) { return *this; } 150 if (this == &rhs) { return *this; } 105 151 106 // Copy base class data 152 // Copy base class data 107 // 153 // 108 G4CSGSolid::operator=(rhs); 154 G4CSGSolid::operator=(rhs); 109 155 110 // Copy data 156 // Copy data 111 // 157 // 112 halfCarTolerance = rhs.halfCarTolerance; << 113 fDx1 = rhs.fDx1; fDx2 = rhs.fDx2; 158 fDx1 = rhs.fDx1; fDx2 = rhs.fDx2; 114 fDy1 = rhs.fDy1; fDy2 = rhs.fDy2; 159 fDy1 = rhs.fDy1; fDy2 = rhs.fDy2; 115 fDz = rhs.fDz; 160 fDz = rhs.fDz; 116 fHx = rhs.fHx; fHy = rhs.fHy; << 117 for (G4int i=0; i<4; ++i) { fPlanes[i] = rh << 118 161 119 return *this; 162 return *this; 120 } 163 } 121 164 122 ////////////////////////////////////////////// << 165 //////////////////////////////////////////////////////////////////////////// 123 // 166 // 124 // Set all parameters, as for constructor - se << 125 << 126 void G4Trd::SetAllParameters(G4double pdx1, G4 << 127 G4double pdy1, G4 << 128 { << 129 // Reset data of the base class << 130 fCubicVolume = 0.; << 131 fSurfaceArea = 0.; << 132 fRebuildPolyhedron = true; << 133 << 134 // Set parameters << 135 fDx1 = pdx1; fDx2 = pdx2; << 136 fDy1 = pdy1; fDy2 = pdy2; << 137 fDz = pdz; << 138 << 139 CheckParameters(); << 140 MakePlanes(); << 141 } << 142 << 143 ////////////////////////////////////////////// << 144 // 167 // 145 // Check dimensions << 146 168 147 void G4Trd::CheckParameters() << 169 void G4Trd::SetAllParameters ( G4double pdx1, G4double pdx2, G4double pdy1, >> 170 G4double pdy2, G4double pdz ) 148 { 171 { 149 G4double dmin = 2*kCarTolerance; << 172 CheckAndSetAllParameters (pdx1, pdx2, pdy1, pdy2, pdz); 150 if ((fDx1 < 0 || fDx2 < 0 || fDy1 < 0 || fDy << 151 (fDx1 < dmin && fDx2 < dmin) || << 152 (fDy1 < dmin && fDy2 < dmin)) << 153 { << 154 std::ostringstream message; << 155 message << "Invalid (too small or negative << 156 << GetName() << 157 << "\n X - " << fDx1 << ", " << f << 158 << "\n Y - " << fDy1 << ", " << f << 159 << "\n Z - " << fDz; << 160 G4Exception("G4Trd::CheckParameters()", "G << 161 FatalException, message); << 162 } << 163 } 173 } 164 174 165 ////////////////////////////////////////////// << 166 // << 167 // Set side planes << 168 << 169 void G4Trd::MakePlanes() << 170 { << 171 G4double dx = fDx1 - fDx2; << 172 G4double dy = fDy1 - fDy2; << 173 G4double dz = 2*fDz; << 174 fHx = std::sqrt(dy*dy + dz*dz); << 175 fHy = std::sqrt(dx*dx + dz*dz); << 176 175 177 // Set X planes at -Y & +Y << 176 ///////////////////////////////////////////////////////////////////////// 178 // << 179 fPlanes[0].a = 0.; << 180 fPlanes[0].b = -dz/fHx; << 181 fPlanes[0].c = dy/fHx; << 182 fPlanes[0].d = fPlanes[0].b*fDy1 + fPlanes[0 << 183 << 184 fPlanes[1].a = fPlanes[0].a; << 185 fPlanes[1].b = -fPlanes[0].b; << 186 fPlanes[1].c = fPlanes[0].c; << 187 fPlanes[1].d = fPlanes[0].d; << 188 << 189 // Set Y planes at -X & +X << 190 // << 191 fPlanes[2].a = -dz/fHy; << 192 fPlanes[2].b = 0.; << 193 fPlanes[2].c = dx/fHy; << 194 fPlanes[2].d = fPlanes[2].a*fDx1 + fPlanes[2 << 195 << 196 fPlanes[3].a = -fPlanes[2].a; << 197 fPlanes[3].b = fPlanes[2].b; << 198 fPlanes[3].c = fPlanes[2].c; << 199 fPlanes[3].d = fPlanes[2].d; << 200 } << 201 << 202 ////////////////////////////////////////////// << 203 // << 204 // Get volume << 205 << 206 G4double G4Trd::GetCubicVolume() << 207 { << 208 if (fCubicVolume == 0.) << 209 { << 210 fCubicVolume = 2*fDz*( (fDx1+fDx2)*(fDy1+f << 211 (fDx2-fDx1)*(fDy2-f << 212 } << 213 return fCubicVolume; << 214 } << 215 << 216 ////////////////////////////////////////////// << 217 // << 218 // Get surface area << 219 << 220 G4double G4Trd::GetSurfaceArea() << 221 { << 222 if (fSurfaceArea == 0.) << 223 { << 224 fSurfaceArea = << 225 4*(fDx1*fDy1 + fDx2*fDy2) + 2*(fDx1+fDx2 << 226 } << 227 return fSurfaceArea; << 228 } << 229 << 230 ////////////////////////////////////////////// << 231 // 177 // 232 // Dispatch to parameterisation for replicatio 178 // Dispatch to parameterisation for replication mechanism dimension 233 // computation & modification << 179 // computation & modification. 234 180 235 void G4Trd::ComputeDimensions( G4VPVPara 181 void G4Trd::ComputeDimensions( G4VPVParameterisation* p, 236 const G4int n, 182 const G4int n, 237 const G4VPhysic 183 const G4VPhysicalVolume* pRep ) 238 { 184 { 239 p->ComputeDimensions(*this,n,pRep); 185 p->ComputeDimensions(*this,n,pRep); 240 } 186 } 241 187 242 ////////////////////////////////////////////// << 243 // << 244 // Get bounding box << 245 188 246 void G4Trd::BoundingLimits(G4ThreeVector& pMin << 189 /////////////////////////////////////////////////////////////////////////// 247 { << 248 G4double dx1 = GetXHalfLength1(); << 249 G4double dx2 = GetXHalfLength2(); << 250 G4double dy1 = GetYHalfLength1(); << 251 G4double dy2 = GetYHalfLength2(); << 252 G4double dz = GetZHalfLength(); << 253 << 254 G4double xmax = std::max(dx1,dx2); << 255 G4double ymax = std::max(dy1,dy2); << 256 pMin.set(-xmax,-ymax,-dz); << 257 pMax.set( xmax, ymax, dz); << 258 << 259 // Check correctness of the bounding box << 260 // << 261 if (pMin.x() >= pMax.x() || pMin.y() >= pMax << 262 { << 263 std::ostringstream message; << 264 message << "Bad bounding box (min >= max) << 265 << GetName() << " !" << 266 << "\npMin = " << pMin << 267 << "\npMax = " << pMax; << 268 G4Exception("G4Trd::BoundingLimits()", "Ge << 269 DumpInfo(); << 270 } << 271 } << 272 << 273 ////////////////////////////////////////////// << 274 // 190 // 275 // Calculate extent under transform and specif 191 // Calculate extent under transform and specified limit 276 192 277 G4bool G4Trd::CalculateExtent( const EAxis pAx 193 G4bool G4Trd::CalculateExtent( const EAxis pAxis, 278 const G4VoxelLi 194 const G4VoxelLimits& pVoxelLimit, 279 const G4AffineT 195 const G4AffineTransform& pTransform, 280 G4double& 196 G4double& pMin, G4double& pMax ) const 281 { 197 { 282 G4ThreeVector bmin, bmax; << 198 if (!pTransform.IsRotated()) 283 G4bool exist; << 284 << 285 // Check bounding box (bbox) << 286 // << 287 BoundingLimits(bmin,bmax); << 288 G4BoundingEnvelope bbox(bmin,bmax); << 289 #ifdef G4BBOX_EXTENT << 290 return bbox.CalculateExtent(pAxis,pVoxelLimi << 291 #endif << 292 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox << 293 { 199 { 294 return exist = pMin < pMax; << 200 // Special case handling for unrotated solids >> 201 // Compute x/y/z mins and maxs respecting limits, with early returns >> 202 // if outside limits. Then switch() on pAxis >> 203 >> 204 G4double xoffset,xMin,xMax; >> 205 G4double yoffset,yMin,yMax; >> 206 G4double zoffset,zMin,zMax; >> 207 >> 208 zoffset=pTransform.NetTranslation().z(); >> 209 zMin=zoffset-fDz; >> 210 zMax=zoffset+fDz; >> 211 if (pVoxelLimit.IsZLimited()) >> 212 { >> 213 if ( (zMin>pVoxelLimit.GetMaxZExtent()+kCarTolerance) >> 214 || (zMax<pVoxelLimit.GetMinZExtent()-kCarTolerance) ) >> 215 { >> 216 return false; >> 217 } >> 218 else >> 219 { >> 220 if (zMin<pVoxelLimit.GetMinZExtent()) >> 221 { >> 222 zMin=pVoxelLimit.GetMinZExtent(); >> 223 } >> 224 if (zMax>pVoxelLimit.GetMaxZExtent()) >> 225 { >> 226 zMax=pVoxelLimit.GetMaxZExtent(); >> 227 } >> 228 } >> 229 } >> 230 xoffset=pTransform.NetTranslation().x(); >> 231 if (fDx2 >= fDx1) >> 232 { >> 233 xMax = xoffset+(fDx1+fDx2)/2+(zMax-zoffset)*(fDx2-fDx1)/(2*fDz) ; >> 234 xMin = 2*xoffset - xMax ; >> 235 } >> 236 else >> 237 { >> 238 xMax = xoffset+(fDx1+fDx2)/2+(zMin-zoffset)*(fDx2-fDx1)/(2*fDz) ; >> 239 xMin = 2*xoffset - xMax ; >> 240 } >> 241 if (pVoxelLimit.IsXLimited()) >> 242 { >> 243 if ( (xMin>pVoxelLimit.GetMaxXExtent()+kCarTolerance) >> 244 || (xMax<pVoxelLimit.GetMinXExtent()-kCarTolerance) ) >> 245 { >> 246 return false; >> 247 } >> 248 else >> 249 { >> 250 if (xMin<pVoxelLimit.GetMinXExtent()) >> 251 { >> 252 xMin=pVoxelLimit.GetMinXExtent(); >> 253 } >> 254 if (xMax>pVoxelLimit.GetMaxXExtent()) >> 255 { >> 256 xMax=pVoxelLimit.GetMaxXExtent(); >> 257 } >> 258 } >> 259 } >> 260 yoffset= pTransform.NetTranslation().y() ; >> 261 if(fDy2 >= fDy1) >> 262 { >> 263 yMax = yoffset+(fDy2+fDy1)/2+(zMax-zoffset)*(fDy2-fDy1)/(2*fDz) ; >> 264 yMin = 2*yoffset - yMax ; >> 265 } >> 266 else >> 267 { >> 268 yMax = yoffset+(fDy2+fDy1)/2+(zMin-zoffset)*(fDy2-fDy1)/(2*fDz) ; >> 269 yMin = 2*yoffset - yMax ; >> 270 } >> 271 if (pVoxelLimit.IsYLimited()) >> 272 { >> 273 if ( (yMin>pVoxelLimit.GetMaxYExtent()+kCarTolerance) >> 274 || (yMax<pVoxelLimit.GetMinYExtent()-kCarTolerance) ) >> 275 { >> 276 return false; >> 277 } >> 278 else >> 279 { >> 280 if (yMin<pVoxelLimit.GetMinYExtent()) >> 281 { >> 282 yMin=pVoxelLimit.GetMinYExtent(); >> 283 } >> 284 if (yMax>pVoxelLimit.GetMaxYExtent()) >> 285 { >> 286 yMax=pVoxelLimit.GetMaxYExtent(); >> 287 } >> 288 } >> 289 } >> 290 >> 291 switch (pAxis) >> 292 { >> 293 case kXAxis: >> 294 pMin=xMin; >> 295 pMax=xMax; >> 296 break; >> 297 case kYAxis: >> 298 pMin=yMin; >> 299 pMax=yMax; >> 300 break; >> 301 case kZAxis: >> 302 pMin=zMin; >> 303 pMax=zMax; >> 304 break; >> 305 default: >> 306 break; >> 307 } >> 308 >> 309 // Add 2*Tolerance to avoid precision troubles ? >> 310 // >> 311 pMin-=kCarTolerance; >> 312 pMax+=kCarTolerance; >> 313 >> 314 return true; 295 } 315 } >> 316 else >> 317 { >> 318 // General rotated case - create and clip mesh to boundaries 296 319 297 // Set bounding envelope (benv) and calculat << 320 G4bool existsAfterClip=false; 298 // << 321 G4ThreeVectorList *vertices; 299 G4double dx1 = GetXHalfLength1(); << 322 300 G4double dx2 = GetXHalfLength2(); << 323 pMin=+kInfinity; 301 G4double dy1 = GetYHalfLength1(); << 324 pMax=-kInfinity; 302 G4double dy2 = GetYHalfLength2(); << 325 303 G4double dz = GetZHalfLength(); << 326 // Calculate rotated vertex coordinates 304 << 327 // 305 G4ThreeVectorList baseA(4), baseB(4); << 328 vertices=CreateRotatedVertices(pTransform); 306 baseA[0].set(-dx1,-dy1,-dz); << 329 ClipCrossSection(vertices,0,pVoxelLimit,pAxis,pMin,pMax); 307 baseA[1].set( dx1,-dy1,-dz); << 330 ClipCrossSection(vertices,4,pVoxelLimit,pAxis,pMin,pMax); 308 baseA[2].set( dx1, dy1,-dz); << 331 ClipBetweenSections(vertices,0,pVoxelLimit,pAxis,pMin,pMax); 309 baseA[3].set(-dx1, dy1,-dz); << 332 310 baseB[0].set(-dx2,-dy2, dz); << 333 if (pMin!=kInfinity||pMax!=-kInfinity) 311 baseB[1].set( dx2,-dy2, dz); << 334 { 312 baseB[2].set( dx2, dy2, dz); << 335 existsAfterClip=true; 313 baseB[3].set(-dx2, dy2, dz); << 336 314 << 337 // Add 2*tolerance to avoid precision troubles 315 std::vector<const G4ThreeVectorList *> polyg << 338 // 316 polygons[0] = &baseA; << 339 pMin-=kCarTolerance; 317 polygons[1] = &baseB; << 340 pMax+=kCarTolerance; 318 << 341 319 G4BoundingEnvelope benv(bmin,bmax,polygons); << 342 } 320 exist = benv.CalculateExtent(pAxis,pVoxelLim << 343 else 321 return exist; << 344 { >> 345 // Check for case where completely enveloping clipping volume >> 346 // If point inside then we are confident that the solid completely >> 347 // envelopes the clipping volume. Hence set min/max extents according >> 348 // to clipping volume extents along the specified axis. >> 349 >> 350 G4ThreeVector clipCentre( >> 351 (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5, >> 352 (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5, >> 353 (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5); >> 354 >> 355 if (Inside(pTransform.Inverse().TransformPoint(clipCentre))!=kOutside) >> 356 { >> 357 existsAfterClip=true; >> 358 pMin=pVoxelLimit.GetMinExtent(pAxis); >> 359 pMax=pVoxelLimit.GetMaxExtent(pAxis); >> 360 } >> 361 } >> 362 delete vertices; >> 363 return existsAfterClip; >> 364 } 322 } 365 } 323 366 324 ////////////////////////////////////////////// << 367 /////////////////////////////////////////////////////////////////// 325 // 368 // 326 // Return whether point inside/outside/on surf 369 // Return whether point inside/outside/on surface, using tolerance 327 370 328 EInside G4Trd::Inside( const G4ThreeVector& p 371 EInside G4Trd::Inside( const G4ThreeVector& p ) const 329 { << 372 { 330 G4double dx = fPlanes[3].a*std::abs(p.x())+f << 373 EInside in=kOutside; 331 G4double dy = fPlanes[1].b*std::abs(p.y())+f << 374 G4double x,y,zbase1,zbase2; 332 G4double dxy = std::max(dx,dy); << 375 >> 376 if (std::fabs(p.z())<=fDz-kCarTolerance/2) >> 377 { >> 378 zbase1=p.z()+fDz; // Dist from -ve z plane >> 379 zbase2=fDz-p.z(); // Dist from +ve z plane 333 380 334 G4double dz = std::abs(p.z())-fDz; << 381 // Check whether inside x tolerance 335 G4double dist = std::max(dz,dxy); << 382 // >> 383 x=0.5*(fDx2*zbase1+fDx1*zbase2)/fDz - kCarTolerance/2; >> 384 if (std::fabs(p.x())<=x) >> 385 { >> 386 y=0.5*((fDy2*zbase1+fDy1*zbase2))/fDz - kCarTolerance/2; >> 387 if (std::fabs(p.y())<=y) >> 388 { >> 389 in=kInside; >> 390 } >> 391 else if (std::fabs(p.y())<=y+kCarTolerance) >> 392 { >> 393 in=kSurface; >> 394 } >> 395 } >> 396 else if (std::fabs(p.x())<=x+kCarTolerance) >> 397 { >> 398 // y = y half width of shape at z of point + tolerant boundary >> 399 // >> 400 y=0.5*((fDy2*zbase1+fDy1*zbase2))/fDz + kCarTolerance/2; >> 401 if (std::fabs(p.y())<=y) >> 402 { >> 403 in=kSurface; >> 404 } >> 405 } >> 406 } >> 407 else if (std::fabs(p.z())<=fDz+kCarTolerance/2) >> 408 { >> 409 // Only need to check outer tolerant boundaries >> 410 // >> 411 zbase1=p.z()+fDz; // Dist from -ve z plane >> 412 zbase2=fDz-p.z(); // Dist from +ve z plane 336 413 337 return (dist > halfCarTolerance) ? kOutside << 414 // x = x half width of shape at z of point plus tolerance 338 ((dist > -halfCarTolerance) ? kSurface : k << 415 // >> 416 x=0.5*(fDx2*zbase1+fDx1*zbase2)/fDz + kCarTolerance/2; >> 417 if (std::fabs(p.x())<=x) >> 418 { >> 419 // y = y half width of shape at z of point >> 420 // >> 421 y=0.5*((fDy2*zbase1+fDy1*zbase2))/fDz + kCarTolerance/2; >> 422 if (std::fabs(p.y())<=y) in=kSurface; >> 423 } >> 424 } >> 425 return in; 339 } 426 } 340 427 341 ////////////////////////////////////////////// 428 ////////////////////////////////////////////////////////////////////////// 342 // 429 // 343 // Determine side where point is, and return c << 430 // Calculate side nearest to p, and return normal >> 431 // If two sides are equidistant, normal of first side (x/y/z) >> 432 // encountered returned 344 433 345 G4ThreeVector G4Trd::SurfaceNormal( const G4Th 434 G4ThreeVector G4Trd::SurfaceNormal( const G4ThreeVector& p ) const 346 { 435 { 347 G4int nsurf = 0; // number of surfaces where << 436 G4ThreeVector norm, sumnorm(0.,0.,0.); 348 << 437 G4int noSurfaces = 0; 349 // Check Z faces << 438 G4double z = 2.0*fDz, tanx, secx, newpx, widx; 350 // << 439 G4double tany, secy, newpy, widy; 351 G4double nz = 0; << 440 G4double distx, disty, distz, fcos; 352 G4double dz = std::abs(p.z()) - fDz; << 441 G4double delta = 0.5*kCarTolerance; 353 if (std::abs(dz) <= halfCarTolerance) << 442 354 { << 443 tanx = (fDx2 - fDx1)/z; 355 nz = (p.z() < 0) ? -1 : 1; << 444 secx = std::sqrt(1.0+tanx*tanx); 356 ++nsurf; << 445 newpx = std::fabs(p.x())-p.z()*tanx; >> 446 widx = fDx2 - fDz*tanx; >> 447 >> 448 tany = (fDy2 - fDy1)/z; >> 449 secy = std::sqrt(1.0+tany*tany); >> 450 newpy = std::fabs(p.y())-p.z()*tany; >> 451 widy = fDy2 - fDz*tany; >> 452 >> 453 distx = std::fabs(newpx-widx)/secx; // perp. distance to x side >> 454 disty = std::fabs(newpy-widy)/secy; // to y side >> 455 distz = std::fabs(std::fabs(p.z())-fDz); // to z side >> 456 >> 457 fcos = 1.0/secx; >> 458 G4ThreeVector nX = G4ThreeVector( fcos,0,-tanx*fcos); >> 459 G4ThreeVector nmX = G4ThreeVector(-fcos,0,-tanx*fcos); >> 460 >> 461 fcos = 1.0/secy; >> 462 G4ThreeVector nY = G4ThreeVector(0, fcos,-tany*fcos); >> 463 G4ThreeVector nmY = G4ThreeVector(0,-fcos,-tany*fcos); >> 464 G4ThreeVector nZ = G4ThreeVector( 0, 0, 1.0); >> 465 >> 466 if (distx <= delta) >> 467 { >> 468 noSurfaces ++; >> 469 if ( p.x() >= 0.) sumnorm += nX; >> 470 else sumnorm += nmX; >> 471 } >> 472 if (disty <= delta) >> 473 { >> 474 noSurfaces ++; >> 475 if ( p.y() >= 0.) sumnorm += nY; >> 476 else sumnorm += nmY; >> 477 } >> 478 if (distz <= delta) >> 479 { >> 480 noSurfaces ++; >> 481 if ( p.z() >= 0.) sumnorm += nZ; >> 482 else sumnorm -= nZ; 357 } 483 } 358 << 484 if ( noSurfaces == 0 ) 359 // Check Y faces << 360 // << 361 G4double ny = 0; << 362 G4double dy1 = fPlanes[0].b*p.y(); << 363 G4double dy2 = fPlanes[0].c*p.z() + fPlanes[ << 364 if (std::abs(dy2 + dy1) <= halfCarTolerance) << 365 { 485 { 366 ny += fPlanes[0].b; << 367 nz += fPlanes[0].c; << 368 ++nsurf; << 369 } << 370 if (std::abs(dy2 - dy1) <= halfCarTolerance) << 371 { << 372 ny += fPlanes[1].b; << 373 nz += fPlanes[1].c; << 374 ++nsurf; << 375 } << 376 << 377 // Check X faces << 378 // << 379 G4double nx = 0; << 380 G4double dx1 = fPlanes[2].a*p.x(); << 381 G4double dx2 = fPlanes[2].c*p.z() + fPlanes[ << 382 if (std::abs(dx2 + dx1) <= halfCarTolerance) << 383 { << 384 nx += fPlanes[2].a; << 385 nz += fPlanes[2].c; << 386 ++nsurf; << 387 } << 388 if (std::abs(dx2 - dx1) <= halfCarTolerance) << 389 { << 390 nx += fPlanes[3].a; << 391 nz += fPlanes[3].c; << 392 ++nsurf; << 393 } << 394 << 395 // Return normal << 396 // << 397 if (nsurf == 1) return {nx,ny,nz}; << 398 else if (nsurf != 0) return G4ThreeVector(nx << 399 else << 400 { << 401 // Point is not on the surface << 402 // << 403 #ifdef G4CSGDEBUG 486 #ifdef G4CSGDEBUG 404 std::ostringstream message; << 487 G4Exception("G4Trd::SurfaceNormal(p)", "GeomSolids1002", JustWarning, 405 G4long oldprc = message.precision(16); << 488 "Point p is not on surface !?" ); 406 message << "Point p is not on surface (!?) << 489 #endif 407 << GetName() << G4endl; << 490 norm = ApproxSurfaceNormal(p); 408 message << "Position:\n"; << 409 message << " p.x() = " << p.x()/mm << " << 410 message << " p.y() = " << p.y()/mm << " << 411 message << " p.z() = " << p.z()/mm << " << 412 G4cout.precision(oldprc) ; << 413 G4Exception("G4Trd::SurfaceNormal(p)", "Ge << 414 JustWarning, message ); << 415 DumpInfo(); << 416 #endif << 417 return ApproxSurfaceNormal(p); << 418 } 491 } >> 492 else if ( noSurfaces == 1 ) norm = sumnorm; >> 493 else norm = sumnorm.unit(); >> 494 return norm; 419 } 495 } 420 496 421 ////////////////////////////////////////////// << 497 >> 498 ///////////////////////////////////////////////////////////////////////////// 422 // 499 // 423 // Algorithm for SurfaceNormal() following the 500 // Algorithm for SurfaceNormal() following the original specification 424 // for points not on the surface 501 // for points not on the surface 425 502 426 G4ThreeVector G4Trd::ApproxSurfaceNormal( cons 503 G4ThreeVector G4Trd::ApproxSurfaceNormal( const G4ThreeVector& p ) const 427 { 504 { 428 G4double dist = -DBL_MAX; << 505 G4ThreeVector norm; 429 G4int iside = 0; << 506 G4double z,tanx,secx,newpx,widx; 430 for (G4int i=0; i<4; ++i) << 507 G4double tany,secy,newpy,widy; 431 { << 508 G4double distx,disty,distz,fcos; 432 G4double d = fPlanes[i].a*p.x() + << 509 433 fPlanes[i].b*p.y() + << 510 z=2.0*fDz; 434 fPlanes[i].c*p.z() + fPlanes[ << 511 435 if (d > dist) { dist = d; iside = i; } << 512 tanx=(fDx2-fDx1)/z; 436 } << 513 secx=std::sqrt(1.0+tanx*tanx); >> 514 newpx=std::fabs(p.x())-p.z()*tanx; >> 515 widx=fDx2-fDz*tanx; >> 516 >> 517 tany=(fDy2-fDy1)/z; >> 518 secy=std::sqrt(1.0+tany*tany); >> 519 newpy=std::fabs(p.y())-p.z()*tany; >> 520 widy=fDy2-fDz*tany; >> 521 >> 522 distx=std::fabs(newpx-widx)/secx; // perpendicular distance to x side >> 523 disty=std::fabs(newpy-widy)/secy; // to y side >> 524 distz=std::fabs(std::fabs(p.z())-fDz); // to z side 437 525 438 G4double distz = std::abs(p.z()) - fDz; << 526 // find closest side 439 if (dist > distz) << 527 // 440 return { fPlanes[iside].a, fPlanes[iside]. << 528 if (distx<=disty) >> 529 { >> 530 if (distx<=distz) >> 531 { >> 532 // Closest to X >> 533 // >> 534 fcos=1.0/secx; >> 535 // normal=(+/-std::cos(ang),0,-std::sin(ang)) >> 536 if (p.x()>=0) >> 537 norm=G4ThreeVector(fcos,0,-tanx*fcos); >> 538 else >> 539 norm=G4ThreeVector(-fcos,0,-tanx*fcos); >> 540 } >> 541 else >> 542 { >> 543 // Closest to Z >> 544 // >> 545 if (p.z()>=0) >> 546 norm=G4ThreeVector(0,0,1); >> 547 else >> 548 norm=G4ThreeVector(0,0,-1); >> 549 } >> 550 } 441 else 551 else 442 return { 0, 0, (G4double)((p.z() < 0) ? -1 << 552 { >> 553 if (disty<=distz) >> 554 { >> 555 // Closest to Y >> 556 // >> 557 fcos=1.0/secy; >> 558 if (p.y()>=0) >> 559 norm=G4ThreeVector(0,fcos,-tany*fcos); >> 560 else >> 561 norm=G4ThreeVector(0,-fcos,-tany*fcos); >> 562 } >> 563 else >> 564 { >> 565 // Closest to Z >> 566 // >> 567 if (p.z()>=0) >> 568 norm=G4ThreeVector(0,0,1); >> 569 else >> 570 norm=G4ThreeVector(0,0,-1); >> 571 } >> 572 } >> 573 return norm; 443 } 574 } 444 575 445 ////////////////////////////////////////////// << 576 //////////////////////////////////////////////////////////////////////////// 446 // 577 // 447 // Calculate distance to shape from outside 578 // Calculate distance to shape from outside 448 // - return kInfinity if no intersection << 579 // - return kInfinity if no intersection >> 580 // >> 581 // ALGORITHM: >> 582 // For each component, calculate pair of minimum and maximum intersection >> 583 // values for which the particle is in the extent of the shape >> 584 // - The smallest (MAX minimum) allowed distance of the pairs is intersect >> 585 // - Z plane intersectin uses tolerance >> 586 // - XZ YZ planes use logic & *SLIGHTLY INCORRECT* tolerance >> 587 // (this saves at least 1 sqrt, 1 multiply and 1 divide... in applicable >> 588 // cases) >> 589 // - Note: XZ and YZ planes each divide space into four regions, >> 590 // characterised by ss1 ss2 >> 591 // NOTE: >> 592 // >> 593 // `Inside' safe - meaningful answers given if point is inside the exact >> 594 // shape. >> 595 >> 596 G4double G4Trd::DistanceToIn( const G4ThreeVector& p, >> 597 const G4ThreeVector& v ) const >> 598 { >> 599 G4double snxt = kInfinity ; // snxt = default return value >> 600 G4double smin,smax; >> 601 G4double s1,s2,tanxz,tanyz,ds1,ds2; >> 602 G4double ss1,ss2,sn1=0.,sn2=0.,Dist; 449 603 450 G4double G4Trd::DistanceToIn(const G4ThreeVect << 604 if ( v.z() ) // Calculate valid z intersect range 451 const G4ThreeVect << 605 { 452 { << 606 if ( v.z() > 0 ) // Calculate smax: must be +ve or no intersection. 453 // Z intersections << 607 { 454 // << 608 Dist = fDz - p.z() ; // to plane at +dz 455 if ((std::abs(p.z()) - fDz) >= -halfCarToler << 456 return kInfinity; << 457 G4double invz = (-v.z() == 0) ? DBL_MAX : -1 << 458 G4double dz = (invz < 0) ? fDz : -fDz; << 459 G4double tzmin = (p.z() + dz)*invz; << 460 G4double tzmax = (p.z() - dz)*invz; << 461 609 462 // Y intersections << 610 if (Dist >= 0.5*kCarTolerance) 463 // << 611 { 464 G4double tmin0 = tzmin, tmax0 = tzmax; << 612 smax = Dist/v.z() ; 465 G4double ya = fPlanes[0].b*v.y(), yb = fPlan << 613 smin = -(fDz + p.z())/v.z() ; 466 G4double yc = fPlanes[0].b*p.y(), yd = fPlan << 614 } 467 G4double cos0 = yb + ya; << 615 else return snxt ; 468 G4double dis0 = yd + yc; << 616 } 469 if (dis0 >= -halfCarTolerance) << 617 else // v.z <0 >> 618 { >> 619 Dist=fDz+p.z(); // plane at -dz >> 620 >> 621 if ( Dist >= 0.5*kCarTolerance ) >> 622 { >> 623 smax = -Dist/v.z() ; >> 624 smin = (fDz - p.z())/v.z() ; >> 625 } >> 626 else return snxt ; >> 627 } >> 628 if (smin < 0 ) smin = 0 ; >> 629 } >> 630 else // v.z=0 470 { 631 { 471 if (cos0 >= 0) return kInfinity; << 632 if (std::fabs(p.z()) >= fDz ) return snxt ; // Outside & no intersect 472 G4double tmp = -dis0/cos0; << 633 else 473 if (tmin0 < tmp) tmin0 = tmp; << 634 { >> 635 smin = 0 ; // Always inside z range >> 636 smax = kInfinity; >> 637 } 474 } 638 } 475 else if (cos0 > 0) << 639 >> 640 // Calculate x intersection range >> 641 // >> 642 // Calc half width at p.z, and components towards planes >> 643 >> 644 tanxz = (fDx2 - fDx1)*0.5/fDz ; >> 645 s1 = 0.5*(fDx1+fDx2) + tanxz*p.z() ; // x half width at p.z >> 646 ds1 = v.x() - tanxz*v.z() ; // Components of v towards faces at +-x >> 647 ds2 = v.x() + tanxz*v.z() ; >> 648 ss1 = s1 - p.x() ; // -delta x to +ve plane >> 649 // -ve when outside >> 650 ss2 = -s1 - p.x() ; // -delta x to -ve plane >> 651 // +ve when outside >> 652 >> 653 if (ss1 < 0 && ss2 <= 0 ) 476 { 654 { 477 G4double tmp = -dis0/cos0; << 655 if (ds1 < 0) // In +ve coord Area 478 if (tmax0 > tmp) tmax0 = tmp; << 656 { >> 657 sn1 = ss1/ds1 ; >> 658 >> 659 if ( ds2 < 0 ) sn2 = ss2/ds2 ; >> 660 else sn2 = kInfinity ; >> 661 } >> 662 else return snxt ; 479 } 663 } >> 664 else if ( ss1 >= 0 && ss2 > 0 ) >> 665 { >> 666 if ( ds2 > 0 ) // In -ve coord Area >> 667 { >> 668 sn1 = ss2/ds2 ; 480 669 481 G4double tmin1 = tmin0, tmax1 = tmax0; << 670 if (ds1 > 0) sn2 = ss1/ds1 ; 482 G4double cos1 = yb - ya; << 671 else sn2 = kInfinity; 483 G4double dis1 = yd - yc; << 672 484 if (dis1 >= -halfCarTolerance) << 673 } >> 674 else return snxt ; >> 675 } >> 676 else if (ss1 >= 0 && ss2 <= 0 ) 485 { 677 { 486 if (cos1 >= 0) return kInfinity; << 678 // Inside Area - calculate leaving distance 487 G4double tmp = -dis1/cos1; << 679 // *Don't* use exact distance to side for tolerance 488 if (tmin1 < tmp) tmin1 = tmp; << 680 // = ss1*std::cos(ang xz) >> 681 // = ss1/std::sqrt(1.0+tanxz*tanxz) >> 682 sn1 = 0 ; >> 683 >> 684 if ( ds1 > 0 ) >> 685 { >> 686 if (ss1 > 0.5*kCarTolerance) sn2 = ss1/ds1 ; // Leave +ve side extent >> 687 else return snxt ; // Leave immediately by +ve >> 688 } >> 689 else sn2 = kInfinity ; >> 690 >> 691 if ( ds2 < 0 ) >> 692 { >> 693 if ( ss2 < -0.5*kCarTolerance ) >> 694 { >> 695 Dist = ss2/ds2 ; // Leave -ve side extent >> 696 if ( Dist < sn2 ) sn2 = Dist ; >> 697 } >> 698 else return snxt ; >> 699 } 489 } 700 } 490 else if (cos1 > 0) << 701 else if (ss1 < 0 && ss2 > 0 ) 491 { 702 { 492 G4double tmp = -dis1/cos1; << 703 // Within +/- plane cross-over areas (not on boundaries ss1||ss2==0) 493 if (tmax1 > tmp) tmax1 = tmp; << 704 >> 705 if ( ds1 >= 0 || ds2 <= 0 ) >> 706 { >> 707 return snxt ; >> 708 } >> 709 else // Will intersect & stay inside >> 710 { >> 711 sn1 = ss1/ds1 ; >> 712 Dist = ss2/ds2 ; >> 713 if (Dist > sn1 ) sn1 = Dist ; >> 714 sn2 = kInfinity ; >> 715 } 494 } 716 } 495 717 496 // X intersections << 718 // Reduce allowed range of distances as appropriate 497 // << 719 498 G4double tmin2 = tmin1, tmax2 = tmax1; << 720 if ( sn1 > smin ) smin = sn1 ; 499 G4double xa = fPlanes[2].a*v.x(), xb = fPlan << 721 if ( sn2 < smax ) smax = sn2 ; 500 G4double xc = fPlanes[2].a*p.x(), xd = fPlan << 722 501 G4double cos2 = xb + xa; << 723 // Check for incompatible ranges (eg z intersects between 50 ->100 and x 502 G4double dis2 = xd + xc; << 724 // only 10-40 -> no intersection) 503 if (dis2 >= -halfCarTolerance) << 725 >> 726 if ( smax < smin ) return snxt ; >> 727 >> 728 // Calculate valid y intersection range >> 729 // (repeat of x intersection code) >> 730 >> 731 tanyz = (fDy2-fDy1)*0.5/fDz ; >> 732 s2 = 0.5*(fDy1+fDy2) + tanyz*p.z() ; // y half width at p.z >> 733 ds1 = v.y() - tanyz*v.z() ; // Components of v towards faces at +-y >> 734 ds2 = v.y() + tanyz*v.z() ; >> 735 ss1 = s2 - p.y() ; // -delta y to +ve plane >> 736 ss2 = -s2 - p.y() ; // -delta y to -ve plane >> 737 >> 738 if ( ss1 < 0 && ss2 <= 0 ) 504 { 739 { 505 if (cos2 >= 0) return kInfinity; << 740 if (ds1 < 0 ) // In +ve coord Area 506 G4double tmp = -dis2/cos2; << 741 { 507 if (tmin2 < tmp) tmin2 = tmp; << 742 sn1 = ss1/ds1 ; >> 743 if ( ds2 < 0 ) sn2 = ss2/ds2 ; >> 744 else sn2 = kInfinity ; >> 745 } >> 746 else return snxt ; 508 } 747 } 509 else if (cos2 > 0) << 748 else if ( ss1 >= 0 && ss2 > 0 ) 510 { 749 { 511 G4double tmp = -dis2/cos2; << 750 if ( ds2 > 0 ) // In -ve coord Area 512 if (tmax2 > tmp) tmax2 = tmp; << 751 { >> 752 sn1 = ss2/ds2 ; >> 753 if ( ds1 > 0 ) sn2 = ss1/ds1 ; >> 754 else sn2 = kInfinity ; >> 755 } >> 756 else return snxt ; 513 } 757 } 514 << 758 else if (ss1 >= 0 && ss2 <= 0 ) 515 G4double tmin3 = tmin2, tmax3 = tmax2; << 516 G4double cos3 = xb - xa; << 517 G4double dis3 = xd - xc; << 518 if (dis3 >= -halfCarTolerance) << 519 { 759 { 520 if (cos3 >= 0) return kInfinity; << 760 // Inside Area - calculate leaving distance 521 G4double tmp = -dis3/cos3; << 761 // *Don't* use exact distance to side for tolerance 522 if (tmin3 < tmp) tmin3 = tmp; << 762 // = ss1*std::cos(ang yz) >> 763 // = ss1/std::sqrt(1.0+tanyz*tanyz) >> 764 sn1 = 0 ; >> 765 >> 766 if ( ds1 > 0 ) >> 767 { >> 768 if (ss1 > 0.5*kCarTolerance) sn2 = ss1/ds1 ; // Leave +ve side extent >> 769 else return snxt ; // Leave immediately by +ve >> 770 } >> 771 else sn2 = kInfinity ; >> 772 >> 773 if ( ds2 < 0 ) >> 774 { >> 775 if ( ss2 < -0.5*kCarTolerance ) >> 776 { >> 777 Dist = ss2/ds2 ; // Leave -ve side extent >> 778 if (Dist < sn2) sn2=Dist; >> 779 } >> 780 else return snxt ; >> 781 } 523 } 782 } 524 else if (cos3 > 0) << 783 else if (ss1 < 0 && ss2 > 0 ) 525 { 784 { 526 G4double tmp = -dis3/cos3; << 785 // Within +/- plane cross-over areas (not on boundaries ss1||ss2==0) 527 if (tmax3 > tmp) tmax3 = tmp; << 786 >> 787 if (ds1 >= 0 || ds2 <= 0 ) >> 788 { >> 789 return snxt ; >> 790 } >> 791 else // Will intersect & stay inside >> 792 { >> 793 sn1 = ss1/ds1 ; >> 794 Dist = ss2/ds2 ; >> 795 if (Dist > sn1 ) sn1 = Dist ; >> 796 sn2 = kInfinity ; >> 797 } 528 } 798 } >> 799 >> 800 // Reduce allowed range of distances as appropriate 529 801 530 // Find distance << 802 if ( sn1 > smin) smin = sn1 ; 531 // << 803 if ( sn2 < smax) smax = sn2 ; 532 G4double tmin = tmin3, tmax = tmax3; << 804 533 if (tmax <= tmin + halfCarTolerance) return << 805 // Check for incompatible ranges (eg x intersects between 50 ->100 and y 534 return (tmin < halfCarTolerance ) ? 0. : tmi << 806 // only 10-40 -> no intersection). Set snxt if ok >> 807 >> 808 if ( smax > smin ) snxt = smin ; >> 809 if (snxt < 0.5*kCarTolerance ) snxt = 0.0 ; >> 810 >> 811 return snxt ; 535 } 812 } 536 813 537 ////////////////////////////////////////////// << 814 ///////////////////////////////////////////////////////////////////////// 538 // 815 // 539 // Calculate exact shortest distance to any bo << 816 // Approximate distance to shape 540 // This is the best fast estimation of the sho << 817 // Calculate perpendicular distances to z/x/y surfaces, return largest 541 // - returns 0 if point is inside << 818 // which is the most fast estimation of shortest distance to Trd >> 819 // - Safe underestimate >> 820 // - If point within exact shape, return 0 542 821 543 G4double G4Trd::DistanceToIn( const G4ThreeVec 822 G4double G4Trd::DistanceToIn( const G4ThreeVector& p ) const 544 { 823 { 545 G4double dx = fPlanes[3].a*std::abs(p.x())+f << 824 G4double safe=0.0; 546 G4double dy = fPlanes[1].b*std::abs(p.y())+f << 825 G4double tanxz,distx,safx; 547 G4double dxy = std::max(dx,dy); << 826 G4double tanyz,disty,safy; >> 827 G4double zbase; >> 828 >> 829 safe=std::fabs(p.z())-fDz; >> 830 if (safe<0) safe=0; // Also used to ensure x/y distances >> 831 // POSITIVE 548 832 549 G4double dz = std::abs(p.z())-fDz; << 833 zbase=fDz+p.z(); 550 G4double dist = std::max(dz,dxy); << 551 834 552 return (dist > 0) ? dist : 0.; << 835 // Find distance along x direction to closest x plane >> 836 // >> 837 tanxz=(fDx2-fDx1)*0.5/fDz; >> 838 // widx=fDx1+tanxz*(fDz+p.z()); // x width at p.z >> 839 // distx=std::fabs(p.x())-widx; // distance to plane >> 840 distx=std::fabs(p.x())-(fDx1+tanxz*zbase); >> 841 if (distx>safe) >> 842 { >> 843 safx=distx/std::sqrt(1.0+tanxz*tanxz); // vector Dist=Dist*std::cos(ang) >> 844 if (safx>safe) safe=safx; >> 845 } >> 846 >> 847 // Find distance along y direction to slanted wall >> 848 tanyz=(fDy2-fDy1)*0.5/fDz; >> 849 // widy=fDy1+tanyz*(fDz+p.z()); // y width at p.z >> 850 // disty=std::fabs(p.y())-widy; // distance to plane >> 851 disty=std::fabs(p.y())-(fDy1+tanyz*zbase); >> 852 if (disty>safe) >> 853 { >> 854 safy=disty/std::sqrt(1.0+tanyz*tanyz); // distance along vector >> 855 if (safy>safe) safe=safy; >> 856 } >> 857 return safe; 553 } 858 } 554 859 555 ////////////////////////////////////////////// << 860 //////////////////////////////////////////////////////////////////////// 556 // 861 // 557 // Calculate distance to surface of shape from << 862 // Calcluate distance to surface of shape from inside 558 // find normal at exit point, if required << 863 // Calculate distance to x/y/z planes - smallest is exiting distance 559 // - when leaving the surface, return 0 << 864 // - z planes have std. check for tolerance 560 << 865 // - xz yz planes have check based on distance || to x or y axis 561 G4double G4Trd::DistanceToOut(const G4ThreeVec << 866 // (not corrected for slope of planes) 562 const G4bool cal << 867 // ?BUG? If v.z==0 are there cases when snside not set???? 563 G4bool* va << 868 >> 869 G4double G4Trd::DistanceToOut( const G4ThreeVector& p, >> 870 const G4ThreeVector& v, >> 871 const G4bool calcNorm, >> 872 G4bool *validNorm, >> 873 G4ThreeVector *n ) const 564 { 874 { 565 // Z intersections << 875 ESide side = kUndefined, snside = kUndefined; 566 // << 876 G4double snxt,pdist; 567 if ((std::abs(p.z()) - fDz) >= -halfCarToler << 877 G4double central,ss1,ss2,ds1,ds2,sn=0.,sn2=0.; >> 878 G4double tanxz=0.,cosxz=0.,tanyz=0.,cosyz=0.; >> 879 >> 880 if (calcNorm) *validNorm=true; // All normals are valid >> 881 >> 882 // Calculate z plane intersection >> 883 if (v.z()>0) 568 { 884 { 569 if (calcNorm) << 885 pdist=fDz-p.z(); >> 886 if (pdist>kCarTolerance/2) 570 { 887 { 571 *validNorm = true; << 888 snxt=pdist/v.z(); 572 n->set(0, 0, (p.z() < 0) ? -1 : 1); << 889 side=kPZ; >> 890 } >> 891 else >> 892 { >> 893 if (calcNorm) >> 894 { >> 895 *n=G4ThreeVector(0,0,1); >> 896 } >> 897 return snxt=0; >> 898 } >> 899 } >> 900 else if (v.z()<0) >> 901 { >> 902 pdist=fDz+p.z(); >> 903 if (pdist>kCarTolerance/2) >> 904 { >> 905 snxt=-pdist/v.z(); >> 906 side=kMZ; >> 907 } >> 908 else >> 909 { >> 910 if (calcNorm) >> 911 { >> 912 *n=G4ThreeVector(0,0,-1); >> 913 } >> 914 return snxt=0; 573 } 915 } 574 return 0; << 575 } 916 } 576 G4double vz = v.z(); << 917 else 577 G4double tmax = (vz == 0) ? DBL_MAX : (std:: << 918 { 578 G4int iside = (vz < 0) ? -4 : -2; // little << 919 snxt=kInfinity; >> 920 } 579 921 580 // Y intersections << 581 // 922 // 582 G4int i = 0; << 923 // Calculate x intersection 583 for ( ; i<2; ++i) << 924 // 584 { << 925 tanxz=(fDx2-fDx1)*0.5/fDz; 585 G4double cosa = fPlanes[i].b*v.y() + fPlan << 926 central=0.5*(fDx1+fDx2); 586 if (cosa > 0) << 927 >> 928 // +ve plane (1) >> 929 // >> 930 ss1=central+tanxz*p.z()-p.x(); // distance || x axis to plane >> 931 // (+ve if point inside) >> 932 ds1=v.x()-tanxz*v.z(); // component towards plane at +x >> 933 // (-ve if +ve -> -ve direction) >> 934 // -ve plane (2) >> 935 // >> 936 ss2=-tanxz*p.z()-p.x()-central; //distance || x axis to plane >> 937 // (-ve if point inside) >> 938 ds2=tanxz*v.z()+v.x(); // component towards plane at -x >> 939 >> 940 if (ss1>0&&ss2<0) >> 941 { >> 942 // Normal case - entirely inside region >> 943 if (ds1<=0&&ds2<0) >> 944 { >> 945 if (ss2<-kCarTolerance/2) >> 946 { >> 947 sn=ss2/ds2; // Leave by -ve side >> 948 snside=kMX; >> 949 } >> 950 else >> 951 { >> 952 sn=0; // Leave immediately by -ve side >> 953 snside=kMX; >> 954 } >> 955 } >> 956 else if (ds1>0&&ds2>=0) 587 { 957 { 588 G4double dist = fPlanes[i].b*p.y()+fPlan << 958 if (ss1>kCarTolerance/2) 589 if (dist >= -halfCarTolerance) << 590 { 959 { 591 if (calcNorm) << 960 sn=ss1/ds1; // Leave by +ve side >> 961 snside=kPX; >> 962 } >> 963 else >> 964 { >> 965 sn=0; // Leave immediately by +ve side >> 966 snside=kPX; >> 967 } >> 968 } >> 969 else if (ds1>0&&ds2<0) >> 970 { >> 971 if (ss1>kCarTolerance/2) >> 972 { >> 973 // sn=ss1/ds1; // Leave by +ve side >> 974 if (ss2<-kCarTolerance/2) 592 { 975 { 593 *validNorm = true; << 976 sn=ss1/ds1; // Leave by +ve side 594 n->set(0, fPlanes[i].b, fPlanes[i].c << 977 sn2=ss2/ds2; >> 978 if (sn2<sn) >> 979 { >> 980 sn=sn2; >> 981 snside=kMX; >> 982 } >> 983 else >> 984 { >> 985 snside=kPX; >> 986 } 595 } 987 } 596 return 0; << 988 else >> 989 { >> 990 sn=0; // Leave immediately by -ve >> 991 snside=kMX; >> 992 } >> 993 } >> 994 else >> 995 { >> 996 sn=0; // Leave immediately by +ve side >> 997 snside=kPX; >> 998 } >> 999 } >> 1000 else >> 1001 { >> 1002 // Must be || to both >> 1003 // >> 1004 sn=kInfinity; // Don't leave by either side >> 1005 } >> 1006 } >> 1007 else if (ss1<=0&&ss2<0) >> 1008 { >> 1009 // Outside, in +ve Area >> 1010 >> 1011 if (ds1>0) >> 1012 { >> 1013 sn=0; // Away from shape >> 1014 // Left by +ve side >> 1015 snside=kPX; >> 1016 } >> 1017 else >> 1018 { >> 1019 if (ds2<0) >> 1020 { >> 1021 // Ignore +ve plane and use -ve plane intersect >> 1022 // >> 1023 sn=ss2/ds2; // Leave by -ve side >> 1024 snside=kMX; >> 1025 } >> 1026 else >> 1027 { >> 1028 // Must be || to both -> exit determined by other axes >> 1029 // >> 1030 sn=kInfinity; // Don't leave by either side 597 } 1031 } 598 G4double tmp = -dist/cosa; << 599 if (tmax > tmp) { tmax = tmp; iside = i; << 600 } 1032 } 601 } 1033 } >> 1034 else if (ss1>0&&ss2>=0) >> 1035 { >> 1036 // Outside, in -ve Area 602 1037 603 // X intersections << 1038 if (ds2<0) 604 // << 1039 { 605 for ( ; i<4; ++i) << 1040 sn=0; // away from shape >> 1041 // Left by -ve side >> 1042 snside=kMX; >> 1043 } >> 1044 else >> 1045 { >> 1046 if (ds1>0) >> 1047 { >> 1048 // Ignore +ve plane and use -ve plane intersect >> 1049 // >> 1050 sn=ss1/ds1; // Leave by +ve side >> 1051 snside=kPX; >> 1052 } >> 1053 else >> 1054 { >> 1055 // Must be || to both -> exit determined by other axes >> 1056 // >> 1057 sn=kInfinity; // Don't leave by either side >> 1058 } >> 1059 } >> 1060 } >> 1061 >> 1062 // Update minimum exit distance >> 1063 >> 1064 if (sn<snxt) 606 { 1065 { 607 G4double cosa = fPlanes[i].a*v.x()+fPlanes << 1066 snxt=sn; 608 if (cosa > 0) << 1067 side=snside; >> 1068 } >> 1069 if (snxt>0) >> 1070 { >> 1071 // Calculate y intersection >> 1072 >> 1073 tanyz=(fDy2-fDy1)*0.5/fDz; >> 1074 central=0.5*(fDy1+fDy2); >> 1075 >> 1076 // +ve plane (1) >> 1077 // >> 1078 ss1=central+tanyz*p.z()-p.y(); // distance || y axis to plane >> 1079 // (+ve if point inside) >> 1080 ds1=v.y()-tanyz*v.z(); // component towards +ve plane >> 1081 // (-ve if +ve -> -ve direction) >> 1082 // -ve plane (2) >> 1083 // >> 1084 ss2=-tanyz*p.z()-p.y()-central; // distance || y axis to plane >> 1085 // (-ve if point inside) >> 1086 ds2=tanyz*v.z()+v.y(); // component towards -ve plane >> 1087 >> 1088 if (ss1>0&&ss2<0) >> 1089 { >> 1090 // Normal case - entirely inside region >> 1091 >> 1092 if (ds1<=0&&ds2<0) >> 1093 { >> 1094 if (ss2<-kCarTolerance/2) >> 1095 { >> 1096 sn=ss2/ds2; // Leave by -ve side >> 1097 snside=kMY; >> 1098 } >> 1099 else >> 1100 { >> 1101 sn=0; // Leave immediately by -ve side >> 1102 snside=kMY; >> 1103 } >> 1104 } >> 1105 else if (ds1>0&&ds2>=0) >> 1106 { >> 1107 if (ss1>kCarTolerance/2) >> 1108 { >> 1109 sn=ss1/ds1; // Leave by +ve side >> 1110 snside=kPY; >> 1111 } >> 1112 else >> 1113 { >> 1114 sn=0; // Leave immediately by +ve side >> 1115 snside=kPY; >> 1116 } >> 1117 } >> 1118 else if (ds1>0&&ds2<0) >> 1119 { >> 1120 if (ss1>kCarTolerance/2) >> 1121 { >> 1122 // sn=ss1/ds1; // Leave by +ve side >> 1123 if (ss2<-kCarTolerance/2) >> 1124 { >> 1125 sn=ss1/ds1; // Leave by +ve side >> 1126 sn2=ss2/ds2; >> 1127 if (sn2<sn) >> 1128 { >> 1129 sn=sn2; >> 1130 snside=kMY; >> 1131 } >> 1132 else >> 1133 { >> 1134 snside=kPY; >> 1135 } >> 1136 } >> 1137 else >> 1138 { >> 1139 sn=0; // Leave immediately by -ve >> 1140 snside=kMY; >> 1141 } >> 1142 } >> 1143 else >> 1144 { >> 1145 sn=0; // Leave immediately by +ve side >> 1146 snside=kPY; >> 1147 } >> 1148 } >> 1149 else >> 1150 { >> 1151 // Must be || to both >> 1152 // >> 1153 sn=kInfinity; // Don't leave by either side >> 1154 } >> 1155 } >> 1156 else if (ss1<=0&&ss2<0) 609 { 1157 { 610 G4double dist = fPlanes[i].a*p.x()+fPlan << 1158 // Outside, in +ve Area 611 if (dist >= -halfCarTolerance) << 1159 >> 1160 if (ds1>0) 612 { 1161 { 613 if (calcNorm) << 1162 sn=0; // Away from shape >> 1163 // Left by +ve side >> 1164 snside=kPY; >> 1165 } >> 1166 else >> 1167 { >> 1168 if (ds2<0) >> 1169 { >> 1170 // Ignore +ve plane and use -ve plane intersect >> 1171 // >> 1172 sn=ss2/ds2; // Leave by -ve side >> 1173 snside=kMY; >> 1174 } >> 1175 else 614 { 1176 { 615 *validNorm = true; << 1177 // Must be || to both -> exit determined by other axes 616 n->set(fPlanes[i].a, fPlanes[i].b, << 1178 // >> 1179 sn=kInfinity; // Don't leave by either side 617 } 1180 } 618 return 0; << 619 } 1181 } 620 G4double tmp = -dist/cosa; << 1182 } 621 if (tmax > tmp) { tmax = tmp; iside = i; << 1183 else if (ss1>0&&ss2>=0) >> 1184 { >> 1185 // Outside, in -ve Area >> 1186 if (ds2<0) >> 1187 { >> 1188 sn=0; // away from shape >> 1189 // Left by -ve side >> 1190 snside=kMY; >> 1191 } >> 1192 else >> 1193 { >> 1194 if (ds1>0) >> 1195 { >> 1196 // Ignore +ve plane and use -ve plane intersect >> 1197 // >> 1198 sn=ss1/ds1; // Leave by +ve side >> 1199 snside=kPY; >> 1200 } >> 1201 else >> 1202 { >> 1203 // Must be || to both -> exit determined by other axes >> 1204 // >> 1205 sn=kInfinity; // Don't leave by either side >> 1206 } >> 1207 } >> 1208 } >> 1209 >> 1210 // Update minimum exit distance >> 1211 >> 1212 if (sn<snxt) >> 1213 { >> 1214 snxt=sn; >> 1215 side=snside; 622 } 1216 } 623 } 1217 } 624 1218 625 // Set normal, if required, and return dista << 626 // << 627 if (calcNorm) 1219 if (calcNorm) 628 { 1220 { 629 *validNorm = true; << 1221 switch (side) 630 if (iside < 0) << 1222 { 631 n->set(0, 0, iside + 3); // (-4+3)=-1, ( << 1223 case kPX: 632 else << 1224 cosxz=1.0/std::sqrt(1.0+tanxz*tanxz); 633 n->set(fPlanes[iside].a, fPlanes[iside]. << 1225 *n=G4ThreeVector(cosxz,0,-tanxz*cosxz); >> 1226 break; >> 1227 case kMX: >> 1228 cosxz=-1.0/std::sqrt(1.0+tanxz*tanxz); >> 1229 *n=G4ThreeVector(cosxz,0,tanxz*cosxz); >> 1230 break; >> 1231 case kPY: >> 1232 cosyz=1.0/std::sqrt(1.0+tanyz*tanyz); >> 1233 *n=G4ThreeVector(0,cosyz,-tanyz*cosyz); >> 1234 break; >> 1235 case kMY: >> 1236 cosyz=-1.0/std::sqrt(1.0+tanyz*tanyz); >> 1237 *n=G4ThreeVector(0,cosyz,tanyz*cosyz); >> 1238 break; >> 1239 case kPZ: >> 1240 *n=G4ThreeVector(0,0,1); >> 1241 break; >> 1242 case kMZ: >> 1243 *n=G4ThreeVector(0,0,-1); >> 1244 break; >> 1245 default: >> 1246 DumpInfo(); >> 1247 G4Exception("G4Trd::DistanceToOut(p,v,..)", >> 1248 "GeomSolids1002", JustWarning, >> 1249 "Undefined side for valid surface normal to solid."); >> 1250 break; >> 1251 } 634 } 1252 } 635 return tmax; << 1253 return snxt; 636 } 1254 } 637 1255 638 ////////////////////////////////////////////// << 1256 /////////////////////////////////////////////////////////////////////////// 639 // 1257 // 640 // Calculate exact shortest distance to any bo 1258 // Calculate exact shortest distance to any boundary from inside 641 // - returns 0 if point is outside << 1259 // - Returns 0 is point outside 642 1260 643 G4double G4Trd::DistanceToOut( const G4ThreeVe 1261 G4double G4Trd::DistanceToOut( const G4ThreeVector& p ) const 644 { 1262 { >> 1263 G4double safe=0.0; >> 1264 G4double tanxz,xdist,saf1; >> 1265 G4double tanyz,ydist,saf2; >> 1266 G4double zbase; >> 1267 645 #ifdef G4CSGDEBUG 1268 #ifdef G4CSGDEBUG 646 if( Inside(p) == kOutside ) 1269 if( Inside(p) == kOutside ) 647 { 1270 { 648 std::ostringstream message; << 1271 G4int oldprc = G4cout.precision(16) ; 649 G4long oldprc = message.precision(16); << 1272 G4cout << G4endl ; 650 message << "Point p is outside (!?) of sol << 1273 DumpInfo(); 651 message << "Position:\n"; << 1274 G4cout << "Position:" << G4endl << G4endl ; 652 message << " p.x() = " << p.x()/mm << " << 1275 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ; 653 message << " p.y() = " << p.y()/mm << " << 1276 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ; 654 message << " p.z() = " << p.z()/mm << " << 1277 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ; 655 G4cout.precision(oldprc); << 1278 G4cout.precision(oldprc) ; 656 G4Exception("G4Trd::DistanceToOut(p)", "Ge << 1279 G4Exception("G4Trd::DistanceToOut(p)", "GeomSolids1002", JustWarning, 657 JustWarning, message ); << 1280 "Point p is outside !?" ); 658 DumpInfo(); << 659 } 1281 } 660 #endif 1282 #endif 661 G4double dx = fPlanes[3].a*std::abs(p.x())+f << 662 G4double dy = fPlanes[1].b*std::abs(p.y())+f << 663 G4double dxy = std::max(dx,dy); << 664 1283 665 G4double dz = std::abs(p.z())-fDz; << 1284 safe=fDz-std::fabs(p.z()); // z perpendicular Dist 666 G4double dist = std::max(dz,dxy); << 1285 >> 1286 zbase=fDz+p.z(); 667 1287 668 return (dist < 0) ? -dist : 0.; << 1288 // xdist = distance perpendicular to z axis to closest x plane from p >> 1289 // = (x half width of shape at p.z) - std::fabs(p.x) >> 1290 // >> 1291 tanxz=(fDx2-fDx1)*0.5/fDz; >> 1292 xdist=fDx1+tanxz*zbase-std::fabs(p.x()); >> 1293 saf1=xdist/std::sqrt(1.0+tanxz*tanxz); // x*std::cos(ang_xz) = >> 1294 // shortest (perpendicular) >> 1295 // distance to plane >> 1296 tanyz=(fDy2-fDy1)*0.5/fDz; >> 1297 ydist=fDy1+tanyz*zbase-std::fabs(p.y()); >> 1298 saf2=ydist/std::sqrt(1.0+tanyz*tanyz); >> 1299 >> 1300 // Return minimum x/y/z distance >> 1301 // >> 1302 if (safe>saf1) safe=saf1; >> 1303 if (safe>saf2) safe=saf2; >> 1304 >> 1305 if (safe<0) safe=0; >> 1306 return safe; 669 } 1307 } 670 1308 671 ////////////////////////////////////////////// << 1309 //////////////////////////////////////////////////////////////////////////// 672 // 1310 // 673 // GetEntityType << 1311 // Create a List containing the transformed vertices 674 << 1312 // Ordering [0-3] -fDz cross section 675 G4GeometryType G4Trd::GetEntityType() const << 1313 // [4-7] +fDz cross section such that [0] is below [4], 676 { << 1314 // [1] below [5] etc. 677 return {"G4Trd"}; << 1315 // Note: >> 1316 // Caller has deletion resposibility >> 1317 >> 1318 G4ThreeVectorList* >> 1319 G4Trd::CreateRotatedVertices( const G4AffineTransform& pTransform ) const >> 1320 { >> 1321 G4ThreeVectorList *vertices; >> 1322 vertices=new G4ThreeVectorList(); >> 1323 if (vertices) >> 1324 { >> 1325 vertices->reserve(8); >> 1326 G4ThreeVector vertex0(-fDx1,-fDy1,-fDz); >> 1327 G4ThreeVector vertex1(fDx1,-fDy1,-fDz); >> 1328 G4ThreeVector vertex2(fDx1,fDy1,-fDz); >> 1329 G4ThreeVector vertex3(-fDx1,fDy1,-fDz); >> 1330 G4ThreeVector vertex4(-fDx2,-fDy2,fDz); >> 1331 G4ThreeVector vertex5(fDx2,-fDy2,fDz); >> 1332 G4ThreeVector vertex6(fDx2,fDy2,fDz); >> 1333 G4ThreeVector vertex7(-fDx2,fDy2,fDz); >> 1334 >> 1335 vertices->push_back(pTransform.TransformPoint(vertex0)); >> 1336 vertices->push_back(pTransform.TransformPoint(vertex1)); >> 1337 vertices->push_back(pTransform.TransformPoint(vertex2)); >> 1338 vertices->push_back(pTransform.TransformPoint(vertex3)); >> 1339 vertices->push_back(pTransform.TransformPoint(vertex4)); >> 1340 vertices->push_back(pTransform.TransformPoint(vertex5)); >> 1341 vertices->push_back(pTransform.TransformPoint(vertex6)); >> 1342 vertices->push_back(pTransform.TransformPoint(vertex7)); >> 1343 } >> 1344 else >> 1345 { >> 1346 DumpInfo(); >> 1347 G4Exception("G4Trd::CreateRotatedVertices()", >> 1348 "GeomSolids0003", FatalException, >> 1349 "Error in allocation of vertices. Out of memory !"); >> 1350 } >> 1351 return vertices; 678 } 1352 } 679 1353 680 ////////////////////////////////////////////// 1354 ////////////////////////////////////////////////////////////////////////// 681 // 1355 // 682 // IsFaceted << 1356 // GetEntityType 683 1357 684 G4bool G4Trd::IsFaceted() const << 1358 G4GeometryType G4Trd::GetEntityType() const 685 { 1359 { 686 return true; << 1360 return G4String("G4Trd"); 687 } 1361 } 688 1362 689 ////////////////////////////////////////////// 1363 ////////////////////////////////////////////////////////////////////////// 690 // 1364 // 691 // Make a clone of the object 1365 // Make a clone of the object 692 // 1366 // 693 G4VSolid* G4Trd::Clone() const 1367 G4VSolid* G4Trd::Clone() const 694 { 1368 { 695 return new G4Trd(*this); 1369 return new G4Trd(*this); 696 } 1370 } 697 1371 698 ////////////////////////////////////////////// 1372 ////////////////////////////////////////////////////////////////////////// 699 // 1373 // 700 // Stream object contents to an output stream 1374 // Stream object contents to an output stream 701 1375 702 std::ostream& G4Trd::StreamInfo( std::ostream& 1376 std::ostream& G4Trd::StreamInfo( std::ostream& os ) const 703 { 1377 { 704 G4long oldprc = os.precision(16); << 1378 G4int oldprc = os.precision(16); 705 os << "------------------------------------- 1379 os << "-----------------------------------------------------------\n" 706 << " *** Dump for solid - " << GetName 1380 << " *** Dump for solid - " << GetName() << " ***\n" 707 << " ================================= 1381 << " ===================================================\n" 708 << " Solid type: G4Trd\n" 1382 << " Solid type: G4Trd\n" 709 << " Parameters: \n" 1383 << " Parameters: \n" 710 << " half length X, surface -dZ: " << 1384 << " half length X, surface -dZ: " << fDx1/mm << " mm \n" 711 << " half length X, surface +dZ: " << 1385 << " half length X, surface +dZ: " << fDx2/mm << " mm \n" 712 << " half length Y, surface -dZ: " << 1386 << " half length Y, surface -dZ: " << fDy1/mm << " mm \n" 713 << " half length Y, surface +dZ: " << 1387 << " half length Y, surface +dZ: " << fDy2/mm << " mm \n" 714 << " half length Z : " << << 1388 << " half length Z : " << fDz/mm << " mm \n" 715 << "------------------------------------- 1389 << "-----------------------------------------------------------\n"; 716 os.precision(oldprc); 1390 os.precision(oldprc); 717 1391 718 return os; 1392 return os; 719 } 1393 } 720 1394 721 ////////////////////////////////////////////// << 1395 >> 1396 //////////////////////////////////////////////////////////////////////// 722 // 1397 // 723 // Return a point randomly and uniformly selec << 1398 // GetPointOnSurface >> 1399 // >> 1400 // Return a point (G4ThreeVector) randomly and uniformly >> 1401 // selected on the solid surface 724 1402 725 G4ThreeVector G4Trd::GetPointOnSurface() const 1403 G4ThreeVector G4Trd::GetPointOnSurface() const 726 { 1404 { 727 // Set areas << 1405 G4double px, py, pz, tgX, tgY, secX, secY, select, sumS, tmp; 728 // << 1406 G4double Sxy1, Sxy2, Sxy, Sxz, Syz; 729 G4double sxz = (fDx1 + fDx2)*fHx; << 730 G4double syz = (fDy1 + fDy2)*fHy; << 731 G4double ssurf[6] = { 4.*fDx1*fDy1, sxz, sxz << 732 ssurf[1] += ssurf[0]; << 733 ssurf[2] += ssurf[1]; << 734 ssurf[3] += ssurf[2]; << 735 ssurf[4] += ssurf[3]; << 736 ssurf[5] += ssurf[4]; << 737 1407 738 // Select face << 1408 tgX = 0.5*(fDx2-fDx1)/fDz; 739 // << 1409 secX = std::sqrt(1+tgX*tgX); 740 G4double select = ssurf[5]*G4QuickRand(); << 1410 tgY = 0.5*(fDy2-fDy1)/fDz; 741 G4int k = 5; << 1411 secY = std::sqrt(1+tgY*tgY); 742 k -= (G4int)(select <= ssurf[4]); << 1412 743 k -= (G4int)(select <= ssurf[3]); << 1413 // calculate 0.25 of side surfaces, sumS is 0.25 of total surface 744 k -= (G4int)(select <= ssurf[2]); << 1414 745 k -= (G4int)(select <= ssurf[1]); << 1415 Sxy1 = fDx1*fDy1; 746 k -= (G4int)(select <= ssurf[0]); << 1416 Sxy2 = fDx2*fDy2; 747 << 1417 Sxy = Sxy1 + Sxy2; 748 // Generate point on selected surface << 1418 Sxz = (fDx1 + fDx2)*fDz*secY; 749 // << 1419 Syz = (fDy1 + fDy2)*fDz*secX; 750 G4double u = G4QuickRand(); << 1420 sumS = Sxy + Sxz + Syz; 751 G4double v = G4QuickRand(); << 1421 752 switch(k) << 1422 select = sumS*G4UniformRand(); >> 1423 >> 1424 if( select < Sxy ) // Sxy1 or Sxy2 753 { 1425 { 754 case 0: // base at -Z << 1426 if( select < Sxy1 ) 755 { << 756 return { (2.*u - 1.)*fDx1, (2.*v - 1.)*f << 757 } << 758 case 1: // X face at -Y << 759 { << 760 if (u + v > 1.) { u = 1. - u; v = 1. - v << 761 G4ThreeVector p0(-fDx1,-fDy1,-fDz); << 762 G4ThreeVector p1( fDx2,-fDy2, fDz); << 763 return (select <= ssurf[0] + fDx1*fHx) ? << 764 (1. - u - v)*p0 + u*p1 + v*G4ThreeVect << 765 (1. - u - v)*p0 + u*p1 + v*G4ThreeVect << 766 } << 767 case 2: // X face at +Y << 768 { 1427 { 769 if (u + v > 1.) { u = 1. - u; v = 1. - v << 1428 pz = -fDz; 770 G4ThreeVector p0( fDx1, fDy1,-fDz); << 1429 px = -fDx1 + 2*fDx1*G4UniformRand(); 771 G4ThreeVector p1(-fDx2, fDy2, fDz); << 1430 py = -fDy1 + 2*fDy1*G4UniformRand(); 772 return (select <= ssurf[1] + fDx1*fHx) ? << 773 (1. - u - v)*p0 + u*p1 + v*G4ThreeVect << 774 (1. - u - v)*p0 + u*p1 + v*G4ThreeVect << 775 } 1431 } 776 case 3: // Y face at -X << 1432 else 777 { 1433 { 778 if (u + v > 1.) { u = 1. - u; v = 1. - v << 1434 pz = fDz; 779 G4ThreeVector p0(-fDx1, fDy1,-fDz); << 1435 px = -fDx2 + 2*fDx2*G4UniformRand(); 780 G4ThreeVector p1(-fDx2,-fDy2, fDz); << 1436 py = -fDy2 + 2*fDy2*G4UniformRand(); 781 return (select <= ssurf[2] + fDy1*fHy) ? << 782 (1. - u - v)*p0 + u*p1 + v*G4ThreeVect << 783 (1. - u - v)*p0 + u*p1 + v*G4ThreeVect << 784 } << 785 case 4: // Y face at +X << 786 { << 787 if (u + v > 1.) { u = 1. - u; v = 1. - v << 788 G4ThreeVector p0( fDx1,-fDy1,-fDz); << 789 G4ThreeVector p1( fDx2, fDy2, fDz); << 790 return (select <= ssurf[3] + fDy1*fHy) ? << 791 (1. - u - v)*p0 + u*p1 + v*G4ThreeVect << 792 (1. - u - v)*p0 + u*p1 + v*G4ThreeVect << 793 } << 794 case 5: // base at +Z << 795 { << 796 return { (2.*u - 1.)*fDx2, (2.*v - 1.)*f << 797 } 1437 } 798 } 1438 } 799 return {0., 0., 0.}; << 1439 else if ( ( select - Sxy ) < Sxz ) // Sxz >> 1440 { >> 1441 pz = -fDz + 2*fDz*G4UniformRand(); >> 1442 tmp = fDx1 + (pz + fDz)*tgX; >> 1443 px = -tmp + 2*tmp*G4UniformRand(); >> 1444 tmp = fDy1 + (pz + fDz)*tgY; >> 1445 >> 1446 if(G4UniformRand() > 0.5) { py = tmp; } >> 1447 else { py = -tmp; } >> 1448 } >> 1449 else // Syz >> 1450 { >> 1451 pz = -fDz + 2*fDz*G4UniformRand(); >> 1452 tmp = fDy1 + (pz + fDz)*tgY; >> 1453 py = -tmp + 2*tmp*G4UniformRand(); >> 1454 tmp = fDx1 + (pz + fDz)*tgX; >> 1455 >> 1456 if(G4UniformRand() > 0.5) { px = tmp; } >> 1457 else { px = -tmp; } >> 1458 } >> 1459 return G4ThreeVector(px,py,pz); 800 } 1460 } 801 1461 802 ////////////////////////////////////////////// << 1462 /////////////////////////////////////////////////////////////////////// 803 // 1463 // 804 // Methods for visualisation 1464 // Methods for visualisation 805 1465 806 void G4Trd::DescribeYourselfTo ( G4VGraphicsSc 1466 void G4Trd::DescribeYourselfTo ( G4VGraphicsScene& scene ) const 807 { 1467 { 808 scene.AddSolid (*this); 1468 scene.AddSolid (*this); 809 } 1469 } 810 1470 811 G4Polyhedron* G4Trd::CreatePolyhedron () const 1471 G4Polyhedron* G4Trd::CreatePolyhedron () const 812 { 1472 { 813 return new G4PolyhedronTrd2 (fDx1, fDx2, fDy 1473 return new G4PolyhedronTrd2 (fDx1, fDx2, fDy1, fDy2, fDz); 814 } 1474 } 815 1475 816 #endif << 1476 G4NURBS* G4Trd::CreateNURBS () const >> 1477 { >> 1478 // return new G4NURBSbox (fDx, fDy, fDz); >> 1479 return 0; >> 1480 } 817 1481