| 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 G4PolyconeSide, the face <<
27 // one conical side of a polycone <<
28 // 26 //
29 // Author: David C. Williams (davidw@scipp.ucs << 27 // $Id: G4PolyconeSide.cc 100427 2016-10-21 12:58:28Z gcosmo $
>> 28 //
>> 29 //
>> 30 // --------------------------------------------------------------------
>> 31 // GEANT 4 class source file
>> 32 //
>> 33 //
>> 34 // G4PolyconeSide.cc
>> 35 //
>> 36 // Implementation of the face representing one conical side of a polycone
>> 37 //
30 // ------------------------------------------- 38 // --------------------------------------------------------------------
31 39
32 #include "G4PolyconeSide.hh" 40 #include "G4PolyconeSide.hh"
33 #include "meshdefs.hh" 41 #include "meshdefs.hh"
34 #include "G4PhysicalConstants.hh" 42 #include "G4PhysicalConstants.hh"
35 #include "G4IntersectingCone.hh" 43 #include "G4IntersectingCone.hh"
36 #include "G4ClippablePolygon.hh" 44 #include "G4ClippablePolygon.hh"
37 #include "G4AffineTransform.hh" 45 #include "G4AffineTransform.hh"
38 #include "G4SolidExtentList.hh" 46 #include "G4SolidExtentList.hh"
39 #include "G4GeometryTolerance.hh" 47 #include "G4GeometryTolerance.hh"
40 48
41 #include "Randomize.hh" 49 #include "Randomize.hh"
42 50
43 // This new field helps to use the class G4PlS 51 // This new field helps to use the class G4PlSideManager.
44 // 52 //
45 G4PlSideManager G4PolyconeSide::subInstanceMan 53 G4PlSideManager G4PolyconeSide::subInstanceManager;
46 54
47 // This macro changes the references to fields 55 // This macro changes the references to fields that are now encapsulated
48 // in the class G4PlSideData. 56 // in the class G4PlSideData.
49 // 57 //
50 #define G4MT_pcphix ((subInstanceManager.offse << 58 #define G4MT_pcphi ((subInstanceManager.offset[instanceID]).fPhi)
51 #define G4MT_pcphiy ((subInstanceManager.offse <<
52 #define G4MT_pcphiz ((subInstanceManager.offse <<
53 #define G4MT_pcphik ((subInstanceManager.offse <<
54 59
55 // Returns the private data instance manager. 60 // Returns the private data instance manager.
56 // 61 //
57 const G4PlSideManager& G4PolyconeSide::GetSubI 62 const G4PlSideManager& G4PolyconeSide::GetSubInstanceManager()
58 { 63 {
59 return subInstanceManager; 64 return subInstanceManager;
60 } 65 }
61 66
>> 67 //
62 // Constructor 68 // Constructor
63 // 69 //
64 // Values for r1,z1 and r2,z2 should be specif 70 // Values for r1,z1 and r2,z2 should be specified in clockwise
65 // order in (r,z). 71 // order in (r,z).
66 // 72 //
67 G4PolyconeSide::G4PolyconeSide( const G4Polyco << 73 G4PolyconeSide::G4PolyconeSide( const G4PolyconeSideRZ *prevRZ,
68 const G4Polyco << 74 const G4PolyconeSideRZ *tail,
69 const G4Polyco << 75 const G4PolyconeSideRZ *head,
70 const G4Polyco << 76 const G4PolyconeSideRZ *nextRZ,
71 G4double 77 G4double thePhiStart,
72 G4double 78 G4double theDeltaPhi,
73 G4bool t 79 G4bool thePhiIsOpen,
74 G4bool i 80 G4bool isAllBehind )
>> 81 : ncorners(0), corners(0)
75 { 82 {
>> 83
76 instanceID = subInstanceManager.CreateSubIns 84 instanceID = subInstanceManager.CreateSubInstance();
77 85
78 kCarTolerance = G4GeometryTolerance::GetInst 86 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
79 G4MT_pcphix = 0.0; G4MT_pcphiy = 0.0; G4MT_p << 87 fSurfaceArea = 0.0;
>> 88 G4MT_pcphi.first = G4ThreeVector(0,0,0);
>> 89 G4MT_pcphi.second= 0.0;
80 90
81 // 91 //
82 // Record values 92 // Record values
83 // 93 //
84 r[0] = tail->r; z[0] = tail->z; 94 r[0] = tail->r; z[0] = tail->z;
85 r[1] = head->r; z[1] = head->z; 95 r[1] = head->r; z[1] = head->z;
86 96
87 phiIsOpen = thePhiIsOpen; 97 phiIsOpen = thePhiIsOpen;
88 if (phiIsOpen) 98 if (phiIsOpen)
89 { 99 {
90 deltaPhi = theDeltaPhi; 100 deltaPhi = theDeltaPhi;
91 startPhi = thePhiStart; 101 startPhi = thePhiStart;
92 102
93 // 103 //
94 // Set phi values to our conventions 104 // Set phi values to our conventions
95 // 105 //
96 while (deltaPhi < 0.0) // Loop checking 106 while (deltaPhi < 0.0) // Loop checking, 13.08.2015, G.Cosmo
97 deltaPhi += twopi; 107 deltaPhi += twopi;
98 while (startPhi < 0.0) // Loop checking 108 while (startPhi < 0.0) // Loop checking, 13.08.2015, G.Cosmo
99 startPhi += twopi; 109 startPhi += twopi;
100 110
101 // 111 //
102 // Calculate corner coordinates 112 // Calculate corner coordinates
103 // 113 //
104 ncorners = 4; 114 ncorners = 4;
105 corners = new G4ThreeVector[ncorners]; 115 corners = new G4ThreeVector[ncorners];
106 116
107 corners[0] = G4ThreeVector( tail->r*std::c 117 corners[0] = G4ThreeVector( tail->r*std::cos(startPhi),
108 tail->r*std::s 118 tail->r*std::sin(startPhi), tail->z );
109 corners[1] = G4ThreeVector( head->r*std::c 119 corners[1] = G4ThreeVector( head->r*std::cos(startPhi),
110 head->r*std::s 120 head->r*std::sin(startPhi), head->z );
111 corners[2] = G4ThreeVector( tail->r*std::c 121 corners[2] = G4ThreeVector( tail->r*std::cos(startPhi+deltaPhi),
112 tail->r*std::s 122 tail->r*std::sin(startPhi+deltaPhi), tail->z );
113 corners[3] = G4ThreeVector( head->r*std::c 123 corners[3] = G4ThreeVector( head->r*std::cos(startPhi+deltaPhi),
114 head->r*std::s 124 head->r*std::sin(startPhi+deltaPhi), head->z );
115 } 125 }
116 else 126 else
117 { 127 {
118 deltaPhi = twopi; 128 deltaPhi = twopi;
119 startPhi = 0.0; 129 startPhi = 0.0;
120 } 130 }
121 131
122 allBehind = isAllBehind; 132 allBehind = isAllBehind;
123 133
124 // 134 //
125 // Make our intersecting cone 135 // Make our intersecting cone
126 // 136 //
127 cone = new G4IntersectingCone( r, z ); 137 cone = new G4IntersectingCone( r, z );
128 138
129 // 139 //
130 // Calculate vectors in r,z space 140 // Calculate vectors in r,z space
131 // 141 //
132 rS = r[1]-r[0]; zS = z[1]-z[0]; 142 rS = r[1]-r[0]; zS = z[1]-z[0];
133 length = std::sqrt( rS*rS + zS*zS); 143 length = std::sqrt( rS*rS + zS*zS);
134 rS /= length; zS /= length; 144 rS /= length; zS /= length;
135 145
136 rNorm = +zS; 146 rNorm = +zS;
137 zNorm = -rS; 147 zNorm = -rS;
138 148
139 G4double lAdj; 149 G4double lAdj;
140 150
141 prevRS = r[0]-prevRZ->r; 151 prevRS = r[0]-prevRZ->r;
142 prevZS = z[0]-prevRZ->z; 152 prevZS = z[0]-prevRZ->z;
143 lAdj = std::sqrt( prevRS*prevRS + prevZS*pre 153 lAdj = std::sqrt( prevRS*prevRS + prevZS*prevZS );
144 prevRS /= lAdj; 154 prevRS /= lAdj;
145 prevZS /= lAdj; 155 prevZS /= lAdj;
146 156
147 rNormEdge[0] = rNorm + prevZS; 157 rNormEdge[0] = rNorm + prevZS;
148 zNormEdge[0] = zNorm - prevRS; 158 zNormEdge[0] = zNorm - prevRS;
149 lAdj = std::sqrt( rNormEdge[0]*rNormEdge[0] 159 lAdj = std::sqrt( rNormEdge[0]*rNormEdge[0] + zNormEdge[0]*zNormEdge[0] );
150 rNormEdge[0] /= lAdj; 160 rNormEdge[0] /= lAdj;
151 zNormEdge[0] /= lAdj; 161 zNormEdge[0] /= lAdj;
152 162
153 nextRS = nextRZ->r-r[1]; 163 nextRS = nextRZ->r-r[1];
154 nextZS = nextRZ->z-z[1]; 164 nextZS = nextRZ->z-z[1];
155 lAdj = std::sqrt( nextRS*nextRS + nextZS*nex 165 lAdj = std::sqrt( nextRS*nextRS + nextZS*nextZS );
156 nextRS /= lAdj; 166 nextRS /= lAdj;
157 nextZS /= lAdj; 167 nextZS /= lAdj;
158 168
159 rNormEdge[1] = rNorm + nextZS; 169 rNormEdge[1] = rNorm + nextZS;
160 zNormEdge[1] = zNorm - nextRS; 170 zNormEdge[1] = zNorm - nextRS;
161 lAdj = std::sqrt( rNormEdge[1]*rNormEdge[1] 171 lAdj = std::sqrt( rNormEdge[1]*rNormEdge[1] + zNormEdge[1]*zNormEdge[1] );
162 rNormEdge[1] /= lAdj; 172 rNormEdge[1] /= lAdj;
163 zNormEdge[1] /= lAdj; 173 zNormEdge[1] /= lAdj;
164 } 174 }
165 175
>> 176
>> 177 //
166 // Fake default constructor - sets only member 178 // Fake default constructor - sets only member data and allocates memory
167 // for usage restri 179 // for usage restricted to object persistency.
168 // 180 //
169 G4PolyconeSide::G4PolyconeSide( __void__& ) 181 G4PolyconeSide::G4PolyconeSide( __void__& )
170 : startPhi(0.), deltaPhi(0.), << 182 : startPhi(0.), deltaPhi(0.), phiIsOpen(false), allBehind(false),
171 rNorm(0.), zNorm(0.), rS(0.), zS(0.), leng << 183 cone(0), rNorm(0.), zNorm(0.), rS(0.), zS(0.), length(0.),
172 prevRS(0.), prevZS(0.), nextRS(0.), nextZS << 184 prevRS(0.), prevZS(0.), nextRS(0.), nextZS(0.), ncorners(0), corners(0),
173 kCarTolerance(0.), instanceID(0) << 185 kCarTolerance(0.), fSurfaceArea(0.), instanceID(0)
174 { 186 {
175 r[0] = r[1] = 0.; 187 r[0] = r[1] = 0.;
176 z[0] = z[1] = 0.; 188 z[0] = z[1] = 0.;
177 rNormEdge[0]= rNormEdge[1] = 0.; 189 rNormEdge[0]= rNormEdge[1] = 0.;
178 zNormEdge[0]= zNormEdge[1] = 0.; 190 zNormEdge[0]= zNormEdge[1] = 0.;
179 } 191 }
180 192
>> 193
>> 194 //
181 // Destructor 195 // Destructor
182 // 196 //
183 G4PolyconeSide::~G4PolyconeSide() 197 G4PolyconeSide::~G4PolyconeSide()
184 { 198 {
185 delete cone; 199 delete cone;
186 if (phiIsOpen) { delete [] corners; } 200 if (phiIsOpen) { delete [] corners; }
187 } 201 }
188 202
>> 203
>> 204 //
189 // Copy constructor 205 // Copy constructor
190 // 206 //
191 G4PolyconeSide::G4PolyconeSide( const G4Polyco << 207 G4PolyconeSide::G4PolyconeSide( const G4PolyconeSide &source )
>> 208 : G4VCSGface(), ncorners(0), corners(0)
192 { 209 {
193 instanceID = subInstanceManager.CreateSubIns 210 instanceID = subInstanceManager.CreateSubInstance();
194 211
195 CopyStuff( source ); 212 CopyStuff( source );
196 } 213 }
197 214
>> 215
>> 216 //
198 // Assignment operator 217 // Assignment operator
199 // 218 //
200 G4PolyconeSide& G4PolyconeSide::operator=( con << 219 G4PolyconeSide& G4PolyconeSide::operator=( const G4PolyconeSide &source )
201 { 220 {
202 if (this == &source) { return *this; } 221 if (this == &source) { return *this; }
203 222
204 delete cone; 223 delete cone;
205 if (phiIsOpen) { delete [] corners; } 224 if (phiIsOpen) { delete [] corners; }
206 225
207 CopyStuff( source ); 226 CopyStuff( source );
208 227
209 return *this; 228 return *this;
210 } 229 }
211 230
>> 231
>> 232 //
212 // CopyStuff 233 // CopyStuff
213 // 234 //
214 void G4PolyconeSide::CopyStuff( const G4Polyco << 235 void G4PolyconeSide::CopyStuff( const G4PolyconeSide &source )
215 { 236 {
216 r[0] = source.r[0]; 237 r[0] = source.r[0];
217 r[1] = source.r[1]; 238 r[1] = source.r[1];
218 z[0] = source.z[0]; 239 z[0] = source.z[0];
219 z[1] = source.z[1]; 240 z[1] = source.z[1];
220 241
221 startPhi = source.startPhi; 242 startPhi = source.startPhi;
222 deltaPhi = source.deltaPhi; 243 deltaPhi = source.deltaPhi;
223 phiIsOpen = source.phiIsOpen; 244 phiIsOpen = source.phiIsOpen;
224 allBehind = source.allBehind; 245 allBehind = source.allBehind;
225 246
226 kCarTolerance = source.kCarTolerance; 247 kCarTolerance = source.kCarTolerance;
227 fSurfaceArea = source.fSurfaceArea; 248 fSurfaceArea = source.fSurfaceArea;
228 249
229 cone = new G4IntersectingCone( *source.co 250 cone = new G4IntersectingCone( *source.cone );
230 251
231 rNorm = source.rNorm; 252 rNorm = source.rNorm;
232 zNorm = source.zNorm; 253 zNorm = source.zNorm;
233 rS = source.rS; 254 rS = source.rS;
234 zS = source.zS; 255 zS = source.zS;
235 length = source.length; 256 length = source.length;
236 prevRS = source.prevRS; 257 prevRS = source.prevRS;
237 prevZS = source.prevZS; 258 prevZS = source.prevZS;
238 nextRS = source.nextRS; 259 nextRS = source.nextRS;
239 nextZS = source.nextZS; 260 nextZS = source.nextZS;
240 261
241 rNormEdge[0] = source.rNormEdge[0]; 262 rNormEdge[0] = source.rNormEdge[0];
242 rNormEdge[1] = source.rNormEdge[1]; 263 rNormEdge[1] = source.rNormEdge[1];
243 zNormEdge[0] = source.zNormEdge[0]; 264 zNormEdge[0] = source.zNormEdge[0];
244 zNormEdge[1] = source.zNormEdge[1]; 265 zNormEdge[1] = source.zNormEdge[1];
245 266
246 if (phiIsOpen) 267 if (phiIsOpen)
247 { 268 {
248 ncorners = 4; 269 ncorners = 4;
249 corners = new G4ThreeVector[ncorners]; 270 corners = new G4ThreeVector[ncorners];
250 271
251 corners[0] = source.corners[0]; 272 corners[0] = source.corners[0];
252 corners[1] = source.corners[1]; 273 corners[1] = source.corners[1];
253 corners[2] = source.corners[2]; 274 corners[2] = source.corners[2];
254 corners[3] = source.corners[3]; 275 corners[3] = source.corners[3];
255 } 276 }
256 } 277 }
257 278
>> 279
>> 280 //
258 // Intersect 281 // Intersect
259 // 282 //
260 G4bool G4PolyconeSide::Intersect( const G4Thre << 283 G4bool G4PolyconeSide::Intersect( const G4ThreeVector &p,
261 const G4Thre << 284 const G4ThreeVector &v,
262 G4bool 285 G4bool outgoing,
263 G4doub 286 G4double surfTolerance,
264 G4doub << 287 G4double &distance,
265 G4doub << 288 G4double &distFromSurface,
266 G4Thre << 289 G4ThreeVector &normal,
267 G4bool << 290 G4bool &isAllBehind )
268 { 291 {
269 G4double s1=0., s2=0.; << 292 G4double s1, s2;
270 G4double normSign = outgoing ? +1 : -1; 293 G4double normSign = outgoing ? +1 : -1;
271 294
272 isAllBehind = allBehind; 295 isAllBehind = allBehind;
273 296
274 // 297 //
275 // Check for two possible intersections 298 // Check for two possible intersections
276 // 299 //
277 G4int nside = cone->LineHitsCone( p, v, &s1, 300 G4int nside = cone->LineHitsCone( p, v, &s1, &s2 );
278 if (nside == 0) return false; 301 if (nside == 0) return false;
279 302
280 // 303 //
281 // Check the first side first, since it is ( 304 // Check the first side first, since it is (supposed to be) closest
282 // 305 //
283 G4ThreeVector hit = p + s1*v; 306 G4ThreeVector hit = p + s1*v;
284 307
285 if (PointOnCone( hit, normSign, p, v, normal 308 if (PointOnCone( hit, normSign, p, v, normal ))
286 { 309 {
287 // 310 //
288 // Good intersection! What about the norma 311 // Good intersection! What about the normal?
289 // 312 //
290 if (normSign*v.dot(normal) > 0) 313 if (normSign*v.dot(normal) > 0)
291 { 314 {
292 // 315 //
293 // We have a valid intersection, but it 316 // We have a valid intersection, but it could very easily
294 // be behind the point. To decide if we 317 // be behind the point. To decide if we tolerate this,
295 // we have to see if the point p is on t 318 // we have to see if the point p is on the surface near
296 // the intersecting point. 319 // the intersecting point.
297 // 320 //
298 // What does it mean exactly for the poi 321 // What does it mean exactly for the point p to be "near"
299 // the intersection? It means that if we 322 // the intersection? It means that if we draw a line from
300 // p to the hit, the line remains entire 323 // p to the hit, the line remains entirely within the
301 // tolerance bounds of the cone. To test 324 // tolerance bounds of the cone. To test this, we can
302 // ask if the normal is correct near p. 325 // ask if the normal is correct near p.
303 // 326 //
304 G4double pr = p.perp(); 327 G4double pr = p.perp();
305 if (pr < DBL_MIN) pr = DBL_MIN; 328 if (pr < DBL_MIN) pr = DBL_MIN;
306 G4ThreeVector pNormal( rNorm*p.x()/pr, r 329 G4ThreeVector pNormal( rNorm*p.x()/pr, rNorm*p.y()/pr, zNorm );
307 if (normSign*v.dot(pNormal) > 0) 330 if (normSign*v.dot(pNormal) > 0)
308 { 331 {
309 // 332 //
310 // p and intersection in same hemisphe 333 // p and intersection in same hemisphere
311 // 334 //
312 G4double distOutside2; 335 G4double distOutside2;
313 distFromSurface = -normSign*DistanceAw 336 distFromSurface = -normSign*DistanceAway( p, false, distOutside2 );
314 if (distOutside2 < surfTolerance*surfT 337 if (distOutside2 < surfTolerance*surfTolerance)
315 { 338 {
316 if (distFromSurface > -surfTolerance 339 if (distFromSurface > -surfTolerance)
317 { 340 {
318 // 341 //
319 // We are just inside or away from 342 // We are just inside or away from the
320 // surface. Accept *any* value of 343 // surface. Accept *any* value of distance.
321 // 344 //
322 distance = s1; 345 distance = s1;
323 return true; 346 return true;
324 } 347 }
325 } 348 }
326 } 349 }
327 else 350 else
328 distFromSurface = s1; 351 distFromSurface = s1;
329 352
330 // 353 //
331 // Accept positive distances 354 // Accept positive distances
332 // 355 //
333 if (s1 > 0) 356 if (s1 > 0)
334 { 357 {
335 distance = s1; 358 distance = s1;
336 return true; 359 return true;
337 } 360 }
338 } 361 }
339 } 362 }
340 363
341 if (nside==1) return false; 364 if (nside==1) return false;
342 365
343 // 366 //
344 // Well, try the second hit 367 // Well, try the second hit
345 // 368 //
346 hit = p + s2*v; 369 hit = p + s2*v;
347 370
348 if (PointOnCone( hit, normSign, p, v, normal 371 if (PointOnCone( hit, normSign, p, v, normal ))
349 { 372 {
350 // 373 //
351 // Good intersection! What about the norma 374 // Good intersection! What about the normal?
352 // 375 //
353 if (normSign*v.dot(normal) > 0) 376 if (normSign*v.dot(normal) > 0)
354 { 377 {
355 G4double pr = p.perp(); 378 G4double pr = p.perp();
356 if (pr < DBL_MIN) pr = DBL_MIN; 379 if (pr < DBL_MIN) pr = DBL_MIN;
357 G4ThreeVector pNormal( rNorm*p.x()/pr, r 380 G4ThreeVector pNormal( rNorm*p.x()/pr, rNorm*p.y()/pr, zNorm );
358 if (normSign*v.dot(pNormal) > 0) 381 if (normSign*v.dot(pNormal) > 0)
359 { 382 {
360 G4double distOutside2; 383 G4double distOutside2;
361 distFromSurface = -normSign*DistanceAw 384 distFromSurface = -normSign*DistanceAway( p, false, distOutside2 );
362 if (distOutside2 < surfTolerance*surfT 385 if (distOutside2 < surfTolerance*surfTolerance)
363 { 386 {
364 if (distFromSurface > -surfTolerance 387 if (distFromSurface > -surfTolerance)
365 { 388 {
366 distance = s2; 389 distance = s2;
367 return true; 390 return true;
368 } 391 }
369 } 392 }
370 } 393 }
371 else 394 else
372 distFromSurface = s2; 395 distFromSurface = s2;
373 396
374 if (s2 > 0) 397 if (s2 > 0)
375 { 398 {
376 distance = s2; 399 distance = s2;
377 return true; 400 return true;
378 } 401 }
379 } 402 }
380 } 403 }
381 404
382 // 405 //
383 // Better luck next time 406 // Better luck next time
384 // 407 //
385 return false; 408 return false;
386 } 409 }
387 410
>> 411
>> 412 //
388 // Distance 413 // Distance
389 // 414 //
390 G4double G4PolyconeSide::Distance( const G4Thr << 415 G4double G4PolyconeSide::Distance( const G4ThreeVector &p, G4bool outgoing )
391 { 416 {
392 G4double normSign = outgoing ? -1 : +1; 417 G4double normSign = outgoing ? -1 : +1;
393 G4double distFrom, distOut2; 418 G4double distFrom, distOut2;
394 419
395 // 420 //
396 // We have two tries for each hemisphere. Tr 421 // We have two tries for each hemisphere. Try the closest first.
397 // 422 //
398 distFrom = normSign*DistanceAway( p, false, 423 distFrom = normSign*DistanceAway( p, false, distOut2 );
399 if (distFrom > -0.5*kCarTolerance ) 424 if (distFrom > -0.5*kCarTolerance )
400 { 425 {
401 // 426 //
402 // Good answer 427 // Good answer
403 // 428 //
404 if (distOut2 > 0) 429 if (distOut2 > 0)
405 return std::sqrt( distFrom*distFrom + di 430 return std::sqrt( distFrom*distFrom + distOut2 );
406 else 431 else
407 return std::fabs(distFrom); 432 return std::fabs(distFrom);
408 } 433 }
409 434
410 // 435 //
411 // Try second side. 436 // Try second side.
412 // 437 //
413 distFrom = normSign*DistanceAway( p, true, 438 distFrom = normSign*DistanceAway( p, true, distOut2 );
414 if (distFrom > -0.5*kCarTolerance) 439 if (distFrom > -0.5*kCarTolerance)
415 { 440 {
416 441
417 if (distOut2 > 0) 442 if (distOut2 > 0)
418 return std::sqrt( distFrom*distFrom + di 443 return std::sqrt( distFrom*distFrom + distOut2 );
419 else 444 else
420 return std::fabs(distFrom); 445 return std::fabs(distFrom);
421 } 446 }
422 447
423 return kInfinity; 448 return kInfinity;
424 } 449 }
425 450
>> 451
>> 452 //
426 // Inside 453 // Inside
427 // 454 //
428 EInside G4PolyconeSide::Inside( const G4ThreeV << 455 EInside G4PolyconeSide::Inside( const G4ThreeVector &p,
429 G4double 456 G4double tolerance,
430 G4double << 457 G4double *bestDistance )
431 { 458 {
432 G4double distFrom, distOut2, dist2; 459 G4double distFrom, distOut2, dist2;
433 G4double edgeRZnorm; 460 G4double edgeRZnorm;
434 461
435 distFrom = DistanceAway( p, distOut2, &edge 462 distFrom = DistanceAway( p, distOut2, &edgeRZnorm );
436 dist2 = distFrom*distFrom + distOut2; 463 dist2 = distFrom*distFrom + distOut2;
437 464
438 *bestDistance = std::sqrt( dist2); 465 *bestDistance = std::sqrt( dist2);
439 466
440 // Okay then, inside or out? 467 // Okay then, inside or out?
441 // 468 //
442 if ( (std::fabs(edgeRZnorm) < tolerance) 469 if ( (std::fabs(edgeRZnorm) < tolerance)
443 && (distOut2< tolerance*tolerance) ) 470 && (distOut2< tolerance*tolerance) )
444 return kSurface; 471 return kSurface;
445 else if (edgeRZnorm < 0) 472 else if (edgeRZnorm < 0)
446 return kInside; 473 return kInside;
447 else 474 else
448 return kOutside; 475 return kOutside;
449 } 476 }
450 477
>> 478
>> 479 //
451 // Normal 480 // Normal
452 // 481 //
453 G4ThreeVector G4PolyconeSide::Normal( const G4 << 482 G4ThreeVector G4PolyconeSide::Normal( const G4ThreeVector &p,
454 G4 << 483 G4double *bestDistance )
455 { 484 {
456 if (p == G4ThreeVector(0.,0.,0.)) { return 485 if (p == G4ThreeVector(0.,0.,0.)) { return p; }
457 486
458 G4double dFrom, dOut2; 487 G4double dFrom, dOut2;
459 488
460 dFrom = DistanceAway( p, false, dOut2 ); 489 dFrom = DistanceAway( p, false, dOut2 );
461 490
462 *bestDistance = std::sqrt( dFrom*dFrom + dOu 491 *bestDistance = std::sqrt( dFrom*dFrom + dOut2 );
463 492
464 G4double rds = p.perp(); 493 G4double rds = p.perp();
465 if (rds!=0.) { return {rNorm*p.x()/rds,rNorm << 494 if (rds!=0.) { return G4ThreeVector(rNorm*p.x()/rds,rNorm*p.y()/rds,zNorm); }
466 return G4ThreeVector( 0.,0., zNorm ).unit(); 495 return G4ThreeVector( 0.,0., zNorm ).unit();
467 } 496 }
468 497
>> 498
>> 499 //
469 // Extent 500 // Extent
470 // 501 //
471 G4double G4PolyconeSide::Extent( const G4Three 502 G4double G4PolyconeSide::Extent( const G4ThreeVector axis )
472 { 503 {
473 if (axis.perp2() < DBL_MIN) 504 if (axis.perp2() < DBL_MIN)
474 { 505 {
475 // 506 //
476 // Special case 507 // Special case
477 // 508 //
478 return axis.z() < 0 ? -cone->ZLo() : cone- 509 return axis.z() < 0 ? -cone->ZLo() : cone->ZHi();
479 } 510 }
480 511
481 // 512 //
482 // Is the axis pointing inside our phi gap? 513 // Is the axis pointing inside our phi gap?
483 // 514 //
484 if (phiIsOpen) 515 if (phiIsOpen)
485 { 516 {
486 G4double phi = GetPhi(axis); 517 G4double phi = GetPhi(axis);
487 while( phi < startPhi ) // Loop checkin 518 while( phi < startPhi ) // Loop checking, 13.08.2015, G.Cosmo
488 phi += twopi; 519 phi += twopi;
489 520
490 if (phi > deltaPhi+startPhi) 521 if (phi > deltaPhi+startPhi)
491 { 522 {
492 // 523 //
493 // Yeah, looks so. Make four three vecto 524 // Yeah, looks so. Make four three vectors defining the phi
494 // opening 525 // opening
495 // 526 //
496 G4double cosP = std::cos(startPhi), sinP 527 G4double cosP = std::cos(startPhi), sinP = std::sin(startPhi);
497 G4ThreeVector a( r[0]*cosP, r[0]*sinP, z 528 G4ThreeVector a( r[0]*cosP, r[0]*sinP, z[0] );
498 G4ThreeVector b( r[1]*cosP, r[1]*sinP, z 529 G4ThreeVector b( r[1]*cosP, r[1]*sinP, z[1] );
499 cosP = std::cos(startPhi+deltaPhi); sinP 530 cosP = std::cos(startPhi+deltaPhi); sinP = std::sin(startPhi+deltaPhi);
500 G4ThreeVector c( r[0]*cosP, r[0]*sinP, z 531 G4ThreeVector c( r[0]*cosP, r[0]*sinP, z[0] );
501 G4ThreeVector d( r[1]*cosP, r[1]*sinP, z 532 G4ThreeVector d( r[1]*cosP, r[1]*sinP, z[1] );
502 533
503 G4double ad = axis.dot(a), 534 G4double ad = axis.dot(a),
504 bd = axis.dot(b), << 535 bd = axis.dot(b),
505 cd = axis.dot(c), << 536 cd = axis.dot(c),
506 dd = axis.dot(d); << 537 dd = axis.dot(d);
507 538
508 if (bd > ad) ad = bd; 539 if (bd > ad) ad = bd;
509 if (cd > ad) ad = cd; 540 if (cd > ad) ad = cd;
510 if (dd > ad) ad = dd; 541 if (dd > ad) ad = dd;
511 542
512 return ad; 543 return ad;
513 } 544 }
514 } 545 }
515 546
516 // 547 //
517 // Check either end 548 // Check either end
518 // 549 //
519 G4double aPerp = axis.perp(); 550 G4double aPerp = axis.perp();
520 551
521 G4double a = aPerp*r[0] + axis.z()*z[0]; 552 G4double a = aPerp*r[0] + axis.z()*z[0];
522 G4double b = aPerp*r[1] + axis.z()*z[1]; 553 G4double b = aPerp*r[1] + axis.z()*z[1];
523 554
524 if (b > a) a = b; 555 if (b > a) a = b;
525 556
526 return a; 557 return a;
527 } 558 }
528 559
>> 560
>> 561 //
529 // CalculateExtent 562 // CalculateExtent
530 // 563 //
531 // See notes in G4VCSGface 564 // See notes in G4VCSGface
532 // 565 //
533 void G4PolyconeSide::CalculateExtent( const EA 566 void G4PolyconeSide::CalculateExtent( const EAxis axis,
534 const G4 << 567 const G4VoxelLimits &voxelLimit,
535 const G4 << 568 const G4AffineTransform &transform,
536 G4 << 569 G4SolidExtentList &extentList )
537 { 570 {
538 G4ClippablePolygon polygon; 571 G4ClippablePolygon polygon;
539 572
540 // 573 //
541 // Here we will approximate (ala G4Cons) and 574 // Here we will approximate (ala G4Cons) and divide our conical section
542 // into segments, like G4Polyhedra. When doi 575 // into segments, like G4Polyhedra. When doing so, the radius
543 // is extented far enough such that the segm 576 // is extented far enough such that the segments always lie
544 // just outside the surface of the conical s 577 // just outside the surface of the conical section we are
545 // approximating. 578 // approximating.
546 // 579 //
547 580
548 // 581 //
549 // Choose phi size of our segment(s) based o 582 // Choose phi size of our segment(s) based on constants as
550 // defined in meshdefs.hh 583 // defined in meshdefs.hh
551 // 584 //
552 G4int numPhi = (G4int)(deltaPhi/kMeshAngleDe 585 G4int numPhi = (G4int)(deltaPhi/kMeshAngleDefault) + 1;
553 if (numPhi < kMinMeshSections) 586 if (numPhi < kMinMeshSections)
554 numPhi = kMinMeshSections; 587 numPhi = kMinMeshSections;
555 else if (numPhi > kMaxMeshSections) 588 else if (numPhi > kMaxMeshSections)
556 numPhi = kMaxMeshSections; 589 numPhi = kMaxMeshSections;
557 590
558 G4double sigPhi = deltaPhi/numPhi; 591 G4double sigPhi = deltaPhi/numPhi;
559 592
560 // 593 //
561 // Determine radius factor to keep segments 594 // Determine radius factor to keep segments outside
562 // 595 //
563 G4double rFudge = 1.0/std::cos(0.5*sigPhi); 596 G4double rFudge = 1.0/std::cos(0.5*sigPhi);
564 597
565 // 598 //
566 // Decide which radius to use on each end of 599 // Decide which radius to use on each end of the side,
567 // and whether a transition mesh is required 600 // and whether a transition mesh is required
568 // 601 //
569 // {r0,z0} - Beginning of this side 602 // {r0,z0} - Beginning of this side
570 // {r1,z1} - Ending of this side 603 // {r1,z1} - Ending of this side
571 // {r2,z0} - Beginning of transition piece 604 // {r2,z0} - Beginning of transition piece connecting previous
572 // side (and ends at beginning of 605 // side (and ends at beginning of this side)
573 // 606 //
574 // So, order is 2 --> 0 --> 1. 607 // So, order is 2 --> 0 --> 1.
575 // ------- 608 // -------
576 // 609 //
577 // r2 < 0 indicates that no transition piece 610 // r2 < 0 indicates that no transition piece is required
578 // 611 //
579 G4double r0, r1, r2, z0, z1; 612 G4double r0, r1, r2, z0, z1;
580 613
581 r2 = -1; // By default: no transition piece 614 r2 = -1; // By default: no transition piece
582 615
583 if (rNorm < -DBL_MIN) 616 if (rNorm < -DBL_MIN)
584 { 617 {
585 // 618 //
586 // This side faces *inward*, and so our me 619 // This side faces *inward*, and so our mesh has
587 // the same radius 620 // the same radius
588 // 621 //
589 r1 = r[1]; 622 r1 = r[1];
590 z1 = z[1]; 623 z1 = z[1];
591 z0 = z[0]; 624 z0 = z[0];
592 r0 = r[0]; 625 r0 = r[0];
593 626
594 r2 = -1; 627 r2 = -1;
595 628
596 if (prevZS > DBL_MIN) 629 if (prevZS > DBL_MIN)
597 { 630 {
598 // 631 //
599 // The previous side is facing outwards 632 // The previous side is facing outwards
600 // 633 //
601 if ( prevRS*zS - prevZS*rS > 0 ) 634 if ( prevRS*zS - prevZS*rS > 0 )
602 { 635 {
603 // 636 //
604 // Transition was convex: build transi 637 // Transition was convex: build transition piece
605 // 638 //
606 if (r[0] > DBL_MIN) r2 = r[0]*rFudge; 639 if (r[0] > DBL_MIN) r2 = r[0]*rFudge;
607 } 640 }
608 else 641 else
609 { 642 {
610 // 643 //
611 // Transition was concave: short this 644 // Transition was concave: short this side
612 // 645 //
613 FindLineIntersect( z0, r0, zS, rS, 646 FindLineIntersect( z0, r0, zS, rS,
614 z0, r0*rFudge, prev 647 z0, r0*rFudge, prevZS, prevRS*rFudge, z0, r0 );
615 } 648 }
616 } 649 }
617 650
618 if ( nextZS > DBL_MIN && (rS*nextZS - zS*n 651 if ( nextZS > DBL_MIN && (rS*nextZS - zS*nextRS < 0) )
619 { 652 {
620 // 653 //
621 // The next side is facing outwards, for 654 // The next side is facing outwards, forming a
622 // concave transition: short this side 655 // concave transition: short this side
623 // 656 //
624 FindLineIntersect( z1, r1, zS, rS, 657 FindLineIntersect( z1, r1, zS, rS,
625 z1, r1*rFudge, nextZS 658 z1, r1*rFudge, nextZS, nextRS*rFudge, z1, r1 );
626 } 659 }
627 } 660 }
628 else if (rNorm > DBL_MIN) 661 else if (rNorm > DBL_MIN)
629 { 662 {
630 // 663 //
631 // This side faces *outward* and is given 664 // This side faces *outward* and is given a boost to
632 // it radius 665 // it radius
633 // 666 //
634 r0 = r[0]*rFudge; 667 r0 = r[0]*rFudge;
635 z0 = z[0]; 668 z0 = z[0];
636 r1 = r[1]*rFudge; 669 r1 = r[1]*rFudge;
637 z1 = z[1]; 670 z1 = z[1];
638 671
639 if (prevZS < -DBL_MIN) 672 if (prevZS < -DBL_MIN)
640 { 673 {
641 // 674 //
642 // The previous side is facing inwards 675 // The previous side is facing inwards
643 // 676 //
644 if ( prevRS*zS - prevZS*rS > 0 ) 677 if ( prevRS*zS - prevZS*rS > 0 )
645 { 678 {
646 // 679 //
647 // Transition was convex: build transi 680 // Transition was convex: build transition piece
648 // 681 //
649 if (r[0] > DBL_MIN) r2 = r[0]; 682 if (r[0] > DBL_MIN) r2 = r[0];
650 } 683 }
651 else 684 else
652 { 685 {
653 // 686 //
654 // Transition was concave: short this 687 // Transition was concave: short this side
655 // 688 //
656 FindLineIntersect( z0, r0, zS, rS*rFud 689 FindLineIntersect( z0, r0, zS, rS*rFudge,
657 z0, r[0], prevZS, p 690 z0, r[0], prevZS, prevRS, z0, r0 );
658 } 691 }
659 } 692 }
660 693
661 if ( nextZS < -DBL_MIN && (rS*nextZS - zS* 694 if ( nextZS < -DBL_MIN && (rS*nextZS - zS*nextRS < 0) )
662 { 695 {
663 // 696 //
664 // The next side is facing inwards, form 697 // The next side is facing inwards, forming a
665 // concave transition: short this side 698 // concave transition: short this side
666 // 699 //
667 FindLineIntersect( z1, r1, zS, rS*rFudge 700 FindLineIntersect( z1, r1, zS, rS*rFudge,
668 z1, r[1], nextZS, nex 701 z1, r[1], nextZS, nextRS, z1, r1 );
669 } 702 }
670 } 703 }
671 else 704 else
672 { 705 {
673 // 706 //
674 // This side is perpendicular to the z axi 707 // This side is perpendicular to the z axis (is a disk)
675 // 708 //
676 // Whether or not r0 needs a rFudge factor 709 // Whether or not r0 needs a rFudge factor depends
677 // on the normal of the previous edge. Sim 710 // on the normal of the previous edge. Similar with r1
678 // and the next edge. No transition piece 711 // and the next edge. No transition piece is required.
679 // 712 //
680 r0 = r[0]; 713 r0 = r[0];
681 r1 = r[1]; 714 r1 = r[1];
682 z0 = z[0]; 715 z0 = z[0];
683 z1 = z[1]; 716 z1 = z[1];
684 717
685 if (prevZS > DBL_MIN) r0 *= rFudge; 718 if (prevZS > DBL_MIN) r0 *= rFudge;
686 if (nextZS > DBL_MIN) r1 *= rFudge; 719 if (nextZS > DBL_MIN) r1 *= rFudge;
687 } 720 }
688 721
689 // 722 //
690 // Loop 723 // Loop
691 // 724 //
692 G4double phi = startPhi, 725 G4double phi = startPhi,
693 cosPhi = std::cos(phi), 726 cosPhi = std::cos(phi),
694 sinPhi = std::sin(phi); 727 sinPhi = std::sin(phi);
695 728
696 G4ThreeVector v0( r0*cosPhi, r0*sinPhi, z0 ) 729 G4ThreeVector v0( r0*cosPhi, r0*sinPhi, z0 ),
697 v1( r1*cosPhi, r1*sinPhi, 730 v1( r1*cosPhi, r1*sinPhi, z1 ),
698 v2, w0, w1, w2; 731 v2, w0, w1, w2;
699 transform.ApplyPointTransform( v0 ); 732 transform.ApplyPointTransform( v0 );
700 transform.ApplyPointTransform( v1 ); 733 transform.ApplyPointTransform( v1 );
701 734
702 if (r2 >= 0) 735 if (r2 >= 0)
703 { 736 {
704 v2 = G4ThreeVector( r2*cosPhi, r2*sinPhi, 737 v2 = G4ThreeVector( r2*cosPhi, r2*sinPhi, z0 );
705 transform.ApplyPointTransform( v2 ); 738 transform.ApplyPointTransform( v2 );
706 } 739 }
707 740
708 do // Loop checking, 13.08.2015, G.Cosmo 741 do // Loop checking, 13.08.2015, G.Cosmo
709 { 742 {
710 phi += sigPhi; 743 phi += sigPhi;
711 if (numPhi == 1) phi = startPhi+deltaPhi; 744 if (numPhi == 1) phi = startPhi+deltaPhi; // Try to avoid roundoff
712 cosPhi = std::cos(phi), 745 cosPhi = std::cos(phi),
713 sinPhi = std::sin(phi); 746 sinPhi = std::sin(phi);
714 747
715 w0 = G4ThreeVector( r0*cosPhi, r0*sinPhi, 748 w0 = G4ThreeVector( r0*cosPhi, r0*sinPhi, z0 );
716 w1 = G4ThreeVector( r1*cosPhi, r1*sinPhi, 749 w1 = G4ThreeVector( r1*cosPhi, r1*sinPhi, z1 );
717 transform.ApplyPointTransform( w0 ); 750 transform.ApplyPointTransform( w0 );
718 transform.ApplyPointTransform( w1 ); 751 transform.ApplyPointTransform( w1 );
719 752
720 G4ThreeVector deltaV = r0 > r1 ? w0-v0 : w 753 G4ThreeVector deltaV = r0 > r1 ? w0-v0 : w1-v1;
721 754
722 // 755 //
723 // Build polygon, taking special care to k 756 // Build polygon, taking special care to keep the vertices
724 // in order 757 // in order
725 // 758 //
726 polygon.ClearAllVertices(); 759 polygon.ClearAllVertices();
727 760
728 polygon.AddVertexInOrder( v0 ); 761 polygon.AddVertexInOrder( v0 );
729 polygon.AddVertexInOrder( v1 ); 762 polygon.AddVertexInOrder( v1 );
730 polygon.AddVertexInOrder( w1 ); 763 polygon.AddVertexInOrder( w1 );
731 polygon.AddVertexInOrder( w0 ); 764 polygon.AddVertexInOrder( w0 );
732 765
733 // 766 //
734 // Get extent 767 // Get extent
735 // 768 //
736 if (polygon.PartialClip( voxelLimit, axis 769 if (polygon.PartialClip( voxelLimit, axis ))
737 { 770 {
738 // 771 //
739 // Get dot product of normal with target 772 // Get dot product of normal with target axis
740 // 773 //
741 polygon.SetNormal( deltaV.cross(v1-v0).u 774 polygon.SetNormal( deltaV.cross(v1-v0).unit() );
742 775
743 extentList.AddSurface( polygon ); 776 extentList.AddSurface( polygon );
744 } 777 }
745 778
746 if (r2 >= 0) 779 if (r2 >= 0)
747 { 780 {
748 // 781 //
749 // Repeat, for transition piece 782 // Repeat, for transition piece
750 // 783 //
751 w2 = G4ThreeVector( r2*cosPhi, r2*sinPhi 784 w2 = G4ThreeVector( r2*cosPhi, r2*sinPhi, z0 );
752 transform.ApplyPointTransform( w2 ); 785 transform.ApplyPointTransform( w2 );
753 786
754 polygon.ClearAllVertices(); 787 polygon.ClearAllVertices();
755 788
756 polygon.AddVertexInOrder( v2 ); 789 polygon.AddVertexInOrder( v2 );
757 polygon.AddVertexInOrder( v0 ); 790 polygon.AddVertexInOrder( v0 );
758 polygon.AddVertexInOrder( w0 ); 791 polygon.AddVertexInOrder( w0 );
759 polygon.AddVertexInOrder( w2 ); 792 polygon.AddVertexInOrder( w2 );
760 793
761 if (polygon.PartialClip( voxelLimit, axi 794 if (polygon.PartialClip( voxelLimit, axis ))
762 { 795 {
763 polygon.SetNormal( deltaV.cross(v0-v2) 796 polygon.SetNormal( deltaV.cross(v0-v2).unit() );
764 797
765 extentList.AddSurface( polygon ); 798 extentList.AddSurface( polygon );
766 } 799 }
767 800
768 v2 = w2; 801 v2 = w2;
769 } 802 }
770 803
771 // 804 //
772 // Next vertex 805 // Next vertex
773 // 806 //
774 v0 = w0; 807 v0 = w0;
775 v1 = w1; 808 v1 = w1;
776 } while( --numPhi > 0 ); 809 } while( --numPhi > 0 );
777 810
778 // 811 //
779 // We are almost done. But, it is important 812 // We are almost done. But, it is important that we leave no
780 // gaps in the surface of our solid. By usin 813 // gaps in the surface of our solid. By using rFudge, however,
781 // we've done exactly that, if we have a phi 814 // we've done exactly that, if we have a phi segment.
782 // Add two additional faces if necessary 815 // Add two additional faces if necessary
783 // 816 //
784 if (phiIsOpen && rNorm > DBL_MIN) 817 if (phiIsOpen && rNorm > DBL_MIN)
785 { 818 {
786 cosPhi = std::cos(startPhi); 819 cosPhi = std::cos(startPhi);
787 sinPhi = std::sin(startPhi); 820 sinPhi = std::sin(startPhi);
788 821
789 G4ThreeVector a0( r[0]*cosPhi, r[0]*sinPhi 822 G4ThreeVector a0( r[0]*cosPhi, r[0]*sinPhi, z[0] ),
790 a1( r[1]*cosPhi, r[1]*sinPhi 823 a1( r[1]*cosPhi, r[1]*sinPhi, z[1] ),
791 b0( r0*cosPhi, r0*sinPhi, z[ 824 b0( r0*cosPhi, r0*sinPhi, z[0] ),
792 b1( r1*cosPhi, r1*sinPhi, z[ 825 b1( r1*cosPhi, r1*sinPhi, z[1] );
793 826
794 transform.ApplyPointTransform( a0 ); 827 transform.ApplyPointTransform( a0 );
795 transform.ApplyPointTransform( a1 ); 828 transform.ApplyPointTransform( a1 );
796 transform.ApplyPointTransform( b0 ); 829 transform.ApplyPointTransform( b0 );
797 transform.ApplyPointTransform( b1 ); 830 transform.ApplyPointTransform( b1 );
798 831
799 polygon.ClearAllVertices(); 832 polygon.ClearAllVertices();
800 833
801 polygon.AddVertexInOrder( a0 ); 834 polygon.AddVertexInOrder( a0 );
802 polygon.AddVertexInOrder( a1 ); 835 polygon.AddVertexInOrder( a1 );
803 polygon.AddVertexInOrder( b0 ); 836 polygon.AddVertexInOrder( b0 );
804 polygon.AddVertexInOrder( b1 ); 837 polygon.AddVertexInOrder( b1 );
805 838
806 if (polygon.PartialClip( voxelLimit , axis 839 if (polygon.PartialClip( voxelLimit , axis))
807 { 840 {
808 G4ThreeVector normal( sinPhi, -cosPhi, 0 841 G4ThreeVector normal( sinPhi, -cosPhi, 0 );
809 polygon.SetNormal( transform.TransformAx 842 polygon.SetNormal( transform.TransformAxis( normal ) );
810 843
811 extentList.AddSurface( polygon ); 844 extentList.AddSurface( polygon );
812 } 845 }
813 846
814 cosPhi = std::cos(startPhi+deltaPhi); 847 cosPhi = std::cos(startPhi+deltaPhi);
815 sinPhi = std::sin(startPhi+deltaPhi); 848 sinPhi = std::sin(startPhi+deltaPhi);
816 849
817 a0 = G4ThreeVector( r[0]*cosPhi, r[0]*sinP 850 a0 = G4ThreeVector( r[0]*cosPhi, r[0]*sinPhi, z[0] ),
818 a1 = G4ThreeVector( r[1]*cosPhi, r[1]*sinP 851 a1 = G4ThreeVector( r[1]*cosPhi, r[1]*sinPhi, z[1] ),
819 b0 = G4ThreeVector( r0*cosPhi, r0*sinPhi, 852 b0 = G4ThreeVector( r0*cosPhi, r0*sinPhi, z[0] ),
820 b1 = G4ThreeVector( r1*cosPhi, r1*sinPhi, 853 b1 = G4ThreeVector( r1*cosPhi, r1*sinPhi, z[1] );
821 transform.ApplyPointTransform( a0 ); 854 transform.ApplyPointTransform( a0 );
822 transform.ApplyPointTransform( a1 ); 855 transform.ApplyPointTransform( a1 );
823 transform.ApplyPointTransform( b0 ); 856 transform.ApplyPointTransform( b0 );
824 transform.ApplyPointTransform( b1 ); 857 transform.ApplyPointTransform( b1 );
825 858
826 polygon.ClearAllVertices(); 859 polygon.ClearAllVertices();
827 860
828 polygon.AddVertexInOrder( a0 ); 861 polygon.AddVertexInOrder( a0 );
829 polygon.AddVertexInOrder( a1 ); 862 polygon.AddVertexInOrder( a1 );
830 polygon.AddVertexInOrder( b0 ); 863 polygon.AddVertexInOrder( b0 );
831 polygon.AddVertexInOrder( b1 ); 864 polygon.AddVertexInOrder( b1 );
832 865
833 if (polygon.PartialClip( voxelLimit, axis 866 if (polygon.PartialClip( voxelLimit, axis ))
834 { 867 {
835 G4ThreeVector normal( -sinPhi, cosPhi, 0 868 G4ThreeVector normal( -sinPhi, cosPhi, 0 );
836 polygon.SetNormal( transform.TransformAx 869 polygon.SetNormal( transform.TransformAxis( normal ) );
837 870
838 extentList.AddSurface( polygon ); 871 extentList.AddSurface( polygon );
839 } 872 }
840 } 873 }
841 874
842 return; 875 return;
843 } 876 }
844 877
>> 878
>> 879 //
845 // GetPhi 880 // GetPhi
846 // 881 //
847 // Calculate Phi for a given 3-vector (point), 882 // Calculate Phi for a given 3-vector (point), if not already cached for the
848 // same point, in the attempt to avoid consecu 883 // same point, in the attempt to avoid consecutive computation of the same
849 // quantity 884 // quantity
850 // 885 //
851 G4double G4PolyconeSide::GetPhi( const G4Three 886 G4double G4PolyconeSide::GetPhi( const G4ThreeVector& p )
852 { 887 {
853 G4double val=0.; 888 G4double val=0.;
854 G4ThreeVector vphi(G4MT_pcphix, G4MT_pcphiy, <<
855 889
856 if (vphi != p) << 890 if (G4MT_pcphi.first != p)
857 { 891 {
858 val = p.phi(); 892 val = p.phi();
859 G4MT_pcphix = p.x(); G4MT_pcphiy = p.y(); << 893 G4MT_pcphi.first = p;
860 G4MT_pcphik = val; << 894 G4MT_pcphi.second = val;
861 } 895 }
862 else 896 else
863 { 897 {
864 val = G4MT_pcphik; << 898 val = G4MT_pcphi.second;
865 } 899 }
866 return val; 900 return val;
867 } 901 }
868 902
>> 903
>> 904 //
869 // DistanceAway 905 // DistanceAway
870 // 906 //
871 // Calculate distance of a point from our coni 907 // Calculate distance of a point from our conical surface, including the effect
872 // of any phi segmentation 908 // of any phi segmentation
873 // 909 //
874 // Arguments: 910 // Arguments:
875 // p - (in) Point to check 911 // p - (in) Point to check
876 // opposite - (in) If true, check opposi 912 // opposite - (in) If true, check opposite hemisphere (see below)
877 // distOutside - (out) Additional distance 913 // distOutside - (out) Additional distance outside the edges of the surface
878 // edgeRZnorm - (out) if negative, point i 914 // edgeRZnorm - (out) if negative, point is inside
879 // 915 //
880 // return value = distance from the conical p 916 // return value = distance from the conical plane, if extrapolated beyond edges,
881 // signed by whether the point 917 // signed by whether the point is in inside or outside the shape
882 // 918 //
883 // Notes: 919 // Notes:
884 // * There are two answers, depending on whic 920 // * There are two answers, depending on which hemisphere is considered.
885 // 921 //
886 G4double G4PolyconeSide::DistanceAway( const G << 922 G4double G4PolyconeSide::DistanceAway( const G4ThreeVector &p,
887 G 923 G4bool opposite,
888 G << 924 G4double &distOutside2,
889 G << 925 G4double *edgeRZnorm )
890 { 926 {
891 // 927 //
892 // Convert our point to r and z 928 // Convert our point to r and z
893 // 929 //
894 G4double rx = p.perp(), zx = p.z(); 930 G4double rx = p.perp(), zx = p.z();
895 931
896 // 932 //
897 // Change sign of r if opposite says we shou 933 // Change sign of r if opposite says we should
898 // 934 //
899 if (opposite) rx = -rx; 935 if (opposite) rx = -rx;
900 936
901 // 937 //
902 // Calculate return value 938 // Calculate return value
903 // 939 //
904 G4double deltaR = rx - r[0], deltaZ = zx - 940 G4double deltaR = rx - r[0], deltaZ = zx - z[0];
905 G4double answer = deltaR*rNorm + deltaZ*zNor 941 G4double answer = deltaR*rNorm + deltaZ*zNorm;
906 942
907 // 943 //
908 // Are we off the surface in r,z space? 944 // Are we off the surface in r,z space?
909 // 945 //
910 G4double q = deltaR*rS + deltaZ*zS; 946 G4double q = deltaR*rS + deltaZ*zS;
911 if (q < 0) 947 if (q < 0)
912 { 948 {
913 distOutside2 = q*q; 949 distOutside2 = q*q;
914 if (edgeRZnorm != nullptr) << 950 if (edgeRZnorm) *edgeRZnorm = deltaR*rNormEdge[0] + deltaZ*zNormEdge[0];
915 *edgeRZnorm = deltaR*rNormEdge[0] + delt <<
916 } 951 }
917 else if (q > length) 952 else if (q > length)
918 { 953 {
919 distOutside2 = sqr( q-length ); 954 distOutside2 = sqr( q-length );
920 if (edgeRZnorm != nullptr) << 955 if (edgeRZnorm)
921 { 956 {
922 deltaR = rx - r[1]; 957 deltaR = rx - r[1];
923 deltaZ = zx - z[1]; 958 deltaZ = zx - z[1];
924 *edgeRZnorm = deltaR*rNormEdge[1] + delt 959 *edgeRZnorm = deltaR*rNormEdge[1] + deltaZ*zNormEdge[1];
925 } 960 }
926 } 961 }
927 else 962 else
928 { 963 {
929 distOutside2 = 0.; << 964 distOutside2 = 0;
930 if (edgeRZnorm != nullptr) *edgeRZnorm = a << 965 if (edgeRZnorm) *edgeRZnorm = answer;
931 } 966 }
932 967
933 if (phiIsOpen) 968 if (phiIsOpen)
934 { 969 {
935 // 970 //
936 // Finally, check phi 971 // Finally, check phi
937 // 972 //
938 G4double phi = GetPhi(p); 973 G4double phi = GetPhi(p);
939 while( phi < startPhi ) // Loop checkin 974 while( phi < startPhi ) // Loop checking, 13.08.2015, G.Cosmo
940 phi += twopi; 975 phi += twopi;
941 976
942 if (phi > startPhi+deltaPhi) 977 if (phi > startPhi+deltaPhi)
943 { 978 {
944 // 979 //
945 // Oops. Are we closer to the start phi 980 // Oops. Are we closer to the start phi or end phi?
946 // 981 //
947 G4double d1 = phi-startPhi-deltaPhi; 982 G4double d1 = phi-startPhi-deltaPhi;
948 while( phi > startPhi ) // Loop check 983 while( phi > startPhi ) // Loop checking, 13.08.2015, G.Cosmo
949 phi -= twopi; 984 phi -= twopi;
950 G4double d2 = startPhi-phi; 985 G4double d2 = startPhi-phi;
951 986
952 if (d2 < d1) d1 = d2; 987 if (d2 < d1) d1 = d2;
953 988
954 // 989 //
955 // Add result to our distance 990 // Add result to our distance
956 // 991 //
957 G4double dist = d1*rx; 992 G4double dist = d1*rx;
958 993
959 distOutside2 += dist*dist; 994 distOutside2 += dist*dist;
960 if (edgeRZnorm != nullptr) << 995 if (edgeRZnorm)
961 { 996 {
962 *edgeRZnorm = std::max(std::fabs(*edge 997 *edgeRZnorm = std::max(std::fabs(*edgeRZnorm),std::fabs(dist));
963 } 998 }
964 } 999 }
965 } 1000 }
966 1001
967 return answer; 1002 return answer;
968 } 1003 }
969 1004
>> 1005
>> 1006 //
970 // DistanceAway 1007 // DistanceAway
971 // 1008 //
972 // Special version of DistanceAway for Inside. 1009 // Special version of DistanceAway for Inside.
973 // Opposite parameter is not used, instead use 1010 // Opposite parameter is not used, instead use sign of rx for choosing the side
974 // 1011 //
975 G4double G4PolyconeSide::DistanceAway( const G << 1012 G4double G4PolyconeSide::DistanceAway( const G4ThreeVector &p,
976 G << 1013 G4double &distOutside2,
977 G << 1014 G4double *edgeRZnorm )
978 { 1015 {
979 // 1016 //
980 // Convert our point to r and z 1017 // Convert our point to r and z
981 // 1018 //
982 G4double rx = p.perp(), zx = p.z(); 1019 G4double rx = p.perp(), zx = p.z();
983 1020
984 // 1021 //
985 // Change sign of r if we should 1022 // Change sign of r if we should
986 // 1023 //
987 G4int part = 1; 1024 G4int part = 1;
988 if (rx < 0) part = -1; 1025 if (rx < 0) part = -1;
989 1026
990 // 1027 //
991 // Calculate return value 1028 // Calculate return value
992 // 1029 //
993 G4double deltaR = rx - r[0]*part, deltaZ = z << 1030 G4double deltaR = rx - r[0]*part, deltaZ = zx - z[0];
994 G4double answer = deltaR*rNorm*part + deltaZ 1031 G4double answer = deltaR*rNorm*part + deltaZ*zNorm;
995 1032
996 // 1033 //
997 // Are we off the surface in r,z space? 1034 // Are we off the surface in r,z space?
998 // 1035 //
999 G4double q = deltaR*rS*part + deltaZ*zS; 1036 G4double q = deltaR*rS*part + deltaZ*zS;
1000 if (q < 0) 1037 if (q < 0)
1001 { 1038 {
1002 distOutside2 = q*q; 1039 distOutside2 = q*q;
1003 if (edgeRZnorm != nullptr) << 1040 if (edgeRZnorm)
1004 { 1041 {
1005 *edgeRZnorm = deltaR*rNormEdge[0]*part 1042 *edgeRZnorm = deltaR*rNormEdge[0]*part + deltaZ*zNormEdge[0];
1006 } 1043 }
1007 } 1044 }
1008 else if (q > length) 1045 else if (q > length)
1009 { 1046 {
1010 distOutside2 = sqr( q-length ); 1047 distOutside2 = sqr( q-length );
1011 if (edgeRZnorm != nullptr) << 1048 if (edgeRZnorm)
1012 { 1049 {
1013 deltaR = rx - r[1]*part; 1050 deltaR = rx - r[1]*part;
1014 deltaZ = zx - z[1]; 1051 deltaZ = zx - z[1];
1015 *edgeRZnorm = deltaR*rNormEdge[1]*part 1052 *edgeRZnorm = deltaR*rNormEdge[1]*part + deltaZ*zNormEdge[1];
1016 } 1053 }
1017 } 1054 }
1018 else 1055 else
1019 { 1056 {
1020 distOutside2 = 0.; << 1057 distOutside2 = 0;
1021 if (edgeRZnorm != nullptr) *edgeRZnorm = << 1058 if (edgeRZnorm) *edgeRZnorm = answer;
1022 } 1059 }
1023 1060
1024 if (phiIsOpen) 1061 if (phiIsOpen)
1025 { 1062 {
1026 // 1063 //
1027 // Finally, check phi 1064 // Finally, check phi
1028 // 1065 //
1029 G4double phi = GetPhi(p); 1066 G4double phi = GetPhi(p);
1030 while( phi < startPhi ) // Loop checki 1067 while( phi < startPhi ) // Loop checking, 13.08.2015, G.Cosmo
1031 phi += twopi; 1068 phi += twopi;
1032 1069
1033 if (phi > startPhi+deltaPhi) 1070 if (phi > startPhi+deltaPhi)
1034 { 1071 {
1035 // 1072 //
1036 // Oops. Are we closer to the start phi 1073 // Oops. Are we closer to the start phi or end phi?
1037 // 1074 //
1038 G4double d1 = phi-startPhi-deltaPhi; 1075 G4double d1 = phi-startPhi-deltaPhi;
1039 while( phi > startPhi ) // Loop chec 1076 while( phi > startPhi ) // Loop checking, 13.08.2015, G.Cosmo
1040 phi -= twopi; 1077 phi -= twopi;
1041 G4double d2 = startPhi-phi; 1078 G4double d2 = startPhi-phi;
1042 1079
1043 if (d2 < d1) d1 = d2; 1080 if (d2 < d1) d1 = d2;
1044 1081
1045 // 1082 //
1046 // Add result to our distance 1083 // Add result to our distance
1047 // 1084 //
1048 G4double dist = d1*rx*part; 1085 G4double dist = d1*rx*part;
1049 1086
1050 distOutside2 += dist*dist; 1087 distOutside2 += dist*dist;
1051 if (edgeRZnorm != nullptr) << 1088 if (edgeRZnorm)
1052 { 1089 {
1053 *edgeRZnorm = std::max(std::fabs(*edg 1090 *edgeRZnorm = std::max(std::fabs(*edgeRZnorm),std::fabs(dist));
1054 } 1091 }
1055 } 1092 }
1056 } 1093 }
1057 1094
1058 return answer; 1095 return answer;
1059 } 1096 }
1060 1097
>> 1098
>> 1099 //
1061 // PointOnCone 1100 // PointOnCone
1062 // 1101 //
1063 // Decide if a point is on a cone and return 1102 // Decide if a point is on a cone and return normal if it is
1064 // 1103 //
1065 G4bool G4PolyconeSide::PointOnCone( const G4T << 1104 G4bool G4PolyconeSide::PointOnCone( const G4ThreeVector &hit,
1066 G4d 1105 G4double normSign,
1067 const G4T << 1106 const G4ThreeVector &p,
1068 const G4T << 1107 const G4ThreeVector &v,
1069 G4T << 1108 G4ThreeVector &normal )
1070 { 1109 {
1071 G4double rx = hit.perp(); 1110 G4double rx = hit.perp();
1072 // 1111 //
1073 // Check radial/z extent, as appropriate 1112 // Check radial/z extent, as appropriate
1074 // 1113 //
1075 if (!cone->HitOn( rx, hit.z() )) return fal 1114 if (!cone->HitOn( rx, hit.z() )) return false;
1076 1115
1077 if (phiIsOpen) 1116 if (phiIsOpen)
1078 { 1117 {
1079 G4double phiTolerant = 2.0*kCarTolerance/ 1118 G4double phiTolerant = 2.0*kCarTolerance/(rx+kCarTolerance);
1080 // 1119 //
1081 // Check phi segment. Here we have to be 1120 // Check phi segment. Here we have to be careful
1082 // to use the standard method consistent 1121 // to use the standard method consistent with
1083 // PolyPhiFace. See PolyPhiFace::InsideEd 1122 // PolyPhiFace. See PolyPhiFace::InsideEdgesExact
1084 // 1123 //
1085 G4double phi = GetPhi(hit); 1124 G4double phi = GetPhi(hit);
1086 while( phi < startPhi-phiTolerant ) // 1125 while( phi < startPhi-phiTolerant ) // Loop checking, 13.08.2015, G.Cosmo
1087 phi += twopi; 1126 phi += twopi;
1088 1127
1089 if (phi > startPhi+deltaPhi+phiTolerant) 1128 if (phi > startPhi+deltaPhi+phiTolerant) return false;
1090 1129
1091 if (phi > startPhi+deltaPhi-phiTolerant) 1130 if (phi > startPhi+deltaPhi-phiTolerant)
1092 { 1131 {
1093 // 1132 //
1094 // Exact treatment 1133 // Exact treatment
1095 // 1134 //
1096 G4ThreeVector qx = p + v; 1135 G4ThreeVector qx = p + v;
1097 G4ThreeVector qa = qx - corners[2], 1136 G4ThreeVector qa = qx - corners[2],
1098 qb = qx - corners[3]; 1137 qb = qx - corners[3];
1099 G4ThreeVector qacb = qa.cross(qb); 1138 G4ThreeVector qacb = qa.cross(qb);
1100 1139
1101 if (normSign*qacb.dot(v) < 0) return fa 1140 if (normSign*qacb.dot(v) < 0) return false;
1102 } 1141 }
1103 else if (phi < phiTolerant) 1142 else if (phi < phiTolerant)
1104 { 1143 {
1105 G4ThreeVector qx = p + v; 1144 G4ThreeVector qx = p + v;
1106 G4ThreeVector qa = qx - corners[1], 1145 G4ThreeVector qa = qx - corners[1],
1107 qb = qx - corners[0]; 1146 qb = qx - corners[0];
1108 G4ThreeVector qacb = qa.cross(qb); 1147 G4ThreeVector qacb = qa.cross(qb);
1109 1148
1110 if (normSign*qacb.dot(v) < 0) return fa 1149 if (normSign*qacb.dot(v) < 0) return false;
1111 } 1150 }
1112 } 1151 }
1113 1152
1114 // 1153 //
1115 // We have a good hit! Calculate normal 1154 // We have a good hit! Calculate normal
1116 // 1155 //
1117 if (rx < DBL_MIN) 1156 if (rx < DBL_MIN)
1118 normal = G4ThreeVector( 0, 0, zNorm < 0 ? 1157 normal = G4ThreeVector( 0, 0, zNorm < 0 ? -1 : 1 );
1119 else 1158 else
1120 normal = G4ThreeVector( rNorm*hit.x()/rx, 1159 normal = G4ThreeVector( rNorm*hit.x()/rx, rNorm*hit.y()/rx, zNorm );
1121 return true; 1160 return true;
1122 } 1161 }
1123 1162
>> 1163
>> 1164 //
1124 // FindLineIntersect 1165 // FindLineIntersect
1125 // 1166 //
1126 // Decide the point at which two 2-dimensiona 1167 // Decide the point at which two 2-dimensional lines intersect
1127 // 1168 //
1128 // Equation of line: x = x1 + s*tx1 1169 // Equation of line: x = x1 + s*tx1
1129 // y = y1 + s*ty1 1170 // y = y1 + s*ty1
1130 // 1171 //
1131 // It is assumed that the lines are *not* par 1172 // It is assumed that the lines are *not* parallel
1132 // 1173 //
1133 void G4PolyconeSide::FindLineIntersect( G4dou 1174 void G4PolyconeSide::FindLineIntersect( G4double x1, G4double y1,
1134 G4dou 1175 G4double tx1, G4double ty1,
1135 G4dou 1176 G4double x2, G4double y2,
1136 G4dou 1177 G4double tx2, G4double ty2,
1137 G4dou << 1178 G4double &x, G4double &y )
1138 { 1179 {
1139 // 1180 //
1140 // The solution is a simple linear equation 1181 // The solution is a simple linear equation
1141 // 1182 //
1142 G4double deter = tx1*ty2 - tx2*ty1; 1183 G4double deter = tx1*ty2 - tx2*ty1;
1143 1184
1144 G4double s1 = ((x2-x1)*ty2 - tx2*(y2-y1))/d 1185 G4double s1 = ((x2-x1)*ty2 - tx2*(y2-y1))/deter;
1145 G4double s2 = ((x2-x1)*ty1 - tx1*(y2-y1))/d 1186 G4double s2 = ((x2-x1)*ty1 - tx1*(y2-y1))/deter;
1146 1187
1147 // 1188 //
1148 // We want the answer to not depend on whic 1189 // We want the answer to not depend on which order the
1149 // lines were specified. Take average. 1190 // lines were specified. Take average.
1150 // 1191 //
1151 x = 0.5*( x1+s1*tx1 + x2+s2*tx2 ); 1192 x = 0.5*( x1+s1*tx1 + x2+s2*tx2 );
1152 y = 0.5*( y1+s1*ty1 + y2+s2*ty2 ); 1193 y = 0.5*( y1+s1*ty1 + y2+s2*ty2 );
1153 } 1194 }
1154 1195
>> 1196
>> 1197 //
1155 // Calculate surface area for GetPointOnSurfa 1198 // Calculate surface area for GetPointOnSurface()
1156 // 1199 //
1157 G4double G4PolyconeSide::SurfaceArea() 1200 G4double G4PolyconeSide::SurfaceArea()
1158 { 1201 {
1159 if(fSurfaceArea==0.) << 1202 if(fSurfaceArea==0)
1160 { 1203 {
1161 fSurfaceArea = (r[0]+r[1])* std::sqrt(sqr 1204 fSurfaceArea = (r[0]+r[1])* std::sqrt(sqr(r[0]-r[1])+sqr(z[0]-z[1]));
1162 fSurfaceArea *= 0.5*(deltaPhi); 1205 fSurfaceArea *= 0.5*(deltaPhi);
1163 } 1206 }
1164 return fSurfaceArea; 1207 return fSurfaceArea;
1165 } 1208 }
1166 1209
>> 1210
>> 1211 //
1167 // GetPointOnFace 1212 // GetPointOnFace
1168 // 1213 //
1169 G4ThreeVector G4PolyconeSide::GetPointOnFace( 1214 G4ThreeVector G4PolyconeSide::GetPointOnFace()
1170 { 1215 {
1171 G4double x,y,zz; 1216 G4double x,y,zz;
1172 G4double rr,phi,dz,dr; 1217 G4double rr,phi,dz,dr;
1173 dr=r[1]-r[0];dz=z[1]-z[0]; 1218 dr=r[1]-r[0];dz=z[1]-z[0];
1174 phi=startPhi+deltaPhi*G4UniformRand(); 1219 phi=startPhi+deltaPhi*G4UniformRand();
1175 rr=r[0]+dr*G4UniformRand(); 1220 rr=r[0]+dr*G4UniformRand();
1176 1221
1177 x=rr*std::cos(phi); 1222 x=rr*std::cos(phi);
1178 y=rr*std::sin(phi); 1223 y=rr*std::sin(phi);
1179 1224
1180 // PolyconeSide has a Ring Form 1225 // PolyconeSide has a Ring Form
1181 // 1226 //
1182 if (dz==0.) 1227 if (dz==0.)
1183 { 1228 {
1184 zz=z[0]; 1229 zz=z[0];
1185 } 1230 }
1186 else 1231 else
1187 { 1232 {
1188 if(dr==0.) // PolyconeSide has a Tube Fo 1233 if(dr==0.) // PolyconeSide has a Tube Form
1189 { 1234 {
1190 zz = z[0]+dz*G4UniformRand(); 1235 zz = z[0]+dz*G4UniformRand();
1191 } 1236 }
1192 else 1237 else
1193 { 1238 {
1194 zz = z[0]+(rr-r[0])*dz/dr; 1239 zz = z[0]+(rr-r[0])*dz/dr;
1195 } 1240 }
1196 } 1241 }
1197 1242
1198 return {x,y,zz}; << 1243 return G4ThreeVector(x,y,zz);
1199 } 1244 }
1200 1245