The following is observed for Belle II simulation with Geant version geant4-10-06-patch-01 When performing backward extrapolation of a track, Geant calls the function G4InterpolationDriver<T>::AdvanceChordLimited(). At the first step, this function calls Base::GetEquationOfMotion()->RightHandSide(y, fdydx); Which calls G4ErrorMag_UsualEqRhs::EvaluateRhsGivenB(). The equation of motion G4ErrorMag_UsualEqRhs reverses the direction of magnetic field: B-> -B for G4ErrorMode_PropBackwards. Thus, the values fdydx are calculated using the reversed field. The stack trace at this point is: #0 G4ErrorMag_UsualEqRhs::EvaluateRhsGivenB (this=0x2d2cd1e0, y=0x7fffffff3150, B=0x7fffffff3040, dydx=0x2ba2eb18) at source/geometry/magneticfield/src/G4ErrorMag_UsualEqRhs.cc:53 #1 0x00007fffe484db3d in G4EquationOfMotion::RightHandSide (this=0x2d2cd1e0, y=0x7fffffff3150, dydx=0x2ba2eb18) at source/geometry/magneticfield/include/G4EquationOfMotion.icc:71 #2 0x00007fffe485bb92 in G4InterpolationDriver<G4DormandPrince745>::AdvanceChordLimited (this=0x2ba2ea80, track=..., hstep=100, epsStep=0.0001, chordDistance=0.25) at source/geometry/magneticfield/include/G4InterpolationDriver.icc:218 #3 0x00007fffe484c3f4 in G4BFieldIntegrationDriver::AdvanceChordLimited ( this=0x2c832770, yCurrent=..., stepMax=100, epsStep=0.0001, chordDistance=0.25) at source/geometry/magneticfield/src/G4BFieldIntegrationDriver.cc:104 #4 0x00007fffe48f2773 in G4ChordFinder::AdvanceChordLimited (this=0x2c831af0, yCurrent=..., stepInitial=100, epsStep_Relative=0.0001) at source/geometry/magneticfield/include/G4ChordFinder.icc:95 #5 0x00007fffe48f040a in G4PropagatorInField::ComputeStep (this=0x1a300db0, pFieldTrack=..., CurrentProposedStepLength=100, currentSafety=@0x7fffffff3de0: 0, pPhysVol=0x16586cc0, canRelaxDeltaChord=false) at source/geometry/navigation/src/G4PropagatorInField.cc:318 #6 0x00007fffdc2a2ee0 in G4Transportation::AlongStepGetPhysicalInteractionLength (this=0x2b63a5a0, track=..., currentMinimumStep=100, currentSafety=@0x7fffffff3de0: 0, selection=0x1a2f965c) at source/processes/transportation/src/G4Transportation.cc:358 #7 0x00007fffe40a05ca in G4VProcess::AlongStepGPIL (this=0x2b63a5a0, track=..., previousStepSize=0, currentMinimumStep=100, proposedSafety=@0x7fffffff3de0: 0, selection=0x1a2f965c) at source/processes/management/include/G4VProcess.hh:488 #8 0x00007fffe409ed43 in G4SteppingManager::DefinePhysicalStepLength ( this=0x1a2f94d0) at source/tracking/src/G4SteppingManager2.cc:242 #9 0x00007fffe409bbf0 in G4SteppingManager::Stepping (this=0x1a2f94d0) at source/tracking/src/G4SteppingManager.cc:179 #10 0x00007fffd978516c in G4ErrorPropagator::MakeOneStep (this=0x2ee39af0, currentTS_FREE=0x7fffffff4be0) at source/error_propagation/src/G4ErrorPropagator.cc:375 #11 0x00007fffd9784b9d in G4ErrorPropagator::PropagateOneStep ( this=0x2ee39af0, currentTS=0x7fffffff4be0) at source/error_propagation/src/G4ErrorPropagator.cc:243 Then, the step is calculated by OneGoodStep(). This function calls G4DormandPrince745::Stepper(), which call RightHandSide(). This time, however, the equation of motion G4Mag_UsualEqRhs is used. This equation does not have any special condition to reverse the magnetic field direction. As the result, the magnetic field used to calculate the step is not consistent with the field used previously. The part of stack trace that differs from the previous one is: #0 G4Mag_UsualEqRhs::EvaluateRhsGivenB (this=0x2c831b50, y=0x7fffffff2d90, B=0x7fffffff2c60, dydx=0x2c832098) at source/geometry/magneticfield/src/G4Mag_UsualEqRhs.cc:50 #1 0x00007fffe484db3d in G4EquationOfMotion::RightHandSide (this=0x2c831b50, y=0x7fffffff2d90, dydx=0x2c832098) at source/geometry/magneticfield/include/G4EquationOfMotion.icc:71 #2 0x00007fffe484dbc7 in G4MagIntegratorStepper::RightHandSide ( this=0x2c832070, y=0x7fffffff2d90, dydx=0x2c832098) at source/geometry/magneticfield/include/G4MagIntegratorStepper.icc:105 #3 0x00007fffe4865477 in G4DormandPrince745::Stepper (this=0x2c832070, yInput=0x7fffffff3150, dydx=0x2ba2eb18, hstep=100, yOut=0x7fffffff3000, yErr=0x7fffffff3060) at source/geometry/magneticfield/src/G4DormandPrince745.cc:144 #4 0x00007fffe4865162 in G4DormandPrince745::Stepper (this=0x2c832070, yInput=0x7fffffff3150, dydx=0x2ba2eb18, hstep=100, yOutput=0x7fffffff3000, yError=0x7fffffff3060, dydxOutput=0x7fffffff2fa0) at source/geometry/magneticfield/src/G4DormandPrince745.cc:73 #5 0x00007fffe485e3e3 in G4InterpolationDriver<G4DormandPrince745>::OneGoodStep (this=0x2ba2ea80, it=..., y=..., dydx=..., hstep=@0x7fffffff3148: 100, epsStep=0.0001, curveLength=0) at source/geometry/magneticfield/include/G4InterpolationDriver.icc:455 #6 0x00007fffe485be01 in G4InterpolationDriver<G4DormandPrince745>::AdvanceChordLimited (this=0x2ba2ea80, track=..., hstep=100, epsStep=0.0001, chordDistance=0.25) at source/geometry/magneticfield/include/G4InterpolationDriver.icc:251
Hello John, Kirill, in CMS for backward propagation we historically are using G4CrassicalRK4. Will it have the same problem? Vladimir
Indeed what you report is a problem. Could you please clarify how the field propagation was set up - i.e. how the equation of motion, stepper and driver were created ?
Sorry for late reply. I am not sure if the information I am providing below is sufficient; I am not usually working on this, I just found the bug. Please ask for additional information if necessary. First, G4ErrorFreeTrajState is created using the measured data: G4ErrorFreeTrajState g4eState(nameG4e, positionG4e, momentumG4e, covarianceG4e); Then, the extrapolation is initialized with the following function: if (m_G4RunMgr) { m_G4RunMgr->SetUserAction((G4UserTrackingAction*)NULL); m_G4RunMgr->SetUserAction((G4UserSteppingAction*)NULL); } if (mode == G4ErrorMode_PropBackwards) { if (m_StdStepper) { m_StdStepper->SetEquationOfMotion(m_BackwardEquationOfMotion); } } if (m_Propagator == NULL) m_Propagator = new G4ErrorPropagator(); m_Propagator->SetStepN(0); G4ErrorPropagatorData::GetErrorPropagatorData()->SetState(G4ErrorState_Propagating); where equations were previously initialized to m_StdStepper = const_cast<G4MagIntegratorStepper*>(fieldManager->GetChordFinder()->GetIntegrationDriver()->GetStepper()); m_ForwardEquationOfMotion = m_StdStepper->GetEquationOfMotion(); m_BackwardEquationOfMotion = new G4ErrorMag_UsualEqRhs(m_MagneticField); We can see that G4ErrorMag_UsualEqRhs is set explicitly; perhaps it is not a correct thing to do that anymore? Then, it may be a Belle II rather than Geant bug.
Thanks for providing the key part of your user code. In Geant4 10.6 we changed the default integration method to a new driver 'G4BFieldIntegrationDriver' which uses the 'G4InterpolationDriver<G4DormandPrince745>' integration driver (and another method for long steps, whose angle is beyond 2 pi). This integration method creates multiple instances of the 'stepper' in order to maintain a set of smaller integration intervals with which to interpolate. As a result what you obtain in the call m_StdStepper = const_cast<G4MagIntegratorStepper*>(fieldManager->GetChordFinder()->GetIntegrationDriver()->GetStepper()); is only one of these steppers. To confirm that this is the problem, you can add the following method if you have a local installation of Geant4. In G4InterpolationDriver.hh add the method: public: void SetEquationOfMotion(G4EquationOfMotion* equation) override; and in G4InterpolationDriver.icc its implementation: template <class T> void G4InterpolationDriver<T>::SetEquationOfMotion(G4EquationOfMotion* equation) { // for (StepperIterator is = fSteppers.begin(); is <= fLastStepper; ++is) // { // is->stepper->SetEquationOfMotion(equation); // } G4Exception("G4InterpolationDriver: SetEquationOfMotion is incompatible with this Driver ", "GeomField0003", FatalException, "please use a different type of Driver."); } It will be more reliable to use the most basic way to construct the equation of motion, and pass it to a stepper (and driver) which you construct. In this case you will be in control of the type of driver. Using the existing / old method m_StdStepper->SetEquationOfMotion(m_BackwardEquationOfMotion); could only be fully relied upon in this case. The chain to control your creation of equation stepper, driver can be created in the following order: 1. Equation: G4Mag_UsualEqRhs or G4ErrorMag_UsualRhs 2. Stepper: G4ClassicalRK4 or G4DormandPrince745 auto stepper = new G4DormandPrince745( fEquation ); 3. Driver: G4IntegrationDriver<stepper_type> or the old G4Mag_IntDriver G4VIntegrationDriver driver = new G4IntegrationDriver<G4DormandPrince745>( stepMinimum, stepper, 6 ); 4. G4ChordFinder G4ChordFinder chordFinder = new G4ChordFinder( driver ); Once you have ensured that the driver is of your type, then the call to change the equation will be reliable without the proposed change above. ( Else one could also control the equation from the G4FieldManager's ConfigureForTrack method. This could simply switch between a 'forward driver' and a 'backward driver' using the flag comparison G4ErrorPropagatorData::GetErrorPropagatorData()->GetMode() == G4ErrorMode_PropBackwards) which is currently made in G4ErrorMag_UsualRhs. )
The recipe should work. Can you please let us know if you have tried it, and if it worked?
