1 | // |
---|
2 | // ******************************************************************** |
---|
3 | // * License and Disclaimer * |
---|
4 | // * * |
---|
5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
7 | // * conditions of the Geant4 Software License, included in the file * |
---|
8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
9 | // * include a list of copyright holders. * |
---|
10 | // * * |
---|
11 | // * Neither the authors of this software system, nor their employing * |
---|
12 | // * institutes,nor the agencies providing financial support for this * |
---|
13 | // * work make any representation or warranty, express or implied, * |
---|
14 | // * regarding this software system or assume any liability for its * |
---|
15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
16 | // * for the full disclaimer and the limitation of liability. * |
---|
17 | // * * |
---|
18 | // * This code implementation is the result of the scientific and * |
---|
19 | // * technical work of the GEANT4 collaboration. * |
---|
20 | // * By using, copying, modifying or distributing the software (or * |
---|
21 | // * any work based on the software) you agree to acknowledge its * |
---|
22 | // * use in resulting scientific publications, and indicate your * |
---|
23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
24 | // ******************************************************************** |
---|
25 | // |
---|
26 | // $Id: G4BREPSolidPolyhedra.cc,v 1.35 2008/01/22 16:04:58 tnikitin Exp $ |
---|
27 | // GEANT4 tag $Name: HEAD $ |
---|
28 | // |
---|
29 | // ---------------------------------------------------------------------- |
---|
30 | // GEANT 4 class source file |
---|
31 | // |
---|
32 | // G4BREPSolidPolyhedra.cc |
---|
33 | // |
---|
34 | // ---------------------------------------------------------------------- |
---|
35 | // The polygonal solid G4BREPSolidPolyhedra is a shape defined by an inner |
---|
36 | // and outer polygonal surface and two planes perpendicular to the Z axis. |
---|
37 | // Each polygonal surface is created by linking a series of polygons created |
---|
38 | // at different planes perpendicular to the Z-axis. All these polygons all |
---|
39 | // have the same number of sides (sides) and are defined at the same Z planes |
---|
40 | // for both inner and outer polygonal surfaces. |
---|
41 | // ---------------------------------------------------------------------- |
---|
42 | // |
---|
43 | // History |
---|
44 | // ------- |
---|
45 | // |
---|
46 | // Bug-fix #266 by R.Chytracek: |
---|
47 | // - The situation when phi1 = 0 dphi1 = 2*pi and all RMINs = 0.0 is handled |
---|
48 | // now. In this case the inner planes are not created. The fix goes even |
---|
49 | // further this means it considers more than 2 z-planes and inner planes |
---|
50 | // are not created whenever two consecutive RMINs are = 0.0 . |
---|
51 | // |
---|
52 | // Corrections by S.Giani: |
---|
53 | // - Xaxis now corresponds to phi=0 |
---|
54 | // - partial angle = phiTotal / Nsides |
---|
55 | // - end planes exact boundary calculation for phiTotal < 2pi |
---|
56 | // (also including case with RMIN=RMAX) |
---|
57 | // - Xaxis now properly rotated to compute correct scope of vertixes |
---|
58 | // - corrected surface orientation for outer faces parallel to Z |
---|
59 | // - completed explicit setting of the orientation for all faces |
---|
60 | // - some comparison between doubles avoided by using tolerances |
---|
61 | // - visualisation parameters made consistent with the use made by |
---|
62 | // constructor of the input arguments (i.e. circumscribed radius). |
---|
63 | // ---------------------------------------------------------------------- |
---|
64 | |
---|
65 | #include "G4BREPSolidPolyhedra.hh" |
---|
66 | #include "G4FPlane.hh" |
---|
67 | |
---|
68 | #include <sstream> |
---|
69 | |
---|
70 | G4BREPSolidPolyhedra::G4BREPSolidPolyhedra(const G4String& name, |
---|
71 | G4double start_angle, |
---|
72 | G4double opening_angle, |
---|
73 | G4int sides, |
---|
74 | G4int num_z_planes, |
---|
75 | G4double z_start, |
---|
76 | G4double z_values[], |
---|
77 | G4double RMIN[], |
---|
78 | G4double RMAX[] ) |
---|
79 | : G4BREPSolid(name) |
---|
80 | { |
---|
81 | G4int sections = num_z_planes - 1; |
---|
82 | |
---|
83 | if( opening_angle >= 2*pi-perMillion ) |
---|
84 | { |
---|
85 | nb_of_surfaces = 2*(sections * sides) + 2; |
---|
86 | } |
---|
87 | else |
---|
88 | { |
---|
89 | nb_of_surfaces = 2*(sections * sides) + 4; |
---|
90 | } |
---|
91 | |
---|
92 | //SurfaceVec = new G4Surface*[nb_of_surfaces]; |
---|
93 | G4int MaxNbOfSurfaces = nb_of_surfaces; |
---|
94 | G4Surface** MaxSurfaceVec = new G4Surface*[MaxNbOfSurfaces]; |
---|
95 | |
---|
96 | G4Vector3D Axis(0,0,1); |
---|
97 | G4Vector3D XAxis(1,0,0); |
---|
98 | G4Vector3D TmpAxis; |
---|
99 | G4Point3D Origin(0,0,z_start); |
---|
100 | G4Point3D LocalOrigin(0,0,z_start); |
---|
101 | G4double Length; |
---|
102 | G4int Count = 0 ; |
---|
103 | G4double PartAngle = (opening_angle)/sides; |
---|
104 | |
---|
105 | |
---|
106 | /////////////////////////////////////////////////// |
---|
107 | // Preconditions check |
---|
108 | |
---|
109 | // Detecting minimal required number of sides |
---|
110 | // |
---|
111 | if( sides < 3 ) |
---|
112 | { |
---|
113 | G4Exception( "G4BREPSolidPolyhedra::G4BREPSolidPolyhedra()", |
---|
114 | "InvalidSetup", FatalException, |
---|
115 | "The solid must have at least 3 sides!" ); |
---|
116 | } |
---|
117 | |
---|
118 | // Detecting minimal required number of z-sections |
---|
119 | // |
---|
120 | if( num_z_planes < 2 ) |
---|
121 | { |
---|
122 | G4Exception( "G4BREPSolidPolyhedra::G4BREPSolidPolyhedra()", |
---|
123 | "InvalidSetup", FatalException, |
---|
124 | "The solid must have at least 2 z-sections!" ); |
---|
125 | } |
---|
126 | |
---|
127 | // Detect invalid configurations at the ends of polyhedra which |
---|
128 | // would not lead to a valid solid creation and likely to a crash |
---|
129 | // |
---|
130 | if( z_values[0] == z_values[1] |
---|
131 | || z_values[sections-1] == z_values[sections] ) |
---|
132 | { |
---|
133 | G4Exception( "G4BREPSolidPolyhedra::G4BREPSolidPolyhedra()", |
---|
134 | "InvalidSetup", FatalException, |
---|
135 | "The solid must have the first 2 and the last 2 z-values different!" ); |
---|
136 | } |
---|
137 | |
---|
138 | // Find out how the z-values sequence is ordered |
---|
139 | // |
---|
140 | G4bool increasing; |
---|
141 | if( z_values[0] < z_values[1] ) |
---|
142 | { |
---|
143 | increasing = true; |
---|
144 | } |
---|
145 | else |
---|
146 | { |
---|
147 | increasing = false; |
---|
148 | } |
---|
149 | |
---|
150 | // Detecting polyhedra teeth. |
---|
151 | // It's forbidden to specify unordered, e.g. non-increasing or |
---|
152 | // non-decreasing sequence of z-values. It may be provided by a |
---|
153 | // specific solid in a future. |
---|
154 | // |
---|
155 | for( G4int idx = 0; idx < sections; idx++ ) |
---|
156 | { |
---|
157 | if( ( z_values[idx] > z_values[idx+1] && increasing ) || |
---|
158 | ( z_values[idx] < z_values[idx+1] && !increasing ) ) |
---|
159 | { |
---|
160 | // ERROR! Invalid sequence of z-values |
---|
161 | // |
---|
162 | std::ostringstream msgstr; |
---|
163 | msgstr << G4endl |
---|
164 | << "ERROR: unordered, non-increasing or non-decreasing sequence" |
---|
165 | << G4endl |
---|
166 | << " of z_values detected!" |
---|
167 | << G4endl |
---|
168 | << " Check z_values with indexes: " |
---|
169 | << idx << " " << (idx+1) << "." << G4endl << std::ends; |
---|
170 | G4String message = msgstr.str(); |
---|
171 | G4Exception( "G4BREPSolidPolyhedra::G4BREPSolidPolyhedra()", |
---|
172 | "InvalidSetup", FatalException, message ); |
---|
173 | } |
---|
174 | } |
---|
175 | |
---|
176 | /////////////////////////////////////////////////// |
---|
177 | #ifdef G4_EXPERIMENTAL_CODE |
---|
178 | // There is one problem when sequence of z values is not increasing in a |
---|
179 | // regular way, in other words, it's not purely increasing or decreasing |
---|
180 | // Irregular sequence can be provided in order to define a polyhedra having |
---|
181 | // teeth as shown on the picture bellow |
---|
182 | // In this sequence can happen the following: |
---|
183 | // z[a-1] > z[a] < z[a+1] && z[a+1] >= z[a-1] |
---|
184 | // One has to check the RMAX and RMIN values due to the possible |
---|
185 | // intersections. |
---|
186 | // |
---|
187 | // 1 2 3 |
---|
188 | // ___ ___ ____ |
---|
189 | // 00/ 00/ _ 000/ |
---|
190 | // 0/ 0/ |0 00| |
---|
191 | // V___ V__+0 00+-- |
---|
192 | // 0000 00000 00000 |
---|
193 | // ---- ----- ----- |
---|
194 | // ------------------------------------ z-axis |
---|
195 | // |
---|
196 | // |
---|
197 | // NOTE: This picture doesn't show all the possible configurations of |
---|
198 | // a polyhedra having teeth when looking at its profile. |
---|
199 | // The picture shows only one half of the polyhedra's profile |
---|
200 | ///////////////////////////////////////////////////////////////////////// |
---|
201 | |
---|
202 | // Experimental code! Not recommended for production, it's incomplete! |
---|
203 | // The task is to identify invalid combination of z, RMIN and RMAX values |
---|
204 | // in the case of toothydra :-) |
---|
205 | // |
---|
206 | G4int toothIdx; |
---|
207 | |
---|
208 | for( G4int idx = 1; idx < sections+1; idx++ ) |
---|
209 | { |
---|
210 | if( z_values[idx-1] > z_values[idx] ) |
---|
211 | { |
---|
212 | G4double toothdist = std::fabs( z_values[idx-1] - z_values[idx] ); |
---|
213 | G4double aftertoothdist = std::fabs( z_values[idx+1] - z_values[idx] ); |
---|
214 | if( toothdist > aftertoothdist ) |
---|
215 | { |
---|
216 | // Check for possible intersection |
---|
217 | // |
---|
218 | if( RMAX[idx-1] < RMAX[idx+1] || RMIN[idx-1] > RMIN[idx+1] ) |
---|
219 | { |
---|
220 | // ERROR! The surface conflict! |
---|
221 | // |
---|
222 | std::ostringstream msgstr; |
---|
223 | msgstr << G4endl |
---|
224 | << "ERROR: unordered sequence of z_values detected with" |
---|
225 | << G4endl |
---|
226 | << " conflicting RMAX or RMIN values!" |
---|
227 | << G4endl |
---|
228 | << " Check z_values with indexes: " |
---|
229 | << (idx-1) << " " << idx << " " << (idx+1) << "." |
---|
230 | << G4endl << std::ends; |
---|
231 | G4Exception( "G4BREPSolidPolyhedra::G4BREPSolidPolyhedra()", |
---|
232 | "InvalidSetup", FatalException, msgstr.str() ); |
---|
233 | } |
---|
234 | } |
---|
235 | } |
---|
236 | } |
---|
237 | #endif // G4_EXPERIMENTAL_CODE |
---|
238 | /////////////////////////////////////////////////// |
---|
239 | |
---|
240 | for(G4int a=0;a<sections;a++) |
---|
241 | { |
---|
242 | Length = z_values[a+1] - z_values[a]; |
---|
243 | |
---|
244 | if( Length != 0.0 ) |
---|
245 | { |
---|
246 | TmpAxis= XAxis; |
---|
247 | TmpAxis.rotateZ(start_angle); |
---|
248 | |
---|
249 | // L. Broglia: Be careful in the construction of the planes, |
---|
250 | // see G4FPlane |
---|
251 | // |
---|
252 | for( G4int b = 0; b < sides; b++ ) |
---|
253 | { |
---|
254 | // Create inner side by calculation of points for the planar surface |
---|
255 | // boundary. The order of the points gives the surface sense -> changed |
---|
256 | // to explicit sense set-up by R. Chytracek, 12/02/2002 |
---|
257 | // We must check if a pair of two consecutive RMINs is not = 0.0, |
---|
258 | // this means no inner plane exists! |
---|
259 | // |
---|
260 | if( RMIN[a] != 0.0 ) |
---|
261 | { |
---|
262 | if( RMIN[a+1] != 0.0 ) |
---|
263 | { |
---|
264 | // Standard case |
---|
265 | // |
---|
266 | MaxSurfaceVec[Count] = |
---|
267 | CreateTrapezoidalSurface( RMIN[a], RMIN[a+1], |
---|
268 | LocalOrigin, Length, |
---|
269 | TmpAxis, PartAngle, EInverse ); |
---|
270 | } |
---|
271 | else |
---|
272 | { |
---|
273 | // The special case of r1 > r2 where we end at the |
---|
274 | // point (0,0,z[a+1]) |
---|
275 | // |
---|
276 | MaxSurfaceVec[Count] = |
---|
277 | CreateTriangularSurface( RMIN[a], RMIN[a+1], |
---|
278 | LocalOrigin, Length, |
---|
279 | TmpAxis, PartAngle, EInverse ); |
---|
280 | } |
---|
281 | } |
---|
282 | else if( RMIN[a+1] != 0.0 ) |
---|
283 | { |
---|
284 | // The special case of r1 < r2 where we start at the |
---|
285 | // point ( 0,0,z[a]) |
---|
286 | // |
---|
287 | MaxSurfaceVec[Count] = |
---|
288 | CreateTriangularSurface( RMIN[a], RMIN[a+1], LocalOrigin, Length, |
---|
289 | TmpAxis, PartAngle, EInverse ); |
---|
290 | } |
---|
291 | else |
---|
292 | { |
---|
293 | // Insert nothing into the vector of sufaces, we'll replicate |
---|
294 | // the vector anyway later |
---|
295 | // |
---|
296 | MaxSurfaceVec[Count] = 0; |
---|
297 | |
---|
298 | // We need to reduce the number of planes by 1, |
---|
299 | // one we have just skipped |
---|
300 | // |
---|
301 | nb_of_surfaces--; |
---|
302 | } |
---|
303 | |
---|
304 | if( MaxSurfaceVec[Count] != 0 ) |
---|
305 | { |
---|
306 | // Rotate axis back for the other surface point calculation |
---|
307 | // only in the case any of the Create* methods above have been |
---|
308 | // called because they modify the passed in TmpAxis |
---|
309 | // |
---|
310 | TmpAxis.rotateZ(-PartAngle); |
---|
311 | } |
---|
312 | |
---|
313 | Count++; |
---|
314 | |
---|
315 | // Create outer side |
---|
316 | |
---|
317 | if( RMAX[a] != 0.0 ) |
---|
318 | { |
---|
319 | if( RMAX[a+1] != 0.0 ) |
---|
320 | { |
---|
321 | // Standard case |
---|
322 | // |
---|
323 | MaxSurfaceVec[Count] = |
---|
324 | CreateTrapezoidalSurface( RMAX[a], RMAX[a+1], |
---|
325 | LocalOrigin, Length, |
---|
326 | TmpAxis, PartAngle, ENormal ); |
---|
327 | } |
---|
328 | else |
---|
329 | { |
---|
330 | // The special case of r1 > r2 where we end |
---|
331 | // at the point (0,0,z[a+1]) |
---|
332 | // |
---|
333 | MaxSurfaceVec[Count] = |
---|
334 | CreateTriangularSurface( RMAX[a], RMAX[a+1], |
---|
335 | LocalOrigin, Length, |
---|
336 | TmpAxis, PartAngle, ENormal ); |
---|
337 | } |
---|
338 | } |
---|
339 | else if( RMAX[a+1] != 0.0 ) |
---|
340 | { |
---|
341 | // The special case of r1 < r2 where we start |
---|
342 | // at the point ( 0,0,z[a]) |
---|
343 | // |
---|
344 | MaxSurfaceVec[Count] = |
---|
345 | CreateTriangularSurface( RMAX[a], RMAX[a+1], LocalOrigin, Length, |
---|
346 | TmpAxis, PartAngle, ENormal ); |
---|
347 | } |
---|
348 | else |
---|
349 | { |
---|
350 | // Two consecutive RMAX values can't be zero as |
---|
351 | // it's against the definition of BREP polyhedra |
---|
352 | // |
---|
353 | G4Exception( "G4BREPSolidPolyhedra::G4BREPSolidPolyhedra()", |
---|
354 | "InvalidSetup", FatalException, |
---|
355 | "Two consecutive RMAX values cannot be zero!" ); |
---|
356 | } |
---|
357 | |
---|
358 | Count++; |
---|
359 | } // End of for loop over sides |
---|
360 | } |
---|
361 | else |
---|
362 | { |
---|
363 | // Create planar surfaces perpendicular to z-axis |
---|
364 | |
---|
365 | ESurfaceSense OuterSurfSense, InnerSurfSense; |
---|
366 | |
---|
367 | if( RMAX[a] != RMAX[a+1] && RMIN[a] != RMIN[a+1] ) |
---|
368 | { |
---|
369 | // We're about to create a planar surface perpendicular to z-axis |
---|
370 | // We can have the 8 following configurations here: |
---|
371 | // |
---|
372 | // 1. 2. 3. 4. |
---|
373 | // --+ +-- --+ +-- |
---|
374 | // xx|-> <-|xx xx| |xx |
---|
375 | // xx+-- --+xx --+ +-- |
---|
376 | // xxxxx xxxxx | | |
---|
377 | // xxxxx xxxxx +-- --+ |
---|
378 | // xx+-- --+xx |xx xx| |
---|
379 | // xx|-> <-|xx +-- --+ |
---|
380 | // --+ +-- |
---|
381 | // -------------------------- Z axis |
---|
382 | // |
---|
383 | ////////////////////////////////////////////////////////////// |
---|
384 | ////////////////////////////////////////////////////////////// |
---|
385 | // |
---|
386 | // 5. 6. 7. 8. |
---|
387 | // --+ +-- --+ +-- |
---|
388 | // xx|-> <-|xx xx|-> <-|xx |
---|
389 | // --+-- --+-- xx+-- --+xx |
---|
390 | // <-|xx xx|-> xxxxx xxxxx |
---|
391 | // +-- --+ --+xx xx+-- |
---|
392 | // <-|xx xx|-> |
---|
393 | // +-- --+ |
---|
394 | // -------------------------- Z axis |
---|
395 | // |
---|
396 | // NOTE: The pictures shows only one half of polyhedra! |
---|
397 | // The arrows show the expected surface normal direction. |
---|
398 | // The configuration No. 3 and 4 are not valid solids! |
---|
399 | |
---|
400 | // Eliminate the invalid cases 3 and 4. |
---|
401 | // At this point is guaranteed that each RMIN[i] < RMAX[i] |
---|
402 | // where i in in interval 0 < i < num_z_planes-1. So: |
---|
403 | // |
---|
404 | if( RMIN[a] > RMAX[a+1] || RMAX[a] < RMIN[a+1] ) |
---|
405 | { |
---|
406 | std::stringstream s; |
---|
407 | s << G4endl << "The values of RMIN[" << a << "] & RMAX[" << a+1 |
---|
408 | << "] or RMAX[" << a << "] & RMIN[" << a+1 << "] " |
---|
409 | << "generate an invalid configuration of solid: " |
---|
410 | << name.c_str() << "!" << G4endl << std::ends; |
---|
411 | G4String message = s.str(); |
---|
412 | G4Exception( "G4BREPSolidPolyhedra::G4BREPSolidPolyhedra()", |
---|
413 | "InvalidSetup", FatalException, message ); |
---|
414 | } |
---|
415 | |
---|
416 | // We need to clasify all the cases in order to figure out |
---|
417 | // the planar surface sense |
---|
418 | // |
---|
419 | if( RMAX[a] > RMAX[a+1] ) |
---|
420 | { |
---|
421 | // Cases 1, 5, 7 |
---|
422 | // |
---|
423 | if( RMIN[a] < RMIN[a+1] ) |
---|
424 | { |
---|
425 | // Case 1 |
---|
426 | OuterSurfSense = EInverse; |
---|
427 | InnerSurfSense = EInverse; |
---|
428 | } |
---|
429 | else if( RMAX[a+1] != RMIN[a]) |
---|
430 | { |
---|
431 | // Case 7 |
---|
432 | OuterSurfSense = EInverse; |
---|
433 | InnerSurfSense = ENormal; |
---|
434 | } |
---|
435 | else |
---|
436 | { |
---|
437 | // Case 5 |
---|
438 | OuterSurfSense = EInverse; |
---|
439 | InnerSurfSense = ENormal; |
---|
440 | } |
---|
441 | } |
---|
442 | else |
---|
443 | { |
---|
444 | // Cases 2, 6, 8 |
---|
445 | if( RMIN[a] > RMIN[a+1] ) |
---|
446 | { |
---|
447 | // Case 2 |
---|
448 | OuterSurfSense = ENormal; |
---|
449 | InnerSurfSense = ENormal; |
---|
450 | } |
---|
451 | else if( RMIN[a+1] != RMAX[a] ) |
---|
452 | { |
---|
453 | // Case 8 |
---|
454 | OuterSurfSense = ENormal; |
---|
455 | InnerSurfSense = EInverse; |
---|
456 | } |
---|
457 | else |
---|
458 | { |
---|
459 | // Case 6 |
---|
460 | OuterSurfSense = ENormal; |
---|
461 | InnerSurfSense = EInverse; |
---|
462 | } |
---|
463 | } |
---|
464 | |
---|
465 | TmpAxis= XAxis; |
---|
466 | TmpAxis.rotateZ(start_angle); |
---|
467 | |
---|
468 | // Compute the outer planar surface |
---|
469 | // |
---|
470 | MaxSurfaceVec[Count] = |
---|
471 | ComputePlanarSurface( RMAX[a], RMAX[a+1], LocalOrigin, TmpAxis, |
---|
472 | sides, PartAngle, OuterSurfSense ); |
---|
473 | if( MaxSurfaceVec[Count] == 0 ) |
---|
474 | { |
---|
475 | // No surface was created |
---|
476 | // |
---|
477 | nb_of_surfaces--; |
---|
478 | } |
---|
479 | Count++; |
---|
480 | |
---|
481 | TmpAxis= XAxis; |
---|
482 | TmpAxis.rotateZ(start_angle); |
---|
483 | |
---|
484 | // Compute the inner planar surface |
---|
485 | // |
---|
486 | MaxSurfaceVec[Count] = |
---|
487 | ComputePlanarSurface( RMIN[a], RMIN[a+1], LocalOrigin, TmpAxis, |
---|
488 | sides, PartAngle, InnerSurfSense ); |
---|
489 | if( MaxSurfaceVec[Count] == 0 ) |
---|
490 | { |
---|
491 | // No surface was created |
---|
492 | // |
---|
493 | nb_of_surfaces--; |
---|
494 | } |
---|
495 | Count++; |
---|
496 | |
---|
497 | // Since we can create here at maximum 2 surfaces |
---|
498 | // we need to reflect this in the total |
---|
499 | // |
---|
500 | nb_of_surfaces -= (2*(sides-1)); |
---|
501 | } |
---|
502 | else |
---|
503 | { |
---|
504 | // The case where only one of the radius values has changed |
---|
505 | // |
---|
506 | // RMAX RMIN |
---|
507 | // change change |
---|
508 | // |
---|
509 | // 1 2 3 4 |
---|
510 | // --+ +-- ----- ----- |
---|
511 | // 00|-> <-|00 00000 00000 |
---|
512 | // 00+-- --+00 --+00 00+-- |
---|
513 | // 00000 00000 <-|00 00|-> |
---|
514 | // +-- --+ |
---|
515 | // --------------------------- Z axis |
---|
516 | // |
---|
517 | // NOTE: The picture shows only one half of polyhedra! |
---|
518 | |
---|
519 | G4double R1, R2; |
---|
520 | ESurfaceSense SurfSense; |
---|
521 | |
---|
522 | // The case by case classification |
---|
523 | // |
---|
524 | if( RMAX[a] != RMAX[a+1] ) |
---|
525 | { |
---|
526 | // Cases 1, 2 |
---|
527 | // |
---|
528 | R1 = RMAX[a]; |
---|
529 | R2 = RMAX[a+1]; |
---|
530 | if( R1 > R2 ) |
---|
531 | { |
---|
532 | // Case 1 |
---|
533 | // |
---|
534 | SurfSense = EInverse; |
---|
535 | } |
---|
536 | else |
---|
537 | { |
---|
538 | // Case 2 |
---|
539 | // |
---|
540 | SurfSense = ENormal; |
---|
541 | } |
---|
542 | } |
---|
543 | else if(RMIN[a] != RMIN[a+1]) |
---|
544 | { |
---|
545 | // Cases 3, 4 |
---|
546 | // |
---|
547 | R1 = RMIN[a]; |
---|
548 | R2 = RMIN[a+1]; |
---|
549 | if( R1 > R2 ) |
---|
550 | { |
---|
551 | // Case 3 |
---|
552 | // |
---|
553 | SurfSense = ENormal; |
---|
554 | } |
---|
555 | else |
---|
556 | { |
---|
557 | // Case 4 |
---|
558 | // |
---|
559 | SurfSense = EInverse; |
---|
560 | } |
---|
561 | } |
---|
562 | else |
---|
563 | { |
---|
564 | G4cerr << "ERROR - G4BREPSolidPolyhedra::G4BREPSolidPolyhedra()" |
---|
565 | << G4endl |
---|
566 | << " Error in construction." |
---|
567 | << G4endl |
---|
568 | << " Exactly the same z, rmin and rmax given for" |
---|
569 | << G4endl |
---|
570 | << " consecutive indices, " << a << " and " << a+1 |
---|
571 | << G4endl; |
---|
572 | G4Exception( "G4BREPSolidPolyhedra::G4BREPSolidPolyhedra()", |
---|
573 | "Notification", JustWarning, "Construction Error!" ); |
---|
574 | continue; |
---|
575 | } |
---|
576 | TmpAxis= XAxis; |
---|
577 | TmpAxis.rotateZ(start_angle); |
---|
578 | |
---|
579 | MaxSurfaceVec[Count] = |
---|
580 | ComputePlanarSurface( R1, R2, LocalOrigin, TmpAxis, |
---|
581 | sides, PartAngle, SurfSense ); |
---|
582 | if( MaxSurfaceVec[Count] == 0 ) |
---|
583 | { |
---|
584 | // No surface was created |
---|
585 | // |
---|
586 | nb_of_surfaces--; |
---|
587 | } |
---|
588 | Count++; |
---|
589 | |
---|
590 | // Since we can create here at maximum 1 surface |
---|
591 | // we need to reflect this in the total |
---|
592 | // |
---|
593 | nb_of_surfaces -= ((2*sides) - 1); |
---|
594 | } |
---|
595 | } // End of if( Length != 0.0 ) |
---|
596 | |
---|
597 | LocalOrigin = LocalOrigin + (Length*Axis); |
---|
598 | |
---|
599 | } // End of for loop over z sections |
---|
600 | |
---|
601 | if(opening_angle >= 2*pi-perMillion) |
---|
602 | { |
---|
603 | // Create the end planes for the configuration where delta phi >= 2*PI |
---|
604 | |
---|
605 | TmpAxis = XAxis; |
---|
606 | TmpAxis.rotateZ(start_angle); |
---|
607 | |
---|
608 | MaxSurfaceVec[Count] = |
---|
609 | ComputePlanarSurface( RMIN[0], RMAX[0], Origin, TmpAxis, |
---|
610 | sides, PartAngle, ENormal ); |
---|
611 | |
---|
612 | if( MaxSurfaceVec[Count] == 0 ) |
---|
613 | { |
---|
614 | // No surface was created |
---|
615 | // |
---|
616 | nb_of_surfaces--; |
---|
617 | } |
---|
618 | Count++; |
---|
619 | |
---|
620 | // Reset plane axis |
---|
621 | // |
---|
622 | TmpAxis = XAxis; |
---|
623 | TmpAxis.rotateZ(start_angle); |
---|
624 | |
---|
625 | MaxSurfaceVec[Count] = |
---|
626 | ComputePlanarSurface( RMIN[sections], RMAX[sections], |
---|
627 | LocalOrigin, TmpAxis, |
---|
628 | sides, PartAngle, EInverse ); |
---|
629 | |
---|
630 | if( MaxSurfaceVec[Count] == 0 ) |
---|
631 | { |
---|
632 | // No surface was created |
---|
633 | // |
---|
634 | nb_of_surfaces--; |
---|
635 | } |
---|
636 | Count++; |
---|
637 | } |
---|
638 | else |
---|
639 | { |
---|
640 | // If delta phi < 2*PI then create a single boundary |
---|
641 | // (case with RMIN=0 included) |
---|
642 | |
---|
643 | // Create the lateral planars |
---|
644 | // |
---|
645 | TmpAxis = XAxis; |
---|
646 | G4Vector3D TmpAxis2 = XAxis; |
---|
647 | TmpAxis.rotateZ(start_angle); |
---|
648 | TmpAxis2.rotateZ(start_angle); |
---|
649 | TmpAxis2.rotateZ(start_angle); |
---|
650 | |
---|
651 | LocalOrigin = Origin; |
---|
652 | G4int points = sections*2+2; |
---|
653 | G4int PointCount = 0; |
---|
654 | |
---|
655 | G4Point3DVector GapPointList(points); |
---|
656 | G4Point3DVector GapPointList2(points); |
---|
657 | |
---|
658 | for(G4int d=0;d<sections+1;d++) |
---|
659 | { |
---|
660 | GapPointList[PointCount] = LocalOrigin + (RMAX[d]*TmpAxis); |
---|
661 | GapPointList[points-1-PointCount] = LocalOrigin + (RMIN[d]*TmpAxis); |
---|
662 | |
---|
663 | GapPointList2[PointCount] = LocalOrigin + (RMAX[d]*TmpAxis2); |
---|
664 | GapPointList2[points-1-PointCount] = LocalOrigin + (RMIN[d]*TmpAxis2); |
---|
665 | |
---|
666 | PointCount++; |
---|
667 | |
---|
668 | Length = z_values[d+1] - z_values[d]; |
---|
669 | LocalOrigin = LocalOrigin+(Length*Axis); |
---|
670 | } |
---|
671 | |
---|
672 | // Add the lateral planars to the surfaces list and set/reverse sense |
---|
673 | // |
---|
674 | MaxSurfaceVec[Count++] = new G4FPlane( &GapPointList, 0, ENormal ); |
---|
675 | MaxSurfaceVec[Count++] = new G4FPlane( &GapPointList2, 0, EInverse ); |
---|
676 | |
---|
677 | TmpAxis = XAxis; |
---|
678 | TmpAxis.rotateZ(start_angle); |
---|
679 | TmpAxis.rotateZ(opening_angle); |
---|
680 | |
---|
681 | // Create end planes |
---|
682 | // |
---|
683 | G4Point3DVector EndPointList ((sides+1)*2); |
---|
684 | G4Point3DVector EndPointList2((sides+1)*2); |
---|
685 | |
---|
686 | for(G4int c=0;c<sides+1;c++) |
---|
687 | { |
---|
688 | // outer polylines for origin end and opposite side |
---|
689 | EndPointList[c] = Origin + (RMAX[0] * TmpAxis); |
---|
690 | EndPointList[(sides+1)*2-1-c] = Origin + (RMIN[0] * TmpAxis); |
---|
691 | EndPointList2[c] = LocalOrigin + (RMAX[sections] * TmpAxis); |
---|
692 | EndPointList2[(sides+1)*2-1-c] = LocalOrigin + (RMIN[sections] * TmpAxis); |
---|
693 | TmpAxis.rotateZ(-PartAngle); |
---|
694 | } |
---|
695 | |
---|
696 | // Add the end planes to the surfaces list |
---|
697 | // Note the surface sense in this case is reversed |
---|
698 | // It's because here we have created the end planes in reversed order |
---|
699 | // than it's done by ComputePlanarSurface() method |
---|
700 | // |
---|
701 | if(RMAX[0]-RMIN[0] >= perMillion) |
---|
702 | { |
---|
703 | MaxSurfaceVec[Count] = new G4FPlane( &EndPointList, 0, EInverse ); |
---|
704 | } |
---|
705 | else |
---|
706 | { |
---|
707 | MaxSurfaceVec[Count] = 0; |
---|
708 | nb_of_surfaces--; |
---|
709 | } |
---|
710 | |
---|
711 | Count++; |
---|
712 | |
---|
713 | if(RMAX[sections]-RMIN[sections] >= perMillion) |
---|
714 | { |
---|
715 | MaxSurfaceVec[Count] = new G4FPlane( &EndPointList2, 0, ENormal ); |
---|
716 | } |
---|
717 | else |
---|
718 | { |
---|
719 | MaxSurfaceVec[Count] = 0; |
---|
720 | nb_of_surfaces--; |
---|
721 | } |
---|
722 | } |
---|
723 | |
---|
724 | // Now let's replicate the relevant surfaces into |
---|
725 | // G4BREPSolid's vector of surfaces |
---|
726 | // |
---|
727 | SurfaceVec = new G4Surface*[nb_of_surfaces]; |
---|
728 | G4int sf = 0; G4int zeroCount = 0; |
---|
729 | for( G4int srf = 0; srf < MaxNbOfSurfaces; srf++ ) |
---|
730 | { |
---|
731 | if( MaxSurfaceVec[srf] != 0 ) |
---|
732 | { |
---|
733 | if( sf < nb_of_surfaces ) |
---|
734 | { |
---|
735 | SurfaceVec[sf] = MaxSurfaceVec[srf]; |
---|
736 | } |
---|
737 | sf++; |
---|
738 | } |
---|
739 | else |
---|
740 | { |
---|
741 | zeroCount++; |
---|
742 | } |
---|
743 | } |
---|
744 | |
---|
745 | if( sf != nb_of_surfaces ) |
---|
746 | { |
---|
747 | G4cerr << "ERROR - G4BREPSolidPolyhedra::G4BREPSolidPolyhedra()" << G4endl |
---|
748 | << " Bad number of surfaces!" << G4endl |
---|
749 | << " sf : " << sf |
---|
750 | << " nb_of_surfaces: " << nb_of_surfaces |
---|
751 | << " Count : " << Count |
---|
752 | << G4endl; |
---|
753 | G4Exception( "G4BREPSolidPolyhedra::G4BREPSolidPolyhedra()", |
---|
754 | "FatalError", FatalException, |
---|
755 | "INTERNAL ERROR: Going bananas!" ); |
---|
756 | } |
---|
757 | |
---|
758 | // Clean up the temporary vector of surfaces |
---|
759 | // |
---|
760 | delete [] MaxSurfaceVec; |
---|
761 | |
---|
762 | // Store the original parameters, to be used in visualisation |
---|
763 | // Note radii are not scaled because this BREP uses the radius of the |
---|
764 | // circumscribed circle and also graphics_reps/G4Polyhedron uses the |
---|
765 | // radius of the circumscribed circle. |
---|
766 | |
---|
767 | // Save contructor parameters |
---|
768 | // |
---|
769 | constructorParams.start_angle = start_angle; |
---|
770 | constructorParams.opening_angle = opening_angle; |
---|
771 | constructorParams.sides = sides; |
---|
772 | constructorParams.num_z_planes = num_z_planes; |
---|
773 | constructorParams.z_start = z_start; |
---|
774 | constructorParams.z_values = 0; |
---|
775 | constructorParams.RMIN = 0; |
---|
776 | constructorParams.RMAX = 0; |
---|
777 | |
---|
778 | if( num_z_planes > 0 ) |
---|
779 | { |
---|
780 | constructorParams.z_values = new G4double[num_z_planes]; |
---|
781 | constructorParams.RMIN = new G4double[num_z_planes]; |
---|
782 | constructorParams.RMAX = new G4double[num_z_planes]; |
---|
783 | for( G4int idx = 0; idx < num_z_planes; idx++ ) |
---|
784 | { |
---|
785 | constructorParams.z_values[idx] = z_values[idx]; |
---|
786 | constructorParams.RMIN[idx] = RMIN[idx]; |
---|
787 | constructorParams.RMAX[idx] = RMAX[idx]; |
---|
788 | } |
---|
789 | } |
---|
790 | |
---|
791 | // z_values[0] should be equal to z_start, for consistency |
---|
792 | // with what the constructor does. |
---|
793 | // Otherwise the z_values that are shifted by (z_values[0] - z_start) , |
---|
794 | // because z_values are only used in the form |
---|
795 | // length = z_values[d+1] - z_values[d]; // JA Apr 2, 97 |
---|
796 | |
---|
797 | if( z_values[0] != z_start ) |
---|
798 | { |
---|
799 | G4cerr << "ERROR - G4BREPSolidPolyhedra::G4BREPSolidPolyhedra()" << G4endl |
---|
800 | << " Wrong solid parameters: " |
---|
801 | << " z_values[0]= " << z_values[0] << " is not equal to " |
---|
802 | << " z_start= " << z_start << "." << G4endl; |
---|
803 | G4Exception( "G4BREPSolidPolyhedra::G4BREPSolidPolyhedra()", |
---|
804 | "Notification", JustWarning, |
---|
805 | "Construction Error. z_values[0] must be equal to z_start!" ); |
---|
806 | constructorParams.z_values[0]= z_start; |
---|
807 | } |
---|
808 | |
---|
809 | active=1; |
---|
810 | Initialize(); |
---|
811 | } |
---|
812 | |
---|
813 | G4BREPSolidPolyhedra::G4BREPSolidPolyhedra( __void__& a ) |
---|
814 | : G4BREPSolid(a) |
---|
815 | { |
---|
816 | } |
---|
817 | |
---|
818 | G4BREPSolidPolyhedra::~G4BREPSolidPolyhedra() |
---|
819 | { |
---|
820 | if( constructorParams.num_z_planes > 0 ) |
---|
821 | { |
---|
822 | delete [] constructorParams.z_values; |
---|
823 | delete [] constructorParams.RMIN; |
---|
824 | delete [] constructorParams.RMAX; |
---|
825 | } |
---|
826 | } |
---|
827 | |
---|
828 | void G4BREPSolidPolyhedra::Initialize() |
---|
829 | { |
---|
830 | // Calc bounding box for solids and surfaces |
---|
831 | // Convert concave planes to convex |
---|
832 | // |
---|
833 | ShortestDistance=1000000; |
---|
834 | CheckSurfaceNormals(); |
---|
835 | if(!Box || !AxisBox) |
---|
836 | IsConvex(); |
---|
837 | |
---|
838 | CalcBBoxes(); |
---|
839 | } |
---|
840 | |
---|
841 | void G4BREPSolidPolyhedra::Reset() const |
---|
842 | { |
---|
843 | Active(1); |
---|
844 | ((G4BREPSolidPolyhedra*)this)->intersectionDistance=kInfinity; |
---|
845 | StartInside(0); |
---|
846 | for(register G4int a=0;a<nb_of_surfaces;a++) |
---|
847 | SurfaceVec[a]->Reset(); |
---|
848 | ShortestDistance = kInfinity; |
---|
849 | } |
---|
850 | |
---|
851 | EInside G4BREPSolidPolyhedra::Inside(register const G4ThreeVector& Pt) const |
---|
852 | { |
---|
853 | // This function find if the point Pt is inside, |
---|
854 | // outside or on the surface of the solid |
---|
855 | |
---|
856 | const G4double sqrHalfTolerance = kCarTolerance*kCarTolerance*0.25; |
---|
857 | |
---|
858 | G4Vector3D v(1, 0, 0.01); |
---|
859 | G4Vector3D Pttmp(Pt); |
---|
860 | G4Vector3D Vtmp(v); |
---|
861 | G4Ray r(Pttmp, Vtmp); |
---|
862 | |
---|
863 | // Check if point is inside the Polyhedra bounding box |
---|
864 | // |
---|
865 | if( !GetBBox()->Inside(Pttmp) ) |
---|
866 | { |
---|
867 | return kOutside; |
---|
868 | } |
---|
869 | |
---|
870 | // Set the surfaces to active again |
---|
871 | // |
---|
872 | Reset(); |
---|
873 | |
---|
874 | // Test if the bounding box of each surface is intersected |
---|
875 | // by the ray. If not, the surface is deactivated. |
---|
876 | // |
---|
877 | TestSurfaceBBoxes(r); |
---|
878 | |
---|
879 | G4int hits=0, samehit=0; |
---|
880 | |
---|
881 | for(G4int a=0; a < nb_of_surfaces; a++) |
---|
882 | { |
---|
883 | G4Surface* surface = SurfaceVec[a]; |
---|
884 | |
---|
885 | if(surface->IsActive()) |
---|
886 | { |
---|
887 | // count the number of intersections. |
---|
888 | // if this number is odd, the start of the ray is |
---|
889 | // inside the volume bounded by the surfaces, so |
---|
890 | // increment the number of intersection by 1 if the |
---|
891 | // point is not on the surface and if this intersection |
---|
892 | // was not found before |
---|
893 | // |
---|
894 | if( (surface->Intersect(r)) & 1 ) |
---|
895 | { |
---|
896 | // test if the point is on the surface |
---|
897 | // |
---|
898 | if(surface->GetDistance() < sqrHalfTolerance) |
---|
899 | { |
---|
900 | return kSurface; |
---|
901 | } |
---|
902 | // test if this intersection was found before |
---|
903 | // |
---|
904 | for(G4int i=0; i<a; i++) |
---|
905 | { |
---|
906 | if(surface->GetDistance() == SurfaceVec[i]->GetDistance()) |
---|
907 | { |
---|
908 | samehit++; |
---|
909 | break; |
---|
910 | } |
---|
911 | } |
---|
912 | |
---|
913 | // count the number of surfaces intersected by the ray |
---|
914 | // |
---|
915 | if(!samehit) |
---|
916 | { |
---|
917 | hits++; |
---|
918 | } |
---|
919 | } |
---|
920 | } |
---|
921 | } |
---|
922 | |
---|
923 | // if the number of surfaces intersected is odd, |
---|
924 | // the point is inside the solid |
---|
925 | // |
---|
926 | return ( (hits&1) ? kInside : kOutside ); |
---|
927 | } |
---|
928 | |
---|
929 | G4ThreeVector |
---|
930 | G4BREPSolidPolyhedra::SurfaceNormal(const G4ThreeVector& Pt) const |
---|
931 | { |
---|
932 | // This function calculates the normal of the closest surface |
---|
933 | // to the given point |
---|
934 | // Note : the sense of the normal depends on the sense of the surface |
---|
935 | |
---|
936 | G4int iplane; |
---|
937 | G4bool normflag = false; |
---|
938 | const G4double sqrHalfTolerance = kCarTolerance*kCarTolerance*0.25; |
---|
939 | |
---|
940 | // Determine if the point is on the surface |
---|
941 | // |
---|
942 | G4double minDist = kInfinity; |
---|
943 | G4int normPlane = 0; |
---|
944 | for(iplane = 0; iplane < nb_of_surfaces; iplane++) |
---|
945 | { |
---|
946 | G4double dist = std::fabs(SurfaceVec[iplane]->HowNear(Pt)); |
---|
947 | if( minDist > dist ) |
---|
948 | { |
---|
949 | minDist = dist; |
---|
950 | normPlane = iplane; |
---|
951 | } |
---|
952 | if( dist < sqrHalfTolerance) |
---|
953 | { |
---|
954 | // the point is on this surface, so take this as the |
---|
955 | // the surface to consider for computing the normal |
---|
956 | // |
---|
957 | normflag = true; |
---|
958 | break; |
---|
959 | } |
---|
960 | } |
---|
961 | |
---|
962 | // Calculate the normal at this point, if the point is on the |
---|
963 | // surface, otherwise compute the normal to the closest surface |
---|
964 | // |
---|
965 | if ( normflag ) // point on surface |
---|
966 | { |
---|
967 | G4ThreeVector norm = SurfaceVec[iplane]->SurfaceNormal(Pt); |
---|
968 | return norm.unit(); |
---|
969 | } |
---|
970 | else // point not on surface |
---|
971 | { |
---|
972 | G4FPlane* nPlane = (G4FPlane*)(SurfaceVec[normPlane]); |
---|
973 | G4ThreeVector hitPt = nPlane->GetSrfPoint(); |
---|
974 | G4ThreeVector hitNorm = nPlane->SurfaceNormal(hitPt); |
---|
975 | return hitNorm.unit(); |
---|
976 | } |
---|
977 | } |
---|
978 | |
---|
979 | G4double G4BREPSolidPolyhedra::DistanceToIn(const G4ThreeVector& Pt) const |
---|
980 | { |
---|
981 | // Calculates the shortest distance ("safety") from a point |
---|
982 | // outside the solid to any boundary of this solid. |
---|
983 | // Return 0 if the point is already inside. |
---|
984 | |
---|
985 | G4double *dists = new G4double[nb_of_surfaces]; |
---|
986 | G4int a; |
---|
987 | |
---|
988 | // Set the surfaces to active again |
---|
989 | // |
---|
990 | Reset(); |
---|
991 | |
---|
992 | // compute the shortest distance of the point to each surfaces |
---|
993 | // Be careful : it's a signed value |
---|
994 | // |
---|
995 | for(a=0; a< nb_of_surfaces; a++) |
---|
996 | dists[a] = SurfaceVec[a]->HowNear(Pt); |
---|
997 | |
---|
998 | G4double Dist = kInfinity; |
---|
999 | |
---|
1000 | // if dists[] is positive, the point is outside |
---|
1001 | // so take the shortest of the shortest positive distances |
---|
1002 | // dists[] can be equal to 0 : point on a surface |
---|
1003 | // ( Problem with the G4FPlane : there is no inside and no outside... |
---|
1004 | // So, to test if the point is inside to return 0, utilize the Inside |
---|
1005 | // function. But I don`t know if it is really needed because dToIn is |
---|
1006 | // called only if the point is outside ) |
---|
1007 | // |
---|
1008 | for(a = 0; a < nb_of_surfaces; a++) |
---|
1009 | if( std::fabs(Dist) > std::fabs(dists[a]) ) |
---|
1010 | //if( dists[a] >= 0) |
---|
1011 | Dist = dists[a]; |
---|
1012 | |
---|
1013 | delete[] dists; |
---|
1014 | |
---|
1015 | if(Dist == kInfinity) |
---|
1016 | { |
---|
1017 | // the point is inside the solid or on a surface |
---|
1018 | // |
---|
1019 | return 0; |
---|
1020 | } |
---|
1021 | else |
---|
1022 | { |
---|
1023 | return std::fabs(Dist); |
---|
1024 | } |
---|
1025 | } |
---|
1026 | |
---|
1027 | G4double |
---|
1028 | G4BREPSolidPolyhedra::DistanceToIn(register const G4ThreeVector& Pt, |
---|
1029 | register const G4ThreeVector& V) const |
---|
1030 | { |
---|
1031 | // Calculates the distance from a point outside the solid |
---|
1032 | // to the solid`s boundary along a specified direction vector. |
---|
1033 | // |
---|
1034 | // Note : Intersections with boundaries less than the |
---|
1035 | // tolerance must be ignored if the direction |
---|
1036 | // is away from the boundary |
---|
1037 | |
---|
1038 | G4int a; |
---|
1039 | |
---|
1040 | // Set the surfaces to active again |
---|
1041 | // |
---|
1042 | Reset(); |
---|
1043 | |
---|
1044 | const G4double sqrHalfTolerance = kCarTolerance*kCarTolerance*0.25; |
---|
1045 | G4Vector3D Pttmp(Pt); |
---|
1046 | G4Vector3D Vtmp(V); |
---|
1047 | G4Ray r(Pttmp, Vtmp); |
---|
1048 | |
---|
1049 | // Test if the bounding box of each surface is intersected |
---|
1050 | // by the ray. If not, the surface become deactive. |
---|
1051 | // |
---|
1052 | TestSurfaceBBoxes(r); |
---|
1053 | |
---|
1054 | ShortestDistance = kInfinity; |
---|
1055 | |
---|
1056 | for(a=0; a< nb_of_surfaces; a++) |
---|
1057 | { |
---|
1058 | if( SurfaceVec[a]->IsActive() ) |
---|
1059 | { |
---|
1060 | // test if the ray intersect the surface |
---|
1061 | // |
---|
1062 | if( SurfaceVec[a]->Intersect(r) ) |
---|
1063 | { |
---|
1064 | G4double surfDistance = SurfaceVec[a]->GetDistance(); |
---|
1065 | |
---|
1066 | // if more than 1 surface is intersected, |
---|
1067 | // take the nearest one |
---|
1068 | // |
---|
1069 | if( surfDistance < ShortestDistance ) |
---|
1070 | { |
---|
1071 | if( surfDistance > sqrHalfTolerance ) |
---|
1072 | { |
---|
1073 | ShortestDistance = surfDistance; |
---|
1074 | } |
---|
1075 | else |
---|
1076 | { |
---|
1077 | // the point is within the boundary |
---|
1078 | // ignore it if the direction is away from the boundary |
---|
1079 | // |
---|
1080 | G4Vector3D Norm = SurfaceVec[a]->SurfaceNormal(Pttmp); |
---|
1081 | |
---|
1082 | if( (Norm * Vtmp) < 0 ) |
---|
1083 | { |
---|
1084 | ShortestDistance = surfDistance; |
---|
1085 | // ShortestDistance = surfDistance==0 |
---|
1086 | // ? sqrHalfTolerance |
---|
1087 | // : surfDistance; |
---|
1088 | } |
---|
1089 | } |
---|
1090 | } |
---|
1091 | } |
---|
1092 | } |
---|
1093 | } |
---|
1094 | |
---|
1095 | // Be careful ! |
---|
1096 | // SurfaceVec->Distance is in fact the squared distance |
---|
1097 | // |
---|
1098 | if(ShortestDistance != kInfinity) |
---|
1099 | { |
---|
1100 | return std::sqrt(ShortestDistance); |
---|
1101 | } |
---|
1102 | else // no intersection |
---|
1103 | { |
---|
1104 | return kInfinity; |
---|
1105 | } |
---|
1106 | } |
---|
1107 | |
---|
1108 | G4double |
---|
1109 | G4BREPSolidPolyhedra::DistanceToOut(register const G4ThreeVector& Pt, |
---|
1110 | register const G4ThreeVector& V, |
---|
1111 | const G4bool, |
---|
1112 | G4bool *validNorm, |
---|
1113 | G4ThreeVector * ) const |
---|
1114 | { |
---|
1115 | // Calculates the distance from a point inside the solid |
---|
1116 | // to the solid`s boundary along a specified direction vector. |
---|
1117 | // Return 0 if the point is already outside (even number of |
---|
1118 | // intersections greater than the tolerance). |
---|
1119 | // |
---|
1120 | // Note : If the shortest distance to a boundary is less |
---|
1121 | // than the tolerance, it is ignored. This allows |
---|
1122 | // for a point within a tolerant boundary to leave |
---|
1123 | // immediately |
---|
1124 | |
---|
1125 | G4int parity = 0; |
---|
1126 | |
---|
1127 | // Set the surfaces to active again |
---|
1128 | // |
---|
1129 | Reset(); |
---|
1130 | |
---|
1131 | const G4double sqrHalfTolerance = kCarTolerance*kCarTolerance*0.25; |
---|
1132 | G4Vector3D Ptv = Pt; |
---|
1133 | G4int a; |
---|
1134 | |
---|
1135 | // I don`t understand this line |
---|
1136 | // |
---|
1137 | if(validNorm) |
---|
1138 | *validNorm=false; |
---|
1139 | |
---|
1140 | G4Vector3D Pttmp(Pt); |
---|
1141 | G4Vector3D Vtmp(V); |
---|
1142 | |
---|
1143 | G4Ray r(Pttmp, Vtmp); |
---|
1144 | |
---|
1145 | // Test if the bounding box of each surface is intersected |
---|
1146 | // by the ray. If not, the surface become deactive. |
---|
1147 | // |
---|
1148 | TestSurfaceBBoxes(r); |
---|
1149 | |
---|
1150 | ShortestDistance = kInfinity; // this is actually the square of the distance |
---|
1151 | |
---|
1152 | for(a=0; a< nb_of_surfaces; a++) |
---|
1153 | { |
---|
1154 | G4double surfDistance = SurfaceVec[a]->GetDistance(); |
---|
1155 | |
---|
1156 | if(SurfaceVec[a]->IsActive()) |
---|
1157 | { |
---|
1158 | G4int intersects = SurfaceVec[a]->Intersect(r); |
---|
1159 | |
---|
1160 | // test if the ray intersects the surface |
---|
1161 | // |
---|
1162 | if( intersects != 0 ) |
---|
1163 | { |
---|
1164 | parity += 1; |
---|
1165 | |
---|
1166 | // if more than 1 surface is intersected, take the nearest one |
---|
1167 | // |
---|
1168 | if( surfDistance < ShortestDistance ) |
---|
1169 | { |
---|
1170 | if( surfDistance > sqrHalfTolerance ) |
---|
1171 | { |
---|
1172 | ShortestDistance = surfDistance; |
---|
1173 | } |
---|
1174 | else |
---|
1175 | { |
---|
1176 | // the point is within the boundary: ignore it |
---|
1177 | // |
---|
1178 | parity -= 1; |
---|
1179 | } |
---|
1180 | } |
---|
1181 | } |
---|
1182 | } |
---|
1183 | } |
---|
1184 | |
---|
1185 | G4double distance = 0.; |
---|
1186 | |
---|
1187 | // Be careful ! |
---|
1188 | // SurfaceVec->Distance is in fact the squared distance |
---|
1189 | // |
---|
1190 | // This condition was changed in order to give not zero answer |
---|
1191 | // when particle is passing the border of two Touching Surfaces |
---|
1192 | // and the distance to this surfaces is not zero. |
---|
1193 | // parity is for the points on the boundary, |
---|
1194 | // parity is counting only surfDistance<kCarTolerance/2. |
---|
1195 | // |
---|
1196 | // if((ShortestDistance != kInfinity) && (parity&1)) |
---|
1197 | // |
---|
1198 | // |
---|
1199 | if((ShortestDistance != kInfinity) || (parity&1)) |
---|
1200 | { |
---|
1201 | distance = std::sqrt(ShortestDistance); |
---|
1202 | } |
---|
1203 | |
---|
1204 | return distance; |
---|
1205 | } |
---|
1206 | |
---|
1207 | G4double G4BREPSolidPolyhedra::DistanceToOut(const G4ThreeVector& Pt) const |
---|
1208 | { |
---|
1209 | // Calculates the shortest distance ("safety") from a point |
---|
1210 | // inside the solid to any boundary of this solid. |
---|
1211 | // Return 0 if the point is already outside. |
---|
1212 | |
---|
1213 | G4double *dists = new G4double[nb_of_surfaces]; |
---|
1214 | G4int a; |
---|
1215 | |
---|
1216 | // Set the surfaces to active again |
---|
1217 | // |
---|
1218 | Reset(); |
---|
1219 | |
---|
1220 | // calculate the shortest distance of the point to each surfaces |
---|
1221 | // Be careful : it's a signed value |
---|
1222 | // |
---|
1223 | for(a=0; a< nb_of_surfaces; a++) |
---|
1224 | { |
---|
1225 | dists[a] = SurfaceVec[a]->HowNear(Pt); |
---|
1226 | } |
---|
1227 | |
---|
1228 | G4double Dist = kInfinity; |
---|
1229 | |
---|
1230 | // if dists[] is negative, the point is inside |
---|
1231 | // so take the shortest of the shortest negative distances |
---|
1232 | // dists[] can be equal to 0 : point on a surface |
---|
1233 | // ( Problem with the G4FPlane : there is no inside and no outside... |
---|
1234 | // So, to test if the point is outside to return 0, utilize the Inside |
---|
1235 | // function. But I don`t know if it is really needed because dToOut is |
---|
1236 | // called only if the point is inside ) |
---|
1237 | |
---|
1238 | for(a = 0; a < nb_of_surfaces; a++) |
---|
1239 | { |
---|
1240 | if( std::fabs(Dist) > std::fabs(dists[a]) ) |
---|
1241 | { |
---|
1242 | //if( dists[a] <= 0) |
---|
1243 | Dist = dists[a]; |
---|
1244 | } |
---|
1245 | } |
---|
1246 | |
---|
1247 | delete[] dists; |
---|
1248 | |
---|
1249 | if(Dist == kInfinity) |
---|
1250 | { |
---|
1251 | // the point is ouside the solid or on a surface |
---|
1252 | // |
---|
1253 | return 0; |
---|
1254 | } |
---|
1255 | else |
---|
1256 | { |
---|
1257 | // return Dist; |
---|
1258 | return std::fabs(Dist); |
---|
1259 | } |
---|
1260 | } |
---|
1261 | |
---|
1262 | std::ostream& G4BREPSolidPolyhedra::StreamInfo(std::ostream& os) const |
---|
1263 | { |
---|
1264 | |
---|
1265 | // Streams solid contents to output stream. |
---|
1266 | |
---|
1267 | G4BREPSolid::StreamInfo( os ) |
---|
1268 | << "\n start_angle: " << constructorParams.start_angle |
---|
1269 | << "\n opening_angle: " << constructorParams.opening_angle |
---|
1270 | << "\n sides: " << constructorParams.sides |
---|
1271 | << "\n num_z_planes: " << constructorParams.num_z_planes |
---|
1272 | << "\n z_start: " << constructorParams.z_start |
---|
1273 | << "\n z_values: "; |
---|
1274 | G4int idx; |
---|
1275 | for( idx = 0; idx < constructorParams.num_z_planes; idx++ ) |
---|
1276 | { |
---|
1277 | os << constructorParams.z_values[idx] << " "; |
---|
1278 | } |
---|
1279 | os << "\n RMIN: "; |
---|
1280 | for( idx = 0; idx < constructorParams.num_z_planes; idx++ ) |
---|
1281 | { |
---|
1282 | os << constructorParams.RMIN[idx] << " "; |
---|
1283 | } |
---|
1284 | os << "\n RMAX: "; |
---|
1285 | for( idx = 0; idx < constructorParams.num_z_planes; idx++ ) |
---|
1286 | { |
---|
1287 | os << constructorParams.RMAX[idx] << " "; |
---|
1288 | } |
---|
1289 | os << "\n-----------------------------------------------------------\n"; |
---|
1290 | |
---|
1291 | return os; |
---|
1292 | } |
---|
1293 | |
---|
1294 | G4Surface* |
---|
1295 | G4BREPSolidPolyhedra::CreateTrapezoidalSurface( G4double r1, |
---|
1296 | G4double r2, |
---|
1297 | const G4Point3D& origin, |
---|
1298 | G4double distance, |
---|
1299 | G4Vector3D& xAxis, |
---|
1300 | G4double partAngle, |
---|
1301 | ESurfaceSense sense ) |
---|
1302 | { |
---|
1303 | // The surface to be returned |
---|
1304 | // |
---|
1305 | G4Surface* trapsrf = 0; |
---|
1306 | G4Point3DVector PointList(4); |
---|
1307 | G4Vector3D zAxis(0,0,1); |
---|
1308 | |
---|
1309 | PointList[0] = origin + ( r1 * xAxis); |
---|
1310 | PointList[3] = origin + ( distance * zAxis) + (r2 * xAxis); |
---|
1311 | |
---|
1312 | xAxis.rotateZ( partAngle ); |
---|
1313 | |
---|
1314 | PointList[2] = origin + ( distance * zAxis) + (r2 * xAxis); |
---|
1315 | PointList[1] = origin + ( r1 * xAxis); |
---|
1316 | |
---|
1317 | // Return the planar trapezoidal surface |
---|
1318 | // |
---|
1319 | trapsrf = new G4FPlane( &PointList, 0, sense ); |
---|
1320 | |
---|
1321 | return trapsrf; |
---|
1322 | } |
---|
1323 | |
---|
1324 | G4Surface* |
---|
1325 | G4BREPSolidPolyhedra::CreateTriangularSurface( G4double r1, |
---|
1326 | G4double r2, |
---|
1327 | const G4Point3D& origin, |
---|
1328 | G4double distance, |
---|
1329 | G4Vector3D& xAxis, |
---|
1330 | G4double partAngle, |
---|
1331 | ESurfaceSense sense ) |
---|
1332 | { |
---|
1333 | // The surface to be returned |
---|
1334 | // |
---|
1335 | G4Surface* trapsrf = 0; |
---|
1336 | G4Point3DVector PointList(3); |
---|
1337 | G4Vector3D zAxis(0,0,1); |
---|
1338 | |
---|
1339 | PointList[0] = origin + ( r1 * xAxis); |
---|
1340 | PointList[2] = origin + ( distance * zAxis) + (r2 * xAxis); |
---|
1341 | |
---|
1342 | xAxis.rotateZ( partAngle ); |
---|
1343 | |
---|
1344 | if( r1 < r2 ) |
---|
1345 | { |
---|
1346 | PointList[1] = origin + ( distance * zAxis) + (r2 * xAxis); |
---|
1347 | } |
---|
1348 | else |
---|
1349 | { |
---|
1350 | PointList[1] = origin + ( r1 * xAxis); |
---|
1351 | } |
---|
1352 | |
---|
1353 | // Return the planar trapezoidal surface |
---|
1354 | // |
---|
1355 | trapsrf = new G4FPlane( &PointList, 0, sense ); |
---|
1356 | |
---|
1357 | return trapsrf; |
---|
1358 | } |
---|
1359 | |
---|
1360 | G4Surface* |
---|
1361 | G4BREPSolidPolyhedra::ComputePlanarSurface( G4double r1, |
---|
1362 | G4double r2, |
---|
1363 | const G4Point3D& origin, |
---|
1364 | G4Vector3D& xAxis, |
---|
1365 | G4int sides, |
---|
1366 | G4double partAngle, |
---|
1367 | ESurfaceSense sense ) |
---|
1368 | { |
---|
1369 | // This method can be called only when r1 != r2, |
---|
1370 | // otherwise it returns 0 which means that no surface can be |
---|
1371 | // created out of the given radius pair. |
---|
1372 | // This method requires the xAxis to be pre-rotated properly. |
---|
1373 | |
---|
1374 | G4Point3DVector OuterPointList( sides ); |
---|
1375 | G4Point3DVector InnerPointList( sides ); |
---|
1376 | |
---|
1377 | G4double rIn, rOut; |
---|
1378 | G4Surface* planarSrf = 0; |
---|
1379 | |
---|
1380 | if( r1 < r2 ) |
---|
1381 | { |
---|
1382 | rIn = r1; |
---|
1383 | rOut = r2; |
---|
1384 | } |
---|
1385 | else if( r1 > r2 ) |
---|
1386 | { |
---|
1387 | rIn = r2; |
---|
1388 | rOut = r1; |
---|
1389 | } |
---|
1390 | else |
---|
1391 | { |
---|
1392 | // Invalid precondition, the radius values are r1 == r2, |
---|
1393 | // which means we can create only polyline but no surface |
---|
1394 | // |
---|
1395 | return 0; |
---|
1396 | } |
---|
1397 | |
---|
1398 | for( G4int pidx = 0; pidx < sides; pidx++ ) |
---|
1399 | { |
---|
1400 | // Outer polyline |
---|
1401 | // |
---|
1402 | OuterPointList[pidx] = origin + ( rOut * xAxis); |
---|
1403 | // Inner polyline |
---|
1404 | // |
---|
1405 | InnerPointList[pidx] = origin + ( rIn * xAxis); |
---|
1406 | xAxis.rotateZ( partAngle ); |
---|
1407 | } |
---|
1408 | |
---|
1409 | if( rIn != 0.0 && rOut != 0.0 ) |
---|
1410 | { |
---|
1411 | // Standard case |
---|
1412 | // |
---|
1413 | planarSrf = new G4FPlane( &OuterPointList, &InnerPointList, sense ); |
---|
1414 | } |
---|
1415 | else if( rOut != 0.0 ) |
---|
1416 | { |
---|
1417 | // Special case where inner radius is zero so no polyline |
---|
1418 | // is actually created |
---|
1419 | // |
---|
1420 | planarSrf = new G4FPlane( &OuterPointList, 0, sense ); |
---|
1421 | } |
---|
1422 | else |
---|
1423 | { |
---|
1424 | // No surface being created |
---|
1425 | // This should not happen as filtered out by precondition check above |
---|
1426 | } |
---|
1427 | |
---|
1428 | return planarSrf; |
---|
1429 | } |
---|
1430 | |
---|
1431 | // In graphics_reps: |
---|
1432 | |
---|
1433 | #include "G4Polyhedron.hh" |
---|
1434 | |
---|
1435 | G4Polyhedron* G4BREPSolidPolyhedra::CreatePolyhedron() const |
---|
1436 | { |
---|
1437 | return new G4PolyhedronPgon( constructorParams.start_angle, |
---|
1438 | constructorParams.opening_angle, |
---|
1439 | constructorParams.sides, |
---|
1440 | constructorParams.num_z_planes, |
---|
1441 | constructorParams.z_values, |
---|
1442 | constructorParams.RMIN, |
---|
1443 | constructorParams.RMAX); |
---|
1444 | } |
---|