136{
137
138
139
140 if(CurrentProposedStepLength<kCarTolerance)
141 {
142 return kInfinity;
143 }
144
145
146
147 if (fpTrajectoryFilter)
148 {
150 }
151
152
153
155 G4double TruePathLength = CurrentProposedStepLength;
159 G4bool first_substep =
true;
160
162 fParticleIsLooping = false;
163
164
165
166
167
169
170 fSetFieldMgr= false;
171
173
174
175
177
180
181
182
183
184 if( CurrentProposedStepLength >= fLargestAcceptableStep )
185 {
189
190 G4double trialProposedStep = 1.e2 * ( 10.0 * cm +
192 GetSolid()->DistanceToOut(StartPointA, VelocityUnit) );
193 CurrentProposedStepLength= std::min( trialProposedStep,
194 fLargestAcceptableStep );
195 }
196 epsilon = fCurrentFieldMgr->
GetDeltaOneStep() / CurrentProposedStepLength;
197
200 if( epsilon < epsilonMin ) epsilon = epsilonMin;
201 if( epsilon > epsilonMax ) epsilon = epsilonMax;
203
204
205
206
207
208
209 if( fNoZeroStep > fActionThreshold_NoZeroSteps )
210 {
212
213 stepTrial= fFull_CurveLen_of_LastAttempt;
214 if( (stepTrial <= 0.0) && (fLast_ProposedStepLength > 0.0) )
215 stepTrial= fLast_ProposedStepLength;
216
218 if( (fNoZeroStep < fSevereActionThreshold_NoZeroSteps)
219 && (stepTrial > 100.0*fZeroStepThreshold) )
220 {
221
222
223 decreaseFactor= 0.25;
224 }
225 else
226 {
227
228
229
230
231 if( stepTrial > 100.0*fZeroStepThreshold )
232 decreaseFactor = 0.35;
233 else if( stepTrial > 30.0*fZeroStepThreshold )
234 decreaseFactor= 0.5;
235 else if( stepTrial > 10.0*fZeroStepThreshold )
236 decreaseFactor= 0.75;
237 else
238 decreaseFactor= 0.9;
239 }
240 stepTrial *= decreaseFactor;
241
242#ifdef G4DEBUG_FIELD
243 G4cout <<
" G4PropagatorInField::ComputeStep(): " <<
G4endl
244 << " Decreasing step - "
245 << " decreaseFactor= " << std::setw(8) << decreaseFactor
246 << " stepTrial = " << std::setw(18) << stepTrial << " "
247 <<
" fZeroStepThreshold = " << fZeroStepThreshold <<
G4endl;
249 stepTrial, pFieldTrack);
250#endif
251 if( stepTrial == 0.0 )
252 {
253 std::ostringstream message;
254 message << "Particle abandoned due to lack of progress in field."
256 <<
" Properties : " << pFieldTrack <<
G4endl
257 <<
" Attempting a zero step = " << stepTrial <<
G4endl
258 << " while attempting to progress after " << fNoZeroStep
259 << " trial steps. Will abandon step.";
260 G4Exception(
"G4PropagatorInField::ComputeStep()",
"GeomNav1002",
262 fParticleIsLooping= true;
263 return 0;
264 }
265 if( stepTrial < CurrentProposedStepLength )
266 CurrentProposedStepLength = stepTrial;
267 }
268 fLast_ProposedStepLength = CurrentProposedStepLength;
269
270 G4int do_loop_count = 0;
271 do
272 {
275
276 if( !first_substep) {
278 }
279
280
281
282 h_TrialStepSize = CurrentProposedStepLength - StepTaken;
283
284
285
287 CurrentState,
288 h_TrialStepSize,
289 fEpsilonStep,
290 fPreviousSftOrigin,
291 fPreviousSafety
292 );
293
294
295 fFull_CurveLen_of_LastAttempt = s_length_taken;
296
300
301
303 NewSafety, LinearStepLength,
304 InterSectionPointE );
305
306
307
308 if( first_substep ) {
309 currentSafety = NewSafety;
310 }
311
312 if( intersects )
313 {
315
316
317
318 G4bool recalculatedEndPt=
false;
319
320 G4bool found_intersection = fIntersectionLocator->
321 EstimateIntersectionPoint( SubStepStartState, CurrentState,
322 InterSectionPointE, IntersectPointVelct_G,
323 recalculatedEndPt,fPreviousSafety,fPreviousSftOrigin);
324 intersects = intersects && found_intersection;
325 if( found_intersection ) {
326 End_PointAndTangent= IntersectPointVelct_G;
327 StepTaken = TruePathLength = IntersectPointVelct_G.
GetCurveLength()
329 } else {
330
331 if( recalculatedEndPt ){
332 CurrentState= IntersectPointVelct_G;
333 }
334 }
335 }
336 if( !intersects )
337 {
338 StepTaken += s_length_taken;
339
340 if (fpTrajectoryFilter) {
342 }
343 }
344 first_substep = false;
345
346#ifdef G4DEBUG_FIELD
347 if( fNoZeroStep > fActionThreshold_NoZeroSteps ) {
349 CurrentState, CurrentProposedStepLength,
350 NewSafety, do_loop_count, pPhysVol );
351 }
352 if( (fVerboseLevel > 1) && (do_loop_count > fMax_loop_count-10 )) {
353 if( do_loop_count == fMax_loop_count-9 ){
354 G4cout <<
" G4PropagatorInField::ComputeStep(): " <<
G4endl
355 <<
" Difficult track - taking many sub steps." <<
G4endl;
356 }
357 printStatus( SubStepStartState, CurrentState, CurrentProposedStepLength,
358 NewSafety, do_loop_count, pPhysVol );
359 }
360#endif
361
362 do_loop_count++;
363
364 } while( (!intersects )
365 && (StepTaken + kCarTolerance < CurrentProposedStepLength)
366 && ( do_loop_count < fMax_loop_count ) );
367
368 if( do_loop_count >= fMax_loop_count )
369 {
370 fParticleIsLooping = true;
371
372 if ( fVerboseLevel > 0 )
373 {
374 G4cout <<
" G4PropagateInField::ComputeStep(): " <<
G4endl
375 << " Killing looping particle "
376
377 << " after " << do_loop_count << " field substeps "
378 << " totaling " << StepTaken / mm << " mm " ;
379 if( pPhysVol )
381 else
382 G4cout <<
" in unknown or null volume. " ;
384 }
385 }
386
387 if( !intersects )
388 {
389
390
391
392 End_PointAndTangent = CurrentState;
393 TruePathLength = StepTaken;
394 }
395
396
397
398 pFieldTrack = End_PointAndTangent;
399
400#ifdef G4VERBOSE
401
402
404 - End_PointAndTangent.
GetCurveLength()) > 3.e-4 * TruePathLength )
405 {
406 std::ostringstream message;
407 message << "Curve length mis-match between original state "
408 <<
"and proposed endpoint of propagation." <<
G4endl
409 << " The curve length of the endpoint should be: "
411 << " and it is instead: "
413 << " A difference of: "
416 <<
" Original state = " << OriginalState <<
G4endl
417 << " Proposed state = " << End_PointAndTangent;
420 }
421#endif
422
423
424
425
426
427 if( ( (TruePathLength < fZeroStepThreshold)
428 && ( TruePathLength+kCarTolerance < CurrentProposedStepLength )
429 )
430 || ( TruePathLength < 0.5*kCarTolerance )
431 )
432 {
433 fNoZeroStep++;
434 }
435 else{
436 fNoZeroStep = 0;
437 }
438
439 if( fNoZeroStep > fAbandonThreshold_NoZeroSteps )
440 {
441 fParticleIsLooping = true;
442 std::ostringstream message;
443 message <<
"Particle is stuck; it will be killed." <<
G4endl
444 << " Zero progress for " << fNoZeroStep << " attempted steps."
446 << " Proposed Step is " << CurrentProposedStepLength
447 <<
" but Step Taken is "<< fFull_CurveLen_of_LastAttempt <<
G4endl
448 << " For Particle with Charge = " << fCharge
449 << " Momentum = "<< fInitialMomentumModulus
450 <<
" Mass = " << fMass <<
G4endl;
451 if( pPhysVol )
452 message <<
" in volume " << pPhysVol->
GetName() ;
453 else
454 message << " in unknown or null volume. " ;
457 fNoZeroStep = 0;
458 }
459
460 return TruePathLength;
461}
void SetChargeMomentumMass(G4double pCharge, G4double pMomentum, G4double pMass)
G4double AdvanceChordLimited(G4FieldTrack &yCurrent, G4double stepInitial, G4double epsStep_Relative, const G4ThreeVector latestSafetyOrigin, G4double lasestSafetyRadius)
G4double GetMinimumEpsilonStep() const
G4double GetDeltaOneStep() const
const G4ThreeVector & GetMomentumDir() const
G4double GetCurveLength() const
G4ThreeVector GetPosition() const
virtual void LocateGlobalPointWithinVolume(const G4ThreeVector &position)
G4VPhysicalVolume * GetWorldVolume() const
void printStatus(const G4FieldTrack &startFT, const G4FieldTrack ¤tFT, G4double requestStep, G4double safety, G4int step, G4VPhysicalVolume *startVolume)
G4FieldManager * FindAndSetFieldManager(G4VPhysicalVolume *pCurrentPhysVol)
void PrintStepLengthDiagnostic(G4double currentProposedStepLength, G4double decreaseFactor, G4double stepTrial, const G4FieldTrack &aFieldTrack)
G4ChordFinder * GetChordFinder()
G4bool IntersectChord(const G4ThreeVector &StartPointA, const G4ThreeVector &EndPointB, G4double &NewSafety, G4double &LinearStepLength, G4ThreeVector &IntersectionPoint)
void SetEpsilonStep(G4double newEps)
void CreateNewTrajectorySegment()
virtual void TakeIntermediatePoint(G4ThreeVector newPoint)=0
G4LogicalVolume * GetLogicalVolume() const
const G4String & GetName() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)