[1316] | 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 | // SBTrun |
---|
| 28 | // |
---|
| 29 | // Implementation of the batch solid test |
---|
| 30 | // |
---|
| 31 | #include "SBTrun.hh" |
---|
| 32 | |
---|
| 33 | #include "Randomize.hh" |
---|
| 34 | #include "G4VSolid.hh" |
---|
| 35 | |
---|
| 36 | #include "SBTVisManager.hh" |
---|
| 37 | #include "G4Polyline.hh" |
---|
| 38 | #include "G4Circle.hh" |
---|
| 39 | #include "G4Color.hh" |
---|
| 40 | #include "G4VisAttributes.hh" |
---|
| 41 | #include "G4GeometryTolerance.hh" |
---|
| 42 | #include "G4SolidStore.hh" |
---|
| 43 | |
---|
| 44 | #include <iomanip> |
---|
| 45 | #include <sstream> |
---|
| 46 | #include <time.h> |
---|
| 47 | |
---|
| 48 | #define DEBUG 0 |
---|
| 49 | |
---|
| 50 | // |
---|
| 51 | // Constructor |
---|
| 52 | // |
---|
| 53 | SBTrun::SBTrun() |
---|
| 54 | { |
---|
| 55 | // |
---|
| 56 | // Defaults |
---|
| 57 | // |
---|
| 58 | SetDefaults(); |
---|
| 59 | |
---|
| 60 | // |
---|
| 61 | // Zero error list |
---|
| 62 | // |
---|
| 63 | errorList = 0; |
---|
| 64 | CurrentSolid = ""; |
---|
| 65 | } |
---|
| 66 | |
---|
| 67 | |
---|
| 68 | // |
---|
| 69 | // Destructor |
---|
| 70 | // |
---|
| 71 | SBTrun::~SBTrun() |
---|
| 72 | { |
---|
| 73 | ClearErrors(); |
---|
| 74 | } |
---|
| 75 | |
---|
| 76 | |
---|
| 77 | // |
---|
| 78 | // SetDefaults |
---|
| 79 | // |
---|
| 80 | // Set defaults values |
---|
| 81 | // |
---|
| 82 | void SBTrun::SetDefaults() |
---|
| 83 | { |
---|
| 84 | target = G4ThreeVector( 0, 0, 0 ); |
---|
| 85 | widths = G4ThreeVector( 1*m, 1*m, 1*m ); |
---|
| 86 | grids = G4ThreeVector( 0, 0, 0 ); |
---|
| 87 | |
---|
| 88 | maxPoints = 10000; |
---|
| 89 | maxErrors = 100; |
---|
| 90 | } |
---|
| 91 | |
---|
| 92 | |
---|
| 93 | // |
---|
| 94 | // GetRandomPoint |
---|
| 95 | // |
---|
| 96 | // Return a random point in three dimensions using the current |
---|
| 97 | // specs |
---|
| 98 | // |
---|
| 99 | G4ThreeVector SBTrun::GetRandomPoint() const { |
---|
| 100 | G4double dx = widths.x()*GaussianRandom(10*m/widths.x()), |
---|
| 101 | dy = widths.y()*GaussianRandom(10*m/widths.y()), |
---|
| 102 | dz = widths.z()*GaussianRandom(10*m/widths.z()); |
---|
| 103 | |
---|
| 104 | |
---|
| 105 | if (grids.x() > 0) dx = grids.x()*G4int( dx/grids.x()+1 ); |
---|
| 106 | if (grids.y() > 0) dy = grids.y()*G4int( dy/grids.y()+1 ); |
---|
| 107 | if (grids.z() > 0) dz = grids.z()*G4int( dz/grids.z()+1 ); |
---|
| 108 | |
---|
| 109 | G4ThreeVector randvec( dx, dy, dz ); |
---|
| 110 | |
---|
| 111 | return target + randvec; |
---|
| 112 | } |
---|
| 113 | |
---|
| 114 | |
---|
| 115 | // |
---|
| 116 | // GaussianRandom |
---|
| 117 | // |
---|
| 118 | // Return Gaussian random number of unit width |
---|
| 119 | // |
---|
| 120 | // A classic, slow, but remarkably effective algorithm. Certainly good |
---|
| 121 | // enough for our purposes. |
---|
| 122 | // |
---|
| 123 | G4double SBTrun::GaussianRandom(const G4double cutoff) const { |
---|
| 124 | if (cutoff <= 0) G4Exception( "Illegal cutoff" ); |
---|
| 125 | |
---|
| 126 | G4double answer; |
---|
| 127 | do { |
---|
| 128 | answer = -3.0; |
---|
| 129 | for( G4int j = 0; j < 6; j++ ) answer += G4UniformRand(); |
---|
| 130 | answer *= std::sqrt(2.0); |
---|
| 131 | } while( std::fabs(answer) > cutoff ); |
---|
| 132 | |
---|
| 133 | return(answer); |
---|
| 134 | } |
---|
| 135 | |
---|
| 136 | |
---|
| 137 | // |
---|
| 138 | // RunTest |
---|
| 139 | // |
---|
| 140 | // Do your stuff! |
---|
| 141 | // |
---|
| 142 | void SBTrun::RunTest( const G4VSolid *testVolume, std::ostream &logger ) |
---|
| 143 | { |
---|
| 144 | // |
---|
| 145 | // Clear error list |
---|
| 146 | // |
---|
| 147 | |
---|
| 148 | ClearErrors(); |
---|
| 149 | |
---|
| 150 | // |
---|
| 151 | // Output test parameters |
---|
| 152 | // |
---|
| 153 | time_t now = time(0); |
---|
| 154 | time(&now); |
---|
| 155 | G4String dateTime(ctime(&now)); |
---|
| 156 | |
---|
| 157 | { |
---|
| 158 | /* |
---|
| 159 | MEDERNACH Emmanuel |
---|
| 160 | aug 2000 |
---|
| 161 | |
---|
| 162 | We try to retrieve the solid being tested |
---|
| 163 | to write in the log file |
---|
| 164 | */ |
---|
| 165 | |
---|
| 166 | G4SolidStore* SolidStore = G4SolidStore::GetInstance(); |
---|
| 167 | G4cout << ((*SolidStore)[0])->GetName() << G4endl ; |
---|
| 168 | } |
---|
| 169 | |
---|
| 170 | |
---|
| 171 | logger << "% SBT logged output " << dateTime; |
---|
| 172 | logger << "% " << CurrentSolid << G4endl; |
---|
| 173 | logger << "% target = " << target << G4endl; |
---|
| 174 | logger << "% widths = " << widths << G4endl; |
---|
| 175 | logger << "% grids = " << grids << G4endl; |
---|
| 176 | logger << "% maxPoints = " << maxPoints << G4endl; |
---|
| 177 | logger << "% maxErrors = " << maxErrors << G4endl; |
---|
| 178 | |
---|
| 179 | // |
---|
| 180 | // Setup lists |
---|
| 181 | // |
---|
| 182 | SBTrunPointList inside(100); |
---|
| 183 | SBTrunPointList outside(100); |
---|
| 184 | SBTrunPointList surface(100); |
---|
| 185 | |
---|
| 186 | // |
---|
| 187 | // Set iostream precision to 20 digits |
---|
| 188 | // |
---|
| 189 | logger << std::setprecision(20); |
---|
| 190 | |
---|
| 191 | // |
---|
| 192 | // Set clock |
---|
| 193 | // |
---|
| 194 | // clock_t start = clock(); |
---|
| 195 | |
---|
| 196 | // |
---|
| 197 | // Loop over points |
---|
| 198 | // |
---|
| 199 | G4int nIn = 0, nOut = 0, nSurf = 0; |
---|
| 200 | |
---|
| 201 | G4int nPoint = 0; |
---|
| 202 | G4int nError = 0; |
---|
| 203 | |
---|
| 204 | std::setprecision(20); |
---|
| 205 | |
---|
| 206 | for(;;) { |
---|
| 207 | // |
---|
| 208 | // Generate a random point |
---|
| 209 | // |
---|
| 210 | G4ThreeVector point = GetRandomPoint(); |
---|
| 211 | |
---|
| 212 | #if DEBUG |
---|
| 213 | G4cerr << std::setprecision(32); |
---|
| 214 | G4cerr << "Current point is : " << point* (1/cm) << G4endl ; |
---|
| 215 | #endif |
---|
| 216 | // |
---|
| 217 | // Catagorize, test, and store |
---|
| 218 | // |
---|
| 219 | EInside catagory = testVolume->Inside( point ); |
---|
| 220 | |
---|
| 221 | |
---|
| 222 | switch( catagory ) { |
---|
| 223 | case kOutside: |
---|
| 224 | //G4cerr << "kOutside " << G4endl; |
---|
| 225 | |
---|
| 226 | nOut++; |
---|
| 227 | TestOutsidePoint( testVolume, &nError, &inside, &outside, point, logger ); |
---|
| 228 | outside.AddPoint( point ); |
---|
| 229 | break; |
---|
| 230 | |
---|
| 231 | case kInside: |
---|
| 232 | //G4cerr << "kInside " << G4endl; |
---|
| 233 | |
---|
| 234 | nIn++; |
---|
| 235 | TestInsidePoint( testVolume, &nError, &outside, point, logger ); |
---|
| 236 | inside.AddPoint( point ); |
---|
| 237 | break; |
---|
| 238 | |
---|
| 239 | case kSurface: |
---|
| 240 | //G4cerr << "kSurface " << G4endl; |
---|
| 241 | |
---|
| 242 | nSurf++; |
---|
| 243 | // TestSurfacePoint( testVolume, &nError, point, logger ); |
---|
| 244 | surface.AddPoint( point ); |
---|
| 245 | break; |
---|
| 246 | } |
---|
| 247 | |
---|
| 248 | // |
---|
| 249 | // Return early if we have accumulated enough errors |
---|
| 250 | // |
---|
| 251 | if (nError >= maxErrors) { |
---|
| 252 | logger << "% End of test (maximum number errors) "; |
---|
| 253 | break; |
---|
| 254 | } |
---|
| 255 | |
---|
| 256 | if (++nPoint >= maxPoints) { |
---|
| 257 | logger << "% End of test (maximum number points) "; |
---|
| 258 | break; |
---|
| 259 | } |
---|
| 260 | |
---|
| 261 | } |
---|
| 262 | |
---|
| 263 | time(&now); |
---|
| 264 | G4String dateTime2(ctime(&now)); |
---|
| 265 | logger << dateTime2; |
---|
| 266 | |
---|
| 267 | logger << "% Statistics: points=" << nPoint << " errors=" << CountErrors() |
---|
| 268 | << " errors reported=" << nError << G4endl; |
---|
| 269 | |
---|
| 270 | logger << "% inside=" << nIn << " outside=" << nOut |
---|
| 271 | << " surface=" << nSurf << G4endl; |
---|
| 272 | logger << "% cpu time=" << clock()/CLOCKS_PER_SEC << G4endl; |
---|
| 273 | |
---|
| 274 | logger << "%(End of file)" << G4endl; |
---|
| 275 | } |
---|
| 276 | |
---|
| 277 | |
---|
| 278 | // |
---|
| 279 | // DrawError |
---|
| 280 | // |
---|
| 281 | // Recover previously logged error and display it |
---|
| 282 | // |
---|
| 283 | G4int SBTrun::DrawError( const G4VSolid *testVolume, std::istream &logger, |
---|
| 284 | const G4int errorIndex, SBTVisManager *visManager ) const |
---|
| 285 | { |
---|
| 286 | G4ThreeVector p, v; |
---|
| 287 | |
---|
| 288 | // Now required for drawing |
---|
| 289 | G4Transform3D objectTransformation; |
---|
| 290 | |
---|
| 291 | // |
---|
| 292 | // This draws the target solid |
---|
| 293 | // |
---|
| 294 | G4VisAttributes redStuff( G4Color(1,0,0) ); |
---|
| 295 | visManager->Draw( *testVolume, redStuff, objectTransformation ); |
---|
| 296 | |
---|
| 297 | if ( errorIndex > 0 ) { |
---|
| 298 | |
---|
| 299 | // |
---|
| 300 | // Recover information from log file |
---|
| 301 | // |
---|
| 302 | G4int error = GetLoggedPV( logger, errorIndex, p, v ); |
---|
| 303 | if (error) return error; |
---|
| 304 | |
---|
| 305 | G4cout << "DrawError: p=" << p << ", v=" << v << G4endl; |
---|
| 306 | |
---|
| 307 | // |
---|
| 308 | // Draw away |
---|
| 309 | // |
---|
| 310 | // visManager->ClearView(); |
---|
| 311 | |
---|
| 312 | // |
---|
| 313 | // This draws the trajectory |
---|
| 314 | // |
---|
| 315 | G4VisAttributes blueStuff( G4Color(0,0,1) ); |
---|
| 316 | |
---|
| 317 | G4Polyline polyline; |
---|
| 318 | polyline.push_back( p ); |
---|
| 319 | polyline.push_back( p + 4*v*m ); |
---|
| 320 | polyline.SetVisAttributes( blueStuff ); |
---|
| 321 | visManager->Draw( polyline, objectTransformation); |
---|
| 322 | |
---|
| 323 | // |
---|
| 324 | // This draws the initial point p |
---|
| 325 | // |
---|
| 326 | G4Circle circle(p); |
---|
| 327 | circle.SetWorldSize( 5*cm ); |
---|
| 328 | circle.SetVisAttributes( blueStuff ); |
---|
| 329 | visManager->Draw( circle, objectTransformation); |
---|
| 330 | } |
---|
| 331 | |
---|
| 332 | // visManager->Show(); |
---|
| 333 | |
---|
| 334 | return 0; |
---|
| 335 | } |
---|
| 336 | |
---|
| 337 | |
---|
| 338 | // |
---|
| 339 | // DebugInside |
---|
| 340 | // |
---|
| 341 | // Recover previously logged error and invoke G4VSolid::Inside |
---|
| 342 | // |
---|
| 343 | G4int SBTrun::DebugInside( const G4VSolid *testVolume, std::istream &logger, const G4int errorIndex ) const |
---|
| 344 | { |
---|
| 345 | G4ThreeVector p, v; |
---|
| 346 | |
---|
| 347 | // |
---|
| 348 | // Recover information from log file |
---|
| 349 | // |
---|
| 350 | G4int error = GetLoggedPV( logger, errorIndex, p, v ); |
---|
| 351 | if (error) return error; |
---|
| 352 | |
---|
| 353 | // |
---|
| 354 | // Call |
---|
| 355 | // |
---|
| 356 | G4cout << "testVolume->Inside(p): " << testVolume->Inside( p ) << G4endl; |
---|
| 357 | return 0; |
---|
| 358 | } |
---|
| 359 | |
---|
| 360 | |
---|
| 361 | // |
---|
| 362 | // DebugToInP |
---|
| 363 | // |
---|
| 364 | // Recover previously logged error and invoke G4VSolid::DistanceToIn(p) |
---|
| 365 | // |
---|
| 366 | G4int SBTrun::DebugToInP( const G4VSolid *testVolume, std::istream &logger, const G4int errorIndex ) const |
---|
| 367 | { |
---|
| 368 | G4ThreeVector p, v; |
---|
| 369 | |
---|
| 370 | // |
---|
| 371 | // Recover information from log file |
---|
| 372 | // |
---|
| 373 | G4int error = GetLoggedPV( logger, errorIndex, p, v ); |
---|
| 374 | if (error) return error; |
---|
| 375 | |
---|
| 376 | // |
---|
| 377 | // Call |
---|
| 378 | // |
---|
| 379 | G4cout << "testVolume->DistanceToIn(p): " << testVolume->DistanceToIn( p ) << G4endl; |
---|
| 380 | return 0; |
---|
| 381 | } |
---|
| 382 | |
---|
| 383 | |
---|
| 384 | // |
---|
| 385 | // DebugToInPV |
---|
| 386 | // |
---|
| 387 | // Recover previously logged error and invoke G4VSolid::DistanceToIn(p,v) |
---|
| 388 | // |
---|
| 389 | G4int SBTrun::DebugToInPV( const G4VSolid *testVolume, std::istream &logger, const G4int errorIndex ) const |
---|
| 390 | { |
---|
| 391 | G4ThreeVector p, v; |
---|
| 392 | |
---|
| 393 | // |
---|
| 394 | // Recover information from log file |
---|
| 395 | // |
---|
| 396 | G4int error = GetLoggedPV( logger, errorIndex, p, v ); |
---|
| 397 | if (error) return error; |
---|
| 398 | |
---|
| 399 | // |
---|
| 400 | // Call |
---|
| 401 | // |
---|
| 402 | G4double answer = testVolume->DistanceToIn( p, v ); |
---|
| 403 | G4cout << "testVolume->DistanceToIn(p,v): " << answer << G4endl; |
---|
| 404 | |
---|
| 405 | p += answer*v; |
---|
| 406 | |
---|
| 407 | G4cout << "testVolume->Inside(p+=answer*v):" << p << " " << testVolume->Inside(p) << G4endl; |
---|
| 408 | return 0; |
---|
| 409 | } |
---|
| 410 | |
---|
| 411 | |
---|
| 412 | // |
---|
| 413 | // DebugToOutP |
---|
| 414 | // |
---|
| 415 | // Recover previously logged error and invoke G4VSolid::DistanceToOut(p) |
---|
| 416 | // |
---|
| 417 | G4int SBTrun::DebugToOutP( const G4VSolid *testVolume, std::istream &logger, const G4int errorIndex ) const |
---|
| 418 | { |
---|
| 419 | G4ThreeVector p, v; |
---|
| 420 | |
---|
| 421 | // |
---|
| 422 | // Recover information from log file |
---|
| 423 | // |
---|
| 424 | G4int error = GetLoggedPV( logger, errorIndex, p, v ); |
---|
| 425 | if (error) return error; |
---|
| 426 | |
---|
| 427 | // |
---|
| 428 | // Call |
---|
| 429 | // |
---|
| 430 | G4cout << "testVolume->DistanceToOut(p): " << testVolume->DistanceToOut( p ) << G4endl; |
---|
| 431 | return 0; |
---|
| 432 | } |
---|
| 433 | |
---|
| 434 | |
---|
| 435 | // |
---|
| 436 | // DebugToOutPV |
---|
| 437 | // |
---|
| 438 | // Recover previously logged error and invoke G4VSolid::DistanceToOut(p,v) |
---|
| 439 | // |
---|
| 440 | G4int SBTrun::DebugToOutPV( const G4VSolid *testVolume, std::istream &logger, const G4int errorIndex ) const |
---|
| 441 | { |
---|
| 442 | G4ThreeVector p, v; |
---|
| 443 | |
---|
| 444 | // |
---|
| 445 | // Recover information from log file |
---|
| 446 | // |
---|
| 447 | G4int error = GetLoggedPV( logger, errorIndex, p, v ); |
---|
| 448 | if (error) return error; |
---|
| 449 | |
---|
| 450 | // |
---|
| 451 | // Call |
---|
| 452 | // |
---|
| 453 | G4bool validNorm; |
---|
| 454 | G4ThreeVector norm; |
---|
| 455 | G4double answer = testVolume->DistanceToOut( p, v ); |
---|
| 456 | G4cout << "testVolume->DistanceToOut( p, v ): " << answer << " validNorm: " << validNorm << G4endl; |
---|
| 457 | |
---|
| 458 | p += answer*v; |
---|
| 459 | G4cout << "testVolume->Inside(p += answer*v): " << testVolume->Inside(p) << G4endl; |
---|
| 460 | |
---|
| 461 | return 0; |
---|
| 462 | } |
---|
| 463 | |
---|
| 464 | // |
---|
| 465 | // DebugSurfNorm |
---|
| 466 | // |
---|
| 467 | // Recover previously logged error and invoke G4VSolid::SurfaceNormal(p) |
---|
| 468 | // |
---|
| 469 | G4int SBTrun::DebugSurfNorm( const G4VSolid *testVolume, std::istream &logger, const G4int errorIndex ) const |
---|
| 470 | { |
---|
| 471 | G4ThreeVector p, v; |
---|
| 472 | |
---|
| 473 | // |
---|
| 474 | // Recover information from log file |
---|
| 475 | // |
---|
| 476 | G4int error = GetLoggedPV( logger, errorIndex, p, v ); |
---|
| 477 | if (error) return error; |
---|
| 478 | |
---|
| 479 | // |
---|
| 480 | // Call |
---|
| 481 | // |
---|
| 482 | G4ThreeVector norm = testVolume->SurfaceNormal(p); |
---|
| 483 | G4cout << "testVolume->SurfaceNormal(p): " << norm << G4endl; |
---|
| 484 | |
---|
| 485 | G4double answer = norm.dot(v); |
---|
| 486 | G4cout << "norm.dot(v): " << answer << G4endl; |
---|
| 487 | |
---|
| 488 | return 0; |
---|
| 489 | } |
---|
| 490 | |
---|
| 491 | // |
---|
| 492 | // -------------------------------------- |
---|
| 493 | // Tests |
---|
| 494 | // |
---|
| 495 | |
---|
| 496 | // |
---|
| 497 | // TestOutsidePoint |
---|
| 498 | // |
---|
| 499 | void SBTrun::TestOutsidePoint( const G4VSolid *testVolume, G4int *nError, |
---|
| 500 | const SBTrunPointList *inside, const SBTrunPointList *outside, |
---|
| 501 | const G4ThreeVector point, std::ostream &logger ) |
---|
| 502 | { |
---|
| 503 | G4int i, n = inside->NumPoints(); |
---|
| 504 | |
---|
| 505 | G4double safeDistance = testVolume->DistanceToIn( point ); |
---|
| 506 | |
---|
| 507 | if (safeDistance <= 0.0) { |
---|
| 508 | ReportError( nError, point, 0, safeDistance,"T0: DistanceToIn(p) <= 0", logger ); |
---|
| 509 | return; |
---|
| 510 | } |
---|
| 511 | |
---|
| 512 | for( i=0; i < n; i++ ) { |
---|
| 513 | G4ThreeVector vr = (*inside)[i] - point; |
---|
| 514 | G4ThreeVector v = vr.unit(); |
---|
| 515 | |
---|
| 516 | G4double dist = testVolume->DistanceToIn( point, v ); |
---|
| 517 | if (dist <= 0) { |
---|
| 518 | ReportError( nError, point, v, safeDistance, "T0: DistanceToIn(p,v) <= 0", logger ); |
---|
| 519 | continue; |
---|
| 520 | } |
---|
| 521 | if (dist == kInfinity) { |
---|
| 522 | ReportError( nError, point, v, safeDistance, "T0: DistanceToIn(p,v) == kInfinity", logger ); |
---|
| 523 | continue; |
---|
| 524 | } |
---|
| 525 | if (dist < safeDistance-1E-10) { |
---|
| 526 | ReportError( nError, point, v, safeDistance, "T0: DistanceToIn(p,v) < DistanceToIn(p)", logger ); |
---|
| 527 | continue; |
---|
| 528 | } |
---|
| 529 | |
---|
| 530 | G4ThreeVector p = point + dist*v; |
---|
| 531 | |
---|
| 532 | EInside insideOrNot = testVolume->Inside( p ); |
---|
| 533 | if (insideOrNot == kOutside) { |
---|
| 534 | ReportError( nError, point, v, safeDistance, "T0: DistanceToIn(p,v) undershoots", logger ); |
---|
| 535 | continue; |
---|
| 536 | } |
---|
| 537 | if (insideOrNot == kInside) { |
---|
| 538 | ReportError( nError, point, v, safeDistance, "TO: DistanceToIn(p,v) overshoots", logger ); |
---|
| 539 | continue; |
---|
| 540 | } |
---|
| 541 | |
---|
| 542 | dist = testVolume->DistanceToIn( p ); |
---|
| 543 | |
---|
| 544 | //if (dist != 0) { |
---|
| 545 | if (dist > G4GeometryTolerance::GetInstance()->GetSurfaceTolerance()) { |
---|
| 546 | ReportError( nError, p, v, safeDistance, "T02: DistanceToIn(p) should be zero", logger ); |
---|
| 547 | // logger << "Dist != 0 : " << dist << endl; |
---|
| 548 | continue; |
---|
| 549 | } |
---|
| 550 | |
---|
| 551 | dist = testVolume->DistanceToOut( p ); |
---|
| 552 | //if (dist != 0) { |
---|
| 553 | if (dist > G4GeometryTolerance::GetInstance()->GetSurfaceTolerance()) { |
---|
| 554 | ReportError( nError, p, v, safeDistance, "T02: DistanceToOut(p) should be zero", logger ); |
---|
| 555 | continue; |
---|
| 556 | } |
---|
| 557 | |
---|
| 558 | dist = testVolume->DistanceToIn( p, v ); |
---|
| 559 | safeDistance = testVolume->DistanceToIn( p ); |
---|
| 560 | // |
---|
| 561 | // Beware! We might expect dist to be precisely zero, but this may not |
---|
| 562 | // be true at corners due to roundoff of the calculation of p = point + dist*v. |
---|
| 563 | // It should, however, *not* be infinity. |
---|
| 564 | // |
---|
| 565 | if (dist == kInfinity) { |
---|
| 566 | ReportError( nError, p, v, safeDistance, "T02: DistanceToIn(p,v) == kInfinity", logger ); |
---|
| 567 | continue; |
---|
| 568 | } |
---|
| 569 | |
---|
| 570 | G4bool validNorm; |
---|
| 571 | G4ThreeVector norm; |
---|
| 572 | |
---|
| 573 | dist = testVolume->DistanceToOut( p, v, true, &validNorm, &norm ); |
---|
| 574 | if (dist == 0) continue; |
---|
| 575 | |
---|
| 576 | if (dist == kInfinity) { |
---|
| 577 | ReportError( nError, p, v, safeDistance, "T02: DistanceToOut(p,v) == kInfinity", logger ); |
---|
| 578 | continue; |
---|
| 579 | } |
---|
| 580 | else if (dist < 0) { |
---|
| 581 | ReportError( nError, p, v, safeDistance, "T02: DistanceToOut(p,v) < 0", logger ); |
---|
| 582 | continue; |
---|
| 583 | } |
---|
| 584 | |
---|
| 585 | if (validNorm) { |
---|
| 586 | if (norm.dot(v) < 0) { |
---|
| 587 | ReportError( nError, p, v, safeDistance, "T02: Outgoing normal incorrect", logger ); |
---|
| 588 | continue; |
---|
| 589 | } |
---|
| 590 | } |
---|
| 591 | |
---|
| 592 | G4ThreeVector norm1 = testVolume->SurfaceNormal( p ); |
---|
| 593 | if (norm1.dot(v) > 0) { |
---|
| 594 | ReportError( nError, p, v, safeDistance, "T02: Ingoing surfaceNormal is incorrect", logger ); |
---|
| 595 | } |
---|
| 596 | |
---|
| 597 | |
---|
| 598 | G4ThreeVector p2 = p + v*dist; |
---|
| 599 | |
---|
| 600 | insideOrNot = testVolume->Inside(p2); |
---|
| 601 | if (insideOrNot == kInside) { |
---|
| 602 | ReportError( nError, p, v, safeDistance, "T02: DistanceToOut(p,v) undershoots", logger ); |
---|
| 603 | continue; |
---|
| 604 | } |
---|
| 605 | if (insideOrNot == kOutside) { |
---|
| 606 | ReportError( nError, p, v, safeDistance, "TO2: DistanceToOut(p,v) overshoots", logger ); |
---|
| 607 | continue; |
---|
| 608 | } |
---|
| 609 | |
---|
| 610 | G4ThreeVector norm2 = testVolume->SurfaceNormal( p2 ); |
---|
| 611 | if (norm2.dot(v) < 0) { |
---|
| 612 | if (testVolume->DistanceToIn(p2,v) != 0) |
---|
| 613 | ReportError( nError, p2, v, safeDistance, "T02: Outgoing surfaceNormal is incorrect", logger ); |
---|
| 614 | } |
---|
| 615 | if (validNorm) { |
---|
| 616 | if (norm.dot(norm2) < 0.0) { |
---|
| 617 | ReportError( nError, p2, v, safeDistance, "T02: SurfaceNormal and DistanceToOut disagree on normal", logger ); |
---|
| 618 | } |
---|
| 619 | } |
---|
| 620 | |
---|
| 621 | if (validNorm) { |
---|
| 622 | dist = testVolume->DistanceToIn(p2,v); |
---|
| 623 | if (dist == 0) { |
---|
| 624 | // |
---|
| 625 | // We may have grazed a corner, which is a problem of design. |
---|
| 626 | // Check distance out |
---|
| 627 | // |
---|
| 628 | if (testVolume->DistanceToOut(p2,v) != 0) { |
---|
| 629 | ReportError( nError, p, v, safeDistance, "TO2: DistanceToOut incorrectly returns validNorm==true (line of sight)(c)", logger ); |
---|
| 630 | continue; |
---|
| 631 | } |
---|
| 632 | } |
---|
| 633 | else if (dist != kInfinity) { |
---|
| 634 | ReportError( nError, p, v, safeDistance, "TO2: DistanceToOut incorrectly returns validNorm==true (line of sight)", logger ); |
---|
| 635 | continue; |
---|
| 636 | } |
---|
| 637 | |
---|
| 638 | G4int j; |
---|
| 639 | for( j=0; j < n; j++ ) { |
---|
| 640 | G4ThreeVector p2top = (*inside)[j] - p2; |
---|
| 641 | |
---|
| 642 | if (p2top.dot(norm) > 0) { |
---|
| 643 | ReportError( nError, p, v,safeDistance, |
---|
| 644 | "T02: DistanceToOut incorrectly returns validNorm==true (horizon)", logger ); |
---|
| 645 | continue; |
---|
| 646 | } |
---|
| 647 | } |
---|
| 648 | } // if valid normal |
---|
| 649 | } // Loop over inside points |
---|
| 650 | |
---|
| 651 | n = outside->NumPoints(); |
---|
| 652 | |
---|
| 653 | for( i=0; i < n; i++ ) { |
---|
| 654 | G4ThreeVector vr = (*outside)[i] - point; |
---|
| 655 | if (vr.mag2() < DBL_MIN) continue; |
---|
| 656 | |
---|
| 657 | G4ThreeVector v = vr.unit(); |
---|
| 658 | |
---|
| 659 | #if DEBUG |
---|
| 660 | G4cerr << std::setprecision(32); |
---|
| 661 | G4cerr << "point = " << point << G4endl ; |
---|
| 662 | G4cerr << "v = " << v << G4endl ; |
---|
| 663 | #endif |
---|
| 664 | |
---|
| 665 | G4double dist = testVolume->DistanceToIn( point, v ); |
---|
| 666 | |
---|
| 667 | #if DEBUG |
---|
| 668 | G4cerr << "finished ..\n"; |
---|
| 669 | #endif |
---|
| 670 | |
---|
| 671 | if (dist <= 0) { |
---|
| 672 | ReportError( nError, point, v, safeDistance, "T03: DistanceToIn(p,v) <= 0", logger ); |
---|
| 673 | continue; |
---|
| 674 | } |
---|
| 675 | if (dist == kInfinity) { |
---|
| 676 | //G4cout << "dist == kInfinity" << G4endl ; |
---|
| 677 | continue; |
---|
| 678 | } |
---|
| 679 | if (dist < safeDistance-1E-10) { |
---|
| 680 | ReportError( nError, point, v, safeDistance, "T03: DistanceToIn(p,v) < DistanceToIn(p)", logger ); |
---|
| 681 | continue; |
---|
| 682 | } |
---|
| 683 | G4ThreeVector p = point + dist*v; |
---|
| 684 | |
---|
| 685 | EInside insideOrNot = testVolume->Inside( p ); |
---|
| 686 | if (insideOrNot == kOutside) { |
---|
| 687 | ReportError( nError, point, v, safeDistance, "T03: DistanceToIn(p,v) undershoots", logger ); |
---|
| 688 | continue; |
---|
| 689 | } |
---|
| 690 | if (insideOrNot == kInside) { |
---|
| 691 | ReportError( nError, point, v, safeDistance, "TO3: DistanceToIn(p,v) overshoots", logger ); |
---|
| 692 | continue; |
---|
| 693 | } |
---|
| 694 | } // Loop over outside points |
---|
| 695 | } |
---|
| 696 | |
---|
| 697 | // |
---|
| 698 | // TestInsidePoint |
---|
| 699 | // |
---|
| 700 | void SBTrun::TestInsidePoint( const G4VSolid *testVolume, G4int *nError, |
---|
| 701 | const SBTrunPointList *outside, const G4ThreeVector point, std::ostream &logger ) |
---|
| 702 | { |
---|
| 703 | G4int i, n = outside->NumPoints(); |
---|
| 704 | |
---|
| 705 | G4double safeDistance = testVolume->DistanceToOut( point ); |
---|
| 706 | if (safeDistance <= 0.0) { |
---|
| 707 | ReportError( nError, point, 0, safeDistance, "TI: DistanceToOut(p) <= 0", logger ); |
---|
| 708 | return; |
---|
| 709 | } |
---|
| 710 | |
---|
| 711 | G4VisExtent visExtent(testVolume->GetExtent()); |
---|
| 712 | if (point.x() < visExtent.GetXmin() || |
---|
| 713 | point.x() > visExtent.GetXmax() || |
---|
| 714 | point.y() < visExtent.GetYmin() || |
---|
| 715 | point.y() > visExtent.GetYmax() || |
---|
| 716 | point.z() < visExtent.GetZmin() || |
---|
| 717 | point.z() > visExtent.GetZmax() ) { |
---|
| 718 | ReportError( nError, point, 0, safeDistance, "TI: Point is outside G4VisExtent", logger ); |
---|
| 719 | } |
---|
| 720 | |
---|
| 721 | for( i=0; i < n; i++ ) { |
---|
| 722 | G4ThreeVector vr = (*outside)[i] - point; |
---|
| 723 | G4ThreeVector v = vr.unit(); |
---|
| 724 | |
---|
| 725 | G4bool validNorm; |
---|
| 726 | G4ThreeVector norm; |
---|
| 727 | |
---|
| 728 | G4double dist = testVolume->DistanceToOut( point, v, true, &validNorm, &norm ); |
---|
| 729 | G4double NormalDist ; |
---|
| 730 | |
---|
| 731 | NormalDist = testVolume->DistanceToOut( point ); |
---|
| 732 | |
---|
| 733 | if (dist <= 0) { |
---|
| 734 | ReportError( nError, point, v, safeDistance, "TI: DistanceToOut(p,v) <= 0 Normal Dist = ", logger ); |
---|
| 735 | continue; |
---|
| 736 | } |
---|
| 737 | if (dist == kInfinity) { |
---|
| 738 | ReportError( nError, point, v, safeDistance, "TI: DistanceToOut(p,v) == kInfinity", logger ); |
---|
| 739 | continue; |
---|
| 740 | } |
---|
| 741 | if (dist < safeDistance-1E-10) { |
---|
| 742 | ReportError( nError, point, v, safeDistance, "TI: DistanceToOut(p,v) < DistanceToIn(p)", logger ); |
---|
| 743 | continue; |
---|
| 744 | } |
---|
| 745 | |
---|
| 746 | if (validNorm) { |
---|
| 747 | if (norm.dot(v) < 0) { |
---|
| 748 | ReportError( nError, point, v, safeDistance, "TI: Outgoing normal incorrect", logger ); |
---|
| 749 | continue; |
---|
| 750 | } |
---|
| 751 | } |
---|
| 752 | |
---|
| 753 | G4ThreeVector p = point + v*dist; |
---|
| 754 | |
---|
| 755 | EInside insideOrNot = testVolume->Inside(p); |
---|
| 756 | if (insideOrNot == kInside) { |
---|
| 757 | ReportError( nError, point, v, safeDistance, "TI: DistanceToOut(p,v) undershoots", logger ); |
---|
| 758 | continue; |
---|
| 759 | } |
---|
| 760 | if (insideOrNot == kOutside) { |
---|
| 761 | ReportError( nError, point, v, safeDistance, "TI: DistanceToOut(p,v) overshoots", logger ); |
---|
| 762 | continue; |
---|
| 763 | } |
---|
| 764 | |
---|
| 765 | G4ThreeVector norm1 = testVolume->SurfaceNormal( p ); |
---|
| 766 | if (norm1.dot(v) < 0) { |
---|
| 767 | if (testVolume->DistanceToIn(p,v) != 0) |
---|
| 768 | ReportError( nError, p, v, safeDistance, "TI: SurfaceNormal is incorrect", logger ); |
---|
| 769 | } |
---|
| 770 | } |
---|
| 771 | } |
---|
| 772 | |
---|
| 773 | |
---|
| 774 | // |
---|
| 775 | // ReportError |
---|
| 776 | // |
---|
| 777 | // Report the specified error message, but only if it has not been reported a zillion |
---|
| 778 | // times already. |
---|
| 779 | // |
---|
| 780 | void SBTrun::ReportError( G4int *nError, const G4ThreeVector p, |
---|
| 781 | const G4ThreeVector v, G4double distance, |
---|
| 782 | const G4String comment, std::ostream &logger ) |
---|
| 783 | { |
---|
| 784 | // |
---|
| 785 | // Have we encountered this error message before? |
---|
| 786 | // |
---|
| 787 | SBTrunErrorList *last=0, *errors = errorList; |
---|
| 788 | while( errors ) { |
---|
| 789 | // |
---|
| 790 | // Note: below is an expensive comparison. This could be replaced with something |
---|
| 791 | // faster if we really wanted, like a hash. But is it worth the trouble? Nah. |
---|
| 792 | // |
---|
| 793 | if (errors->message == comment) { |
---|
| 794 | // |
---|
| 795 | // Yup. Increment count, and return now if the count is too high |
---|
| 796 | // |
---|
| 797 | if ( ++errors->nUsed > 5 ) return; |
---|
| 798 | break; |
---|
| 799 | } |
---|
| 800 | last = errors; |
---|
| 801 | errors = errors->next; |
---|
| 802 | } |
---|
| 803 | |
---|
| 804 | if (errors == 0) { |
---|
| 805 | // |
---|
| 806 | // New error: add it the end of our list |
---|
| 807 | // |
---|
| 808 | errors = new SBTrunErrorList; |
---|
| 809 | errors->message = comment; |
---|
| 810 | errors->nUsed = 1; |
---|
| 811 | errors->next = 0; |
---|
| 812 | if (errorList) |
---|
| 813 | last->next = errors; |
---|
| 814 | else |
---|
| 815 | errorList = errors; |
---|
| 816 | } |
---|
| 817 | |
---|
| 818 | // |
---|
| 819 | // Output the message |
---|
| 820 | // |
---|
| 821 | logger << "% " << comment; |
---|
| 822 | if (errors->nUsed == 5) logger << " (any further such errors suppressed)"; |
---|
| 823 | logger << " DistanceToIn = " << distance ; |
---|
| 824 | logger << G4endl; |
---|
| 825 | |
---|
| 826 | logger << ++(*nError) << " " << p.x() << " " << p.y() << " " << p.z() |
---|
| 827 | << " " << v.x() << " " << v.y() << " " << v.z() << G4endl; |
---|
| 828 | } |
---|
| 829 | |
---|
| 830 | |
---|
| 831 | // |
---|
| 832 | // ClearErrors |
---|
| 833 | // |
---|
| 834 | // Reset list of errors (and clear memory) |
---|
| 835 | // |
---|
| 836 | void SBTrun::ClearErrors() |
---|
| 837 | { |
---|
| 838 | SBTrunErrorList *here, *next; |
---|
| 839 | |
---|
| 840 | here = errorList; |
---|
| 841 | while( here ) { |
---|
| 842 | next = here->next; |
---|
| 843 | delete here; |
---|
| 844 | here = next; |
---|
| 845 | } |
---|
| 846 | errorList = 0; |
---|
| 847 | } |
---|
| 848 | |
---|
| 849 | // |
---|
| 850 | // CountErrors |
---|
| 851 | // |
---|
| 852 | // Count up all errors |
---|
| 853 | // |
---|
| 854 | G4int SBTrun::CountErrors() const |
---|
| 855 | { |
---|
| 856 | SBTrunErrorList *here; |
---|
| 857 | G4int answer = 0; |
---|
| 858 | |
---|
| 859 | here = errorList; |
---|
| 860 | while( here ) { |
---|
| 861 | answer += here->nUsed; |
---|
| 862 | here = here->next; |
---|
| 863 | } |
---|
| 864 | |
---|
| 865 | return answer; |
---|
| 866 | } |
---|
| 867 | |
---|
| 868 | |
---|
| 869 | // |
---|
| 870 | // GetLoggedPV |
---|
| 871 | // |
---|
| 872 | // Get the p and v vectors stored in a test3 log file |
---|
| 873 | // |
---|
| 874 | G4int SBTrun::GetLoggedPV( std::istream &logger, const G4int errorIndex, |
---|
| 875 | G4ThreeVector &p, G4ThreeVector &v ) const |
---|
| 876 | { |
---|
| 877 | // |
---|
| 878 | // Search for the requested error index, skipping comments along the way |
---|
| 879 | // |
---|
| 880 | for(;;) { |
---|
| 881 | while( logger.peek() == '%' ) logger.ignore( 99999, '\n' ); |
---|
| 882 | |
---|
| 883 | if (logger.peek() == EOF) return 1; |
---|
| 884 | |
---|
| 885 | G4int thisIndex; |
---|
| 886 | |
---|
| 887 | logger >> thisIndex; |
---|
| 888 | if (thisIndex == errorIndex) break; |
---|
| 889 | |
---|
| 890 | logger.ignore( 99999, '\n' ); |
---|
| 891 | } |
---|
| 892 | |
---|
| 893 | // |
---|
| 894 | // Extract the vectors |
---|
| 895 | // |
---|
| 896 | // We should probably have some error checking here... |
---|
| 897 | // |
---|
| 898 | G4double x, y, z; |
---|
| 899 | |
---|
| 900 | logger >> x >> y >> z; |
---|
| 901 | p = G4ThreeVector(x*mm,y*mm,z*mm); |
---|
| 902 | |
---|
| 903 | logger >> x >> y >> z; |
---|
| 904 | v = G4ThreeVector(x*mm,y*mm,z*mm); |
---|
| 905 | |
---|
| 906 | return 0; |
---|
| 907 | } |
---|
| 908 | |
---|
| 909 | G4String SBTrun::CurrentSolid = "" ; |
---|
| 910 | |
---|
| 911 | G4String SBTrun::GetCurrentSolid(void) |
---|
| 912 | { |
---|
| 913 | return CurrentSolid; |
---|
| 914 | } |
---|
| 915 | |
---|
| 916 | void SBTrun::SetCurrentSolid(G4String SolidName) |
---|
| 917 | { |
---|
| 918 | CurrentSolid = SolidName; |
---|
| 919 | } |
---|
| 920 | |
---|
| 921 | // |
---|
| 922 | // -------------------------------------- |
---|
| 923 | // SBTrunPointList stuff |
---|
| 924 | // |
---|
| 925 | |
---|
| 926 | SBTrunPointList::SBTrunPointList( G4int size ) |
---|
| 927 | { |
---|
| 928 | pointList = new G4ThreeVector[size]; |
---|
| 929 | maxPoints = size; |
---|
| 930 | numPoints = 0; |
---|
| 931 | } |
---|
| 932 | |
---|
| 933 | SBTrunPointList::~SBTrunPointList() |
---|
| 934 | { |
---|
| 935 | delete [] pointList; |
---|
| 936 | } |
---|
| 937 | |
---|
| 938 | void SBTrunPointList::AddPoint( G4ThreeVector newPoint ) |
---|
| 939 | { |
---|
| 940 | if (numPoints < maxPoints) { |
---|
| 941 | // |
---|
| 942 | // Since we still have space, append the point on the |
---|
| 943 | // end of the list |
---|
| 944 | // |
---|
| 945 | pointList[numPoints] = newPoint; |
---|
| 946 | numPoints++; |
---|
| 947 | } |
---|
| 948 | else { |
---|
| 949 | // |
---|
| 950 | // Our list is filled up, so replace a random |
---|
| 951 | // entry |
---|
| 952 | // |
---|
| 953 | G4int irand = (G4int)((G4UniformRand()*( (G4double)maxPoints ))); |
---|
| 954 | pointList[irand] = newPoint; |
---|
| 955 | } |
---|
| 956 | } |
---|
| 957 | |
---|
| 958 | |
---|
| 959 | |
---|