GEANT321/fiface/flufin.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.3  1997/06/20 18:32:55  japost
   6 * A summary of the problem buf fix and its original report:
   7 *
   8 * ------------------------------------------------------------------------------
   9 * From: Shawn McKee - (313) 764-4395 <mckee@pooh.physics.lsa.umich.edu>
  10 * Subject: RE: Problem with FLUKA interface
  11 * Date: Thu, 12 Jun 97 18:09:31 METDST
  12 *
  13 *
  14 * I have seen this error with FLUKA.  I traced it to
  15 * the following statement:
  16 *
  17 *      IF (RNDEVT .GE. SINE/FSIG) THEN
  18 *
  19 * Apparently FSIG (for this event) was not initialized
  20 * properly (it is 0.0).  My fix was to insert a check
  21 * on FSIG in the FLUFIN routine source.
  22 *
  23 * ------------------------------------------------------------------------------
  24 * From: pniemine@estwm0.wm.estec.esa.nl
  25 * Subject: Problem with FLUKA interface - cern.heplib #5667
  26 * Date: Thu, 12 Jun 97 17:21:14 METDST
  27 *
  28 * In trying to run a GEANT program with FLUKA interface in our VMS-AXP (Alpha)
  29 * cluster, I get the following error message (originating from FLUFIN) once
  30 * a hadronic interaction is to be simulated:
  31 *
  32 * arithmetic trap, floating/decimal divide by zero
  33 * ------------------------------------------------------------------------------
  34 *
  35 * Revision 1.2  1996/02/27 14:12:57  ravndal
  36 * Correction for overstopped antiprotons
  37 *
  38 * Revision 1.1.1.1  1995/10/24 10:19:53  cernlib
  39 * Geant
  40 *
  41 *
  42 #include "geant321/pilot.h"
  43 *CMZ :  3.21/02 19/05/94  13.35.12  by  S.Ravndal
  44 *FCA :          05/01/99  09:47:25  by  Federico Carminati
  45 *               Added Xi- to the list of particles that are forced
  46 *               to decay
  47 *-- Author :
  48       SUBROUTINE FLUFIN
  49 #include "geant321/gcbank.inc"
  50 #include "geant321/gccuts.inc"
  51 #include "geant321/gcjloc.inc"
  52 #include "geant321/gcflag.inc"
  53 #include "geant321/gckine.inc"
  54 #include "geant321/gcking.inc"
  55 #include "geant321/gcmate.inc"
  56 #include "geant321/gcphys.inc"
  57 #include "geant321/gctrak.inc"
  58 #include "geant321/gsecti.inc"
  59 #include "geant321/gctmed.inc"
  60 #include "geant321/gcunit.inc"
  61 #include "geant321/dimpar.inc"
  62
  63 #if !defined(CERNLIB_SINGLE)
  64 #include "geant321/finuct.inc"
  65 #endif
  66 #include "geant321/finuc.inc"
  67       REAL RNDM(1)
  68 #if !defined(CERNLIB_SINGLE)
  69       DOUBLE PRECISION AOCMBM, AMSS , ZTAR, RHO , ZLIN, ZLEL, ZLRAD,
  70      +ZUL
  71 #endif
  72       COMMON / FKMAPA / AOCMBM (MXXMDF), AMSS (MXXMDF), ZTAR   (MXXMDF),
  73      +                  RHO    (MXXMDF), ZLIN (MXXMDF), ZLEL   (MXXMDF),
  74      +                  ZLRAD  (MXXMDF), ZUL  (MXXMDF), MEDIUM (MXXRGN),
  75      +                  MULFLG (MXXMDF),IFCOMP(MXXMDF), MSSNUM (MXXMDF),
  76      +                  NREGS, NMATF, MTBSNM
  77 #if !defined(CERNLIB_SINGLE)
  78 #include "geant321/part2t.inc"
  79 #endif
  80 #include "geant321/part2.inc"
  81 #if !defined(CERNLIB_SINGLE)
  82 #include "geant321/comcont.inc"
  83 #endif
  84 #include "geant321/comcon.inc"
  85 #if !defined(CERNLIB_SINGLE)
  86 #include "geant321/fheavyt.inc"
  87 #endif
  88 #include "geant321/fheavy.inc"
  89 #include "geant321/paprop.inc"
  90 #if !defined(CERNLIB_SINGLE)
  91 #include "geant321/papropt.inc"
  92 #endif
  93 #include "geant321/gfkdis.inc"
  94 #if !defined(CERNLIB_SINGLE)
  95       DOUBLE PRECISION POO,EKE,TXI,TYI,TZI,AMM,WE,ONE,PGEANT,DMOD
  96 #endif
  97       PARAMETER (ONE=1)
  98       DIMENSION IGTOFL(49),IFLTOG(39),IHVTOG(6),ZSAMP(50)
  99       DATA IGTOFL / 0, 0, 0, 0, 0, 0,23,13,14,12, 15,16, 8, 1, 2,19, 0,
 100      +17,21,22, 20, 34, 36, 38, 9,18, 31, 32, 33, 35, 37, 39, 17*0/
 101
 102       DATA IFLTOG /14,15, 3, 2, 4, 4, 1,13,25, 5, 6,10, 8, 9,11,12,18,
 103      +26,16,21, 19,20, 7, 7*0, 27, 28, 29, 22, 30, 23, 31, 24, 32/
 104       DATA IHVTOG /13,14,45,46,49,47/
 105 *
 106       NP = 0
 107       NPHEAV = 0
 108 *
 109 *    Stopped particles:
 110 *    o also xi- are forced to decay
 111 *    o Neutral particles are sent to GHSTOP
 112 *    o pi+ and K+/K- are forced to decay
 113 *    o pi-, antiprotons and antineutrons are sent to FLUKA
 114 *      for annihilation (not here but later in this routine)
 115       IF ((IGF.EQ.1).OR.
 116      +      (GEKIN.EQ.0..AND.ITRTYP.EQ.3.AND.IPART.NE.25)) THEN
 117          CALL GHEISH
 118          IGF = 0
 119          GOTO 999
 120       ELSE IF (GEKIN.EQ.0..AND.
 121      +   (IPART.EQ.8.OR.IPART.EQ.12.OR.IPART.EQ.11.OR.IPART.EQ.23)) THEN
 122          CALL GDECAY
 123          NMEC=NMEC+1
 124          LMEC(NMEC)=5
 125          ISTOP=1
 126          GOTO 999
 127       ENDIF
 128 *
 129       IF (IFINIT(5) .EQ. 0) CALL FLINIT
 130       INT=0
 131       IJ=IGTOFL(IPART)
 132       IF(IJ.EQ.0) GOTO 110
 133       NMEC = NMEC + 1
 134       EKE = GEKIN
 135       TXI = VECT(4)
 136       TYI = VECT(5)
 137       TZI = VECT(6)
 138       DMOD = ONE/SQRT(TXI**2+TYI**2+TZI**2)
 139       TXI = TXI*DMOD
 140       TYI = TYI*DMOD
 141       TZI = TZI*DMOD
 142       WE  = 1.
 143       JMA = LQ(JMATE-NMAT)
 144       NCOMP = ABS (Q(JMA+11))
 145       AMM = Q(JMA+6)
 146       JMIXT = LQ(JMA-5)
 147
 148 *    Antiprotons, antineutrons and pi- are sent to
 149 *    eventv for annihilation
 150       IF (GEKIN.EQ.0..AND.
 151      +         (IPART.EQ.15.OR.IPART.EQ.9.OR.IPART.EQ.25)) THEN
 152          IF(NCOMP.LE.1) THEN
 153             AMSS(1) = Q(JMA+6)
 154             ZTAR(1) = Q(JMA+7)
 155             MSSNUM(1) = 0
 156             RHO(1) = Q(JMA+8)
 157          ELSE
 158             ZSAMP(1) = 0.
 159             DO 10 I=1,NCOMP
 160                ZSAMP(I+1) = ZSAMP(I) + Q(JMIXT+NCOMP+I)
 161    10       CONTINUE
 162             CALL GRNDM(RNDM,1)
 163             ZCONT=ZSAMP(NCOMP+1)*RNDM(1)
 164             DO 20 I=1,NCOMP
 165                IF(ZCONT.LE.ZSAMP(I+1)) GO TO 30
 166    20       CONTINUE
 167             I = NCOMP
 168    30       CONTINUE
 169             AMSS(1)   = Q(JMIXT+I)
 170             MSSNUM(1) = 0
 171             ZTAR(1)   = Q(JMIXT+NCOMP+I)
 172             RHO(1)    = Q(JMIXT+2*NCOMP+I)*DENS
 173          END IF
 174          EKE = 1E-9
 175          POO=SQRT(EKE*(EKE+2*AM(IJ)))
 176          CALL EVENTV(IJ,POO,EKE,TXI,TYI,TZI,WE,1)
 177          NMEC=NMEC+1
 178          LMEC(NMEC)=17
 179          GOTO 80
 180       ELSE IF (GEKIN.LE.CUTHAD .AND. ITRTYP.EQ.4) THEN
 181          DESTEP = DESTEP + GEKIN
 182          GEKIN  = 0.
 183          GETOT  = AMASS
 184          VECT(7) = 0.
 185          ISTOP = 1
 186          LMEC(NMEC)=30
 187          GO TO 110
 188       ENDIF
 189 *
 190       CALL GRNDM(RNDM,1)
 191       RNDEVT=RNDM(1)
 192
 193       IF (FSIG.LE.0.) THEN
 194          SRAT = 1.E7
 195       ELSE
 196          SRAT = SINE/FSIG
 197       ENDIF
 198       IF ( RNDEVT .GE. SRAT) THEN
 199          IF (GEKIN .GT. 0.02) THEN
 200             POO=SQRT(EKE*(EKE+2*AM(IJ)))
 201          ELSE
 202             GO TO 110
 203          END IF
 204          INT=1
 205          LMEC(NMEC)=13
 206          IF(NCOMP.LE.1) THEN
 207             CALL NUCREL(IJ,POO,EKE,TXI,TYI,TZI,AMM,WE)
 208          ELSE
 209             CALL GRNDM(RNDM,1)
 210             RCONT=ELXNOR*RNDM(1)
 211             DO 40  I=1,NCOMP
 212                IF(RCONT.LE.CABELX(I)) GO TO 50
 213    40       CONTINUE
 214             I=NCOMP
 215    50       CONTINUE
 216             CALL NUCREL(IJ,POO,EKE,TXI,TYI,TZI,ONE*Q(JMIXT+I),WE)
 217          END IF
 218       ELSE
 219          LMEC(NMEC)=20
 220          IF (IHADR.EQ.2) THEN
 221             ISTOP = 2
 222             DESTEP = DESTEP + GETOT
 223             GO TO 110
 224          ENDIF
 225          IF (GEKIN .GT. 0.02) THEN
 226             POO=SQRT(EKE*(EKE+2*AM(IJ)))
 227          ELSE
 228             IF ((IJ.EQ.2 .OR. IJ.EQ.9 .OR. IJ.EQ.14 .OR. IJ.EQ.16)
 229      +            .AND. GEKIN .GT. 0.0) THEN
 230                POO=SQRT(EKE*(EKE+2*AM(IJ)))
 231             ELSE
 232                NMEC=NMEC-1
 233                GO TO 110
 234             END IF
 235          END IF
 236          INT=2
 237          IF(NCOMP.LE.1) THEN
 238             AMSS(1) = Q(JMA+6)
 239             ZTAR(1) = Q(JMA+7)
 240             MSSNUM(1) = 0
 241             RHO(1) = Q(JMA+8)
 242          ELSE
 243             CALL GRNDM(RNDM,1)
 244             RCONT=ANXNOR*RNDM(1)
 245             DO 60  I=1,NCOMP
 246                IF(RCONT.LE.CABINX(I)) GO TO 70
 247    60       CONTINUE
 248             I=NCOMP
 249    70       CONTINUE
 250             AMSS(1)   = Q(JMIXT+I)
 251             MSSNUM(1) = 0
 252             ZTAR(1)   = Q(JMIXT+NCOMP+I)
 253             RHO(1)    = Q(JMIXT+2*NCOMP+I)*DENS
 254          END IF
 255          CALL EVENTV(IJ,POO,EKE,TXI,TYI,TZI,WE,1)
 256       END IF
 257 *
 258    80 IF(NP.EQ.1.AND.NPHEAV.EQ.0.AND.KPART(1).EQ.IJ) THEN
 259          VECT(4)=CXR(1)
 260          VECT(5)=CYR(1)
 261          VECT(6)=CZR(1)
 262          VECT(7)=SQRT(TKI(1)*(TKI(1)+2*AMASS))
 263          GETOT=TKI(1)+AMASS
 264          GEKIN=TKI(1)
 265       ELSE
 266          ISTOP=1
 267          NSTAK1 = MIN(NP,MXGKIN-NGKINE)
 268          IF(NP.GT.NSTAK1) THEN
 269             WRITE(CHMAIL,10000) NP-NSTAK1
 270             CALL GMAIL(0,0)
 271          ENDIF
 272          DO 90  K=1,NSTAK1
 273             NGKINE = NGKINE + 1
 274             IF (KPART(K) .EQ. 24 .OR. KPART(K) .EQ. 25) THEN
 275                KPART(K) = 19
 276                CALL GRNDM(RNDM,1)
 277                IF (RNDM(1) .GT. 0.5) KPART(K) = 12
 278             END IF
 279             IGEPAR = IFLTOG(KPART(K))
 280             JPA = LQ(JPART-IGEPAR)
 281             AGEMAS = Q(JPA+7)
 282             PGEANT = SQRT(TKI(K)*(TKI(K)+2*AGEMAS))
 283             GKIN(1,NGKINE)=CXR(K)*PGEANT
 284             GKIN(2,NGKINE)=CYR(K)*PGEANT
 285             GKIN(3,NGKINE)=CZR(K)*PGEANT
 286             GKIN(4,NGKINE)=TKI(K)+AGEMAS
 287             GKIN(5,NGKINE)=IGEPAR
 288             TOFD(NGKINE)=0.0
 289             GPOS(1,NGKINE) = VECT(1)
 290             GPOS(2,NGKINE) = VECT(2)
 291             GPOS(3,NGKINE) = VECT(3)
 292    90    CONTINUE
 293 *
 294          NSTAK2 = MIN(NPHEAV,MXGKIN-NGKINE)
 295          IF(NPHEAV.GT.NSTAK2) THEN
 296             WRITE(CHMAIL,10100) NPHEAV-NSTAK2
 297             CALL GMAIL(0,0)
 298          ENDIF
 299          DO 100 K=1,NSTAK2
 300             NGKINE = NGKINE + 1
 301             IGEPAR = IHVTOG(KHEAVY(K))
 302             JPA = LQ(JPART-IGEPAR)
 303             AGEMAS = Q(JPA+7)
 304             PGEANT = SQRT(TKHEAV(K)*(TKHEAV(K)+2*AGEMAS))
 305             GKIN(1,NGKINE)=CXHEAV(K)*PGEANT
 306             GKIN(2,NGKINE)=CYHEAV(K)*PGEANT
 307             GKIN(3,NGKINE)=CZHEAV(K)*PGEANT
 308             GKIN(4,NGKINE)=TKHEAV(K)+AGEMAS
 309             GKIN(5,NGKINE)=IGEPAR
 310             TOFD(NGKINE)=0.0
 311             GPOS(1,NGKINE) = VECT(1)
 312             GPOS(2,NGKINE) = VECT(2)
 313             GPOS(3,NGKINE) = VECT(3)
 314   100    CONTINUE
 315 *
 316          KCASE=NAMEC(12)
 317       END IF
 318   110 CONTINUE
 319       ZINTHA = GARNDM(DUMMY)
 320       SLHADR = SLENG
 321       STEPHA = 1.0E10
 322 10000 FORMAT(' **** FLUFIN: Stack overflow, ',I6,' particles lost')
 323 10100 FORMAT(' **** FLUFIN: Stack overflow, ',I6,
 324      +' heavy particles lost')
 325   999 END