GEANT321/gheisha/casx0.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:21:05  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/02 29/03/94  15.41.40  by  S.Giani
  11 *-- Author :
  12       SUBROUTINE CASX0(K,INT,NFL)
  13 C
  14 C *** CASCADE OF XI0 ***
  15 C *** NVE 20-JAN-1989 CERN GENEVA ***
  16 C
  17 C XI0  UNDERGOES INTERACTION WITH NUCLEON WITHIN NUCLEUS.
  18 C CHECK IF ENERGETICALLY POSSIBLE TO PRODUCE PIONS/KAONS.
  19 C IF NOT, ASSUME NUCLEAR EXCITATION OCCURS, DEGRADE INPUT PARTICLE
  20 C IN ENERGY AND NO OTHER PARTICLES ARE PRODUCED.
  21 C IF REACTION IS POSSIBLE FIND CORRECT NUMBER OF PIONS/PROTONS/
  22 C NEUTRONS PRODUCED USING AN INTERPOLATION TO MULTIPLICITY DATA.
  23 C REPLACE SOME PIONS OR PROTONS/NEUTRONS BY KAONS OR STRANGE BARYONS
  24 C ACCORDING TO AVERAGE MULTIPLICITY PER INELASTIC REACTIONS.
  25 C
  26 #include "geant321/mxgkgh.inc"
  27 #include "geant321/s_consts.inc"
  28 #include "geant321/s_curpar.inc"
  29 #include "geant321/s_result.inc"
  30 #include "geant321/s_prntfl.inc"
  31 #include "geant321/s_kginit.inc"
  32 #include "geant321/limits.inc"
  33 C
  34       REAL N
  35       DIMENSION PMUL(2,1200),ANORM(2,60),CECH(10),IIPA(12,2),B(2)
  36       DIMENSION RNDM(1)
  37       SAVE PMUL,ANORM
  38       DATA CECH/0.50,0.45,0.40,0.35,0.30,0.25,0.06,0.04,0.005,0./
  39 C --- ARRAY IIPA DENOTES THE STRANGENESS AND CHARGE EXCHAGE REACTIONS ---
  40 C XI0 P --> S+ S0,  XI0 P --> S0 S+
  41 C XI0 P --> S+ L0,  XI0 P --> L0 S+
  42 C XI0 P --> P XI0
  43 C XI0 N --> S0 S0
  44 C XI0 N --> L0 L0
  45 C XI0 N --> XI- P,  XI0 N --> P XI-
  46 C XI0 N --> S+ S-,  XI0 N --> S- S+
  47 C XI0 N --> N XI0
  48       DATA IIPA/20,21,20,18,14, 21,18,27,14,20,22,16,
  49      *          21,20,18,20,26, 21,18,14,27,22,20,26/
  50       DATA B/0.7,0.7/,C/1.25/
  51 C
  52 C --- INITIALIZATION INDICATED BY KGINIT(20) ---
  53       IF (KGINIT(20) .NE. 0) GO TO 10
  54       KGINIT(20)=1
  55 C
  56 C --- INITIALIZE PMUL AND ANORM ARRAYS ---
  57       DO 9000 J=1,1200
  58       DO 9001 I=1,2
  59       PMUL(I,J)=0.0
  60       IF (J .LE. 60) ANORM(I,J)=0.0
  61  9001 CONTINUE
  62  9000 CONTINUE
  63 C
  64 C *** COMPUTE NORMALIZATION CONSTANTS ***
  65 C
  66 C --- FOR P TARGET ---
  67       L=0
  68       DO 1 NP1=1,20
  69       NP=NP1-1
  70       NMM1=NP1-2
  71       IF (NMM1 .LE. 0) NMM1=1
  72       NPP1=NP1+1
  73       DO 1 NM1=NMM1,NPP1
  74       NM=NM1-1
  75       DO 1 NZ1=1,20
  76       NZ=NZ1-1
  77       L=L+1
  78       IF (L .GT. 1200) GO TO 1
  79       NT=NP+NM+NZ
  80       IF ((NT .LE. 0) .OR. (NT .GT. 60)) GO TO 1
  81       PMUL(1,L)=PMLTPC(NP,NM,NZ,NT,B(1),C)
  82       ANORM(1,NT)=ANORM(1,NT)+PMUL(1,L)
  83  1    CONTINUE
  84 C --- FOR N TARGET ---
  85       L=0
  86       DO 2 NP1=1,20
  87       NP=NP1-1
  88       NMM1=NP1-1
  89       IF (NMM1 .LE. 0) NMM1=1
  90       NPP1=NP1+2
  91       DO 2 NM1=NMM1,NPP1
  92       NM=NM1-1
  93       DO 2 NZ1=1,20
  94       NZ=NZ1-1
  95       L=L+1
  96       IF (L .GT. 1200) GO TO 2
  97       NT=NP+NM+NZ
  98       IF ((NT .LE. 0) .OR. (NT .GT. 60)) GO TO 2
  99       PMUL(2,L)=PMLTPC(NP,NM,NZ,NT,B(2),C)
 100       ANORM(2,NT)=ANORM(2,NT)+PMUL(2,L)
 101  2    CONTINUE
 102 C
 103       DO 3 I=1,60
 104       IF (ANORM(1,I) .GT. 0.) ANORM(1,I)=1./ANORM(1,I)
 105       IF (ANORM(2,I) .GT. 0.) ANORM(2,I)=1./ANORM(2,I)
 106  3    CONTINUE
 107 C
 108       IF (.NOT. NPRT(10)) GO TO 10
 109 C
 110       WRITE(NEWBCD,2001)
 111  2001 FORMAT('0*CASX0* TABLES FOR MULT. DATA XI0 INDUCED REACTION',
 112      $ ' FOR DEFINITION OF NUMBERS SEE FORTRAN CODING')
 113       DO 4 NFL=1,2
 114       WRITE(NEWBCD,2002) NFL
 115  2002 FORMAT(' *CASX0* TARGET PARTICLE FLAG',2X,I5)
 116       WRITE(NEWBCD,2003) (ANORM(NFL,I),I=1,60)
 117       WRITE(NEWBCD,2003) (PMUL(NFL,I),I=1,1200)
 118  2003 FORMAT(1H ,10E12.4)
 119  4    CONTINUE
 120 C
 121 C --- SELECT TARGET NUCLEON ---
 122  10   CONTINUE
 123       NFL=2
 124       CALL GRNDM(RNDM,1)
 125       IF (RNDM(1) .LT. (ZNO2/ATNO2)) NFL=1
 126       TARMAS=RMASS(14)
 127       IF (NFL .EQ. 2) TARMAS=RMASS(16)
 128       S=AMASQ+TARMAS**2+2.0*TARMAS*EN
 129       RS=SQRT(S)
 130       ENP(8)=AMASQ+TARMAS**2+2.0*TARMAS*ENP(6)
 131       ENP(9)=SQRT(ENP(8))
 132       EAB=RS-TARMAS-RMASS(26)
 133 C
 134 C --- RESET STRANGENESS FIXING FLAG ---
 135       NVEFIX=0
 136 C
 137 C *** ELASTIC SCATTERING ***
 138       NP=0
 139       NM=0
 140       NZ=0
 141       N=0.
 142       IPA(1)=26
 143       IPA(2)=14
 144       IF (NFL .EQ. 2) IPA(2)=16
 145 C
 146       IF (INT .EQ. 2) GO TO 20
 147 C
 148 C *** INTRODUCE CHARGE AND STRANGENESS EXCHANGE REACTIONS ***
 149       IPLAB=IFIX(P*2.5)+1
 150       IF (IPLAB .GT. 10) IPLAB=10
 151       CALL GRNDM(RNDM,1)
 152       IF (RNDM(1) .GT. (CECH(IPLAB)/ATNO2**0.42)) GO TO 120
 153       CALL GRNDM(RNDM,1)
 154       RAN=RNDM(1)
 155       IRN=IFIX(RAN*5.)+1
 156       IF (NFL .EQ. 2) IRN=5+IFIX(RAN*7.)+1
 157       IF (NFL .EQ. 1) IRN=MAX(IRN,5)
 158       IF (NFL .EQ. 2) IRN=MAX(IRN,12)
 159       IPA(1)=IIPA(IRN,1)
 160       IPA(2)=IIPA(IRN,2)
 161       GO TO 120
 162 C
 163 C --- CHECK IF ENERGETICALLY POSSIBLE TO PRODUCE ONE EXTRA PION ---
 164  20   CONTINUE
 165       IF (EAB .LE. RMASS(7)) GO TO 55
 166 C
 167 C --- NO. OF TOTAL PARTICLES VS SQRT(S)-MP-MSM ---
 168       ALEAB=LOG(EAB)
 169       N=3.62567+0.665843*ALEAB+0.336514*ALEAB*ALEAB
 170      * +0.117712*ALEAB*ALEAB*ALEAB+0.0136912*ALEAB*ALEAB*ALEAB*ALEAB
 171       N=N-2.
 172 C
 173 C --- NORMALIZATION CONSTANT FOR  KNO-DISTRIBUTION ---
 174       ANPN=0.
 175       DO 21 NT=1,60
 176       TEST=-(PI/4.0)*(NT/N)**2
 177       IF (TEST .LT. EXPXL) TEST=EXPXL
 178       IF (TEST .GT. EXPXU) TEST=EXPXU
 179       DUM1=PI*NT/(2.0*N*N)
 180       DUM2=ABS(DUM1)
 181       DUM3=EXP(TEST)
 182       ADDNVE=0.0
 183       IF (DUM2 .GE. 1.0) ADDNVE=DUM1*DUM3
 184       IF ((DUM2 .LT. 1.0) .AND. (DUM3 .GE. 1.0E-10)) ADDNVE=DUM1*DUM3
 185       ANPN=ANPN+ADDNVE
 186  21   CONTINUE
 187       ANPN=1./ANPN
 188 C
 189 C --- CHECK FOR TARGET NUCLEON TYPE ---
 190       CALL GRNDM(RNDM,1)
 191       RAN=RNDM(1)
 192       EXCS=0.
 193       GO TO (30,40),NFL
 194 C
 195 C --- PROTON TARGET ---
 196  30   CONTINUE
 197       L=0
 198       DO 31 NP1=1,20
 199       NP=NP1-1
 200       NMM1=NP1-2
 201       IF (NMM1 .LE. 0) NMM1=1
 202       NPP1=NP1+1
 203       DO 31 NM1=NMM1,NPP1
 204       NM=NM1-1
 205       DO 31 NZ1=1,20
 206       NZ=NZ1-1
 207       L=L+1
 208       IF (L .GT. 1200) GO TO 31
 209       NT=NP+NM+NZ
 210       IF ((NT .LE. 0) .OR. (NT .GT. 60)) GO TO 31
 211       TEST=-(PI/4.0)*(NT/N)**2
 212       IF (TEST .LT. EXPXL) TEST=EXPXL
 213       IF (TEST .GT. EXPXU) TEST=EXPXU
 214       DUM1=ANPN*PI*NT*PMUL(1,L)*ANORM(1,NT)/(2.0*N*N)
 215       DUM2=ABS(DUM1)
 216       DUM3=EXP(TEST)
 217       ADDNVE=0.0
 218       IF (DUM2 .GE. 1.0) ADDNVE=DUM1*DUM3
 219       IF ((DUM2 .LT. 1.0) .AND. (DUM3 .GE. 1.0E-10)) ADDNVE=DUM1*DUM3
 220       EXCS=EXCS+ADDNVE
 221       IF (RAN .LT. EXCS) GO TO 100
 222    31 CONTINUE
 223       GO TO 80
 224 C
 225 C --- NEUTRON TARGET ---
 226  40   CONTINUE
 227       L=0
 228       DO 41 NP1=1,20
 229       NP=NP1-1
 230       NMM1=NP1-1
 231       IF (NMM1 .LE. 0) NMM1=1
 232       NPP1=NP1+2
 233       DO 41 NM1=NMM1,NPP1
 234       NM=NM1-1
 235       DO 41 NZ1=1,20
 236       NZ=NZ1-1
 237       L=L+1
 238       IF (L .GT. 1200) GO TO 41
 239       NT=NP+NM+NZ
 240       IF ((NT .LE. 0) .OR. (NT .GT. 60)) GO TO 41
 241       TEST=-(PI/4.0)*(NT/N)**2
 242       IF (TEST .LT. EXPXL) TEST=EXPXL
 243       IF (TEST .GT. EXPXU) TEST=EXPXU
 244       DUM1=ANPN*PI*NT*PMUL(2,L)*ANORM(2,NT)/(2.0*N*N)
 245       DUM2=ABS(DUM1)
 246       DUM3=EXP(TEST)
 247       ADDNVE=0.0
 248       IF (DUM2 .GE. 1.0) ADDNVE=DUM1*DUM3
 249       IF ((DUM2 .LT. 1.0) .AND. (DUM3 .GE. 1.0E-10)) ADDNVE=DUM1*DUM3
 250       EXCS=EXCS+ADDNVE
 251       IF (RAN .LT. EXCS) GO TO 100
 252    41 CONTINUE
 253       GO TO 80
 254 C
 255  50   CONTINUE
 256       IF (NPRT(4)) WRITE(NEWBCD,1003) EAB,N,NFL,NP,NM,NZ
 257  1003 FORMAT(' *CASX0* XI0 -INDUCED CASCADE,',
 258      $ ' AVAIL. ENERGY',2X,F8.4,
 259      $ 2X,'<NTOT>',2X,F8.4,2X,'FROM',4(2X,I3),2X,'PARTICLES')
 260       IF (INT .EQ. 1) CALL TWOB(27,NFL,N)
 261       IF (INT .EQ. 2) CALL GENXPT(27,NFL,N)
 262       GO TO 9999
 263 C
 264 C *** ENERGETICALLY NOT POSSIBLE TO PRODUCE ONE EXTRA PION ***
 265  55   CONTINUE
 266       IF (NPRT(4)) WRITE(NEWBCD,1001)
 267  1001 FORMAT('0*CASX0* CASCADE ENERGETICALLY NOT POSSIBLE',
 268      $ ' CONTINUE WITH QUASI-ELASTIC SCATTERING')
 269       GO TO 53
 270 C
 271 C *** EXCLUSIVE REACTION NOT FOUND ***
 272  80   CONTINUE
 273       IF (NPRT(4)) WRITE(NEWBCD,1004) RS,N
 274  1004 FORMAT(' *CASX0* XI0 -INDUCED CASCADE,',
 275      $ ' EXCLUSIVE REACTION NOT FOUND',
 276      $ ' TRY ELASTIC SCATTERING  AVAIL. ENERGY',2X,F8.4,2X,
 277      $ '<NTOT>',2X,F8.4)
 278 C
 279  53   CONTINUE
 280       INT=1
 281       NP=0
 282       NM=0
 283       NZ=0
 284       IPA(1)=26
 285       IPA(2)=14
 286       IF (NFL .EQ. 2) IPA(2)=16
 287       GO TO 120
 288 C
 289 C *** INELASTIC INTERACTION HAS OCCURRED ***
 290 C *** NUMBER OF SECONDARY MESONS DETERMINED BY KNO DISTRIBUTION ***
 291  100  CONTINUE
 292       DO 101 I=1,60
 293       IPA(I)=0
 294  101  CONTINUE
 295 C
 296       IF (INT .LE. 0) GO TO 131
 297 C
 298 C --- TAKE TARGET NUCLEON TYPE INTO ACCOUNT ---
 299       GO TO (102,112),NFL
 300 C
 301 C --- PROTON TARGET ---
 302  102  CONTINUE
 303 C --- CHECK FOR TOTAL CHARGE OF FINAL STATE MESONS TO DETERMINE ---
 304 C --- THE KIND OF BARYONS TO BE PRODUCED TAKING INTO ACCOUNT    ---
 305 C --- CHARGE AND STRANGENESS CONSERVATION                       ---
 306       NCHT=NP-NM
 307       IF (NCHT .LT. 1) GO TO 103
 308       IF (NCHT .EQ. 1) GO TO 104
 309       IF (NCHT .GT. 1) GO TO 105
 310 C
 311  103  CONTINUE
 312 C --- XI0 P ---
 313       IPA(1)=26
 314       IPA(2)=14
 315       IF (NCHT .EQ. 0) GO TO 120
 316 C --- CHARGE MISMATCH ==> TAKE A S+ AND CORRECT THE STRANGENESS ---
 317 C --- BY REPLACING A PI- BY K- ---
 318 C --- S+ P ---
 319       IPA(1)=20
 320       IPA(2)=14
 321       NVEFIX=1
 322       GO TO 120
 323 C
 324  104  CONTINUE
 325 C --- XI0 N ---
 326       IPA(1)=26
 327       IPA(2)=16
 328       CALL GRNDM(RNDM,1)
 329       IF (RNDM(1) .LT. 0.5) GO TO 120
 330 C --- XI- P ---
 331       IPA(1)=27
 332       IPA(2)=14
 333       GO TO 120
 334 C
 335  105  CONTINUE
 336 C --- XI- N ---
 337       IPA(1)=27
 338       IPA(2)=16
 339       GO TO 120
 340 C
 341 C --- NEUTRON TARGET ---
 342  112  CONTINUE
 343 C --- CHECK FOR TOTAL CHARGE OF FINAL STATE MESONS TO DETERMINE ---
 344 C --- THE KIND OF BARYONS TO BE PRODUCED TAKING INTO ACCOUNT    ---
 345 C --- CHARGE AND STRANGENESS CONSERVATION                       ---
 346       NCHT=NP-NM
 347       IF (NCHT .LT. 0) GO TO 113
 348       IF (NCHT .EQ. 0) GO TO 114
 349       IF (NCHT .GT. 0) GO TO 115
 350 C
 351  113  CONTINUE
 352 C --- XI0 P ---
 353       IPA(1)=26
 354       IPA(2)=14
 355       IF (NCHT .EQ. -1) GO TO 120
 356 C --- CHARGE MISMATCH ==> TAKE A S+ AND CORRECT THE STRANGENESS ---
 357 C --- BY REPLACING A PI- BY K- ---
 358 C --- S+ P ---
 359       IPA(1)=20
 360       IPA(2)=14
 361       NVEFIX=1
 362       GO TO 120
 363 C
 364  114  CONTINUE
 365 C --- XI0 N ---
 366       IPA(1)=26
 367       IPA(2)=16
 368       CALL GRNDM(RNDM,1)
 369       IF (RNDM(1) .LT. 0.5) GO TO 120
 370 C --- XI- P ---
 371       IPA(1)=27
 372       IPA(2)=14
 373       GO TO 120
 374 C
 375  115  CONTINUE
 376 C --- XI- N ---
 377       IPA(1)=27
 378       IPA(2)=16
 379 C
 380 C --- TAKE PIONS FOR ALL SECONDARY MESONS ---
 381 C --- REPLACE PI BY K IN CASE OF STRANGENESS TO BE FIXED ---
 382  120  CONTINUE
 383       NT=2
 384 C
 385       IF (NP .EQ. 0) GO TO 122
 386 C
 387 C --- PI+ ---
 388       DO 121 I=1,NP
 389       NT=NT+1
 390       IPA(NT)=7
 391  121  CONTINUE
 392 C
 393  122  CONTINUE
 394       IF (NM .EQ. 0) GO TO 124
 395 C
 396 C --- PI- ---
 397       DO 123 I=1,NM
 398       NT=NT+1
 399       IPA(NT)=9
 400       IF (NVEFIX .GE. 1) IPA(NT)=13
 401       IF (NPRT(4) .AND. (NVEFIX .GE. 1)) PRINT 3000
 402  3000 FORMAT(' *CASX0* K- INTRODUCED')
 403       NVEFIX=NVEFIX-1
 404  123  CONTINUE
 405 C
 406  124  CONTINUE
 407       IF (NZ .EQ. 0) GO TO 130
 408 C
 409 C --- PI0 ---
 410       DO 125 I=1,NZ
 411       NT=NT+1
 412       IPA(NT)=8
 413  125  CONTINUE
 414 C
 415 C --- ALL SECONDARY PARTICLES HAVE BEEN DEFINED ---
 416 C --- NOW GO FOR MOMENTA AND X VALUES ---
 417  130  CONTINUE
 418       IF (NPRT(4)) WRITE(NEWBCD,2004) NT,(IPA(I),I=1,60)
 419  2004 FORMAT(' *CASX0* ',I3,' PARTICLES PRODUCED. MASS INDEX ARRAY : '/
 420      $ 3(1H ,20(I3,1X)/))
 421       GO TO 50
 422 C
 423  131  CONTINUE
 424       IF (NPRT(4)) WRITE(NEWBCD,2005)
 425  2005 FORMAT(' *CASX0* NO PARTICLES PRODUCED')
 426 C
 427  9999 CONTINUE
 428       END