| 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 G4PolyPhiFace, the face t <<
27 // polyhedra at its phi opening. <<
28 // 26 //
29 // Author: David C. Williams (davidw@scipp.ucs << 27 // $Id: G4PolyPhiFace.cc 95301 2016-02-04 11:25:33Z gcosmo $
>> 28 //
>> 29 //
>> 30 // --------------------------------------------------------------------
>> 31 // GEANT 4 class source file
>> 32 //
>> 33 //
>> 34 // G4PolyPhiFace.cc
>> 35 //
>> 36 // Implementation of the face that bounds a polycone or polyhedra at
>> 37 // its phi opening.
>> 38 //
30 // ------------------------------------------- 39 // --------------------------------------------------------------------
31 40
32 #include "G4PolyPhiFace.hh" 41 #include "G4PolyPhiFace.hh"
33 #include "G4ClippablePolygon.hh" 42 #include "G4ClippablePolygon.hh"
34 #include "G4ReduciblePolygon.hh" 43 #include "G4ReduciblePolygon.hh"
35 #include "G4AffineTransform.hh" 44 #include "G4AffineTransform.hh"
36 #include "G4SolidExtentList.hh" 45 #include "G4SolidExtentList.hh"
37 #include "G4GeometryTolerance.hh" 46 #include "G4GeometryTolerance.hh"
38 47
39 #include "Randomize.hh" 48 #include "Randomize.hh"
40 #include "G4TwoVector.hh" 49 #include "G4TwoVector.hh"
41 50
>> 51 //
42 // Constructor 52 // Constructor
43 // 53 //
44 // Points r,z should be supplied in clockwise 54 // Points r,z should be supplied in clockwise order in r,z. For example:
45 // 55 //
46 // [1]---------[2] ^ R 56 // [1]---------[2] ^ R
47 // | | | 57 // | | |
48 // | | +--> 58 // | | +--> z
49 // [0]---------[3] 59 // [0]---------[3]
50 // 60 //
51 G4PolyPhiFace::G4PolyPhiFace( const G4Reducibl << 61 G4PolyPhiFace::G4PolyPhiFace( const G4ReduciblePolygon *rz,
52 G4double p 62 G4double phi,
53 G4double d 63 G4double deltaPhi,
54 G4double p 64 G4double phiOther )
>> 65 : fSurfaceArea(0.), triangles(0)
55 { 66 {
56 kCarTolerance = G4GeometryTolerance::GetInst 67 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
57 68
58 numEdges = rz->NumVertices(); 69 numEdges = rz->NumVertices();
59 << 70
60 rMin = rz->Amin(); 71 rMin = rz->Amin();
61 rMax = rz->Amax(); 72 rMax = rz->Amax();
62 zMin = rz->Bmin(); 73 zMin = rz->Bmin();
63 zMax = rz->Bmax(); 74 zMax = rz->Bmax();
64 75
65 // 76 //
66 // Is this the "starting" phi edge of the tw 77 // Is this the "starting" phi edge of the two?
67 // 78 //
68 G4bool start = (phiOther > phi); 79 G4bool start = (phiOther > phi);
69 80
70 // 81 //
71 // Build radial vector 82 // Build radial vector
72 // 83 //
73 radial = G4ThreeVector( std::cos(phi), std:: 84 radial = G4ThreeVector( std::cos(phi), std::sin(phi), 0.0 );
74 85
75 // 86 //
76 // Build normal 87 // Build normal
77 // 88 //
78 G4double zSign = start ? 1 : -1; 89 G4double zSign = start ? 1 : -1;
79 normal = G4ThreeVector( zSign*radial.y(), -z 90 normal = G4ThreeVector( zSign*radial.y(), -zSign*radial.x(), 0 );
80 91
81 // 92 //
82 // Is allBehind? 93 // Is allBehind?
83 // 94 //
84 allBehind = (zSign*(std::cos(phiOther)*radia 95 allBehind = (zSign*(std::cos(phiOther)*radial.y() - std::sin(phiOther)*radial.x()) < 0);
85 96
86 // 97 //
87 // Adjacent edges 98 // Adjacent edges
88 // 99 //
89 G4double midPhi = phi + (start ? +0.5 : -0.5 100 G4double midPhi = phi + (start ? +0.5 : -0.5)*deltaPhi;
90 G4double cosMid = std::cos(midPhi), 101 G4double cosMid = std::cos(midPhi),
91 sinMid = std::sin(midPhi); 102 sinMid = std::sin(midPhi);
>> 103
92 // 104 //
93 // Allocate corners 105 // Allocate corners
94 // 106 //
95 const std::size_t maxEdges = numEdges>0 ? nu << 107 corners = new G4PolyPhiFaceVertex[numEdges];
96 corners = new G4PolyPhiFaceVertex[maxEdges]; <<
97 // 108 //
98 // Fill them 109 // Fill them
99 // 110 //
100 G4ReduciblePolygonIterator iterRZ(rz); 111 G4ReduciblePolygonIterator iterRZ(rz);
101 112
102 G4PolyPhiFaceVertex* corn = corners; << 113 G4PolyPhiFaceVertex *corn = corners;
103 G4PolyPhiFaceVertex* helper = corners; << 114 G4PolyPhiFaceVertex *helper=corners;
104 115
105 iterRZ.Begin(); 116 iterRZ.Begin();
106 do // Loop checking, 13.08.2015, G.Cosmo 117 do // Loop checking, 13.08.2015, G.Cosmo
107 { 118 {
108 corn->r = iterRZ.GetA(); 119 corn->r = iterRZ.GetA();
109 corn->z = iterRZ.GetB(); 120 corn->z = iterRZ.GetB();
110 corn->x = corn->r*radial.x(); 121 corn->x = corn->r*radial.x();
111 corn->y = corn->r*radial.y(); 122 corn->y = corn->r*radial.y();
112 123
113 // Add pointer on prev corner 124 // Add pointer on prev corner
114 // 125 //
115 if( corn == corners ) 126 if( corn == corners )
116 { corn->prev = corners+maxEdges-1;} << 127 { corn->prev = corners+numEdges-1;}
117 else 128 else
118 { corn->prev = helper; } 129 { corn->prev = helper; }
119 130
120 // Add pointer on next corner 131 // Add pointer on next corner
121 // 132 //
122 if( corn < corners+maxEdges-1 ) << 133 if( corn < corners+numEdges-1 )
123 { corn->next = corn+1;} 134 { corn->next = corn+1;}
124 else 135 else
125 { corn->next = corners; } 136 { corn->next = corners; }
126 137
127 helper = corn; 138 helper = corn;
128 } while( ++corn, iterRZ.Next() ); 139 } while( ++corn, iterRZ.Next() );
129 140
130 // 141 //
131 // Allocate edges 142 // Allocate edges
132 // 143 //
133 edges = new G4PolyPhiFaceEdge[maxEdges]; << 144 edges = new G4PolyPhiFaceEdge[numEdges];
134 145
135 // 146 //
136 // Fill them 147 // Fill them
137 // 148 //
138 G4double rFact = std::cos(0.5*deltaPhi); 149 G4double rFact = std::cos(0.5*deltaPhi);
139 G4double rFactNormalize = 1.0/std::sqrt(1.0+ 150 G4double rFactNormalize = 1.0/std::sqrt(1.0+rFact*rFact);
140 151
141 G4PolyPhiFaceVertex* prev = corners+maxEdges << 152 G4PolyPhiFaceVertex *prev = corners+numEdges-1,
142 * here = corners; << 153 *here = corners;
143 G4PolyPhiFaceEdge* edge = edges; << 154 G4PolyPhiFaceEdge *edge = edges;
144 do // Loop checking, 13.08.2015, G.Cosmo 155 do // Loop checking, 13.08.2015, G.Cosmo
145 { 156 {
146 G4ThreeVector sideNorm; 157 G4ThreeVector sideNorm;
147 158
148 edge->v0 = prev; 159 edge->v0 = prev;
149 edge->v1 = here; 160 edge->v1 = here;
150 161
151 G4double dr = here->r - prev->r, 162 G4double dr = here->r - prev->r,
152 dz = here->z - prev->z; 163 dz = here->z - prev->z;
153 164
154 edge->length = std::sqrt( dr*dr + dz*dz ); 165 edge->length = std::sqrt( dr*dr + dz*dz );
155 166
156 edge->tr = dr/edge->length; 167 edge->tr = dr/edge->length;
157 edge->tz = dz/edge->length; 168 edge->tz = dz/edge->length;
158 169
159 if ((here->r < DBL_MIN) && (prev->r < DBL_ 170 if ((here->r < DBL_MIN) && (prev->r < DBL_MIN))
160 { 171 {
161 // 172 //
162 // Sigh! Always exceptions! 173 // Sigh! Always exceptions!
163 // This edge runs at r==0, so its adjoin 174 // This edge runs at r==0, so its adjoing surface is not a
164 // PolyconeSide or PolyhedraSide, but th 175 // PolyconeSide or PolyhedraSide, but the opposite PolyPhiFace.
165 // 176 //
166 G4double zSignOther = start ? -1 : 1; 177 G4double zSignOther = start ? -1 : 1;
167 sideNorm = G4ThreeVector( zSignOther*st 178 sideNorm = G4ThreeVector( zSignOther*std::sin(phiOther),
168 -zSignOther*st 179 -zSignOther*std::cos(phiOther), 0 );
169 } 180 }
170 else 181 else
171 { 182 {
172 sideNorm = G4ThreeVector( edge->tz*cosMi 183 sideNorm = G4ThreeVector( edge->tz*cosMid,
173 edge->tz*sinMi 184 edge->tz*sinMid,
174 -edge->tr*rFact 185 -edge->tr*rFact );
175 sideNorm *= rFactNormalize; 186 sideNorm *= rFactNormalize;
176 } 187 }
177 sideNorm += normal; 188 sideNorm += normal;
178 189
179 edge->norm3D = sideNorm.unit(); 190 edge->norm3D = sideNorm.unit();
180 } while( edge++, prev=here, ++here < corners << 191 } while( edge++, prev=here, ++here < corners+numEdges );
181 192
182 // 193 //
183 // Go back and fill in corner "normals" 194 // Go back and fill in corner "normals"
184 // 195 //
185 G4PolyPhiFaceEdge* prevEdge = edges+maxEdges << 196 G4PolyPhiFaceEdge *prevEdge = edges+numEdges-1;
186 edge = edges; 197 edge = edges;
187 do // Loop checking, 13.08.2015, G.Cosmo 198 do // Loop checking, 13.08.2015, G.Cosmo
188 { 199 {
189 // 200 //
190 // Calculate vertex 2D normals (on the phi 201 // Calculate vertex 2D normals (on the phi surface)
191 // 202 //
192 G4double rPart = prevEdge->tr + edge->tr; 203 G4double rPart = prevEdge->tr + edge->tr;
193 G4double zPart = prevEdge->tz + edge->tz; 204 G4double zPart = prevEdge->tz + edge->tz;
194 G4double norm = std::sqrt( rPart*rPart + z 205 G4double norm = std::sqrt( rPart*rPart + zPart*zPart );
195 G4double rNorm = +zPart/norm; 206 G4double rNorm = +zPart/norm;
196 G4double zNorm = -rPart/norm; 207 G4double zNorm = -rPart/norm;
197 208
198 edge->v0->rNorm = rNorm; 209 edge->v0->rNorm = rNorm;
199 edge->v0->zNorm = zNorm; 210 edge->v0->zNorm = zNorm;
200 211
201 // 212 //
202 // Calculate the 3D normals. 213 // Calculate the 3D normals.
203 // 214 //
204 // Find the vector perpendicular to the z 215 // Find the vector perpendicular to the z axis
205 // that defines the plane that contains th 216 // that defines the plane that contains the vertex normal
206 // 217 //
207 G4ThreeVector xyVector; 218 G4ThreeVector xyVector;
208 219
209 if (edge->v0->r < DBL_MIN) 220 if (edge->v0->r < DBL_MIN)
210 { 221 {
211 // 222 //
212 // This is a vertex at r==0, which is a 223 // This is a vertex at r==0, which is a special
213 // case. The normal we will construct la 224 // case. The normal we will construct lays in the
214 // plane at the center of the phi openin 225 // plane at the center of the phi opening.
215 // 226 //
216 // We also know that rNorm < 0 227 // We also know that rNorm < 0
217 // 228 //
218 G4double zSignOther = start ? -1 : 1; 229 G4double zSignOther = start ? -1 : 1;
219 G4ThreeVector normalOther( zSignOther*s 230 G4ThreeVector normalOther( zSignOther*std::sin(phiOther),
220 -zSignOther*s 231 -zSignOther*std::cos(phiOther), 0 );
221 232
222 xyVector = - normal - normalOther; 233 xyVector = - normal - normalOther;
223 } 234 }
224 else 235 else
225 { 236 {
226 // 237 //
227 // This is a vertex at r > 0. The plane 238 // This is a vertex at r > 0. The plane
228 // is the average of the normal and the 239 // is the average of the normal and the
229 // normal of the adjacent phi face 240 // normal of the adjacent phi face
230 // 241 //
231 xyVector = G4ThreeVector( cosMid, sinMid 242 xyVector = G4ThreeVector( cosMid, sinMid, 0 );
232 if (rNorm < 0) 243 if (rNorm < 0)
233 xyVector -= normal; 244 xyVector -= normal;
234 else 245 else
235 xyVector += normal; 246 xyVector += normal;
236 } 247 }
237 248
238 // 249 //
239 // Combine it with the r/z direction from 250 // Combine it with the r/z direction from the face
240 // 251 //
241 edge->v0->norm3D = rNorm*xyVector.unit() + 252 edge->v0->norm3D = rNorm*xyVector.unit() + G4ThreeVector( 0, 0, zNorm );
242 } while( prevEdge=edge, ++edge < edges+maxE << 253 } while( prevEdge=edge, ++edge < edges+numEdges );
243 254
244 // 255 //
245 // Build point on surface 256 // Build point on surface
246 // 257 //
247 G4double rAve = 0.5*(rMax-rMin), 258 G4double rAve = 0.5*(rMax-rMin),
248 zAve = 0.5*(zMax-zMin); 259 zAve = 0.5*(zMax-zMin);
249 surface = G4ThreeVector( rAve*radial.x(), rA 260 surface = G4ThreeVector( rAve*radial.x(), rAve*radial.y(), zAve );
250 } 261 }
251 262
>> 263
>> 264 //
252 // Diagnose 265 // Diagnose
253 // 266 //
254 // Throw an exception if something is found in 267 // Throw an exception if something is found inconsistent with
255 // the solid. 268 // the solid.
256 // 269 //
257 // For debugging purposes only 270 // For debugging purposes only
258 // 271 //
259 void G4PolyPhiFace::Diagnose( G4VSolid* owner << 272 void G4PolyPhiFace::Diagnose( G4VSolid *owner )
260 { 273 {
261 G4PolyPhiFaceVertex* corner = corners; << 274 G4PolyPhiFaceVertex *corner = corners;
262 do // Loop checking, 13.08.2015, G.Cosmo 275 do // Loop checking, 13.08.2015, G.Cosmo
263 { 276 {
264 G4ThreeVector test(corner->x, corner->y, c 277 G4ThreeVector test(corner->x, corner->y, corner->z);
265 test -= 1E-6*corner->norm3D; 278 test -= 1E-6*corner->norm3D;
266 279
267 if (owner->Inside(test) != kInside) 280 if (owner->Inside(test) != kInside)
268 G4Exception( "G4PolyPhiFace::Diagnose()" 281 G4Exception( "G4PolyPhiFace::Diagnose()", "GeomSolids0002",
269 FatalException, "Bad vertex 282 FatalException, "Bad vertex normal found." );
270 } while( ++corner < corners+numEdges ); 283 } while( ++corner < corners+numEdges );
271 } 284 }
272 285
>> 286
>> 287 //
273 // Fake default constructor - sets only member 288 // Fake default constructor - sets only member data and allocates memory
274 // for usage restri 289 // for usage restricted to object persistency.
275 // 290 //
276 G4PolyPhiFace::G4PolyPhiFace( __void__&) 291 G4PolyPhiFace::G4PolyPhiFace( __void__&)
277 : numEdges(0), rMin(0.), rMax(0.), zMin(0.), << 292 : numEdges(0), edges(0), corners(0), rMin(0.), rMax(0.), zMin(0.), zMax(0.),
>> 293 allBehind(false), kCarTolerance(0.), fSurfaceArea(0.), triangles(0)
278 { 294 {
279 } 295 }
280 296
281 297
282 // 298 //
283 // Destructor 299 // Destructor
284 // 300 //
285 G4PolyPhiFace::~G4PolyPhiFace() 301 G4PolyPhiFace::~G4PolyPhiFace()
286 { 302 {
287 delete [] edges; 303 delete [] edges;
288 delete [] corners; 304 delete [] corners;
289 } 305 }
290 306
291 307
292 // 308 //
293 // Copy constructor 309 // Copy constructor
294 // 310 //
295 G4PolyPhiFace::G4PolyPhiFace( const G4PolyPhiF << 311 G4PolyPhiFace::G4PolyPhiFace( const G4PolyPhiFace &source )
>> 312 : G4VCSGface()
296 { 313 {
297 CopyStuff( source ); 314 CopyStuff( source );
298 } 315 }
299 316
300 317
301 // 318 //
302 // Assignment operator 319 // Assignment operator
303 // 320 //
304 G4PolyPhiFace& G4PolyPhiFace::operator=( const << 321 G4PolyPhiFace& G4PolyPhiFace::operator=( const G4PolyPhiFace &source )
305 { 322 {
306 if (this == &source) { return *this; } 323 if (this == &source) { return *this; }
307 324
308 delete [] edges; 325 delete [] edges;
309 delete [] corners; 326 delete [] corners;
310 327
311 CopyStuff( source ); 328 CopyStuff( source );
312 329
313 return *this; 330 return *this;
314 } 331 }
315 332
316 333
317 // 334 //
318 // CopyStuff (protected) 335 // CopyStuff (protected)
319 // 336 //
320 void G4PolyPhiFace::CopyStuff( const G4PolyPhi << 337 void G4PolyPhiFace::CopyStuff( const G4PolyPhiFace &source )
321 { 338 {
322 // 339 //
323 // The simple stuff 340 // The simple stuff
324 // 341 //
325 numEdges = source.numEdges; 342 numEdges = source.numEdges;
326 normal = source.normal; 343 normal = source.normal;
327 radial = source.radial; 344 radial = source.radial;
328 surface = source.surface; 345 surface = source.surface;
329 rMin = source.rMin; 346 rMin = source.rMin;
330 rMax = source.rMax; 347 rMax = source.rMax;
331 zMin = source.zMin; 348 zMin = source.zMin;
332 zMax = source.zMax; 349 zMax = source.zMax;
333 allBehind = source.allBehind; 350 allBehind = source.allBehind;
334 triangles = nullptr; << 351 triangles = 0;
335 352
336 kCarTolerance = source.kCarTolerance; 353 kCarTolerance = source.kCarTolerance;
337 fSurfaceArea = source.fSurfaceArea; 354 fSurfaceArea = source.fSurfaceArea;
338 355
339 const std::size_t maxEdges = (numEdges > 0) <<
340 <<
341 // 356 //
342 // Corner dynamic array 357 // Corner dynamic array
343 // 358 //
344 corners = new G4PolyPhiFaceVertex[maxEdges]; << 359 corners = new G4PolyPhiFaceVertex[numEdges];
345 G4PolyPhiFaceVertex *corn = corners, 360 G4PolyPhiFaceVertex *corn = corners,
346 *sourceCorn = source.cor 361 *sourceCorn = source.corners;
347 do // Loop checking, 13.08.2015, G.Cosmo 362 do // Loop checking, 13.08.2015, G.Cosmo
348 { 363 {
349 *corn = *sourceCorn; 364 *corn = *sourceCorn;
350 } while( ++sourceCorn, ++corn < corners+maxE << 365 } while( ++sourceCorn, ++corn < corners+numEdges );
351 366
352 // 367 //
353 // Edge dynamic array 368 // Edge dynamic array
354 // 369 //
355 edges = new G4PolyPhiFaceEdge[maxEdges]; << 370 edges = new G4PolyPhiFaceEdge[numEdges];
356 371
357 G4PolyPhiFaceVertex* prev = corners+maxEdges << 372 G4PolyPhiFaceVertex *prev = corners+numEdges-1,
358 * here = corners; << 373 *here = corners;
359 G4PolyPhiFaceEdge* edge = edges, << 374 G4PolyPhiFaceEdge *edge = edges,
360 * sourceEdge = source.edg << 375 *sourceEdge = source.edges;
361 do // Loop checking, 13.08.2015, G.Cosmo 376 do // Loop checking, 13.08.2015, G.Cosmo
362 { 377 {
363 *edge = *sourceEdge; 378 *edge = *sourceEdge;
364 edge->v0 = prev; 379 edge->v0 = prev;
365 edge->v1 = here; 380 edge->v1 = here;
366 } while( ++sourceEdge, ++edge, prev=here, ++ << 381 } while( ++sourceEdge, ++edge, prev=here, ++here < corners+numEdges );
367 } 382 }
368 383
>> 384
>> 385 //
369 // Intersect 386 // Intersect
370 // 387 //
371 G4bool G4PolyPhiFace::Intersect( const G4Three << 388 G4bool G4PolyPhiFace::Intersect( const G4ThreeVector &p,
372 const G4Three << 389 const G4ThreeVector &v,
373 G4bool 390 G4bool outgoing,
374 G4doubl 391 G4double surfTolerance,
375 G4doubl << 392 G4double &distance,
376 G4doubl << 393 G4double &distFromSurface,
377 G4Three << 394 G4ThreeVector &aNormal,
378 G4bool& << 395 G4bool &isAllBehind )
379 { 396 {
380 G4double normSign = outgoing ? +1 : -1; 397 G4double normSign = outgoing ? +1 : -1;
381 398
382 // 399 //
383 // These don't change 400 // These don't change
384 // 401 //
385 isAllBehind = allBehind; 402 isAllBehind = allBehind;
386 aNormal = normal; 403 aNormal = normal;
387 404
388 // 405 //
389 // Correct normal? Here we have straight sid 406 // Correct normal? Here we have straight sides, and can safely ignore
390 // intersections where the dot product with 407 // intersections where the dot product with the normal is zero.
391 // 408 //
392 G4double dotProd = normSign*normal.dot(v); 409 G4double dotProd = normSign*normal.dot(v);
393 410
394 if (dotProd <= 0) return false; 411 if (dotProd <= 0) return false;
395 412
396 // 413 //
397 // Calculate distance to surface. If the sid 414 // Calculate distance to surface. If the side is too far
398 // behind the point, we must reject it. 415 // behind the point, we must reject it.
399 // 416 //
400 G4ThreeVector ps = p - surface; 417 G4ThreeVector ps = p - surface;
401 distFromSurface = -normSign*ps.dot(normal); 418 distFromSurface = -normSign*ps.dot(normal);
402 419
403 if (distFromSurface < -surfTolerance) return 420 if (distFromSurface < -surfTolerance) return false;
404 421
405 // 422 //
406 // Calculate precise distance to intersectio 423 // Calculate precise distance to intersection with the side
407 // (along the trajectory, not normal to the 424 // (along the trajectory, not normal to the surface)
408 // 425 //
409 distance = distFromSurface/dotProd; 426 distance = distFromSurface/dotProd;
410 427
411 // 428 //
412 // Calculate intersection point in r,z 429 // Calculate intersection point in r,z
413 // 430 //
414 G4ThreeVector ip = p + distance*v; 431 G4ThreeVector ip = p + distance*v;
415 432
416 G4double r = radial.dot(ip); 433 G4double r = radial.dot(ip);
417 434
418 // 435 //
419 // And is it inside the r/z extent? 436 // And is it inside the r/z extent?
420 // 437 //
421 return InsideEdgesExact( r, ip.z(), normSign 438 return InsideEdgesExact( r, ip.z(), normSign, p, v );
422 } 439 }
423 440
>> 441
>> 442 //
424 // Distance 443 // Distance
425 // 444 //
426 G4double G4PolyPhiFace::Distance( const G4Thre << 445 G4double G4PolyPhiFace::Distance( const G4ThreeVector &p, G4bool outgoing )
427 { 446 {
428 G4double normSign = outgoing ? +1 : -1; 447 G4double normSign = outgoing ? +1 : -1;
429 // 448 //
430 // Correct normal? 449 // Correct normal?
431 // 450 //
432 G4ThreeVector ps = p - surface; 451 G4ThreeVector ps = p - surface;
433 G4double distPhi = -normSign*normal.dot(ps); 452 G4double distPhi = -normSign*normal.dot(ps);
434 453
435 if (distPhi < -0.5*kCarTolerance) 454 if (distPhi < -0.5*kCarTolerance)
436 return kInfinity; 455 return kInfinity;
437 else if (distPhi < 0) 456 else if (distPhi < 0)
438 distPhi = 0.0; 457 distPhi = 0.0;
439 458
440 // 459 //
441 // Calculate projected point in r,z 460 // Calculate projected point in r,z
442 // 461 //
443 G4double r = radial.dot(p); 462 G4double r = radial.dot(p);
444 463
445 // 464 //
446 // Are we inside the face? 465 // Are we inside the face?
447 // 466 //
448 G4double distRZ2; 467 G4double distRZ2;
449 468
450 if (InsideEdges( r, p.z(), &distRZ2, nullptr << 469 if (InsideEdges( r, p.z(), &distRZ2, 0 ))
451 { 470 {
452 // 471 //
453 // Yup, answer is just distPhi 472 // Yup, answer is just distPhi
454 // 473 //
455 return distPhi; 474 return distPhi;
456 } 475 }
457 else 476 else
458 { 477 {
459 // 478 //
460 // Nope. Penalize by distance out 479 // Nope. Penalize by distance out
461 // 480 //
462 return std::sqrt( distPhi*distPhi + distRZ 481 return std::sqrt( distPhi*distPhi + distRZ2 );
463 } 482 }
464 } 483 }
>> 484
465 485
>> 486 //
466 // Inside 487 // Inside
467 // 488 //
468 EInside G4PolyPhiFace::Inside( const G4ThreeVe << 489 EInside G4PolyPhiFace::Inside( const G4ThreeVector &p,
469 G4double 490 G4double tolerance,
470 G4double* << 491 G4double *bestDistance )
471 { 492 {
472 // 493 //
473 // Get distance along phi, which if negative 494 // Get distance along phi, which if negative means the point
474 // is nominally inside the shape. 495 // is nominally inside the shape.
475 // 496 //
476 G4ThreeVector ps = p - surface; 497 G4ThreeVector ps = p - surface;
477 G4double distPhi = normal.dot(ps); 498 G4double distPhi = normal.dot(ps);
478 499
479 // 500 //
480 // Calculate projected point in r,z 501 // Calculate projected point in r,z
481 // 502 //
482 G4double r = radial.dot(p); 503 G4double r = radial.dot(p);
483 504
484 // 505 //
485 // Are we inside the face? 506 // Are we inside the face?
486 // 507 //
487 G4double distRZ2; 508 G4double distRZ2;
488 G4PolyPhiFaceVertex* base3Dnorm = nullptr; << 509 G4PolyPhiFaceVertex *base3Dnorm;
489 G4ThreeVector* head3Dnorm = nullptr; << 510 G4ThreeVector *head3Dnorm;
490 511
491 if (InsideEdges( r, p.z(), &distRZ2, &base3D 512 if (InsideEdges( r, p.z(), &distRZ2, &base3Dnorm, &head3Dnorm ))
492 { 513 {
493 // 514 //
494 // Looks like we're inside. Distance is di 515 // Looks like we're inside. Distance is distance in phi.
495 // 516 //
496 *bestDistance = std::fabs(distPhi); 517 *bestDistance = std::fabs(distPhi);
497 518
498 // 519 //
499 // Use distPhi to decide fate 520 // Use distPhi to decide fate
500 // 521 //
501 if (distPhi < -tolerance) return kInside; 522 if (distPhi < -tolerance) return kInside;
502 if (distPhi < tolerance) return kSurface; 523 if (distPhi < tolerance) return kSurface;
503 return kOutside; 524 return kOutside;
504 } 525 }
505 else 526 else
506 { 527 {
507 // 528 //
508 // We're outside the extent of the face, 529 // We're outside the extent of the face,
509 // so the distance is penalized by distanc 530 // so the distance is penalized by distance from edges in RZ
510 // 531 //
511 *bestDistance = std::sqrt( distPhi*distPhi 532 *bestDistance = std::sqrt( distPhi*distPhi + distRZ2 );
512 533
513 // 534 //
514 // Use edge normal to decide fate 535 // Use edge normal to decide fate
515 // 536 //
516 G4ThreeVector cc( base3Dnorm->r*radial.x() 537 G4ThreeVector cc( base3Dnorm->r*radial.x(),
517 base3Dnorm->r*radial.y() << 538 base3Dnorm->r*radial.y(),
518 base3Dnorm->z ); << 539 base3Dnorm->z );
519 cc = p - cc; 540 cc = p - cc;
520 G4double normDist = head3Dnorm->dot(cc); 541 G4double normDist = head3Dnorm->dot(cc);
521 if ( distRZ2 > tolerance*tolerance ) 542 if ( distRZ2 > tolerance*tolerance )
522 { 543 {
523 // 544 //
524 // We're far enough away that kSurface i 545 // We're far enough away that kSurface is not possible
525 // 546 //
526 return normDist < 0 ? kInside : kOutside 547 return normDist < 0 ? kInside : kOutside;
527 } 548 }
528 549
529 if (normDist < -tolerance) return kInside; 550 if (normDist < -tolerance) return kInside;
530 if (normDist < tolerance) return kSurface 551 if (normDist < tolerance) return kSurface;
531 return kOutside; 552 return kOutside;
532 } 553 }
533 } 554 }
534 555
>> 556
>> 557 //
535 // Normal 558 // Normal
536 // 559 //
537 // This virtual member is simple for our plane 560 // This virtual member is simple for our planer shape,
538 // which has only one normal 561 // which has only one normal
539 // 562 //
540 G4ThreeVector G4PolyPhiFace::Normal( const G4T << 563 G4ThreeVector G4PolyPhiFace::Normal( const G4ThreeVector &p,
541 G4d << 564 G4double *bestDistance )
542 { 565 {
543 // 566 //
544 // Get distance along phi, which if negative 567 // Get distance along phi, which if negative means the point
545 // is nominally inside the shape. 568 // is nominally inside the shape.
546 // 569 //
547 G4double distPhi = normal.dot(p); 570 G4double distPhi = normal.dot(p);
548 571
549 // 572 //
550 // Calculate projected point in r,z 573 // Calculate projected point in r,z
551 // 574 //
552 G4double r = radial.dot(p); 575 G4double r = radial.dot(p);
553 576
554 // 577 //
555 // Are we inside the face? 578 // Are we inside the face?
556 // 579 //
557 G4double distRZ2; 580 G4double distRZ2;
558 581
559 if (InsideEdges( r, p.z(), &distRZ2, nullptr << 582 if (InsideEdges( r, p.z(), &distRZ2, 0 ))
560 { 583 {
561 // 584 //
562 // Yup, answer is just distPhi 585 // Yup, answer is just distPhi
563 // 586 //
564 *bestDistance = std::fabs(distPhi); 587 *bestDistance = std::fabs(distPhi);
565 } 588 }
566 else 589 else
567 { 590 {
568 // 591 //
569 // Nope. Penalize by distance out 592 // Nope. Penalize by distance out
570 // 593 //
571 *bestDistance = std::sqrt( distPhi*distPhi 594 *bestDistance = std::sqrt( distPhi*distPhi + distRZ2 );
572 } 595 }
573 596
574 return normal; 597 return normal;
575 } 598 }
576 599
>> 600
>> 601 //
577 // Extent 602 // Extent
578 // 603 //
579 // This actually isn't needed by polycone or p 604 // This actually isn't needed by polycone or polyhedra...
580 // 605 //
581 G4double G4PolyPhiFace::Extent( const G4ThreeV 606 G4double G4PolyPhiFace::Extent( const G4ThreeVector axis )
582 { 607 {
583 G4double max = -kInfinity; 608 G4double max = -kInfinity;
584 609
585 G4PolyPhiFaceVertex* corner = corners; << 610 G4PolyPhiFaceVertex *corner = corners;
586 do // Loop checking, 13.08.2015, G.Cosmo 611 do // Loop checking, 13.08.2015, G.Cosmo
587 { 612 {
588 G4double here = axis.x()*corner->r*radial. 613 G4double here = axis.x()*corner->r*radial.x()
589 + axis.y()*corner->r*radial.y() 614 + axis.y()*corner->r*radial.y()
590 + axis.z()*corner->z; 615 + axis.z()*corner->z;
591 if (here > max) max = here; 616 if (here > max) max = here;
592 } while( ++corner < corners + numEdges ); 617 } while( ++corner < corners + numEdges );
593 618
594 return max; 619 return max;
595 } 620 }
596 621
>> 622
>> 623 //
597 // CalculateExtent 624 // CalculateExtent
598 // 625 //
599 // See notes in G4VCSGface 626 // See notes in G4VCSGface
600 // 627 //
601 void G4PolyPhiFace::CalculateExtent( const EAx 628 void G4PolyPhiFace::CalculateExtent( const EAxis axis,
602 const G4V << 629 const G4VoxelLimits &voxelLimit,
603 const G4A << 630 const G4AffineTransform &transform,
604 G4S << 631 G4SolidExtentList &extentList )
605 { 632 {
606 // 633 //
607 // Construct a (sometimes big) clippable pol 634 // Construct a (sometimes big) clippable polygon,
608 // 635 //
609 // Perform the necessary transformations whi 636 // Perform the necessary transformations while doing so
610 // 637 //
611 G4ClippablePolygon polygon; 638 G4ClippablePolygon polygon;
612 639
613 G4PolyPhiFaceVertex* corner = corners; << 640 G4PolyPhiFaceVertex *corner = corners;
614 do // Loop checking, 13.08.2015, G.Cosmo 641 do // Loop checking, 13.08.2015, G.Cosmo
615 { 642 {
616 G4ThreeVector point( 0, 0, corner->z ); 643 G4ThreeVector point( 0, 0, corner->z );
617 point += radial*corner->r; 644 point += radial*corner->r;
618 645
619 polygon.AddVertexInOrder( transform.Transf 646 polygon.AddVertexInOrder( transform.TransformPoint( point ) );
620 } while( ++corner < corners + numEdges ); 647 } while( ++corner < corners + numEdges );
621 648
622 // 649 //
623 // Clip away 650 // Clip away
624 // 651 //
625 if (polygon.PartialClip( voxelLimit, axis )) 652 if (polygon.PartialClip( voxelLimit, axis ))
626 { 653 {
627 // 654 //
628 // Add it to the list 655 // Add it to the list
629 // 656 //
630 polygon.SetNormal( transform.TransformAxis 657 polygon.SetNormal( transform.TransformAxis(normal) );
631 extentList.AddSurface( polygon ); 658 extentList.AddSurface( polygon );
632 } 659 }
633 } 660 }
634 661
>> 662
>> 663 //
>> 664 //-------------------------------------------------------
>> 665
>> 666
>> 667 //
635 // InsideEdgesExact 668 // InsideEdgesExact
636 // 669 //
637 // Decide if the point in r,z is inside the ed 670 // Decide if the point in r,z is inside the edges of our face,
638 // **but** do so consistently with other faces 671 // **but** do so consistently with other faces.
639 // 672 //
640 // This routine has functionality similar to I 673 // This routine has functionality similar to InsideEdges, but uses
641 // an algorithm to decide if a trajectory fall 674 // an algorithm to decide if a trajectory falls inside or outside the
642 // face that uses only the trajectory p,v valu 675 // face that uses only the trajectory p,v values and the three dimensional
643 // points representing the edges of the polygo 676 // points representing the edges of the polygon. The objective is to plug up
644 // any leaks between touching G4PolyPhiFaces ( 677 // any leaks between touching G4PolyPhiFaces (at r==0) and any other face
645 // that uses the same convention. 678 // that uses the same convention.
646 // 679 //
647 // See: "Computational Geometry in C (Second E 680 // See: "Computational Geometry in C (Second Edition)"
648 // http://cs.smith.edu/~orourke/ 681 // http://cs.smith.edu/~orourke/
649 // 682 //
650 G4bool G4PolyPhiFace::InsideEdgesExact( G4doub 683 G4bool G4PolyPhiFace::InsideEdgesExact( G4double r, G4double z,
651 G4doub 684 G4double normSign,
652 const G4Thre << 685 const G4ThreeVector &p,
653 const G4Thre << 686 const G4ThreeVector &v )
654 { 687 {
655 // 688 //
656 // Quick check of extent 689 // Quick check of extent
657 // 690 //
658 if ( (r < rMin-kCarTolerance) 691 if ( (r < rMin-kCarTolerance)
659 || (r > rMax+kCarTolerance) ) return false 692 || (r > rMax+kCarTolerance) ) return false;
660 693
661 if ( (z < zMin-kCarTolerance) 694 if ( (z < zMin-kCarTolerance)
662 || (z > zMax+kCarTolerance) ) return false 695 || (z > zMax+kCarTolerance) ) return false;
663 696
664 // 697 //
665 // Exact check: loop over all vertices 698 // Exact check: loop over all vertices
666 // 699 //
667 G4double qx = p.x() + v.x(), 700 G4double qx = p.x() + v.x(),
668 qy = p.y() + v.y(), 701 qy = p.y() + v.y(),
669 qz = p.z() + v.z(); 702 qz = p.z() + v.z();
670 703
671 G4int answer = 0; 704 G4int answer = 0;
672 G4PolyPhiFaceVertex *corn = corners, 705 G4PolyPhiFaceVertex *corn = corners,
673 *prev = corners+numEdges 706 *prev = corners+numEdges-1;
674 707
675 G4double cornZ, prevZ; 708 G4double cornZ, prevZ;
676 709
677 prevZ = ExactZOrder( z, qx, qy, qz, v, normS 710 prevZ = ExactZOrder( z, qx, qy, qz, v, normSign, prev );
678 do // Loop checking, 13.08.2015, G.Cosmo 711 do // Loop checking, 13.08.2015, G.Cosmo
679 { 712 {
680 // 713 //
681 // Get z order of this vertex, and compare 714 // Get z order of this vertex, and compare to previous vertex
682 // 715 //
683 cornZ = ExactZOrder( z, qx, qy, qz, v, nor 716 cornZ = ExactZOrder( z, qx, qy, qz, v, normSign, corn );
684 717
685 if (cornZ < 0) 718 if (cornZ < 0)
686 { 719 {
687 if (prevZ < 0) continue; 720 if (prevZ < 0) continue;
688 } 721 }
689 else if (cornZ > 0) 722 else if (cornZ > 0)
690 { 723 {
691 if (prevZ > 0) continue; 724 if (prevZ > 0) continue;
692 } 725 }
693 else 726 else
694 { 727 {
695 // 728 //
696 // By chance, we overlap exactly (within 729 // By chance, we overlap exactly (within precision) with
697 // the current vertex. Continue if the s 730 // the current vertex. Continue if the same happened previously
698 // (e.g. the previous vertex had the sam 731 // (e.g. the previous vertex had the same z value)
699 // 732 //
700 if (prevZ == 0) continue; 733 if (prevZ == 0) continue;
701 734
702 // 735 //
703 // Otherwise, to decide what to do, we n 736 // Otherwise, to decide what to do, we need to know what is
704 // coming up next. Specifically, we need 737 // coming up next. Specifically, we need to find the next vertex
705 // with a non-zero z order. 738 // with a non-zero z order.
706 // 739 //
707 // One might worry about infinite loops, 740 // One might worry about infinite loops, but the above conditional
708 // should prevent it 741 // should prevent it
709 // 742 //
710 G4PolyPhiFaceVertex *next = corn; 743 G4PolyPhiFaceVertex *next = corn;
711 G4double nextZ; 744 G4double nextZ;
712 do // Loop checking, 13.08.2015, G.Co 745 do // Loop checking, 13.08.2015, G.Cosmo
713 { 746 {
714 ++next; << 747 next++;
715 if (next == corners+numEdges) next = c 748 if (next == corners+numEdges) next = corners;
716 749
717 nextZ = ExactZOrder( z, qx, qy, qz, v, 750 nextZ = ExactZOrder( z, qx, qy, qz, v, normSign, next );
718 } while( nextZ == 0 ); 751 } while( nextZ == 0 );
719 752
720 // 753 //
721 // If we won't be changing direction, go 754 // If we won't be changing direction, go to the next vertex
722 // 755 //
723 if (nextZ*prevZ < 0) continue; 756 if (nextZ*prevZ < 0) continue;
724 } 757 }
725 758
726 759
727 // 760 //
728 // We overlap in z with the side of the fa 761 // We overlap in z with the side of the face that stretches from
729 // vertex "prev" to "corn". On which side 762 // vertex "prev" to "corn". On which side (left or right) do
730 // we lay with respect to this segment? 763 // we lay with respect to this segment?
731 // 764 //
732 G4ThreeVector qa( qx - prev->x, qy - prev- 765 G4ThreeVector qa( qx - prev->x, qy - prev->y, qz - prev->z ),
733 qb( qx - corn->x, qy - corn- 766 qb( qx - corn->x, qy - corn->y, qz - corn->z );
734 767
735 G4double aboveOrBelow = normSign*qa.cross( 768 G4double aboveOrBelow = normSign*qa.cross(qb).dot(v);
736 769
737 if (aboveOrBelow > 0) 770 if (aboveOrBelow > 0)
738 ++answer; << 771 answer++;
739 else if (aboveOrBelow < 0) 772 else if (aboveOrBelow < 0)
740 --answer; << 773 answer--;
741 else 774 else
742 { 775 {
743 // 776 //
744 // A precisely zero answer here means we 777 // A precisely zero answer here means we exactly
745 // intersect (within roundoff) the edge 778 // intersect (within roundoff) the edge of the face.
746 // Return true in this case. 779 // Return true in this case.
747 // 780 //
748 return true; 781 return true;
749 } 782 }
750 } while( prevZ = cornZ, prev=corn, ++corn < 783 } while( prevZ = cornZ, prev=corn, ++corn < corners+numEdges );
751 784
752 return answer!=0; 785 return answer!=0;
753 } 786 }
754 << 787
755 // InsideEdges (don't care about distance) << 788
>> 789 //
>> 790 // InsideEdges (don't care aboud distance)
756 // 791 //
757 // Decide if the point in r,z is inside the ed 792 // Decide if the point in r,z is inside the edges of our face
758 // 793 //
759 // This routine can be made a zillion times qu 794 // This routine can be made a zillion times quicker by implementing
760 // better code, for example: 795 // better code, for example:
761 // 796 //
762 // int pnpoly(int npol, float *xp, float *y 797 // int pnpoly(int npol, float *xp, float *yp, float x, float y)
763 // { 798 // {
764 // int i, j, c = 0; 799 // int i, j, c = 0;
765 // for (i = 0, j = npol-1; i < npol; j = 800 // for (i = 0, j = npol-1; i < npol; j = i++) {
766 // if ((((yp[i]<=y) && (y<yp[j])) || 801 // if ((((yp[i]<=y) && (y<yp[j])) ||
767 // ((yp[j]<=y) && (y<yp[i]))) && 802 // ((yp[j]<=y) && (y<yp[i]))) &&
768 // (x < (xp[j] - xp[i]) * (y - yp[i 803 // (x < (xp[j] - xp[i]) * (y - yp[i]) / (yp[j] - yp[i]) + xp[i]))
769 // 804 //
770 // c = !c; 805 // c = !c;
771 // } 806 // }
772 // return c; 807 // return c;
773 // } 808 // }
774 // 809 //
775 // See "Point in Polyon Strategies", Eric Hain 810 // See "Point in Polyon Strategies", Eric Haines [Graphic Gems IV] pp. 24-46
776 // 811 //
777 // The algorithm below is rather unique, but i << 812 // My algorithm below is rather unique, but is based on code needed to
778 // calculate the distance to the shape. Left h << 813 // calculate the distance to the shape. I left it in here because ...
>> 814 // well ... to test it better.
779 // 815 //
780 G4bool G4PolyPhiFace::InsideEdges( G4double r, 816 G4bool G4PolyPhiFace::InsideEdges( G4double r, G4double z )
781 { 817 {
782 // 818 //
783 // Quick check of extent 819 // Quick check of extent
784 // 820 //
785 if ( r < rMin || r > rMax ) return false; 821 if ( r < rMin || r > rMax ) return false;
786 if ( z < zMin || z > zMax ) return false; 822 if ( z < zMin || z > zMax ) return false;
787 823
788 // 824 //
789 // More thorough check 825 // More thorough check
790 // 826 //
791 G4double notUsed; 827 G4double notUsed;
792 828
793 return InsideEdges( r, z, ¬Used, nullptr << 829 return InsideEdges( r, z, ¬Used, 0 );
794 } 830 }
795 831
>> 832
>> 833 //
796 // InsideEdges (care about distance) 834 // InsideEdges (care about distance)
797 // 835 //
798 // Decide if the point in r,z is inside the ed 836 // Decide if the point in r,z is inside the edges of our face
799 // 837 //
800 G4bool G4PolyPhiFace::InsideEdges( G4double r, 838 G4bool G4PolyPhiFace::InsideEdges( G4double r, G4double z,
801 G4double* b << 839 G4double *bestDist2,
802 G4PolyPhiFa << 840 G4PolyPhiFaceVertex **base3Dnorm,
803 G4ThreeVect << 841 G4ThreeVector **head3Dnorm )
804 { 842 {
805 G4double bestDistance2 = kInfinity; 843 G4double bestDistance2 = kInfinity;
806 G4bool answer = false; << 844 G4bool answer = 0;
807 845
808 G4PolyPhiFaceEdge* edge = edges; << 846 G4PolyPhiFaceEdge *edge = edges;
809 do // Loop checking, 13.08.2015, G.Cosmo 847 do // Loop checking, 13.08.2015, G.Cosmo
810 { 848 {
811 G4PolyPhiFaceVertex* testMe = nullptr; << 849 G4PolyPhiFaceVertex *testMe;
812 G4PolyPhiFaceVertex* v0_vertex = edge->v0; <<
813 G4PolyPhiFaceVertex* v1_vertex = edge->v1; <<
814 // 850 //
815 // Get distance perpendicular to the edge 851 // Get distance perpendicular to the edge
816 // 852 //
817 G4double dr = (r-v0_vertex->r), dz = (z-v0 << 853 G4double dr = (r-edge->v0->r), dz = (z-edge->v0->z);
818 854
819 G4double distOut = dr*edge->tz - dz*edge-> 855 G4double distOut = dr*edge->tz - dz*edge->tr;
820 G4double distance2 = distOut*distOut; 856 G4double distance2 = distOut*distOut;
821 if (distance2 > bestDistance2) continue; 857 if (distance2 > bestDistance2) continue; // No hope!
822 858
823 // 859 //
824 // Check to see if normal intersects edge 860 // Check to see if normal intersects edge within the edge's boundary
825 // 861 //
826 G4double q = dr*edge->tr + dz*edge->tz; 862 G4double q = dr*edge->tr + dz*edge->tz;
827 863
828 // 864 //
829 // If it doesn't, penalize distance2 appro 865 // If it doesn't, penalize distance2 appropriately
830 // 866 //
831 if (q < 0) 867 if (q < 0)
832 { 868 {
833 distance2 += q*q; 869 distance2 += q*q;
834 testMe = v0_vertex; << 870 testMe = edge->v0;
835 } 871 }
836 else if (q > edge->length) 872 else if (q > edge->length)
837 { 873 {
838 G4double s2 = q-edge->length; 874 G4double s2 = q-edge->length;
839 distance2 += s2*s2; 875 distance2 += s2*s2;
840 testMe = v1_vertex; << 876 testMe = edge->v1;
>> 877 }
>> 878 else
>> 879 {
>> 880 testMe = 0;
841 } 881 }
842 882
843 // 883 //
844 // Closest edge so far? 884 // Closest edge so far?
845 // 885 //
846 if (distance2 < bestDistance2) 886 if (distance2 < bestDistance2)
847 { 887 {
848 bestDistance2 = distance2; 888 bestDistance2 = distance2;
849 if (testMe != nullptr) << 889 if (testMe)
850 { 890 {
851 G4double distNorm = dr*testMe->rNorm + 891 G4double distNorm = dr*testMe->rNorm + dz*testMe->zNorm;
852 answer = (distNorm <= 0); 892 answer = (distNorm <= 0);
853 if (base3Dnorm != nullptr) << 893 if (base3Dnorm)
854 { 894 {
855 *base3Dnorm = testMe; 895 *base3Dnorm = testMe;
856 *head3Dnorm = &testMe->norm3D; 896 *head3Dnorm = &testMe->norm3D;
857 } 897 }
858 } 898 }
859 else 899 else
860 { 900 {
861 answer = (distOut <= 0); 901 answer = (distOut <= 0);
862 if (base3Dnorm != nullptr) << 902 if (base3Dnorm)
863 { 903 {
864 *base3Dnorm = v0_vertex; << 904 *base3Dnorm = edge->v0;
865 *head3Dnorm = &edge->norm3D; 905 *head3Dnorm = &edge->norm3D;
866 } 906 }
867 } 907 }
868 } 908 }
869 } while( ++edge < edges + numEdges ); 909 } while( ++edge < edges + numEdges );
870 910
871 *bestDist2 = bestDistance2; 911 *bestDist2 = bestDistance2;
872 return answer; 912 return answer;
873 } 913 }
874 914
>> 915 //
875 // Calculation of Surface Area of a Triangle 916 // Calculation of Surface Area of a Triangle
876 // In the same time Random Point in Triangle i 917 // In the same time Random Point in Triangle is given
877 // 918 //
878 G4double G4PolyPhiFace::SurfaceTriangle( const << 919 G4double G4PolyPhiFace::SurfaceTriangle( G4ThreeVector p1,
879 const << 920 G4ThreeVector p2,
880 const << 921 G4ThreeVector p3,
881 G4Thr << 922 G4ThreeVector *p4 )
882 { 923 {
883 G4ThreeVector v, w; 924 G4ThreeVector v, w;
884 925
885 v = p3 - p1; 926 v = p3 - p1;
886 w = p1 - p2; 927 w = p1 - p2;
887 G4double lambda1 = G4UniformRand(); 928 G4double lambda1 = G4UniformRand();
888 G4double lambda2 = lambda1*G4UniformRand(); 929 G4double lambda2 = lambda1*G4UniformRand();
889 930
890 *p4=p2 + lambda1*w + lambda2*v; 931 *p4=p2 + lambda1*w + lambda2*v;
891 return 0.5*(v.cross(w)).mag(); 932 return 0.5*(v.cross(w)).mag();
892 } 933 }
893 934
>> 935 //
894 // Compute surface area 936 // Compute surface area
895 // 937 //
896 G4double G4PolyPhiFace::SurfaceArea() 938 G4double G4PolyPhiFace::SurfaceArea()
897 { 939 {
898 if ( fSurfaceArea==0. ) { Triangulate(); } 940 if ( fSurfaceArea==0. ) { Triangulate(); }
899 return fSurfaceArea; 941 return fSurfaceArea;
900 } 942 }
901 943
>> 944 //
902 // Return random point on face 945 // Return random point on face
903 // 946 //
904 G4ThreeVector G4PolyPhiFace::GetPointOnFace() 947 G4ThreeVector G4PolyPhiFace::GetPointOnFace()
905 { 948 {
906 Triangulate(); 949 Triangulate();
907 return surface_point; 950 return surface_point;
908 } 951 }
909 952
910 // 953 //
911 // Auxiliary Functions used for Finding the Po 954 // Auxiliary Functions used for Finding the PointOnFace using Triangulation
912 // 955 //
913 956
>> 957 //
914 // Calculation of 2*Area of Triangle with Sign 958 // Calculation of 2*Area of Triangle with Sign
915 // 959 //
916 G4double G4PolyPhiFace::Area2( const G4TwoVect << 960 G4double G4PolyPhiFace::Area2( G4TwoVector a,
917 const G4TwoVect << 961 G4TwoVector b,
918 const G4TwoVect << 962 G4TwoVector c )
919 { 963 {
920 return ((b.x()-a.x())*(c.y()-a.y())- 964 return ((b.x()-a.x())*(c.y()-a.y())-
921 (c.x()-a.x())*(b.y()-a.y())); 965 (c.x()-a.x())*(b.y()-a.y()));
922 } 966 }
923 967
>> 968 //
924 // Boolean function for sign of Surface 969 // Boolean function for sign of Surface
925 // 970 //
926 G4bool G4PolyPhiFace::Left( const G4TwoVector& << 971 G4bool G4PolyPhiFace::Left( G4TwoVector a,
927 const G4TwoVector& << 972 G4TwoVector b,
928 const G4TwoVector& << 973 G4TwoVector c )
929 { 974 {
930 return Area2(a,b,c)>0; 975 return Area2(a,b,c)>0;
931 } 976 }
932 977
>> 978 //
933 // Boolean function for sign of Surface 979 // Boolean function for sign of Surface
934 // 980 //
935 G4bool G4PolyPhiFace::LeftOn( const G4TwoVecto << 981 G4bool G4PolyPhiFace::LeftOn( G4TwoVector a,
936 const G4TwoVecto << 982 G4TwoVector b,
937 const G4TwoVecto << 983 G4TwoVector c )
938 { 984 {
939 return Area2(a,b,c)>=0; 985 return Area2(a,b,c)>=0;
940 } 986 }
941 987
>> 988 //
942 // Boolean function for sign of Surface 989 // Boolean function for sign of Surface
943 // 990 //
944 G4bool G4PolyPhiFace::Collinear( const G4TwoVe << 991 G4bool G4PolyPhiFace::Collinear( G4TwoVector a,
945 const G4TwoVe << 992 G4TwoVector b,
946 const G4TwoVe << 993 G4TwoVector c )
947 { 994 {
948 return Area2(a,b,c)==0; 995 return Area2(a,b,c)==0;
949 } 996 }
950 997
>> 998 //
951 // Boolean function for finding "Proper" Inter 999 // Boolean function for finding "Proper" Intersection
952 // That means Intersection of two lines segmen 1000 // That means Intersection of two lines segments (a,b) and (c,d)
953 // 1001 //
954 G4bool G4PolyPhiFace::IntersectProp( const G4T << 1002 G4bool G4PolyPhiFace::IntersectProp( G4TwoVector a,
955 const G4T << 1003 G4TwoVector b,
956 const G4T << 1004 G4TwoVector c, G4TwoVector d )
957 { 1005 {
958 if( Collinear(a,b,c) || Collinear(a,b,d)|| 1006 if( Collinear(a,b,c) || Collinear(a,b,d)||
959 Collinear(c,d,a) || Collinear(c,d,b) ) 1007 Collinear(c,d,a) || Collinear(c,d,b) ) { return false; }
960 1008
961 G4bool Positive; 1009 G4bool Positive;
962 Positive = !(Left(a,b,c))^!(Left(a,b,d)); 1010 Positive = !(Left(a,b,c))^!(Left(a,b,d));
963 return Positive && ((!Left(c,d,a)^!Left(c,d, << 1011 return Positive && (!Left(c,d,a)^!Left(c,d,b));
964 } 1012 }
965 1013
>> 1014 //
966 // Boolean function for determining if Point c 1015 // Boolean function for determining if Point c is between a and b
967 // For the tree points(a,b,c) on the same line 1016 // For the tree points(a,b,c) on the same line
968 // 1017 //
969 G4bool G4PolyPhiFace::Between( const G4TwoVect << 1018 G4bool G4PolyPhiFace::Between( G4TwoVector a, G4TwoVector b, G4TwoVector c )
970 { 1019 {
971 if( !Collinear(a,b,c) ) { return false; } 1020 if( !Collinear(a,b,c) ) { return false; }
972 1021
973 if(a.x()!=b.x()) 1022 if(a.x()!=b.x())
974 { 1023 {
975 return ((a.x()<=c.x())&&(c.x()<=b.x()))|| 1024 return ((a.x()<=c.x())&&(c.x()<=b.x()))||
976 ((a.x()>=c.x())&&(c.x()>=b.x())); 1025 ((a.x()>=c.x())&&(c.x()>=b.x()));
977 } 1026 }
978 else 1027 else
979 { 1028 {
980 return ((a.y()<=c.y())&&(c.y()<=b.y()))|| 1029 return ((a.y()<=c.y())&&(c.y()<=b.y()))||
981 ((a.y()>=c.y())&&(c.y()>=b.y())); 1030 ((a.y()>=c.y())&&(c.y()>=b.y()));
982 } 1031 }
983 } 1032 }
984 1033
>> 1034 //
985 // Boolean function for finding Intersection " 1035 // Boolean function for finding Intersection "Proper" or not
986 // Between two line segments (a,b) and (c,d) 1036 // Between two line segments (a,b) and (c,d)
987 // 1037 //
988 G4bool G4PolyPhiFace::Intersect( const G4TwoVe << 1038 G4bool G4PolyPhiFace::Intersect( G4TwoVector a,
989 const G4TwoVe << 1039 G4TwoVector b,
990 const G4TwoVe << 1040 G4TwoVector c, G4TwoVector d )
991 { 1041 {
992 if( IntersectProp(a,b,c,d) ) 1042 if( IntersectProp(a,b,c,d) )
993 { return true; } 1043 { return true; }
994 else if( Between(a,b,c)|| 1044 else if( Between(a,b,c)||
995 Between(a,b,d)|| 1045 Between(a,b,d)||
996 Between(c,d,a)|| 1046 Between(c,d,a)||
997 Between(c,d,b) ) 1047 Between(c,d,b) )
998 { return true; } 1048 { return true; }
999 else 1049 else
1000 { return false; } 1050 { return false; }
1001 } 1051 }
1002 1052
>> 1053 //
1003 // Boolean Diagonalie help to determine 1054 // Boolean Diagonalie help to determine
1004 // if diagonal s of segment (a,b) is convex o 1055 // if diagonal s of segment (a,b) is convex or reflex
1005 // 1056 //
1006 G4bool G4PolyPhiFace::Diagonalie( G4PolyPhiFa << 1057 G4bool G4PolyPhiFace::Diagonalie( G4PolyPhiFaceVertex *a,
1007 G4PolyPhiFa << 1058 G4PolyPhiFaceVertex *b )
1008 { 1059 {
1009 G4PolyPhiFaceVertex *corner = triangles; 1060 G4PolyPhiFaceVertex *corner = triangles;
1010 G4PolyPhiFaceVertex *corner_next=triangle 1061 G4PolyPhiFaceVertex *corner_next=triangles;
1011 1062
1012 // For each Edge (corner,corner_next) 1063 // For each Edge (corner,corner_next)
1013 do // Loop checking, 13.08.2015, G.Cosmo 1064 do // Loop checking, 13.08.2015, G.Cosmo
1014 { 1065 {
1015 corner_next=corner->next; 1066 corner_next=corner->next;
1016 1067
1017 // Skip edges incident to a of b 1068 // Skip edges incident to a of b
1018 // 1069 //
1019 if( (corner!=a)&&(corner_next!=a) 1070 if( (corner!=a)&&(corner_next!=a)
1020 &&(corner!=b)&&(corner_next!=b) ) 1071 &&(corner!=b)&&(corner_next!=b) )
1021 { 1072 {
1022 G4TwoVector rz1,rz2,rz3,rz4; 1073 G4TwoVector rz1,rz2,rz3,rz4;
1023 rz1 = G4TwoVector(a->r,a->z); 1074 rz1 = G4TwoVector(a->r,a->z);
1024 rz2 = G4TwoVector(b->r,b->z); 1075 rz2 = G4TwoVector(b->r,b->z);
1025 rz3 = G4TwoVector(corner->r,corner->z) 1076 rz3 = G4TwoVector(corner->r,corner->z);
1026 rz4 = G4TwoVector(corner_next->r,corne 1077 rz4 = G4TwoVector(corner_next->r,corner_next->z);
1027 if( Intersect(rz1,rz2,rz3,rz4) ) { ret 1078 if( Intersect(rz1,rz2,rz3,rz4) ) { return false; }
1028 } 1079 }
1029 corner=corner->next; 1080 corner=corner->next;
1030 1081
1031 } while( corner != triangles ); 1082 } while( corner != triangles );
1032 1083
1033 return true; 1084 return true;
1034 } 1085 }
1035 1086
>> 1087 //
1036 // Boolean function that determine if b is In 1088 // Boolean function that determine if b is Inside Cone (a0,a,a1)
1037 // being a the center of the Cone 1089 // being a the center of the Cone
1038 // 1090 //
1039 G4bool G4PolyPhiFace::InCone( G4PolyPhiFaceVe << 1091 G4bool G4PolyPhiFace::InCone( G4PolyPhiFaceVertex *a, G4PolyPhiFaceVertex *b )
1040 { 1092 {
1041 // a0,a and a1 are consecutive vertices 1093 // a0,a and a1 are consecutive vertices
1042 // 1094 //
1043 G4PolyPhiFaceVertex *a0,*a1; 1095 G4PolyPhiFaceVertex *a0,*a1;
1044 a1=a->next; 1096 a1=a->next;
1045 a0=a->prev; 1097 a0=a->prev;
1046 1098
1047 G4TwoVector arz,arz0,arz1,brz; 1099 G4TwoVector arz,arz0,arz1,brz;
1048 arz=G4TwoVector(a->r,a->z);arz0=G4TwoVector 1100 arz=G4TwoVector(a->r,a->z);arz0=G4TwoVector(a0->r,a0->z);
1049 arz1=G4TwoVector(a1->r,a1->z);brz=G4TwoVect 1101 arz1=G4TwoVector(a1->r,a1->z);brz=G4TwoVector(b->r,b->z);
1050 1102
1051 1103
1052 if(LeftOn(arz,arz1,arz0)) // If a is conve 1104 if(LeftOn(arz,arz1,arz0)) // If a is convex vertex
1053 { 1105 {
1054 return Left(arz,brz,arz0)&&Left(brz,arz,a 1106 return Left(arz,brz,arz0)&&Left(brz,arz,arz1);
1055 } 1107 }
1056 else // Else a is ref 1108 else // Else a is reflex
1057 { 1109 {
1058 return !( LeftOn(arz,brz,arz1)&&LeftOn(br 1110 return !( LeftOn(arz,brz,arz1)&&LeftOn(brz,arz,arz0));
1059 } 1111 }
1060 } 1112 }
1061 1113
>> 1114 //
1062 // Boolean function finding if Diagonal is po 1115 // Boolean function finding if Diagonal is possible
1063 // inside Polycone or PolyHedra 1116 // inside Polycone or PolyHedra
1064 // 1117 //
1065 G4bool G4PolyPhiFace::Diagonal( G4PolyPhiFace << 1118 G4bool G4PolyPhiFace::Diagonal( G4PolyPhiFaceVertex *a, G4PolyPhiFaceVertex *b )
1066 { 1119 {
1067 return InCone(a,b) && InCone(b,a) && Diagon 1120 return InCone(a,b) && InCone(b,a) && Diagonalie(a,b);
1068 } 1121 }
1069 1122
>> 1123 //
1070 // Initialisation for Triangulisation by ear 1124 // Initialisation for Triangulisation by ear tips
1071 // For details see "Computational Geometry in 1125 // For details see "Computational Geometry in C" by Joseph O'Rourke
1072 // 1126 //
1073 void G4PolyPhiFace::EarInit() 1127 void G4PolyPhiFace::EarInit()
1074 { 1128 {
1075 G4PolyPhiFaceVertex* corner = triangles; << 1129 G4PolyPhiFaceVertex *corner = triangles;
1076 G4PolyPhiFaceVertex* c_prev, *c_next; << 1130 G4PolyPhiFaceVertex *c_prev,*c_next;
1077 1131
1078 do // Loop checking, 13.08.2015, G.Cosmo 1132 do // Loop checking, 13.08.2015, G.Cosmo
1079 { 1133 {
1080 // We need to determine three consecutiv 1134 // We need to determine three consecutive vertices
1081 // 1135 //
1082 c_next=corner->next; 1136 c_next=corner->next;
1083 c_prev=corner->prev; 1137 c_prev=corner->prev;
1084 1138
1085 // Calculation of ears 1139 // Calculation of ears
1086 // 1140 //
1087 corner->ear=Diagonal(c_prev,c_next); 1141 corner->ear=Diagonal(c_prev,c_next);
1088 corner=corner->next; 1142 corner=corner->next;
1089 1143
1090 } while( corner!=triangles ); 1144 } while( corner!=triangles );
1091 } 1145 }
1092 1146
>> 1147 //
1093 // Triangulisation by ear tips for Polycone o 1148 // Triangulisation by ear tips for Polycone or Polyhedra
1094 // For details see "Computational Geometry in 1149 // For details see "Computational Geometry in C" by Joseph O'Rourke
1095 // 1150 //
1096 void G4PolyPhiFace::Triangulate() 1151 void G4PolyPhiFace::Triangulate()
1097 { 1152 {
1098 // The copy of Polycone is made and this co 1153 // The copy of Polycone is made and this copy is reordered in order to
1099 // have a list of triangles. This list is u 1154 // have a list of triangles. This list is used for GetPointOnFace().
1100 1155
1101 const std::size_t maxEdges = (numEdges > 0) << 1156 G4PolyPhiFaceVertex *tri_help = new G4PolyPhiFaceVertex[numEdges];
1102 auto tri_help = new G4PolyPhiFaceVertex[max <<
1103 triangles = tri_help; 1157 triangles = tri_help;
1104 G4PolyPhiFaceVertex* triang = triangles; << 1158 G4PolyPhiFaceVertex *triang = triangles;
1105 1159
1106 std::vector<G4double> areas; 1160 std::vector<G4double> areas;
1107 std::vector<G4ThreeVector> points; 1161 std::vector<G4ThreeVector> points;
1108 G4double area=0.; 1162 G4double area=0.;
1109 G4PolyPhiFaceVertex *v0,*v1,*v2,*v3,*v4; 1163 G4PolyPhiFaceVertex *v0,*v1,*v2,*v3,*v4;
1110 v2=triangles; 1164 v2=triangles;
1111 1165
1112 // Make copy for prev/next for triang=corne 1166 // Make copy for prev/next for triang=corners
1113 // 1167 //
1114 G4PolyPhiFaceVertex* helper = corners; << 1168 G4PolyPhiFaceVertex *helper = corners;
1115 G4PolyPhiFaceVertex* helper2 = corners; << 1169 G4PolyPhiFaceVertex *helper2 = corners;
1116 do // Loop checking, 13.08.2015, G.Cosmo 1170 do // Loop checking, 13.08.2015, G.Cosmo
1117 { 1171 {
1118 triang->r = helper->r; 1172 triang->r = helper->r;
1119 triang->z = helper->z; 1173 triang->z = helper->z;
1120 triang->x = helper->x; 1174 triang->x = helper->x;
1121 triang->y = helper->y; << 1175 triang->y= helper->y;
1122 1176
1123 // add pointer on prev corner 1177 // add pointer on prev corner
1124 // 1178 //
1125 if( helper==corners ) 1179 if( helper==corners )
1126 { triang->prev=triangles+maxEdges-1; } << 1180 { triang->prev=triangles+numEdges-1; }
1127 else 1181 else
1128 { triang->prev=helper2; } 1182 { triang->prev=helper2; }
1129 1183
1130 // add pointer on next corner 1184 // add pointer on next corner
1131 // 1185 //
1132 if( helper<corners+maxEdges-1 ) << 1186 if( helper<corners+numEdges-1 )
1133 { triang->next=triang+1; } 1187 { triang->next=triang+1; }
1134 else 1188 else
1135 { triang->next=triangles; } 1189 { triang->next=triangles; }
1136 helper2=triang; 1190 helper2=triang;
1137 helper=helper->next; 1191 helper=helper->next;
1138 triang=triang->next; 1192 triang=triang->next;
1139 1193
1140 } while( helper!=corners ); 1194 } while( helper!=corners );
1141 1195
1142 EarInit(); 1196 EarInit();
1143 1197
1144 G4int n=numEdges; 1198 G4int n=numEdges;
1145 G4int i=0; 1199 G4int i=0;
1146 G4ThreeVector p1,p2,p3,p4; 1200 G4ThreeVector p1,p2,p3,p4;
1147 const G4int max_n_loops=numEdges*10000; // 1201 const G4int max_n_loops=numEdges*10000; // protection against infinite loop
1148 1202
1149 // Each step of outer loop removes one ear 1203 // Each step of outer loop removes one ear
1150 // 1204 //
1151 while(n>3) // Loop checking, 13.08.201 1205 while(n>3) // Loop checking, 13.08.2015, G.Cosmo
1152 { // Inner loop searches for 1206 { // Inner loop searches for one ear
1153 v2=triangles; 1207 v2=triangles;
1154 do // Loop checking, 13.08.2015, G.Cos 1208 do // Loop checking, 13.08.2015, G.Cosmo
1155 { 1209 {
1156 if(v2->ear) // Ear found. Fill variabl 1210 if(v2->ear) // Ear found. Fill variables
1157 { 1211 {
1158 // (v1,v3) is diagonal 1212 // (v1,v3) is diagonal
1159 // 1213 //
1160 v3=v2->next; v4=v3->next; 1214 v3=v2->next; v4=v3->next;
1161 v1=v2->prev; v0=v1->prev; 1215 v1=v2->prev; v0=v1->prev;
1162 1216
1163 // Calculate areas and points 1217 // Calculate areas and points
1164 1218
1165 p1=G4ThreeVector((v2)->x,(v2)->y,(v2) 1219 p1=G4ThreeVector((v2)->x,(v2)->y,(v2)->z);
1166 p2=G4ThreeVector((v1)->x,(v1)->y,(v1) 1220 p2=G4ThreeVector((v1)->x,(v1)->y,(v1)->z);
1167 p3=G4ThreeVector((v3)->x,(v3)->y,(v3) 1221 p3=G4ThreeVector((v3)->x,(v3)->y,(v3)->z);
1168 1222
1169 G4double result1 = SurfaceTriangle(p1 1223 G4double result1 = SurfaceTriangle(p1,p2,p3,&p4 );
1170 points.push_back(p4); 1224 points.push_back(p4);
1171 areas.push_back(result1); 1225 areas.push_back(result1);
1172 area=area+result1; 1226 area=area+result1;
1173 1227
1174 // Update earity of diagonal endpoint 1228 // Update earity of diagonal endpoints
1175 // 1229 //
1176 v1->ear=Diagonal(v0,v3); 1230 v1->ear=Diagonal(v0,v3);
1177 v3->ear=Diagonal(v1,v4); 1231 v3->ear=Diagonal(v1,v4);
1178 1232
1179 // Cut off the ear v2 1233 // Cut off the ear v2
1180 // Has to be done for a copy and not 1234 // Has to be done for a copy and not for real PolyPhiFace
1181 // 1235 //
1182 v1->next=v3; 1236 v1->next=v3;
1183 v3->prev=v1; 1237 v3->prev=v1;
1184 triangles=v3; // In case the head was 1238 triangles=v3; // In case the head was v2
1185 --n; << 1239 n--;
1186 1240
1187 break; // out of inner loop 1241 break; // out of inner loop
1188 } // end if ear found 1242 } // end if ear found
1189 1243
1190 v2=v2->next; 1244 v2=v2->next;
1191 1245
1192 } while( v2!=triangles ); 1246 } while( v2!=triangles );
1193 1247
1194 ++i; << 1248 i++;
1195 if(i>=max_n_loops) 1249 if(i>=max_n_loops)
1196 { 1250 {
1197 G4Exception( "G4PolyPhiFace::Triangulat 1251 G4Exception( "G4PolyPhiFace::Triangulation()",
1198 "GeomSolids0003", FatalExc 1252 "GeomSolids0003", FatalException,
1199 "Maximum number of steps i 1253 "Maximum number of steps is reached for triangulation!" );
1200 } 1254 }
1201 } // end outer while loop 1255 } // end outer while loop
1202 1256
1203 if(v2->next != nullptr) << 1257 if(v2->next)
1204 { 1258 {
1205 // add last triangle 1259 // add last triangle
1206 // 1260 //
1207 v2=v2->next; 1261 v2=v2->next;
1208 p1=G4ThreeVector((v2)->x,(v2)->y,(v2)->z 1262 p1=G4ThreeVector((v2)->x,(v2)->y,(v2)->z);
1209 p2=G4ThreeVector((v2->next)->x,(v2->next 1263 p2=G4ThreeVector((v2->next)->x,(v2->next)->y,(v2->next)->z);
1210 p3=G4ThreeVector((v2->prev)->x,(v2->prev 1264 p3=G4ThreeVector((v2->prev)->x,(v2->prev)->y,(v2->prev)->z);
1211 G4double result1 = SurfaceTriangle(p1,p2 1265 G4double result1 = SurfaceTriangle(p1,p2,p3,&p4 );
1212 points.push_back(p4); 1266 points.push_back(p4);
1213 areas.push_back(result1); 1267 areas.push_back(result1);
1214 area=area+result1; 1268 area=area+result1;
1215 } 1269 }
1216 1270
1217 // Surface Area is stored 1271 // Surface Area is stored
1218 // 1272 //
1219 fSurfaceArea = area; 1273 fSurfaceArea = area;
1220 1274
1221 // Second Step: choose randomly one surface 1275 // Second Step: choose randomly one surface
1222 // 1276 //
1223 G4double chose = area*G4UniformRand(); 1277 G4double chose = area*G4UniformRand();
1224 1278
1225 // Third Step: Get a point on choosen surfa 1279 // Third Step: Get a point on choosen surface
1226 // 1280 //
1227 G4double Achose1, Achose2; 1281 G4double Achose1, Achose2;
1228 Achose1=0; Achose2=0.; 1282 Achose1=0; Achose2=0.;
1229 i=0; 1283 i=0;
1230 do // Loop checking, 13.08.2015, G.Cosm 1284 do // Loop checking, 13.08.2015, G.Cosmo
1231 { 1285 {
1232 Achose2+=areas[i]; 1286 Achose2+=areas[i];
1233 if(chose>=Achose1 && chose<Achose2) 1287 if(chose>=Achose1 && chose<Achose2)
1234 { 1288 {
1235 G4ThreeVector point; 1289 G4ThreeVector point;
1236 point=points[i]; << 1290 point=points[i] ;
1237 surface_point=point; << 1291 surface_point=point;
1238 break; << 1292 break;
1239 } 1293 }
1240 ++i; Achose1=Achose2; << 1294 i++; Achose1=Achose2;
1241 } while( i<numEdges-2 ); 1295 } while( i<numEdges-2 );
1242 1296
1243 delete [] tri_help; 1297 delete [] tri_help;
1244 tri_help = nullptr; << 1298 tri_help = 0;
1245 } 1299 }
1246 1300