Problem 1990

Summary: Rotation in G4CascadeInterface::createBullet() leads to irreproducible results
Product: Geant4 Reporter: Kirill Chilikin <chilikin>
Component: processes/hadronic/models/cascadeAssignee: dennis.herbert.wright
Status: RESOLVED FIXED    
Severity: normal    
Priority: P4    
Version: 10.1   
Hardware: PC   
OS: Linux   

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.