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