GEANT321/fluka/fkflav.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:20:06  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/02 29/03/94  15.41.45  by  S.Giani
  11 *-- Author :
  12 *=== flavor ===========================================================*
  13 *
  14       SUBROUTINE FKFLAV(IT,LT,LL,E0,IV,RE,KFR1,KFR2,ISU,BETE,KFA1,KFA2,
  15      &                  KFA3,KFA4,IOPT)
  16
  17 #include "geant321/dblprc.inc"
  18 #include "geant321/dimpar.inc"
  19 #include "geant321/iounit.inc"
  20 *
  21 *----------------------------------------------------------------------*
  22 *                                                                      *
  23 *     New version from A. Ferrari ( 22 august 1990 ): it is almost the *
  24 *     the same as before, but with a few corrections important for isu *
  25 *     2 and 3 and also for isu=4 if it switchs to 100 continue         *
  26 *                                                                      *
  27 *----------------------------------------------------------------------*
  28 #include "geant321/part.inc"
  29 *
  30 *   The following are the masses of the quarks: the d quark mass is
  31 *   assumed to be the same as the u one. They are quite different from
  32 *   the last values from the particle data group, but any change
  33 *   can imply a change also in the BET parameter in common INPDAT
  34       PARAMETER ( UQUARM = 0.3D+00 )
  35       PARAMETER ( SQUARM = 0.5D+00 )
  36       PARAMETER ( CQUARM = 2.1D+00 )
  37       PARAMETER ( BQUARM = 5.0D+00 )
  38 *
  39       DIMENSION RE(*),KFR1(*),KFR2(*),IV(*)
  40       REAL RNDM(2)
  41 *
  42 C     CHOICE OF THE QUARK FLAVOUR
  43       IF (LT.EQ.1)WRITE(LUNOUT,288)IT,LT,LL,E0,ISU,BETE,KFA1,KFA2
  44   288 FORMAT(3I5,E12.4,I5,E12.4,2I5
  45      *,' FLAVOR IT,LT,LL,E0,ISU,BETE,KFA1,KFA2')
  46       I=IT
  47       J=IT-1
  48       IVA=1
  49       IVX=IV(I)
  50       IF (I .LE. 1) THEN
  51          IF (IOPT.EQ.2) THEN
  52             KX1 = KFA1
  53             KX2 = KFA2
  54          ELSE IF (IOPT.EQ.3.AND.LL.EQ.1) THEN
  55             KX1=KFA2
  56             KX2=0
  57          ELSE IF (IOPT.EQ.4 .AND. KFA1.LE.6 .AND. LL.EQ.1) THEN
  58             KX1=KFA2
  59             KX2=KFA3
  60          ELSE IF (IOPT.EQ.4 .AND. KFA1.GT.6 .AND. LL.EQ.0) THEN
  61             KX1=KFA2
  62             KX2=KFA3
  63          ELSE IF (IOPT.EQ.5 .AND. LL.EQ.0) THEN
  64             KX1=KFA3
  65             KX2=KFA4
  66          ELSE IF (IOPT.EQ.5 .AND. LL.EQ.1) THEN
  67             KX1=KFA1
  68             KX2=KFA2
  69          ELSE
  70             KX1=KFA1
  71             KX2=0
  72          END IF
  73          RX = E0
  74       ELSE
  75          KX1=KFR1(J)
  76          KX2=KFR2(J)
  77          RX =RE(J)
  78       END IF
  79       IF (KX1.GT.0.AND.KX2.GT.0) THEN
  80          BET=10.D+00
  81       ELSE
  82          BET=BETE
  83       END IF
  84       CALL GRNDM(RNDM,2)
  85       Z1=RNDM(1)
  86       Z2=RNDM(2)
  87       IF(ISU.EQ.4) GO TO 300
  88       IF(ISU.EQ.3) GO TO 200
  89       IF(ISU.EQ.2) GO TO 100
  90 C     U FLAVOUR
  91          KF1=1
  92          KF2=1
  93       GO TO 20
  94 C     U/D FLAVOURS
  95   100 CONTINUE
  96          IF (KX1.EQ.1.OR.KX1.EQ.7) THEN
  97             IIAA=1
  98          ELSE IF (KX1.EQ.2.OR.KX1.EQ.8) THEN
  99             IIAA=2
 100          ELSE
 101             IIAA=0
 102          END IF
 103   110 CONTINUE
 104          IF (IIAA .EQ. 1) THEN
 105             PD=0.6666666666666667D+00
 106             PU=0.3333333333333333D+00
 107          ELSE IF (IIAA .EQ. 2) THEN
 108             PU=0.6666666666666667D+00
 109             PD=0.3333333333333333D+00
 110          ELSE
 111             PU=0.5D+00
 112             PD=0.5D+00
 113          END IF
 114          PS=0.D+00
 115          PC=0.D+00
 116          IF (Z1 .LE. PD) THEN
 117             KF1=2
 118          ELSE
 119             KF1=1
 120          END IF
 121          IF (Z2 .LE. PD) THEN
 122             KF2=2
 123          ELSE
 124             KF2=1
 125          END IF
 126       GO TO 20
 127 C     U/D/S FLAVOURS
 128   200 CONTINUE
 129          IF (KX1.EQ.1.OR.KX1.EQ.7) THEN
 130             IIAA=1
 131          ELSE IF (KX1.EQ.2.OR.KX1.EQ.8) THEN
 132             IIAA=2
 133          ELSE
 134             IIAA=0
 135          END IF
 136   210 CONTINUE
 137          IF (RX .LE. 1.019D+00 ) GO TO 110
 138          X1=RX
 139          X2=UQUARM
 140          PU=BETA(X1,X2,BET)
 141          X2=SQUARM
 142          PS=BETA(X1,X2,BET)
 143          PTOT=2.D+00*PU+PS
 144          PU1=PU/PTOT
 145          PS =PS/PTOT
 146          PC=0.D+00
 147          IF (IIAA .EQ. 1) THEN
 148             PU=0.6666666666666667D+00*PU1
 149             PD=2.D+00*PU1-PU
 150          ELSE IF (IIAA .EQ. 2) THEN
 151             PD=0.6666666666666667D+00*PU1
 152             PU=2.D+00*PU1-PD
 153          ELSE
 154             PU=PU1
 155             PD=PU
 156          END IF
 157          IF (Z1 .LE. PU) THEN
 158             KF1 = 1
 159          ELSE IF ( Z1 .LE. PU + PD ) THEN
 160             KF1 = 2
 161          ELSE
 162             KF1 = 3
 163          END IF
 164          IF (Z2 .LE. PU) THEN
 165             KF2 = 1
 166          ELSE IF ( Z2 .LE. PU + PD ) THEN
 167             KF2 = 2
 168          ELSE
 169             KF2 = 3
 170          END IF
 171       GO TO 20
 172 C     U/D/S/C FLAVOUR
 173   300 CONTINUE
 174          GO TO (11,12,13,14,14,11,12,13,14,14),IVX
 175    11    CONTINUE
 176             IF (KX1.EQ.4.OR.KX1.EQ.10) THEN
 177                GM=AM(129)
 178             ELSE
 179                GM=AM(127)
 180             END IF
 181             IF (KX1.EQ.1.OR.KX1.EQ.7) THEN
 182                IIAA=1
 183             ELSE IF (KX1.EQ.2.OR.KX1.EQ.8) THEN
 184                IIAA=2
 185             ELSE
 186                IIAA=0
 187             END IF
 188          GO TO 15
 189    12    CONTINUE
 190             IF (KX1.EQ.4.OR.KX1.EQ.10) THEN
 191                GM=AM(170)
 192             ELSE
 193                GM=AM(127)
 194             END IF
 195             IIAA=0
 196          GO TO 15
 197    13    CONTINUE
 198 *  |  |  +-------------------------------------------------------------*
 199 *  |  |  | The following if replaces the cards:
 200 *  |  |  |      GM=3.85D0
 201 *  |  |  |      IF(KX1.EQ.4.AND.KX2.EQ.4) GM=4.89D0
 202 *  |  |  |      IF(KX1.NE.4.AND.KX2.NE.4) GM=2.770D0
 203 *  |  |  |      IF(KX1.EQ.10.AND.KX2.EQ.10) GM=4.89D0
 204 *  |  |  |      IF(KX1.NE.10.AND.KX2.NE.10) GM=2.770D0
 205 *  |  |  | It is completely equivalent except for the combination
 206 *  |  |  | 4-4 which now gives GM = 4.89, while in the original
 207 *  |  |  | coding gave GM = 2.77, because of the last condition
 208 *  |  |  | always overrides the first one (it seems to be a mistake)
 209             IF (KX1 .EQ. 4 .OR. KX1 .EQ. 10 ) THEN
 210                IF ( KX2 .EQ. KX1 ) THEN
 211                   GM = AM(170)
 212                ELSE
 213                   GM = AM(169)
 214                END IF
 215             ELSE IF (KX2 .EQ. 4 .OR. KX2 .EQ. 10 ) THEN
 216                GM = AM(169)
 217             ELSE
 218                GM = AM(166)
 219             END IF
 220             IIAA=0
 221          GO TO 15
 222    14    CONTINUE
 223 *  |  |  +-------------------------------------------------------------*
 224 *  |  |  | The following if replaces the cards:
 225 *  |  |  |      GM=3.684D0
 226 *  |  |  |      IF(KX1.NE.4.AND.KX2.NE.4) GM=2.140D0
 227 *  |  |  |      IF(KX1.NE.10.AND.KX2.NE.10) GM=2.140D0
 228 *  |  |  | It is equivalent: only for the combinations
 229 *  |  |  | 4-10 and 10-4 we get GM = 3.684
 230 *  |  |  | It is not clear if it is correct since 4-x,x-4 (x.ne.10),
 231 *  |  |  | 10-x,x-10 (x.ne.4) give GM = 2.14
 232             IF ((KX1.EQ.4.AND.KX2.EQ.10).OR.(KX1.EQ.10.AND.KX2.EQ.4))
 233      &         THEN
 234                GM = AM(129)
 235             ELSE
 236                GM = AM(127)
 237             END IF
 238             IF (IVX.EQ.4.OR.IVX.EQ.9) THEN
 239                KAXI=KX1
 240             ELSE IF (IVX.EQ.5.OR.IVX.EQ.10) THEN
 241                KAXI=KX2
 242             ELSE
 243 *  |  |  |   Kaxi = 0 added for completeness, maybe it is useless
 244                KAXI=0
 245             END IF
 246             IF (KAXI.EQ.1.OR.KAXI.EQ.7) THEN
 247                IIAA=1
 248             ELSE IF (KAXI.EQ.2.OR.KAXI.EQ.8) THEN
 249                IIAA=2
 250             ELSE
 251                IIAA=0
 252             END IF
 253          GO TO 15
 254    15    CONTINUE
 255          IF (RX .LE. GM) GO TO 200
 256          X1=RX
 257          X2=UQUARM
 258          PU=BETA(X1,X2,BET)
 259          X2=SQUARM
 260          PS=BETA(X1,X2,BET)
 261          X2=CQUARM
 262          PC=BETA(X1,X2,BET)
 263          PTOT=2.D+00*PU+PS+PC
 264          PU1=PU/PTOT
 265          PS=PS/PTOT
 266          PC=PC/PTOT
 267          IF (IIAA .EQ. 1) THEN
 268             PU=0.6666666666666667D+00*PU1
 269             PD=2.D+00*PU1-PU
 270          ELSE IF (IIAA .EQ. 2) THEN
 271             PD=0.6666666666666667D+00*PU1
 272             PU=2.D+00*PU1-PD
 273          ELSE
 274             PU=PU1
 275             PD=PU
 276          END IF
 277          IF (Z1 .LE. PU) THEN
 278             KF1 = 1
 279          ELSE IF ( Z1 .LE. PU + PD ) THEN
 280             KF1 = 2
 281          ELSE IF ( Z1 .LE. PU + PD + PS ) THEN
 282             KF1 = 3
 283          ELSE
 284             KF1 = 4
 285          END IF
 286          IF (Z2 .LE. PU) THEN
 287             KF2 = 1
 288          ELSE IF ( Z2 .LE. PU + PD ) THEN
 289             KF2 = 2
 290          ELSE IF ( Z2 .LE. PU + PD + PS ) THEN
 291             KF2 = 3
 292          ELSE
 293             KF2 = 4
 294          END IF
 295       GO TO 20
 296    20 CONTINUE
 297 C*****CHOICE OF THE QUARKFLAVOURS IN DEPENDENCE OF THE VERTEX IV
 298       IVX=IV(I)
 299       GO TO (1,2,3,4,5,1,2,3,4,5),IVX
 300     1 CONTINUE
 301          IF (LL.EQ.1) THEN
 302             KFR1(I)=KF1+6
 303          ELSE
 304             KFR1(I)=KF1
 305          END IF
 306          KFR2(I)=0
 307       GO TO 30
 308     2 CONTINUE
 309          IF (LL.EQ.1) THEN
 310             KFR1(I)=KF1
 311             KFR2(I)=KF2
 312          ELSE
 313             KFR1(I)=KF1+6
 314             KFR2(I)=KF2+6
 315          END IF
 316       GO TO 30
 317     3 CONTINUE
 318          KFR2(I)=0
 319          IF (LL.EQ.1) THEN
 320             KFR1(I)=KF1+6
 321          ELSE
 322             KFR1(I)=KF1
 323          END IF
 324       GO TO 30
 325     4 CONTINUE
 326          IF (LL.EQ.1) THEN
 327             KFR1(I)=KF1
 328          ELSE
 329             KFR1(I)=KF1+6
 330          END IF
 331          KFR2(I)=KX2
 332       GO TO 30
 333     5 CONTINUE
 334          KFR1(I)=KX1
 335          IF (LL.EQ.1) THEN
 336             KFR2(I)=KF2
 337          ELSE
 338             KFR2(I)=KF2+6
 339          END IF
 340       GO TO 30
 341    30 CONTINUE
 342       IF(LT.EQ.0) GO TO 80
 343       WRITE(LUNOUT,6)PU,PD,PS,PC,KX1,KX2
 344     6 FORMAT(1H0,' FLAVOR PU,PD,PS,PC,KX1,KX2',4F8.4,2I5)
 345       IF(I.EQ.1) GO TO 40
 346       WRITE(LUNOUT,60)IV(I),LL,KFR1(J),KFR2(J),KFR1(I),KFR2(I)
 347       GO TO 50
 348    40 WRITE(LUNOUT,70)IV(I),LL,KFA1,KFA2,KFR1(I),KFR2(I)
 349    50 CONTINUE
 350    60 FORMAT(1H0,22HIV,LL,Q1A,Q2A,Q1N,Q2N=,6I3)
 351    70 FORMAT(1H0,'IV(I),LL,KFA1,KFA2,KFR1(I),KFR2(I)=',6I3 )
 352    80 CONTINUE
 353       RETURN
 354       END