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