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