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 #include <iomanip>
00034
00035 #include "G4ChordFinder.hh"
00036 #include "G4SystemOfUnits.hh"
00037 #include "G4MagneticField.hh"
00038 #include "G4Mag_UsualEqRhs.hh"
00039 #include "G4ClassicalRK4.hh"
00040
00041
00042
00043
00044 G4ChordFinder::G4ChordFinder(G4MagInt_Driver* pIntegrationDriver)
00045 : fDefaultDeltaChord( 0.25 * mm ),
00046 fDeltaChord( fDefaultDeltaChord ),
00047 fFirstFraction(0.999), fFractionLast(1.00), fFractionNextEstimate(0.98),
00048 fMultipleRadius(15.0),
00049 fStatsVerbose(0),
00050 fDriversStepper(0),
00051 fAllocatedStepper(false),
00052 fEquation(0),
00053 fTotalNoTrials_FNC(0), fNoCalls_FNC(0), fmaxTrials_FNC(0)
00054 {
00055
00056 fIntgrDriver= pIntegrationDriver;
00057 fAllocatedStepper= false;
00058
00059 fLastStepEstimate_Unconstrained = DBL_MAX;
00060
00061 SetFractions_Last_Next( fFractionLast, fFractionNextEstimate);
00062
00063 }
00064
00065
00066
00067
00068 G4ChordFinder::G4ChordFinder( G4MagneticField* theMagField,
00069 G4double stepMinimum,
00070 G4MagIntegratorStepper* pItsStepper )
00071 : fDefaultDeltaChord( 0.25 * mm ),
00072 fDeltaChord( fDefaultDeltaChord ),
00073 fFirstFraction(0.999), fFractionLast(1.00), fFractionNextEstimate(0.98),
00074 fMultipleRadius(15.0),
00075 fStatsVerbose(0),
00076 fDriversStepper(0),
00077 fAllocatedStepper(false),
00078 fEquation(0),
00079 fTotalNoTrials_FNC(0), fNoCalls_FNC(0), fmaxTrials_FNC(0)
00080 {
00081
00082
00083
00084 G4Mag_EqRhs *pEquation = new G4Mag_UsualEqRhs(theMagField);
00085 fEquation = pEquation;
00086 fLastStepEstimate_Unconstrained = DBL_MAX;
00087
00088
00089 SetFractions_Last_Next( fFractionLast, fFractionNextEstimate);
00090
00091
00092
00093
00094
00095 if( pItsStepper == 0 )
00096 {
00097 pItsStepper = fDriversStepper = new G4ClassicalRK4(pEquation);
00098 fAllocatedStepper= true;
00099 }
00100 else
00101 {
00102 fAllocatedStepper= false;
00103 }
00104 fIntgrDriver = new G4MagInt_Driver(stepMinimum, pItsStepper,
00105 pItsStepper->GetNumberOfVariables() );
00106 }
00107
00108
00109
00110
00111 G4ChordFinder::~G4ChordFinder()
00112 {
00113 delete fEquation;
00114 if( fAllocatedStepper)
00115 {
00116 delete fDriversStepper;
00117 }
00118 delete fIntgrDriver;
00119
00120 if( fStatsVerbose ) { PrintStatistics(); }
00121 }
00122
00123
00124
00125
00126 void
00127 G4ChordFinder::SetFractions_Last_Next( G4double fractLast, G4double fractNext )
00128 {
00129
00130 if( fractLast == -1.0 ) fractLast = 1.0;
00131 if( fractNext == -1.0 ) fractNext = 0.98;
00132
00133
00134
00135
00136 if( fStatsVerbose )
00137 {
00138 G4cout << " ChordFnd> Trying to set fractions: "
00139 << " first " << fFirstFraction
00140 << " last " << fractLast
00141 << " next " << fractNext
00142 << " and multiple " << fMultipleRadius
00143 << G4endl;
00144 }
00145
00146 if( (fractLast > 0.0) && (fractLast <=1.0) )
00147 {
00148 fFractionLast= fractLast;
00149 }
00150 else
00151 {
00152 G4cerr << "G4ChordFinder::SetFractions_Last_Next: Invalid "
00153 << " fraction Last = " << fractLast
00154 << " must be 0 < fractionLast <= 1 " << G4endl;
00155 }
00156 if( (fractNext > 0.0) && (fractNext <1.0) )
00157 {
00158 fFractionNextEstimate = fractNext;
00159 }
00160 else
00161 {
00162 G4cerr << "G4ChordFinder:: SetFractions_Last_Next: Invalid "
00163 << " fraction Next = " << fractNext
00164 << " must be 0 < fractionNext < 1 " << G4endl;
00165 }
00166 }
00167
00168
00169
00170
00171 G4double
00172 G4ChordFinder::AdvanceChordLimited( G4FieldTrack& yCurrent,
00173 G4double stepMax,
00174 G4double epsStep,
00175 const G4ThreeVector latestSafetyOrigin,
00176 G4double latestSafetyRadius )
00177 {
00178 G4double stepPossible;
00179 G4double dyErr;
00180 G4FieldTrack yEnd( yCurrent);
00181 G4double startCurveLen= yCurrent.GetCurveLength();
00182 G4double nextStep;
00183
00184 stepPossible= FindNextChord(yCurrent, stepMax, yEnd, dyErr, epsStep,
00185 &nextStep, latestSafetyOrigin, latestSafetyRadius
00186 );
00187
00188
00189 G4bool good_advance;
00190
00191 if ( dyErr < epsStep * stepPossible )
00192 {
00193
00194
00195 yCurrent = yEnd;
00196 good_advance = true;
00197 }
00198 else
00199 {
00200
00201
00202 good_advance = fIntgrDriver->AccurateAdvance(yCurrent, stepPossible,
00203 epsStep, nextStep);
00204 if ( ! good_advance )
00205 {
00206
00207 stepPossible= yCurrent.GetCurveLength()-startCurveLen;
00208 }
00209 }
00210 return stepPossible;
00211 }
00212
00213
00214
00215
00216 G4double
00217 G4ChordFinder::FindNextChord( const G4FieldTrack& yStart,
00218 G4double stepMax,
00219 G4FieldTrack& yEnd,
00220 G4double& dyErrPos,
00221 G4double epsStep,
00222 G4double* pStepForAccuracy,
00223 const G4ThreeVector,
00224 G4double
00225 )
00226 {
00227
00228
00229 G4FieldTrack yCurrent= yStart;
00230 G4double stepTrial, stepForAccuracy;
00231 G4double dydx[G4FieldTrack::ncompSVEC];
00232
00233
00234
00235
00236
00237 G4bool validEndPoint= false;
00238 G4double dChordStep, lastStepLength;
00239
00240 fIntgrDriver-> GetDerivatives( yCurrent, dydx );
00241
00242 G4int noTrials=0;
00243 const G4double safetyFactor= fFirstFraction;
00244
00245 stepTrial = std::min( stepMax, safetyFactor*fLastStepEstimate_Unconstrained );
00246
00247 G4double newStepEst_Uncons= 0.0;
00248 do
00249 {
00250 G4double stepForChord;
00251 yCurrent = yStart;
00252
00253
00254 fIntgrDriver->QuickAdvance( yCurrent, dydx, stepTrial,
00255 dChordStep, dyErrPos);
00256
00257
00258
00259 validEndPoint = AcceptableMissDist(dChordStep);
00260
00261 lastStepLength = stepTrial;
00262
00263
00264 stepForChord = NewStep(stepTrial, dChordStep, newStepEst_Uncons );
00265
00266 if( ! validEndPoint )
00267 {
00268 if( stepTrial<=0.0 )
00269 {
00270 stepTrial = stepForChord;
00271 }
00272 else if (stepForChord <= stepTrial)
00273 {
00274
00275 stepTrial = std::min( stepForChord, fFractionLast * stepTrial);
00276 }
00277 else
00278 {
00279 stepTrial *= 0.1;
00280 }
00281 }
00282 noTrials++;
00283 }
00284 while( ! validEndPoint );
00285
00286 if( newStepEst_Uncons > 0.0 )
00287 {
00288 fLastStepEstimate_Unconstrained= newStepEst_Uncons;
00289 }
00290
00291 AccumulateStatistics( noTrials );
00292
00293 if( pStepForAccuracy )
00294 {
00295
00296
00297 G4double dyErr_relative = dyErrPos/(epsStep*lastStepLength);
00298 if( dyErr_relative > 1.0 )
00299 {
00300 stepForAccuracy = fIntgrDriver->ComputeNewStepSize( dyErr_relative,
00301 lastStepLength );
00302 }
00303 else
00304 {
00305 stepForAccuracy = 0.0;
00306 }
00307 *pStepForAccuracy = stepForAccuracy;
00308 }
00309
00310 #ifdef TEST_CHORD_PRINT
00311 static int dbg=0;
00312 if( dbg )
00313 {
00314 G4cout << "ChordF/FindNextChord: NoTrials= " << noTrials
00315 << " StepForGoodChord=" << std::setw(10) << stepTrial << G4endl;
00316 }
00317 #endif
00318 yEnd= yCurrent;
00319 return stepTrial;
00320 }
00321
00322
00323
00324
00325 G4double G4ChordFinder::NewStep(G4double stepTrialOld,
00326 G4double dChordStep,
00327 G4double& stepEstimate_Unconstrained )
00328 {
00329
00330
00331
00332
00333 G4double stepTrial;
00334
00335 #if 1
00336
00337 if (dChordStep > 0.0)
00338 {
00339 stepEstimate_Unconstrained =
00340 stepTrialOld*std::sqrt( fDeltaChord / dChordStep );
00341 stepTrial = fFractionNextEstimate * stepEstimate_Unconstrained;
00342 }
00343 else
00344 {
00345
00346 stepTrial = stepTrialOld * 2.;
00347 }
00348
00349 if( stepTrial <= 0.001 * stepTrialOld)
00350 {
00351 if ( dChordStep > 1000.0 * fDeltaChord )
00352 {
00353 stepTrial= stepTrialOld * 0.03;
00354 }
00355 else
00356 {
00357 if ( dChordStep > 100. * fDeltaChord )
00358 {
00359 stepTrial= stepTrialOld * 0.1;
00360 }
00361 else
00362 {
00363 stepTrial= stepTrialOld * 0.5;
00364 }
00365 }
00366 }
00367 else if (stepTrial > 1000.0 * stepTrialOld)
00368 {
00369 stepTrial= 1000.0 * stepTrialOld;
00370 }
00371
00372 if( stepTrial == 0.0 )
00373 {
00374 stepTrial= 0.000001;
00375 }
00376
00377 #else
00378
00379 if ( dChordStep > 1000. * fDeltaChord )
00380 {
00381 stepTrial= stepTrialOld * 0.03;
00382 }
00383 else
00384 {
00385 if ( dChordStep > 100. * fDeltaChord )
00386 {
00387 stepTrial= stepTrialOld * 0.1;
00388 }
00389 else
00390 {
00391 stepTrial= stepTrialOld * 0.5;
00392 }
00393 }
00394
00395 #endif
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406 return stepTrial;
00407 }
00408
00409
00410
00411
00412 G4FieldTrack
00413 G4ChordFinder::ApproxCurvePointS( const G4FieldTrack& CurveA_PointVelocity,
00414 const G4FieldTrack& CurveB_PointVelocity,
00415 const G4FieldTrack& ApproxCurveV,
00416 const G4ThreeVector& CurrentE_Point,
00417 const G4ThreeVector& CurrentF_Point,
00418 const G4ThreeVector& PointG,
00419 G4bool first, G4double eps_step)
00420 {
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430 G4FieldTrack EndPoint(CurveA_PointVelocity);
00431 if(!first){EndPoint= ApproxCurveV;}
00432
00433 G4ThreeVector Point_A,Point_B;
00434 Point_A=CurveA_PointVelocity.GetPosition();
00435 Point_B=CurveB_PointVelocity.GetPosition();
00436
00437 G4double xa,xb,xc,ya,yb,yc;
00438
00439
00440
00441 if(first)
00442 {
00443 xa=0.;
00444 ya=(PointG-Point_A).mag();
00445 xb=(Point_A-CurrentF_Point).mag();
00446 yb=-(PointG-CurrentF_Point).mag();
00447 xc=(Point_A-Point_B).mag();
00448 yc=-(CurrentE_Point-Point_B).mag();
00449 }
00450 else
00451 {
00452 xa=0.;
00453 ya=(Point_A-CurrentE_Point).mag();
00454 xb=(Point_A-CurrentF_Point).mag();
00455 yb=(PointG-CurrentF_Point).mag();
00456 xc=(Point_A-Point_B).mag();
00457 yc=-(Point_B-PointG).mag();
00458 if(xb==0.)
00459 {
00460 EndPoint=
00461 ApproxCurvePointV(CurveA_PointVelocity, CurveB_PointVelocity,
00462 CurrentE_Point, eps_step);
00463 return EndPoint;
00464 }
00465 }
00466
00467 const G4double tolerance= 1.e-12;
00468 if(std::abs(ya)<=tolerance||std::abs(yc)<=tolerance)
00469 {
00470 ;
00471
00472 }
00473 else
00474 {
00475 G4double test_step = InvParabolic(xa,ya,xb,yb,xc,yc);
00476 G4double curve;
00477 if(first)
00478 {
00479 curve=std::abs(EndPoint.GetCurveLength()
00480 -ApproxCurveV.GetCurveLength());
00481 }
00482 else
00483 {
00484 test_step=(test_step-xb);
00485 curve=std::abs(EndPoint.GetCurveLength()
00486 -CurveB_PointVelocity.GetCurveLength());
00487 xb=(CurrentF_Point-Point_B).mag();
00488 }
00489
00490 if(test_step<=0) { test_step=0.1*xb; }
00491 if(test_step>=xb) { test_step=0.5*xb; }
00492 if(test_step>=curve){ test_step=0.5*curve; }
00493
00494 if(curve*(1.+eps_step)<xb)
00495 {
00496 test_step=0.5*curve;
00497 }
00498
00499 fIntgrDriver->AccurateAdvance(EndPoint,test_step, eps_step);
00500
00501 #ifdef G4DEBUG_FIELD
00502
00503
00504 G4cout << "G4ChordFinder::ApproxCurvePointS() - test-step ShF = "
00505 << test_step << " EndPoint = " << EndPoint << G4endl;
00506
00507
00508
00509 G4FieldTrack TestTrack( CurveA_PointVelocity);
00510 TestTrack = ApproxCurvePointV( CurveA_PointVelocity,
00511 CurveB_PointVelocity,
00512 CurrentE_Point, eps_step );
00513 G4cout.precision(14);
00514 G4cout << "G4ChordFinder::BrentApprox = " << EndPoint << G4endl;
00515 G4cout << "G4ChordFinder::LinearApprox= " << TestTrack << G4endl;
00516 #endif
00517 }
00518 return EndPoint;
00519 }
00520
00521
00522
00523
00524 G4FieldTrack G4ChordFinder::
00525 ApproxCurvePointV( const G4FieldTrack& CurveA_PointVelocity,
00526 const G4FieldTrack& CurveB_PointVelocity,
00527 const G4ThreeVector& CurrentE_Point,
00528 G4double eps_step)
00529 {
00530
00531
00532
00533 G4FieldTrack Current_PointVelocity = CurveA_PointVelocity;
00534
00535 G4ThreeVector CurveA_Point= CurveA_PointVelocity.GetPosition();
00536 G4ThreeVector CurveB_Point= CurveB_PointVelocity.GetPosition();
00537
00538 G4ThreeVector ChordAB_Vector= CurveB_Point - CurveA_Point;
00539 G4ThreeVector ChordAE_Vector= CurrentE_Point - CurveA_Point;
00540
00541 G4double ABdist= ChordAB_Vector.mag();
00542 G4double curve_length;
00543 G4double AE_fraction;
00544
00545 curve_length= CurveB_PointVelocity.GetCurveLength()
00546 - CurveA_PointVelocity.GetCurveLength();
00547
00548 G4double integrationInaccuracyLimit= std::max( perMillion, 0.5*eps_step );
00549 if( curve_length < ABdist * (1. - integrationInaccuracyLimit) )
00550 {
00551 #ifdef G4DEBUG_FIELD
00552 G4cerr << " Warning in G4ChordFinder::ApproxCurvePoint: "
00553 << G4endl
00554 << " The two points are further apart than the curve length "
00555 << G4endl
00556 << " Dist = " << ABdist
00557 << " curve length = " << curve_length
00558 << " relativeDiff = " << (curve_length-ABdist)/ABdist
00559 << G4endl;
00560 if( curve_length < ABdist * (1. - 10*eps_step) )
00561 {
00562 std::ostringstream message;
00563 message << "Unphysical curve length." << G4endl
00564 << "The size of the above difference exceeds allowed limits."
00565 << G4endl
00566 << "Aborting.";
00567 G4Exception("G4ChordFinder::ApproxCurvePointV()", "GeomField0003",
00568 FatalException, message);
00569 }
00570 #endif
00571
00572
00573
00574 }
00575
00576 G4double new_st_length;
00577
00578 if ( ABdist > 0.0 )
00579 {
00580 AE_fraction = ChordAE_Vector.mag() / ABdist;
00581 }
00582 else
00583 {
00584 AE_fraction = 0.5;
00585 #ifdef G4DEBUG_FIELD
00586 G4cout << "Warning in G4ChordFinder::ApproxCurvePointV():"
00587 << " A and B are the same point!" << G4endl
00588 << " Chord AB length = " << ChordAE_Vector.mag() << G4endl
00589 << G4endl;
00590 #endif
00591 }
00592
00593 if( (AE_fraction> 1.0 + perMillion) || (AE_fraction< 0.) )
00594 {
00595 #ifdef G4DEBUG_FIELD
00596 G4cerr << " G4ChordFinder::ApproxCurvePointV() - Warning:"
00597 << " Anomalous condition:AE > AB or AE/AB <= 0 " << G4endl
00598 << " AE_fraction = " << AE_fraction << G4endl
00599 << " Chord AE length = " << ChordAE_Vector.mag() << G4endl
00600 << " Chord AB length = " << ABdist << G4endl << G4endl;
00601 G4cerr << " OK if this condition occurs after a recalculation of 'B'"
00602 << G4endl << " Otherwise it is an error. " << G4endl ;
00603 #endif
00604
00605
00606
00607
00608
00609 AE_fraction = 0.5;
00610 }
00611
00612 new_st_length= AE_fraction * curve_length;
00613
00614 if ( AE_fraction > 0.0 )
00615 {
00616 fIntgrDriver->AccurateAdvance(Current_PointVelocity,
00617 new_st_length, eps_step );
00618
00619
00620 }
00621
00622
00623
00624
00625 G4cout.precision(14);
00626
00627 return Current_PointVelocity;
00628 }
00629
00630
00631
00632
00633 void
00634 G4ChordFinder::PrintStatistics()
00635 {
00636
00637
00638 G4cout << "G4ChordFinder statistics report: " << G4endl;
00639 G4cout
00640 << " No trials: " << fTotalNoTrials_FNC
00641 << " No Calls: " << fNoCalls_FNC
00642 << " Max-trial: " << fmaxTrials_FNC
00643 << G4endl;
00644 G4cout
00645 << " Parameters: "
00646 << " fFirstFraction " << fFirstFraction
00647 << " fFractionLast " << fFractionLast
00648 << " fFractionNextEstimate " << fFractionNextEstimate
00649 << G4endl;
00650 }
00651
00652
00653
00654
00655 void G4ChordFinder::TestChordPrint( G4int noTrials,
00656 G4int lastStepTrial,
00657 G4double dChordStep,
00658 G4double nextStepTrial )
00659 {
00660 G4int oldprec= G4cout.precision(5);
00661 G4cout << " ChF/fnc: notrial " << std::setw( 3) << noTrials
00662 << " this_step= " << std::setw(10) << lastStepTrial;
00663 if( std::fabs( (dChordStep / fDeltaChord) - 1.0 ) < 0.001 )
00664 {
00665 G4cout.precision(8);
00666 }
00667 else
00668 {
00669 G4cout.precision(6);
00670 }
00671 G4cout << " dChordStep= " << std::setw(12) << dChordStep;
00672 if( dChordStep > fDeltaChord ) { G4cout << " d+"; }
00673 else { G4cout << " d-"; }
00674 G4cout.precision(5);
00675 G4cout << " new_step= " << std::setw(10)
00676 << fLastStepEstimate_Unconstrained
00677 << " new_step_constr= " << std::setw(10)
00678 << lastStepTrial << G4endl;
00679 G4cout << " nextStepTrial = " << std::setw(10) << nextStepTrial << G4endl;
00680 G4cout.precision(oldprec);
00681 }