Problem 1990 - Rotation in G4CascadeInterface::createBullet() leads to irreproducible results
Summary: Rotation in G4CascadeInterface::createBullet() leads to irreproducible results
Status: RESOLVED FIXED
Alias: None
Product: Geant4
Classification: Unclassified
Component: processes/hadronic/models/cascade (show other problems)
Version: 10.1
Hardware: PC Linux
: P4 normal
Assignee: dennis.herbert.wright
URL:
Depends on:
Blocks:
 
Reported: 2017-07-05 11:18 CEST by Kirill Chilikin
Modified: 2017-11-27 21:20 CET (History)
0 users

See Also:


Attachments

Note You need to log in before you can comment on or make changes to this problem.
Description Kirill Chilikin 2017-07-05 11:18:37 CEST
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.
Comment 1 dennis.herbert.wright 2017-11-27 21:20:43 CET
Thanks for reporting this.  The rotation has been removed.