GEANT321/gheisha/caspim.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:21:00  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 CASPIM(K,INT,NFL)
  13 C
  14 C *** CASCADE OF PI- ***
  15 C *** NVE 04-MAY-1988 CERN GENEVA ***
  16 C
  17 C ORIGIN : H.FESEFELDT 13-SEP-1987
  18 C
  19 C PI-  UNDERGOES INTERACTION WITH NUCLEON WITHIN NUCLEUS.
  20 C CHECK IF ENERGETICALLY POSSIBLE TO PRODUCE PIONS/KAONS.
  21 C IF NOT ASSUME NUCLEAR EXCITATION OCCURS AND INPUT PARTICLE
  22 C IS DEGRADED IN ENERGY.    NO OTHER PARTICLES PRODUCED.
  23 C IF REACTION IS POSSIBLE FIND CORRECT NUMBER OF PIONS/PROTONS/
  24 C NEUTRONS PRODUCED USING AN INTERPOLATION TO MULTIPLICITY DATA.
  25 C REPLACE SOME PIONS OR PROTONS/NEUTRONS BY KAONS OR STRANGE BARYONS
  26 C ACCORDING TO AVERAGE MULTIPLICITY PER INELASTIC REACTIONS.
  27 C
  28 #include "geant321/mxgkgh.inc"
  29 #include "geant321/s_consts.inc"
  30 #include "geant321/s_curpar.inc"
  31 #include "geant321/s_result.inc"
  32 #include "geant321/s_prntfl.inc"
  33 #include "geant321/limits.inc"
  34 #include "geant321/s_kginit.inc"
  35 C
  36       REAL N
  37       DIMENSION PMUL(2,1200),ANORM(2,60),SUPP(10),CECH(10),B(2)
  38       DIMENSION RNDM(1)
  39       SAVE PMUL,ANORM
  40       DATA SUPP/0.,0.4,0.55,0.65,0.75,0.82,0.86,0.90,0.94,0.98/
  41       DATA CECH/1.,0.95,0.79,0.32,0.19,0.16,0.14,0.12,0.10,0.08/
  42       DATA B/0.7,0.7/,C/1.25/
  43 C
  44 C --- INITIALIZATION INDICATED BY KGINIT(16) ---
  45       IF (KGINIT(16) .NE. 0) GO TO 10
  46       KGINIT(16)=1
  47 C
  48 C --- INITIALIZE PMUL AND ANORM ARRAYS ---
  49       DO 9000 J=1,1200
  50       DO 9001 I=1,2
  51       PMUL(I,J)=0.0
  52       IF (J .LE. 60) ANORM(I,J)=0.0
  53  9001 CONTINUE
  54  9000 CONTINUE
  55 C
  56 C *** COMPUTATION OF NORMALIZATION CONSTANTS ***
  57 C
  58 C --- P TARGET ---
  59       L=0
  60       DO 1100 NP1=1,20
  61       NP=NP1-1
  62       NMM1=NP1-1
  63       IF (NMM1 .LE. 1) NMM1=1
  64       NPP1=NP1+1
  65 C
  66       DO 1101 NM1=NMM1,NPP1
  67       NM=NM1-1
  68 C
  69       DO 1102 NZ1=1,20
  70       NZ=NZ1-1
  71       L=L+1
  72       IF (L .GT. 1200) GOTO 1199
  73       NT=NP+NM+NZ
  74       IF (NT .LE. 0) GO TO 1102
  75       IF (NT .GT. 60) GO TO 1102
  76       PMUL(1,L)=PMLTPC(NP,NM,NZ,NT,B(1),C)
  77       ANORM(1,NT)=ANORM(1,NT)+PMUL(1,L)
  78  1102 CONTINUE
  79 C
  80  1101 CONTINUE
  81 C
  82  1100 CONTINUE
  83 C
  84  1199 CONTINUE
  85 C
  86 C --- N TARGET ---
  87       L=0
  88       DO 1200 NP1=1,20
  89       NP=NP1-1
  90       NPP1=NP1+2
  91 C
  92       DO 1201 NM1=NP1,NPP1
  93       NM=NM1-1
  94 C
  95       DO 1202 NZ1=1,20
  96       NZ=NZ1-1
  97       L=L+1
  98       IF (L .GT. 1200) GO TO 1299
  99       NT=NP+NM+NZ
 100       IF (NT .LE. 0) GO TO 1202
 101       IF (NT .GT. 60) GO TO 1202
 102       PMUL(2,L)=PMLTPC(NP,NM,NZ,NT,B(2),C)
 103       ANORM(2,NT)=ANORM(2,NT)+PMUL(2,L)
 104  1202 CONTINUE
 105 C
 106  1201 CONTINUE
 107 C
 108  1200 CONTINUE
 109 C
 110  1299 CONTINUE
 111 C
 112       DO 3 I=1,60
 113       IF (ANORM(1,I) .GT. 0.0) ANORM(1,I)=1.0/ANORM(1,I)
 114       IF (ANORM(2,I) .GT. 0.0) ANORM(2,I)=1.0/ANORM(2,I)
 115     3 CONTINUE
 116 C
 117       IF (.NOT. NPRT(10)) GO TO 10
 118       WRITE(NEWBCD,2001)
 119       DO 4 NFL=1,2
 120       WRITE(NEWBCD,2002) NFL
 121       WRITE(NEWBCD,2003) (ANORM(NFL,I),I=1,60)
 122       WRITE(NEWBCD,2003) (PMUL(NFL,I),I=1,1200)
 123     4 CONTINUE
 124 C
 125 C --- CHOOSE PROTON OR NEUTRON AS TARGET ---
 126  10   CONTINUE
 127       NFL=2
 128       CALL GRNDM(RNDM,1)
 129       IF (RNDM(1) .LT. ZNO2/ATNO2) NFL=1
 130       TARMAS=RMASS(14)
 131       IF (NFL .EQ. 2) TARMAS=RMASS(16)
 132       S=AMASQ+TARMAS**2+2.0*TARMAS*EN
 133       RS=SQRT(S)
 134       ENP(8)=AMASQ+TARMAS**2+2.0*TARMAS*ENP(6)
 135       ENP(9)=SQRT(ENP(8))
 136       EAB=RS-TARMAS-RMASS(9)
 137 C
 138 C --- ELASTIC SCATTERING ---
 139       NP=0
 140       NM=0
 141       NZ=0
 142       N=0.0
 143       IPA(1)=9
 144       IPA(2)=14
 145       IF (NFL .EQ. 2) IPA(2)=16
 146       IF (INT .EQ. 2) GOTO 20
 147       GOTO 100
 148 C
 149 C --- CHECK IF ENERGETICALLY POSSIBLE TO PRODUCE ONE EXTRA PION IN REACT.
 150  20   CONTINUE
 151       IF (EAB .LE. RMASS(9)) GO TO 55
 152 C
 153 C --- SUPPRESSION OF HIGH MULTIPLICITY EVENTS AT LOW MOMENTUM ---
 154       IEAB=IFIX(EAB*5.0)+1
 155       IF (IEAB .GT. 10) GO TO 22
 156       CALL GRNDM(RNDM,1)
 157       IF (RNDM(1) .LT. SUPP(IEAB)) GO TO 22
 158 C
 159 C --- CHARGE EXCHANGE REACTION (IS INCLUDED IN INELASTIC CROSS SECTION)
 160       IPLAB=IFIX(P*5.0)+1
 161       IF (IPLAB .GT. 10) IPLAB=10
 162       CALL GRNDM(RNDM,1)
 163       IF (RNDM(1) .GT. CECH(IPLAB)) GO TO 23
 164 C
 165       IF (NFL .EQ. 1) GOTO 24
 166 C
 167 C --- N TARGET ---
 168       INT=1
 169       IPA(1)=9
 170       IPA(2)=16
 171       GO TO 100
 172 C
 173 C --- P TARGET ---
 174  24   CONTINUE
 175       IPA(1)=8
 176       IPA(2)=16
 177       GO TO 100
 178 C
 179  23   CONTINUE
 180       N=1.0
 181 C
 182       IF (NFL .EQ. 1) GO TO 26
 183 C
 184 C --- N TARGET ---
 185       DUM=-(1+B(2))**2/(2.0*C**2)
 186       IF (DUM .LT. EXPXL) DUM=EXPXL
 187       IF (DUM .GT. EXPXU) DUM=EXPXU
 188       W0=EXP(DUM)
 189       DUM=-(-1+B(2))**2/(2.0*C**2)
 190       IF (DUM .LT. EXPXL) DUM=EXPXL
 191       IF (DUM .GT. EXPXU) DUM=EXPXU
 192       WM=EXP(DUM)
 193       CALL GRNDM(RNDM,1)
 194       RAN=RNDM(1)
 195       NP=0
 196       NM=0
 197       NZ=1
 198       IF (RAN .LT. W0/(W0+WM)) GO TO 50
 199       NP=0
 200       NM=1
 201       NZ=0
 202       GO TO 50
 203 C
 204 C --- P TARGET ---
 205  26   CONTINUE
 206       DUM=-(1+B(1))**2/(2.0*C**2)
 207       IF (DUM .LT. EXPXL) DUM=EXPXL
 208       IF (DUM .GT. EXPXU) DUM=EXPXU
 209       W0=EXP(DUM)
 210       WP=EXP(DUM)
 211       DUM=-(-1+B(1))**2/(2.0*C**2)
 212       IF (DUM .LT. EXPXL) DUM=EXPXL
 213       IF (DUM .GT. EXPXU) DUM=EXPXU
 214       WM=EXP(DUM)
 215       WP=WP*10.
 216       WT=W0+WP+WM
 217       WP=W0+WP
 218       CALL GRNDM(RNDM,1)
 219       RAN=RNDM(1)
 220       NP=0
 221       NM=0
 222       NZ=1
 223       IF (RAN .LT. W0/WT) GO TO 50
 224       NP=1
 225       NM=0
 226       NZ=0
 227       IF (RAN .LT. WP/WT) GO TO 50
 228       NP=0
 229       NM=1
 230       NZ=0
 231       GOTO 50
 232 C
 233  22   CONTINUE
 234       ALEAB=LOG(EAB)
 235 C
 236 C --- NO. OF TOTAL PARTICLES VS SQRT(S)-2*MP ---
 237       N=3.62567+0.665843*ALEAB+0.336514*ALEAB*ALEAB
 238      $ +0.117712*ALEAB*ALEAB*ALEAB+0.0136912*ALEAB*ALEAB*ALEAB*ALEAB
 239       N=N-2.0
 240 C
 241 C --- NORMALIZATION CONSTANT FOR KNO-DISTRIBUTION ---
 242       ANPN=0.0
 243       DO 21 NT=1,60
 244       TEST=-(PI/4.0)*(NT/N)**2
 245       IF (TEST .LT. EXPXL) TEST=EXPXL
 246       IF (TEST .GT. EXPXU) TEST=EXPXU
 247       DUM1=PI*NT/(2.0*N*N)
 248       DUM2=ABS(DUM1)
 249       DUM3=EXP(TEST)
 250       ADDNVE=0.0
 251       IF (DUM2 .GE. 1.0) ADDNVE=DUM1*DUM3
 252       IF ((DUM2 .LT. 1.0) .AND. (DUM3 .GE. 1.0E-10)) ADDNVE=DUM1*DUM3
 253       ANPN=ANPN+ADDNVE
 254    21 CONTINUE
 255       ANPN=1.0/ANPN
 256 C
 257       CALL GRNDM(RNDM,1)
 258       RAN=RNDM(1)
 259       EXCS=0.0
 260       IF (NFL .EQ. 2) GO TO 40
 261 C
 262 C --- P TARGET ---
 263       L=0
 264       DO 310 NP1=1,20
 265       NP=NP1-1
 266       NMM1=NP1-1
 267       IF (NMM1 .LE. 1) NMM1=1
 268       NPP1=NP1+1
 269 C
 270       DO 311 NM1=NMM1,NPP1
 271       NM=NM1-1
 272 C
 273       DO 312 NZ1=1,20
 274       NZ=NZ1-1
 275       L=L+1
 276       IF (L .GT. 1200) GO TO 80
 277       NT=NP+NM+NZ
 278       IF (NT .LE. 0) GO TO 312
 279       IF (NT .GT. 60) GO TO 312
 280       TEST=-(PI/4.0)*(NT/N)**2
 281       IF (TEST .LT. EXPXL) TEST=EXPXL
 282       IF (TEST .GT. EXPXU) TEST=EXPXU
 283       DUM1=ANPN*PI*NT*PMUL(1,L)*ANORM(1,NT)/(2.0*N*N)
 284       DUM2=ABS(DUM1)
 285       DUM3=EXP(TEST)
 286       ADDNVE=0.0
 287       IF (DUM2 .GE. 1.0) ADDNVE=DUM1*DUM3
 288       IF ((DUM2 .LT. 1.0) .AND. (DUM3 .GE. 1.0E-10)) ADDNVE=DUM1*DUM3
 289       EXCS=EXCS+ADDNVE
 290       IF (RAN .LT. EXCS) GOTO 50
 291  312  CONTINUE
 292 C
 293  311  CONTINUE
 294 C
 295  310  CONTINUE
 296       GOTO 80
 297 C
 298 C --- N TARGET ---
 299  40   CONTINUE
 300       L=0
 301       DO 410 NP1=1,20
 302       NP=NP1-1
 303       NPP1=NP1+2
 304 C
 305       DO 411 NM1=NP1,NPP1
 306       NM=NM1-1
 307 C
 308       DO 412 NZ1=1,20
 309       NZ=NZ1-1
 310       L=L+1
 311       IF (L .GT. 1200) GO TO 80
 312       NT=NP+NM+NZ
 313       IF (NT .LE. 0) GO TO 412
 314       IF (NT .GT. 60) GO TO 412
 315       TEST=-(PI/4.0)*(NT/N)**2
 316       IF (TEST .LT. EXPXL) TEST=EXPXL
 317       IF (TEST .GT. EXPXU) TEST=EXPXU
 318       DUM1=ANPN*PI*NT*PMUL(2,L)*ANORM(2,NT)/(2.0*N*N)
 319       DUM2=ABS(DUM1)
 320       DUM3=EXP(TEST)
 321       ADDNVE=0.0
 322       IF (DUM2 .GE. 1.0) ADDNVE=DUM1*DUM3
 323       IF ((DUM2 .LT. 1.0) .AND. (DUM3 .GE. 1.0E-10)) ADDNVE=DUM1*DUM3
 324       EXCS=EXCS+ADDNVE
 325       IF (RAN .LT. EXCS) GOTO 50
 326  412  CONTINUE
 327 C
 328  411  CONTINUE
 329 C
 330  410  CONTINUE
 331       GO TO 80
 332 C
 333  50   CONTINUE
 334       IF (NFL .EQ. 2) GO TO 65
 335 C
 336 C --- P TARGET ---
 337       IF (NP .EQ. NM) GO TO 61
 338       IF (NP .EQ. 1+NM) GO TO 63
 339       IPA(1)=8
 340       IPA(2)=14
 341       GO TO 100
 342 C
 343  61   CONTINUE
 344       CALL GRNDM(RNDM,1)
 345       IF (RNDM(1) .LT. 0.75) GO TO 62
 346       IPA(1)=8
 347       IPA(2)=16
 348       GO TO 100
 349 C
 350  62   CONTINUE
 351       IPA(1)=9
 352       IPA(2)=14
 353       GO TO 100
 354 C
 355  63   CONTINUE
 356       IPA(1)=9
 357       IPA(2)=16
 358       GO TO 100
 359 C
 360 C --- N TARGET ---
 361  65   CONTINUE
 362       IF (NP .EQ. -1+NM) GO TO 66
 363       IF (NP .EQ. NM) GO TO 68
 364       IPA(1)=8
 365       IPA(2)=16
 366       GO TO 100
 367 C
 368  66   CONTINUE
 369       CALL GRNDM(RNDM,1)
 370       IF (RNDM(1) .LT. 0.50) GO TO 67
 371       IPA(1)=8
 372       IPA(2)=16
 373       GO TO 100
 374 C
 375  67   CONTINUE
 376       IPA(1)=9
 377       IPA(2)=14
 378       GO TO 100
 379 C
 380  68   CONTINUE
 381       IPA(1)=9
 382       IPA(2)=16
 383       GO TO 100
 384 C
 385  70   CONTINUE
 386       IF (NPRT(4)) WRITE(NEWBCD,1003) EAB,N,NFL,NP,NM,NZ
 387       CALL STPAIR
 388       IF (INT .EQ. 1) CALL TWOB(9,NFL,N)
 389       IF (INT .EQ. 2) CALL GENXPT(9,NFL,N)
 390       GO TO 9999
 391 C
 392 C --- ENERGETICALLY NOT POSSIBLE TO PRODUCE CASCADE-PARTICLES ---
 393 C --- CONTINUE WITH QUASI-ELASTIC SCATTERING ---
 394  55   CONTINUE
 395       IF (NPRT(4)) WRITE(NEWBCD,1001)
 396       GO TO 53
 397 C
 398 C --- EXCLUSIVE REACTION NOT FOUND ---
 399  80   CONTINUE
 400       IF (NPRT(4)) WRITE(NEWBCD,1004) RS,N
 401 C
 402  53   CONTINUE
 403       INT=1
 404       NP=0
 405       NM=0
 406       NZ=0
 407       N=0.0
 408       IPA(1)=9
 409       IPA(2)=14
 410       IF (NFL .EQ. 2) IPA(2)=16
 411 C
 412  100  CONTINUE
 413       DO 101 I=3,60
 414       IPA(I)=0
 415  101  CONTINUE
 416       IF (INT .LE. 0) GO TO 131
 417 C
 418  120  CONTINUE
 419       NT=2
 420       IF (NP .EQ. 0) GO TO 122
 421       DO 121 I=1,NP
 422       NT=NT+1
 423       IPA(NT)=7
 424  121  CONTINUE
 425 C
 426  122  CONTINUE
 427       IF (NM .EQ. 0) GO TO 124
 428       DO 123 I=1,NM
 429       NT=NT+1
 430       IPA(NT)=9
 431  123  CONTINUE
 432 C
 433  124  CONTINUE
 434       IF (NZ .EQ. 0) GO TO 130
 435       DO 125 I=1,NZ
 436       NT=NT+1
 437       IPA(NT)=8
 438  125  CONTINUE
 439 C
 440  130  CONTINUE
 441       IF (NPRT(4)) WRITE(NEWBCD,2004) NT,(IPA(I),I=1,20)
 442       IF (IPA(1) .EQ. 7) NP=NP+1
 443       IF (IPA(1) .EQ. 8) NZ=NZ+1
 444       IF (IPA(1) .EQ. 9) NM=NM+1
 445       GO TO 70
 446 C
 447  131  CONTINUE
 448       IF (NPRT(4)) WRITE(NEWBCD,2005)
 449 C
 450 1001  FORMAT('0*CASPIM* CASCADE ENERGETICALLY NOT POSSIBLE',
 451      $ ' CONTINUE WITH QUASI-ELASTIC SCATTERING')
 452 1003  FORMAT(' *CASPIM* PION- -INDUCED CASCADE, AVAIL. ENERGY',2X,F8.4,
 453      $ 2X,'<NTOT>',2X,F8.4,2X,'FROM',4(2X,I3),2X,'PARTICLES')
 454 1004  FORMAT(' *CASPIM* PION- -INDUCED CASCADE, EXCLUSIVE REACTION',
 455      $ ' NOT FOUND TRY ELASTIC SCATTERING  AVAIL. ENERGY',2X,F8.4,2X,
 456      * '<NTOT>',2X,F8.4)
 457 2001  FORMAT('0*CASPIM* TABLES FOR MULTIPLICITY DATA PION- INDUCED',
 458      $ 'REACTION FOR DEFINITION OF NUMBERS SEE FORTRAN CODING')
 459 2002  FORMAT(' *CASPIM* TARGET PARTICLE FLAG',2X,I5)
 460 2003  FORMAT(1H ,10E12.4)
 461 2004  FORMAT(' *CASPIM* ',I3,2X,'PARTICLES, MASS INDEX ARRAY',2X,20I4)
 462 2005  FORMAT(' *CASPIM* NO PARTICLES PRODUCED')
 463 C
 464  9999 CONTINUE
 465       END