00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037 #include "G4GeomTestSegment.hh"
00038
00039 #include "G4VSolid.hh"
00040 #include "G4GeomTestLogger.hh"
00041 #include "G4GeometryTolerance.hh"
00042
00043
00044
00045
00046 G4GeomTestSegment::G4GeomTestSegment( const G4VSolid *theSolid,
00047 const G4ThreeVector &theP,
00048 const G4ThreeVector &theV,
00049 G4GeomTestLogger *logger )
00050 : solid(theSolid),
00051 p0(theP),
00052 v(theV)
00053 {
00054 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
00055 FindPoints(logger);
00056 }
00057
00058
00059
00060
00061
00062 const G4VSolid * G4GeomTestSegment::GetSolid() const { return solid; }
00063
00064 const G4ThreeVector &G4GeomTestSegment::GetP() const { return p0; }
00065
00066 const G4ThreeVector &G4GeomTestSegment::GetV() const { return v; }
00067
00068
00069
00070
00071
00072 G4int G4GeomTestSegment::GetNumberPoints() const
00073 {
00074 return points.size();
00075 }
00076
00077
00078
00079
00080
00081 const G4GeomTestPoint &G4GeomTestSegment::GetPoint( G4int i ) const
00082 {
00083 return points[i];
00084 }
00085
00086
00087
00088
00089
00090
00091
00092 void G4GeomTestSegment::FindPoints( G4GeomTestLogger *logger )
00093 {
00094 FindSomePoints( logger, true );
00095 FindSomePoints( logger, false );
00096
00097 PatchInconsistencies( logger );
00098 }
00099
00100
00101
00102
00103
00104
00105
00106
00107 void G4GeomTestSegment::PatchInconsistencies( G4GeomTestLogger *logger )
00108 {
00109 if (points.size() == 0) return;
00110
00111
00112
00113
00114 std::sort( points.begin(), points.end() );
00115
00116
00117
00118
00119 std::vector<G4GeomTestPoint>::iterator curr = points.begin();
00120 do {
00121
00122 std::vector<G4GeomTestPoint>::iterator next = curr + 1;
00123
00124
00125
00126
00127 while (next != points.end() &&
00128 next->GetDistance()-curr->GetDistance() < kCarTolerance) {
00129
00130
00131
00132 if (next->Entering() != curr->Entering()) {
00133 curr = next + 1;
00134 next = curr + 1;
00135
00136 break;
00137 }
00138
00139
00140
00141
00142 next = points.erase(next);
00143 curr = next - 1;
00144
00145 }
00146
00147 if (curr == points.end()) break;
00148
00149
00150
00151
00152
00153 if (!curr->Entering()) {
00154
00155
00156
00157
00158 G4double ds = curr->GetDistance();
00159 G4ThreeVector p = p0 + ds*v;
00160
00161 G4ThreeVector p1 = p - 10*kCarTolerance*v;
00162
00163 if (solid->Inside(p1) == kOutside||(solid->Inside(p1)== kSurface)) {
00164
00165
00166
00167
00168 curr = points.insert(curr, G4GeomTestPoint( p, ds, true ) );
00169 ++curr;
00170 break;
00171 }
00172
00173
00174
00175
00176 if (curr != points.begin()) {
00177 std::vector<G4GeomTestPoint>::iterator prev = curr - 1;
00178
00179 ds = prev->GetDistance();
00180 p = p0 + ds*v;
00181
00182 p1 = p + 10*kCarTolerance*v;
00183 if ((solid->Inside(p1) == kOutside)||(solid->Inside(p1)== kSurface)) {
00184
00185
00186
00187
00188
00189 curr = points.insert(curr, G4GeomTestPoint( p, ds, true ) );
00190 ++curr;
00191 break;
00192 }
00193 }
00194
00195
00196
00197
00198 logger->SolidProblem( solid, "Spurious exiting intersection point", p );
00199 curr = points.erase(curr);
00200 break;
00201 }
00202
00203
00204
00205
00206
00207 if (next == points.end() || next->Entering() ) {
00208
00209
00210
00211 G4double ds = curr->GetDistance();
00212 G4ThreeVector p = p0 + ds*v;
00213 G4ThreeVector p1 = p + 10*kCarTolerance*v;
00214
00215 if (solid->Inside(p1) == kOutside) {
00216
00217
00218
00219
00220 curr = points.insert(next, G4GeomTestPoint( p, ds, false ) );
00221 break;
00222 }
00223
00224 if (next != points.end()) {
00225
00226
00227
00228 ds = next->GetDistance();
00229 p = p0 + ds*v;
00230 p1 = p - 10*kCarTolerance*v;
00231 if (solid->Inside(p1) == kOutside) {
00232
00233
00234
00235
00236 curr = points.insert(next, G4GeomTestPoint( p, ds, false ) );
00237 break;
00238 }
00239 }
00240
00241
00242
00243
00244 logger->SolidProblem( solid, "Spurious entering intersection point", p );
00245 curr = points.erase(curr);
00246 }
00247
00248 if(curr!=points.end()){curr = next + 1;}
00249
00250 } while( curr != points.end() );
00251
00252
00253
00254
00255 if (points.size()&1)
00256 logger->SolidProblem( solid,
00257 "Solid has odd number of intersection points", p0 );
00258
00259 return;
00260 }
00261
00262
00263
00264
00265
00266 void G4GeomTestSegment::FindSomePoints( G4GeomTestLogger *logger,
00267 G4bool forward )
00268 {
00269 G4double sign = forward ? +1 : -1;
00270
00271 G4ThreeVector p(p0);
00272 G4ThreeVector vSearch(sign*v);
00273 G4double ds(0);
00274 G4bool entering;
00275 G4double vSurfN;
00276
00277
00278
00279
00280 G4double dist;
00281 switch(solid->Inside(p)) {
00282 case kInside:
00283 dist = solid->DistanceToOut(p,vSearch);
00284 if (dist >= kInfinity) {
00285 logger->SolidProblem( solid,
00286 "DistanceToOut(p,v) = kInfinity for point inside", p );
00287 return;
00288 }
00289 ds += sign*dist;
00290 entering = false;
00291 break;
00292 case kOutside:
00293 dist = solid->DistanceToIn(p,vSearch);
00294 if (dist >= kInfinity) return;
00295 ds += sign*dist;
00296 entering = true;
00297 break;
00298 case kSurface:
00299 vSurfN=v.dot(solid->SurfaceNormal(p));
00300 if(std::fabs(vSurfN)<kCarTolerance)vSurfN=0;
00301 entering = (vSurfN < 0);
00302 break;
00303 default:
00304 logger->SolidProblem( solid,
00305 "Inside returns illegal enumerated value", p );
00306 return;
00307 }
00308
00309
00310
00311
00312
00313
00314 G4int nzero=0;
00315
00316 for(;;) {
00317
00318
00319
00320 p = p0 + ds*v;
00321
00322 if (nzero > 2) {
00323
00324
00325
00326
00327 G4double push = 1E-6;
00328 ds += sign*push;
00329 for(;;) {
00330 p = p0 + ds*v;
00331 EInside inside = solid->Inside(p);
00332 if (inside == kInside) {
00333 entering = true;
00334 break;
00335 }
00336 else if (inside == kOutside) {
00337 entering = false;
00338 break;
00339 }
00340
00341 push = 2*push;
00342 if (push > 0.1) {
00343 logger->SolidProblem( solid,
00344 "Push fails to fix geometry inconsistency", p );
00345 return;
00346 }
00347 ds += sign*push;
00348 }
00349 }
00350 else {
00351
00352
00353
00354
00355 points.push_back( G4GeomTestPoint( p, ds, entering==forward ) );
00356
00357 }
00358
00359
00360
00361
00362 if (entering) {
00363 dist = solid->DistanceToOut(p,vSearch);
00364 if (dist >= kInfinity) {
00365 logger->SolidProblem( solid,
00366 "DistanceToOut(p,v) = kInfinity for point inside", p );
00367 return;
00368 }
00369
00370 if ( (dist > kCarTolerance)
00371 && (solid->Inside(p + (dist*0.999)*vSearch) == kOutside) ) {
00372
00373
00374
00375 if (solid->Inside(p) == kOutside)
00376 logger->SolidProblem(solid,
00377 "Entering point is outside (possible roundoff error)",p);
00378 else
00379 logger->SolidProblem(solid,
00380 "DistanceToOut(p,v) brings trajectory well outside solid",p);
00381 return;
00382 }
00383
00384 entering = false;
00385 }
00386 else {
00387 dist = solid->DistanceToIn(p,vSearch);
00388 if (dist >= kInfinity) return;
00389
00390 if ( (dist > kCarTolerance)
00391 && (solid->Inside(p + (dist*0.999)*vSearch) == kInside) ) {
00392
00393
00394
00395 if (solid->Inside(p) == kInside)
00396 logger->SolidProblem(solid,
00397 "Exiting point is inside (possible roundoff error)", p);
00398 else
00399 logger->SolidProblem(solid,
00400 "DistanceToIn(p,v) brings trajectory well inside solid", p);
00401 return;
00402 }
00403
00404 entering = true;
00405 }
00406
00407
00408
00409
00410 if (dist <= 0) {
00411 nzero++;
00412 }
00413 else {
00414 nzero=0;
00415 ds += sign*dist;
00416 }
00417 }
00418 }