Fix for bug #76597: floating point exception in the test of Pythia 6.4 with Omega...

diff --git a/PYTHIA6/pythia6.4.21/pythia-6.4.21.f b/PYTHIA6/pythia6.4.21/pythia-6.4.21.f
index b0c1bb9610babebbfaf4a048128d77b2cb1ae83d..521a846659eed98cec1ee493c1b3dd3d6bf5a4cf 100644 (file)
--- a/PYTHIA6/pythia6.4.21/pythia-6.4.21.f
+++ b/PYTHIA6/pythia6.4.21/pythia-6.4.21.f
@@ -63836,7 +63836,11 @@ C...(Borrow slot N+2 for temporary direction.)
       DO 850 I=IC1+1,IC2-1
         IF((K(I,1).EQ.1.OR.K(I,1).EQ.2).AND.
      &  KCHG(PYCOMP(K(I,2)),2).NE.0) THEN
-          FRAC1=FOUR(IC2,I)/(FOUR(IC1,I)+FOUR(IC2,I))
+          IF (ABS(FOUR(IC1,I)+FOUR(IC2,I)).GT.0.D0) THEN
+             FRAC1=FOUR(IC2,I)/(FOUR(IC1,I)+FOUR(IC2,I))
+          ELSE
+             FRAC1 = 1.D0
+          ENDIF
           DO 840 J=1,4
             P(N+2,J)=P(N+2,J)+FRAC1*P(I,J)
   840     CONTINUE
@@ -64042,8 +64046,13 @@ C...Calculate Pg, a part of which will be added to Phad later. (EN)
           ALPHA=1D0
           BETA=1D0
         ELSE
-          ALPHA=FOUR(NSAV+1,I2)/FOUR(I1,I2)
-          BETA=FOUR(NSAV+1,I1)/FOUR(I1,I2)
+           IF (ABS(FOUR(I1,I2)).GT.0.D0) THEN
+              ALPHA=FOUR(NSAV+1,I2)/FOUR(I1,I2)
+              BETA=FOUR(NSAV+1,I1)/FOUR(I1,I2)
+           ELSE
+              ALPHA=1D0
+              BETA=1D0
+           ENDIF
         ENDIF
         DO 1030 J=1,4
           PG(J)=ALPHA*P(I1,J)+BETA*P(I2,J)