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