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: G4BREPSolidPCone.cc,v 1.38 2006/06/29 18:41:24 gunter Exp $ |
---|
27 | // GEANT4 tag $Name: HEAD $ |
---|
28 | // |
---|
29 | // ---------------------------------------------------------------------- |
---|
30 | // GEANT 4 class source file |
---|
31 | // |
---|
32 | // G4BREPSolidPCone.cc |
---|
33 | // |
---|
34 | // ---------------------------------------------------------------------- |
---|
35 | // The polyconical solid G4BREPSolidPCone is a shape defined by a set of |
---|
36 | // inner and outer conical or cylindrical surface sections and two planes |
---|
37 | // perpendicular to the Z axis. Each conical surface is defined by its |
---|
38 | // radius at two different planes perpendicular to the Z-axis. Inner and |
---|
39 | // outer conical surfaces are defined using common Z planes. |
---|
40 | // ---------------------------------------------------------------------- |
---|
41 | |
---|
42 | #include "G4Types.hh" |
---|
43 | #include <sstream> |
---|
44 | |
---|
45 | #include "G4BREPSolidPCone.hh" |
---|
46 | #include "G4FCylindricalSurface.hh" |
---|
47 | #include "G4FConicalSurface.hh" |
---|
48 | #include "G4CircularCurve.hh" |
---|
49 | #include "G4FPlane.hh" |
---|
50 | |
---|
51 | typedef enum |
---|
52 | { |
---|
53 | EInverse = 0, |
---|
54 | ENormal = 1 |
---|
55 | } ESurfaceSense; |
---|
56 | |
---|
57 | G4BREPSolidPCone::G4BREPSolidPCone(const G4String& name, |
---|
58 | G4double start_angle, |
---|
59 | G4double opening_angle, |
---|
60 | G4int num_z_planes, // sections, |
---|
61 | G4double z_start, |
---|
62 | G4double z_values[], |
---|
63 | G4double RMIN[], |
---|
64 | G4double RMAX[] ) |
---|
65 | : G4BREPSolid(name) |
---|
66 | { |
---|
67 | G4int sections= num_z_planes-1; |
---|
68 | nb_of_surfaces = 2*sections+2; |
---|
69 | SurfaceVec = new G4Surface*[nb_of_surfaces]; |
---|
70 | G4ThreeVector Axis(0,0,1); |
---|
71 | G4ThreeVector Origin(0,0,z_start); |
---|
72 | G4double Length; |
---|
73 | G4ThreeVector LocalOrigin(0,0,z_start); |
---|
74 | G4int a, b = 0; |
---|
75 | |
---|
76 | G4ThreeVector PlaneAxis(0, 0, 1); |
---|
77 | G4ThreeVector PlaneDir (0, 1, 0); |
---|
78 | |
---|
79 | /////////////////////////////////////////////////// |
---|
80 | // Test delta phi |
---|
81 | |
---|
82 | // At the moment (11/03/2002) the phi section is not implemented |
---|
83 | // so we take a G4 application down if there is a request for such |
---|
84 | // a configuration |
---|
85 | // |
---|
86 | if( opening_angle < 2*pi-perMillion ) |
---|
87 | { |
---|
88 | G4Exception("G4BREPSolidPCone::G4BREPSolidPCone()", |
---|
89 | "NotImplemented", FatalException, |
---|
90 | "Sorry, phi-section not supported yet, try to use G4Polycone instead!"); |
---|
91 | } |
---|
92 | |
---|
93 | /////////////////////////////////////////////////// |
---|
94 | // Test the validity of the R values |
---|
95 | |
---|
96 | // RMIN[0] and RMIN[num_z_planes-1] cannot be = 0 |
---|
97 | // when RMIN[0] or RMIN[num_z_planes-1] are = 0 |
---|
98 | // |
---|
99 | if( ((RMIN[0] == 0) && (RMAX[0] == 0)) || |
---|
100 | ((RMIN[num_z_planes-1] == 0) && (RMAX[num_z_planes-1] == 0)) ) |
---|
101 | G4Exception("G4BREPSolidPCone::G4BREPSolidPCone()", |
---|
102 | "InvalidSetup", FatalException, |
---|
103 | "RMIN at the extremities cannot be 0 when RMAX=0 !"); |
---|
104 | |
---|
105 | // only RMAX[0] and RMAX[num_z_planes-1] can be = 0 |
---|
106 | // |
---|
107 | for(a = 1; a < num_z_planes-1; a++) |
---|
108 | if (RMAX[a] == 0) |
---|
109 | G4Exception("G4BREPSolidPCone::G4BREPSolidPCone()", |
---|
110 | "InvalidSetup", FatalException, |
---|
111 | "RMAX inside the solid cannot be 0 !"); |
---|
112 | |
---|
113 | // RMAX[a] must be greater than RMIN[a] |
---|
114 | // |
---|
115 | for(a = 1; a < num_z_planes-1; a++) |
---|
116 | if (RMIN[a] >= RMAX[a]) |
---|
117 | G4Exception("G4BREPSolidPCone::G4BREPSolidPCone()", |
---|
118 | "InvalidSetup", FatalException, |
---|
119 | "RMAX must be greater than RMIN in the middle Z planes !"); |
---|
120 | |
---|
121 | if( (RMIN[num_z_planes-1] > RMAX[num_z_planes-1] ) |
---|
122 | || (RMIN[0] > RMAX[0]) ) |
---|
123 | G4Exception("G4BREPSolidPCone::G4BREPSolidPCone()", |
---|
124 | "InvalidSetup", FatalException, |
---|
125 | "RMAX must be greater or equal than RMIN at the ends !"); |
---|
126 | |
---|
127 | // Create surfaces |
---|
128 | |
---|
129 | for( a = 0; a < sections; a++) |
---|
130 | { |
---|
131 | // Surface length |
---|
132 | // |
---|
133 | Length = z_values[a+1] - z_values[a]; |
---|
134 | |
---|
135 | if( Length == 0 ) |
---|
136 | { |
---|
137 | // We need to create planar surface(s) |
---|
138 | // |
---|
139 | if( RMAX[a] != RMAX[a+1] && RMIN[a] != RMIN[a+1] ) |
---|
140 | { |
---|
141 | // We can have the 8 following configurations here: |
---|
142 | // |
---|
143 | // 1. 2. 3. 4. |
---|
144 | // --+ +-- --+ +-- |
---|
145 | // xx|-> <-|xx xx| |xx |
---|
146 | // xx+-- --+xx --+ +-- |
---|
147 | // xxxxx xxxxx | | |
---|
148 | // xxxxx xxxxx +-- --+ |
---|
149 | // xx+-- --+xx |xx xx| |
---|
150 | // xx|-> <-|xx +-- --+ |
---|
151 | // --+ +-- |
---|
152 | // -------------------------- Z axis |
---|
153 | // |
---|
154 | ////////////////////////////////////////////////////////////// |
---|
155 | ////////////////////////////////////////////////////////////// |
---|
156 | // |
---|
157 | // 5. 6. 7. 8. |
---|
158 | // --+ +-- --+ +-- |
---|
159 | // xx|-> <-|xx xx|-> <-|xx |
---|
160 | // --+-- --+-- xx+-- --+xx |
---|
161 | // <-|xx xx|-> xxxxx xxxxx |
---|
162 | // +-- --+ --+xx xx+-- |
---|
163 | // <-|xx xx|-> |
---|
164 | // +-- --+ |
---|
165 | // -------------------------- Z axis |
---|
166 | // |
---|
167 | // NOTE: The pictures show only one half of polycone! |
---|
168 | // The arrows show the expected surface normal direction. |
---|
169 | // The configuration No. 3 and 4 are not valid solids! |
---|
170 | |
---|
171 | // Eliminate the invalid cases 3 and 4. |
---|
172 | // At this point is guaranteed that each RMIN[i] < RMAX[i] |
---|
173 | // where i in in interval 0 < i < num_z_planes-1. So: |
---|
174 | // |
---|
175 | if( RMIN[a] > RMAX[a+1] || RMAX[a] < RMIN[a+1] ) |
---|
176 | { |
---|
177 | std::ostringstream os; |
---|
178 | os << "The values of RMIN[" << a |
---|
179 | << "] & RMAX[" << a+1 << "] or RMAX[" << a |
---|
180 | << "] & RMIN[" << a+1 << "] " |
---|
181 | << "generate an invalid configuration for solid: " |
---|
182 | << name.c_str() << "!" << G4endl; |
---|
183 | G4String message = os.str(); |
---|
184 | G4Exception("G4BREPSolidPCone::G4BREPSolidPCone()", |
---|
185 | "InvalidSetup", FatalException, message ); |
---|
186 | } |
---|
187 | |
---|
188 | ESurfaceSense UpSurfSense, LowSurfSense; |
---|
189 | |
---|
190 | // We need to classify all the cases in order to figure out |
---|
191 | // the planar surface sense |
---|
192 | // |
---|
193 | if( RMAX[a] > RMAX[a+1] ) |
---|
194 | { |
---|
195 | // Cases 1, 5, 7 |
---|
196 | // |
---|
197 | if( RMIN[a] < RMIN[a+1] ) |
---|
198 | { |
---|
199 | // Case 1 |
---|
200 | // |
---|
201 | UpSurfSense = ENormal; |
---|
202 | LowSurfSense = ENormal; |
---|
203 | } |
---|
204 | else if( RMAX[a+1] != RMIN[a]) |
---|
205 | { |
---|
206 | // Case 7 |
---|
207 | // |
---|
208 | UpSurfSense = ENormal; |
---|
209 | LowSurfSense = EInverse; |
---|
210 | } |
---|
211 | else |
---|
212 | { |
---|
213 | // Case 5 |
---|
214 | // |
---|
215 | UpSurfSense = ENormal; |
---|
216 | LowSurfSense = EInverse; |
---|
217 | } |
---|
218 | } |
---|
219 | else |
---|
220 | { |
---|
221 | // Cases 2, 6, 8 |
---|
222 | // |
---|
223 | if( RMIN[a] > RMIN[a+1] ) |
---|
224 | { |
---|
225 | // Case 2 |
---|
226 | // |
---|
227 | UpSurfSense = EInverse; |
---|
228 | LowSurfSense = EInverse; |
---|
229 | } |
---|
230 | else if( RMIN[a+1] != RMAX[a] ) |
---|
231 | { |
---|
232 | // Case 8 |
---|
233 | // |
---|
234 | UpSurfSense = EInverse; |
---|
235 | LowSurfSense = ENormal; |
---|
236 | } |
---|
237 | else |
---|
238 | { |
---|
239 | // Case 6 |
---|
240 | UpSurfSense = EInverse; |
---|
241 | LowSurfSense = ENormal; |
---|
242 | } |
---|
243 | } |
---|
244 | |
---|
245 | SurfaceVec[b] = |
---|
246 | ComputePlanarSurface( RMAX[a], RMAX[a+1], LocalOrigin, PlaneAxis, |
---|
247 | PlaneDir, UpSurfSense ); |
---|
248 | //SurfaceVec[b]->SetSameSense( UpSurfSense ); |
---|
249 | b++; |
---|
250 | |
---|
251 | SurfaceVec[b] = |
---|
252 | ComputePlanarSurface( RMIN[a], RMIN[a+1], LocalOrigin, PlaneAxis, |
---|
253 | PlaneDir, LowSurfSense ); |
---|
254 | //SurfaceVec[b]->SetSameSense( LowSurfSense ); |
---|
255 | b++; |
---|
256 | } |
---|
257 | else |
---|
258 | { |
---|
259 | // The original code creating single planar surface |
---|
260 | // in case where only either RMAX or RMIN have changed |
---|
261 | // at the same Z value; e.g.: |
---|
262 | // |
---|
263 | // RMAX RMIN |
---|
264 | // change change |
---|
265 | // |
---|
266 | // 1 2 3 4 |
---|
267 | // --+ +-- ----- ----- |
---|
268 | // 00|-> <-|00 00000 00000 |
---|
269 | // 00+-- --+00 --+00 00+-- |
---|
270 | // 00000 00000 <-|00 00|-> |
---|
271 | // +-- --+ |
---|
272 | // --------------------------- Z axis |
---|
273 | // |
---|
274 | // NOTE: The picture shows only one half of polycone! |
---|
275 | |
---|
276 | G4double R1, R2; |
---|
277 | ESurfaceSense SurfSense; |
---|
278 | |
---|
279 | // test where is the plane surface |
---|
280 | // if( RMAX[a] != RMAX[a+1] ) |
---|
281 | // { |
---|
282 | // R1 = RMAX[a]; |
---|
283 | // R2 = RMAX[a+1]; |
---|
284 | // } |
---|
285 | // else if(RMIN[a] != RMIN[a+1]) |
---|
286 | // { |
---|
287 | // R1 = RMIN[a]; |
---|
288 | // R2 = RMIN[a+1]; |
---|
289 | // } |
---|
290 | // else |
---|
291 | // { |
---|
292 | // G4cerr << "Error in construction of G4BREPSolidPCone. \n" |
---|
293 | // << "Exactly the same z, rmin and rmax given for \n" |
---|
294 | // << "consecutive indices, " << a << " and " << a+1 << G4endl; |
---|
295 | // continue; |
---|
296 | // } |
---|
297 | if( RMAX[a] != RMAX[a+1] ) |
---|
298 | { |
---|
299 | // Cases 1, 2 |
---|
300 | // |
---|
301 | R1 = RMAX[a]; |
---|
302 | R2 = RMAX[a+1]; |
---|
303 | if( R1 > R2 ) |
---|
304 | { |
---|
305 | // Case 1 |
---|
306 | // |
---|
307 | SurfSense = ENormal; |
---|
308 | } |
---|
309 | else |
---|
310 | { |
---|
311 | // Case 2 |
---|
312 | // |
---|
313 | SurfSense = EInverse; |
---|
314 | } |
---|
315 | } |
---|
316 | else if(RMIN[a] != RMIN[a+1]) |
---|
317 | { |
---|
318 | // Cases 1, 2 |
---|
319 | // |
---|
320 | R1 = RMIN[a]; |
---|
321 | R2 = RMIN[a+1]; |
---|
322 | if( R1 > R2 ) |
---|
323 | { |
---|
324 | // Case 3 |
---|
325 | // |
---|
326 | SurfSense = EInverse; |
---|
327 | } |
---|
328 | else |
---|
329 | { |
---|
330 | // Case 4 |
---|
331 | // |
---|
332 | SurfSense = ENormal; |
---|
333 | } |
---|
334 | } |
---|
335 | else |
---|
336 | { |
---|
337 | G4cerr << "Error in construction of G4BREPSolidPCone. \n" |
---|
338 | << "Exactly the same z, rmin and rmax given for \n" |
---|
339 | << "consecutive indices, " << a << " and " << a+1 << G4endl; |
---|
340 | continue; |
---|
341 | } |
---|
342 | |
---|
343 | SurfaceVec[b] = |
---|
344 | ComputePlanarSurface( R1, R2, LocalOrigin, PlaneAxis, |
---|
345 | PlaneDir, SurfSense ); |
---|
346 | //SurfaceVec[b]->SetSameSense( SurfSense ); |
---|
347 | b++; |
---|
348 | |
---|
349 | // SurfaceVec[b]->SetSameSense(1); |
---|
350 | nb_of_surfaces--; |
---|
351 | } |
---|
352 | } |
---|
353 | else |
---|
354 | { |
---|
355 | // The surface to create is conical or cylindrical |
---|
356 | |
---|
357 | // Inner PCone |
---|
358 | // |
---|
359 | if(RMIN[a] != RMIN[a+1]) |
---|
360 | { |
---|
361 | // Create cone |
---|
362 | // |
---|
363 | if(RMIN[a] > RMIN[a+1]) |
---|
364 | { |
---|
365 | G4Vector3D ConeOrigin = G4Vector3D(LocalOrigin) ; |
---|
366 | SurfaceVec[b] = new G4FConicalSurface(ConeOrigin, Axis, Length, |
---|
367 | RMIN[a+1], RMIN[a]); |
---|
368 | SurfaceVec[b]->SetSameSense(0); |
---|
369 | } |
---|
370 | else |
---|
371 | { |
---|
372 | G4Vector3D Axis2 = G4Vector3D( -1*Axis ); |
---|
373 | G4Vector3D LocalOrigin2 = G4Vector3D( LocalOrigin + (Length*Axis) ); |
---|
374 | G4Vector3D ConeOrigin = LocalOrigin2; |
---|
375 | SurfaceVec[b] = new G4FConicalSurface(ConeOrigin, Axis2, Length, |
---|
376 | RMIN[a], RMIN[a+1]); |
---|
377 | SurfaceVec[b]->SetSameSense(0); |
---|
378 | } |
---|
379 | b++; |
---|
380 | } |
---|
381 | else |
---|
382 | { |
---|
383 | if( RMIN[a] == 0 ) |
---|
384 | { |
---|
385 | // Do not create any surface and decrease nb_of_surfaces |
---|
386 | // |
---|
387 | nb_of_surfaces--; |
---|
388 | } |
---|
389 | else |
---|
390 | { |
---|
391 | // Create cylinder |
---|
392 | // |
---|
393 | G4Vector3D CylOrigin = G4Vector3D( LocalOrigin ); |
---|
394 | SurfaceVec[b] = new G4FCylindricalSurface(CylOrigin, Axis, |
---|
395 | RMIN[a], Length ); |
---|
396 | SurfaceVec[b]->SetSameSense(0); |
---|
397 | b++; |
---|
398 | } |
---|
399 | } |
---|
400 | |
---|
401 | // Outer PCone |
---|
402 | // |
---|
403 | if(RMAX[a] != RMAX[a+1]) |
---|
404 | { |
---|
405 | // Create cone |
---|
406 | // |
---|
407 | if(RMAX[a] > RMAX[a+1]) |
---|
408 | { |
---|
409 | G4Vector3D ConeOrigin = G4Vector3D( LocalOrigin ); |
---|
410 | SurfaceVec[b] = new G4FConicalSurface(ConeOrigin, Axis, Length, RMAX[a+1], RMAX[a]); |
---|
411 | SurfaceVec[b]->SetSameSense(1); |
---|
412 | } |
---|
413 | else |
---|
414 | { |
---|
415 | G4Vector3D Axis2 = G4Vector3D( -1*Axis ); |
---|
416 | G4Vector3D LocalOrigin2 = G4Vector3D( LocalOrigin + (Length*Axis) ); |
---|
417 | G4Vector3D ConeOrigin = LocalOrigin2; |
---|
418 | SurfaceVec[b] = new G4FConicalSurface(ConeOrigin, Axis2, Length, |
---|
419 | RMAX[a], RMAX[a+1]); |
---|
420 | SurfaceVec[b]->SetSameSense(1); |
---|
421 | } |
---|
422 | b++; |
---|
423 | } |
---|
424 | else |
---|
425 | { |
---|
426 | // Create cylinder |
---|
427 | // |
---|
428 | G4Vector3D CylOrigin = G4Vector3D( LocalOrigin ); |
---|
429 | |
---|
430 | if (RMAX[a] == 0) |
---|
431 | { |
---|
432 | // Do not create any surface and decrease nb_of_surfaces |
---|
433 | // |
---|
434 | nb_of_surfaces--; |
---|
435 | } |
---|
436 | else |
---|
437 | { |
---|
438 | // Create cylinder |
---|
439 | // |
---|
440 | G4Vector3D CylOrigin = G4Vector3D( LocalOrigin ); |
---|
441 | SurfaceVec[b] = new G4FCylindricalSurface(CylOrigin, Axis, |
---|
442 | RMAX[a], Length ); |
---|
443 | SurfaceVec[b]->SetSameSense(1); |
---|
444 | b++; |
---|
445 | } |
---|
446 | } |
---|
447 | } |
---|
448 | |
---|
449 | // Move surface origin to next section |
---|
450 | // |
---|
451 | LocalOrigin = LocalOrigin + (Length*Axis); |
---|
452 | } |
---|
453 | |
---|
454 | /////////////////////////////////////////////////// |
---|
455 | // Create two end planes |
---|
456 | |
---|
457 | G4CurveVector cv; |
---|
458 | G4CircularCurve* tmp; |
---|
459 | |
---|
460 | if(RMIN[0] < RMAX[0]) |
---|
461 | { |
---|
462 | // Create start G4Plane & boundaries |
---|
463 | // |
---|
464 | G4Point3D ArcStart1a = G4Point3D( Origin + (RMIN[0]*PlaneDir) ); |
---|
465 | G4Point3D ArcStart1b = G4Point3D( Origin + (RMAX[0]*PlaneDir) ); |
---|
466 | |
---|
467 | if( RMIN[0] > 0.0 ) |
---|
468 | { |
---|
469 | tmp = new G4CircularCurve; |
---|
470 | tmp->Init(G4Axis2Placement3D(PlaneDir, PlaneAxis, Origin), RMIN[0]); |
---|
471 | tmp->SetBounds(ArcStart1a, ArcStart1a); |
---|
472 | tmp->SetSameSense(0); |
---|
473 | cv.push_back(tmp); |
---|
474 | } |
---|
475 | |
---|
476 | tmp = new G4CircularCurve; |
---|
477 | tmp->Init(G4Axis2Placement3D(PlaneDir, PlaneAxis, Origin), RMAX[0]); |
---|
478 | tmp->SetBounds(ArcStart1b, ArcStart1b); |
---|
479 | tmp->SetSameSense(1); |
---|
480 | cv.push_back(tmp); |
---|
481 | |
---|
482 | SurfaceVec[nb_of_surfaces-2] = new G4FPlane(PlaneDir, -PlaneAxis, Origin); |
---|
483 | SurfaceVec[nb_of_surfaces-2]->SetBoundaries(&cv); |
---|
484 | SurfaceVec[nb_of_surfaces-2]->SetSameSense(0); |
---|
485 | } |
---|
486 | else |
---|
487 | { |
---|
488 | // RMIN[0] == RMAX[0] so no surface is needed, it is a line! |
---|
489 | // |
---|
490 | nb_of_surfaces--; |
---|
491 | } |
---|
492 | |
---|
493 | if(RMIN[sections] < RMAX[sections]) |
---|
494 | { |
---|
495 | // Create end G4Plane & boundaries |
---|
496 | // |
---|
497 | G4Point3D ArcStart2a = G4Point3D( LocalOrigin+(RMIN[sections]*PlaneDir) ); |
---|
498 | G4Point3D ArcStart2b = G4Point3D( LocalOrigin+(RMAX[sections]*PlaneDir) ); |
---|
499 | |
---|
500 | cv.clear(); |
---|
501 | |
---|
502 | if( RMIN[sections] > 0.0 ) |
---|
503 | { |
---|
504 | tmp = new G4CircularCurve; |
---|
505 | tmp->Init(G4Axis2Placement3D(PlaneDir, PlaneAxis, LocalOrigin), |
---|
506 | RMIN[sections]); |
---|
507 | tmp->SetBounds(ArcStart2a, ArcStart2a); |
---|
508 | tmp->SetSameSense(0); |
---|
509 | cv.push_back(tmp); |
---|
510 | } |
---|
511 | |
---|
512 | tmp = new G4CircularCurve; |
---|
513 | tmp->Init(G4Axis2Placement3D(PlaneDir, PlaneAxis, LocalOrigin), |
---|
514 | RMAX[sections]); |
---|
515 | tmp->SetBounds(ArcStart2b, ArcStart2b); |
---|
516 | tmp->SetSameSense(1); |
---|
517 | cv.push_back(tmp); |
---|
518 | |
---|
519 | SurfaceVec[nb_of_surfaces-1] = new G4FPlane(PlaneDir, PlaneAxis, |
---|
520 | LocalOrigin); |
---|
521 | SurfaceVec[nb_of_surfaces-1]->SetBoundaries(&cv); |
---|
522 | |
---|
523 | // set sense of the surface |
---|
524 | // |
---|
525 | SurfaceVec[nb_of_surfaces-1]->SetSameSense(0); |
---|
526 | } |
---|
527 | else |
---|
528 | { |
---|
529 | // RMIN[0] == RMAX[0] so no surface is needed, it is a line! |
---|
530 | // |
---|
531 | nb_of_surfaces--; |
---|
532 | } |
---|
533 | |
---|
534 | // Save contructor parameters |
---|
535 | // |
---|
536 | constructorParams.start_angle = start_angle; |
---|
537 | constructorParams.opening_angle = opening_angle; |
---|
538 | constructorParams.num_z_planes = num_z_planes; |
---|
539 | constructorParams.z_start = z_start; |
---|
540 | constructorParams.z_values = 0; |
---|
541 | constructorParams.RMIN = 0; |
---|
542 | constructorParams.RMAX = 0; |
---|
543 | |
---|
544 | if( num_z_planes > 0 ) |
---|
545 | { |
---|
546 | constructorParams.z_values = new G4double[num_z_planes]; |
---|
547 | constructorParams.RMIN = new G4double[num_z_planes]; |
---|
548 | constructorParams.RMAX = new G4double[num_z_planes]; |
---|
549 | for( G4int idx = 0; idx < num_z_planes; idx++ ) |
---|
550 | { |
---|
551 | constructorParams.z_values[idx] = z_values[idx]; |
---|
552 | constructorParams.RMIN[idx] = RMIN[idx]; |
---|
553 | constructorParams.RMAX[idx] = RMAX[idx]; |
---|
554 | } |
---|
555 | } |
---|
556 | |
---|
557 | active=1; |
---|
558 | Initialize(); |
---|
559 | } |
---|
560 | |
---|
561 | G4BREPSolidPCone::G4BREPSolidPCone( __void__& a ) |
---|
562 | : G4BREPSolid(a) |
---|
563 | { |
---|
564 | } |
---|
565 | |
---|
566 | G4BREPSolidPCone::~G4BREPSolidPCone() |
---|
567 | { |
---|
568 | if( constructorParams.num_z_planes > 0 ) |
---|
569 | { |
---|
570 | delete [] constructorParams.z_values; |
---|
571 | delete [] constructorParams.RMIN; |
---|
572 | delete [] constructorParams.RMAX; |
---|
573 | } |
---|
574 | } |
---|
575 | |
---|
576 | void G4BREPSolidPCone::Initialize() |
---|
577 | { |
---|
578 | // Computes the bounding box for solids and surfaces |
---|
579 | // Converts concave planes to convex |
---|
580 | // |
---|
581 | ShortestDistance=1000000; |
---|
582 | CheckSurfaceNormals(); |
---|
583 | |
---|
584 | if(!Box || !AxisBox) |
---|
585 | IsConvex(); |
---|
586 | |
---|
587 | CalcBBoxes(); |
---|
588 | } |
---|
589 | |
---|
590 | EInside G4BREPSolidPCone::Inside(register const G4ThreeVector& Pt) const |
---|
591 | { |
---|
592 | // This function find if the point Pt is inside, |
---|
593 | // outside or on the surface of the solid |
---|
594 | |
---|
595 | G4Vector3D v(1, 0, 0.01); |
---|
596 | //G4Vector3D v(1, 0, 0); // This will miss the planar surface perp. to Z axis |
---|
597 | //G4Vector3D v(0, 0, 1); // This works, however considered as hack not a fix |
---|
598 | G4Vector3D Pttmp(Pt); |
---|
599 | G4Vector3D Vtmp(v); |
---|
600 | G4Ray r(Pttmp, Vtmp); |
---|
601 | |
---|
602 | // Check if point is inside the PCone bounding box |
---|
603 | // |
---|
604 | if( !GetBBox()->Inside(Pttmp) ) |
---|
605 | { |
---|
606 | return kOutside; |
---|
607 | } |
---|
608 | |
---|
609 | // Set the surfaces to active again |
---|
610 | // |
---|
611 | Reset(); |
---|
612 | |
---|
613 | // Test if the bounding box of each surface is intersected |
---|
614 | // by the ray. If not, the surface is deactivated. |
---|
615 | // |
---|
616 | TestSurfaceBBoxes(r); |
---|
617 | |
---|
618 | G4double dist = kInfinity; |
---|
619 | G4bool isIntersected = false; |
---|
620 | G4int WhichSurface = 0; |
---|
621 | |
---|
622 | const G4double sqrHalfTolerance = kCarTolerance*kCarTolerance*0.25; |
---|
623 | |
---|
624 | // Chech if the point is on the surface, otherwise |
---|
625 | // find the nearest intersected suface. If there are not intersections the |
---|
626 | // point is outside |
---|
627 | |
---|
628 | for(G4int a=0; a < nb_of_surfaces; a++) |
---|
629 | { |
---|
630 | G4Surface* surface = SurfaceVec[a]; |
---|
631 | |
---|
632 | if( surface->IsActive() ) |
---|
633 | { |
---|
634 | G4double hownear = surface->HowNear(Pt); |
---|
635 | |
---|
636 | if( std::fabs( hownear ) < sqrHalfTolerance ) |
---|
637 | { |
---|
638 | return kSurface; |
---|
639 | } |
---|
640 | |
---|
641 | if( surface->Intersect(r) ) |
---|
642 | { |
---|
643 | isIntersected = true; |
---|
644 | hownear = surface->GetDistance(); |
---|
645 | |
---|
646 | if ( std::fabs( hownear ) < dist ) |
---|
647 | { |
---|
648 | dist = hownear; |
---|
649 | WhichSurface = a; |
---|
650 | } |
---|
651 | } |
---|
652 | } |
---|
653 | } |
---|
654 | |
---|
655 | if ( !isIntersected ) |
---|
656 | { |
---|
657 | return kOutside; |
---|
658 | } |
---|
659 | |
---|
660 | // Find the point of intersection on the surface and the normal |
---|
661 | // !!!! be carefull the distance is std::sqrt(dist) !!!! |
---|
662 | |
---|
663 | dist = std::sqrt( dist ); |
---|
664 | G4Vector3D IntersectionPoint = Pttmp + dist*Vtmp; |
---|
665 | G4Vector3D Normal = |
---|
666 | SurfaceVec[WhichSurface]->SurfaceNormal( IntersectionPoint ); |
---|
667 | |
---|
668 | G4double dot = Normal*Vtmp; |
---|
669 | |
---|
670 | return ( (dot > 0) ? kInside : kOutside ); |
---|
671 | } |
---|
672 | |
---|
673 | G4ThreeVector |
---|
674 | G4BREPSolidPCone::SurfaceNormal(const G4ThreeVector& Pt) const |
---|
675 | { |
---|
676 | // This function calculates the normal of the closest surface |
---|
677 | // at a point on the surface |
---|
678 | // Note : the sense of the normal depends on the sense of the surface |
---|
679 | |
---|
680 | G4int isurf; |
---|
681 | G4bool normflag = false; |
---|
682 | |
---|
683 | const G4double sqrHalfTolerance = kCarTolerance*kCarTolerance*0.25; |
---|
684 | |
---|
685 | // Determine if the point is on the surface |
---|
686 | // |
---|
687 | G4double minDist = kInfinity; |
---|
688 | G4int normSurface = 0; |
---|
689 | for(isurf = 0; isurf < nb_of_surfaces; isurf++) |
---|
690 | { |
---|
691 | G4double dist = std::fabs(SurfaceVec[isurf]->HowNear(Pt)); |
---|
692 | if( minDist > dist ) |
---|
693 | { |
---|
694 | minDist = dist; |
---|
695 | normSurface = isurf; |
---|
696 | } |
---|
697 | if( dist < sqrHalfTolerance) |
---|
698 | { |
---|
699 | // the point is on this surface |
---|
700 | // |
---|
701 | normflag = true; |
---|
702 | break; |
---|
703 | } |
---|
704 | } |
---|
705 | |
---|
706 | // Calculate the normal at this point, if the point is on the |
---|
707 | // surface, otherwise compute the normal to the closest surface |
---|
708 | // |
---|
709 | if ( normflag ) // point on surface |
---|
710 | { |
---|
711 | G4ThreeVector norm = SurfaceVec[isurf]->SurfaceNormal(Pt); |
---|
712 | return norm.unit(); |
---|
713 | } |
---|
714 | else // point not on surface |
---|
715 | { |
---|
716 | G4Surface* nSurface = SurfaceVec[normSurface]; |
---|
717 | G4ThreeVector hitPt = nSurface->GetClosestHit(); |
---|
718 | G4ThreeVector hitNorm = nSurface->SurfaceNormal(hitPt); |
---|
719 | return hitNorm.unit(); |
---|
720 | } |
---|
721 | } |
---|
722 | |
---|
723 | G4double G4BREPSolidPCone::DistanceToIn(const G4ThreeVector& Pt) const |
---|
724 | { |
---|
725 | // Calculates the shortest distance ("safety") from a point |
---|
726 | // outside the solid to any boundary of this solid. |
---|
727 | // Return 0 if the point is already inside. |
---|
728 | |
---|
729 | G4double *dists = new G4double[nb_of_surfaces]; |
---|
730 | G4int a; |
---|
731 | |
---|
732 | // Set the surfaces to active again |
---|
733 | // |
---|
734 | Reset(); |
---|
735 | |
---|
736 | // Compute the shortest distance of the point to each surfaces |
---|
737 | // Be careful : it's a signed value |
---|
738 | // |
---|
739 | for(a=0; a< nb_of_surfaces; a++) |
---|
740 | dists[a] = SurfaceVec[a]->HowNear(Pt); |
---|
741 | |
---|
742 | G4double Dist = kInfinity; |
---|
743 | |
---|
744 | // if dists[] is positive, the point is outside |
---|
745 | // so take the shortest of the shortest positive distances |
---|
746 | // dists[] can be equal to 0 : point on a surface |
---|
747 | // ( Problem with the G4FPlane : there is no inside and no outside... |
---|
748 | // So, to test if the point is inside to return 0, utilize the Inside |
---|
749 | // function. But I don`t know if it is really needed because dToIn is |
---|
750 | // called only if the point is outside ) |
---|
751 | // |
---|
752 | for(a = 0; a < nb_of_surfaces; a++) |
---|
753 | if( std::fabs(Dist) > std::fabs(dists[a]) ) |
---|
754 | //if( dists[a] >= 0) |
---|
755 | Dist = dists[a]; |
---|
756 | |
---|
757 | delete[] dists; |
---|
758 | |
---|
759 | if(Dist == kInfinity) |
---|
760 | // the point is inside the solid or on a surface |
---|
761 | // |
---|
762 | return 0; |
---|
763 | else |
---|
764 | return std::fabs(Dist); |
---|
765 | } |
---|
766 | |
---|
767 | G4double G4BREPSolidPCone::DistanceToIn(register const G4ThreeVector& Pt, |
---|
768 | register const G4ThreeVector& V) const |
---|
769 | { |
---|
770 | // Calculates the distance from a point outside the solid |
---|
771 | // to the solid`s boundary along a specified direction vector. |
---|
772 | // |
---|
773 | // Note : Intersections with boundaries less than the |
---|
774 | // tolerance must be ignored if the direction |
---|
775 | // is away from the boundary |
---|
776 | |
---|
777 | G4int a; |
---|
778 | |
---|
779 | // Set the surfaces to active again |
---|
780 | // |
---|
781 | Reset(); |
---|
782 | |
---|
783 | G4double sqrHalfTolerance = kCarTolerance*kCarTolerance*0.25; |
---|
784 | G4Vector3D Pttmp(Pt); |
---|
785 | G4Vector3D Vtmp(V); |
---|
786 | G4Ray r(Pttmp, Vtmp); |
---|
787 | |
---|
788 | // Test if the bounding box of each surface is intersected |
---|
789 | // by the ray. If not, the surface become deactive. |
---|
790 | // |
---|
791 | TestSurfaceBBoxes(r); |
---|
792 | |
---|
793 | ShortestDistance = kInfinity; |
---|
794 | |
---|
795 | for(a=0; a< nb_of_surfaces; a++) |
---|
796 | { |
---|
797 | if(SurfaceVec[a]->IsActive()) |
---|
798 | { |
---|
799 | // test if the ray intersect the surface |
---|
800 | // |
---|
801 | G4Vector3D Norm = SurfaceVec[a]->SurfaceNormal(Pttmp); |
---|
802 | G4double hownear = SurfaceVec[a]->HowNear(Pt); |
---|
803 | |
---|
804 | if( (Norm * Vtmp) < 0 && std::fabs(hownear) < sqrHalfTolerance ) |
---|
805 | return 0; |
---|
806 | |
---|
807 | if( (SurfaceVec[a]->Intersect(r)) ) |
---|
808 | { |
---|
809 | // if more than 1 surface is intersected, |
---|
810 | // take the nearest one |
---|
811 | G4double distance = SurfaceVec[a]->GetDistance(); |
---|
812 | |
---|
813 | if( distance < ShortestDistance ) |
---|
814 | { |
---|
815 | if( distance > sqrHalfTolerance ) |
---|
816 | { |
---|
817 | ShortestDistance = distance; |
---|
818 | } |
---|
819 | } |
---|
820 | } |
---|
821 | } |
---|
822 | } |
---|
823 | |
---|
824 | // Be careful ! |
---|
825 | // SurfaceVec->Distance is in fact the squared distance |
---|
826 | // |
---|
827 | if(ShortestDistance != kInfinity) |
---|
828 | return( std::sqrt(ShortestDistance) ); |
---|
829 | else |
---|
830 | // no intersection |
---|
831 | // |
---|
832 | return kInfinity; |
---|
833 | } |
---|
834 | |
---|
835 | G4double G4BREPSolidPCone::DistanceToOut(register const G4ThreeVector& Pt, |
---|
836 | register const G4ThreeVector& V, |
---|
837 | const G4bool, |
---|
838 | G4bool *validNorm, |
---|
839 | G4ThreeVector *) const |
---|
840 | { |
---|
841 | // Calculates the distance from a point inside the solid |
---|
842 | // to the solid`s boundary along a specified direction vector. |
---|
843 | // Return 0 if the point is already outside. |
---|
844 | // |
---|
845 | // Note : If the shortest distance to a boundary is less |
---|
846 | // than the tolerance, it is ignored. This allows |
---|
847 | // for a point within a tolerant boundary to leave |
---|
848 | // immediately |
---|
849 | |
---|
850 | // Set the surfaces to active again |
---|
851 | // |
---|
852 | Reset(); |
---|
853 | |
---|
854 | const G4double sqrHalfTolerance = kCarTolerance*kCarTolerance*0.25; |
---|
855 | G4Vector3D Ptv = G4Vector3D( Pt ); |
---|
856 | G4int a; |
---|
857 | |
---|
858 | // I don`t understand this line |
---|
859 | // |
---|
860 | if(validNorm) |
---|
861 | *validNorm=false; |
---|
862 | |
---|
863 | G4Vector3D Pttmp(Pt); |
---|
864 | G4Vector3D Vtmp(V); |
---|
865 | |
---|
866 | G4Ray r(Pttmp, Vtmp); |
---|
867 | |
---|
868 | // Test if the bounding box of each surface is intersected |
---|
869 | // by the ray. If not, the surface become deactive. |
---|
870 | // |
---|
871 | TestSurfaceBBoxes(r); |
---|
872 | |
---|
873 | ShortestDistance = kInfinity; |
---|
874 | |
---|
875 | for(a=0; a< nb_of_surfaces; a++) |
---|
876 | { |
---|
877 | if( SurfaceVec[a]->IsActive() ) |
---|
878 | { |
---|
879 | G4Vector3D Norm = SurfaceVec[a]->SurfaceNormal(Pttmp); |
---|
880 | G4double hownear = SurfaceVec[a]->HowNear(Pt); |
---|
881 | |
---|
882 | if( (Norm * Vtmp) > 0 && std::fabs( hownear ) < sqrHalfTolerance ) |
---|
883 | return 0; |
---|
884 | |
---|
885 | // test if the ray intersect the surface |
---|
886 | // |
---|
887 | if( SurfaceVec[a]->Intersect(r) ) |
---|
888 | { |
---|
889 | // if more than 1 surface is intersected, |
---|
890 | // take the nearest one |
---|
891 | // |
---|
892 | G4double distance = SurfaceVec[a]->GetDistance(); |
---|
893 | |
---|
894 | if( distance < ShortestDistance ) |
---|
895 | { |
---|
896 | if( distance > sqrHalfTolerance ) |
---|
897 | { |
---|
898 | ShortestDistance = distance; |
---|
899 | } |
---|
900 | } |
---|
901 | } |
---|
902 | } |
---|
903 | } |
---|
904 | |
---|
905 | // Be careful ! |
---|
906 | // SurfaceVec->Distance is in fact the squared distance |
---|
907 | // |
---|
908 | if(ShortestDistance != kInfinity) |
---|
909 | return std::sqrt(ShortestDistance); |
---|
910 | else |
---|
911 | // if no intersection is founded, the point is outside |
---|
912 | // |
---|
913 | return 0; |
---|
914 | } |
---|
915 | |
---|
916 | G4double G4BREPSolidPCone::DistanceToOut(const G4ThreeVector& Pt) const |
---|
917 | { |
---|
918 | // Calculates the shortest distance ("safety") from a point |
---|
919 | // inside the solid to any boundary of this solid. |
---|
920 | // Return 0 if the point is already outside. |
---|
921 | |
---|
922 | G4double *dists = new G4double[nb_of_surfaces]; |
---|
923 | G4int a; |
---|
924 | |
---|
925 | // Set the surfaces to active again |
---|
926 | Reset(); |
---|
927 | |
---|
928 | // calcul of the shortest distance of the point to each surfaces |
---|
929 | // Be careful : it's a signed value |
---|
930 | // |
---|
931 | for(a=0; a< nb_of_surfaces; a++) |
---|
932 | dists[a] = SurfaceVec[a]->HowNear(Pt); |
---|
933 | |
---|
934 | G4double Dist = kInfinity; |
---|
935 | |
---|
936 | // if dists[] is negative, the point is inside |
---|
937 | // so take the shortest of the shortest negative distances |
---|
938 | // dists[] can be equal to 0 : point on a surface |
---|
939 | // ( Problem with the G4FPlane : there is no inside and no outside... |
---|
940 | // So, to test if the point is outside to return 0, utilize the Inside |
---|
941 | // function. But I don`t know if it is really needed because dToOut is |
---|
942 | // called only if the point is inside ) |
---|
943 | |
---|
944 | for(a = 0; a < nb_of_surfaces; a++) |
---|
945 | if( std::fabs(Dist) > std::fabs(dists[a]) ) |
---|
946 | //if( dists[a] <= 0) |
---|
947 | Dist = dists[a]; |
---|
948 | |
---|
949 | delete[] dists; |
---|
950 | |
---|
951 | if(Dist == kInfinity) |
---|
952 | // the point is ouside the solid or on a surface |
---|
953 | // |
---|
954 | return 0; |
---|
955 | else |
---|
956 | return std::fabs(Dist); |
---|
957 | } |
---|
958 | |
---|
959 | std::ostream& G4BREPSolidPCone::StreamInfo(std::ostream& os) const |
---|
960 | { |
---|
961 | // Streams solid contents to output stream. |
---|
962 | |
---|
963 | G4BREPSolid::StreamInfo( os ) |
---|
964 | << "\n start_angle: " << constructorParams.start_angle |
---|
965 | << "\n opening_angle: " << constructorParams.opening_angle |
---|
966 | << "\n num_z_planes: " << constructorParams.num_z_planes |
---|
967 | << "\n z_start: " << constructorParams.z_start |
---|
968 | << "\n z_values: "; |
---|
969 | G4int idx; |
---|
970 | for( idx = 0; idx < constructorParams.num_z_planes; idx++ ) |
---|
971 | { |
---|
972 | os << constructorParams.z_values[idx] << " "; |
---|
973 | } |
---|
974 | os << "\n RMIN: "; |
---|
975 | for( idx = 0; idx < constructorParams.num_z_planes; idx++ ) |
---|
976 | { |
---|
977 | os << constructorParams.RMIN[idx] << " "; |
---|
978 | } |
---|
979 | os << "\n RMAX: "; |
---|
980 | for( idx = 0; idx < constructorParams.num_z_planes; idx++ ) |
---|
981 | { |
---|
982 | os << constructorParams.RMAX[idx] << " "; |
---|
983 | } |
---|
984 | os << "\n-----------------------------------------------------------\n"; |
---|
985 | |
---|
986 | return os; |
---|
987 | } |
---|
988 | |
---|
989 | void G4BREPSolidPCone::Reset() const |
---|
990 | { |
---|
991 | Active(1); |
---|
992 | ((G4BREPSolidPCone*)this)->intersectionDistance=kInfinity; |
---|
993 | StartInside(0); |
---|
994 | for(register int a=0;a<nb_of_surfaces;a++) |
---|
995 | SurfaceVec[a]->Reset(); |
---|
996 | ShortestDistance = kInfinity; |
---|
997 | } |
---|
998 | |
---|
999 | G4Surface* |
---|
1000 | G4BREPSolidPCone::ComputePlanarSurface( G4double r1, |
---|
1001 | G4double r2, |
---|
1002 | G4ThreeVector& origin, |
---|
1003 | G4ThreeVector& planeAxis, |
---|
1004 | G4ThreeVector& planeDirection, |
---|
1005 | G4int surfSense ) |
---|
1006 | { |
---|
1007 | // The planar surface to return |
---|
1008 | G4Surface* planarFace = 0; |
---|
1009 | |
---|
1010 | G4CurveVector cv1; |
---|
1011 | G4CircularCurve *tmp1, *tmp2; |
---|
1012 | |
---|
1013 | // Create plane surface |
---|
1014 | G4Point3D ArcStart1 = G4Point3D( origin + (r1 * planeDirection) ); |
---|
1015 | G4Point3D ArcStart2 = G4Point3D( origin + (r2 * planeDirection) ); |
---|
1016 | |
---|
1017 | if(r1 != 0) |
---|
1018 | { |
---|
1019 | tmp1 = new G4CircularCurve; |
---|
1020 | tmp1->Init(G4Axis2Placement3D( planeDirection, planeAxis, origin), r1); |
---|
1021 | tmp1->SetBounds(ArcStart1, ArcStart1); |
---|
1022 | |
---|
1023 | if( surfSense ) |
---|
1024 | tmp1->SetSameSense(1); |
---|
1025 | else |
---|
1026 | tmp1->SetSameSense(0); |
---|
1027 | |
---|
1028 | cv1.push_back(tmp1); |
---|
1029 | } |
---|
1030 | |
---|
1031 | if(r2 != 0) |
---|
1032 | { |
---|
1033 | tmp2 = new G4CircularCurve; |
---|
1034 | tmp2->Init(G4Axis2Placement3D( planeDirection, planeAxis, origin), r2); |
---|
1035 | tmp2->SetBounds(ArcStart2, ArcStart2); |
---|
1036 | |
---|
1037 | if( surfSense ) |
---|
1038 | tmp2->SetSameSense(0); |
---|
1039 | else |
---|
1040 | tmp2->SetSameSense(1); |
---|
1041 | |
---|
1042 | cv1.push_back(tmp2); |
---|
1043 | } |
---|
1044 | |
---|
1045 | planarFace = new G4FPlane( planeDirection, planeAxis, origin, surfSense ); |
---|
1046 | |
---|
1047 | planarFace->SetBoundaries(&cv1); |
---|
1048 | |
---|
1049 | return planarFace; |
---|
1050 | } |
---|
1051 | |
---|
1052 | // In graphics_reps: |
---|
1053 | |
---|
1054 | #include "G4Polyhedron.hh" |
---|
1055 | |
---|
1056 | G4Polyhedron* G4BREPSolidPCone::CreatePolyhedron() const |
---|
1057 | { |
---|
1058 | return new G4PolyhedronPcon( constructorParams.start_angle, |
---|
1059 | constructorParams.opening_angle, |
---|
1060 | constructorParams.num_z_planes, |
---|
1061 | constructorParams.z_values, |
---|
1062 | constructorParams.RMIN, |
---|
1063 | constructorParams.RMAX ); |
---|
1064 | } |
---|