source: trunk/source/geometry/navigation/src/G4GeomTestSegment.cc@ 1350

Last change on this file since 1350 was 1347, checked in by garnier, 15 years ago

geant4 tag 9.4
File size: 10.8 KB
Line 
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//
27// $Id: G4GeomTestSegment.cc,v 1.13 2010/08/20 09:03:54 gcosmo Exp $
28// GEANT4 tag $Name: geant4-09-04-ref-00 $
29//
30// --------------------------------------------------------------------
31// GEANT 4 class source file
32//
33// G4GeomTestSegment
34//
35// Author: D.C.Williams, UCSC (davidw@scipp.ucsc.edu)
36// --------------------------------------------------------------------
37
38#include "G4GeomTestSegment.hh"
39
40#include "G4VSolid.hh"
41#include "G4GeomTestLogger.hh"
42#include "G4GeometryTolerance.hh"
43
44//
45// Constructor
46//
47G4GeomTestSegment::G4GeomTestSegment( const G4VSolid *theSolid,
48 const G4ThreeVector &theP,
49 const G4ThreeVector &theV,
50 G4GeomTestLogger *logger )
51 : solid(theSolid),
52 p0(theP),
53 v(theV)
54{
55 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
56 FindPoints(logger);
57}
58
59
60//
61// Simple accessors
62//
63const G4VSolid * G4GeomTestSegment::GetSolid() const { return solid; }
64
65const G4ThreeVector &G4GeomTestSegment::GetP() const { return p0; }
66
67const G4ThreeVector &G4GeomTestSegment::GetV() const { return v; }
68
69
70//
71// Return number points
72//
73G4int G4GeomTestSegment::GetNumberPoints() const
74{
75 return points.size();
76}
77
78
79//
80// Return ith point
81//
82const G4GeomTestPoint &G4GeomTestSegment::GetPoint( G4int i ) const
83{
84 return points[i];
85}
86
87
88//
89// FindPoints
90//
91// Do the dirty work
92//
93void G4GeomTestSegment::FindPoints( G4GeomTestLogger *logger )
94{
95 FindSomePoints( logger, true );
96 FindSomePoints( logger, false );
97
98 PatchInconsistencies( logger );
99}
100
101
102//
103// PatchInconsistencies
104//
105// Search for inconsistancies in the hit list and patch
106// them up, if possible
107//
108void G4GeomTestSegment::PatchInconsistencies( G4GeomTestLogger *logger )
109{
110 if (points.size() == 0) return;
111
112 //
113 // Sort
114 //
115 std::sort( points.begin(), points.end() );
116
117 //
118 // Loop over entering/leaving pairs
119 //
120 std::vector<G4GeomTestPoint>::iterator curr = points.begin();
121 do {
122
123 std::vector<G4GeomTestPoint>::iterator next = curr + 1;
124
125 //
126 // Is the next point close by?
127 //
128 while (next != points.end() &&
129 next->GetDistance()-curr->GetDistance() < kCarTolerance) {
130 //
131 // Yup. If we have an entering/leaving pair, all is okay
132 //
133 if (next->Entering() != curr->Entering()) {
134 curr = next + 1;
135 next = curr + 1;
136
137 break;
138 }
139
140 //
141 // Otherwise, they are duplicate, and just delete one
142 //
143 next = points.erase(next);
144 curr = next - 1;
145
146 }
147
148 if (curr == points.end()) break;
149
150 //
151 // The next point should be entering...
152 //
153
154 if (!curr->Entering()) {
155 //
156 // Point is out of sequence:
157 // Test solid just before this point
158 //
159 G4double s = curr->GetDistance();
160 G4ThreeVector p = p0 + s*v;
161
162 G4ThreeVector p1 = p - 10*kCarTolerance*v;
163
164 if (solid->Inside(p1) == kOutside||(solid->Inside(p1)== kSurface)) {
165 //
166 // We are missing an entrance point near the current
167 // point. Add one.
168 //
169 curr = points.insert(curr, G4GeomTestPoint( p, s, true ) );
170 ++curr;
171 break;
172 }
173
174 //
175 // Test solid just after last point
176 //
177 if (curr != points.begin()) {
178 std::vector<G4GeomTestPoint>::iterator prev = curr - 1;
179
180 s = prev->GetDistance();
181 p = p0 + s*v;
182
183 p1 = p + 10*kCarTolerance*v;
184 if ((solid->Inside(p1) == kOutside)||(solid->Inside(p1)== kSurface)) {
185
186 //
187 // We are missing an entrance point near the previous
188 // point. Add one.
189 //
190 curr = points.insert(curr, G4GeomTestPoint( p, s, true ) );
191 ++curr;
192 break;
193 }
194 }
195
196 //
197 // Oh oh. No solution. Delete the current point and complain.
198 //
199 logger->SolidProblem( solid, "Spurious exiting intersection point", p );
200 curr = points.erase(curr);
201 break;
202 }
203
204 //
205 // The following point should be leaving
206 //
207
208 if (next == points.end() || next->Entering() ) {
209 //
210 // But is missing. Check immediately after this point.
211 //
212 G4double s = curr->GetDistance();
213 G4ThreeVector p = p0 + s*v;
214 G4ThreeVector p1 = p + 10*kCarTolerance*v;
215
216 if (solid->Inside(p1) == kOutside) {
217 //
218 // We are missing an exit point near the current
219 // point. Add one.
220 //
221 curr = points.insert(next, G4GeomTestPoint( p, s, false ) );
222 break;
223 }
224
225 if (next != points.end()) {
226 //
227 // Check just before next point
228 //
229 s = next->GetDistance();
230 p = p0 + s*v;
231 p1 = p - 10*kCarTolerance*v;
232 if (solid->Inside(p1) == kOutside) {
233 //
234 // We are missing an exit point before the next
235 // point. Add one.
236 //
237 curr = points.insert(next, G4GeomTestPoint( p, s, false ) );
238 break;
239 }
240 }
241
242 //
243 // Oh oh. Hanging entering point. Delete and complain.
244 //
245 logger->SolidProblem( solid, "Spurious entering intersection point", p );
246 curr = points.erase(curr);
247 }
248
249 if(curr!=points.end()){curr = next + 1;}
250
251 } while( curr != points.end() );
252
253 //
254 // Double check
255 //
256 if (points.size()&1)
257 logger->SolidProblem( solid,
258 "Solid has odd number of intersection points", p0 );
259
260 return;
261}
262
263
264//
265// Find intersections either in front or behind the point
266//
267void G4GeomTestSegment::FindSomePoints( G4GeomTestLogger *logger,
268 G4bool forward )
269{
270 G4double sign = forward ? +1 : -1;
271
272 G4ThreeVector p(p0);
273 G4ThreeVector vSearch(sign*v);
274 G4double s(0);
275 G4bool entering;
276 G4double vSurfN;
277
278 // Look for nearest intersection point in the specified
279 // direction and return if there isn't one
280 //
281 G4double dist;
282 switch(solid->Inside(p)) {
283 case kInside:
284 dist = solid->DistanceToOut(p,vSearch);
285 if (dist >= kInfinity) {
286 logger->SolidProblem( solid,
287 "DistanceToOut(p,v) = kInfinity for point inside", p );
288 return;
289 }
290 s += sign*dist;
291 entering = false;
292 break;
293 case kOutside:
294 dist = solid->DistanceToIn(p,vSearch);
295 if (dist >= kInfinity) return;
296 s += sign*dist;
297 entering = true;
298 break;
299 case kSurface:
300 vSurfN=v.dot(solid->SurfaceNormal(p));
301 if(std::fabs(vSurfN)<kCarTolerance)vSurfN=0;
302 entering = (vSurfN < 0);
303 break;
304 default:
305 logger->SolidProblem( solid,
306 "Inside returns illegal enumerated value", p );
307 return;
308 }
309
310 //
311 // Loop
312 //
313 // nzero = the number of consecutive zeros
314 //
315 G4int nzero=0;
316
317 for(;;) {
318 //
319 // Locate point
320 //
321 p = p0 + s*v;
322
323 if (nzero > 2) {
324 //
325 // Oops. In a loop. Probably along a spherical or cylindrical surface.
326 // Let's give the tool a little help with a push
327 //
328 G4double push = 1E-6;
329 s += sign*push;
330 for(;;) {
331 p = p0 + s*v;
332 EInside inside = solid->Inside(p);
333 if (inside == kInside) {
334 entering = true;
335 break;
336 }
337 else if (inside == kOutside) {
338 entering = false;
339 break;
340 }
341
342 push = 2*push;
343 if (push > 0.1) {
344 logger->SolidProblem( solid,
345 "Push fails to fix geometry inconsistency", p );
346 return;
347 }
348 s += sign*push;
349 }
350 }
351 else {
352
353 //
354 // Record point
355 //
356 points.push_back( G4GeomTestPoint( p, s, entering==forward ) );
357
358 }
359
360 //
361 // Find next intersection
362 //
363 if (entering) {
364 dist = solid->DistanceToOut(p,vSearch);
365 if (dist >= kInfinity) {
366 logger->SolidProblem( solid,
367 "DistanceToOut(p,v) = kInfinity for point inside", p );
368 return;
369 }
370
371 if ( (dist > kCarTolerance)
372 && (solid->Inside(p + (dist*0.999)*vSearch) == kOutside) ) {
373 //
374 // This shouldn't have happened
375 //
376 if (solid->Inside(p) == kOutside)
377 logger->SolidProblem(solid,
378 "Entering point is outside (possible roundoff error)",p);
379 else
380 logger->SolidProblem(solid,
381 "DistanceToOut(p,v) brings trajectory well outside solid",p);
382 return;
383 }
384
385 entering = false;
386 }
387 else {
388 dist = solid->DistanceToIn(p,vSearch);
389 if (dist >= kInfinity) return;
390
391 if ( (dist > kCarTolerance)
392 && (solid->Inside(p + (dist*0.999)*vSearch) == kInside) ) {
393 //
394 // This shouldn't have happened
395 //
396 if (solid->Inside(p) == kInside)
397 logger->SolidProblem(solid,
398 "Exiting point is inside (possible roundoff error)", p);
399 else
400 logger->SolidProblem(solid,
401 "DistanceToIn(p,v) brings trajectory well inside solid", p);
402 return;
403 }
404
405 entering = true;
406 }
407
408 //
409 // Update s
410 //
411 if (dist <= 0) {
412 nzero++;
413 }
414 else {
415 nzero=0;
416 s += sign*dist;
417 }
418 }
419}
Note: See TracBrowser for help on using the repository browser.