Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // Class G4VIntersectionLocator implementation 27 // 28 // 27.10.08 - John Apostolakis, Tatiana Nikiti 29 // ------------------------------------------- 30 31 #include <iomanip> 32 #include <sstream> 33 34 #include "globals.hh" 35 #include "G4ios.hh" 36 #include "G4AutoDelete.hh" 37 #include "G4SystemOfUnits.hh" 38 #include "G4VIntersectionLocator.hh" 39 #include "G4TouchableHandle.hh" 40 #include "G4GeometryTolerance.hh" 41 42 ////////////////////////////////////////////// 43 // 44 // Constructor 45 // 46 G4VIntersectionLocator::G4VIntersectionLocator 47 : fiNavigator(theNavigator) 48 { 49 kCarTolerance = G4GeometryTolerance::GetInst 50 51 if( fiNavigator->GetExternalNavigation() == 52 { 53 fHelpingNavigator = new G4Navigator(); 54 } 55 else // Must clone the navigator, together w 56 { 57 fHelpingNavigator = fiNavigator->Clone(); 58 } 59 } 60 61 ////////////////////////////////////////////// 62 // 63 // Destructor. 64 // 65 G4VIntersectionLocator::~G4VIntersectionLocato 66 { 67 delete fHelpingNavigator; 68 delete fpTouchable; 69 } 70 71 ////////////////////////////////////////////// 72 // 73 // Dump status of propagator to cout (old meth 74 // 75 void 76 G4VIntersectionLocator::printStatus( const G4F 77 const G4F 78 G4d 79 G4d 80 G4i 81 { 82 std::ostringstream os; 83 printStatus( StartFT,CurrentFT,requestStep,s 84 G4cout << os.str(); 85 } 86 87 ////////////////////////////////////////////// 88 // 89 // Dumps status of propagator. 90 // 91 void 92 G4VIntersectionLocator::printStatus( const G4F 93 const G4F 94 G4d 95 G4d 96 G4i 97 std 98 G4i 99 { 100 // const G4int verboseLevel= fVerboseLevel; 101 const G4ThreeVector StartPosition = St 102 const G4ThreeVector StartUnitVelocity = St 103 const G4ThreeVector CurrentPosition = Cu 104 const G4ThreeVector CurrentUnitVelocity = Cu 105 106 G4double step_len = CurrentFT.GetCurveLength 107 G4long oldprc; // cout/cerr precision setti 108 109 if( ((stepNo == 0) && (verboseLevel <3)) || 110 { 111 oldprc = os.precision(4); 112 os << std::setw( 6) << " " 113 << std::setw( 25) << " Current Posi 114 << G4endl; 115 os << std::setw( 5) << "Step#" 116 << std::setw(10) << " s " << " " 117 << std::setw(10) << "X(mm)" << " " 118 << std::setw(10) << "Y(mm)" << " " 119 << std::setw(10) << "Z(mm)" << " " 120 << std::setw( 7) << " N_x " << " " 121 << std::setw( 7) << " N_y " << " " 122 << std::setw( 7) << " N_z " << " " 123 os << std::setw( 7) << " Delta|N|" << " " 124 << std::setw( 9) << "StepLen" << " 125 << std::setw(12) << "StartSafety" < 126 << std::setw( 9) << "PhsStep" << " 127 os << G4endl; 128 os.precision(oldprc); 129 } 130 if((stepNo == 0) && (verboseLevel <=3)) 131 { 132 // Recurse to print the start values 133 // 134 printStatus( StartFT, StartFT, -1.0, safet 135 } 136 if( verboseLevel <= 3 ) 137 { 138 if( stepNo >= 0) 139 { 140 os << std::setw( 4) << stepNo << " "; 141 } 142 else 143 { 144 os << std::setw( 5) << "Start" ; 145 } 146 oldprc = os.precision(8); 147 os << std::setw(10) << CurrentFT.GetCurveL 148 os << std::setw(10) << CurrentPosition.x() 149 << std::setw(10) << CurrentPosition 150 << std::setw(10) << CurrentPosition 151 os.precision(4); 152 os << std::setw( 7) << CurrentUnitVelocity 153 << std::setw( 7) << CurrentUnitVelo 154 << std::setw( 7) << CurrentUnitVelo 155 os.precision(3); 156 os << std::setw( 7) 157 << CurrentFT.GetMomentum().mag()- S 158 << " "; 159 os << std::setw( 9) << step_len << " "; 160 os << std::setw(12) << safety << " "; 161 if( requestStep != -1.0 ) 162 { 163 os << std::setw( 9) << requestStep << " 164 } 165 else 166 { 167 os << std::setw( 9) << "Init/NotKnown" < 168 } 169 os << G4endl; 170 os.precision(oldprc); 171 } 172 else // if( verboseLevel > 3 ) 173 { 174 // Multi-line output 175 176 os << "Step taken was " << step_len 177 << " out of PhysicalStep= " << req 178 os << "Final safety is: " << safety << G4e 179 os << "Chord length = " << (CurrentPositio 180 << G4endl; 181 os << G4endl; 182 } 183 } 184 185 ////////////////////////////////////////////// 186 // 187 // ReEstimateEndPoint. 188 // 189 G4FieldTrack G4VIntersectionLocator:: 190 ReEstimateEndpoint( const G4FieldTrack& Curren 191 const G4FieldTrack& Estima 192 G4double , // l 193 G4double 194 #ifdef G4DEBUG_FIELD 195 curveDist 196 #endif 197 ) 198 { 199 G4FieldTrack newEndPoint( CurrentStateA ); 200 auto integrDriver = GetChordFinderFor()->Get 201 202 G4FieldTrack retEndPoint( CurrentStateA ); 203 G4bool goodAdvance; 204 G4int itrial = 0; 205 const G4int no_trials = 20; 206 207 208 G4double endCurveLen= EstimatedEndStateB.Get 209 210 do // Loop checking, 07.10.2016, JA 211 { 212 G4double currentCurveLen = newEndPoint.Get 213 G4double advanceLength = endCurveLen - cur 214 if (std::abs(advanceLength)<kCarTolerance) 215 { 216 goodAdvance=true; 217 } 218 else 219 { 220 goodAdvance = integrDriver->AccurateAdva 221 222 } 223 } 224 while( !goodAdvance && (++itrial < no_trials 225 226 if( goodAdvance ) 227 { 228 retEndPoint = newEndPoint; 229 } 230 else 231 { 232 retEndPoint = EstimatedEndStateB; // Could 233 } 234 235 // All the work is done 236 // below are some diagnostics only -- befor 237 // 238 const G4String MethodName("G4VIntersectionLo 239 240 #ifdef G4VERBOSE 241 G4int latest_good_trials = 0; 242 if( itrial > 1) 243 { 244 if( fVerboseLevel > 0 ) 245 { 246 G4cout << MethodName << " called - goodA 247 << " trials = " << itrial 248 << " previous good= " << latest_g 249 << G4endl; 250 } 251 latest_good_trials = 0; 252 } 253 else 254 { 255 ++latest_good_trials; 256 } 257 #endif 258 259 #ifdef G4DEBUG_FIELD 260 G4double lengthDone = newEndPoint.GetCurveLe 261 - CurrentStateA.GetCurve 262 if( !goodAdvance ) 263 { 264 if( fVerboseLevel >= 3 ) 265 { 266 G4cout << MethodName << "> AccurateAdvan 267 G4cout << " in " << itrial << " integrat 268 G4cout << " It went only " << lengthDone 269 << " -- a difference of " << curv 270 G4cout << " ReEstimateEndpoint> Reset en 271 << G4endl; 272 } 273 } 274 G4double linearDist = ( EstimatedEndStateB.G 275 - CurrentStateA.GetPosit 276 static G4int noInaccuracyWarnings = 0; 277 G4int maxNoWarnings = 10; 278 if ( (noInaccuracyWarnings < maxNoWarnings 279 || (fVerboseLevel > 1) ) 280 { 281 G4ThreeVector move = newEndPoint.GetPosi 282 - EstimatedEndStateB. 283 std::ostringstream message; 284 message.precision(12); 285 message << " Integration inaccuracy requ 286 << " an adjustment in the step's 287 << " Two mid-points are furthe 288 << " curve length difference" 289 << " Dist = " << linearD 290 << " curve length = " << curveDi 291 message << " Correction applied is " << 292 << " Old Estimated B position= 293 << EstimatedEndStateB.GetPositio 294 << " Recalculated Position= 295 << newEndPoint.GetPosition() << 296 << " Change ( new - old ) = 297 G4Exception("G4VIntersectionLocator::ReE 298 "GeomNav1002", JustWarning, 299 } 300 /* 301 #else 302 // Statistics on the RMS value of the correc 303 304 static G4ThreadLocal G4int noCorrections = 0 305 ++noCorrections; 306 if( goodAdvance ) 307 { 308 static G4ThreadLocal G4double sumCorrectio 309 sumCorrectionsSq += (EstimatedEndStateB.Ge 310 newEndPoint.GetPositi 311 } 312 */ 313 #endif 314 315 return retEndPoint; 316 } 317 318 319 ////////////////////////////////////////////// 320 // 321 // ReEstimateEndPoint. 322 // 323 // The following values are returned in curv 324 // 0: Normal - no problem 325 // 1: Unexpected co-incidence - milder m 326 // 2: Real mixup - EndB is NOT past Star 327 // ( ie. StartA's curve-lengh is bi 328 329 330 G4bool G4VIntersectionLocator:: 331 CheckAndReEstimateEndpoint( const G4FieldTrack 332 const G4FieldTrack 333 G4FieldTrack 334 G4int& 335 { 336 G4double linDistSq, curveDist; 337 338 G4bool recalculated = false; 339 curveError= 0; 340 341 linDistSq = ( EstimatedEndB.GetPosition() 342 - CurrentStartA.GetPosition() ). 343 curveDist = EstimatedEndB.GetCurveLength() 344 - CurrentStartA.GetCurveLength(); 345 if( (curveDist>=0.0) 346 && (curveDist*curveDist *(1.0+2.0*fiEpsil 347 { 348 G4FieldTrack newEndPointFT = EstimatedEndB 349 350 if (curveDist>0.0) 351 { 352 // Re-integrate to obtain a new B 353 RevisedEndPoint = ReEstimateEndpoint( C 354 E 355 l 356 c 357 recalculated = true; 358 } 359 else 360 { 361 // Zero length -> no advance! 362 newEndPointFT = CurrentStartA; 363 recalculated = true; 364 curveError = 1; // Unexpected co-incid 365 366 G4Exception("G4MultiLevelLocator::Estim 367 "GeomNav1002", JustWarning, 368 "A & B are at equal distance in 2nd 369 } 370 } 371 372 // Sanity check 373 // 374 if( curveDist < 0.0 ) 375 { 376 curveError = 2; // Real mixup 377 } 378 return recalculated; 379 } 380 381 ////////////////////////////////////////////// 382 // 383 // Method for finding SurfaceNormal of Interse 384 // 385 G4ThreeVector G4VIntersectionLocator:: 386 GetLocalSurfaceNormal(const G4ThreeVector& Cur 387 { 388 G4ThreeVector Normal(G4ThreeVector(0.0,0.0,0 389 G4VPhysicalVolume* located; 390 391 validNormal = false; 392 fHelpingNavigator->SetWorldVolume(GetNavigat 393 located = fHelpingNavigator->LocateGlobalPoi 394 395 delete fpTouchable; 396 fpTouchable = fHelpingNavigator->CreateTouch 397 398 // To check if we can use GetGlobalExitNorma 399 // 400 G4ThreeVector localPosition = fpTouchable->G 401 ->GetTopTransform().TransformP 402 403 // Issue: in the case of coincident surfaces 404 // which side you are located onto (c 405 // TO-DO: use direction (of chord) to identi 406 407 if( located != nullptr) 408 { 409 G4LogicalVolume* pLogical= located->GetLog 410 G4VSolid* pSolid; 411 412 if( (pLogical != nullptr) && ( (pSolid=pLo 413 { 414 if ( ( pSolid->Inside(localPosition)==kS 415 || ( pSolid->DistanceToOut(localPositi 416 { 417 Normal = pSolid->SurfaceNormal(localPo 418 validNormal = true; 419 420 #ifdef G4DEBUG_FIELD 421 if( std::fabs(Normal.mag2() - 1.0 ) > 422 { 423 G4cerr << "PROBLEM in G4VIntersectio 424 << G4endl; 425 G4cerr << " Normal is not unit - ma 426 G4cerr << " at trial local point " 427 G4cerr << " Solid is " << *pSolid 428 } 429 #endif 430 } 431 } 432 } 433 return Normal; 434 } 435 436 ////////////////////////////////////////////// 437 // 438 // Adjustment of Found Intersection 439 // 440 G4bool G4VIntersectionLocator:: 441 AdjustmentOfFoundIntersection( const G4ThreeVe 442 const G4ThreeVe 443 const G4ThreeVe 444 const G4ThreeVe 445 const G4bool 446 G4ThreeVe 447 G4double& 448 G4double& 449 G4ThreeVe 450 { 451 G4double dist,lambda; 452 G4ThreeVector Normal, NewPoint, Point_G; 453 G4bool goodAdjust = false, Intersects_FP = f 454 455 // Get SurfaceNormal of Intersecting Solid 456 // 457 Normal = GetGlobalSurfaceNormal(CurrentE_Poi 458 if(!validNormal) { return false; } 459 460 // Intersection between Line and Plane 461 // 462 G4double n_d_m = Normal.dot(MomentumDir); 463 if ( std::abs(n_d_m)>kCarTolerance ) 464 { 465 #ifdef G4VERBOSE 466 if ( fVerboseLevel>1 ) 467 { 468 G4Exception("G4VIntersectionLocator::Adj 469 "GeomNav0003", JustWarning, 470 "No intersection. Parallels 471 } 472 #endif 473 lambda =- Normal.dot(CurrentF_Point-Curren 474 475 // New candidate for Intersection 476 // 477 NewPoint = CurrentF_Point+lambda*MomentumD 478 479 // Distance from CurrentF to Calculated In 480 // 481 dist = std::abs(lambda); 482 483 if ( dist<kCarTolerance*0.001 ) { return 484 485 // Calculation of new intersection point o 486 // 487 if ( IntersectAF ) // First part interse 488 { 489 G4double stepLengthFP; 490 G4ThreeVector Point_P = CurrentA_Point; 491 GetNavigatorFor()->LocateGlobalPointWith 492 Intersects_FP = IntersectChord( Point_P, 493 fPreviou 494 stepLeng 495 496 } 497 else // Second part intersects 498 { 499 G4double stepLengthFP; 500 GetNavigatorFor()->LocateGlobalPointWith 501 Intersects_FP = IntersectChord( CurrentF 502 fPreviou 503 stepLeng 504 } 505 if ( Intersects_FP ) 506 { 507 goodAdjust = true; 508 IntersectionPoint = Point_G; 509 } 510 } 511 512 return goodAdjust; 513 } 514 515 ////////////////////////////////////////////// 516 // 517 // GetSurfaceNormal. 518 // 519 G4ThreeVector G4VIntersectionLocator:: 520 GetSurfaceNormal(const G4ThreeVector& CurrentI 521 G4bool& validNormal) 522 { 523 G4ThreeVector NormalAtEntry; // ( -10. , -10 524 525 G4ThreeVector NormalAtEntryLast, NormalAtEnt 526 G4bool validNormalLast; 527 528 // Relies on a call to Navigator::ComputeSte 529 // this call 530 // 531 NormalAtEntryLast = GetLastSurfaceNormal( Cu 532 // May return valid=false in cases, includ 533 // - if the candidate volume was not foun 534 // - a replica was involved -- determined 535 // (This list is not complete.) 536 537 #ifdef G4DEBUG_FIELD 538 if ( validNormalLast 539 && ( std::fabs(NormalAtEntryLast.mag2() - 1 540 { 541 std::ostringstream message; 542 message << "PROBLEM: Normal is not unit - 543 << NormalAtEntryLast.mag() << G4en 544 message << " at trial intersection point 545 message << " Obtained from Get *Last* Su 546 G4Exception("G4VIntersectionLocator::GetSu 547 "GeomNav1002", JustWarning, me 548 } 549 #endif 550 551 if( validNormalLast ) 552 { 553 NormalAtEntry = NormalAtEntryLast; 554 } 555 validNormal = validNormalLast; 556 557 return NormalAtEntry; 558 } 559 560 ////////////////////////////////////////////// 561 // 562 // GetGlobalSurfaceNormal. 563 // 564 G4ThreeVector G4VIntersectionLocator:: 565 GetGlobalSurfaceNormal(const G4ThreeVector& Cu 566 G4bool& validNorm 567 { 568 G4ThreeVector localNormal = GetLocalSurfaceN 569 G4AffineTransform localToGlobal = / 570 fHelpingNavigator->GetLocalToGlobalTrans 571 G4ThreeVector globalNormal = localToGlobal.T 572 573 #ifdef G4DEBUG_FIELD 574 if( validNormal && ( std::fabs(globalNormal. 575 { 576 std::ostringstream message; 577 message << "****************************** 578 << G4endl; 579 message << " Bad Normal in G4VIntersection 580 << G4endl; 581 message << " * Constituents: " << G4endl; 582 message << " Local Normal= " << localN 583 message << " Transform: " << G4endl 584 << " Net Translation= " << lo 585 << G4endl 586 << " Net Rotation = " << lo 587 << G4endl; 588 message << " * Result: " << G4endl; 589 message << " Global Normal= " << local 590 message << "****************************** 591 G4Exception("G4VIntersectionLocator::GetGl 592 "GeomNav1002", JustWarning, me 593 } 594 #endif 595 596 return globalNormal; 597 } 598 599 ////////////////////////////////////////////// 600 // 601 // GetLastSurfaceNormal. 602 // 603 G4ThreeVector G4VIntersectionLocator:: 604 GetLastSurfaceNormal( const G4ThreeVector& int 605 G4bool& normalIsVa 606 { 607 G4ThreeVector normalVec; 608 G4bool validNorm; 609 normalVec = fiNavigator->GetGlobalExitNormal 610 normalIsValid = validNorm; 611 612 return normalVec; 613 } 614 615 ////////////////////////////////////////////// 616 // 617 // ReportTrialStep. 618 // 619 void G4VIntersectionLocator::ReportTrialStep( 620 const 621 const 622 const 623 const 624 625 { 626 G4double ABchord_length = ChordAB_v.mag(); 627 G4double MomDir_dot_Norm = NewMomentumDir.do 628 G4double MomDir_dot_ABchord; 629 MomDir_dot_ABchord = (1.0 / ABchord_length) 630 631 std::ostringstream outStream; 632 outStream << std::setw(6) << " Step# " 633 << std::setw(17) << " |ChordEF|(mag)" << " 634 << std::setw(18) << " uMomentum.Normal" << 635 << std::setw(18) << " uMomentum.ABdir " << 636 << std::setw(16) << " AB-dist " << 637 << " Chord Vector (EF) " 638 << G4endl; 639 outStream.precision(7); 640 outStream << " " << std::setw(5) << step_no 641 << " " << std::setw(18) << ChordEF_v.mag() 642 << " " << std::setw(18) << MomDir_dot_Norm 643 << " " << std::setw(18) << MomDir_dot_ABch 644 << " " << std::setw(12) << ABchord_length 645 << " " << ChordEF_v 646 << G4endl; 647 outStream << " MomentumDir= " << " " << NewM 648 << " Normal at Entry E= " << NormalAtEntry 649 << " AB chord = " << ChordAB_v 650 << G4endl; 651 G4cout << outStream.str(); 652 653 if( ( std::fabs(NormalAtEntry.mag2() - 1.0) 654 { 655 std::ostringstream message; 656 message << "Normal is not unit - mag= " << 657 << " ValidNormalAtE = " << 658 G4Exception("G4VIntersectionLocator::Repor 659 "GeomNav1002", JustWarning, me 660 } 661 return; 662 } 663 664 ////////////////////////////////////////////// 665 // 666 // LocateGlobalPointWithinVolumeAndCheck. 667 // 668 // Locate point using navigator: updates state 669 // By default, it assumes that the point is in 670 // and returns true. 671 // In check mode, checks whether the point is 672 // If it is inside, it returns true 673 // If not, issues a warning and returns fals 674 // 675 G4bool G4VIntersectionLocator:: 676 LocateGlobalPointWithinVolumeAndCheck( const G 677 { 678 G4bool good = true; 679 G4Navigator* nav = GetNavigatorFor(); 680 const G4String 681 MethodName("G4VIntersectionLocator::LocateGl 682 683 if( fCheckMode ) 684 { 685 G4bool navCheck= nav->IsCheckModeActive(); 686 nav->CheckMode(true); 687 688 // Identify the current volume 689 690 G4TouchableHandle startTH= nav->CreateTouc 691 G4VPhysicalVolume* motherPhys = startTH->G 692 G4VSolid* motherSolid = startTH-> 693 G4AffineTransform transform = nav->GetGlob 694 G4int motherCopyNo = motherPhys->GetCopyNo 695 696 // Let's check that the point is inside th 697 G4ThreeVector localPosition = transform.T 698 EInside inMother = motherSolid->Ins 699 if( inMother != kInside ) 700 { 701 std::ostringstream message; 702 message << "Position located " 703 << ( inMother == kSurface ? " on 704 << "expected volume" << G4endl 705 << " Safety (from Outside) = " 706 << motherSolid->DistanceToIn(loc 707 G4Exception("G4VIntersectionLocator::Loc 708 "GeomNav1002", JustWarning, 709 } 710 711 // 1. Simple next step - quick relocation 712 // nav->LocateGlobalPointWithinVolume( pos 713 714 // 2. Full relocation - to cross-check ans 715 G4VPhysicalVolume* nextPhysical= nav->Loca 716 if( (nextPhysical != motherPhys) 717 || (nextPhysical->GetCopyNo() != mothe 718 ) 719 { 720 G4Exception("G4VIntersectionLocator::Loc 721 "GeomNav1002", JustWarning, 722 "Position located outside ex 723 } 724 nav->CheckMode(navCheck); // Recover orig 725 } 726 else 727 { 728 nav->LocateGlobalPointWithinVolume( positi 729 } 730 return good; 731 } 732 733 ////////////////////////////////////////////// 734 // 735 // LocateGlobalPointWithinVolumeCheckAndReport 736 // 737 void G4VIntersectionLocator:: 738 LocateGlobalPointWithinVolumeCheckAndReport( c 739 c 740 G 741 { 742 // Save value of Check mode first 743 G4bool oldCheck = GetCheckMode(); 744 745 G4bool ok = LocateGlobalPointWithinVolumeAnd 746 if( !ok ) 747 { 748 std::ostringstream message; 749 message << "Failed point location." << G4e 750 << " Code Location info: " << Co 751 G4Exception("G4VIntersectionLocator::Locat 752 "GeomNav1002", JustWarning, me 753 } 754 755 SetCheckMode( oldCheck ); 756 } 757 758 ////////////////////////////////////////////// 759 // 760 // ReportReversedPoints. 761 // 762 void G4VIntersectionLocator:: 763 ReportReversedPoints( std::ostringstream& msg, 764 const G4FieldTrack& Star 765 const G4FieldTrack& EndP 766 G4double NewSafety 767 const G4FieldTrack& A_Pt 768 const G4FieldTrack& B_Pt 769 const G4FieldTrack& SubS 770 const G4ThreeVector& E_P 771 const G4FieldTrack& Appr 772 G4int substep_no, 773 { 774 // Expect that 'msg' can hold the name of t 775 776 // FieldTrack 'points' A and B have been ta 777 // Whereas A should be before B, it is foun 778 G4int verboseLevel= 5; 779 G4double curveDist = B_PtVel.GetCurveLength 780 G4VIntersectionLocator::printStatus( A_PtVe 781 -1.0, NewSafety, s 782 msg << "Error in advancing propagation." << 783 << " The final curve point is NOT fur 784 << " than the original!" << G4endl 785 << " Going *backwards* from len(A) = 786 << " to len(B) = " << B_PtVel.GetCurve 787 << " Curve distance is " << curveD 788 << G4endl 789 << " Point A' (start) is " << A_Pt 790 << " Point B' (end) is " << B_Pt 791 msg << " fEpsStep= " << epsStep << G4e 792 793 G4long oldprc = msg.precision(20); 794 msg << " In full precision, the position, m 795 << " ... are: " << G4endl; 796 msg << " Point A[0] (Curve start) is " << 797 << " Point S (Sub start) is " << 798 << " Point A' (Current start) is " << 799 << " Point E (Trial Point) is " << 800 << " Point F (Intersection) is " << 801 << " Point B' (Current end) is " << 802 << " Point B[0] (Curve end) is " << 803 << G4endl 804 << " LocateIntersection parameters are 805 << " Substep no (total) = " << su 806 << " Substep no = " << su 807 msg.precision(oldprc); 808 } 809 810 ////////////////////////////////////////////// 811 // 812 // ReportProgress. 813 // 814 void G4VIntersectionLocator::ReportProgress( s 815 const G4Fi 816 const G4Fi 817 G4in 818 const G4Fi 819 const G4Fi 820 G4do 821 G4in 822 823 { 824 oss << "ReportProgress: Current status of in 825 if( depth > 0 ) { oss << " Depth= " << depth 826 oss << " Substep no = " << substep_no << G4e 827 G4int verboseLevel = 5; 828 G4double safetyPrev = -1.0; // Add as argum 829 830 printStatus( StartPointVel, EndPointVel, -1. 831 oss, verboseLevel); 832 oss << " * Start and end-point of requested 833 oss << " ** State of point A: "; 834 printStatus( A_PtVel, A_PtVel, -1.0, safetyP 835 oss, verboseLevel); 836 oss << " ** State of point B: "; 837 printStatus( A_PtVel, B_PtVel, -1.0, safetyL 838 oss, verboseLevel); 839 } 840 841 ////////////////////////////////////////////// 842 // 843 // ReportImmediateHit. 844 // 845 void 846 G4VIntersectionLocator::ReportImmediateHit( co 847 co 848 co 849 850 un 851 { 852 static G4ThreadLocal unsigned int occurredOn 853 static G4ThreadLocal G4ThreeVector* ptrLast 854 if( ptrLast == nullptr ) 855 { 856 ptrLast= new G4ThreeVector( DBL_MAX, DBL_ 857 G4AutoDelete::Register(ptrLast); 858 } 859 G4ThreeVector &lastStart= *ptrLast; 860 861 if( (TrialPoint - StartPosition).mag2() < to 862 { 863 static G4ThreadLocal unsigned int numUnmo 864 static G4ThreadLocal unsigned int numStil 865 866 G4cout << "Intersection F == start A in " 867 G4cout << "Start Point: " << StartPositio 868 G4cout << " Start-Trial: " << TrialPoint 869 G4cout << " Start-last: " << StartPositio 870 871 if( (StartPosition - lastStart).mag() < t 872 { 873 // We are at position of last 'Start' 874 ++numUnmoved; 875 ++numStill; 876 G4cout << " { Unmoved: " << " still#= 877 << " total # = " << numUnmoved 878 } 879 else 880 { 881 numStill = 0; 882 } 883 G4cout << " Occurred: " << ++occurredOnTo 884 G4cout << " out of total calls= " << num 885 G4cout << G4endl; 886 lastStart = StartPosition; 887 } 888 } // End of ReportImmediateHit() 889