GEANT321/giface/gheish.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:21:15  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/02 29/03/94  15.41.38  by  S.Giani
  11 *-- Author :
  12       SUBROUTINE GHEISH
  13 C
  14 C *** MAIN STEERING FOR HADRON SHOWER DEVELOPMENT ***
  15 C *** NVE 15-JUN-1988 CERN GENEVA ***
  16 C
  17 C CALLED BY : GUHADR (USER ROUTINE)
  18 C ORIGIN : F.CARMINATI, H.FESEFELDT
  19 C                       ROUTINES : CALIM  16-SEP-1987
  20 C                                  SETRES 19-AUG-1985
  21 C                                  INTACT 06-OCT-1987
  22 C
  23 #include "geant321/gcbank.inc"
  24 #include "geant321/gcjloc.inc"
  25 #include "geant321/gccuts.inc"
  26 #include "geant321/gcflag.inc"
  27 #include "geant321/gckine.inc"
  28 #include "geant321/gcking.inc"
  29 #include "geant321/gcmate.inc"
  30 #include "geant321/gcphys.inc"
  31 #include "geant321/gctmed.inc"
  32 #include "geant321/gctrak.inc"
  33 #include "geant321/gsecti.inc"
  34 #include "geant321/gcunit.inc"
  35 C --- GHEISHA COMMONS ---
  36 #include "geant321/mxgkgh.inc"
  37 #include "geant321/s_blankp.inc"
  38 #include "geant321/s_consts.inc"
  39 #include "geant321/s_event.inc"
  40 #include "geant321/s_prntfl.inc"
  41 C
  42 C --- "NEVENT" CHANGED TO "KEVENT" IN COMMON /CURPAR/ DUE TO CLASH ---
  43 C --- WITH VARIABLE "NEVENT" IN GEANT COMMON ---
  44 C
  45       PARAMETER (MXGKCU=MXGKGH)
  46       COMMON /CURPAR /WEIGHT(10),DDELTN,IFILE,IRUN,NEVT,KEVENT,SHFLAG,
  47      $                ITHST,ITTOT,ITLST,IFRND,TOFCUT,CMOM(5),CENG(5),
  48      $                RS,S,ENP(10),NP,NM,NN,NR,NO,NZ,IPA(MXGKCU),
  49      $                ATNO2,ZNO2
  50 C
  51 C --- "IPART" CHANGED TO "KPART" IN COMMON /RESULT/ DUE TO CLASH ---
  52 C --- WITH VARIABLE "IPART" IN GEANT COMMON ---
  53 C
  54       COMMON /RESULT/ XEND,YEND,ZEND,RCA,RCE,AMAS,NCH,TOF,PX,PY,PZ,
  55      $                USERW,INTCT,P,EN,EK,AMASQ,DELTN,ITK,NTK,KPART,IND,
  56      $                LCALO,ICEL,SINL,COSL,SINP,COSP,
  57      $                XOLD,YOLD,ZOLD,POLD,PXOLD,PYOLD,PZOLD,
  58      $                XSCAT,YSCAT,ZSCAT,PSCAT,PXSCAT,PYSCAT,PZSCAT
  59                       REAL NCH,INTCT
  60 C
  61 C --- "ABSL(21)" CHANGED TO "ABSLTH(21)" IN COMMON /MAT/ DUE TO CLASH ---
  62 C --- WITH VARIABLE "ABSL" IN GEANT COMMON ---
  63 C
  64       COMMON /MAT/ LMAT,
  65      $             DEN(21),RADLTH(21),ATNO(21),ZNO(21),ABSLTH(21),
  66      $             CDEN(21),MDEN(21),X0DEN(21),X1DEN(21),RION(21),
  67      $             MATID(21),MATID1(21,24),PARMAT(21,10),
  68      $             IFRAT,IFRAC(21),FRAC1(21,10),DEN1(21,10),
  69      $             ATNO1(21,10),ZNO1(21,10)
  70 C
  71       DIMENSION IPELOS(35)
  72       SAVE IDEOL
  73 C
  74 C --- TRANSFER GEANT CUT-OFFS INTO GHEISHA VALUES ---
  75       DIMENSION CUTS(5)
  76       EQUIVALENCE (CUTS(1),CUTGAM)
  77       DIMENSION RNDM(1)
  78 C
  79 #include "geant321/pcodim.inc"
  80 #include "geant321/pcodat.inc"
  81 C
  82 C --- DENOTE STABLE PARTICLES ACCORDING TO GHEISHA CODE ---
  83 C --- STABLE : GAMMA, NEUTRINO, ELECTRON, PROTON AND HEAVY FRAGMENTS ---
  84 C --- WHEN STOPPING THESE PARTICLES ONLY LOOSE THEIR KINETIC ENERGY ---
  85       DATA IPELOS/
  86      $             1,   1,   0,   1,   0,   0,   0,   0,
  87      $             0,   0,   0,   0,   0,   1,   0,   0,
  88      $             0,   0,   0,   0,   0,   0,   0,   0,
  89      $             0,   0,   0,   0,   0,   1,   1,   1,
  90      $             0,   0,   1/
  91 C
  92 C --- LOWERBOUND OF KINETIC ENERGY BIN IN N CROSS-SECTION TABLES ---
  93       DATA TEKLOW /0.0001/
  94 C
  95 C --- KINETIC ENERGY TO SWITCH FROM "CASN" TO "GNSLWD" FOR N CASCADE ---
  96       DATA SWTEKN /0.05/
  97 C
  98       DATA IDEOL/0/
  99 C
 100 C --- INITIALIZE RELEVANT GHEISHA VARIABLES IN CASE NOT DONE ALREADY ---
 101       IF (IFINIT(4) .EQ. 0) CALL GHEINI
 102 C
 103 C --- SET THE INTERACTION MECHANISM TO "HADR" ---
 104       KCASE=NAMEC(12)
 105 C
 106 C --- SET GHEISHA PRINTING FLAGS ACCORDING TO "DEBUG" STEERING CARD --
 107       IF (IDEOL .EQ. IDEBUG) GO TO 9000
 108 C
 109       IF (IDEBUG .NE. 1) GO TO 9001
 110 C
 111 C --- SET SELECTED DEBUGGING FLAGS ---
 112       DO 9002 LL=1,10
 113       IF ((ISWIT(LL) .LE. 100) .OR. (ISWIT(LL) .GT. 110)) GO TO 9002
 114       JJ=ISWIT(LL)-100
 115       NPRT(JJ)=.TRUE.
 116  9002 CONTINUE
 117       GO TO 9000
 118 C
 119 C --- NO DEBUGGING SELECTED ---
 120  9001 CONTINUE
 121       DO 9003 LL=1,10
 122       NPRT(LL)=.FALSE.
 123  9003 CONTINUE
 124       IDEOL=IDEBUG
 125 C
 126  9000 CONTINUE
 127 C
 128 C --- SET THE GHEISHA PARTICLE TYPE TO THE ONE OF GEANT ---
 129       IF(IPART.GT.48) THEN
 130          IF(ISTOP.EQ.0) GOTO 9999
 131          JPA = LQ(JPART-IPART)
 132          AMAS=Q(JPA+7)
 133          NCH =Q(JPA+8)
 134          KPART=-IPART
 135          GOTO 107
 136       ENDIF
 137       NETEST=IKPART(KPART)
 138       IF ((NETEST .EQ. IPART) .OR. (ISTOP .NE. 0)) GO TO 9004
 139 C
 140       PRINT 8881,IPART,KPART,ISTOP
 141  8881 FORMAT(' *GHEISH* IPART,KPART = ',2(I3,1X),' ISTOP = ',I3/
 142      $ ' *GHEISH* ======> PARTICLE TYPES DO NOT MATCH <=======')
 143       STOP
 144 C
 145  9004 CONTINUE
 146       KPART=KIPART(IPART)
 147       KKPART=KPART
 148       AMAS=RMASS(KPART)
 149       NCH=RCHARG(KPART)
 150 C
 151 C --- TRANSPORT THE TRACK NUMBER TO GHEISHA AND INITIALISE SOME NUMBERS
 152  107  NTK=ITRA
 153       INTCT=0.0
 154       NEXT=1
 155       NTOT=0
 156       TOF=0.0
 157 C
 158 C --- FILL RESULT COMMON FOR THIS TRACK WITH GEANT VALUES ---
 159 C --- CALIM CODE ---
 160       XEND=VECT(1)
 161       YEND=VECT(2)
 162       ZEND=VECT(3)
 163       PX=VECT(4)
 164       PY=VECT(5)
 165       PZ=VECT(6)
 166       USERW=UPWGHT
 167 C --- SETRES CODE ---
 168       P=VECT(7)
 169       AMASQ=AMAS*AMAS
 170       EN=SQRT(AMASQ+P*P)
 171       EK=ABS(EN-ABS(AMAS))
 172       ENOLD=EN
 173 C
 174       SINL=0.0
 175       COSL=1.0
 176       SINP=0.0
 177       COSP=1.0
 178 C
 179       IF (ABS(P) .LE. 1.0E-10) GO TO 1
 180       SINL=PZ
 181       COSL=SQRT(ABS(1.0-SINL**2))
 182 C
 183  1    CONTINUE
 184       CALL GRNDM(RNDM,1)
 185       PHI=RNDM(1)*TWPI
 186       IF ((PX .EQ. 0.0) .AND. (PY .EQ. 0.0)) GOTO 3
 187       IF (ABS(PX) .LT. 1.E-10) GOTO 2
 188       PHI=ATAN2(PY,PX)
 189       GOTO 3
 190 C
 191  2    CONTINUE
 192       IF (PY .GT. 0.0) PHI=PI/2.0
 193       IF (PY .LE. 0.0) PHI=3.0*PI/2.0
 194 C
 195  3    CONTINUE
 196       SINP=SIN(PHI)
 197       COSP=COS(PHI)
 198 C
 199 C --- SET GHEISHA INDEX FOR THE CURRENT MEDIUM ALWAYS TO 1 ---
 200       IND=1
 201 C
 202 C --- TRANSFER GLOBAL MATERIAL CONSTANTS FOR CURRENT MEDIUM ---
 203 C --- DETAILED DATA FOR COMPOUNDS IS OBTAINED VIA ROUTINE COMPO ---
 204       ATNO(IND+1)=A
 205       ZNO(IND+1)=Z
 206       DEN(IND+1)=DENS
 207       RADLTH(IND+1)=RADL
 208       ABSLTH(IND+1)=ABSL
 209 C
 210 C --- SETUP PARMAT FOR PHYSICS STEERING ---
 211       PARMAT(IND+1,5)=0.0
 212       PARMAT(IND+1,8)=IPFIS
 213       PARMAT(IND+1,9)=0.0
 214       PARMAT(IND+1,10)=0.0
 215       JTMN=LQ(JTM)
 216       IF (JTMN .LE. 0) GO TO 4
 217       PARMAT(IND+1,5)=Q(JTMN+26)
 218  4    CONTINUE
 219 C
 220 C --- CHECK WHETHER PARTICLE IS STOPPING OR NOT ---
 221       IF (ISTOP .EQ. 0) GO TO 5
 222 C
 223       IF (NPRT(9)) PRINT 1000,KPART
 224  1000 FORMAT(' *GHEISH* STOPPING GHEISHA PARTICLE ',I3)
 225       CALL GHSTOP
 226 C --- IN CASE OF DECAY OF PARTICLE OR USER PARTICLE ==> RETURN ---
 227       IF (LMEC(NMEC) .EQ. 5 .OR. KPART .LT. 0) GO TO 9999
 228 C --- IN CASE OF HAD. INT. WITH GENERATION OF SEC. ==> GO TO 40 ---
 229       IF (IHADR .NE. 2) GO TO 40
 230 C --- ALSO DEPOSIT REST MASS ENERGY FOR IN-STABLE PARTICLES ---
 231       IF (IPELOS(KPART) .EQ. 0) DESTEP=DESTEP+ABS(RMASS(KPART))
 232       GO TO 9999
 233   5   CONTINUE
 234 C
 235 C --- INDICATE LIGHT (<= PI) AND HEAVY PARTICLES (HISTORICALLY) ---
 236 C --- CALIM CODE ---
 237       J=2
 238       TEST=RMASS(7)-0.001
 239       IF (ABS(AMAS) .LT. TEST) J=1
 240 C
 241 C *** DIVISION INTO VARIOUS INTERACTION CHANNELS DENOTED BY "INT" ***
 242 C THE CONVENTION FOR "INT" IS THE FOLLOWING
 243 C
 244 C INT  = -1 REACTION CROSS SECTIONS NOT YET TABULATED/PROGRAMMED
 245 C      =  0 NO INTERACTION
 246 C      =  1 ELEASTIC SCATTERING
 247 C      =  2 INELASTIC SCATTERING
 248 C      =  3 NUCLEAR FISSION WITH INELEASTIC SCATTERING
 249 C      =  4 NEUTRON CAPTURE
 250 C
 251 C --- INTACT CODE ---
 252       KK=ABS(Q(JMA+11))
 253       ALAM1=0.0
 254       CALL GRNDM(RNDM,1)
 255       RAT=RNDM(1)*ALAM
 256       NMEC=NMEC+1
 257       ATNO2=A
 258       ZNO2 =Z
 259 C
 260       DO 6 K=1,KK
 261       IF (KK .LE. 0) GO TO 6
 262 C
 263       IF (KK .EQ. 1) GO TO 7
 264       ATNO2=Q(JMIXT+K)
 265       ZNO2 =Q(JMIXT+K+KK)
 266 C
 267  7    CONTINUE
 268 C
 269 C --- TRY FOR ELASTIC SCATTERING ---
 270       INT=1
 271       LMEC(NMEC)=13
 272       ALAM1=ALAM1+AIEL(K)
 273       IF (RAT .LT. ALAM1) GO TO 8
 274 C
 275 C --- TRY FOR INELASTIC SCATTERING ---
 276       INT=2
 277       LMEC(NMEC)=20
 278       ALAM1=ALAM1+AIIN(K)
 279       IF (RAT .LT. ALAM1) GO TO 8
 280 C
 281 C --- TRY FOR NUCLEAR FISSION WITH INELASTIC SCATTERING ---
 282       INT=3
 283       LMEC(NMEC)=15
 284       ALAM1=ALAM1+AIFI(K)
 285       IF (RAT .LT. ALAM1) GO TO 8
 286 C
 287 C --- TRY FOR NEUTRON CAPTURE ---
 288       INT=4
 289       LMEC(NMEC)=18
 290       ALAM1=ALAM1+AICA(K)
 291       IF (RAT .LT. ALAM1) GO TO 8
 292 C
 293  6    CONTINUE
 294 C --- NO REACTION SELECTED ==> ELASTIC SCATTERING ---
 295       INT=1
 296       LMEC(NMEC)=13
 297 C
 298 C *** TAKE ACTION ACCORDING TO SELECTED REACTION CHANNEL ***
 299 C --- FOLLOWING CODE IS A TRANSLATION OF "CALIM" INTO GEANT JARGON ---
 300 C
 301  8    CONTINUE
 302       IF (NPRT(9)) PRINT 1001,INT
 303  1001 FORMAT(' *GHEISH* INTERACTION TYPE CHOSEN INT = ',I3)
 304 C
 305 C --- IN CASE OF NO INTERACTION OR UNKNOWN CROSS SECTIONS ==> DONE ---
 306       IF (INT .LE. 0) GO TO 40
 307 C
 308 C --- IN CASE OF NON-ELASTIC SCATTERING AND NO GENERATION OF SEC. ---
 309 C --- PARTICLES DEPOSIT TOTAL PARTICLE ENERGY AND RETURN ---
 310       IF ((INT .EQ. 1) .OR. (IHADR .NE. 2)) GO TO 9
 311       ISTOP=2
 312       DESTEP=DESTEP+EN
 313       NGKINE=0
 314       GO TO 9999
 315 C
 316  9    CONTINUE
 317       IF (INT .NE. 4) GO TO 10
 318 C
 319 C --- NEUTRON CAPTURE ---
 320       IF (NPRT(9)) PRINT 2000
 321  2000 FORMAT(' *GHEISH* ROUTINE CAPTUR WILL BE CALLED')
 322       ISTOP=1
 323       CALL CAPTUR(NOPT)
 324       GO TO 40
 325 C
 326  10   CONTINUE
 327       IF (INT .NE. 3) GO TO 11
 328 C --- NUCLEAR FISSION ---
 329       IF (NPRT(9)) PRINT 2001
 330  2001 FORMAT(' *GHEISH* ROUTINE FISSIO WILL BE CALLED')
 331       ISTOP=1
 332       TKIN=FISSIO(EK)
 333       GO TO 40
 334 C
 335  11   CONTINUE
 336 C
 337 C --- ELASTIC AND INELASTIC SCATTERING ---
 338       PV( 1,MXGKPV)=P*PX
 339       PV( 2,MXGKPV)=P*PY
 340       PV( 3,MXGKPV)=P*PZ
 341       PV( 4,MXGKPV)=EN
 342       PV( 5,MXGKPV)=AMAS
 343       PV( 6,MXGKPV)=NCH
 344       PV( 7,MXGKPV)=TOF
 345       PV( 8,MXGKPV)=KPART
 346       PV( 9,MXGKPV)=0.
 347       PV(10,MXGKPV)=USERW
 348 C
 349 C --- ADDITIONAL PARAMETERS TO SIMULATE FERMI MOTION AND EVAPORATION ---
 350       DO 111 JENP=1,10
 351          ENP(JENP)=0.
 352  111  CONTINUE
 353       ENP(5)=EK
 354       ENP(6)=EN
 355       ENP(7)=P
 356 C
 357       IF (INT .NE. 1) GO TO 12
 358 C
 359 C *** ELASTIC SCATTERING PROCESSES ***
 360 C
 361 C --- ONLY NUCLEAR INTERACTIONS FOR HEAVY FRAGMENTS ---
 362       IF ((KPART .GE. 30) .AND. (KPART .LE. 32)) GO TO 35
 363 C
 364 C --- NORMAL ELASTIC SCATTERING FOR LIGHT MEDIA ---
 365       IF (ATNO2 .LT. 1.5) GO TO 35
 366 C
 367 C --- COHERENT ELASTIC SCATTERING FOR HEAVY MEDIA ---
 368       IF (NPRT(9)) PRINT 2002
 369  2002 FORMAT(' *GHEISH* ROUTINE COSCAT WILL BE CALLED')
 370       CALL COSCAT
 371       GO TO 40
 372 C
 373 C *** NON-ELASTIC SCATTERING PROCESSES ***
 374  12   CONTINUE
 375 C
 376 C --- ONLY NUCLEAR INTERACTIONS FOR HEAVY FRAGMENTS ---
 377       IF ((KPART .GE. 30) .AND. (KPART .LE. 32)) GO TO 35
 378 C
 379 C *** USE SOMETIMES NUCLEAR REACTION ROUTINE "NUCREC" FOR LOW ENERGY ***
 380 C *** PROTON AND NEUTRON SCATTERING ***
 381       CALL GRNDM(RNDM,1)
 382       TEST1=RNDM(1)
 383       TEST2=4.5*(EK-0.01)
 384       IF ((KPART .EQ. 14) .AND. (TEST1 .GT. TEST2)) GO TO 85
 385       IF ((KPART .EQ. 16) .AND. (TEST1 .GT. TEST2)) GO TO 86
 386 C
 387 C *** FERMI MOTION AND EVAPORATION ***
 388       TKIN=CINEMA(EK)
 389       PV( 9,MXGKPV)=TKIN
 390       ENP(5)=EK+TKIN
 391 C --- CHECK FOR LOWERBOUND OF EKIN IN CROSS-SECTION TABLES ---
 392       IF (ENP(5) .LE. TEKLOW) ENP(5)=TEKLOW
 393       ENP(6)=ENP(5)+ABS(AMAS)
 394       ENP(7)=(ENP(6)-AMAS)*(ENP(6)+AMAS)
 395       ENP(7)=SQRT(ABS(ENP(7)))
 396       TKIN=FERMI(ENP(5))
 397       ENP(5)=ENP(5)+TKIN
 398 C --- CHECK FOR LOWERBOUND OF EKIN IN CROSS-SECTION TABLES ---
 399       IF (ENP(5) .LE. TEKLOW) ENP(5)=TEKLOW
 400       ENP(6)=ENP(5)+ABS(AMAS)
 401       ENP(7)=(ENP(6)-AMAS)*(ENP(6)+AMAS)
 402       ENP(7)=SQRT(ABS(ENP(7)))
 403       TKIN=EXNU(ENP(5))
 404       ENP(5)=ENP(5)-TKIN
 405 C --- CHECK FOR LOWERBOUND OF EKIN IN CROSS-SECTION TABLES ---
 406       IF (ENP(5) .LE. TEKLOW) ENP(5)=TEKLOW
 407       ENP(6)=ENP(5)+ABS(AMAS)
 408       ENP(7)=(ENP(6)-AMAS)*(ENP(6)+AMAS)
 409       ENP(7)=SQRT(ABS(ENP(7)))
 410 C
 411 C *** IN CASE OF ENERGY ABOVE CUT-OFF LET THE PARTICLE CASCADE ***
 412       TEST=ABS(CHARGE)
 413       IF ((TEST .GT. 1.0E-10) .AND. (ENP(5) .GT. CUTHAD)) GO TO 35
 414       IF ((TEST .LE. 1.0E-10) .AND. (ENP(5) .GT. CUTNEU)) GO TO 35
 415 C
 416 C --- SECOND CHANCE FOR ANTI-BARYONS DUE TO POSSIBLE ANNIHILATION ---
 417       IF ((AMAS .GE. 0.0) .OR. (KPART .LE. 14)) GO TO 13
 418       ANNI=1.3*P
 419       IF (ANNI .GT. 0.4) ANNI=0.4
 420       CALL GRNDM(RNDM,1)
 421       TEST=RNDM(1)
 422       IF (TEST .GT. ANNI) GO TO 35
 423 C
 424 C *** PARTICLE WITH ENERGY BELOW CUT-OFF ***
 425 C --- ==> ONLY NUCLEAR EVAPORATION AND QUASI-ELASTIC SCATTERING ---
 426  13   CONTINUE
 427 C
 428       ISTOP=3
 429 C
 430       IF (NPRT(9)) PRINT 1002,KPART,EK,EN,P,ENP(5),ENP(6),ENP(7)
 431  1002 FORMAT(' *GHEISH* ENERGY BELOW CUT-OFF FOR GHEISHA PARTICLE ',I3/
 432      $ ' EK,EN,P,ENP(5),ENP(6),ENP(7) = ',6(G12.5,1X))
 433 C
 434       IF ((KPART .NE. 14) .AND. (KPART .NE. 16)) GO TO 14
 435       IF (KPART .EQ. 16) GO TO 86
 436 C
 437 C --- SLOW PROTON ---
 438  85   CONTINUE
 439       IF (NPRT(9)) PRINT 2003,EK,KPART
 440  2003 FORMAT(' *GHEISH* ROUTINE NUCREC WILL BE CALLED',
 441      $ ' EK = ',G12.5,' GEV  KPART = ',I3)
 442       CALL NUCREC(NOPT,2)
 443 C
 444       IF (NOPT .NE. 0) GO TO 50
 445 C
 446       IF (NPRT(9)) PRINT 2004,EK,KPART
 447  2004 FORMAT(' *GHEISH* ROUTINE COSCAT WILL BE CALLED',
 448      $ ' EK = ',G12.5,' GEV  KPART = ',I3)
 449       CALL COSCAT
 450       GO TO 40
 451 C
 452 C --- SLOW NEUTRON ---
 453  86   CONTINUE
 454       IF (NPRT(9)) PRINT 2015
 455       NUCFLG=0
 456       CALL GNSLWD(NUCFLG,INT,NFL,TEKLOW)
 457       IF (NUCFLG .NE. 0) GO TO 50
 458       GO TO 40
 459 C
 460 C --- OTHER SLOW PARTICLES ---
 461  14   CONTINUE
 462       IPA(1)=KPART
 463 C --- DECIDE FOR PROTON OR NEUTRON TARGET ---
 464       IPA(2)=16
 465       CALL GRNDM(RNDM,1)
 466       TEST1=RNDM(1)
 467       TEST2=ZNO2/ATNO2
 468       IF (TEST1 .LT. TEST2) IPA(2)=14
 469       AVERN=0.0
 470       NFL=1
 471       IF (IPA(2) .EQ. 16) NFL=2
 472       IPPP=KPART
 473       IF (NPRT(9)) PRINT 2005
 474  2005 FORMAT(' *GHEISH* ROUTINE TWOB WILL BE CALLED')
 475       CALL TWOB(IPPP,NFL,AVERN)
 476       GOTO 40
 477 C
 478 C --- INITIALISATION OF CASCADE QUANTITIES ---
 479  35   CONTINUE
 480 C
 481 C *** CASCADE GENERATION ***
 482 C --- CALCULATE FINAL STATE MULTIPLICITY AND LONGITUDINAL AND ---
 483 C --- TRANSVERSE MOMENTUM DISTRIBUTIONS ---
 484 C
 485 C --- FIXED PARTICLE TYPE TO STEER THE CASCADE ---
 486       KKPART=KPART
 487 C
 488 C --- NO CASCADE FOR LEPTONS ---
 489       IF (KKPART .LE. 6) GO TO 9999
 490 C
 491 C *** WHAT TO DO WITH "NEW PARTICLES" FOR GHEISHA ?????? ***
 492 C --- RETURN FOR THE TIME BEING ---
 493       IF (KKPART .GE. 35) GO TO 9999
 494 C
 495 C --- CASCADE OF HEAVY FRAGMENTS
 496       IF ((KKPART .GE. 30) .AND. (KKPART .LE. 32)) GO TO 390
 497 C
 498 C --- INITIALIZE THE IPA ARRAY ---
 499       CALL VZERO(IPA(1),MXGKCU)
 500 C
 501 C --- CASCADE OF OMEGA - AND OMEGA - BAR ---
 502       IF (KKPART .EQ. 33) GO TO 330
 503       IF (KKPART .EQ. 34) GO TO 331
 504 C
 505       NVEPAR=KKPART-17
 506       IF (NVEPAR .LE. 0) GO TO 15
 507       GO TO (318,319,320,321,322,323,324,325,326,327,328,329),NVEPAR
 508 C
 509  15   CONTINUE
 510       NVEPAR=KKPART-6
 511       GO TO (307,308,309,310,311,312,313,314,315,316,317,318),NVEPAR
 512 C
 513 C --- PI+ CASCADE ---
 514  307  CONTINUE
 515       IF (NPRT(9)) PRINT 2006
 516  2006 FORMAT(' *GHEISH* ROUTINE CASPIP WILL BE CALLED')
 517       CALL CASPIP(J,INT,NFL)
 518       GO TO 40
 519 C
 520 C --- PI0 ==> NO CASCADE ---
 521  308  CONTINUE
 522       GO TO 40
 523 C
 524 C --- PI- CASCADE ---
 525  309  CONTINUE
 526       IF (NPRT(9)) PRINT 2007
 527  2007 FORMAT(' *GHEISH* ROUTINE CASPIM WILL BE CALLED')
 528       CALL CASPIM(J,INT,NFL)
 529       GO TO 40
 530 C
 531 C --- K+ CASCADE ---
 532  310  CONTINUE
 533       IF (NPRT(9)) PRINT 2008
 534  2008 FORMAT(' *GHEISH* ROUTINE CASKP WILL BE CALLED')
 535       CALL CASKP(J,INT,NFL)
 536       GO TO 40
 537 C
 538 C --- K0 CASCADE ---
 539  311  CONTINUE
 540       IF (NPRT(9)) PRINT 2009
 541  2009 FORMAT(' *GHEISH* ROUTINE CASK0 WILL BE CALLED')
 542       CALL CASK0(J,INT,NFL)
 543       GO TO 40
 544 C
 545 C --- K0 BAR CASCADE ---
 546  312  CONTINUE
 547       IF (NPRT(9)) PRINT 2010
 548  2010 FORMAT(' *GHEISH* ROUTINE CASK0B WILL BE CALLED')
 549       CALL CASK0B(J,INT,NFL)
 550       GO TO 40
 551 C
 552 C --- K- CASCADE ---
 553  313  CONTINUE
 554       IF (NPRT(9)) PRINT 2011
 555  2011 FORMAT(' *GHEISH* ROUTINE CASKM WILL BE CALLED')
 556       CALL CASKM(J,INT,NFL)
 557       GO TO 40
 558 C
 559 C --- PROTON CASCADE ---
 560  314  CONTINUE
 561       IF (NPRT(9)) PRINT 2012
 562  2012 FORMAT(' *GHEISH* ROUTINE CASP WILL BE CALLED')
 563       CALL CASP(J,INT,NFL)
 564       GO TO 40
 565 C
 566 C --- PROTON BAR CASCADE ---
 567  315  CONTINUE
 568       IF (NPRT(9)) PRINT 2013
 569  2013 FORMAT(' *GHEISH* ROUTINE CASPB WILL BE CALLED')
 570       CALL CASPB(J,INT,NFL)
 571       GO TO 40
 572 C
 573 C --- NEUTRON CASCADE ---
 574  316  CONTINUE
 575       NUCFLG=0
 576       IF (EK .GT. SWTEKN) THEN
 577          CALL CASN(J,INT,NFL)
 578          IF (NPRT(9)) PRINT 2014
 579  2014 FORMAT(' *GHEISH* ROUTINE CASN WILL BE CALLED')
 580       ELSE
 581          CALL GNSLWD(NUCFLG,INT,NFL,TEKLOW)
 582          IF (NPRT(9)) PRINT 2015
 583  2015 FORMAT(' *GHEISH* ROUTINE GNSLWD WILL BE CALLED')
 584       ENDIF
 585       IF (NUCFLG .NE. 0) GO TO 50
 586       GO TO 40
 587 C
 588 C --- NEUTRON BAR CASCADE ---
 589  317  CONTINUE
 590       IF (NPRT(9)) PRINT 2016
 591  2016 FORMAT(' *GHEISH* ROUTINE CASNB WILL BE CALLED')
 592       CALL CASNB(J,INT,NFL)
 593       GO TO 40
 594 C
 595 C --- LAMBDA CASCADE ---
 596  318  CONTINUE
 597       IF (NPRT(9)) PRINT 2017
 598  2017 FORMAT(' *GHEISH* ROUTINE CASL0 WILL BE CALLED')
 599       CALL CASL0(J,INT,NFL)
 600       GO TO 40
 601 C
 602 C --- LAMBDA BAR CASCADE ---
 603  319  CONTINUE
 604       IF (NPRT(9)) PRINT 2018
 605  2018 FORMAT(' *GHEISH* ROUTINE CASAL0 WILL BE CALLED')
 606       CALL CASAL0(J,INT,NFL)
 607       GO TO 40
 608 C
 609 C --- SIGMA + CASCADE ---
 610  320  CONTINUE
 611       IF (NPRT(9)) PRINT 2019
 612  2019 FORMAT(' *GHEISH* ROUTINE CASSP WILL BE CALLED')
 613       CALL CASSP(J,INT,NFL)
 614       GO TO 40
 615 C
 616 C --- SIGMA 0 ==> NO CASCADE ---
 617  321  CONTINUE
 618       GO TO 40
 619 C
 620 C --- SIGMA - CASCADE ---
 621  322  CONTINUE
 622       IF (NPRT(9)) PRINT 2020
 623  2020 FORMAT(' *GHEISH* ROUTINE CASSM WILL BE CALLED')
 624       CALL CASSM(J,INT,NFL)
 625       GO TO 40
 626 C
 627 C --- SIGMA + BAR CASCADE ---
 628  323  CONTINUE
 629       IF (NPRT(9)) PRINT 2021
 630  2021 FORMAT(' *GHEISH* ROUTINE CASASP WILL BE CALLED')
 631       CALL CASASP(J,INT,NFL)
 632       GO TO 40
 633 C
 634 C --- SIGMA 0 BAR ==> NO CASCADE ---
 635  324  CONTINUE
 636       GO TO 40
 637 C
 638 C --- SIGMA - BAR CASCADE ---
 639  325  CONTINUE
 640       IF (NPRT(9)) PRINT 2022
 641  2022 FORMAT(' *GHEISH* ROUTINE CASASM WILL BE CALLED')
 642       CALL CASASM(J,INT,NFL)
 643       GO TO 40
 644 C
 645 C --- XI 0 CASCADE ---
 646  326  CONTINUE
 647       IF (NPRT(9)) PRINT 2023
 648  2023 FORMAT(' *GHEISH* ROUTINE CASX0 WILL BE CALLED')
 649       CALL CASX0(J,INT,NFL)
 650       GO TO 40
 651 C
 652 C --- XI - CASCADE ---
 653  327  CONTINUE
 654       IF (NPRT(9)) PRINT 2024
 655  2024 FORMAT(' *GHEISH* ROUTINE CASXM WILL BE CALLED')
 656       CALL CASXM(J,INT,NFL)
 657       GO TO 40
 658 C
 659 C --- XI 0 BAR CASCADE ---
 660  328  CONTINUE
 661       IF (NPRT(9)) PRINT 2025
 662  2025 FORMAT(' *GHEISH* ROUTINE CASAX0 WILL BE CALLED')
 663       CALL CASAX0(J,INT,NFL)
 664       GO TO 40
 665 C
 666 C --- XI - BAR CASCADE ---
 667  329  CONTINUE
 668       IF (NPRT(9)) PRINT 2026
 669  2026 FORMAT(' *GHEISH* ROUTINE CASAXM WILL BE CALLED')
 670       CALL CASAXM(J,INT,NFL)
 671       GO TO 40
 672 C
 673 C --- OMEGA - CASCADE ---
 674  330  CONTINUE
 675       IF (NPRT(9)) PRINT 2027
 676  2027 FORMAT(' *GHEISH* ROUTINE CASOM WILL BE CALLED')
 677       CALL CASOM(J,INT,NFL)
 678       GO TO 40
 679 C
 680 C --- OMEGA - BAR CASCADE ---
 681  331  CONTINUE
 682       IF (NPRT(9)) PRINT 2028
 683  2028 FORMAT(' *GHEISH* ROUTINE CASAOM WILL BE CALLED')
 684       CALL CASAOM(J,INT,NFL)
 685       GO TO 40
 686 C
 687 C --- HEAVY FRAGMENT CASCADE ---
 688  390  CONTINUE
 689       IF (NPRT(9)) PRINT 2090
 690  2090 FORMAT(' *GHEISH* ROUTINE CASFRG WILL BE CALLED')
 691       NUCFLG=0
 692       CALL CASFRG(NUCFLG,INT,NFL)
 693       IF (NUCFLG .NE. 0) GO TO 50
 694 C
 695 C *** CHECK WHETHER THERE ARE NEW PARTICLES GENERATED ***
 696  40   CONTINUE
 697       IF ((NTOT .NE. 0) .OR. (KKPART .NE. KPART)) GO TO 50
 698 C
 699 C --- NO SECONDARIES GENERATED AND PARTICLE IS STILL THE SAME ---
 700 C --- ==> COPY EVERYTHING BACK IN THE CURRENT GEANT STACK ---
 701       NGKINE=0
 702       TOFG=TOFG+TOF*0.5E-10
 703 C --- In case of crazy momentum value ==> no change to GEANT stack ---
 704       IF (P .LT. 0.) GO TO 41
 705       VECT(4)=PX
 706       VECT(5)=PY
 707       VECT(6)=PZ
 708       VECT(7)=P
 709       GETOT=EN
 710       GEKIN=EK
 711 C --- CHECK KINETIC ENERGY ---
 712       CALL GEKBIN
 713       EDEP=ABS(ENOLD-EN)
 714       RMASSI=EN-EK
 715       IF (NPRT(9) .AND. (EN .GT. ENOLD))
 716      $ PRINT 8888,EDEP,ENOLD,EN,EK,RMASSI
 717  8888 FORMAT(' *GHEISH* EDEP,ENOLD,EN,EK,M = ',5(G12.5,1X)/
 718      $ ' *GHEISH* =======> EDEP WOULD BE NEGATIVE <========')
 719       IF (ISTOP .EQ. 0) DESTEP=DESTEP+EDEP
 720 C
 721 C --- RE-INITIALIZE THE PROBABILITY FOR HADRONIC INTERACTION ---
 722  41   CONTINUE
 723       CALL GRNDM(RNDM,1)
 724       IF ((RNDM(1) .LE. 0.) .OR. (RNDM(1) .GE. 1.)) GO TO 41
 725       ZINTHA=-LOG(RNDM(1))
 726       SLHADR=SLENG
 727       STEPHA=1.0E10
 728 C
 729       NVEDUM=KIPART(IPART)
 730       IF (NPRT(9)) PRINT 1003,NTOT,IPART,KPART,KKPART,NVEDUM
 731  1003 FORMAT(' *GHEISH* NO SEC. GEN. NTOT,IPART,KPART,KKPART,KIPART = ',
 732      $ 5(I3,1X)/
 733      $ ' CURRENT PARTICLE ON THE STACK AGAIN')
 734       GO TO 9999
 735 C
 736 C *** CURRENT PARTICLE IS NOT THE SAME AS IN THE BEGINNING OR/AND ***
 737 C *** ONE OR MORE SECONDARIES HAVE BEEN GENERATED ***
 738  50   CONTINUE
 739 C
 740       NVEDUM=KIPART(IPART)
 741       IF (NPRT(9)) PRINT 1004,NTOT,IPART,KPART,KKPART,NVEDUM
 742  1004 FORMAT(' *GHEISH* SEC. GEN. NTOT,IPART,KPART,KKPART,KIPART = ',
 743      $ 5(I3,1X))
 744 C
 745 C --- INITIAL PARTICLE TYPE HAS BEEN CHANGED ==> PUT NEW TYPE ON ---
 746 C --- THE GEANT TEMPORARY STACK ---
 747 C
 748 C --- MAKE CHOICE BETWEEN K0 LONG / K0 SHORT ---
 749       IF ((KPART .NE. 11) .AND. (KPART .NE. 12)) GO TO 52
 750       CALL GRNDM(RNDM,1)
 751       KPART=11.5+RNDM(1)
 752 C
 753  52   CONTINUE
 754       ITY=IKPART(KPART)
 755       LNVE=LQ(JPART-ITY)
 756       IF (LNVE .LE. 0) PRINT 1234,NTOT,ITY,LNVE
 757  1234 FORMAT('0*GHEISH* 1234 NTOT,ITY,LNVE = ',3(I10,1X))
 758       IF (LNVE .LE. 0) STOP
 759       IF (ISTOP .EQ. 0) ISTOP=1
 760 C
 761 C --- IN CASE THE NEW PARTICLE IS A NEUTRINO ==> FORGET IT ---
 762       IF (KPART .EQ. 2) GO TO 60
 763 C
 764 C --- PUT PARTICLE ON THE STACK ---
 765       GKIN(1,1)=PX*P
 766       GKIN(2,1)=PY*P
 767       GKIN(3,1)=PZ*P
 768       GKIN(4,1)=SQRT(P*P+RMASS(KPART)**2)
 769       GKIN(5,1)=ITY
 770       TOFD(1)=TOF*0.5E-10
 771       NGKINE = 1
 772       GPOS(1,1) = VECT(1)
 773       GPOS(2,1) = VECT(2)
 774       GPOS(3,1) = VECT(3)
 775 C
 776       IF (NPRT(9)) PRINT 1005,ITY,NGKINE
 777  1005 FORMAT(' *GHEISH* GEANT PART. ',I3,' PUT ONTO STACK AT POS. ',I3)
 778 C
 779 C *** CHECK WHETHER SECONDARIES HAVE BEEN GENERATED AND COPY THEM ***
 780 C *** ALSO ON THE GEANT STACK ***
 781  60   CONTINUE
 782 C
 783 C --- ALL QUANTITIES ARE TAKEN FROM THE GHEISHA STACK WHERE THE ---
 784 C --- CONVENTION IS THE FOLLOWING ---
 785 C
 786 C EVE(INDEX+ 1)= X
 787 C EVE(INDEX+ 2)= Y
 788 C EVE(INDEX+ 3)= Z
 789 C EVE(INDEX+ 4)= NCAL
 790 C EVE(INDEX+ 5)= NCELL
 791 C EVE(INDEX+ 6)= MASS
 792 C EVE(INDEX+ 7)= CHARGE
 793 C EVE(INDEX+ 8)= TOF
 794 C EVE(INDEX+ 9)= PX
 795 C EVE(INDEX+10)= PY
 796 C EVE(INDEX+11)= PZ
 797 C EVE(INDEX+12)= TYPE
 798 C
 799       IF (NTOT .LE. 0) GO TO 9999
 800 C
 801 C --- ONE OR MORE SECONDARIES HAVE BEEN GENERATED ---
 802       DO 61 L=1,NTOT
 803       INDEX=(L-1)*12
 804       JND=EVE(INDEX+12)
 805 C
 806 C --- MAKE CHOICE BETWEEN K0 LONG / K0 SHORT ---
 807       IF ((JND .NE. 11) .AND. (JND .NE. 12)) GO TO 63
 808       CALL GRNDM(RNDM,1)
 809       JND=11.5+RNDM(1)
 810 C
 811 C --- FORGET ABOUT NEUTRINOS ---
 812  63   CONTINUE
 813       IF (JND .EQ. 2) GO TO 61
 814 C
 815 C --- SWITH TO GEANT QUANTITIES ---
 816       ITY=IKPART(JND)
 817       JTY=LQ(JPART-ITY)
 818       IF (JTY .LE. 0) PRINT 1235,NTOT,ITY,JTY
 819  1235 FORMAT('0*GHEISH* 1235 NTOT,ITY,JTY = ',3(I10,1X))
 820       IF (JTY .LE. 0) STOP
 821 *     ITRT=Q(JTY+6)
 822       PLX=EVE(INDEX+9)
 823       PLY=EVE(INDEX+10)
 824       PLZ=EVE(INDEX+11)
 825       ELT=SQRT(PLX*PLX+PLY*PLY+PLZ*PLZ+Q(JTY+7)**2)
 826 C
 827 C --- ADD PARTICLE TO THE STACK IF STACK NOT YET FULL ---
 828       IF (NGKINE .GE. MXGKIN) THEN
 829           WRITE(CHMAIL,1236) NTOT, L
 830  1236     FORMAT(' *** GHEISH: ',I9,' particle produced but only ',
 831      +           I9,' put on the GEANT stack!')
 832           CALL GMAIL(1,1)
 833           GO TO 9999
 834       ENDIF
 835       NGKINE=NGKINE+1
 836       GKIN(1,NGKINE)=PLX
 837       GKIN(2,NGKINE)=PLY
 838       GKIN(3,NGKINE)=PLZ
 839       GKIN(4,NGKINE)=ELT
 840       GKIN(5,NGKINE)=ITY
 841       TOFD(NGKINE)=EVE(INDEX+8)*0.5E-10
 842       GPOS(1,NGKINE) = VECT(1)
 843       GPOS(2,NGKINE) = VECT(2)
 844       GPOS(3,NGKINE) = VECT(3)
 845 C
 846       IF (NPRT(9)) PRINT 1006,ITY,NGKINE,L,(EVE(INDEX+J),J=1,12)
 847  1006 FORMAT(' *GHEISH* GEANT PART. ',I3,' ALSO PUT ONTO STACK AT',
 848      $ ' POS. ',I3/
 849      $ ' EVE(',I2,') = '/12(1H ,12X,G12.5/))
 850 C
 851  61   CONTINUE
 852 C
 853  9999 CONTINUE
 854 C --- LIMIT THE VALUE OF NGKINE IN CASE OF OVERFLOW ---
 855       NGKINE=MIN(NGKINE,MXGKIN)
 856       END