[831] | 1 | // |
---|
| 2 | // ******************************************************************** |
---|
| 3 | // * License and Disclaimer * |
---|
| 4 | // * * |
---|
| 5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
| 6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
| 7 | // * conditions of the Geant4 Software License, included in the file * |
---|
| 8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
| 9 | // * include a list of copyright holders. * |
---|
| 10 | // * * |
---|
| 11 | // * Neither the authors of this software system, nor their employing * |
---|
| 12 | // * institutes,nor the agencies providing financial support for this * |
---|
| 13 | // * work make any representation or warranty, express or implied, * |
---|
| 14 | // * regarding this software system or assume any liability for its * |
---|
| 15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
| 16 | // * for the full disclaimer and the limitation of liability. * |
---|
| 17 | // * * |
---|
| 18 | // * This code implementation is the result of the scientific and * |
---|
| 19 | // * technical work of the GEANT4 collaboration. * |
---|
| 20 | // * By using, copying, modifying or distributing the software (or * |
---|
| 21 | // * any work based on the software) you agree to acknowledge its * |
---|
| 22 | // * use in resulting scientific publications, and indicate your * |
---|
| 23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
| 24 | // ******************************************************************** |
---|
| 25 | // |
---|
| 26 | // |
---|
[850] | 27 | // $Id: G4PolyhedraSide.cc,v 1.15 2008/05/15 11:41:59 gcosmo Exp $ |
---|
[1058] | 28 | // GEANT4 tag $Name: geant4-09-02-ref-02 $ |
---|
[831] | 29 | // |
---|
| 30 | // |
---|
| 31 | // -------------------------------------------------------------------- |
---|
| 32 | // GEANT 4 class source file |
---|
| 33 | // |
---|
| 34 | // |
---|
| 35 | // G4PolyhedraSide.cc |
---|
| 36 | // |
---|
[850] | 37 | // Implementation of the face representing one segmented side of a Polyhedra |
---|
[831] | 38 | // |
---|
| 39 | // -------------------------------------------------------------------- |
---|
| 40 | |
---|
| 41 | #include "G4PolyhedraSide.hh" |
---|
| 42 | #include "G4IntersectingCone.hh" |
---|
| 43 | #include "G4ClippablePolygon.hh" |
---|
| 44 | #include "G4AffineTransform.hh" |
---|
| 45 | #include "G4SolidExtentList.hh" |
---|
| 46 | #include "G4GeometryTolerance.hh" |
---|
| 47 | |
---|
[850] | 48 | #include "Randomize.hh" |
---|
| 49 | |
---|
[831] | 50 | // |
---|
| 51 | // Constructor |
---|
| 52 | // |
---|
| 53 | // Values for r1,z1 and r2,z2 should be specified in clockwise |
---|
| 54 | // order in (r,z). |
---|
| 55 | // |
---|
| 56 | G4PolyhedraSide::G4PolyhedraSide( const G4PolyhedraSideRZ *prevRZ, |
---|
| 57 | const G4PolyhedraSideRZ *tail, |
---|
| 58 | const G4PolyhedraSideRZ *head, |
---|
| 59 | const G4PolyhedraSideRZ *nextRZ, |
---|
| 60 | G4int theNumSide, |
---|
| 61 | G4double thePhiStart, |
---|
| 62 | G4double thePhiTotal, |
---|
| 63 | G4bool thePhiIsOpen, |
---|
| 64 | G4bool isAllBehind ) |
---|
| 65 | { |
---|
| 66 | |
---|
| 67 | kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance(); |
---|
[850] | 68 | fSurfaceArea=0.; |
---|
[831] | 69 | // |
---|
| 70 | // Record values |
---|
| 71 | // |
---|
| 72 | r[0] = tail->r; z[0] = tail->z; |
---|
| 73 | r[1] = head->r; z[1] = head->z; |
---|
| 74 | |
---|
| 75 | G4double phiTotal; |
---|
| 76 | |
---|
| 77 | // |
---|
| 78 | // Set phi to our convention |
---|
| 79 | // |
---|
| 80 | startPhi = thePhiStart; |
---|
| 81 | while (startPhi < 0.0) startPhi += twopi; |
---|
| 82 | |
---|
| 83 | phiIsOpen = thePhiIsOpen; |
---|
| 84 | phiTotal = (phiIsOpen) ? thePhiTotal : twopi; |
---|
| 85 | |
---|
| 86 | allBehind = isAllBehind; |
---|
| 87 | |
---|
| 88 | // |
---|
| 89 | // Make our intersecting cone |
---|
| 90 | // |
---|
| 91 | cone = new G4IntersectingCone( r, z ); |
---|
| 92 | |
---|
| 93 | // |
---|
| 94 | // Construct side plane vector set |
---|
| 95 | // |
---|
| 96 | numSide = theNumSide; |
---|
| 97 | deltaPhi = phiTotal/theNumSide; |
---|
| 98 | endPhi = startPhi+phiTotal; |
---|
| 99 | |
---|
| 100 | vecs = new G4PolyhedraSideVec[numSide]; |
---|
| 101 | |
---|
| 102 | edges = new G4PolyhedraSideEdge[phiIsOpen ? numSide+1 : numSide]; |
---|
| 103 | |
---|
| 104 | // |
---|
| 105 | // ...this is where we start |
---|
| 106 | // |
---|
| 107 | G4double phi = startPhi; |
---|
| 108 | G4ThreeVector a1( r[0]*std::cos(phi), r[0]*std::sin(phi), z[0] ), |
---|
| 109 | b1( r[1]*std::cos(phi), r[1]*std::sin(phi), z[1] ), |
---|
| 110 | c1( prevRZ->r*std::cos(phi), prevRZ->r*std::sin(phi), prevRZ->z ), |
---|
| 111 | d1( nextRZ->r*std::cos(phi), nextRZ->r*std::sin(phi), nextRZ->z ), |
---|
| 112 | a2, b2, c2, d2; |
---|
| 113 | G4PolyhedraSideEdge *edge = edges; |
---|
| 114 | |
---|
| 115 | G4PolyhedraSideVec *vec = vecs; |
---|
| 116 | do |
---|
| 117 | { |
---|
| 118 | // |
---|
| 119 | // ...this is where we are going |
---|
| 120 | // |
---|
| 121 | phi += deltaPhi; |
---|
| 122 | a2 = G4ThreeVector( r[0]*std::cos(phi), r[0]*std::sin(phi), z[0] ); |
---|
| 123 | b2 = G4ThreeVector( r[1]*std::cos(phi), r[1]*std::sin(phi), z[1] ); |
---|
| 124 | c2 = G4ThreeVector( prevRZ->r*std::cos(phi), prevRZ->r*std::sin(phi), prevRZ->z ); |
---|
| 125 | d2 = G4ThreeVector( nextRZ->r*std::cos(phi), nextRZ->r*std::sin(phi), nextRZ->z ); |
---|
| 126 | |
---|
| 127 | G4ThreeVector tt; |
---|
| 128 | |
---|
| 129 | // |
---|
| 130 | // ...build some relevant vectors. |
---|
| 131 | // the point is to sacrifice a little memory with precalcs |
---|
| 132 | // to gain speed |
---|
| 133 | // |
---|
| 134 | vec->center = 0.25*( a1 + a2 + b1 + b2 ); |
---|
| 135 | |
---|
| 136 | tt = b2 + b1 - a2 - a1; |
---|
| 137 | vec->surfRZ = tt.unit(); |
---|
| 138 | if (vec==vecs) lenRZ = 0.25*tt.mag(); |
---|
| 139 | |
---|
| 140 | tt = b2 - b1 + a2 - a1; |
---|
| 141 | vec->surfPhi = tt.unit(); |
---|
| 142 | if (vec==vecs) |
---|
| 143 | { |
---|
| 144 | lenPhi[0] = 0.25*tt.mag(); |
---|
| 145 | tt = b2 - b1; |
---|
| 146 | lenPhi[1] = (0.5*tt.mag()-lenPhi[0])/lenRZ; |
---|
| 147 | } |
---|
| 148 | |
---|
| 149 | tt = vec->surfPhi.cross(vec->surfRZ); |
---|
| 150 | vec->normal = tt.unit(); |
---|
| 151 | |
---|
| 152 | // |
---|
| 153 | // ...edge normals are the average of the normals of |
---|
| 154 | // the two faces they connect. |
---|
| 155 | // |
---|
| 156 | // ...edge normals are necessary if we are to accurately |
---|
| 157 | // decide if a point is "inside" a face. For non-convex |
---|
| 158 | // shapes, it is absolutely necessary to know information |
---|
| 159 | // on adjacent faces to accurate determine this. |
---|
| 160 | // |
---|
| 161 | // ...we don't need them for the phi edges, since that |
---|
| 162 | // information is taken care of internally. The r/z edges, |
---|
| 163 | // however, depend on the adjacent G4PolyhedraSide. |
---|
| 164 | // |
---|
| 165 | G4ThreeVector a12, adj; |
---|
| 166 | |
---|
| 167 | a12 = a2-a1; |
---|
| 168 | |
---|
| 169 | adj = 0.5*(c1+c2-a1-a2); |
---|
| 170 | adj = adj.cross(a12); |
---|
| 171 | adj = adj.unit() + vec->normal; |
---|
| 172 | vec->edgeNorm[0] = adj.unit(); |
---|
| 173 | |
---|
| 174 | a12 = b1-b2; |
---|
| 175 | adj = 0.5*(d1+d2-b1-b2); |
---|
| 176 | adj = adj.cross(a12); |
---|
| 177 | adj = adj.unit() + vec->normal; |
---|
| 178 | vec->edgeNorm[1] = adj.unit(); |
---|
| 179 | |
---|
| 180 | // |
---|
| 181 | // ...the corners are crucial. It is important that |
---|
| 182 | // they are calculated consistently for adjacent |
---|
| 183 | // G4PolyhedraSides, to avoid gaps caused by roundoff. |
---|
| 184 | // |
---|
| 185 | vec->edges[0] = edge; |
---|
| 186 | edge->corner[0] = a1; |
---|
| 187 | edge->corner[1] = b1; |
---|
| 188 | edge++; |
---|
| 189 | vec->edges[1] = edge; |
---|
| 190 | |
---|
| 191 | a1 = a2; |
---|
| 192 | b1 = b2; |
---|
| 193 | c1 = c2; |
---|
| 194 | d1 = d2; |
---|
| 195 | } while( ++vec < vecs+numSide ); |
---|
| 196 | |
---|
| 197 | // |
---|
| 198 | // Clean up hanging edge |
---|
| 199 | // |
---|
| 200 | if (phiIsOpen) |
---|
| 201 | { |
---|
| 202 | edge->corner[0] = a2; |
---|
| 203 | edge->corner[1] = b2; |
---|
| 204 | } |
---|
| 205 | else |
---|
| 206 | { |
---|
| 207 | vecs[numSide-1].edges[1] = edges; |
---|
| 208 | } |
---|
| 209 | |
---|
| 210 | // |
---|
| 211 | // Go back and fill in remaining fields in edges |
---|
| 212 | // |
---|
| 213 | vec = vecs; |
---|
| 214 | G4PolyhedraSideVec *prev = vecs+numSide-1; |
---|
| 215 | do |
---|
| 216 | { |
---|
| 217 | edge = vec->edges[0]; // The edge between prev and vec |
---|
| 218 | |
---|
| 219 | // |
---|
| 220 | // Okay: edge normal is average of normals of adjacent faces |
---|
| 221 | // |
---|
| 222 | G4ThreeVector eNorm = vec->normal + prev->normal; |
---|
| 223 | edge->normal = eNorm.unit(); |
---|
| 224 | |
---|
| 225 | // |
---|
| 226 | // Vertex normal is average of norms of adjacent surfaces (all four) |
---|
| 227 | // However, vec->edgeNorm is unit vector in some direction |
---|
| 228 | // as the sum of normals of adjacent PolyhedraSide with vec. |
---|
| 229 | // The normalization used for this vector should be the same |
---|
| 230 | // for vec and prev. |
---|
| 231 | // |
---|
| 232 | eNorm = vec->edgeNorm[0] + prev->edgeNorm[0]; |
---|
| 233 | edge->cornNorm[0] = eNorm.unit(); |
---|
| 234 | |
---|
| 235 | eNorm = vec->edgeNorm[1] + prev->edgeNorm[1]; |
---|
| 236 | edge->cornNorm[1] = eNorm.unit(); |
---|
| 237 | } while( prev=vec, ++vec < vecs + numSide ); |
---|
| 238 | |
---|
| 239 | if (phiIsOpen) |
---|
| 240 | { |
---|
| 241 | // G4double rFact = std::cos(0.5*deltaPhi); |
---|
| 242 | // |
---|
| 243 | // If phi is open, we need to patch up normals of the |
---|
| 244 | // first and last edges and their corresponding |
---|
| 245 | // vertices. |
---|
| 246 | // |
---|
| 247 | // We use vectors that are in the plane of the |
---|
| 248 | // face. This should be safe. |
---|
| 249 | // |
---|
| 250 | vec = vecs; |
---|
| 251 | |
---|
| 252 | G4ThreeVector normvec = vec->edges[0]->corner[0] |
---|
| 253 | - vec->edges[0]->corner[1]; |
---|
| 254 | normvec = normvec.cross(vec->normal); |
---|
| 255 | if (normvec.dot(vec->surfPhi) > 0) normvec = -normvec; |
---|
| 256 | |
---|
| 257 | vec->edges[0]->normal = normvec.unit(); |
---|
| 258 | |
---|
| 259 | vec->edges[0]->cornNorm[0] = (vec->edges[0]->corner[0] |
---|
| 260 | - vec->center).unit(); |
---|
| 261 | vec->edges[0]->cornNorm[1] = (vec->edges[0]->corner[1] |
---|
| 262 | - vec->center).unit(); |
---|
| 263 | |
---|
| 264 | // |
---|
| 265 | // Repeat for ending phi |
---|
| 266 | // |
---|
| 267 | vec = vecs + numSide - 1; |
---|
| 268 | |
---|
| 269 | normvec = vec->edges[1]->corner[0] - vec->edges[1]->corner[1]; |
---|
| 270 | normvec = normvec.cross(vec->normal); |
---|
| 271 | if (normvec.dot(vec->surfPhi) < 0) normvec = -normvec; |
---|
| 272 | |
---|
| 273 | vec->edges[1]->normal = normvec.unit(); |
---|
| 274 | |
---|
| 275 | vec->edges[1]->cornNorm[0] = (vec->edges[1]->corner[0] |
---|
| 276 | - vec->center).unit(); |
---|
| 277 | vec->edges[1]->cornNorm[1] = (vec->edges[1]->corner[1] |
---|
| 278 | - vec->center).unit(); |
---|
| 279 | } |
---|
| 280 | |
---|
| 281 | // |
---|
| 282 | // edgeNorm is the factor one multiplies the distance along vector phi |
---|
| 283 | // on the surface of one of our sides in order to calculate the distance |
---|
| 284 | // from the edge. (see routine DistanceAway) |
---|
| 285 | // |
---|
| 286 | edgeNorm = 1.0/std::sqrt( 1.0 + lenPhi[1]*lenPhi[1] ); |
---|
| 287 | } |
---|
| 288 | |
---|
| 289 | |
---|
| 290 | // |
---|
| 291 | // Fake default constructor - sets only member data and allocates memory |
---|
| 292 | // for usage restricted to object persistency. |
---|
| 293 | // |
---|
| 294 | G4PolyhedraSide::G4PolyhedraSide( __void__&) |
---|
| 295 | : cone(0), vecs(0), edges(0) |
---|
| 296 | { |
---|
| 297 | } |
---|
| 298 | |
---|
| 299 | |
---|
| 300 | // |
---|
| 301 | // Destructor |
---|
| 302 | // |
---|
| 303 | G4PolyhedraSide::~G4PolyhedraSide() |
---|
| 304 | { |
---|
| 305 | delete cone; |
---|
| 306 | delete [] vecs; |
---|
| 307 | delete [] edges; |
---|
| 308 | } |
---|
| 309 | |
---|
| 310 | |
---|
| 311 | // |
---|
| 312 | // Copy constructor |
---|
| 313 | // |
---|
| 314 | G4PolyhedraSide::G4PolyhedraSide( const G4PolyhedraSide &source ) |
---|
| 315 | : G4VCSGface() |
---|
| 316 | { |
---|
| 317 | CopyStuff( source ); |
---|
| 318 | } |
---|
| 319 | |
---|
| 320 | |
---|
| 321 | // |
---|
| 322 | // Assignment operator |
---|
| 323 | // |
---|
| 324 | G4PolyhedraSide& G4PolyhedraSide::operator=( const G4PolyhedraSide &source ) |
---|
| 325 | { |
---|
| 326 | if (this == &source) return *this; |
---|
| 327 | |
---|
| 328 | delete cone; |
---|
| 329 | delete [] vecs; |
---|
| 330 | delete [] edges; |
---|
| 331 | |
---|
| 332 | CopyStuff( source ); |
---|
| 333 | |
---|
| 334 | return *this; |
---|
| 335 | } |
---|
| 336 | |
---|
| 337 | |
---|
| 338 | // |
---|
| 339 | // CopyStuff |
---|
| 340 | // |
---|
| 341 | void G4PolyhedraSide::CopyStuff( const G4PolyhedraSide &source ) |
---|
| 342 | { |
---|
| 343 | // |
---|
| 344 | // The simple stuff |
---|
| 345 | // |
---|
| 346 | numSide = source.numSide; |
---|
| 347 | r[0] = source.r[0]; |
---|
| 348 | r[1] = source.r[1]; |
---|
| 349 | z[0] = source.z[0]; |
---|
| 350 | z[1] = source.z[1]; |
---|
| 351 | startPhi = source.startPhi; |
---|
| 352 | deltaPhi = source.deltaPhi; |
---|
| 353 | endPhi = source.endPhi; |
---|
| 354 | phiIsOpen = source.phiIsOpen; |
---|
| 355 | allBehind = source.allBehind; |
---|
| 356 | |
---|
| 357 | lenRZ = source.lenRZ; |
---|
| 358 | lenPhi[0] = source.lenPhi[0]; |
---|
| 359 | lenPhi[1] = source.lenPhi[1]; |
---|
| 360 | edgeNorm = source.edgeNorm; |
---|
| 361 | |
---|
| 362 | kCarTolerance = source.kCarTolerance; |
---|
[850] | 363 | fSurfaceArea = source.fSurfaceArea; |
---|
| 364 | |
---|
[831] | 365 | cone = new G4IntersectingCone( *source.cone ); |
---|
| 366 | |
---|
| 367 | // |
---|
| 368 | // Duplicate edges |
---|
| 369 | // |
---|
| 370 | G4int numEdges = phiIsOpen ? numSide+1 : numSide; |
---|
| 371 | edges = new G4PolyhedraSideEdge[numEdges]; |
---|
| 372 | |
---|
| 373 | G4PolyhedraSideEdge *edge = edges, |
---|
| 374 | *sourceEdge = source.edges; |
---|
| 375 | do |
---|
| 376 | { |
---|
| 377 | *edge = *sourceEdge; |
---|
| 378 | } while( ++sourceEdge, ++edge < edges + numEdges); |
---|
| 379 | |
---|
| 380 | // |
---|
| 381 | // Duplicate vecs |
---|
| 382 | // |
---|
| 383 | vecs = new G4PolyhedraSideVec[numSide]; |
---|
| 384 | |
---|
| 385 | G4PolyhedraSideVec *vec = vecs, |
---|
| 386 | *sourceVec = source.vecs; |
---|
| 387 | do |
---|
| 388 | { |
---|
| 389 | *vec = *sourceVec; |
---|
| 390 | vec->edges[0] = edges + (sourceVec->edges[0] - source.edges); |
---|
| 391 | vec->edges[1] = edges + (sourceVec->edges[1] - source.edges); |
---|
| 392 | } while( ++sourceVec, ++vec < vecs + numSide ); |
---|
| 393 | } |
---|
| 394 | |
---|
| 395 | |
---|
| 396 | // |
---|
| 397 | // Intersect |
---|
| 398 | // |
---|
| 399 | // Decide if a line intersects the face. |
---|
| 400 | // |
---|
| 401 | // Arguments: |
---|
| 402 | // p = (in) starting point of line segment |
---|
| 403 | // v = (in) direction of line segment (assumed a unit vector) |
---|
| 404 | // A, B = (in) 2d transform variables (see note top of file) |
---|
| 405 | // normSign = (in) desired sign for dot product with normal (see below) |
---|
| 406 | // surfTolerance = (in) minimum distance from the surface |
---|
| 407 | // vecs = (in) Vector set array |
---|
| 408 | // distance = (out) distance to surface furfilling all requirements |
---|
| 409 | // distFromSurface = (out) distance from the surface |
---|
| 410 | // thisNormal = (out) normal vector of the intersecting surface |
---|
| 411 | // |
---|
| 412 | // Return value: |
---|
| 413 | // true if an intersection is found. Otherwise, output parameters are |
---|
| 414 | // undefined. |
---|
| 415 | // |
---|
| 416 | // Notes: |
---|
| 417 | // * normSign: if we are "inside" the shape and only want to find out how far |
---|
| 418 | // to leave the shape, we only want to consider intersections with surfaces in |
---|
| 419 | // which the trajectory is leaving the shape. Since the normal vectors to the |
---|
| 420 | // surface always point outwards from the inside, this means we want the dot |
---|
| 421 | // product of the trajectory direction v and the normal of the side normals[i] |
---|
| 422 | // to be positive. Thus, we should specify normSign as +1.0. Otherwise, if |
---|
| 423 | // we are outside and want to go in, normSign should be set to -1.0. |
---|
| 424 | // Don't set normSign to zero, or you will get no intersections! |
---|
| 425 | // |
---|
| 426 | // * surfTolerance: see notes on argument "surfTolerance" in routine |
---|
| 427 | // "IntersectSidePlane". |
---|
| 428 | // ----HOWEVER---- We should *not* apply this surface tolerance if the |
---|
| 429 | // starting point is not within phi or z of the surface. Specifically, |
---|
| 430 | // if the starting point p angle in x/y places it on a separate side from the |
---|
| 431 | // intersection or if the starting point p is outside the z bounds of the |
---|
| 432 | // segment, surfTolerance must be ignored or we should *always* accept the |
---|
| 433 | // intersection! |
---|
| 434 | // This is simply because the sides do not have infinite extent. |
---|
| 435 | // |
---|
| 436 | // |
---|
| 437 | G4bool G4PolyhedraSide::Intersect( const G4ThreeVector &p, |
---|
| 438 | const G4ThreeVector &v, |
---|
| 439 | G4bool outgoing, |
---|
| 440 | G4double surfTolerance, |
---|
| 441 | G4double &distance, |
---|
| 442 | G4double &distFromSurface, |
---|
| 443 | G4ThreeVector &normal, |
---|
| 444 | G4bool &isAllBehind ) |
---|
| 445 | { |
---|
| 446 | G4double normSign = outgoing ? +1 : -1; |
---|
| 447 | |
---|
| 448 | // |
---|
| 449 | // ------------------TO BE IMPLEMENTED--------------------- |
---|
| 450 | // Testing the intersection of individual phi faces is |
---|
| 451 | // pretty straight forward. The simple thing therefore is to |
---|
| 452 | // form a loop and check them all in sequence. |
---|
| 453 | // |
---|
| 454 | // But, I worry about one day someone making |
---|
| 455 | // a polygon with a thousands sides. A linear search |
---|
| 456 | // would not be ideal in such a case. |
---|
| 457 | // |
---|
| 458 | // So, it would be nice to be able to quickly decide |
---|
| 459 | // which face would be intersected. One can make a very |
---|
| 460 | // good guess by using the intersection with a cone. |
---|
| 461 | // However, this is only reliable in 99% of the cases. |
---|
| 462 | // |
---|
| 463 | // My solution: make a decent guess as to the one or |
---|
| 464 | // two potential faces might get intersected, and then |
---|
| 465 | // test them. If we have the wrong face, use the test |
---|
| 466 | // to make a better guess. |
---|
| 467 | // |
---|
| 468 | // Since we might have two guesses, form a queue of |
---|
| 469 | // potential intersecting faces. Keep an array of |
---|
| 470 | // already tested faces to avoid doing one more than |
---|
| 471 | // once. |
---|
| 472 | // |
---|
| 473 | // Result: at worst, an iterative search. On average, |
---|
| 474 | // a little more than two tests would be required. |
---|
| 475 | // |
---|
| 476 | G4ThreeVector q = p + v; |
---|
| 477 | |
---|
| 478 | G4int face = 0; |
---|
| 479 | G4PolyhedraSideVec *vec = vecs; |
---|
| 480 | do |
---|
| 481 | { |
---|
| 482 | // |
---|
| 483 | // Correct normal? |
---|
| 484 | // |
---|
| 485 | G4double dotProd = normSign*v.dot(vec->normal); |
---|
| 486 | if (dotProd <= 0) continue; |
---|
| 487 | |
---|
| 488 | // |
---|
| 489 | // Is this face in front of the point along the trajectory? |
---|
| 490 | // |
---|
| 491 | G4ThreeVector delta = p - vec->center; |
---|
| 492 | distFromSurface = -normSign*delta.dot(vec->normal); |
---|
| 493 | |
---|
| 494 | if (distFromSurface < -surfTolerance) continue; |
---|
| 495 | |
---|
| 496 | // |
---|
| 497 | // phi |
---|
| 498 | // c -------- d ^ |
---|
| 499 | // | | | |
---|
| 500 | // a -------- b +---> r/z |
---|
| 501 | // |
---|
| 502 | // |
---|
| 503 | // Do we remain on this particular segment? |
---|
| 504 | // |
---|
| 505 | G4ThreeVector qc = q - vec->edges[1]->corner[0]; |
---|
| 506 | G4ThreeVector qd = q - vec->edges[1]->corner[1]; |
---|
| 507 | |
---|
| 508 | if (normSign*qc.cross(qd).dot(v) < 0) continue; |
---|
| 509 | |
---|
| 510 | G4ThreeVector qa = q - vec->edges[0]->corner[0]; |
---|
| 511 | G4ThreeVector qb = q - vec->edges[0]->corner[1]; |
---|
| 512 | |
---|
| 513 | if (normSign*qa.cross(qb).dot(v) > 0) continue; |
---|
| 514 | |
---|
| 515 | // |
---|
| 516 | // We found the one and only segment we might be intersecting. |
---|
| 517 | // Do we remain within r/z bounds? |
---|
| 518 | // |
---|
| 519 | |
---|
| 520 | if (r[0] > 1/kInfinity && normSign*qa.cross(qc).dot(v) < 0) return false; |
---|
| 521 | if (r[1] > 1/kInfinity && normSign*qb.cross(qd).dot(v) > 0) return false; |
---|
| 522 | |
---|
| 523 | // |
---|
| 524 | // We allow the face to be slightly behind the trajectory |
---|
| 525 | // (surface tolerance) only if the point p is within |
---|
| 526 | // the vicinity of the face |
---|
| 527 | // |
---|
| 528 | if (distFromSurface < 0) |
---|
| 529 | { |
---|
| 530 | G4ThreeVector ps = p - vec->center; |
---|
| 531 | |
---|
| 532 | G4double rz = ps.dot(vec->surfRZ); |
---|
| 533 | if (std::fabs(rz) > lenRZ+surfTolerance) return false; |
---|
| 534 | |
---|
| 535 | G4double pp = ps.dot(vec->surfPhi); |
---|
| 536 | if (std::fabs(pp) > lenPhi[0] + lenPhi[1]*rz + surfTolerance) return false; |
---|
| 537 | } |
---|
| 538 | |
---|
| 539 | |
---|
| 540 | // |
---|
| 541 | // Intersection found. Return answer. |
---|
| 542 | // |
---|
| 543 | distance = distFromSurface/dotProd; |
---|
| 544 | normal = vec->normal; |
---|
| 545 | isAllBehind = allBehind; |
---|
| 546 | return true; |
---|
| 547 | } while( ++vec, ++face < numSide ); |
---|
| 548 | |
---|
| 549 | // |
---|
| 550 | // Oh well. Better luck next time. |
---|
| 551 | // |
---|
| 552 | return false; |
---|
| 553 | } |
---|
| 554 | |
---|
| 555 | |
---|
| 556 | G4double G4PolyhedraSide::Distance( const G4ThreeVector &p, G4bool outgoing ) |
---|
| 557 | { |
---|
| 558 | G4double normSign = outgoing ? -1 : +1; |
---|
| 559 | |
---|
| 560 | // |
---|
| 561 | // Try the closest phi segment first |
---|
| 562 | // |
---|
| 563 | G4int iPhi = ClosestPhiSegment( p.phi() ); |
---|
| 564 | |
---|
| 565 | G4ThreeVector pdotc = p - vecs[iPhi].center; |
---|
| 566 | G4double normDist = pdotc.dot(vecs[iPhi].normal); |
---|
| 567 | |
---|
| 568 | if (normSign*normDist > -0.5*kCarTolerance) |
---|
| 569 | { |
---|
| 570 | return DistanceAway( p, vecs[iPhi], &normDist ); |
---|
| 571 | } |
---|
| 572 | |
---|
| 573 | // |
---|
| 574 | // Now we have an interesting problem... do we try to find the |
---|
| 575 | // closest facing side?? |
---|
| 576 | // |
---|
| 577 | // Considered carefully, the answer is no. We know that if we |
---|
| 578 | // are asking for the distance out, we are supposed to be inside, |
---|
| 579 | // and vice versa. |
---|
| 580 | // |
---|
| 581 | |
---|
| 582 | return kInfinity; |
---|
| 583 | } |
---|
| 584 | |
---|
| 585 | |
---|
| 586 | // |
---|
| 587 | // Inside |
---|
| 588 | // |
---|
| 589 | EInside G4PolyhedraSide::Inside( const G4ThreeVector &p, |
---|
| 590 | G4double tolerance, |
---|
| 591 | G4double *bestDistance ) |
---|
| 592 | { |
---|
| 593 | // |
---|
| 594 | // Which phi segment is closest to this point? |
---|
| 595 | // |
---|
| 596 | G4int iPhi = ClosestPhiSegment( p.phi() ); |
---|
| 597 | |
---|
| 598 | G4double norm; |
---|
| 599 | |
---|
| 600 | // |
---|
| 601 | // Get distance to this segment |
---|
| 602 | // |
---|
| 603 | *bestDistance = DistanceToOneSide( p, vecs[iPhi], &norm ); |
---|
| 604 | |
---|
| 605 | // |
---|
| 606 | // Use distance along normal to decide return value |
---|
| 607 | // |
---|
| 608 | if ( (std::fabs(norm) < tolerance) && (*bestDistance < 2.0*tolerance) ) |
---|
| 609 | return kSurface; |
---|
| 610 | else if (norm < 0) |
---|
| 611 | return kInside; |
---|
| 612 | else |
---|
| 613 | return kOutside; |
---|
| 614 | } |
---|
| 615 | |
---|
| 616 | |
---|
| 617 | // |
---|
| 618 | // Normal |
---|
| 619 | // |
---|
| 620 | G4ThreeVector G4PolyhedraSide::Normal( const G4ThreeVector &p, |
---|
| 621 | G4double *bestDistance ) |
---|
| 622 | { |
---|
| 623 | // |
---|
| 624 | // Which phi segment is closest to this point? |
---|
| 625 | // |
---|
| 626 | G4int iPhi = ClosestPhiSegment( p.phi() ); |
---|
| 627 | |
---|
| 628 | // |
---|
| 629 | // Get distance to this segment |
---|
| 630 | // |
---|
| 631 | G4double norm; |
---|
| 632 | *bestDistance = DistanceToOneSide( p, vecs[iPhi], &norm ); |
---|
| 633 | |
---|
| 634 | return vecs[iPhi].normal; |
---|
| 635 | } |
---|
| 636 | |
---|
| 637 | |
---|
| 638 | // |
---|
| 639 | // Extent |
---|
| 640 | // |
---|
| 641 | G4double G4PolyhedraSide::Extent( const G4ThreeVector axis ) |
---|
| 642 | { |
---|
| 643 | if (axis.perp2() < DBL_MIN) |
---|
| 644 | { |
---|
| 645 | // |
---|
| 646 | // Special case |
---|
| 647 | // |
---|
| 648 | return axis.z() < 0 ? -cone->ZLo() : cone->ZHi(); |
---|
| 649 | } |
---|
| 650 | |
---|
| 651 | G4int iPhi, i1, i2; |
---|
| 652 | G4double best; |
---|
| 653 | G4ThreeVector *list[4]; |
---|
| 654 | |
---|
| 655 | // |
---|
| 656 | // Which phi segment, if any, does the axis belong to |
---|
| 657 | // |
---|
| 658 | iPhi = PhiSegment( axis.phi() ); |
---|
| 659 | |
---|
| 660 | if (iPhi < 0) |
---|
| 661 | { |
---|
| 662 | // |
---|
| 663 | // No phi segment? Check front edge of first side and |
---|
| 664 | // last edge of second side |
---|
| 665 | // |
---|
| 666 | i1 = 0; i2 = numSide-1; |
---|
| 667 | } |
---|
| 668 | else |
---|
| 669 | { |
---|
| 670 | // |
---|
| 671 | // Check all corners of matching phi side |
---|
| 672 | // |
---|
| 673 | i1 = iPhi; i2 = iPhi; |
---|
| 674 | } |
---|
| 675 | |
---|
| 676 | list[0] = vecs[i1].edges[0]->corner; |
---|
| 677 | list[1] = vecs[i1].edges[0]->corner+1; |
---|
| 678 | list[2] = vecs[i2].edges[1]->corner; |
---|
| 679 | list[3] = vecs[i2].edges[1]->corner+1; |
---|
| 680 | |
---|
| 681 | // |
---|
| 682 | // Who's biggest? |
---|
| 683 | // |
---|
| 684 | best = -kInfinity; |
---|
| 685 | G4ThreeVector **vec = list; |
---|
| 686 | do |
---|
| 687 | { |
---|
| 688 | G4double answer = (*vec)->dot(axis); |
---|
| 689 | if (answer > best) best = answer; |
---|
| 690 | } while( ++vec < list+4 ); |
---|
| 691 | |
---|
| 692 | return best; |
---|
| 693 | } |
---|
| 694 | |
---|
| 695 | |
---|
| 696 | // |
---|
| 697 | // CalculateExtent |
---|
| 698 | // |
---|
| 699 | // See notes in G4VCSGface |
---|
| 700 | // |
---|
| 701 | void G4PolyhedraSide::CalculateExtent( const EAxis axis, |
---|
| 702 | const G4VoxelLimits &voxelLimit, |
---|
| 703 | const G4AffineTransform &transform, |
---|
| 704 | G4SolidExtentList &extentList ) |
---|
| 705 | { |
---|
| 706 | // |
---|
| 707 | // Loop over all sides |
---|
| 708 | // |
---|
| 709 | G4PolyhedraSideVec *vec = vecs; |
---|
| 710 | do |
---|
| 711 | { |
---|
| 712 | // |
---|
| 713 | // Fill our polygon with the four corners of |
---|
| 714 | // this side, after the specified transformation |
---|
| 715 | // |
---|
| 716 | G4ClippablePolygon polygon; |
---|
| 717 | |
---|
| 718 | polygon.AddVertexInOrder(transform. |
---|
| 719 | TransformPoint(vec->edges[0]->corner[0])); |
---|
| 720 | polygon.AddVertexInOrder(transform. |
---|
| 721 | TransformPoint(vec->edges[0]->corner[1])); |
---|
| 722 | polygon.AddVertexInOrder(transform. |
---|
| 723 | TransformPoint(vec->edges[1]->corner[1])); |
---|
| 724 | polygon.AddVertexInOrder(transform. |
---|
| 725 | TransformPoint(vec->edges[1]->corner[0])); |
---|
| 726 | |
---|
| 727 | // |
---|
| 728 | // Get extent |
---|
| 729 | // |
---|
| 730 | if (polygon.PartialClip( voxelLimit, axis )) |
---|
| 731 | { |
---|
| 732 | // |
---|
| 733 | // Get dot product of normal along target axis |
---|
| 734 | // |
---|
| 735 | polygon.SetNormal( transform.TransformAxis(vec->normal) ); |
---|
| 736 | |
---|
| 737 | extentList.AddSurface( polygon ); |
---|
| 738 | } |
---|
| 739 | } while( ++vec < vecs+numSide ); |
---|
| 740 | |
---|
| 741 | return; |
---|
| 742 | } |
---|
| 743 | |
---|
| 744 | |
---|
| 745 | // |
---|
| 746 | // IntersectSidePlane |
---|
| 747 | // |
---|
| 748 | // Decide if a line correctly intersects one side plane of our segment. |
---|
| 749 | // It is assumed that the correct side has been chosen, and thus only |
---|
| 750 | // the z bounds (of the entire segment) are checked. |
---|
| 751 | // |
---|
| 752 | // normSign - To be multiplied against normal: |
---|
| 753 | // = +1.0 normal is unchanged |
---|
| 754 | // = -1.0 normal is reversed (now points inward) |
---|
| 755 | // |
---|
| 756 | // Arguments: |
---|
| 757 | // p - (in) Point |
---|
| 758 | // v - (in) Direction |
---|
| 759 | // vec - (in) Description record of the side plane |
---|
| 760 | // normSign - (in) Sign (+/- 1) to apply to normal |
---|
| 761 | // surfTolerance - (in) Surface tolerance (generally > 0, see below) |
---|
| 762 | // distance - (out) Distance along v to intersection |
---|
| 763 | // distFromSurface - (out) Distance from surface normal |
---|
| 764 | // |
---|
| 765 | // Notes: |
---|
| 766 | // surfTolerance - Used to decide if a point is behind the surface, |
---|
| 767 | // a point is allow to be -surfTolerance behind the |
---|
| 768 | // surface (as measured along the normal), but *only* |
---|
| 769 | // if the point is within the r/z bounds + surfTolerance |
---|
| 770 | // of the segment. |
---|
| 771 | // |
---|
| 772 | G4bool G4PolyhedraSide::IntersectSidePlane( const G4ThreeVector &p, |
---|
| 773 | const G4ThreeVector &v, |
---|
| 774 | const G4PolyhedraSideVec vec, |
---|
| 775 | G4double normSign, |
---|
| 776 | G4double surfTolerance, |
---|
| 777 | G4double &distance, |
---|
| 778 | G4double &distFromSurface ) |
---|
| 779 | { |
---|
| 780 | // |
---|
| 781 | // Correct normal? Here we have straight sides, and can safely ignore |
---|
| 782 | // intersections where the dot product with the normal is zero. |
---|
| 783 | // |
---|
| 784 | G4double dotProd = normSign*v.dot(vec.normal); |
---|
| 785 | |
---|
| 786 | if (dotProd <= 0) return false; |
---|
| 787 | |
---|
| 788 | // |
---|
| 789 | // Calculate distance to surface. If the side is too far |
---|
| 790 | // behind the point, we must reject it. |
---|
| 791 | // |
---|
| 792 | G4ThreeVector delta = p - vec.center; |
---|
| 793 | distFromSurface = -normSign*delta.dot(vec.normal); |
---|
| 794 | |
---|
| 795 | if (distFromSurface < -surfTolerance) return false; |
---|
| 796 | |
---|
| 797 | // |
---|
| 798 | // Calculate precise distance to intersection with the side |
---|
| 799 | // (along the trajectory, not normal to the surface) |
---|
| 800 | // |
---|
| 801 | distance = distFromSurface/dotProd; |
---|
| 802 | |
---|
| 803 | // |
---|
| 804 | // Do we fall off the r/z extent of the segment? |
---|
| 805 | // |
---|
| 806 | // Calculate this very, very carefully! Why? |
---|
| 807 | // 1. If a RZ end is at R=0, you can't miss! |
---|
| 808 | // 2. If you just fall off in RZ, the answer must |
---|
| 809 | // be consistent with adjacent G4PolyhedraSide faces. |
---|
| 810 | // (2) implies that only variables used by other G4PolyhedraSide |
---|
| 811 | // faces may be used, which includes only: p, v, and the edge corners. |
---|
| 812 | // It also means that one side is a ">" or "<", which the other |
---|
| 813 | // must be ">=" or "<=". Fortunately, this isn't a new problem. |
---|
| 814 | // The solution below I borrowed from Joseph O'Rourke, |
---|
| 815 | // "Computational Geometry in C (Second Edition)" |
---|
| 816 | // See: http://cs.smith.edu/~orourke/ |
---|
| 817 | // |
---|
| 818 | G4ThreeVector ic = p + distance*v - vec.center; |
---|
| 819 | G4double atRZ = vec.surfRZ.dot(ic); |
---|
| 820 | |
---|
| 821 | if (atRZ < 0) |
---|
| 822 | { |
---|
| 823 | if (r[0]==0) return true; // Can't miss! |
---|
| 824 | |
---|
| 825 | if (atRZ < -lenRZ*1.2) return false; // Forget it! Missed by a mile. |
---|
| 826 | |
---|
| 827 | G4ThreeVector q = p + v; |
---|
| 828 | G4ThreeVector qa = q - vec.edges[0]->corner[0], |
---|
| 829 | qb = q - vec.edges[1]->corner[0]; |
---|
| 830 | G4ThreeVector qacb = qa.cross(qb); |
---|
| 831 | if (normSign*qacb.dot(v) < 0) return false; |
---|
| 832 | |
---|
| 833 | if (distFromSurface < 0) |
---|
| 834 | { |
---|
| 835 | if (atRZ < -lenRZ-surfTolerance) return false; |
---|
| 836 | } |
---|
| 837 | } |
---|
| 838 | else if (atRZ > 0) |
---|
| 839 | { |
---|
| 840 | if (r[1]==0) return true; // Can't miss! |
---|
| 841 | |
---|
| 842 | if (atRZ > lenRZ*1.2) return false; // Missed by a mile |
---|
| 843 | |
---|
| 844 | G4ThreeVector q = p + v; |
---|
| 845 | G4ThreeVector qa = q - vec.edges[0]->corner[1], |
---|
| 846 | qb = q - vec.edges[1]->corner[1]; |
---|
| 847 | G4ThreeVector qacb = qa.cross(qb); |
---|
| 848 | if (normSign*qacb.dot(v) >= 0) return false; |
---|
| 849 | |
---|
| 850 | if (distFromSurface < 0) |
---|
| 851 | { |
---|
| 852 | if (atRZ > lenRZ+surfTolerance) return false; |
---|
| 853 | } |
---|
| 854 | } |
---|
| 855 | |
---|
| 856 | return true; |
---|
| 857 | } |
---|
| 858 | |
---|
| 859 | |
---|
| 860 | // |
---|
| 861 | // LineHitsSegments |
---|
| 862 | // |
---|
| 863 | // Calculate which phi segments a line intersects in three dimensions. |
---|
| 864 | // No check is made as to whether the intersections are within the z bounds of |
---|
| 865 | // the segment. |
---|
| 866 | // |
---|
| 867 | G4int G4PolyhedraSide::LineHitsSegments( const G4ThreeVector &p, |
---|
| 868 | const G4ThreeVector &v, |
---|
| 869 | G4int *i1, G4int *i2 ) |
---|
| 870 | { |
---|
| 871 | G4double s1, s2; |
---|
| 872 | // |
---|
| 873 | // First, decide if and where the line intersects the cone |
---|
| 874 | // |
---|
| 875 | G4int n = cone->LineHitsCone( p, v, &s1, &s2 ); |
---|
| 876 | |
---|
| 877 | if (n==0) return 0; |
---|
| 878 | |
---|
| 879 | // |
---|
| 880 | // Try first intersection. |
---|
| 881 | // |
---|
| 882 | *i1 = PhiSegment( std::atan2( p.y() + s1*v.y(), p.x() + s1*v.x() ) ); |
---|
| 883 | if (n==1) |
---|
| 884 | { |
---|
| 885 | return (*i1 < 0) ? 0 : 1; |
---|
| 886 | } |
---|
| 887 | |
---|
| 888 | // |
---|
| 889 | // Try second intersection |
---|
| 890 | // |
---|
| 891 | *i2 = PhiSegment( std::atan2( p.y() + s2*v.y(), p.x() + s2*v.x() ) ); |
---|
| 892 | if (*i1 == *i2) return 0; |
---|
| 893 | |
---|
| 894 | if (*i1 < 0) |
---|
| 895 | { |
---|
| 896 | if (*i2 < 0) return 0; |
---|
| 897 | *i1 = *i2; |
---|
| 898 | return 1; |
---|
| 899 | } |
---|
| 900 | |
---|
| 901 | if (*i2 < 0) return 1; |
---|
| 902 | |
---|
| 903 | return 2; |
---|
| 904 | } |
---|
| 905 | |
---|
| 906 | |
---|
| 907 | // |
---|
| 908 | // ClosestPhiSegment |
---|
| 909 | // |
---|
| 910 | // Decide which phi segment is closest in phi to the point. |
---|
| 911 | // The result is the same as PhiSegment if there is no phi opening. |
---|
| 912 | // |
---|
| 913 | G4int G4PolyhedraSide::ClosestPhiSegment( G4double phi0 ) |
---|
| 914 | { |
---|
| 915 | G4int iPhi = PhiSegment( phi0 ); |
---|
| 916 | if (iPhi >= 0) return iPhi; |
---|
| 917 | |
---|
| 918 | // |
---|
| 919 | // Boogers! The points falls inside the phi segment. |
---|
| 920 | // Look for the closest point: the start, or end |
---|
| 921 | // |
---|
| 922 | G4double phi = phi0; |
---|
| 923 | |
---|
| 924 | while( phi < startPhi ) phi += twopi; |
---|
| 925 | G4double d1 = phi-endPhi; |
---|
| 926 | |
---|
| 927 | while( phi > startPhi ) phi -= twopi; |
---|
| 928 | G4double d2 = startPhi-phi; |
---|
| 929 | |
---|
| 930 | return (d2 < d1) ? 0 : numSide-1; |
---|
| 931 | } |
---|
| 932 | |
---|
| 933 | |
---|
| 934 | // |
---|
| 935 | // PhiSegment |
---|
| 936 | // |
---|
| 937 | // Decide which phi segment an angle belongs to, counting from zero. |
---|
| 938 | // A value of -1 indicates that the phi value is outside the shape |
---|
| 939 | // (only possible if phiTotal < 360 degrees). |
---|
| 940 | // |
---|
| 941 | G4int G4PolyhedraSide::PhiSegment( G4double phi0 ) |
---|
| 942 | { |
---|
| 943 | // |
---|
| 944 | // How far are we from phiStart? Come up with a positive answer |
---|
| 945 | // that is less than 2*PI |
---|
| 946 | // |
---|
| 947 | G4double phi = phi0 - startPhi; |
---|
| 948 | while( phi < 0 ) phi += twopi; |
---|
| 949 | while( phi > twopi ) phi -= twopi; |
---|
| 950 | |
---|
| 951 | // |
---|
| 952 | // Divide |
---|
| 953 | // |
---|
| 954 | G4int answer = (G4int)(phi/deltaPhi); |
---|
| 955 | |
---|
| 956 | if (answer >= numSide) |
---|
| 957 | { |
---|
| 958 | if (phiIsOpen) |
---|
| 959 | { |
---|
| 960 | return -1; // Looks like we missed |
---|
| 961 | } |
---|
| 962 | else |
---|
| 963 | { |
---|
| 964 | answer = numSide-1; // Probably just roundoff |
---|
| 965 | } |
---|
| 966 | } |
---|
| 967 | |
---|
| 968 | return answer; |
---|
| 969 | } |
---|
| 970 | |
---|
| 971 | |
---|
| 972 | // |
---|
| 973 | // DistanceToOneSide |
---|
| 974 | // |
---|
| 975 | // Arguments: |
---|
| 976 | // p - (in) Point to check |
---|
| 977 | // vec - (in) vector set of this side |
---|
| 978 | // normDist - (out) distance normal to the side or edge, as appropriate, signed |
---|
| 979 | // Return value = total distance from the side |
---|
| 980 | // |
---|
| 981 | G4double G4PolyhedraSide::DistanceToOneSide( const G4ThreeVector &p, |
---|
| 982 | const G4PolyhedraSideVec &vec, |
---|
| 983 | G4double *normDist ) |
---|
| 984 | { |
---|
| 985 | G4ThreeVector pc = p - vec.center; |
---|
| 986 | |
---|
| 987 | // |
---|
| 988 | // Get normal distance |
---|
| 989 | // |
---|
| 990 | *normDist = vec.normal.dot(pc); |
---|
| 991 | |
---|
| 992 | // |
---|
| 993 | // Add edge penalty |
---|
| 994 | // |
---|
| 995 | return DistanceAway( p, vec, normDist ); |
---|
| 996 | } |
---|
| 997 | |
---|
| 998 | |
---|
| 999 | // |
---|
| 1000 | // DistanceAway |
---|
| 1001 | // |
---|
| 1002 | // Add distance from side edges, if necesssary, to total distance, |
---|
| 1003 | // and updates normDist appropriate depending on edge normals. |
---|
| 1004 | // |
---|
| 1005 | G4double G4PolyhedraSide::DistanceAway( const G4ThreeVector &p, |
---|
| 1006 | const G4PolyhedraSideVec &vec, |
---|
| 1007 | G4double *normDist ) |
---|
| 1008 | { |
---|
| 1009 | G4double distOut2; |
---|
| 1010 | G4ThreeVector pc = p - vec.center; |
---|
| 1011 | G4double distFaceNorm = *normDist; |
---|
| 1012 | |
---|
| 1013 | // |
---|
| 1014 | // Okay, are we inside bounds? |
---|
| 1015 | // |
---|
| 1016 | G4double pcDotRZ = pc.dot(vec.surfRZ); |
---|
| 1017 | G4double pcDotPhi = pc.dot(vec.surfPhi); |
---|
| 1018 | |
---|
| 1019 | // |
---|
| 1020 | // Go through all permutations. |
---|
| 1021 | // Phi |
---|
| 1022 | // | | ^ |
---|
| 1023 | // B | H | E | |
---|
| 1024 | // ------[1]------------[3]----- | |
---|
| 1025 | // |XXXXXXXXXXXXXX| +----> RZ |
---|
| 1026 | // C |XXXXXXXXXXXXXX| F |
---|
| 1027 | // |XXXXXXXXXXXXXX| |
---|
| 1028 | // ------[0]------------[2]---- |
---|
| 1029 | // A | G | D |
---|
| 1030 | // | | |
---|
| 1031 | // |
---|
| 1032 | // It's real messy, but at least it's quick |
---|
| 1033 | // |
---|
| 1034 | |
---|
| 1035 | if (pcDotRZ < -lenRZ) |
---|
| 1036 | { |
---|
| 1037 | G4double lenPhiZ = lenPhi[0] - lenRZ*lenPhi[1]; |
---|
| 1038 | G4double distOutZ = pcDotRZ+lenRZ; |
---|
| 1039 | // |
---|
| 1040 | // Below in RZ |
---|
| 1041 | // |
---|
| 1042 | if (pcDotPhi < -lenPhiZ) |
---|
| 1043 | { |
---|
| 1044 | // |
---|
| 1045 | // ...and below in phi. Find distance to point (A) |
---|
| 1046 | // |
---|
| 1047 | G4double distOutPhi = pcDotPhi+lenPhiZ; |
---|
| 1048 | distOut2 = distOutPhi*distOutPhi + distOutZ*distOutZ; |
---|
| 1049 | G4ThreeVector pa = p - vec.edges[0]->corner[0]; |
---|
| 1050 | *normDist = pa.dot(vec.edges[0]->cornNorm[0]); |
---|
| 1051 | } |
---|
| 1052 | else if (pcDotPhi > lenPhiZ) |
---|
| 1053 | { |
---|
| 1054 | // |
---|
| 1055 | // ...and above in phi. Find distance to point (B) |
---|
| 1056 | // |
---|
| 1057 | G4double distOutPhi = pcDotPhi-lenPhiZ; |
---|
| 1058 | distOut2 = distOutPhi*distOutPhi + distOutZ*distOutZ; |
---|
| 1059 | G4ThreeVector pb = p - vec.edges[1]->corner[0]; |
---|
| 1060 | *normDist = pb.dot(vec.edges[1]->cornNorm[0]); |
---|
| 1061 | } |
---|
| 1062 | else |
---|
| 1063 | { |
---|
| 1064 | // |
---|
| 1065 | // ...and inside in phi. Find distance to line (C) |
---|
| 1066 | // |
---|
| 1067 | G4ThreeVector pa = p - vec.edges[0]->corner[0]; |
---|
| 1068 | distOut2 = distOutZ*distOutZ; |
---|
| 1069 | *normDist = pa.dot(vec.edgeNorm[0]); |
---|
| 1070 | } |
---|
| 1071 | } |
---|
| 1072 | else if (pcDotRZ > lenRZ) |
---|
| 1073 | { |
---|
| 1074 | G4double lenPhiZ = lenPhi[0] + lenRZ*lenPhi[1]; |
---|
| 1075 | G4double distOutZ = pcDotRZ-lenRZ; |
---|
| 1076 | // |
---|
| 1077 | // Above in RZ |
---|
| 1078 | // |
---|
| 1079 | if (pcDotPhi < -lenPhiZ) |
---|
| 1080 | { |
---|
| 1081 | // |
---|
| 1082 | // ...and below in phi. Find distance to point (D) |
---|
| 1083 | // |
---|
| 1084 | G4double distOutPhi = pcDotPhi+lenPhiZ; |
---|
| 1085 | distOut2 = distOutPhi*distOutPhi + distOutZ*distOutZ; |
---|
| 1086 | G4ThreeVector pd = p - vec.edges[0]->corner[1]; |
---|
| 1087 | *normDist = pd.dot(vec.edges[0]->cornNorm[1]); |
---|
| 1088 | } |
---|
| 1089 | else if (pcDotPhi > lenPhiZ) |
---|
| 1090 | { |
---|
| 1091 | // |
---|
| 1092 | // ...and above in phi. Find distance to point (E) |
---|
| 1093 | // |
---|
| 1094 | G4double distOutPhi = pcDotPhi-lenPhiZ; |
---|
| 1095 | distOut2 = distOutPhi*distOutPhi + distOutZ*distOutZ; |
---|
| 1096 | G4ThreeVector pe = p - vec.edges[1]->corner[1]; |
---|
| 1097 | *normDist = pe.dot(vec.edges[1]->cornNorm[1]); |
---|
| 1098 | } |
---|
| 1099 | else |
---|
| 1100 | { |
---|
| 1101 | // |
---|
| 1102 | // ...and inside in phi. Find distance to line (F) |
---|
| 1103 | // |
---|
| 1104 | distOut2 = distOutZ*distOutZ; |
---|
| 1105 | G4ThreeVector pd = p - vec.edges[0]->corner[1]; |
---|
| 1106 | *normDist = pd.dot(vec.edgeNorm[1]); |
---|
| 1107 | } |
---|
| 1108 | } |
---|
| 1109 | else |
---|
| 1110 | { |
---|
| 1111 | G4double lenPhiZ = lenPhi[0] + pcDotRZ*lenPhi[1]; |
---|
| 1112 | // |
---|
| 1113 | // We are inside RZ bounds |
---|
| 1114 | // |
---|
| 1115 | if (pcDotPhi < -lenPhiZ) |
---|
| 1116 | { |
---|
| 1117 | // |
---|
| 1118 | // ...and below in phi. Find distance to line (G) |
---|
| 1119 | // |
---|
| 1120 | G4double distOut = edgeNorm*(pcDotPhi+lenPhiZ); |
---|
| 1121 | distOut2 = distOut*distOut; |
---|
| 1122 | G4ThreeVector pd = p - vec.edges[0]->corner[1]; |
---|
| 1123 | *normDist = pd.dot(vec.edges[0]->normal); |
---|
| 1124 | } |
---|
| 1125 | else if (pcDotPhi > lenPhiZ) |
---|
| 1126 | { |
---|
| 1127 | // |
---|
| 1128 | // ...and above in phi. Find distance to line (H) |
---|
| 1129 | // |
---|
| 1130 | G4double distOut = edgeNorm*(pcDotPhi-lenPhiZ); |
---|
| 1131 | distOut2 = distOut*distOut; |
---|
| 1132 | G4ThreeVector pe = p - vec.edges[1]->corner[1]; |
---|
| 1133 | *normDist = pe.dot(vec.edges[1]->normal); |
---|
| 1134 | } |
---|
| 1135 | else |
---|
| 1136 | { |
---|
| 1137 | // |
---|
| 1138 | // Inside bounds! No penalty. |
---|
| 1139 | // |
---|
| 1140 | return std::fabs(distFaceNorm); |
---|
| 1141 | } |
---|
| 1142 | } |
---|
| 1143 | return std::sqrt( distFaceNorm*distFaceNorm + distOut2 ); |
---|
| 1144 | } |
---|
[850] | 1145 | |
---|
| 1146 | |
---|
| 1147 | // |
---|
| 1148 | // Calculation of surface area of a triangle. |
---|
| 1149 | // At the same time a random point in the triangle is given |
---|
| 1150 | // |
---|
| 1151 | G4double G4PolyhedraSide::SurfaceTriangle( G4ThreeVector p1, |
---|
| 1152 | G4ThreeVector p2, |
---|
| 1153 | G4ThreeVector p3, |
---|
| 1154 | G4ThreeVector *p4 ) |
---|
| 1155 | { |
---|
| 1156 | G4ThreeVector v, w; |
---|
| 1157 | |
---|
| 1158 | v = p3 - p1; |
---|
| 1159 | w = p1 - p2; |
---|
| 1160 | G4double lambda1 = G4UniformRand(); |
---|
| 1161 | G4double lambda2 = lambda1*G4UniformRand(); |
---|
| 1162 | |
---|
| 1163 | *p4=p2 + lambda1*w + lambda2*v; |
---|
| 1164 | return 0.5*(v.cross(w)).mag(); |
---|
| 1165 | } |
---|
| 1166 | |
---|
| 1167 | |
---|
| 1168 | // |
---|
| 1169 | // GetPointOnPlane |
---|
| 1170 | // |
---|
| 1171 | // Auxiliary method for GetPointOnSurface() |
---|
| 1172 | // |
---|
| 1173 | G4ThreeVector |
---|
| 1174 | G4PolyhedraSide::GetPointOnPlane( G4ThreeVector p0, G4ThreeVector p1, |
---|
| 1175 | G4ThreeVector p2, G4ThreeVector p3, |
---|
| 1176 | G4double *Area ) |
---|
| 1177 | { |
---|
| 1178 | G4double chose,aOne,aTwo; |
---|
| 1179 | G4ThreeVector point1,point2; |
---|
| 1180 | aOne = SurfaceTriangle(p0,p1,p2,&point1); |
---|
| 1181 | aTwo = SurfaceTriangle(p2,p3,p0,&point2); |
---|
| 1182 | *Area= aOne+aTwo; |
---|
| 1183 | |
---|
| 1184 | chose = G4UniformRand()*(aOne+aTwo); |
---|
| 1185 | if( (chose>=0.) && (chose < aOne) ) |
---|
| 1186 | { |
---|
| 1187 | return (point1); |
---|
| 1188 | } |
---|
| 1189 | return (point2); |
---|
| 1190 | } |
---|
| 1191 | |
---|
| 1192 | |
---|
| 1193 | // |
---|
| 1194 | // SurfaceArea() |
---|
| 1195 | // |
---|
| 1196 | G4double G4PolyhedraSide::SurfaceArea() |
---|
| 1197 | { |
---|
| 1198 | if( fSurfaceArea==0. ) |
---|
| 1199 | { |
---|
| 1200 | // Define the variables |
---|
| 1201 | // |
---|
| 1202 | G4double area,areas; |
---|
| 1203 | G4ThreeVector point1; |
---|
| 1204 | G4ThreeVector v1,v2,v3,v4; |
---|
| 1205 | G4PolyhedraSideVec *vec = vecs; |
---|
| 1206 | areas=0.; |
---|
| 1207 | |
---|
| 1208 | // Do a loop on all SideEdge |
---|
| 1209 | // |
---|
| 1210 | do |
---|
| 1211 | { |
---|
| 1212 | // Define 4points for a Plane or Triangle |
---|
| 1213 | // |
---|
| 1214 | G4ThreeVector v1=vec->edges[0]->corner[0]; |
---|
| 1215 | G4ThreeVector v2=vec->edges[0]->corner[1]; |
---|
| 1216 | G4ThreeVector v3=vec->edges[1]->corner[1]; |
---|
| 1217 | G4ThreeVector v4=vec->edges[1]->corner[0]; |
---|
| 1218 | point1=GetPointOnPlane(v1,v2,v3,v4,&area); |
---|
| 1219 | areas+=area; |
---|
| 1220 | } while( ++vec < vecs + numSide); |
---|
| 1221 | |
---|
| 1222 | fSurfaceArea=areas; |
---|
| 1223 | } |
---|
| 1224 | return fSurfaceArea; |
---|
| 1225 | } |
---|
| 1226 | |
---|
| 1227 | |
---|
| 1228 | // |
---|
| 1229 | // GetPointOnFace() |
---|
| 1230 | // |
---|
| 1231 | G4ThreeVector G4PolyhedraSide::GetPointOnFace() |
---|
| 1232 | { |
---|
| 1233 | // Define the variables |
---|
| 1234 | // |
---|
| 1235 | std::vector<G4double>areas; |
---|
| 1236 | std::vector<G4ThreeVector>points; |
---|
| 1237 | G4double area=0; |
---|
| 1238 | G4double result1; |
---|
| 1239 | G4ThreeVector point1; |
---|
| 1240 | G4ThreeVector v1,v2,v3,v4; |
---|
| 1241 | G4PolyhedraSideVec *vec = vecs; |
---|
| 1242 | |
---|
| 1243 | // Do a loop on all SideEdge |
---|
| 1244 | // |
---|
| 1245 | do |
---|
| 1246 | { |
---|
| 1247 | // Define 4points for a Plane or Triangle |
---|
| 1248 | // |
---|
| 1249 | G4ThreeVector v1=vec->edges[0]->corner[0]; |
---|
| 1250 | G4ThreeVector v2=vec->edges[0]->corner[1]; |
---|
| 1251 | G4ThreeVector v3=vec->edges[1]->corner[1]; |
---|
| 1252 | G4ThreeVector v4=vec->edges[1]->corner[0]; |
---|
| 1253 | point1=GetPointOnPlane(v1,v2,v3,v4,&result1); |
---|
| 1254 | points.push_back(point1); |
---|
| 1255 | areas.push_back(result1); |
---|
| 1256 | area+=result1; |
---|
| 1257 | } while( ++vec < vecs+numSide ); |
---|
| 1258 | |
---|
| 1259 | // Choose randomly one of the surfaces and point on it |
---|
| 1260 | // |
---|
| 1261 | G4double chose = area*G4UniformRand(); |
---|
| 1262 | G4double Achose1,Achose2; |
---|
| 1263 | Achose1=0;Achose2=0.; |
---|
| 1264 | G4int i=0; |
---|
| 1265 | do |
---|
| 1266 | { |
---|
| 1267 | Achose2+=areas[i]; |
---|
| 1268 | if(chose>=Achose1 && chose<Achose2) |
---|
| 1269 | { |
---|
| 1270 | point1=points[i] ; break; |
---|
| 1271 | } |
---|
| 1272 | i++; Achose1=Achose2; |
---|
| 1273 | } while( i<numSide ); |
---|
| 1274 | |
---|
| 1275 | return point1; |
---|
| 1276 | } |
---|