00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042 #include <iomanip>
00043
00044 #include "G4PropagatorInField.hh"
00045 #include "G4ios.hh"
00046 #include "G4SystemOfUnits.hh"
00047 #include "G4ThreeVector.hh"
00048 #include "G4VPhysicalVolume.hh"
00049 #include "G4Navigator.hh"
00050 #include "G4GeometryTolerance.hh"
00051 #include "G4VCurvedTrajectoryFilter.hh"
00052 #include "G4ChordFinder.hh"
00053 #include "G4MultiLevelLocator.hh"
00054
00056
00057
00058
00059 G4PropagatorInField::G4PropagatorInField( G4Navigator *theNavigator,
00060 G4FieldManager *detectorFieldMgr,
00061 G4VIntersectionLocator *vLocator )
00062 :
00063 fMax_loop_count(1000),
00064 fUseSafetyForOptimisation(true),
00065 fZeroStepThreshold( 0.0 ),
00066 fDetectorFieldMgr(detectorFieldMgr),
00067 fpTrajectoryFilter( 0 ),
00068 fNavigator(theNavigator),
00069 fCurrentFieldMgr(detectorFieldMgr),
00070 fSetFieldMgr(false),
00071 fCharge(0.0), fInitialMomentumModulus(0.0), fMass(0.0),
00072 End_PointAndTangent(G4ThreeVector(0.,0.,0.),
00073 G4ThreeVector(0.,0.,0.),0.0,0.0,0.0,0.0,0.0),
00074 fParticleIsLooping(false),
00075 fNoZeroStep(0),
00076 fVerboseLevel(0)
00077 {
00078 if(fDetectorFieldMgr) { fEpsilonStep = fDetectorFieldMgr->GetMaximumEpsilonStep();}
00079 else { fEpsilonStep= 1.0e-5; }
00080 fActionThreshold_NoZeroSteps = 2;
00081 fSevereActionThreshold_NoZeroSteps = 10;
00082 fAbandonThreshold_NoZeroSteps = 50;
00083 fFull_CurveLen_of_LastAttempt = -1;
00084 fLast_ProposedStepLength = -1;
00085 fLargestAcceptableStep = 1000.0 * meter;
00086
00087 fPreviousSftOrigin= G4ThreeVector(0.,0.,0.);
00088 fPreviousSafety= 0.0;
00089 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
00090 fZeroStepThreshold= std::max( 1.0e5 * kCarTolerance, 1.0e-1 * micrometer );
00091
00092 #ifdef G4DEBUG_FIELD
00093 G4cout << " PiF: Zero Step Threshold set to "
00094 << fZeroStepThreshold / millimeter
00095 << " mm." << G4endl;
00096 G4cout << " PiF: Value of kCarTolerance = "
00097 << kCarTolerance / millimeter
00098 << " mm. " << G4endl;
00099 #endif
00100
00101
00102 if(vLocator==0){
00103 fIntersectionLocator= new G4MultiLevelLocator(theNavigator);
00104 fAllocatedLocator=true;
00105 }else{
00106 fIntersectionLocator=vLocator;
00107 fAllocatedLocator=false;
00108 }
00109 RefreshIntersectionLocator();
00110 }
00111
00112 G4PropagatorInField::~G4PropagatorInField()
00113 {
00114 if(fAllocatedLocator)delete fIntersectionLocator;
00115 }
00116
00117
00118 void
00119 G4PropagatorInField::RefreshIntersectionLocator()
00120 {
00121 fIntersectionLocator->SetEpsilonStepFor(fEpsilonStep);
00122 fIntersectionLocator->SetDeltaIntersectionFor(fCurrentFieldMgr->GetDeltaIntersection());
00123 fIntersectionLocator->SetChordFinderFor(GetChordFinder());
00124 fIntersectionLocator->SetSafetyParametersFor( fUseSafetyForOptimisation);
00125 }
00127
00128
00129
00130 G4double
00131 G4PropagatorInField::ComputeStep(
00132 G4FieldTrack& pFieldTrack,
00133 G4double CurrentProposedStepLength,
00134 G4double& currentSafety,
00135 G4VPhysicalVolume* pPhysVol)
00136 {
00137
00138
00139
00140 if(CurrentProposedStepLength<kCarTolerance)
00141 {
00142 return kInfinity;
00143 }
00144
00145
00146
00147 if (fpTrajectoryFilter)
00148 {
00149 fpTrajectoryFilter->CreateNewTrajectorySegment();
00150 }
00151
00152
00153
00154 G4double h_TrialStepSize;
00155 G4double TruePathLength = CurrentProposedStepLength;
00156 G4double StepTaken = 0.0;
00157 G4double s_length_taken, epsilon ;
00158 G4bool intersects;
00159 G4bool first_substep = true;
00160
00161 G4double NewSafety;
00162 fParticleIsLooping = false;
00163
00164
00165
00166
00167
00168 if( !fSetFieldMgr ) fCurrentFieldMgr= FindAndSetFieldManager( pPhysVol );
00169
00170 fSetFieldMgr= false;
00171
00172 GetChordFinder()->SetChargeMomentumMass(fCharge,fInitialMomentumModulus,fMass);
00173
00174
00175
00176 RefreshIntersectionLocator();
00177
00178 G4FieldTrack CurrentState(pFieldTrack);
00179 G4FieldTrack OriginalState = CurrentState;
00180
00181
00182
00183
00184 if( CurrentProposedStepLength >= fLargestAcceptableStep )
00185 {
00186 G4ThreeVector StartPointA, VelocityUnit;
00187 StartPointA = pFieldTrack.GetPosition();
00188 VelocityUnit = pFieldTrack.GetMomentumDir();
00189
00190 G4double trialProposedStep = 1.e2 * ( 10.0 * cm +
00191 fNavigator->GetWorldVolume()->GetLogicalVolume()->
00192 GetSolid()->DistanceToOut(StartPointA, VelocityUnit) );
00193 CurrentProposedStepLength= std::min( trialProposedStep,
00194 fLargestAcceptableStep );
00195 }
00196 epsilon = fCurrentFieldMgr->GetDeltaOneStep() / CurrentProposedStepLength;
00197
00198 G4double epsilonMin= fCurrentFieldMgr->GetMinimumEpsilonStep();
00199 G4double epsilonMax= fCurrentFieldMgr->GetMaximumEpsilonStep();;
00200 if( epsilon < epsilonMin ) epsilon = epsilonMin;
00201 if( epsilon > epsilonMax ) epsilon = epsilonMax;
00202 SetEpsilonStep( epsilon );
00203
00204
00205
00206
00207
00208
00209 if( fNoZeroStep > fActionThreshold_NoZeroSteps )
00210 {
00211 G4double stepTrial;
00212
00213 stepTrial= fFull_CurveLen_of_LastAttempt;
00214 if( (stepTrial <= 0.0) && (fLast_ProposedStepLength > 0.0) )
00215 stepTrial= fLast_ProposedStepLength;
00216
00217 G4double decreaseFactor = 0.9;
00218 if( (fNoZeroStep < fSevereActionThreshold_NoZeroSteps)
00219 && (stepTrial > 100.0*fZeroStepThreshold) )
00220 {
00221
00222
00223 decreaseFactor= 0.25;
00224 }
00225 else
00226 {
00227
00228
00229
00230
00231 if( stepTrial > 100.0*fZeroStepThreshold )
00232 decreaseFactor = 0.35;
00233 else if( stepTrial > 30.0*fZeroStepThreshold )
00234 decreaseFactor= 0.5;
00235 else if( stepTrial > 10.0*fZeroStepThreshold )
00236 decreaseFactor= 0.75;
00237 else
00238 decreaseFactor= 0.9;
00239 }
00240 stepTrial *= decreaseFactor;
00241
00242 #ifdef G4DEBUG_FIELD
00243 G4cout << " G4PropagatorInField::ComputeStep(): " << G4endl
00244 << " Decreasing step - "
00245 << " decreaseFactor= " << std::setw(8) << decreaseFactor
00246 << " stepTrial = " << std::setw(18) << stepTrial << " "
00247 << " fZeroStepThreshold = " << fZeroStepThreshold << G4endl;
00248 PrintStepLengthDiagnostic(CurrentProposedStepLength, decreaseFactor,
00249 stepTrial, pFieldTrack);
00250 #endif
00251 if( stepTrial == 0.0 )
00252 {
00253 std::ostringstream message;
00254 message << "Particle abandoned due to lack of progress in field."
00255 << G4endl
00256 << " Properties : " << pFieldTrack << G4endl
00257 << " Attempting a zero step = " << stepTrial << G4endl
00258 << " while attempting to progress after " << fNoZeroStep
00259 << " trial steps. Will abandon step.";
00260 G4Exception("G4PropagatorInField::ComputeStep()", "GeomNav1002",
00261 JustWarning, message);
00262 fParticleIsLooping= true;
00263 return 0;
00264 }
00265 if( stepTrial < CurrentProposedStepLength )
00266 CurrentProposedStepLength = stepTrial;
00267 }
00268 fLast_ProposedStepLength = CurrentProposedStepLength;
00269
00270 G4int do_loop_count = 0;
00271 do
00272 {
00273 G4FieldTrack SubStepStartState = CurrentState;
00274 G4ThreeVector SubStartPoint = CurrentState.GetPosition();
00275
00276 if( !first_substep) {
00277 fNavigator->LocateGlobalPointWithinVolume( SubStartPoint );
00278 }
00279
00280
00281
00282 h_TrialStepSize = CurrentProposedStepLength - StepTaken;
00283
00284
00285
00286 s_length_taken = GetChordFinder()->AdvanceChordLimited(
00287 CurrentState,
00288 h_TrialStepSize,
00289 fEpsilonStep,
00290 fPreviousSftOrigin,
00291 fPreviousSafety
00292 );
00293
00294
00295 fFull_CurveLen_of_LastAttempt = s_length_taken;
00296
00297 G4ThreeVector EndPointB = CurrentState.GetPosition();
00298 G4ThreeVector InterSectionPointE;
00299 G4double LinearStepLength;
00300
00301
00302 intersects= IntersectChord( SubStartPoint, EndPointB,
00303 NewSafety, LinearStepLength,
00304 InterSectionPointE );
00305
00306
00307
00308 if( first_substep ) {
00309 currentSafety = NewSafety;
00310 }
00311
00312 if( intersects )
00313 {
00314 G4FieldTrack IntersectPointVelct_G(CurrentState);
00315
00316
00317
00318 G4bool recalculatedEndPt= false;
00319
00320 G4bool found_intersection = fIntersectionLocator->
00321 EstimateIntersectionPoint( SubStepStartState, CurrentState,
00322 InterSectionPointE, IntersectPointVelct_G,
00323 recalculatedEndPt,fPreviousSafety,fPreviousSftOrigin);
00324 intersects = intersects && found_intersection;
00325 if( found_intersection ) {
00326 End_PointAndTangent= IntersectPointVelct_G;
00327 StepTaken = TruePathLength = IntersectPointVelct_G.GetCurveLength()
00328 - OriginalState.GetCurveLength();
00329 } else {
00330
00331 if( recalculatedEndPt ){
00332 CurrentState= IntersectPointVelct_G;
00333 }
00334 }
00335 }
00336 if( !intersects )
00337 {
00338 StepTaken += s_length_taken;
00339
00340 if (fpTrajectoryFilter) {
00341 fpTrajectoryFilter->TakeIntermediatePoint(CurrentState.GetPosition());
00342 }
00343 }
00344 first_substep = false;
00345
00346 #ifdef G4DEBUG_FIELD
00347 if( fNoZeroStep > fActionThreshold_NoZeroSteps ) {
00348 printStatus( SubStepStartState,
00349 CurrentState, CurrentProposedStepLength,
00350 NewSafety, do_loop_count, pPhysVol );
00351 }
00352 if( (fVerboseLevel > 1) && (do_loop_count > fMax_loop_count-10 )) {
00353 if( do_loop_count == fMax_loop_count-9 ){
00354 G4cout << " G4PropagatorInField::ComputeStep(): " << G4endl
00355 << " Difficult track - taking many sub steps." << G4endl;
00356 }
00357 printStatus( SubStepStartState, CurrentState, CurrentProposedStepLength,
00358 NewSafety, do_loop_count, pPhysVol );
00359 }
00360 #endif
00361
00362 do_loop_count++;
00363
00364 } while( (!intersects )
00365 && (StepTaken + kCarTolerance < CurrentProposedStepLength)
00366 && ( do_loop_count < fMax_loop_count ) );
00367
00368 if( do_loop_count >= fMax_loop_count )
00369 {
00370 fParticleIsLooping = true;
00371
00372 if ( fVerboseLevel > 0 )
00373 {
00374 G4cout << " G4PropagateInField::ComputeStep(): " << G4endl
00375 << " Killing looping particle "
00376
00377 << " after " << do_loop_count << " field substeps "
00378 << " totaling " << StepTaken / mm << " mm " ;
00379 if( pPhysVol )
00380 G4cout << " in volume " << pPhysVol->GetName() ;
00381 else
00382 G4cout << " in unknown or null volume. " ;
00383 G4cout << G4endl;
00384 }
00385 }
00386
00387 if( !intersects )
00388 {
00389
00390
00391
00392 End_PointAndTangent = CurrentState;
00393 TruePathLength = StepTaken;
00394 }
00395
00396
00397
00398 pFieldTrack = End_PointAndTangent;
00399
00400 #ifdef G4VERBOSE
00401
00402
00403 if( std::fabs(OriginalState.GetCurveLength() + TruePathLength
00404 - End_PointAndTangent.GetCurveLength()) > 3.e-4 * TruePathLength )
00405 {
00406 std::ostringstream message;
00407 message << "Curve length mis-match between original state "
00408 << "and proposed endpoint of propagation." << G4endl
00409 << " The curve length of the endpoint should be: "
00410 << OriginalState.GetCurveLength() + TruePathLength << G4endl
00411 << " and it is instead: "
00412 << End_PointAndTangent.GetCurveLength() << "." << G4endl
00413 << " A difference of: "
00414 << OriginalState.GetCurveLength() + TruePathLength
00415 - End_PointAndTangent.GetCurveLength() << G4endl
00416 << " Original state = " << OriginalState << G4endl
00417 << " Proposed state = " << End_PointAndTangent;
00418 G4Exception("G4PropagatorInField::ComputeStep()",
00419 "GeomNav0003", FatalException, message);
00420 }
00421 #endif
00422
00423
00424
00425
00426
00427 if( ( (TruePathLength < fZeroStepThreshold)
00428 && ( TruePathLength+kCarTolerance < CurrentProposedStepLength )
00429 )
00430 || ( TruePathLength < 0.5*kCarTolerance )
00431 )
00432 {
00433 fNoZeroStep++;
00434 }
00435 else{
00436 fNoZeroStep = 0;
00437 }
00438
00439 if( fNoZeroStep > fAbandonThreshold_NoZeroSteps )
00440 {
00441 fParticleIsLooping = true;
00442 std::ostringstream message;
00443 message << "Particle is stuck; it will be killed." << G4endl
00444 << " Zero progress for " << fNoZeroStep << " attempted steps."
00445 << G4endl
00446 << " Proposed Step is " << CurrentProposedStepLength
00447 << " but Step Taken is "<< fFull_CurveLen_of_LastAttempt << G4endl
00448 << " For Particle with Charge = " << fCharge
00449 << " Momentum = "<< fInitialMomentumModulus
00450 << " Mass = " << fMass << G4endl;
00451 if( pPhysVol )
00452 message << " in volume " << pPhysVol->GetName() ;
00453 else
00454 message << " in unknown or null volume. " ;
00455 G4Exception("G4PropagatorInField::ComputeStep()",
00456 "GeomNav1002", JustWarning, message);
00457 fNoZeroStep = 0;
00458 }
00459
00460 return TruePathLength;
00461 }
00462
00464
00465
00466
00467 void
00468 G4PropagatorInField::printStatus( const G4FieldTrack& StartFT,
00469 const G4FieldTrack& CurrentFT,
00470 G4double requestStep,
00471 G4double safety,
00472 G4int stepNo,
00473 G4VPhysicalVolume* startVolume)
00474 {
00475 const G4int verboseLevel=fVerboseLevel;
00476 const G4ThreeVector StartPosition = StartFT.GetPosition();
00477 const G4ThreeVector StartUnitVelocity = StartFT.GetMomentumDir();
00478 const G4ThreeVector CurrentPosition = CurrentFT.GetPosition();
00479 const G4ThreeVector CurrentUnitVelocity = CurrentFT.GetMomentumDir();
00480
00481 G4double step_len = CurrentFT.GetCurveLength() - StartFT.GetCurveLength();
00482
00483 G4int oldprec;
00484
00485 if( ((stepNo == 0) && (verboseLevel <3)) || (verboseLevel >= 3) )
00486 {
00487 oldprec = G4cout.precision(4);
00488 G4cout << std::setw( 6) << " "
00489 << std::setw( 25) << " Current Position and Direction" << " "
00490 << G4endl;
00491 G4cout << std::setw( 5) << "Step#"
00492 << std::setw(10) << " s " << " "
00493 << std::setw(10) << "X(mm)" << " "
00494 << std::setw(10) << "Y(mm)" << " "
00495 << std::setw(10) << "Z(mm)" << " "
00496 << std::setw( 7) << " N_x " << " "
00497 << std::setw( 7) << " N_y " << " "
00498 << std::setw( 7) << " N_z " << " " ;
00499 G4cout << std::setw( 7) << " Delta|N|" << " "
00500 << std::setw( 9) << "StepLen" << " "
00501 << std::setw(12) << "StartSafety" << " "
00502 << std::setw( 9) << "PhsStep" << " ";
00503 if( startVolume )
00504 { G4cout << std::setw(18) << "NextVolume" << " "; }
00505 G4cout.precision(oldprec);
00506 G4cout << G4endl;
00507 }
00508 if((stepNo == 0) && (verboseLevel <=3))
00509 {
00510
00511
00512 printStatus( StartFT, StartFT, -1.0, safety, -1, startVolume);
00513 }
00514 if( verboseLevel <= 3 )
00515 {
00516 if( stepNo >= 0)
00517 { G4cout << std::setw( 4) << stepNo << " "; }
00518 else
00519 { G4cout << std::setw( 5) << "Start" ; }
00520 oldprec = G4cout.precision(8);
00521 G4cout << std::setw(10) << CurrentFT.GetCurveLength() << " ";
00522 G4cout.precision(8);
00523 G4cout << std::setw(10) << CurrentPosition.x() << " "
00524 << std::setw(10) << CurrentPosition.y() << " "
00525 << std::setw(10) << CurrentPosition.z() << " ";
00526 G4cout.precision(4);
00527 G4cout << std::setw( 7) << CurrentUnitVelocity.x() << " "
00528 << std::setw( 7) << CurrentUnitVelocity.y() << " "
00529 << std::setw( 7) << CurrentUnitVelocity.z() << " ";
00530 G4cout.precision(3);
00531 G4cout << std::setw( 7)
00532 << CurrentFT.GetMomentum().mag()-StartFT.GetMomentum().mag() << " ";
00533 G4cout << std::setw( 9) << step_len << " ";
00534 G4cout << std::setw(12) << safety << " ";
00535 if( requestStep != -1.0 )
00536 { G4cout << std::setw( 9) << requestStep << " "; }
00537 else
00538 { G4cout << std::setw( 9) << "Init/NotKnown" << " "; }
00539 if( startVolume != 0)
00540 { G4cout << std::setw(12) << startVolume->GetName() << " "; }
00541 G4cout.precision(oldprec);
00542 G4cout << G4endl;
00543 }
00544 else
00545 {
00546
00547
00548 G4cout << "Step taken was " << step_len
00549 << " out of PhysicalStep = " << requestStep << G4endl;
00550 G4cout << "Final safety is: " << safety << G4endl;
00551 G4cout << "Chord length = " << (CurrentPosition-StartPosition).mag()
00552 << G4endl;
00553 G4cout << G4endl;
00554 }
00555 }
00556
00558
00559
00560
00561 void
00562 G4PropagatorInField::PrintStepLengthDiagnostic(
00563 G4double CurrentProposedStepLength,
00564 G4double decreaseFactor,
00565 G4double stepTrial,
00566 const G4FieldTrack& )
00567 {
00568 G4int iprec= G4cout.precision(8);
00569 G4cout << " " << std::setw(12) << " PiF: NoZeroStep "
00570 << " " << std::setw(20) << " CurrentProposed len "
00571 << " " << std::setw(18) << " Full_curvelen_last"
00572 << " " << std::setw(18) << " last proposed len "
00573 << " " << std::setw(18) << " decrease factor "
00574 << " " << std::setw(15) << " step trial "
00575 << G4endl;
00576
00577 G4cout << " " << std::setw(10) << fNoZeroStep << " "
00578 << " " << std::setw(20) << CurrentProposedStepLength
00579 << " " << std::setw(18) << fFull_CurveLen_of_LastAttempt
00580 << " " << std::setw(18) << fLast_ProposedStepLength
00581 << " " << std::setw(18) << decreaseFactor
00582 << " " << std::setw(15) << stepTrial
00583 << G4endl;
00584 G4cout.precision( iprec );
00585
00586 }
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597 std::vector<G4ThreeVector>*
00598 G4PropagatorInField::GimmeTrajectoryVectorAndForgetIt() const
00599 {
00600
00601
00602
00603 if (fpTrajectoryFilter)
00604 {
00605 return fpTrajectoryFilter->GimmeThePointsAndForgetThem();
00606 }
00607 else
00608 {
00609 return 0;
00610 }
00611 }
00612
00613 void
00614 G4PropagatorInField::SetTrajectoryFilter(G4VCurvedTrajectoryFilter* filter)
00615 {
00616 fpTrajectoryFilter = filter;
00617 }
00618
00619 void G4PropagatorInField::ClearPropagatorState()
00620 {
00621
00622
00623 fParticleIsLooping= false;
00624 fNoZeroStep= 0;
00625
00626 End_PointAndTangent= G4FieldTrack( G4ThreeVector(0.,0.,0.),
00627 G4ThreeVector(0.,0.,0.),
00628 0.0,0.0,0.0,0.0,0.0);
00629 fFull_CurveLen_of_LastAttempt = -1;
00630 fLast_ProposedStepLength = -1;
00631
00632 fPreviousSftOrigin= G4ThreeVector(0.,0.,0.);
00633 fPreviousSafety= 0.0;
00634 }
00635
00636 G4FieldManager* G4PropagatorInField::
00637 FindAndSetFieldManager( G4VPhysicalVolume* pCurrentPhysicalVolume)
00638 {
00639 G4FieldManager* currentFieldMgr;
00640
00641 currentFieldMgr = fDetectorFieldMgr;
00642 if( pCurrentPhysicalVolume)
00643 {
00644 G4FieldManager *pRegionFieldMgr= 0, *localFieldMgr = 0;
00645 G4LogicalVolume* pLogicalVol= pCurrentPhysicalVolume->GetLogicalVolume();
00646
00647 if( pLogicalVol ) {
00648
00649 G4Region* pRegion= pLogicalVol->GetRegion();
00650 if( pRegion ) {
00651 pRegionFieldMgr= pRegion->GetFieldManager();
00652 if( pRegionFieldMgr )
00653 currentFieldMgr= pRegionFieldMgr;
00654 }
00655
00656
00657 localFieldMgr= pLogicalVol->GetFieldManager();
00658 if ( localFieldMgr )
00659 currentFieldMgr = localFieldMgr;
00660 }
00661 }
00662 fCurrentFieldMgr= currentFieldMgr;
00663
00664
00665 fSetFieldMgr= true;
00666
00667 return currentFieldMgr;
00668 }
00669
00670 G4int G4PropagatorInField::SetVerboseLevel( G4int level )
00671 {
00672 G4int oldval= fVerboseLevel;
00673 fVerboseLevel= level;
00674
00675
00676
00677
00678 G4MagInt_Driver* integrDriver= GetChordFinder()->GetIntegrationDriver();
00679 integrDriver->SetVerboseLevel( fVerboseLevel - 2 );
00680 G4cout << "Set Driver verbosity to " << fVerboseLevel - 2 << G4endl;
00681
00682 return oldval;
00683 }