anchor.diff:  #ifdef PYTHIA
   10  continue
       if(pythia_useRealBeam) then
  -
  -      pIn(1) = yv1(j)*ejfv(j)*c1m6
  -      pIn(2) = yv2(j)*ejfv(j)*c1m6
  -      pIn(3) = sqrt(ejfv(j)**2 - pIn(1)**2 - pIn(2)**2)*c1m3
  +      ! We sample the actual SixTrack tracked beam 1, and sample a colliding
  +      ! particle from beam 2
  +
  +      ! Compute the particle three-momenta in GeV, using division rather
  +      ! than multiplication for scaling to avoid additional numerical error
  +      pIn(1) = yv1(j)*ejfv(j)/c1e6
  +      pIn(2) = yv2(j)*ejfv(j)/c1e6
  +      pIn(3) = sqrt(ejfv(j)**2 - pIn(1)**2 - pIn(2)**2)/c1e3
   
         call scatter_generateParticle(idPro, iElem, j, pGen)
   
  
positive.code:       orbXP = scatter_proList(idPro)%fParams(1)
        orbYP = scatter_proList(idPro)%fParams(2)
      end if
      pVec(1) = (orbXP*e0f)*c1m3
      pVec(2) = (orbYP*e0f)*c1m3
      pVec(3) = -sqrt(e0f**2 - pVec(1)**2 - pVec(2)**2)
  
    case(scatter_proBeamUnCorr)
  
  
negative.code:     sigmaX  = scatter_proList(idPro)%fParams(12)
      sigmaY  = scatter_proList(idPro)%fParams(13)
      normFac = scatter_proList(idPro)%fParams(14)
      ! sRot    = scatter_proList(idPro)%fParams(15)
      ! cRot    = scatter_proList(idPro)%fParams(16)
      ! sOrbXP  = scatter_proList(idPro)%fParams(18)
      ! sOrbYP  = scatter_proList(idPro)%fParams(19)
  
      ! xRot    = (xPos - orbX + sPos*sOrbXP)*cRot - (yPos - orbY + sPos*sOrbXP)*sRot
      ! yRot    = (xPos - orbX + sPos*sOrbYP)*sRot + (yPos - orbY + sPos*sOrbYP)*cRot
  