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