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