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

Last change on this file since 1321 was 1228, checked in by garnier, 16 years ago

update geant4.9.3 tag
File size: 11.5 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.11 2007/11/16 09:39:14 gcosmo Exp $
28// GEANT4 tag $Name: geant4-09-03 $
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 // G4cout<<"Entering Find Some Points vSearch="<<vSearch<<" p="<<p<<G4endl;
278 //
279 // Look for nearest intersection point in the specified
280 // direction and return if there isn't one
281 //
282 G4double dist;
283 switch(solid->Inside(p)) {
284 case kInside:
285 dist = solid->DistanceToOut(p,vSearch);
286 // G4cout<<"Inside DistToOut="<<dist<<G4endl;
287 if (dist >= kInfinity) {
288 logger->SolidProblem( solid,
289 "DistanceToOut(p,v) = kInfinity for point inside", p );
290 return;
291 }
292
293 s += sign*dist;
294 entering = false;
295 break;
296 case kOutside:
297 dist = solid->DistanceToIn(p,vSearch);
298 //G4cout<<"Outside DistToIn="<<dist<<G4endl;
299 if (dist >= kInfinity) return;
300
301 s += sign*dist;
302 entering = true;
303 break;
304 case kSurface:
305 vSurfN=vSearch.dot(solid->SurfaceNormal(p));
306 if(std::abs(vSurfN)<kCarTolerance)vSurfN=0;
307 entering = (vSurfN < 0);
308 //G4cout<<"Surface SurfN="<<solid->SurfaceNormal(p)<<" v.dotN="<<vSurfN<<" entering="<<entering<<G4endl;
309 break;
310 default:
311 logger->SolidProblem( solid,
312 "Inside returns illegal enumerated value", p );
313 return;
314 }
315
316 //
317 // Loop
318 //
319 // nzero = the number of consecutive zeros
320 //
321 G4int nzero=0;
322
323 for(;;) {
324 //
325 // Locate point
326 //
327 p = p0 + s*v;
328
329 if (nzero > 2) {
330 //
331 // Oops. In a loop. Probably along a spherical or cylindrical surface.
332 // Let's give the tool a little help with a push
333 //
334 G4double push = 1E-6;
335 s += sign*push;
336 for(;;) {
337 p = p0 + s*v;
338 EInside inside = solid->Inside(p);
339 if (inside == kInside) {
340 entering = true;
341 break;
342 }
343 else if (inside == kOutside) {
344 entering = false;
345 break;
346 }
347
348 push = 2*push;
349 if (push > 0.1) {
350 logger->SolidProblem( solid,
351 "Push fails to fix geometry inconsistency", p );
352 return;
353 }
354 s += sign*push;
355 }
356 }
357 else {
358
359 //
360 // Record point
361 //
362 points.push_back( G4GeomTestPoint( p, s, entering==forward ) );
363 //G4cout<<"Add point p"<<p<<" s="<<s<<" entering="<<entering<<G4endl;
364 }
365
366 //
367 // Find next intersection
368 //
369 if (entering) {
370 dist = solid->DistanceToOut(p,vSearch);
371 //G4cout<<"if entering distToOut="<<dist<<G4endl;
372 if (dist >= kInfinity) {
373 logger->SolidProblem( solid,
374 "DistanceToOut(p,v) = kInfinity for point inside", p );
375 return;
376 }
377
378 if ( (dist > kCarTolerance)
379 && (solid->Inside(p + (dist*0.999)*vSearch) == kOutside) ) {
380 //
381 // This shouldn't have happened
382 //
383 if (solid->Inside(p) == kOutside)
384 logger->SolidProblem(solid,
385 "Entering point is outside (possible roundoff error)",p);
386 else
387 logger->SolidProblem(solid,
388 "DistanceToOut(p,v) brings trajectory well outside solid",p);
389 return;
390 }
391
392 if(std::abs(dist)<=kCarTolerance){
393 G4double push = 1E-6;
394 s += sign*push;
395 p = p0 + s*v;
396 EInside inside = solid->Inside(p);
397 if (inside == kOutside) {
398 entering = false;
399 break;
400 }
401 }
402
403
404 entering = false;
405 }
406 else {
407 dist = solid->DistanceToIn(p,vSearch);
408 //G4cout<<"if exiting distToIn="<<dist<<G4endl;
409 if (dist >= kInfinity) return;
410
411 if ( (dist > kCarTolerance)
412 && (solid->Inside(p + (dist*0.999)*vSearch) == kInside) ) {
413 //
414 // This shouldn't have happened
415 //
416 if (solid->Inside(p) == kInside)
417 logger->SolidProblem(solid,
418 "Exiting point is inside (possible roundoff error)", p);
419 else
420 logger->SolidProblem(solid,
421 "DistanceToIn(p,v) brings trajectory well inside solid", p);
422 return;
423 }
424
425 entering = true;
426 }
427
428 //
429 // Update s
430 //
431 if (dist <= 0) {
432 nzero++;
433 }
434 else {
435 nzero=0;
436 s += sign*dist;
437 }
438 }
439}
Note: See TracBrowser for help on using the repository browser.