Geant4 Cross Reference |
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 // G4ChordFinder implementation 27 // 28 // Author: J.Apostolakis - Design and implementation - 25.02.1997 29 // ------------------------------------------------------------------- 30 31 #include <iomanip> 32 33 #include "G4ChordFinder.hh" 34 #include "G4SystemOfUnits.hh" 35 #include "G4MagneticField.hh" 36 #include "G4Mag_UsualEqRhs.hh" 37 #include "G4MagIntegratorDriver.hh" 38 // #include "G4ClassicalRK4.hh" 39 // #include "G4CashKarpRKF45.hh" 40 // #include "G4NystromRK4.hh" 41 // #include "G4BogackiShampine23.hh" 42 // #include "G4BogackiShampine45.hh" 43 44 #include "G4DormandPrince745.hh" 45 46 // New templated stepper(s) -- avoid virtual calls to equation rhs 47 #include "G4TDormandPrince45.hh" 48 49 // FSAL type driver / steppers ----- 50 #include "G4FSALIntegrationDriver.hh" 51 #include "G4VFSALIntegrationStepper.hh" 52 #include "G4RK547FEq1.hh" 53 // #include "G4RK547FEq2.hh" 54 // #include "G4RK547FEq3.hh" 55 // #include "G4FSALBogackiShampine45.hh" 56 // #include "G4FSALDormandPrince745.hh" 57 58 // Templated type drivers ----- 59 #include "G4IntegrationDriver.hh" 60 #include "G4InterpolationDriver.hh" 61 62 #include "G4HelixHeum.hh" 63 #include "G4BFieldIntegrationDriver.hh" 64 65 #include "G4QSSDriverCreator.hh" 66 67 #include "G4CachedMagneticField.hh" 68 69 #include <cassert> 70 #include <memory> 71 72 G4bool G4ChordFinder::gVerboseCtor = false; 73 // .......................................................................... 74 75 G4ChordFinder::G4ChordFinder(G4VIntegrationDriver* pIntegrationDriver) 76 : fDefaultDeltaChord(0.25 * mm), fIntgrDriver(pIntegrationDriver) 77 { 78 // Simple constructor -- it does not create equation 79 if( gVerboseCtor ) 80 { 81 G4cout << "G4ChordFinder: Simple constructor -- it uses pre-existing driver." << G4endl; 82 } 83 84 fDeltaChord = fDefaultDeltaChord; // Parameters 85 } 86 87 // .......................................................................... 88 89 G4ChordFinder::G4ChordFinder( G4MagneticField* theMagField, 90 G4double stepMinimum, 91 G4MagIntegratorStepper* pItsStepper, 92 G4int stepperDriverId ) 93 : fDefaultDeltaChord(0.25 * mm) 94 { 95 // Construct the Chord Finder 96 // by creating in inverse order the Driver, the Stepper and EqRhs ... 97 constexpr G4int nVar6 = 6; // Components integrated in Usual Equation of motion 98 99 fDeltaChord = fDefaultDeltaChord; // Parameters 100 101 G4cout << " G4ChordFinder: stepperDriverId: " << stepperDriverId << G4endl; 102 103 G4bool useFSALstepper = (stepperDriverId == kFSALStepperType); // Was 1 104 G4bool useTemplatedStepper= (stepperDriverId == kTemplatedStepperType); // Was 2 105 G4bool useRegularStepper = (stepperDriverId == kRegularStepperType); // Was 3 106 G4bool useBfieldDriver = (stepperDriverId == kBfieldDriverType); // Was 4 107 G4bool useG4QSSDriver = (stepperDriverId == kQss2DriverType) || (stepperDriverId == kQss3DriverType); 108 109 if( stepperDriverId == kQss3DriverType) 110 { 111 stepperDriverId = kQss2DriverType; 112 G4cout << " G4ChordFinder: QSS 3 is currently replaced by QSS 2 driver." << G4endl; 113 } 114 115 using EquationType = G4Mag_UsualEqRhs; 116 117 using TemplatedStepperType = 118 G4TDormandPrince45<EquationType,nVar6>; // 5th order embedded method. High efficiency. 119 const char* TemplatedStepperName = 120 "G4TDormandPrince745 (templated Dormand-Prince45, aka DoPri5): 5th/4th Order 7-stage embedded"; 121 122 using RegularStepperType = 123 G4DormandPrince745; // 5th order embedded method. High efficiency. 124 // G4ClassicalRK4; // The old default 125 // G4CashKarpRKF45; // First embedded method in G4 126 // G4BogackiShampine45; // High efficiency 5th order embedded method 127 // G4NystromRK4; // Nystrom stepper 4th order 128 // G4RK547FEq1; // or 2 or 3 129 const char* RegularStepperName = 130 "G4DormandPrince745 (aka DOPRI5): 5th/4th Order 7-stage embedded"; 131 // "BogackiShampine 45 (Embedded 5th/4th Order, 7-stage)"; 132 // "Nystrom stepper 4th order"; 133 134 using NewFsalStepperType = G4DormandPrince745; // Now works -- 2020.10.08 135 // Was G4RK547FEq1; // or 2 or 3 136 const char* NewFSALStepperName = 137 "G4RK574FEq1> FSAL 4th/5th order 7-stage 'Equilibrium-type' #1."; 138 139 #ifdef G4DEBUG_FIELD 140 static G4bool verboseDebug = true; 141 if( verboseDebug ) 142 { 143 G4cout << "G4ChordFinder 2nd Constructor called. " << G4endl; 144 G4cout << " Arguments: " << G4endl 145 << " - min step = " << stepMinimum << G4endl 146 << " - stepper ptr provided : " 147 << ( pItsStepper==nullptr ? " no " : " yes " ) << G4endl; 148 if( pItsStepper==nullptr ) 149 G4cout << " - stepper/driver Id = " << stepperDriverId << " i.e. " 150 << " useFSAL = " << useFSALstepper 151 << " , useTemplated = " << useTemplatedStepper 152 << " , useRegular = " << useRegularStepper 153 << " , useFSAL = " << useFSALstepper 154 << G4endl; 155 } 156 #endif 157 158 // useHigherStepper = forceHigherEffiencyStepper || useHigherStepper; 159 160 auto pEquation = new G4Mag_UsualEqRhs(theMagField); 161 fEquation = pEquation; 162 163 // G4MagIntegratorStepper* regularStepper = nullptr; 164 // G4VFSALIntegrationStepper* fsalStepper = nullptr; // for FSAL steppers only 165 // G4MagIntegratorStepper* oldFSALStepper = nullptr; 166 167 G4bool errorInStepperCreation = false; 168 169 std::ostringstream message; // In case of failure, load with description ! 170 171 if( pItsStepper != nullptr ) 172 { 173 if( gVerboseCtor ) 174 { 175 G4cout << " G4ChordFinder: Creating G4IntegrationDriver<G4MagIntegratorStepper> with " 176 << " stepMinimum = " << stepMinimum 177 << " numVar= " << pItsStepper->GetNumberOfVariables() << G4endl; 178 } 179 180 // Stepper type is not known - so must use base class G4MagIntegratorStepper 181 if(pItsStepper->isQSS()) 182 { 183 // fIntgrDriver = pItsStepper->build_driver(stepMinimum, true); 184 G4Exception("G4ChordFinder::G4ChordFinder()", 185 "GeomField1001", FatalException, 186 "Cannot provide QSS stepper in constructor. User c-tor with Driver instead."); 187 } 188 else 189 { 190 fIntgrDriver = new G4IntegrationDriver<G4MagIntegratorStepper>( stepMinimum, 191 pItsStepper, pItsStepper->GetNumberOfVariables() ); 192 // Stepper type is not known - so must use base class G4MagIntegratorStepper 193 // Non-interpolating driver used by default. 194 // WAS: fIntgrDriver = pItsStepper->build_driver(stepMinimum); // QSS introduction -- axed 195 } 196 // -- Older: 197 // G4cout << " G4ChordFinder: Creating G4MagInt_Driver with " ... 198 // Type is not known - so must use old class 199 // fIntgrDriver = new G4MagInt_Driver( stepMinimum, pItsStepper, 200 // pItsStepper->GetNumberOfVariables()); 201 } 202 else if ( useTemplatedStepper ) 203 { 204 if( gVerboseCtor ) 205 { 206 G4cout << " G4ChordFinder: Creating Templated Stepper of type> " 207 << TemplatedStepperName << G4endl; 208 } 209 // RegularStepperType* regularStepper = nullptr; // To check the exception 210 auto templatedStepper = new TemplatedStepperType(pEquation); 211 // *** ****************** 212 // 213 // Alternative - for G4NystromRK4: 214 // = new G4NystromRK4(pEquation, 0.1*mm ); 215 fRegularStepperOwned = templatedStepper; 216 if( templatedStepper == nullptr ) 217 { 218 message << "Templated Stepper instantiation FAILED." << G4endl; 219 message << "G4ChordFinder: Attempted to instantiate " 220 << TemplatedStepperName << " type stepper " << G4endl; 221 errorInStepperCreation = true; 222 } 223 else 224 { 225 fIntgrDriver = new G4IntegrationDriver<TemplatedStepperType>( 226 stepMinimum, templatedStepper, nVar6 ); 227 if( gVerboseCtor ) 228 { 229 G4cout << " G4ChordFinder: Using G4IntegrationDriver. " << G4endl; 230 } 231 } 232 233 } 234 else if ( useRegularStepper ) // Plain stepper -- not double ... 235 { 236 auto regularStepper = new RegularStepperType(pEquation); 237 // *** ****************** 238 fRegularStepperOwned = regularStepper; 239 240 if( gVerboseCtor ) 241 { 242 G4cout << " G4ChordFinder: Creating Driver for regular stepper."; 243 } 244 245 if( regularStepper == nullptr ) 246 { 247 message << "Regular Stepper instantiation FAILED." << G4endl; 248 message << "G4ChordFinder: Attempted to instantiate " 249 << RegularStepperName << " type stepper " << G4endl; 250 errorInStepperCreation = true; 251 } 252 else 253 { 254 auto dp5= dynamic_cast<G4DormandPrince745*>(regularStepper); 255 if( dp5 != nullptr ) 256 { 257 fIntgrDriver = new G4InterpolationDriver<G4DormandPrince745>( 258 stepMinimum, dp5, nVar6 ); 259 if( gVerboseCtor ) 260 { 261 G4cout << " Using InterpolationDriver<DoPri5> " << G4endl; 262 } 263 } 264 else 265 { 266 fIntgrDriver = new G4IntegrationDriver<RegularStepperType>( 267 stepMinimum, regularStepper, nVar6 ); 268 if( gVerboseCtor ) 269 { 270 G4cout << " Using IntegrationDriver<DoPri5> " << G4endl; 271 } 272 } 273 } 274 } 275 else if ( useBfieldDriver ) 276 { 277 auto regularStepper = new G4DormandPrince745(pEquation); 278 // *** ****************** 279 // 280 fRegularStepperOwned = regularStepper; 281 282 { 283 using SmallStepDriver = G4InterpolationDriver<G4DormandPrince745>; 284 using LargeStepDriver = G4IntegrationDriver<G4HelixHeum>; 285 286 fLongStepper = std::make_unique<G4HelixHeum>(pEquation); 287 288 fIntgrDriver = new G4BFieldIntegrationDriver( 289 std::make_unique<SmallStepDriver>(stepMinimum, 290 regularStepper, regularStepper->GetNumberOfVariables()), 291 std::make_unique<LargeStepDriver>(stepMinimum, 292 fLongStepper.get(), regularStepper->GetNumberOfVariables()) ); 293 294 if( fIntgrDriver == nullptr) 295 { 296 message << "Using G4BFieldIntegrationDriver with " 297 << RegularStepperName << " type stepper " << G4endl; 298 message << "Driver instantiation FAILED." << G4endl; 299 G4Exception("G4ChordFinder::G4ChordFinder()", 300 "GeomField1001", JustWarning, message); 301 } 302 } 303 } 304 else if( useG4QSSDriver ) 305 { 306 if( stepperDriverId == kQss2DriverType ) 307 { 308 auto qssStepper2 = G4QSSDriverCreator::CreateQss2Stepper(pEquation); 309 if( gVerboseCtor ) 310 { 311 G4cout << "-- Created QSS-2 stepper" << G4endl; 312 } 313 fIntgrDriver = G4QSSDriverCreator::CreateDriver(qssStepper2); 314 } 315 else 316 { 317 auto qssStepper3 = G4QSSDriverCreator::CreateQss3Stepper(pEquation); 318 if( gVerboseCtor ) 319 { 320 G4cout << "-- Created QSS-3 stepper" << G4endl; 321 } 322 fIntgrDriver = G4QSSDriverCreator::CreateDriver(qssStepper3); 323 } 324 if( gVerboseCtor ) 325 { 326 G4cout << "-- G4ChordFinder: Using QSS Driver." << G4endl; 327 } 328 } 329 else 330 { 331 auto fsalStepper= new NewFsalStepperType(pEquation); 332 // *** ****************** 333 fNewFSALStepperOwned = fsalStepper; 334 335 if( fsalStepper == nullptr ) 336 { 337 message << "Stepper instantiation FAILED." << G4endl; 338 message << "Attempted to instantiate " 339 << NewFSALStepperName << " type stepper " << G4endl; 340 G4Exception("G4ChordFinder::G4ChordFinder()", 341 "GeomField1001", JustWarning, message); 342 errorInStepperCreation = true; 343 } 344 else 345 { 346 fIntgrDriver = new 347 G4FSALIntegrationDriver<NewFsalStepperType>(stepMinimum, fsalStepper, 348 fsalStepper->GetNumberOfVariables() ); 349 // ==== Create the driver which knows the class type 350 351 if( fIntgrDriver == nullptr ) 352 { 353 message << "Using G4FSALIntegrationDriver with stepper type: " 354 << NewFSALStepperName << G4endl; 355 message << "Integration Driver instantiation FAILED." << G4endl; 356 G4Exception("G4ChordFinder::G4ChordFinder()", 357 "GeomField1001", JustWarning, message); 358 } 359 } 360 } 361 362 // -- Main work is now done 363 364 // Now check that no error occured, and report it if one did. 365 366 // To test failure to create driver 367 // delete fIntgrDriver; 368 // fIntgrDriver = nullptr; 369 370 // Detect and report Error conditions 371 // 372 if( errorInStepperCreation || (fIntgrDriver == nullptr )) 373 { 374 std::ostringstream errmsg; 375 376 if( errorInStepperCreation ) 377 { 378 errmsg << "ERROR> Failure to create Stepper object." << G4endl 379 << " --------------------------------" << G4endl; 380 } 381 if (fIntgrDriver == nullptr ) 382 { 383 errmsg << "ERROR> Failure to create Integration-Driver object." 384 << G4endl 385 << " -------------------------------------------" 386 << G4endl; 387 } 388 const std::string BoolName[2]= { "False", "True" }; 389 errmsg << " Configuration: (constructor arguments) " << G4endl 390 << " provided Stepper = " << pItsStepper << G4endl 391 << " stepper/driver Id = " << stepperDriverId << " i.e. " 392 << " useTemplated = " << BoolName[useTemplatedStepper] 393 << " useRegular = " << BoolName[useRegularStepper] 394 << " useFSAL = " << BoolName[useFSALstepper] 395 << " using combo BField Driver = " << 396 BoolName[ ! (useFSALstepper||useTemplatedStepper 397 || useRegularStepper ) ] 398 << G4endl; 399 errmsg << message.str(); 400 errmsg << "Aborting."; 401 G4Exception("G4ChordFinder::G4ChordFinder() - constructor 2", 402 "GeomField0003", FatalException, errmsg); 403 } 404 405 assert( ( pItsStepper != nullptr ) 406 || ( fRegularStepperOwned != nullptr ) 407 || ( fNewFSALStepperOwned != nullptr ) 408 || useG4QSSDriver 409 ); 410 assert( fIntgrDriver != nullptr ); 411 } 412 413 // ...................................................................... 414 415 G4ChordFinder::~G4ChordFinder() 416 { 417 delete fEquation; 418 delete fRegularStepperOwned; 419 delete fNewFSALStepperOwned; 420 delete fCachedField; 421 delete fIntgrDriver; 422 } 423 424 // ........................................................................... 425 426 G4FieldTrack 427 G4ChordFinder::ApproxCurvePointS( const G4FieldTrack& CurveA_PointVelocity, 428 const G4FieldTrack& CurveB_PointVelocity, 429 const G4FieldTrack& ApproxCurveV, 430 const G4ThreeVector& CurrentE_Point, 431 const G4ThreeVector& CurrentF_Point, 432 const G4ThreeVector& PointG, 433 G4bool first, G4double eps_step) 434 { 435 // ApproxCurvePointS is 2nd implementation of ApproxCurvePoint. 436 // Use Brent Algorithm (or InvParabolic) when possible. 437 // Given a starting curve point A (CurveA_PointVelocity), curve point B 438 // (CurveB_PointVelocity), a point E which is (generally) not on the curve 439 // and a point F which is on the curve (first approximation), find new 440 // point S on the curve closer to point E. 441 // While advancing towards S utilise 'eps_step' as a measure of the 442 // relative accuracy of each Step. 443 444 G4FieldTrack EndPoint(CurveA_PointVelocity); 445 if(!first) { EndPoint = ApproxCurveV; } 446 447 G4ThreeVector Point_A,Point_B; 448 Point_A=CurveA_PointVelocity.GetPosition(); 449 Point_B=CurveB_PointVelocity.GetPosition(); 450 451 G4double xa,xb,xc,ya,yb,yc; 452 453 // InverseParabolic. AF Intersects (First Part of Curve) 454 455 if(first) 456 { 457 xa=0.; 458 ya=(PointG-Point_A).mag(); 459 xb=(Point_A-CurrentF_Point).mag(); 460 yb=-(PointG-CurrentF_Point).mag(); 461 xc=(Point_A-Point_B).mag(); 462 yc=-(CurrentE_Point-Point_B).mag(); 463 } 464 else 465 { 466 xa=0.; 467 ya=(Point_A-CurrentE_Point).mag(); 468 xb=(Point_A-CurrentF_Point).mag(); 469 yb=(PointG-CurrentF_Point).mag(); 470 xc=(Point_A-Point_B).mag(); 471 yc=-(Point_B-PointG).mag(); 472 if(xb==0.) 473 { 474 EndPoint = ApproxCurvePointV(CurveA_PointVelocity, CurveB_PointVelocity, 475 CurrentE_Point, eps_step); 476 return EndPoint; 477 } 478 } 479 480 const G4double tolerance = 1.e-12; 481 if(std::abs(ya)<=tolerance||std::abs(yc)<=tolerance) 482 { 483 ; // What to do for the moment: return the same point as at start 484 // then PropagatorInField will take care 485 } 486 else 487 { 488 G4double test_step = InvParabolic(xa,ya,xb,yb,xc,yc); 489 G4double curve; 490 if(first) 491 { 492 curve=std::abs(EndPoint.GetCurveLength() 493 -ApproxCurveV.GetCurveLength()); 494 } 495 else 496 { 497 test_step = test_step - xb; 498 curve=std::abs(EndPoint.GetCurveLength() 499 -CurveB_PointVelocity.GetCurveLength()); 500 xb = (CurrentF_Point-Point_B).mag(); 501 } 502 503 if(test_step<=0) { test_step=0.1*xb; } 504 if(test_step>=xb) { test_step=0.5*xb; } 505 if(test_step>=curve){ test_step=0.5*curve; } 506 507 if(curve*(1.+eps_step)<xb) // Similar to ReEstimate Step from 508 { // G4VIntersectionLocator 509 test_step=0.5*curve; 510 } 511 512 fIntgrDriver->AccurateAdvance(EndPoint,test_step, eps_step); 513 514 #ifdef G4DEBUG_FIELD 515 // Printing Brent and Linear Approximation 516 // 517 G4cout << "G4ChordFinder::ApproxCurvePointS() - test-step ShF = " 518 << test_step << " EndPoint = " << EndPoint << G4endl; 519 520 // Test Track 521 // 522 G4FieldTrack TestTrack( CurveA_PointVelocity); 523 TestTrack = ApproxCurvePointV( CurveA_PointVelocity, 524 CurveB_PointVelocity, 525 CurrentE_Point, eps_step ); 526 G4cout.precision(14); 527 G4cout << "G4ChordFinder::BrentApprox = " << EndPoint << G4endl; 528 G4cout << "G4ChordFinder::LinearApprox= " << TestTrack << G4endl; 529 #endif 530 } 531 return EndPoint; 532 } 533 534 535 // ........................................................................... 536 537 G4FieldTrack G4ChordFinder:: 538 ApproxCurvePointV( const G4FieldTrack& CurveA_PointVelocity, 539 const G4FieldTrack& CurveB_PointVelocity, 540 const G4ThreeVector& CurrentE_Point, 541 G4double eps_step) 542 { 543 // If r=|AE|/|AB|, and s=true path lenght (AB) 544 // return the point that is r*s along the curve! 545 546 G4FieldTrack Current_PointVelocity = CurveA_PointVelocity; 547 548 G4ThreeVector CurveA_Point= CurveA_PointVelocity.GetPosition(); 549 G4ThreeVector CurveB_Point= CurveB_PointVelocity.GetPosition(); 550 551 G4ThreeVector ChordAB_Vector= CurveB_Point - CurveA_Point; 552 G4ThreeVector ChordAE_Vector= CurrentE_Point - CurveA_Point; 553 554 G4double ABdist= ChordAB_Vector.mag(); 555 G4double curve_length; // A curve length of AB 556 G4double AE_fraction; 557 558 curve_length= CurveB_PointVelocity.GetCurveLength() 559 - CurveA_PointVelocity.GetCurveLength(); 560 561 G4double integrationInaccuracyLimit= std::max( perMillion, 0.5*eps_step ); 562 if( curve_length < ABdist * (1. - integrationInaccuracyLimit) ) 563 { 564 #ifdef G4DEBUG_FIELD 565 G4cerr << " Warning in G4ChordFinder::ApproxCurvePoint: " 566 << G4endl 567 << " The two points are further apart than the curve length " 568 << G4endl 569 << " Dist = " << ABdist 570 << " curve length = " << curve_length 571 << " relativeDiff = " << (curve_length-ABdist)/ABdist 572 << G4endl; 573 if( curve_length < ABdist * (1. - 10*eps_step) ) 574 { 575 std::ostringstream message; 576 message << "Unphysical curve length." << G4endl 577 << "The size of the above difference exceeds allowed limits." 578 << G4endl 579 << "Aborting."; 580 G4Exception("G4ChordFinder::ApproxCurvePointV()", "GeomField0003", 581 FatalException, message); 582 } 583 #endif 584 // Take default corrective action: adjust the maximum curve length. 585 // NOTE: this case only happens for relatively straight paths. 586 // curve_length = ABdist; 587 } 588 589 G4double new_st_length; 590 591 if ( ABdist > 0.0 ) 592 { 593 AE_fraction = ChordAE_Vector.mag() / ABdist; 594 } 595 else 596 { 597 AE_fraction = 0.5; // Guess .. ?; 598 #ifdef G4DEBUG_FIELD 599 G4cout << "Warning in G4ChordFinder::ApproxCurvePointV():" 600 << " A and B are the same point!" << G4endl 601 << " Chord AB length = " << ChordAE_Vector.mag() << G4endl 602 << G4endl; 603 #endif 604 } 605 606 if( (AE_fraction> 1.0 + perMillion) || (AE_fraction< 0.) ) 607 { 608 #ifdef G4DEBUG_FIELD 609 G4cerr << " G4ChordFinder::ApproxCurvePointV() - Warning:" 610 << " Anomalous condition:AE > AB or AE/AB <= 0 " << G4endl 611 << " AE_fraction = " << AE_fraction << G4endl 612 << " Chord AE length = " << ChordAE_Vector.mag() << G4endl 613 << " Chord AB length = " << ABdist << G4endl << G4endl; 614 G4cerr << " OK if this condition occurs after a recalculation of 'B'" 615 << G4endl << " Otherwise it is an error. " << G4endl ; 616 #endif 617 // This course can now result if B has been re-evaluated, 618 // without E being recomputed (1 July 99). 619 // In this case this is not a "real error" - but it is undesired 620 // and we cope with it by a default corrective action ... 621 // 622 AE_fraction = 0.5; // Default value 623 } 624 625 new_st_length = AE_fraction * curve_length; 626 627 if ( AE_fraction > 0.0 ) 628 { 629 fIntgrDriver->AccurateAdvance(Current_PointVelocity, 630 new_st_length, eps_step ); 631 // 632 // In this case it does not matter if it cannot advance the full distance 633 } 634 635 // If there was a memory of the step_length actually required at the start 636 // of the integration Step, this could be re-used ... 637 638 G4cout.precision(14); 639 640 return Current_PointVelocity; 641 } 642 643 // ........................................................................... 644 645 std::ostream& operator<<( std::ostream& os, const G4ChordFinder& cf) 646 { 647 // Dumping the state of G4ChordFinder 648 os << "State of G4ChordFinder : " << std::endl; 649 os << " delta_chord = " << cf.fDeltaChord; 650 os << " Default d_c = " << cf.fDefaultDeltaChord; 651 652 os << " stats-verbose = " << cf.fStatsVerbose; 653 654 return os; 655 } 656