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