| Summary: | Rotation in G4CascadeInterface::createBullet() leads to irreproducible results | ||
|---|---|---|---|
| Product: | Geant4 | Reporter: | Kirill Chilikin <chilikin> |
| Component: | processes/hadronic/models/cascade | Assignee: | dennis.herbert.wright |
| Status: | RESOLVED FIXED | ||
| Severity: | normal | ||
| Priority: | P4 | ||
| Version: | 10.1 | ||
| Hardware: | PC | ||
| OS: | Linux | ||
Thanks for reporting this. The rotation has been removed. |
In the function G4CascadeInterface::createBullet(), Geant4 tries to determine a rotation to an original frame: G4LorentzVector projectileMomentum = aTrack.Get4Momentum()/GeV; // Rotation/boost to get from z-axis back to original frame bulletInLabFrame = G4LorentzRotation::IDENTITY; // Initialize bulletInLabFrame.rotateZ(-projectileMomentum.phi()); bulletInLabFrame.rotateY(-projectileMomentum.theta()); bulletInLabFrame.invert(); However, the momentum aTrack.Get4Momentum() is already rotated in such a way that its direction is along the z axis. This rotation is performed by G4HadProjectile::Initialise(). This results in the following: A comparison of two simulations of a specific Belle II event is presented below. The random seeds are the same, the difference is that the system is CentOS 6 in the first case and CentOS 7 in the second case. In this event, a hadronic interaction of a proton and an iron nucleus happens, this interaction is handled by the function G4CascadeInterface::ApplyYourself(). The values of projectileMomentum are: CentOS 7: (gdb) p projectileMomentum $196 = {pp = {dx = -8.8817841970012525e-19, dy = 4.4408920985006263e-19, dz = 0.0091146499630734134, static tolerance = 2.22045e-14}, ee = 0.93960956921259553, static tolerance = 2.22045e-14, static metric = 1} CentOS 6: (gdb) p projectileMomentum $156 = {pp = {dx = -1.3322676295501879e-18, dy = 4.4408920985006263e-19, dz = 0.0091146499630734117, static tolerance = 2.22045e-14}, ee = 0.93960956921259553, static tolerance = 2.22045e-14, static metric = 1} The momentum projectileMomentum is already rotated and its direction is along the z axis, just as expected. However, the dx and dz values are not exactly the same, there is a difference of the order of 10^{-18}. This is a numerical-noise level change and might be caused, for example, by minor changes in the system libraries. For the dz component, this difference is many orders of magnitude lower than the value, and does not result in any problems. For the dx component, however, the difference is of the same order as the value itself. This results in a large difference of the values of the rotation angle phi: CentOS 7: (gdb) p -projectileMomentum.phi() $195 = -2.677945044588987 CentOS 6: (gdb) p -projectileMomentum.phi() $155 = -2.819842099193151 After bulletInLabFrame.rotateZ(-projectileMomentum.phi()), the difference in phi propagates into the rotation: CentOS 7: (gdb) n 464 bulletInLabFrame.rotateY(-projectileMomentum.theta()); (gdb) p bulletInLabFrame $199 = {static IDENTITY = { static IDENTITY = <same as static member of an already seen type>, mxx = 1, mxy = 0, mxz = 0, mxt = 0, myx = 0, myy = 1, myz = 0, myt = 0, mzx = 0, mzy = 0, mzz = 1, mzt = 0, mtx = 0, mty = 0, mtz = 0, mtt = 1}, mxx = -0.89442719099991586, mxy = 0.44721359549995809, mxz = 0, mxt = 0, myx = -0.44721359549995809, myy = -0.89442719099991586, myz = -0, myt = -0, mzx = 0, mzy = 0, mzz = 1, mzt = 0, mtx = 0, mty = 0, mtz = 0, mtt = 1} CentOS 6: (gdb) n 464 bulletInLabFrame.rotateY(-projectileMomentum.theta()); (gdb) p bulletInLabFrame $158 = {static IDENTITY = { static IDENTITY = <same as static member of an already seen type>, mxx = 1, mxy = 0, mxz = 0, mxt = 0, myx = 0, myy = 1, myz = 0, myt = 0, mzx = 0, mzy = 0, mzz = 1, mzt = 0, mtx = 0, mty = 0, mtz = 0, mtt = 1}, mxx = -0.94868329805051377, mxy = 0.316227766016838, mxz = 0, mxt = 0, myx = -0.316227766016838, myy = -0.94868329805051377, myz = -0, myt = -0, mzx = 0, mzy = 0, mzz = 1, mzt = 0, mtx = 0, mty = 0, mtz = 0, mtt = 1} Finally, this results in secondary particles of this interaction having different momentum direction.