Good day,
I am simulating UV and visible photons with 10.07. I noticed that even when some photons are never leaving the inner material (total internal reflection), I would see some steps with Velocities (=stepLength/stepTime) corresponding to the outer material.
After pocking around the code, I think I found out why (and I think it is valid for all reflections):
deltaT is calculated in G4Transportation::AlongStepDoIt, using stepData.GetPreStepPoint()->GetVelocity(). At this point, the PreStepPoint has received information from the previous step PostStepPoint (Copied by G4SteppingManager).
Now, what do we have in PostStepPoint->GetVelocity() at the end of the previous step?
A reflection is handled by G4OpBoundaryProcess in a 2-step fashion.
- Step 1: Particle starts in “inner volume”. ProposeVelocity in G4OpBoundaryProcess gives the velocity of the inner volume (G4Track::CalculateVelocityForOpticalPhoton). The momentum and polarization are changed. Despite the reflection, the NextVolume is “outer volume” as G4OpBoundaryProcess cannot change the output of the transportation.
- Step 2: The Particle is in “outer volume”, pointing towards towards “inner volume”. Since the step length is 0, G4OpBoundaryProcess exits without doing much, except call CalculateVelocityForOpticalPhoton, which this time gives the velocity of the “outer volume”. Because Propose was called isVelocityChanged is set to true, we don’t recalculate (but it would not change the outcome). pPostStepPoint receives this velocity.
In the step coming after, although we are in the inner volume, PreStepPoint() has the velocity of the outer volume and the deltaT is calculated incorrectly.
To make sure the issue was not with my code, I show below the output from examples/extended/optical/OpNovice. I modified slightly SteppingVerbose to show the time and velocity used for each step, as well as the volume the photon is in.
G4WT0 > *********************************************************************************************************
G4WT0 > * G4Track Information: Particle = opticalphoton, Track ID = 21, Parent ID = 1
G4WT0 > *********************************************************************************************************
G4WT0 >
G4WT0 > ----------- CalculateVelocityForOpticalPhoton Material: Water, momentum: 2.372646615000221e-06, result: 218.1919701955274
G4WT0 > Step# X(mm) Y(mm) Z(mm) KinE(MeV) dE(MeV) StepLeng TrackLeng GlobalTime NextVolume ProcName
G4WT0 > 0 0.606 0.19 -0.113 2.37e-06 0 0 0 Tank initStep
G4WT0 > G4VDiscreteProcess::PostStepGetPhysicalInteractionLength() - [ OpBoundary]
G4WT0 > Particle type - opticalphoton
G4WT0 > mass: 0[GeV]
G4WT0 > charge: 0[e]
G4WT0 > Direction x: 0.1133339441716328, y: -0.7761738403496036, z: 0.6202496164089497
G4WT0 > Total Momentum = 2.372646615000221e-09[GeV]
G4WT0 > Momentum: 2.689013990034485e-10[GeV], y: -1.841586234957209e-09[GeV], z: 1.47163315282788e-09[GeV]
G4WT0 > Total Energy = 2.372646615000221e-09[GeV]
G4WT0 > Kinetic Energy = 2.372646615000221e-09[GeV]
G4WT0 > MagneticMoment [MeV/T]: 0
G4WT0 > ProperTime = 0[ns]
G4WT0 > in Material Water
G4WT0 > InteractionLength= 1.797693134862316e+307[cm]
G4WT0 > Photon at Boundary!
G4WT0 > thePrePV: Tank
G4WT0 > thePostPV: World
G4WT0 > Old Momentum Direction: (0.1133339441716328,-0.7761738403496036,0.6202496164089497)
G4WT0 > Old Polarization: (-0.09318586388267827,0.6132161051617853,0.784399390070294)
G4WT0 > New Momentum Direction: (0.1423197194964792,0.8598513013574984,0.4903068804292471)
G4WT0 > New Polarization: (-0.09318586388267827,0.6132161051617853,0.784399390070294)
G4WT0 > *** LobeReflection ***
G4WT0 > #Step# X(mm) Y(mm) Z(mm) KinE(MeV) dE(MeV) StepLeng StepTime Velocity TrackLeng Volume NextVolume ProcName
G4WT0 > 1 731 -5e+03 4e+03 2.37e-06 0 6.44e+03 29.5 218 6.44e+03 Tank World Transportation
G4WT0 > G4VDiscreteProcess::PostStepGetPhysicalInteractionLength() - [ OpBoundary]
G4WT0 > Particle type - opticalphoton
G4WT0 > mass: 0[GeV]
G4WT0 > charge: 0[e]
G4WT0 > Direction x: 0.1423197194964792, y: 0.8598513013574984, z: 0.4903068804292471
G4WT0 > Total Momentum = 2.372646615000221e-09[GeV]
G4WT0 > Momentum: 3.376744007111023e-10[GeV], y: 2.040123279569404e-09[GeV], z: 1.163324960161771e-09[GeV]
G4WT0 > Total Energy = 2.372646615000221e-09[GeV]
G4WT0 > Kinetic Energy = 2.372646615000221e-09[GeV]
G4WT0 > MagneticMoment [MeV/T]: 0
G4WT0 > ProperTime = 0[ns]
G4WT0 > in Material Air
G4WT0 > InteractionLength= 1.797693134862316e+307[cm]
G4WT0 > Photon at Boundary!
G4WT0 > thePrePV: World
G4WT0 > thePostPV: Tank
G4WT0 > *** StepTooSmall ***
G4WT0 > #Step# X(mm) Y(mm) Z(mm) KinE(MeV) dE(MeV) StepLeng StepTime Velocity TrackLeng Volume NextVolume ProcName
G4WT0 > 2 731 -5e+03 4e+03 2.37e-06 0 0 0 N/A 6.44e+03 World Tank Transportation
G4WT0 > G4VDiscreteProcess::PostStepGetPhysicalInteractionLength() - [ OpBoundary]
G4WT0 > Particle type - opticalphoton
G4WT0 > mass: 0[GeV]
G4WT0 > charge: 0[e]
G4WT0 > Direction x: 0.1423197194964792, y: 0.8598513013574984, z: 0.4903068804292471
G4WT0 > Total Momentum = 2.372646615000221e-09[GeV]
G4WT0 > Momentum: 3.376744007111023e-10[GeV], y: 2.040123279569404e-09[GeV], z: 1.163324960161771e-09[GeV]
G4WT0 > Total Energy = 2.372646615000221e-09[GeV]
G4WT0 > Kinetic Energy = 2.372646615000221e-09[GeV]
G4WT0 > MagneticMoment [MeV/T]: 0
G4WT0 > ProperTime = 0[ns]
G4WT0 > in Material Water
G4WT0 > InteractionLength= 1.797693134862316e+307[cm]
G4WT0 > Photon at Boundary!
G4WT0 > thePrePV: Tank
G4WT0 > thePostPV: World
G4WT0 > Old Momentum Direction: (0.1423197194964792,0.8598513013574984,0.4903068804292471)
G4WT0 > Old Polarization: (-0.09318586388267827,0.6132161051617853,0.784399390070294)
G4WT0 > New Momentum Direction: (0.1821400361406186,0.8135098843566156,-0.5522922010030378)
G4WT0 > New Polarization: (-0.09318586388267827,0.6132161051617853,0.784399390070294)
G4WT0 > *** LobeReflection ***
G4WT0 > Updating track with new velocity: 218.1919701955274
G4WT0 > #Step# X(mm) Y(mm) Z(mm) KinE(MeV) dE(MeV) StepLeng StepTime Velocity TrackLeng Volume NextVolume ProcName
G4WT0 > 3 1.02e+03 -3.24e+03 5e+03 2.37e-06 0 2.05e+03 6.83 300 8.49e+03 Tank World Transportation
G4WT0 > G4VDiscreteProcess::PostStepGetPhysicalInteractionLength() - [ OpBoundary]
G4WT0 > Particle type - opticalphoton
G4WT0 > mass: 0[GeV]
G4WT0 > charge: 0[e]
G4WT0 > Direction x: 0.1821400361406186, y: 0.8135098843566156, z: -0.5522922010030378
G4WT0 > Total Momentum = 2.372646615000221e-09[GeV]
G4WT0 > Momentum: 4.321539402050567e-10[GeV], y: 1.930171473387945e-09[GeV], z: -1.310394221200879e-09[GeV]
G4WT0 > Total Energy = 2.372646615000221e-09[GeV]
G4WT0 > Kinetic Energy = 2.372646615000221e-09[GeV]
G4WT0 > MagneticMoment [MeV/T]: 0
G4WT0 > ProperTime = 0[ns]
G4WT0 > in Material Air
G4WT0 > InteractionLength= 1.797693134862316e+307[cm]
G4WT0 > Photon at Boundary!
G4WT0 > thePrePV: World
G4WT0 > thePostPV: Tank
G4WT0 > *** StepTooSmall ***
G4WT0 > Updating track with new velocity: 299.792458
G4WT0 > #Step# X(mm) Y(mm) Z(mm) KinE(MeV) dE(MeV) StepLeng StepTime Velocity TrackLeng Volume NextVolume ProcName
G4WT0 > 4 1.02e+03 -3.24e+03 5e+03 2.37e-06 0 0 0 N/A 8.49e+03 World Tank Transportation
G4WT0 > G4VDiscreteProcess::PostStepGetPhysicalInteractionLength() - [ OpBoundary]
G4WT0 > Particle type - opticalphoton
G4WT0 > mass: 0[GeV]
G4WT0 > charge: 0[e]
G4WT0 > Direction x: 0.1821400361406186, y: 0.8135098843566156, z: -0.5522922010030378
G4WT0 > Total Momentum = 2.372646615000221e-09[GeV]
G4WT0 > Momentum: 4.321539402050567e-10[GeV], y: 1.930171473387945e-09[GeV], z: -1.310394221200879e-09[GeV]
G4WT0 > Total Energy = 2.372646615000221e-09[GeV]
G4WT0 > Kinetic Energy = 2.372646615000221e-09[GeV]
G4WT0 > MagneticMoment [MeV/T]: 0
G4WT0 > ProperTime = 0[ns]
G4WT0 > in Material Water
G4WT0 > InteractionLength= 1.797693134862316e+307[cm]
G4WT0 > #Step# X(mm) Y(mm) Z(mm) KinE(MeV) dE(MeV) StepLeng StepTime Velocity TrackLeng Volume NextVolume ProcName
G4WT0 > 5 2.06e+03 1.4e+03 1.85e+03 2.37e-06 2.37e-06 5.7e+03 19 300 1.42e+04 Tank Tank OpAbsorption
I don’t quite know how to fix this issue. In G4 9.6, the problem was not apparent because G4Track::GetVelocity() was not re-calculating the Velocity. During the 0-length step, the velocity was from the inner volume, Proposed by G4OpBoundaryProcess, which was then propagated to the step after… It seems that setting G4Track.useGivenVelocity to “true” would bring back this behavior, but I don’t see (yet) how to implement into my code.
It seems to me that 2 trains of thought collided here: the “Propose” school that would re-calculate the Velocity only when needed, and call “CalculateVelocityForOpticalPhoton” at each step. The problem is that we do that after the deltaT has been computed, using the volume we are IN and not the next volume…
One possibility I see would be to forget all about the “ProposeVelocity” and call track.GetVelocity() instead of stepData.GetPreStepPoint()->GetVelocity() in G4Transportation::AlongStepDoIt.
Let me know what your thoughts are,
Thanks.