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 of G4PolyconeSide, the face << 27 // one conical side of a polycone << 28 // 26 // 29 // Author: David C. Williams (davidw@scipp.ucs << 27 // $Id: G4PolyconeSide.cc,v 1.25 2010-07-12 15:25:37 gcosmo Exp $ >> 28 // GEANT4 tag $Name: geant4-09-04-patch-02 $ >> 29 // >> 30 // >> 31 // -------------------------------------------------------------------- >> 32 // GEANT 4 class source file >> 33 // >> 34 // >> 35 // G4PolyconeSide.cc >> 36 // >> 37 // Implementation of the face representing one conical side of a polycone >> 38 // 30 // ------------------------------------------- 39 // -------------------------------------------------------------------- 31 40 32 #include "G4PolyconeSide.hh" 41 #include "G4PolyconeSide.hh" 33 #include "meshdefs.hh" << 34 #include "G4PhysicalConstants.hh" << 35 #include "G4IntersectingCone.hh" 42 #include "G4IntersectingCone.hh" 36 #include "G4ClippablePolygon.hh" 43 #include "G4ClippablePolygon.hh" 37 #include "G4AffineTransform.hh" 44 #include "G4AffineTransform.hh" >> 45 #include "meshdefs.hh" 38 #include "G4SolidExtentList.hh" 46 #include "G4SolidExtentList.hh" 39 #include "G4GeometryTolerance.hh" 47 #include "G4GeometryTolerance.hh" 40 48 41 #include "Randomize.hh" 49 #include "Randomize.hh" 42 50 43 // This new field helps to use the class G4PlS << 44 // 51 // 45 G4PlSideManager G4PolyconeSide::subInstanceMan << 46 << 47 // This macro changes the references to fields << 48 // in the class G4PlSideData. << 49 // << 50 #define G4MT_pcphix ((subInstanceManager.offse << 51 #define G4MT_pcphiy ((subInstanceManager.offse << 52 #define G4MT_pcphiz ((subInstanceManager.offse << 53 #define G4MT_pcphik ((subInstanceManager.offse << 54 << 55 // Returns the private data instance manager. << 56 // << 57 const G4PlSideManager& G4PolyconeSide::GetSubI << 58 { << 59 return subInstanceManager; << 60 } << 61 << 62 // Constructor 52 // Constructor 63 // 53 // 64 // Values for r1,z1 and r2,z2 should be specif 54 // Values for r1,z1 and r2,z2 should be specified in clockwise 65 // order in (r,z). 55 // order in (r,z). 66 // 56 // 67 G4PolyconeSide::G4PolyconeSide( const G4Polyco << 57 G4PolyconeSide::G4PolyconeSide( const G4PolyconeSideRZ *prevRZ, 68 const G4Polyco << 58 const G4PolyconeSideRZ *tail, 69 const G4Polyco << 59 const G4PolyconeSideRZ *head, 70 const G4Polyco << 60 const G4PolyconeSideRZ *nextRZ, 71 G4double 61 G4double thePhiStart, 72 G4double 62 G4double theDeltaPhi, 73 G4bool t 63 G4bool thePhiIsOpen, 74 G4bool i 64 G4bool isAllBehind ) >> 65 : ncorners(0), corners(0) 75 { 66 { 76 instanceID = subInstanceManager.CreateSubIns << 77 << 78 kCarTolerance = G4GeometryTolerance::GetInst 67 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance(); 79 G4MT_pcphix = 0.0; G4MT_pcphiy = 0.0; G4MT_p << 68 fSurfaceArea = 0.0; >> 69 fPhi.first = G4ThreeVector(0,0,0); >> 70 fPhi.second= 0.0; 80 71 81 // 72 // 82 // Record values 73 // Record values 83 // 74 // 84 r[0] = tail->r; z[0] = tail->z; 75 r[0] = tail->r; z[0] = tail->z; 85 r[1] = head->r; z[1] = head->z; 76 r[1] = head->r; z[1] = head->z; 86 77 87 phiIsOpen = thePhiIsOpen; 78 phiIsOpen = thePhiIsOpen; 88 if (phiIsOpen) 79 if (phiIsOpen) 89 { 80 { 90 deltaPhi = theDeltaPhi; 81 deltaPhi = theDeltaPhi; 91 startPhi = thePhiStart; 82 startPhi = thePhiStart; 92 83 93 // 84 // 94 // Set phi values to our conventions 85 // Set phi values to our conventions 95 // 86 // 96 while (deltaPhi < 0.0) // Loop checking << 87 while (deltaPhi < 0.0) deltaPhi += twopi; 97 deltaPhi += twopi; << 88 while (startPhi < 0.0) startPhi += twopi; 98 while (startPhi < 0.0) // Loop checking << 99 startPhi += twopi; << 100 89 101 // 90 // 102 // Calculate corner coordinates 91 // Calculate corner coordinates 103 // 92 // 104 ncorners = 4; 93 ncorners = 4; 105 corners = new G4ThreeVector[ncorners]; 94 corners = new G4ThreeVector[ncorners]; 106 95 107 corners[0] = G4ThreeVector( tail->r*std::c 96 corners[0] = G4ThreeVector( tail->r*std::cos(startPhi), 108 tail->r*std::s 97 tail->r*std::sin(startPhi), tail->z ); 109 corners[1] = G4ThreeVector( head->r*std::c 98 corners[1] = G4ThreeVector( head->r*std::cos(startPhi), 110 head->r*std::s 99 head->r*std::sin(startPhi), head->z ); 111 corners[2] = G4ThreeVector( tail->r*std::c 100 corners[2] = G4ThreeVector( tail->r*std::cos(startPhi+deltaPhi), 112 tail->r*std::s 101 tail->r*std::sin(startPhi+deltaPhi), tail->z ); 113 corners[3] = G4ThreeVector( head->r*std::c 102 corners[3] = G4ThreeVector( head->r*std::cos(startPhi+deltaPhi), 114 head->r*std::s 103 head->r*std::sin(startPhi+deltaPhi), head->z ); 115 } 104 } 116 else 105 else 117 { 106 { 118 deltaPhi = twopi; 107 deltaPhi = twopi; 119 startPhi = 0.0; 108 startPhi = 0.0; 120 } 109 } 121 110 122 allBehind = isAllBehind; 111 allBehind = isAllBehind; 123 112 124 // 113 // 125 // Make our intersecting cone 114 // Make our intersecting cone 126 // 115 // 127 cone = new G4IntersectingCone( r, z ); 116 cone = new G4IntersectingCone( r, z ); 128 117 129 // 118 // 130 // Calculate vectors in r,z space 119 // Calculate vectors in r,z space 131 // 120 // 132 rS = r[1]-r[0]; zS = z[1]-z[0]; 121 rS = r[1]-r[0]; zS = z[1]-z[0]; 133 length = std::sqrt( rS*rS + zS*zS); 122 length = std::sqrt( rS*rS + zS*zS); 134 rS /= length; zS /= length; 123 rS /= length; zS /= length; 135 124 136 rNorm = +zS; 125 rNorm = +zS; 137 zNorm = -rS; 126 zNorm = -rS; 138 127 139 G4double lAdj; 128 G4double lAdj; 140 129 141 prevRS = r[0]-prevRZ->r; 130 prevRS = r[0]-prevRZ->r; 142 prevZS = z[0]-prevRZ->z; 131 prevZS = z[0]-prevRZ->z; 143 lAdj = std::sqrt( prevRS*prevRS + prevZS*pre 132 lAdj = std::sqrt( prevRS*prevRS + prevZS*prevZS ); 144 prevRS /= lAdj; 133 prevRS /= lAdj; 145 prevZS /= lAdj; 134 prevZS /= lAdj; 146 135 147 rNormEdge[0] = rNorm + prevZS; 136 rNormEdge[0] = rNorm + prevZS; 148 zNormEdge[0] = zNorm - prevRS; 137 zNormEdge[0] = zNorm - prevRS; 149 lAdj = std::sqrt( rNormEdge[0]*rNormEdge[0] 138 lAdj = std::sqrt( rNormEdge[0]*rNormEdge[0] + zNormEdge[0]*zNormEdge[0] ); 150 rNormEdge[0] /= lAdj; 139 rNormEdge[0] /= lAdj; 151 zNormEdge[0] /= lAdj; 140 zNormEdge[0] /= lAdj; 152 141 153 nextRS = nextRZ->r-r[1]; 142 nextRS = nextRZ->r-r[1]; 154 nextZS = nextRZ->z-z[1]; 143 nextZS = nextRZ->z-z[1]; 155 lAdj = std::sqrt( nextRS*nextRS + nextZS*nex 144 lAdj = std::sqrt( nextRS*nextRS + nextZS*nextZS ); 156 nextRS /= lAdj; 145 nextRS /= lAdj; 157 nextZS /= lAdj; 146 nextZS /= lAdj; 158 147 159 rNormEdge[1] = rNorm + nextZS; 148 rNormEdge[1] = rNorm + nextZS; 160 zNormEdge[1] = zNorm - nextRS; 149 zNormEdge[1] = zNorm - nextRS; 161 lAdj = std::sqrt( rNormEdge[1]*rNormEdge[1] 150 lAdj = std::sqrt( rNormEdge[1]*rNormEdge[1] + zNormEdge[1]*zNormEdge[1] ); 162 rNormEdge[1] /= lAdj; 151 rNormEdge[1] /= lAdj; 163 zNormEdge[1] /= lAdj; 152 zNormEdge[1] /= lAdj; 164 } 153 } 165 154 >> 155 >> 156 // 166 // Fake default constructor - sets only member 157 // Fake default constructor - sets only member data and allocates memory 167 // for usage restri 158 // for usage restricted to object persistency. 168 // 159 // 169 G4PolyconeSide::G4PolyconeSide( __void__& ) 160 G4PolyconeSide::G4PolyconeSide( __void__& ) 170 : startPhi(0.), deltaPhi(0.), << 161 : startPhi(0.), deltaPhi(0.), phiIsOpen(false), allBehind(false), 171 rNorm(0.), zNorm(0.), rS(0.), zS(0.), leng << 162 cone(0), rNorm(0.), zNorm(0.), rS(0.), zS(0.), length(0.), 172 prevRS(0.), prevZS(0.), nextRS(0.), nextZS << 163 prevRS(0.), prevZS(0.), nextRS(0.), nextZS(0.), ncorners(0), corners(0), 173 kCarTolerance(0.), instanceID(0) << 164 kCarTolerance(0.), fSurfaceArea(0.) 174 { 165 { 175 r[0] = r[1] = 0.; 166 r[0] = r[1] = 0.; 176 z[0] = z[1] = 0.; 167 z[0] = z[1] = 0.; 177 rNormEdge[0]= rNormEdge[1] = 0.; 168 rNormEdge[0]= rNormEdge[1] = 0.; 178 zNormEdge[0]= zNormEdge[1] = 0.; 169 zNormEdge[0]= zNormEdge[1] = 0.; 179 } 170 } 180 171 >> 172 >> 173 // 181 // Destructor 174 // Destructor 182 // 175 // 183 G4PolyconeSide::~G4PolyconeSide() 176 G4PolyconeSide::~G4PolyconeSide() 184 { 177 { 185 delete cone; 178 delete cone; 186 if (phiIsOpen) { delete [] corners; } 179 if (phiIsOpen) { delete [] corners; } 187 } 180 } 188 181 >> 182 >> 183 // 189 // Copy constructor 184 // Copy constructor 190 // 185 // 191 G4PolyconeSide::G4PolyconeSide( const G4Polyco << 186 G4PolyconeSide::G4PolyconeSide( const G4PolyconeSide &source ) >> 187 : G4VCSGface(), ncorners(0), corners(0) 192 { 188 { 193 instanceID = subInstanceManager.CreateSubIns << 194 << 195 CopyStuff( source ); 189 CopyStuff( source ); 196 } 190 } 197 191 >> 192 >> 193 // 198 // Assignment operator 194 // Assignment operator 199 // 195 // 200 G4PolyconeSide& G4PolyconeSide::operator=( con << 196 G4PolyconeSide& G4PolyconeSide::operator=( const G4PolyconeSide &source ) 201 { 197 { 202 if (this == &source) { return *this; } 198 if (this == &source) { return *this; } 203 199 204 delete cone; 200 delete cone; 205 if (phiIsOpen) { delete [] corners; } 201 if (phiIsOpen) { delete [] corners; } 206 202 207 CopyStuff( source ); 203 CopyStuff( source ); 208 204 209 return *this; 205 return *this; 210 } 206 } 211 207 >> 208 >> 209 // 212 // CopyStuff 210 // CopyStuff 213 // 211 // 214 void G4PolyconeSide::CopyStuff( const G4Polyco << 212 void G4PolyconeSide::CopyStuff( const G4PolyconeSide &source ) 215 { 213 { 216 r[0] = source.r[0]; 214 r[0] = source.r[0]; 217 r[1] = source.r[1]; 215 r[1] = source.r[1]; 218 z[0] = source.z[0]; 216 z[0] = source.z[0]; 219 z[1] = source.z[1]; 217 z[1] = source.z[1]; 220 218 221 startPhi = source.startPhi; 219 startPhi = source.startPhi; 222 deltaPhi = source.deltaPhi; 220 deltaPhi = source.deltaPhi; 223 phiIsOpen = source.phiIsOpen; 221 phiIsOpen = source.phiIsOpen; 224 allBehind = source.allBehind; 222 allBehind = source.allBehind; 225 223 226 kCarTolerance = source.kCarTolerance; 224 kCarTolerance = source.kCarTolerance; 227 fSurfaceArea = source.fSurfaceArea; 225 fSurfaceArea = source.fSurfaceArea; 228 226 229 cone = new G4IntersectingCone( *source.co 227 cone = new G4IntersectingCone( *source.cone ); 230 228 231 rNorm = source.rNorm; 229 rNorm = source.rNorm; 232 zNorm = source.zNorm; 230 zNorm = source.zNorm; 233 rS = source.rS; 231 rS = source.rS; 234 zS = source.zS; 232 zS = source.zS; 235 length = source.length; 233 length = source.length; 236 prevRS = source.prevRS; 234 prevRS = source.prevRS; 237 prevZS = source.prevZS; 235 prevZS = source.prevZS; 238 nextRS = source.nextRS; 236 nextRS = source.nextRS; 239 nextZS = source.nextZS; 237 nextZS = source.nextZS; 240 238 241 rNormEdge[0] = source.rNormEdge[0]; 239 rNormEdge[0] = source.rNormEdge[0]; 242 rNormEdge[1] = source.rNormEdge[1]; 240 rNormEdge[1] = source.rNormEdge[1]; 243 zNormEdge[0] = source.zNormEdge[0]; 241 zNormEdge[0] = source.zNormEdge[0]; 244 zNormEdge[1] = source.zNormEdge[1]; 242 zNormEdge[1] = source.zNormEdge[1]; 245 243 246 if (phiIsOpen) 244 if (phiIsOpen) 247 { 245 { 248 ncorners = 4; 246 ncorners = 4; 249 corners = new G4ThreeVector[ncorners]; 247 corners = new G4ThreeVector[ncorners]; 250 248 251 corners[0] = source.corners[0]; 249 corners[0] = source.corners[0]; 252 corners[1] = source.corners[1]; 250 corners[1] = source.corners[1]; 253 corners[2] = source.corners[2]; 251 corners[2] = source.corners[2]; 254 corners[3] = source.corners[3]; 252 corners[3] = source.corners[3]; 255 } 253 } 256 } 254 } 257 255 >> 256 >> 257 // 258 // Intersect 258 // Intersect 259 // 259 // 260 G4bool G4PolyconeSide::Intersect( const G4Thre << 260 G4bool G4PolyconeSide::Intersect( const G4ThreeVector &p, 261 const G4Thre << 261 const G4ThreeVector &v, 262 G4bool 262 G4bool outgoing, 263 G4doub 263 G4double surfTolerance, 264 G4doub << 264 G4double &distance, 265 G4doub << 265 G4double &distFromSurface, 266 G4Thre << 266 G4ThreeVector &normal, 267 G4bool << 267 G4bool &isAllBehind ) 268 { 268 { 269 G4double s1=0., s2=0.; << 269 G4double s1, s2; 270 G4double normSign = outgoing ? +1 : -1; 270 G4double normSign = outgoing ? +1 : -1; 271 271 272 isAllBehind = allBehind; 272 isAllBehind = allBehind; 273 273 274 // 274 // 275 // Check for two possible intersections 275 // Check for two possible intersections 276 // 276 // 277 G4int nside = cone->LineHitsCone( p, v, &s1, 277 G4int nside = cone->LineHitsCone( p, v, &s1, &s2 ); 278 if (nside == 0) return false; 278 if (nside == 0) return false; 279 279 280 // 280 // 281 // Check the first side first, since it is ( 281 // Check the first side first, since it is (supposed to be) closest 282 // 282 // 283 G4ThreeVector hit = p + s1*v; 283 G4ThreeVector hit = p + s1*v; 284 284 285 if (PointOnCone( hit, normSign, p, v, normal 285 if (PointOnCone( hit, normSign, p, v, normal )) 286 { 286 { 287 // 287 // 288 // Good intersection! What about the norma 288 // Good intersection! What about the normal? 289 // 289 // 290 if (normSign*v.dot(normal) > 0) 290 if (normSign*v.dot(normal) > 0) 291 { 291 { 292 // 292 // 293 // We have a valid intersection, but it 293 // We have a valid intersection, but it could very easily 294 // be behind the point. To decide if we 294 // be behind the point. To decide if we tolerate this, 295 // we have to see if the point p is on t 295 // we have to see if the point p is on the surface near 296 // the intersecting point. 296 // the intersecting point. 297 // 297 // 298 // What does it mean exactly for the poi 298 // What does it mean exactly for the point p to be "near" 299 // the intersection? It means that if we 299 // the intersection? It means that if we draw a line from 300 // p to the hit, the line remains entire 300 // p to the hit, the line remains entirely within the 301 // tolerance bounds of the cone. To test 301 // tolerance bounds of the cone. To test this, we can 302 // ask if the normal is correct near p. 302 // ask if the normal is correct near p. 303 // 303 // 304 G4double pr = p.perp(); 304 G4double pr = p.perp(); 305 if (pr < DBL_MIN) pr = DBL_MIN; 305 if (pr < DBL_MIN) pr = DBL_MIN; 306 G4ThreeVector pNormal( rNorm*p.x()/pr, r 306 G4ThreeVector pNormal( rNorm*p.x()/pr, rNorm*p.y()/pr, zNorm ); 307 if (normSign*v.dot(pNormal) > 0) 307 if (normSign*v.dot(pNormal) > 0) 308 { 308 { 309 // 309 // 310 // p and intersection in same hemisphe 310 // p and intersection in same hemisphere 311 // 311 // 312 G4double distOutside2; 312 G4double distOutside2; 313 distFromSurface = -normSign*DistanceAw 313 distFromSurface = -normSign*DistanceAway( p, false, distOutside2 ); 314 if (distOutside2 < surfTolerance*surfT 314 if (distOutside2 < surfTolerance*surfTolerance) 315 { 315 { 316 if (distFromSurface > -surfTolerance 316 if (distFromSurface > -surfTolerance) 317 { 317 { 318 // 318 // 319 // We are just inside or away from 319 // We are just inside or away from the 320 // surface. Accept *any* value of 320 // surface. Accept *any* value of distance. 321 // 321 // 322 distance = s1; 322 distance = s1; 323 return true; 323 return true; 324 } 324 } 325 } 325 } 326 } 326 } 327 else 327 else 328 distFromSurface = s1; 328 distFromSurface = s1; 329 329 330 // 330 // 331 // Accept positive distances 331 // Accept positive distances 332 // 332 // 333 if (s1 > 0) 333 if (s1 > 0) 334 { 334 { 335 distance = s1; 335 distance = s1; 336 return true; 336 return true; 337 } 337 } 338 } 338 } 339 } 339 } 340 340 341 if (nside==1) return false; 341 if (nside==1) return false; 342 342 343 // 343 // 344 // Well, try the second hit 344 // Well, try the second hit 345 // 345 // 346 hit = p + s2*v; 346 hit = p + s2*v; 347 347 348 if (PointOnCone( hit, normSign, p, v, normal 348 if (PointOnCone( hit, normSign, p, v, normal )) 349 { 349 { 350 // 350 // 351 // Good intersection! What about the norma 351 // Good intersection! What about the normal? 352 // 352 // 353 if (normSign*v.dot(normal) > 0) 353 if (normSign*v.dot(normal) > 0) 354 { 354 { 355 G4double pr = p.perp(); 355 G4double pr = p.perp(); 356 if (pr < DBL_MIN) pr = DBL_MIN; 356 if (pr < DBL_MIN) pr = DBL_MIN; 357 G4ThreeVector pNormal( rNorm*p.x()/pr, r 357 G4ThreeVector pNormal( rNorm*p.x()/pr, rNorm*p.y()/pr, zNorm ); 358 if (normSign*v.dot(pNormal) > 0) 358 if (normSign*v.dot(pNormal) > 0) 359 { 359 { 360 G4double distOutside2; 360 G4double distOutside2; 361 distFromSurface = -normSign*DistanceAw 361 distFromSurface = -normSign*DistanceAway( p, false, distOutside2 ); 362 if (distOutside2 < surfTolerance*surfT 362 if (distOutside2 < surfTolerance*surfTolerance) 363 { 363 { 364 if (distFromSurface > -surfTolerance 364 if (distFromSurface > -surfTolerance) 365 { 365 { 366 distance = s2; 366 distance = s2; 367 return true; 367 return true; 368 } 368 } 369 } 369 } 370 } 370 } 371 else 371 else 372 distFromSurface = s2; 372 distFromSurface = s2; 373 373 374 if (s2 > 0) 374 if (s2 > 0) 375 { 375 { 376 distance = s2; 376 distance = s2; 377 return true; 377 return true; 378 } 378 } 379 } 379 } 380 } 380 } 381 381 382 // 382 // 383 // Better luck next time 383 // Better luck next time 384 // 384 // 385 return false; 385 return false; 386 } 386 } 387 387 388 // Distance << 388 389 // << 389 G4double G4PolyconeSide::Distance( const G4ThreeVector &p, G4bool outgoing ) 390 G4double G4PolyconeSide::Distance( const G4Thr << 391 { 390 { 392 G4double normSign = outgoing ? -1 : +1; 391 G4double normSign = outgoing ? -1 : +1; 393 G4double distFrom, distOut2; 392 G4double distFrom, distOut2; 394 393 395 // 394 // 396 // We have two tries for each hemisphere. Tr 395 // We have two tries for each hemisphere. Try the closest first. 397 // 396 // 398 distFrom = normSign*DistanceAway( p, false, 397 distFrom = normSign*DistanceAway( p, false, distOut2 ); 399 if (distFrom > -0.5*kCarTolerance ) 398 if (distFrom > -0.5*kCarTolerance ) 400 { 399 { 401 // 400 // 402 // Good answer 401 // Good answer 403 // 402 // 404 if (distOut2 > 0) 403 if (distOut2 > 0) 405 return std::sqrt( distFrom*distFrom + di 404 return std::sqrt( distFrom*distFrom + distOut2 ); 406 else 405 else 407 return std::fabs(distFrom); 406 return std::fabs(distFrom); 408 } 407 } 409 408 410 // 409 // 411 // Try second side. 410 // Try second side. 412 // 411 // 413 distFrom = normSign*DistanceAway( p, true, 412 distFrom = normSign*DistanceAway( p, true, distOut2 ); 414 if (distFrom > -0.5*kCarTolerance) 413 if (distFrom > -0.5*kCarTolerance) 415 { 414 { 416 415 417 if (distOut2 > 0) 416 if (distOut2 > 0) 418 return std::sqrt( distFrom*distFrom + di 417 return std::sqrt( distFrom*distFrom + distOut2 ); 419 else 418 else 420 return std::fabs(distFrom); 419 return std::fabs(distFrom); 421 } 420 } 422 421 423 return kInfinity; 422 return kInfinity; 424 } 423 } 425 424 >> 425 >> 426 // 426 // Inside 427 // Inside 427 // 428 // 428 EInside G4PolyconeSide::Inside( const G4ThreeV << 429 EInside G4PolyconeSide::Inside( const G4ThreeVector &p, 429 G4double 430 G4double tolerance, 430 G4double << 431 G4double *bestDistance ) 431 { 432 { 432 G4double distFrom, distOut2, dist2; << 433 // 433 G4double edgeRZnorm; << 434 // Check both sides >> 435 // >> 436 G4double distFrom[2], distOut2[2], dist2[2]; >> 437 G4double edgeRZnorm[2]; 434 438 435 distFrom = DistanceAway( p, distOut2, &edge << 439 distFrom[0] = DistanceAway( p, false, distOut2[0], edgeRZnorm ); 436 dist2 = distFrom*distFrom + distOut2; << 440 distFrom[1] = DistanceAway( p, true, distOut2[1], edgeRZnorm+1 ); 437 << 438 *bestDistance = std::sqrt( dist2); << 439 441 >> 442 dist2[0] = distFrom[0]*distFrom[0] + distOut2[0]; >> 443 dist2[1] = distFrom[1]*distFrom[1] + distOut2[1]; >> 444 >> 445 // >> 446 // Who's closest? >> 447 // >> 448 G4int i = std::fabs(dist2[0]) < std::fabs(dist2[1]) ? 0 : 1; >> 449 >> 450 *bestDistance = std::sqrt( dist2[i] ); >> 451 >> 452 // 440 // Okay then, inside or out? 453 // Okay then, inside or out? 441 // 454 // 442 if ( (std::fabs(edgeRZnorm) < tolerance) << 455 if ( (std::fabs(edgeRZnorm[i]) < tolerance) 443 && (distOut2< tolerance*tolerance) ) << 456 && (distOut2[i] < tolerance*tolerance) ) 444 return kSurface; 457 return kSurface; 445 else if (edgeRZnorm < 0) << 458 else if (edgeRZnorm[i] < 0) 446 return kInside; 459 return kInside; 447 else 460 else 448 return kOutside; 461 return kOutside; 449 } 462 } 450 463 >> 464 >> 465 // 451 // Normal 466 // Normal 452 // 467 // 453 G4ThreeVector G4PolyconeSide::Normal( const G4 << 468 G4ThreeVector G4PolyconeSide::Normal( const G4ThreeVector &p, 454 G4 << 469 G4double *bestDistance ) 455 { 470 { 456 if (p == G4ThreeVector(0.,0.,0.)) { return 471 if (p == G4ThreeVector(0.,0.,0.)) { return p; } 457 472 458 G4double dFrom, dOut2; 473 G4double dFrom, dOut2; 459 474 460 dFrom = DistanceAway( p, false, dOut2 ); 475 dFrom = DistanceAway( p, false, dOut2 ); 461 476 462 *bestDistance = std::sqrt( dFrom*dFrom + dOu 477 *bestDistance = std::sqrt( dFrom*dFrom + dOut2 ); 463 478 464 G4double rds = p.perp(); << 479 G4double rad = p.perp(); 465 if (rds!=0.) { return {rNorm*p.x()/rds,rNorm << 480 if (rad!=0.) { return G4ThreeVector(rNorm*p.x()/rad,rNorm*p.y()/rad,zNorm); } 466 return G4ThreeVector( 0.,0., zNorm ).unit(); 481 return G4ThreeVector( 0.,0., zNorm ).unit(); 467 } 482 } 468 483 >> 484 >> 485 // 469 // Extent 486 // Extent 470 // 487 // 471 G4double G4PolyconeSide::Extent( const G4Three 488 G4double G4PolyconeSide::Extent( const G4ThreeVector axis ) 472 { 489 { 473 if (axis.perp2() < DBL_MIN) 490 if (axis.perp2() < DBL_MIN) 474 { 491 { 475 // 492 // 476 // Special case 493 // Special case 477 // 494 // 478 return axis.z() < 0 ? -cone->ZLo() : cone- 495 return axis.z() < 0 ? -cone->ZLo() : cone->ZHi(); 479 } 496 } 480 497 481 // 498 // 482 // Is the axis pointing inside our phi gap? 499 // Is the axis pointing inside our phi gap? 483 // 500 // 484 if (phiIsOpen) 501 if (phiIsOpen) 485 { 502 { 486 G4double phi = GetPhi(axis); 503 G4double phi = GetPhi(axis); 487 while( phi < startPhi ) // Loop checkin << 504 while( phi < startPhi ) phi += twopi; 488 phi += twopi; << 489 505 490 if (phi > deltaPhi+startPhi) 506 if (phi > deltaPhi+startPhi) 491 { 507 { 492 // 508 // 493 // Yeah, looks so. Make four three vecto 509 // Yeah, looks so. Make four three vectors defining the phi 494 // opening 510 // opening 495 // 511 // 496 G4double cosP = std::cos(startPhi), sinP 512 G4double cosP = std::cos(startPhi), sinP = std::sin(startPhi); 497 G4ThreeVector a( r[0]*cosP, r[0]*sinP, z 513 G4ThreeVector a( r[0]*cosP, r[0]*sinP, z[0] ); 498 G4ThreeVector b( r[1]*cosP, r[1]*sinP, z 514 G4ThreeVector b( r[1]*cosP, r[1]*sinP, z[1] ); 499 cosP = std::cos(startPhi+deltaPhi); sinP 515 cosP = std::cos(startPhi+deltaPhi); sinP = std::sin(startPhi+deltaPhi); 500 G4ThreeVector c( r[0]*cosP, r[0]*sinP, z 516 G4ThreeVector c( r[0]*cosP, r[0]*sinP, z[0] ); 501 G4ThreeVector d( r[1]*cosP, r[1]*sinP, z 517 G4ThreeVector d( r[1]*cosP, r[1]*sinP, z[1] ); 502 518 503 G4double ad = axis.dot(a), 519 G4double ad = axis.dot(a), 504 bd = axis.dot(b), << 520 bd = axis.dot(b), 505 cd = axis.dot(c), << 521 cd = axis.dot(c), 506 dd = axis.dot(d); << 522 dd = axis.dot(d); 507 523 508 if (bd > ad) ad = bd; 524 if (bd > ad) ad = bd; 509 if (cd > ad) ad = cd; 525 if (cd > ad) ad = cd; 510 if (dd > ad) ad = dd; 526 if (dd > ad) ad = dd; 511 527 512 return ad; 528 return ad; 513 } 529 } 514 } 530 } 515 531 516 // 532 // 517 // Check either end 533 // Check either end 518 // 534 // 519 G4double aPerp = axis.perp(); 535 G4double aPerp = axis.perp(); 520 536 521 G4double a = aPerp*r[0] + axis.z()*z[0]; 537 G4double a = aPerp*r[0] + axis.z()*z[0]; 522 G4double b = aPerp*r[1] + axis.z()*z[1]; 538 G4double b = aPerp*r[1] + axis.z()*z[1]; 523 539 524 if (b > a) a = b; 540 if (b > a) a = b; 525 541 526 return a; 542 return a; 527 } 543 } 528 544 >> 545 >> 546 >> 547 // 529 // CalculateExtent 548 // CalculateExtent 530 // 549 // 531 // See notes in G4VCSGface 550 // See notes in G4VCSGface 532 // 551 // 533 void G4PolyconeSide::CalculateExtent( const EA 552 void G4PolyconeSide::CalculateExtent( const EAxis axis, 534 const G4 << 553 const G4VoxelLimits &voxelLimit, 535 const G4 << 554 const G4AffineTransform &transform, 536 G4 << 555 G4SolidExtentList &extentList ) 537 { 556 { 538 G4ClippablePolygon polygon; 557 G4ClippablePolygon polygon; 539 558 540 // 559 // 541 // Here we will approximate (ala G4Cons) and 560 // Here we will approximate (ala G4Cons) and divide our conical section 542 // into segments, like G4Polyhedra. When doi 561 // into segments, like G4Polyhedra. When doing so, the radius 543 // is extented far enough such that the segm 562 // is extented far enough such that the segments always lie 544 // just outside the surface of the conical s 563 // just outside the surface of the conical section we are 545 // approximating. 564 // approximating. 546 // 565 // 547 566 548 // 567 // 549 // Choose phi size of our segment(s) based o 568 // Choose phi size of our segment(s) based on constants as 550 // defined in meshdefs.hh 569 // defined in meshdefs.hh 551 // 570 // 552 G4int numPhi = (G4int)(deltaPhi/kMeshAngleDe 571 G4int numPhi = (G4int)(deltaPhi/kMeshAngleDefault) + 1; 553 if (numPhi < kMinMeshSections) 572 if (numPhi < kMinMeshSections) 554 numPhi = kMinMeshSections; 573 numPhi = kMinMeshSections; 555 else if (numPhi > kMaxMeshSections) 574 else if (numPhi > kMaxMeshSections) 556 numPhi = kMaxMeshSections; 575 numPhi = kMaxMeshSections; 557 576 558 G4double sigPhi = deltaPhi/numPhi; 577 G4double sigPhi = deltaPhi/numPhi; 559 578 560 // 579 // 561 // Determine radius factor to keep segments 580 // Determine radius factor to keep segments outside 562 // 581 // 563 G4double rFudge = 1.0/std::cos(0.5*sigPhi); 582 G4double rFudge = 1.0/std::cos(0.5*sigPhi); 564 583 565 // 584 // 566 // Decide which radius to use on each end of 585 // Decide which radius to use on each end of the side, 567 // and whether a transition mesh is required 586 // and whether a transition mesh is required 568 // 587 // 569 // {r0,z0} - Beginning of this side 588 // {r0,z0} - Beginning of this side 570 // {r1,z1} - Ending of this side 589 // {r1,z1} - Ending of this side 571 // {r2,z0} - Beginning of transition piece 590 // {r2,z0} - Beginning of transition piece connecting previous 572 // side (and ends at beginning of 591 // side (and ends at beginning of this side) 573 // 592 // 574 // So, order is 2 --> 0 --> 1. 593 // So, order is 2 --> 0 --> 1. 575 // ------- 594 // ------- 576 // 595 // 577 // r2 < 0 indicates that no transition piece 596 // r2 < 0 indicates that no transition piece is required 578 // 597 // 579 G4double r0, r1, r2, z0, z1; 598 G4double r0, r1, r2, z0, z1; 580 599 581 r2 = -1; // By default: no transition piece 600 r2 = -1; // By default: no transition piece 582 601 583 if (rNorm < -DBL_MIN) 602 if (rNorm < -DBL_MIN) 584 { 603 { 585 // 604 // 586 // This side faces *inward*, and so our me 605 // This side faces *inward*, and so our mesh has 587 // the same radius 606 // the same radius 588 // 607 // 589 r1 = r[1]; 608 r1 = r[1]; 590 z1 = z[1]; 609 z1 = z[1]; 591 z0 = z[0]; 610 z0 = z[0]; 592 r0 = r[0]; 611 r0 = r[0]; 593 612 594 r2 = -1; 613 r2 = -1; 595 614 596 if (prevZS > DBL_MIN) 615 if (prevZS > DBL_MIN) 597 { 616 { 598 // 617 // 599 // The previous side is facing outwards 618 // The previous side is facing outwards 600 // 619 // 601 if ( prevRS*zS - prevZS*rS > 0 ) 620 if ( prevRS*zS - prevZS*rS > 0 ) 602 { 621 { 603 // 622 // 604 // Transition was convex: build transi 623 // Transition was convex: build transition piece 605 // 624 // 606 if (r[0] > DBL_MIN) r2 = r[0]*rFudge; 625 if (r[0] > DBL_MIN) r2 = r[0]*rFudge; 607 } 626 } 608 else 627 else 609 { 628 { 610 // 629 // 611 // Transition was concave: short this 630 // Transition was concave: short this side 612 // 631 // 613 FindLineIntersect( z0, r0, zS, rS, 632 FindLineIntersect( z0, r0, zS, rS, 614 z0, r0*rFudge, prev 633 z0, r0*rFudge, prevZS, prevRS*rFudge, z0, r0 ); 615 } 634 } 616 } 635 } 617 636 618 if ( nextZS > DBL_MIN && (rS*nextZS - zS*n 637 if ( nextZS > DBL_MIN && (rS*nextZS - zS*nextRS < 0) ) 619 { 638 { 620 // 639 // 621 // The next side is facing outwards, for 640 // The next side is facing outwards, forming a 622 // concave transition: short this side 641 // concave transition: short this side 623 // 642 // 624 FindLineIntersect( z1, r1, zS, rS, 643 FindLineIntersect( z1, r1, zS, rS, 625 z1, r1*rFudge, nextZS 644 z1, r1*rFudge, nextZS, nextRS*rFudge, z1, r1 ); 626 } 645 } 627 } 646 } 628 else if (rNorm > DBL_MIN) 647 else if (rNorm > DBL_MIN) 629 { 648 { 630 // 649 // 631 // This side faces *outward* and is given 650 // This side faces *outward* and is given a boost to 632 // it radius 651 // it radius 633 // 652 // 634 r0 = r[0]*rFudge; 653 r0 = r[0]*rFudge; 635 z0 = z[0]; 654 z0 = z[0]; 636 r1 = r[1]*rFudge; 655 r1 = r[1]*rFudge; 637 z1 = z[1]; 656 z1 = z[1]; 638 657 639 if (prevZS < -DBL_MIN) 658 if (prevZS < -DBL_MIN) 640 { 659 { 641 // 660 // 642 // The previous side is facing inwards 661 // The previous side is facing inwards 643 // 662 // 644 if ( prevRS*zS - prevZS*rS > 0 ) 663 if ( prevRS*zS - prevZS*rS > 0 ) 645 { 664 { 646 // 665 // 647 // Transition was convex: build transi 666 // Transition was convex: build transition piece 648 // 667 // 649 if (r[0] > DBL_MIN) r2 = r[0]; 668 if (r[0] > DBL_MIN) r2 = r[0]; 650 } 669 } 651 else 670 else 652 { 671 { 653 // 672 // 654 // Transition was concave: short this 673 // Transition was concave: short this side 655 // 674 // 656 FindLineIntersect( z0, r0, zS, rS*rFud 675 FindLineIntersect( z0, r0, zS, rS*rFudge, 657 z0, r[0], prevZS, p 676 z0, r[0], prevZS, prevRS, z0, r0 ); 658 } 677 } 659 } 678 } 660 679 661 if ( nextZS < -DBL_MIN && (rS*nextZS - zS* 680 if ( nextZS < -DBL_MIN && (rS*nextZS - zS*nextRS < 0) ) 662 { 681 { 663 // 682 // 664 // The next side is facing inwards, form 683 // The next side is facing inwards, forming a 665 // concave transition: short this side 684 // concave transition: short this side 666 // 685 // 667 FindLineIntersect( z1, r1, zS, rS*rFudge 686 FindLineIntersect( z1, r1, zS, rS*rFudge, 668 z1, r[1], nextZS, nex 687 z1, r[1], nextZS, nextRS, z1, r1 ); 669 } 688 } 670 } 689 } 671 else 690 else 672 { 691 { 673 // 692 // 674 // This side is perpendicular to the z axi 693 // This side is perpendicular to the z axis (is a disk) 675 // 694 // 676 // Whether or not r0 needs a rFudge factor 695 // Whether or not r0 needs a rFudge factor depends 677 // on the normal of the previous edge. Sim 696 // on the normal of the previous edge. Similar with r1 678 // and the next edge. No transition piece 697 // and the next edge. No transition piece is required. 679 // 698 // 680 r0 = r[0]; 699 r0 = r[0]; 681 r1 = r[1]; 700 r1 = r[1]; 682 z0 = z[0]; 701 z0 = z[0]; 683 z1 = z[1]; 702 z1 = z[1]; 684 703 685 if (prevZS > DBL_MIN) r0 *= rFudge; 704 if (prevZS > DBL_MIN) r0 *= rFudge; 686 if (nextZS > DBL_MIN) r1 *= rFudge; 705 if (nextZS > DBL_MIN) r1 *= rFudge; 687 } 706 } 688 707 689 // 708 // 690 // Loop 709 // Loop 691 // 710 // 692 G4double phi = startPhi, 711 G4double phi = startPhi, 693 cosPhi = std::cos(phi), 712 cosPhi = std::cos(phi), 694 sinPhi = std::sin(phi); 713 sinPhi = std::sin(phi); 695 714 696 G4ThreeVector v0( r0*cosPhi, r0*sinPhi, z0 ) 715 G4ThreeVector v0( r0*cosPhi, r0*sinPhi, z0 ), 697 v1( r1*cosPhi, r1*sinPhi, 716 v1( r1*cosPhi, r1*sinPhi, z1 ), 698 v2, w0, w1, w2; 717 v2, w0, w1, w2; 699 transform.ApplyPointTransform( v0 ); 718 transform.ApplyPointTransform( v0 ); 700 transform.ApplyPointTransform( v1 ); 719 transform.ApplyPointTransform( v1 ); 701 720 702 if (r2 >= 0) 721 if (r2 >= 0) 703 { 722 { 704 v2 = G4ThreeVector( r2*cosPhi, r2*sinPhi, 723 v2 = G4ThreeVector( r2*cosPhi, r2*sinPhi, z0 ); 705 transform.ApplyPointTransform( v2 ); 724 transform.ApplyPointTransform( v2 ); 706 } 725 } 707 726 708 do // Loop checking, 13.08.2015, G.Cosmo << 727 do 709 { 728 { 710 phi += sigPhi; 729 phi += sigPhi; 711 if (numPhi == 1) phi = startPhi+deltaPhi; 730 if (numPhi == 1) phi = startPhi+deltaPhi; // Try to avoid roundoff 712 cosPhi = std::cos(phi), 731 cosPhi = std::cos(phi), 713 sinPhi = std::sin(phi); 732 sinPhi = std::sin(phi); 714 733 715 w0 = G4ThreeVector( r0*cosPhi, r0*sinPhi, 734 w0 = G4ThreeVector( r0*cosPhi, r0*sinPhi, z0 ); 716 w1 = G4ThreeVector( r1*cosPhi, r1*sinPhi, 735 w1 = G4ThreeVector( r1*cosPhi, r1*sinPhi, z1 ); 717 transform.ApplyPointTransform( w0 ); 736 transform.ApplyPointTransform( w0 ); 718 transform.ApplyPointTransform( w1 ); 737 transform.ApplyPointTransform( w1 ); 719 738 720 G4ThreeVector deltaV = r0 > r1 ? w0-v0 : w 739 G4ThreeVector deltaV = r0 > r1 ? w0-v0 : w1-v1; 721 740 722 // 741 // 723 // Build polygon, taking special care to k 742 // Build polygon, taking special care to keep the vertices 724 // in order 743 // in order 725 // 744 // 726 polygon.ClearAllVertices(); 745 polygon.ClearAllVertices(); 727 746 728 polygon.AddVertexInOrder( v0 ); 747 polygon.AddVertexInOrder( v0 ); 729 polygon.AddVertexInOrder( v1 ); 748 polygon.AddVertexInOrder( v1 ); 730 polygon.AddVertexInOrder( w1 ); 749 polygon.AddVertexInOrder( w1 ); 731 polygon.AddVertexInOrder( w0 ); 750 polygon.AddVertexInOrder( w0 ); 732 751 733 // 752 // 734 // Get extent 753 // Get extent 735 // 754 // 736 if (polygon.PartialClip( voxelLimit, axis 755 if (polygon.PartialClip( voxelLimit, axis )) 737 { 756 { 738 // 757 // 739 // Get dot product of normal with target 758 // Get dot product of normal with target axis 740 // 759 // 741 polygon.SetNormal( deltaV.cross(v1-v0).u 760 polygon.SetNormal( deltaV.cross(v1-v0).unit() ); 742 761 743 extentList.AddSurface( polygon ); 762 extentList.AddSurface( polygon ); 744 } 763 } 745 764 746 if (r2 >= 0) 765 if (r2 >= 0) 747 { 766 { 748 // 767 // 749 // Repeat, for transition piece 768 // Repeat, for transition piece 750 // 769 // 751 w2 = G4ThreeVector( r2*cosPhi, r2*sinPhi 770 w2 = G4ThreeVector( r2*cosPhi, r2*sinPhi, z0 ); 752 transform.ApplyPointTransform( w2 ); 771 transform.ApplyPointTransform( w2 ); 753 772 754 polygon.ClearAllVertices(); 773 polygon.ClearAllVertices(); 755 774 756 polygon.AddVertexInOrder( v2 ); 775 polygon.AddVertexInOrder( v2 ); 757 polygon.AddVertexInOrder( v0 ); 776 polygon.AddVertexInOrder( v0 ); 758 polygon.AddVertexInOrder( w0 ); 777 polygon.AddVertexInOrder( w0 ); 759 polygon.AddVertexInOrder( w2 ); 778 polygon.AddVertexInOrder( w2 ); 760 779 761 if (polygon.PartialClip( voxelLimit, axi 780 if (polygon.PartialClip( voxelLimit, axis )) 762 { 781 { 763 polygon.SetNormal( deltaV.cross(v0-v2) 782 polygon.SetNormal( deltaV.cross(v0-v2).unit() ); 764 783 765 extentList.AddSurface( polygon ); 784 extentList.AddSurface( polygon ); 766 } 785 } 767 786 768 v2 = w2; 787 v2 = w2; 769 } 788 } 770 789 771 // 790 // 772 // Next vertex 791 // Next vertex 773 // 792 // 774 v0 = w0; 793 v0 = w0; 775 v1 = w1; 794 v1 = w1; 776 } while( --numPhi > 0 ); 795 } while( --numPhi > 0 ); 777 796 778 // 797 // 779 // We are almost done. But, it is important 798 // We are almost done. But, it is important that we leave no 780 // gaps in the surface of our solid. By usin 799 // gaps in the surface of our solid. By using rFudge, however, 781 // we've done exactly that, if we have a phi 800 // we've done exactly that, if we have a phi segment. 782 // Add two additional faces if necessary 801 // Add two additional faces if necessary 783 // 802 // 784 if (phiIsOpen && rNorm > DBL_MIN) 803 if (phiIsOpen && rNorm > DBL_MIN) 785 { 804 { 786 cosPhi = std::cos(startPhi); << 805 G4double cosPhi = std::cos(startPhi), 787 sinPhi = std::sin(startPhi); << 806 sinPhi = std::sin(startPhi); 788 807 789 G4ThreeVector a0( r[0]*cosPhi, r[0]*sinPhi 808 G4ThreeVector a0( r[0]*cosPhi, r[0]*sinPhi, z[0] ), 790 a1( r[1]*cosPhi, r[1]*sinPhi 809 a1( r[1]*cosPhi, r[1]*sinPhi, z[1] ), 791 b0( r0*cosPhi, r0*sinPhi, z[ 810 b0( r0*cosPhi, r0*sinPhi, z[0] ), 792 b1( r1*cosPhi, r1*sinPhi, z[ 811 b1( r1*cosPhi, r1*sinPhi, z[1] ); 793 812 794 transform.ApplyPointTransform( a0 ); 813 transform.ApplyPointTransform( a0 ); 795 transform.ApplyPointTransform( a1 ); 814 transform.ApplyPointTransform( a1 ); 796 transform.ApplyPointTransform( b0 ); 815 transform.ApplyPointTransform( b0 ); 797 transform.ApplyPointTransform( b1 ); 816 transform.ApplyPointTransform( b1 ); 798 817 799 polygon.ClearAllVertices(); 818 polygon.ClearAllVertices(); 800 819 801 polygon.AddVertexInOrder( a0 ); 820 polygon.AddVertexInOrder( a0 ); 802 polygon.AddVertexInOrder( a1 ); 821 polygon.AddVertexInOrder( a1 ); 803 polygon.AddVertexInOrder( b0 ); 822 polygon.AddVertexInOrder( b0 ); 804 polygon.AddVertexInOrder( b1 ); 823 polygon.AddVertexInOrder( b1 ); 805 824 806 if (polygon.PartialClip( voxelLimit , axis 825 if (polygon.PartialClip( voxelLimit , axis)) 807 { 826 { 808 G4ThreeVector normal( sinPhi, -cosPhi, 0 827 G4ThreeVector normal( sinPhi, -cosPhi, 0 ); 809 polygon.SetNormal( transform.TransformAx 828 polygon.SetNormal( transform.TransformAxis( normal ) ); 810 829 811 extentList.AddSurface( polygon ); 830 extentList.AddSurface( polygon ); 812 } 831 } 813 832 814 cosPhi = std::cos(startPhi+deltaPhi); 833 cosPhi = std::cos(startPhi+deltaPhi); 815 sinPhi = std::sin(startPhi+deltaPhi); 834 sinPhi = std::sin(startPhi+deltaPhi); 816 835 817 a0 = G4ThreeVector( r[0]*cosPhi, r[0]*sinP 836 a0 = G4ThreeVector( r[0]*cosPhi, r[0]*sinPhi, z[0] ), 818 a1 = G4ThreeVector( r[1]*cosPhi, r[1]*sinP 837 a1 = G4ThreeVector( r[1]*cosPhi, r[1]*sinPhi, z[1] ), 819 b0 = G4ThreeVector( r0*cosPhi, r0*sinPhi, 838 b0 = G4ThreeVector( r0*cosPhi, r0*sinPhi, z[0] ), 820 b1 = G4ThreeVector( r1*cosPhi, r1*sinPhi, 839 b1 = G4ThreeVector( r1*cosPhi, r1*sinPhi, z[1] ); 821 transform.ApplyPointTransform( a0 ); 840 transform.ApplyPointTransform( a0 ); 822 transform.ApplyPointTransform( a1 ); 841 transform.ApplyPointTransform( a1 ); 823 transform.ApplyPointTransform( b0 ); 842 transform.ApplyPointTransform( b0 ); 824 transform.ApplyPointTransform( b1 ); 843 transform.ApplyPointTransform( b1 ); 825 844 826 polygon.ClearAllVertices(); 845 polygon.ClearAllVertices(); 827 846 828 polygon.AddVertexInOrder( a0 ); 847 polygon.AddVertexInOrder( a0 ); 829 polygon.AddVertexInOrder( a1 ); 848 polygon.AddVertexInOrder( a1 ); 830 polygon.AddVertexInOrder( b0 ); 849 polygon.AddVertexInOrder( b0 ); 831 polygon.AddVertexInOrder( b1 ); 850 polygon.AddVertexInOrder( b1 ); 832 851 833 if (polygon.PartialClip( voxelLimit, axis 852 if (polygon.PartialClip( voxelLimit, axis )) 834 { 853 { 835 G4ThreeVector normal( -sinPhi, cosPhi, 0 854 G4ThreeVector normal( -sinPhi, cosPhi, 0 ); 836 polygon.SetNormal( transform.TransformAx 855 polygon.SetNormal( transform.TransformAxis( normal ) ); 837 856 838 extentList.AddSurface( polygon ); 857 extentList.AddSurface( polygon ); 839 } 858 } 840 } 859 } 841 860 842 return; 861 return; 843 } 862 } 844 863 >> 864 // 845 // GetPhi 865 // GetPhi 846 // 866 // 847 // Calculate Phi for a given 3-vector (point), 867 // Calculate Phi for a given 3-vector (point), if not already cached for the 848 // same point, in the attempt to avoid consecu 868 // same point, in the attempt to avoid consecutive computation of the same 849 // quantity 869 // quantity 850 // 870 // 851 G4double G4PolyconeSide::GetPhi( const G4Three 871 G4double G4PolyconeSide::GetPhi( const G4ThreeVector& p ) 852 { 872 { 853 G4double val=0.; 873 G4double val=0.; 854 G4ThreeVector vphi(G4MT_pcphix, G4MT_pcphiy, << 855 874 856 if (vphi != p) << 875 if (fPhi.first != p) 857 { 876 { 858 val = p.phi(); 877 val = p.phi(); 859 G4MT_pcphix = p.x(); G4MT_pcphiy = p.y(); << 878 fPhi.first = p; 860 G4MT_pcphik = val; << 879 fPhi.second = val; 861 } 880 } 862 else 881 else 863 { 882 { 864 val = G4MT_pcphik; << 883 val = fPhi.second; 865 } 884 } 866 return val; 885 return val; 867 } 886 } 868 887 >> 888 // 869 // DistanceAway 889 // DistanceAway 870 // 890 // 871 // Calculate distance of a point from our coni 891 // Calculate distance of a point from our conical surface, including the effect 872 // of any phi segmentation 892 // of any phi segmentation 873 // 893 // 874 // Arguments: 894 // Arguments: 875 // p - (in) Point to check 895 // p - (in) Point to check 876 // opposite - (in) If true, check opposi 896 // opposite - (in) If true, check opposite hemisphere (see below) 877 // distOutside - (out) Additional distance 897 // distOutside - (out) Additional distance outside the edges of the surface 878 // edgeRZnorm - (out) if negative, point i 898 // edgeRZnorm - (out) if negative, point is inside 879 // 899 // 880 // return value = distance from the conical p 900 // return value = distance from the conical plane, if extrapolated beyond edges, 881 // signed by whether the point 901 // signed by whether the point is in inside or outside the shape 882 // 902 // 883 // Notes: 903 // Notes: 884 // * There are two answers, depending on whic 904 // * There are two answers, depending on which hemisphere is considered. 885 // 905 // 886 G4double G4PolyconeSide::DistanceAway( const G << 906 G4double G4PolyconeSide::DistanceAway( const G4ThreeVector &p, 887 G 907 G4bool opposite, 888 G << 908 G4double &distOutside2, 889 G << 909 G4double *edgeRZnorm ) 890 { 910 { 891 // 911 // 892 // Convert our point to r and z 912 // Convert our point to r and z 893 // 913 // 894 G4double rx = p.perp(), zx = p.z(); 914 G4double rx = p.perp(), zx = p.z(); 895 915 896 // 916 // 897 // Change sign of r if opposite says we shou 917 // Change sign of r if opposite says we should 898 // 918 // 899 if (opposite) rx = -rx; 919 if (opposite) rx = -rx; 900 920 901 // 921 // 902 // Calculate return value 922 // Calculate return value 903 // 923 // 904 G4double deltaR = rx - r[0], deltaZ = zx - 924 G4double deltaR = rx - r[0], deltaZ = zx - z[0]; 905 G4double answer = deltaR*rNorm + deltaZ*zNor 925 G4double answer = deltaR*rNorm + deltaZ*zNorm; 906 926 907 // 927 // 908 // Are we off the surface in r,z space? 928 // Are we off the surface in r,z space? 909 // 929 // 910 G4double q = deltaR*rS + deltaZ*zS; << 930 G4double s = deltaR*rS + deltaZ*zS; 911 if (q < 0) << 931 if (s < 0) 912 { 932 { 913 distOutside2 = q*q; << 933 distOutside2 = s*s; 914 if (edgeRZnorm != nullptr) << 934 if (edgeRZnorm) *edgeRZnorm = deltaR*rNormEdge[0] + deltaZ*zNormEdge[0]; 915 *edgeRZnorm = deltaR*rNormEdge[0] + delt << 916 } 935 } 917 else if (q > length) << 936 else if (s > length) 918 { 937 { 919 distOutside2 = sqr( q-length ); << 938 distOutside2 = sqr( s-length ); 920 if (edgeRZnorm != nullptr) << 939 if (edgeRZnorm) 921 { 940 { 922 deltaR = rx - r[1]; << 941 G4double deltaR = rx - r[1], deltaZ = zx - z[1]; 923 deltaZ = zx - z[1]; << 924 *edgeRZnorm = deltaR*rNormEdge[1] + delt 942 *edgeRZnorm = deltaR*rNormEdge[1] + deltaZ*zNormEdge[1]; 925 } 943 } 926 } 944 } 927 else 945 else 928 { 946 { 929 distOutside2 = 0.; << 947 distOutside2 = 0; 930 if (edgeRZnorm != nullptr) *edgeRZnorm = a << 948 if (edgeRZnorm) *edgeRZnorm = answer; 931 } 949 } 932 950 933 if (phiIsOpen) 951 if (phiIsOpen) 934 { 952 { 935 // 953 // 936 // Finally, check phi 954 // Finally, check phi 937 // 955 // 938 G4double phi = GetPhi(p); 956 G4double phi = GetPhi(p); 939 while( phi < startPhi ) // Loop checkin << 957 while( phi < startPhi ) phi += twopi; 940 phi += twopi; << 941 958 942 if (phi > startPhi+deltaPhi) 959 if (phi > startPhi+deltaPhi) 943 { 960 { 944 // 961 // 945 // Oops. Are we closer to the start phi 962 // Oops. Are we closer to the start phi or end phi? 946 // 963 // 947 G4double d1 = phi-startPhi-deltaPhi; 964 G4double d1 = phi-startPhi-deltaPhi; 948 while( phi > startPhi ) // Loop check << 965 while( phi > startPhi ) phi -= twopi; 949 phi -= twopi; << 950 G4double d2 = startPhi-phi; 966 G4double d2 = startPhi-phi; 951 967 952 if (d2 < d1) d1 = d2; 968 if (d2 < d1) d1 = d2; 953 969 954 // 970 // 955 // Add result to our distance 971 // Add result to our distance 956 // 972 // 957 G4double dist = d1*rx; 973 G4double dist = d1*rx; 958 974 959 distOutside2 += dist*dist; 975 distOutside2 += dist*dist; 960 if (edgeRZnorm != nullptr) << 976 if (edgeRZnorm) 961 { 977 { 962 *edgeRZnorm = std::max(std::fabs(*edge 978 *edgeRZnorm = std::max(std::fabs(*edgeRZnorm),std::fabs(dist)); 963 } 979 } 964 } 980 } 965 } 981 } 966 982 967 return answer; 983 return answer; 968 } 984 } 969 985 970 // DistanceAway << 971 // << 972 // Special version of DistanceAway for Inside. << 973 // Opposite parameter is not used, instead use << 974 // << 975 G4double G4PolyconeSide::DistanceAway( const G << 976 G << 977 G << 978 { << 979 // << 980 // Convert our point to r and z << 981 // << 982 G4double rx = p.perp(), zx = p.z(); << 983 << 984 // << 985 // Change sign of r if we should << 986 // << 987 G4int part = 1; << 988 if (rx < 0) part = -1; << 989 << 990 // << 991 // Calculate return value << 992 // << 993 G4double deltaR = rx - r[0]*part, deltaZ = z << 994 G4double answer = deltaR*rNorm*part + deltaZ << 995 << 996 // << 997 // Are we off the surface in r,z space? << 998 // << 999 G4double q = deltaR*rS*part + deltaZ*zS; << 1000 if (q < 0) << 1001 { << 1002 distOutside2 = q*q; << 1003 if (edgeRZnorm != nullptr) << 1004 { << 1005 *edgeRZnorm = deltaR*rNormEdge[0]*part << 1006 } << 1007 } << 1008 else if (q > length) << 1009 { << 1010 distOutside2 = sqr( q-length ); << 1011 if (edgeRZnorm != nullptr) << 1012 { << 1013 deltaR = rx - r[1]*part; << 1014 deltaZ = zx - z[1]; << 1015 *edgeRZnorm = deltaR*rNormEdge[1]*part << 1016 } << 1017 } << 1018 else << 1019 { << 1020 distOutside2 = 0.; << 1021 if (edgeRZnorm != nullptr) *edgeRZnorm = << 1022 } << 1023 << 1024 if (phiIsOpen) << 1025 { << 1026 // << 1027 // Finally, check phi << 1028 // << 1029 G4double phi = GetPhi(p); << 1030 while( phi < startPhi ) // Loop checki << 1031 phi += twopi; << 1032 << 1033 if (phi > startPhi+deltaPhi) << 1034 { << 1035 // << 1036 // Oops. Are we closer to the start phi << 1037 // << 1038 G4double d1 = phi-startPhi-deltaPhi; << 1039 while( phi > startPhi ) // Loop chec << 1040 phi -= twopi; << 1041 G4double d2 = startPhi-phi; << 1042 << 1043 if (d2 < d1) d1 = d2; << 1044 << 1045 // << 1046 // Add result to our distance << 1047 // << 1048 G4double dist = d1*rx*part; << 1049 << 1050 distOutside2 += dist*dist; << 1051 if (edgeRZnorm != nullptr) << 1052 { << 1053 *edgeRZnorm = std::max(std::fabs(*edg << 1054 } << 1055 } << 1056 } << 1057 << 1058 return answer; << 1059 } << 1060 986 >> 987 // 1061 // PointOnCone 988 // PointOnCone 1062 // 989 // 1063 // Decide if a point is on a cone and return 990 // Decide if a point is on a cone and return normal if it is 1064 // 991 // 1065 G4bool G4PolyconeSide::PointOnCone( const G4T << 992 G4bool G4PolyconeSide::PointOnCone( const G4ThreeVector &hit, 1066 G4d 993 G4double normSign, 1067 const G4T << 994 const G4ThreeVector &p, 1068 const G4T << 995 const G4ThreeVector &v, 1069 G4T << 996 G4ThreeVector &normal ) 1070 { 997 { 1071 G4double rx = hit.perp(); 998 G4double rx = hit.perp(); 1072 // 999 // 1073 // Check radial/z extent, as appropriate 1000 // Check radial/z extent, as appropriate 1074 // 1001 // 1075 if (!cone->HitOn( rx, hit.z() )) return fal 1002 if (!cone->HitOn( rx, hit.z() )) return false; 1076 1003 1077 if (phiIsOpen) 1004 if (phiIsOpen) 1078 { 1005 { 1079 G4double phiTolerant = 2.0*kCarTolerance/ 1006 G4double phiTolerant = 2.0*kCarTolerance/(rx+kCarTolerance); 1080 // 1007 // 1081 // Check phi segment. Here we have to be 1008 // Check phi segment. Here we have to be careful 1082 // to use the standard method consistent 1009 // to use the standard method consistent with 1083 // PolyPhiFace. See PolyPhiFace::InsideEd 1010 // PolyPhiFace. See PolyPhiFace::InsideEdgesExact 1084 // 1011 // 1085 G4double phi = GetPhi(hit); 1012 G4double phi = GetPhi(hit); 1086 while( phi < startPhi-phiTolerant ) // << 1013 while( phi < startPhi-phiTolerant ) phi += twopi; 1087 phi += twopi; << 1088 1014 1089 if (phi > startPhi+deltaPhi+phiTolerant) 1015 if (phi > startPhi+deltaPhi+phiTolerant) return false; 1090 1016 1091 if (phi > startPhi+deltaPhi-phiTolerant) 1017 if (phi > startPhi+deltaPhi-phiTolerant) 1092 { 1018 { 1093 // 1019 // 1094 // Exact treatment 1020 // Exact treatment 1095 // 1021 // 1096 G4ThreeVector qx = p + v; 1022 G4ThreeVector qx = p + v; 1097 G4ThreeVector qa = qx - corners[2], 1023 G4ThreeVector qa = qx - corners[2], 1098 qb = qx - corners[3]; 1024 qb = qx - corners[3]; 1099 G4ThreeVector qacb = qa.cross(qb); 1025 G4ThreeVector qacb = qa.cross(qb); 1100 1026 1101 if (normSign*qacb.dot(v) < 0) return fa 1027 if (normSign*qacb.dot(v) < 0) return false; 1102 } 1028 } 1103 else if (phi < phiTolerant) 1029 else if (phi < phiTolerant) 1104 { 1030 { 1105 G4ThreeVector qx = p + v; 1031 G4ThreeVector qx = p + v; 1106 G4ThreeVector qa = qx - corners[1], 1032 G4ThreeVector qa = qx - corners[1], 1107 qb = qx - corners[0]; 1033 qb = qx - corners[0]; 1108 G4ThreeVector qacb = qa.cross(qb); 1034 G4ThreeVector qacb = qa.cross(qb); 1109 1035 1110 if (normSign*qacb.dot(v) < 0) return fa 1036 if (normSign*qacb.dot(v) < 0) return false; 1111 } 1037 } 1112 } 1038 } 1113 1039 1114 // 1040 // 1115 // We have a good hit! Calculate normal 1041 // We have a good hit! Calculate normal 1116 // 1042 // 1117 if (rx < DBL_MIN) 1043 if (rx < DBL_MIN) 1118 normal = G4ThreeVector( 0, 0, zNorm < 0 ? 1044 normal = G4ThreeVector( 0, 0, zNorm < 0 ? -1 : 1 ); 1119 else 1045 else 1120 normal = G4ThreeVector( rNorm*hit.x()/rx, 1046 normal = G4ThreeVector( rNorm*hit.x()/rx, rNorm*hit.y()/rx, zNorm ); 1121 return true; 1047 return true; 1122 } 1048 } 1123 1049 >> 1050 >> 1051 // 1124 // FindLineIntersect 1052 // FindLineIntersect 1125 // 1053 // 1126 // Decide the point at which two 2-dimensiona 1054 // Decide the point at which two 2-dimensional lines intersect 1127 // 1055 // 1128 // Equation of line: x = x1 + s*tx1 1056 // Equation of line: x = x1 + s*tx1 1129 // y = y1 + s*ty1 1057 // y = y1 + s*ty1 1130 // 1058 // 1131 // It is assumed that the lines are *not* par 1059 // It is assumed that the lines are *not* parallel 1132 // 1060 // 1133 void G4PolyconeSide::FindLineIntersect( G4dou 1061 void G4PolyconeSide::FindLineIntersect( G4double x1, G4double y1, 1134 G4dou 1062 G4double tx1, G4double ty1, 1135 G4dou 1063 G4double x2, G4double y2, 1136 G4dou 1064 G4double tx2, G4double ty2, 1137 G4dou << 1065 G4double &x, G4double &y ) 1138 { 1066 { 1139 // 1067 // 1140 // The solution is a simple linear equation 1068 // The solution is a simple linear equation 1141 // 1069 // 1142 G4double deter = tx1*ty2 - tx2*ty1; 1070 G4double deter = tx1*ty2 - tx2*ty1; 1143 1071 1144 G4double s1 = ((x2-x1)*ty2 - tx2*(y2-y1))/d 1072 G4double s1 = ((x2-x1)*ty2 - tx2*(y2-y1))/deter; 1145 G4double s2 = ((x2-x1)*ty1 - tx1*(y2-y1))/d 1073 G4double s2 = ((x2-x1)*ty1 - tx1*(y2-y1))/deter; 1146 1074 1147 // 1075 // 1148 // We want the answer to not depend on whic 1076 // We want the answer to not depend on which order the 1149 // lines were specified. Take average. 1077 // lines were specified. Take average. 1150 // 1078 // 1151 x = 0.5*( x1+s1*tx1 + x2+s2*tx2 ); 1079 x = 0.5*( x1+s1*tx1 + x2+s2*tx2 ); 1152 y = 0.5*( y1+s1*ty1 + y2+s2*ty2 ); 1080 y = 0.5*( y1+s1*ty1 + y2+s2*ty2 ); 1153 } 1081 } 1154 1082 >> 1083 // 1155 // Calculate surface area for GetPointOnSurfa 1084 // Calculate surface area for GetPointOnSurface() 1156 // 1085 // 1157 G4double G4PolyconeSide::SurfaceArea() 1086 G4double G4PolyconeSide::SurfaceArea() 1158 { 1087 { 1159 if(fSurfaceArea==0.) << 1088 if(fSurfaceArea==0) 1160 { 1089 { 1161 fSurfaceArea = (r[0]+r[1])* std::sqrt(sqr 1090 fSurfaceArea = (r[0]+r[1])* std::sqrt(sqr(r[0]-r[1])+sqr(z[0]-z[1])); 1162 fSurfaceArea *= 0.5*(deltaPhi); 1091 fSurfaceArea *= 0.5*(deltaPhi); 1163 } 1092 } 1164 return fSurfaceArea; 1093 return fSurfaceArea; 1165 } 1094 } 1166 1095 >> 1096 // 1167 // GetPointOnFace 1097 // GetPointOnFace 1168 // 1098 // 1169 G4ThreeVector G4PolyconeSide::GetPointOnFace( 1099 G4ThreeVector G4PolyconeSide::GetPointOnFace() 1170 { 1100 { 1171 G4double x,y,zz; 1101 G4double x,y,zz; 1172 G4double rr,phi,dz,dr; 1102 G4double rr,phi,dz,dr; 1173 dr=r[1]-r[0];dz=z[1]-z[0]; 1103 dr=r[1]-r[0];dz=z[1]-z[0]; 1174 phi=startPhi+deltaPhi*G4UniformRand(); 1104 phi=startPhi+deltaPhi*G4UniformRand(); 1175 rr=r[0]+dr*G4UniformRand(); 1105 rr=r[0]+dr*G4UniformRand(); 1176 1106 1177 x=rr*std::cos(phi); 1107 x=rr*std::cos(phi); 1178 y=rr*std::sin(phi); 1108 y=rr*std::sin(phi); 1179 1109 1180 // PolyconeSide has a Ring Form 1110 // PolyconeSide has a Ring Form 1181 // 1111 // 1182 if (dz==0.) 1112 if (dz==0.) 1183 { 1113 { 1184 zz=z[0]; 1114 zz=z[0]; 1185 } 1115 } 1186 else 1116 else 1187 { 1117 { 1188 if(dr==0.) // PolyconeSide has a Tube Fo 1118 if(dr==0.) // PolyconeSide has a Tube Form 1189 { 1119 { 1190 zz = z[0]+dz*G4UniformRand(); 1120 zz = z[0]+dz*G4UniformRand(); 1191 } 1121 } 1192 else 1122 else 1193 { 1123 { 1194 zz = z[0]+(rr-r[0])*dz/dr; 1124 zz = z[0]+(rr-r[0])*dz/dr; 1195 } 1125 } 1196 } 1126 } 1197 1127 1198 return {x,y,zz}; << 1128 return G4ThreeVector(x,y,zz); 1199 } 1129 } 1200 1130