GEANT321/ghrout/ghetun.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.2  1997/04/29 09:54:58  japost
   6 * Re-dimensioning and filling common arrays to fix bug/error.
   7 *
   8 * More details can be found in this error report, which caused the changes:
   9 *
  10 * Date: Tue, 15 Apr 1997 19:03:49 +0200
  11 * Simone's forward:
  12 * ----------
  13 *
  14 * > From: Dieter Heck <heck@ik3.fzk.de>
  15 * > To: SGIANI@cernvm.cern.ch
  16 * > Subject: GHEISHA_ERROR
  17 * > Date: Tuesday, April 08, 1997 2:09 PM
  18 * >
  19 * >
  20 * > Dear Simone Giani,
  21 * >
  22 * > again I have detected an error within one of the GHEISHA routines.
  23 * > It is the routine GHETUN. This routine contains two arrays keeping the
  24 * > strangenes numbers (SNUM) and baryon numbers (BNUM) of the particles to
  25 * > be treated. Both arrays are dimensioned with 32, while in GHEISHA also
  26 * > OMAGA(-) and ANTI-OMEAG(-) particles are admitted with particle codes 33
  27 * > and 34. In case of treatment of omega particles, the array boundaries
  28 * > are exceeded, which gives an error stop when running with the bounds
  29 * > check option. To correct this error, the arrays must be dimensionend
  30 * > with 34, and the DATA statements have to be filled with the
  31 * > corresponding values:
  32 * >
  33 * >       DIMENSION RNDM(4),SNUM(34),BNUM(34),REDDEC(7)
  34 * >       DATA SNUM/9*0.,1.,0.,0.,-1.,4*0.,-1.,1.,-1.,-1.,-1.,1.,1.,1.,
  35 * >      $          -2.,-2.,2.,2.,3*0.,-3.,3./
  36 * >       DATA BNUM/13*0.,1.,-1.,1.,-1.,1.,-1.,1.,1.,1.,-1.,-1.,-1.,
  37 * >      $          1.,1.,-1.,-1.,2.,3.,4.,1.,-1./
  38 * >
  39 * > Hopefully this correction is helpful for the GHEISHA users.
  40 * >
  41 * > Kind regards,
  42 * >
  43 * >   Dieter Heck
  44 * >
  45 * > =========================================================================
  46 * > From:   Dieter Heck,   Institut fuer Kernphysik 3,
  47 * > Forschungszentrum Karlsruhe,  Postfach 3640,  D-76021 Karlsruhe,  Germany
  48 * > Tel: (49) 7247-82-3777   Fax: (49) 7247-82-4075   e-mail: heck@ik3.fzk.de
  49 * > =========================================================================
  50 *
  51 * Revision 1.1.1.1  1995/10/24 10:21:12  cernlib
  52 * Geant
  53 *
  54 *
  55 #include "geant321/pilot.h"
  56 *CMZ :  3.21/02 29/03/94  15.41.40  by  S.Giani
  57 *-- Author :
  58       SUBROUTINE GHETUN(NT)
  59 C**
  60 C** TUNING OF THE HIGH ENERGY COLLISION MODEL:
  61 C**
  62 C** 1. AVOID THAT PI0 IS LEADING PARTICLE.
  63 C** 2. SOME FINE-TUNING FOR THE NUMBER OF PRODUCED PROTONS AND
  64 C**    NEUTRONS.
  65 C** 3. INTRODUCE A FLAVOUR DEPENDENT CORRECTION FOR SINGLE PARTICLE
  66 C**    SPECTRA.
  67 C** 4. FINE-TUNING OF LEADING PARTICLE SPECTRA AND MOMENTUM
  68 C**    CONSERVATION.
  69 C**
  70 #include "geant321/mxgkgh.inc"
  71 #include "geant321/s_consts.inc"
  72 #include "geant321/s_curpar.inc"
  73 #include "geant321/s_result.inc"
  74 #include "geant321/s_prntfl.inc"
  75 #include "geant321/s_blank.inc"
  76 C
  77 C
  78       DIMENSION RNDM(4),SNUM(34),BNUM(34),REDDEC(7)
  79       DATA SNUM/9*0.,1.,0.,0.,-1.,4*0.,-1.,1.,-1.,-1.,-1.,1.,1.,1.,
  80      $          -2.,-2.,2.,2.,3*0.,-3.,3./
  81       DATA BNUM/13*0.,1.,-1.,1.,-1.,1.,-1.,1.,1.,1.,-1.,-1.,-1.,
  82      $          1.,1.,-1.,-1.,2.,3.,4.,1.,-1./
  83
  84
  85 #if !defined(CERNLIB_SINGLE)
  86       DOUBLE PRECISION AHMF,BHMF
  87 #endif
  88 C**
  89       MX=MXGKPV-10
  90       MX1=MX+1
  91       MX2=MX+2
  92       MX3=MX+3
  93       MX4=MX+4
  94       MX5=MX+5
  95       MX6=MX+6
  96       MX7=MX+7
  97       MX8=MX+8
  98       MX9=MX+9
  99       NT1=NT
 100       IF(NT1.GT.MXGKPV-10) NT1=MXGKPV-10
 101       NT=NT1
 102 C
 103       CALL GRNDM(RNDM,1)
 104       IF(EK.LT.(25.+RNDM(1)*75.)) GOTO 15
 105 C
 106 C**  IF PI0 IS THE HIGHEST MOMENTUM PARTICLE, INTERCHANGE IT WITH A
 107 C**  CHARGED PION.
 108 C
 109       CALL GRNDM(RNDM,4)
 110       REDEN = -0.7 + 0.29*LOG10(EK)
 111       REDAT = 1. - 0.4000*LOG10(ATNO2)
 112       PMAX  = -200.
 113       PMAPIP= -200.
 114       PMAPI0= -200.
 115       PMAPIM= -200.
 116       IPMAX = 0
 117       MAXPIP= 0
 118       MAXPI0= 0
 119       MAXPIM= 0
 120       IF(RNDM(1).GT.(ATNO2/100.-0.28).AND.ABS(NCH).GT.0.5) THEN
 121          DO 46 I=1,NT1
 122             IPHMF=IFIX(PV(8,I)+0.1)
 123             CALL LENGTX(I,PPP)
 124             IF(PPP.GT.PMAX) THEN
 125                PMAX=PPP
 126                IPMAX=I
 127             ENDIF
 128             IF(IPHMF.EQ.7) THEN
 129                IF(PPP.GT.PMAPIP) THEN
 130                   PMAPIP=PPP
 131                   MAXPIP=I
 132                ENDIF
 133             ENDIF
 134             IF(IPHMF.EQ.8) THEN
 135                IF(PPP.GT.PMAPI0) THEN
 136                   PMAPI0=PPP
 137                   MAXPI0=I
 138                ENDIF
 139             ENDIF
 140             IF(IPHMF.EQ.9) THEN
 141                IF(PPP.GT.PMAPIM) THEN
 142                   PMAPIM=PPP
 143                   MAXPIM=I
 144                ENDIF
 145             ENDIF
 146    46    CONTINUE
 147       ENDIF
 148 C**
 149 C**   SOME ADDITIONAL TUNING OF THE NUMBER OF GREY TRACK PARTICLES
 150 C**
 151       IF(NT1.GT.2) THEN
 152       DO 47 I=3,NT1
 153          IPHMF=IFIX(PV(8,I)+0.1)
 154          IF(IPHMF.EQ.14.OR.IPHMF.EQ.16.OR.IPHMF.GE.30) THEN
 155             CALL LENGTX(I,PPP)
 156             IF(PPP.LT.1.5) THEN
 157                IF(RNDM(2).LT.REDEN.OR.RNDM(3).LT.REDAT) THEN
 158                   PV(1,I) = 0.
 159                   PV(2,I) = 0.
 160                   PV(3,I) = 0.
 161                   PV(4,I) = ABS(PV(5,I))
 162                ENDIF
 163             ENDIF
 164          ENDIF
 165    47 CONTINUE
 166       ENDIF
 167 C**
 168       IF(MAXPI0.EQ.0)    GOTO 10
 169       IF(PMAPI0.LT.PMAX) GOTO 10
 170       IF(RNDM(4).LT.PMAPI0/P) THEN
 171       IF(NCH.GT.0.5.AND.MAXPIP.NE.0) THEN
 172          DO 49 J=5,10
 173            PV(J,MX1)=PV(J,MAXPI0)
 174            PV(J,MAXPI0)=PV(J,MAXPIP)
 175            PV(J,MAXPIP)=PV(J,MX1)
 176    49    CONTINUE
 177       ENDIF
 178       IF(NCH.LT.-0.5.AND.MAXPIM.NE.0) THEN
 179          DO 56 J=5,10
 180            PV(J,MX1)=PV(J,MAXPI0)
 181            PV(J,MAXPI0)=PV(J,MAXPIM)
 182            PV(J,MAXPIM)=PV(J,MX1)
 183    56    CONTINUE
 184       ENDIF
 185       ENDIF
 186 C
 187    10 CONTINUE
 188 C**
 189 C** CHECK TOTAL BARYON- NUMBER AND
 190 C** SKIP ZERO MOMENTUM PARTICLES
 191 C**
 192       BNTOT=-BNUM(IPART)-ATNO2
 193       DO 57 I=1,NT1
 194          IPHMF=IFIX(PV(8,I)+0.1)
 195          BNTOT=BNTOT+BNUM(IPHMF)
 196    57 CONTINUE
 197       BNTOT=1.+BNTOT/ATNO2
 198       IF(ATNO2.LT.1.5) BNTOT=0.
 199       CALL GRNDM(RNDM,1)
 200       IF(ATNO2.GT.(75.+RNDM(1)*25.)) BNTOT=0.
 201 C**
 202       II=0
 203       DO 12 I=1,NT1
 204          CALL LENGTX(I,PPP)
 205          IF(PPP.GT.1.E-6) THEN
 206             IPHMF=IFIX(PV(8,I)+0.1)
 207             IF(BNTOT.GT.0.3) THEN
 208             IF(IPHMF.EQ.14.OR.IPHMF.EQ.16.OR.IPHMF.GE.30) THEN
 209                CALL GRNDM(RNDM,1)
 210                IF(RNDM(1).LT.0.5.AND.PPP.LT.1.2) GOTO 12
 211             ENDIF
 212             ENDIF
 213             II=II+1
 214             DO 11 J=1,10
 215                PV(J,II)=PV(J,I)
 216    11       CONTINUE
 217          ENDIF
 218    12 CONTINUE
 219       NT1=II
 220       NT=NT1
 221 C**
 222 C**   EXACT MOMENTUM CONSERVATION AND SOME CORRECTIONS FOR SINGLE
 223 C**   PARTICLE SPECTRA AT HIGH ENERGIES ONLY
 224 C
 225    15 PV(1,MX1) = P*PX
 226       PV(2,MX1) = P*PY
 227       PV(3,MX1) = P*PZ
 228       PV(4,MX1) = EN
 229       PV(5,MX1) = ABS(AMAS)
 230       PV(6,MX1) = NCH
 231       PV(1,MX2) = 0.
 232       PV(2,MX2) = 0.
 233       PV(3,MX2) = 0.
 234       PV(4,MX2) = MP
 235       PV(5,MX2) = MP
 236       PV(6,MX2) = 0.
 237 C
 238       IF(NPRT(4)) THEN
 239          WRITE(NEWBCD,2000)
 240          WRITE(NEWBCD,2001) MX1,(PV(J,MX1),J=1,6)
 241          WRITE(NEWBCD,2001) MX2,(PV(J,MX2),J=1,6)
 242       ENDIF
 243 C
 244       DO 58 J=1,10
 245          PV(J,MX9) = 0.
 246    58 CONTINUE
 247       CALL ADD(MX1,MX2,MX3)
 248       CALL LOR(MX1,MX3,MX4)
 249       CALL LOR(MX2,MX3,MX5)
 250       LEDPAR=0
 251       REDPAR=0.
 252       GESPAR=0.
 253       SNUM1=SNUM(IPART)
 254       IF(IPART.EQ.11.OR.IPART.EQ.12) THEN
 255         CALL GRNDM(RNDM,1)
 256         SNUM1=1.
 257         IF(RNDM(1).LT.0.5) SNUM1=-1.
 258       ENDIF
 259       DO 20 I=1,NT1
 260          IPHMF=IFIX(PV(8,I)+0.1)
 261          IF(IPHMF.LE.6.OR.IPHMF.GT.32) GOTO 19
 262          CALL LENGTX(I,PPP)
 263          IF(PPP.LT.1.E-3) GOTO 19
 264          CALL LOR(I,MX3,MX6)
 265          CALL ANG(MX4,MX6,COST,TETA)
 266          SNUM2=SNUM(IPHMF)
 267          IF(IPHMF.EQ.11.OR.IPHMF.EQ.12) THEN
 268             CALL GRNDM(RNDM,1)
 269             SNUM2=1.
 270             IF(RNDM(1).LT.0.5) SNUM2=-1.
 271          ENDIF
 272          IF(COST.GT.0.) THEN
 273             HFMAS=ABS(AMAS)
 274             REDDEC(1)=ABS(HFMAS    -ABS(PV(5,I)))
 275             REDDEC(2)=ABS(NCH-PV(6,I))
 276             REDDEC(3)=ABS(SNUM1      -SNUM2)
 277             REDDEC(4)=ABS(BNUM(IPART)-BNUM(IPHMF))
 278          ELSE
 279             HFMAS=MP
 280             REDDEC(1)=ABS(HFMAS     -ABS(PV(5,I)))
 281             REDDEC(2)=ABS(ZNO2/ATNO2-PV(6,I))
 282             REDDEC(3)=ABS(SNUM2)
 283             REDDEC(4)=ABS(1.-BNUM(IPHMF))
 284          ENDIF
 285          REDDEC(6)=REDDEC(1)+REDDEC(2)+REDDEC(3)+REDDEC(4)
 286          SBQWGT=REDDEC(6)
 287          IF(HFMAS.LT.0.2) THEN
 288             SBQWGT=(SBQWGT-2.5)*0.10
 289             IF(IPHMF.EQ.8) SBQWGT=0.15
 290          ELSE IF (HFMAS.LT.0.6) THEN
 291             SBQWGT=(SBQWGT-3.0)*0.10
 292          ELSE
 293             SBQWGT=(SBQWGT-2.0)*0.10
 294             IF(IPHMF.EQ.8) SBQWGT=0.15
 295          ENDIF
 296          CALL LENGTX(MX6,PPP)
 297          IF(SBQWGT.GT.0. .AND. PPP.GT.1.E-6) THEN
 298          PLHMF=PPP*COST
 299          PTHMF=PPP*SQRT(1.-COST*COST)
 300          PLHMF=PLHMF*(1.-SBQWGT)
 301          CALL CROSS3(MX4,MX6,MX8)
 302          CALL CROSS3(MX8,MX4,MX8)
 303          CALL LENGTX(MX4,PPP)
 304          PV(1,MX7)=PV(1,MX4)*PLHMF/PPP
 305          PV(2,MX7)=PV(2,MX4)*PLHMF/PPP
 306          PV(3,MX7)=PV(3,MX4)*PLHMF/PPP
 307          CALL LENGTX(MX8,PPP)
 308          PV(1,MX8)=PV(1,MX8)*PTHMF/PPP
 309          PV(2,MX8)=PV(2,MX8)*PTHMF/PPP
 310          PV(3,MX8)=PV(3,MX8)*PTHMF/PPP
 311          CALL ADD3(MX7,MX8,MX6)
 312          CALL LENGTX(MX6,PPP)
 313          AHMF=PPP
 314          BHMF=PV(5,I)
 315          PV(4,MX6)=DSQRT(AHMF**2+BHMF**2)
 316 C
 317          IF(NPRT(4))
 318      $      WRITE(NEWBCD,3001) I,(PV(J,I),J=1,8),SBQWGT
 319 C
 320          CALL LOR(MX6,MX5,I)
 321 C
 322          IF(NPRT(4))
 323      $      WRITE(NEWBCD,3001) I,(PV(J,I),J=1,8),SBQWGT
 324          ENDIF
 325 C
 326          IF(IPHMF.EQ.8) GOTO 19
 327          CALL SUB3(MXGKPV,I,MX8)
 328          CALL LENGTX(MX8,PPP)
 329          REDDEC(5) = PPP/P
 330          REDDEC(7)=REDDEC(5)*2./3. + REDDEC(6)/12.
 331          REDDEC(7) =  1.-REDDEC(7)
 332          IF(REDDEC(7) .LT.  0.) REDDEC(7) =  0.
 333          GESPAR=GESPAR+REDDEC(7)
 334          IF(REDDEC(6).LE.3.75) THEN
 335           IF(REDDEC(7) .GT.  REDPAR) THEN
 336              LEDPAR=I
 337              REDPAR=REDDEC(7)
 338           ENDIF
 339          ENDIF
 340          IF(NPRT(4))
 341      $      WRITE(NEWBCD,2001) I,(PV(J,I),J=1,6),PV(8,I),REDDEC
 342 C
 343   19     CALL ADD3(MX9,I,MX9)
 344 C
 345   20  CONTINUE
 346       IF(NPRT(4))
 347      $   WRITE(NEWBCD,1001) LEDPAR,REDPAR,GESPAR
 348 C**
 349 C** APPLY CORRECTION ON LEADING PARTICLE
 350 C**
 351       CALL GHEPEC(LEDPAR)
 352 C**
 353       RETURN
 354  1001 FORMAT(1H0,'SEARCH FOR LEADING PARTICLE, WEIGHT, TOTAL WEIGHT ',
 355      $ I5,3X,2F10.4)
 356  2000 FORMAT(1H ,'MOMENTUM CONSERVATION AT HIGH ENERGIES: ')
 357  2001 FORMAT(1H ,I3,2X,7F8.3/1H ,5X,7F8.3)
 358  3001 FORMAT(1H ,I3,2X,5F8.3,F5.1,F8.3,F5.1,F8.3)
 359       END