| 1 | SUBROUTINE NIHILA
|
|---|
| 2 |
|
|---|
| 3 | C-----------------------------------------------------------------------
|
|---|
| 4 | C (AN)NIHILA(TION)
|
|---|
| 5 | C
|
|---|
| 6 | C TREATES ANNIHILATION OF ANTINUCLEONS WITH FREE NUCLEONS
|
|---|
| 7 | C MOMENTA CONSERVED IN ALL DIRECTIONS
|
|---|
| 8 | C ENERGY CONSERVED BY MULTPLICATION OF ALL MOMENTA WITH A CORRECTION
|
|---|
| 9 | C FACTOR, CONSERVING MOMENTUM BALANCE
|
|---|
| 10 | C THIS SUBROUTINE IS CALLED FROM BOX60, BOX61, BOX62, AND BOX63
|
|---|
| 11 | C
|
|---|
| 12 | C DESIGN : D. HECK IK3 FZK KARLSRUHE
|
|---|
| 13 | C-----------------------------------------------------------------------
|
|---|
| 14 |
|
|---|
| 15 | IMPLICIT DOUBLE PRECISION (A-H,O-Z)
|
|---|
| 16 | *KEEP,ANNI.
|
|---|
| 17 | COMMON /ANNI/ CAN,CANN
|
|---|
| 18 | DOUBLE PRECISION CAN(50),CANN(50)
|
|---|
| 19 | *KEEP,CONST.
|
|---|
| 20 | COMMON /CONST/ PI,PI2,OB3,TB3,ENEPER
|
|---|
| 21 | DOUBLE PRECISION PI,PI2,OB3,TB3,ENEPER
|
|---|
| 22 | *KEEP,MULT.
|
|---|
| 23 | COMMON /MULT/ EKINL,MSMM,MULTMA,MULTOT
|
|---|
| 24 | DOUBLE PRECISION EKINL
|
|---|
| 25 | INTEGER MSMM,MULTMA(37,13),MULTOT(37,13)
|
|---|
| 26 | *KEEP,PAM.
|
|---|
| 27 | COMMON /PAM/ PAMA,SIGNUM
|
|---|
| 28 | DOUBLE PRECISION PAMA(6000),SIGNUM(6000)
|
|---|
| 29 | *KEEP,PARPAR.
|
|---|
| 30 | COMMON /PARPAR/ CURPAR,SECPAR,PRMPAR,OUTPAR,C,
|
|---|
| 31 | * E00,E00PN,PTOT0,PTOT0N,THICKH,ITYPE,LEVL
|
|---|
| 32 | DOUBLE PRECISION CURPAR(14),SECPAR(14),PRMPAR(14),OUTPAR(14),
|
|---|
| 33 | * C(50),E00,E00PN,PTOT0,PTOT0N,THICKH
|
|---|
| 34 | INTEGER ITYPE,LEVL
|
|---|
| 35 | *KEEP,PARPAE.
|
|---|
| 36 | DOUBLE PRECISION GAMMA,COSTHE,PHI,H,T,X,Y,CHI,BETA,GCM,ECM
|
|---|
| 37 | EQUIVALENCE (CURPAR(2),GAMMA), (CURPAR(3),COSTHE),
|
|---|
| 38 | * (CURPAR(4), PHI ), (CURPAR(5), H ),
|
|---|
| 39 | * (CURPAR(6), T ), (CURPAR(7), X ),
|
|---|
| 40 | * (CURPAR(8), Y ), (CURPAR(9), CHI ),
|
|---|
| 41 | * (CURPAR(10),BETA), (CURPAR(11),GCM ),
|
|---|
| 42 | * (CURPAR(12),ECM )
|
|---|
| 43 | *KEEP,RANDPA.
|
|---|
| 44 | COMMON /RANDPA/ FAC,U1,U2,RD,NSEQ,ISEED,KNOR
|
|---|
| 45 | DOUBLE PRECISION FAC,U1,U2
|
|---|
| 46 | REAL RD(3000)
|
|---|
| 47 | INTEGER ISEED(103,10),NSEQ
|
|---|
| 48 | LOGICAL KNOR
|
|---|
| 49 | *KEEP,RUNPAR.
|
|---|
| 50 | COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
|
|---|
| 51 | * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
|
|---|
| 52 | * MONIOU,MDEBUG,NUCNUC,
|
|---|
| 53 | * CETAPE,
|
|---|
| 54 | * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
|
|---|
| 55 | * N1STTR,MDBASE,
|
|---|
| 56 | * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
|
|---|
| 57 | * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
|
|---|
| 58 | * ,GHEISH,GHESIG
|
|---|
| 59 | COMMON /RUNPAC/ DSN,HOST,USER
|
|---|
| 60 | DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
|
|---|
| 61 | REAL STEPFC
|
|---|
| 62 | INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
|
|---|
| 63 | * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
|
|---|
| 64 | * N1STTR,MDBASE
|
|---|
| 65 | INTEGER CETAPE
|
|---|
| 66 | CHARACTER*79 DSN
|
|---|
| 67 | CHARACTER*20 HOST,USER
|
|---|
| 68 |
|
|---|
| 69 | LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
|
|---|
| 70 | * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
|
|---|
| 71 | * ,GHEISH,GHESIG
|
|---|
| 72 | *KEEP,VKIN.
|
|---|
| 73 | COMMON /VKIN/ BETACM
|
|---|
| 74 | DOUBLE PRECISION BETACM
|
|---|
| 75 | *KEND.
|
|---|
| 76 |
|
|---|
| 77 | DOUBLE PRECISION E(10),PHIPAR(10),PL(10),PTR(10),PTSQ(10)
|
|---|
| 78 | DOUBLE PRECISION PX(10),PY(10)
|
|---|
| 79 | DIMENSION ISEQ(10),NTYP(10)
|
|---|
| 80 | C-----------------------------------------------------------------------
|
|---|
| 81 |
|
|---|
| 82 | IF ( DEBUG ) WRITE(MDEBUG,444) (CURPAR(I),I=1,9)
|
|---|
| 83 | 444 FORMAT(' NIHILA: CURPAR=',1P,9E10.3)
|
|---|
| 84 |
|
|---|
| 85 | IREP = 0
|
|---|
| 86 |
|
|---|
| 87 | C RANDOM DECISION FOR ANNIHILATION WITH PROTON OR NEUTRON
|
|---|
| 88 |
|
|---|
| 89 | 20 CONTINUE
|
|---|
| 90 | CALL RMMAR( RD,2,1 )
|
|---|
| 91 | IREP = IREP + 1
|
|---|
| 92 | C AFTER THE 5TH TRY, QUIT THE ANNIHILATION WITHOUT ANY PION GENERATED
|
|---|
| 93 | IF ( IREP .GT. 5 ) GOTO 999
|
|---|
| 94 | IF ( RD(1) .LE. 0.5 ) THEN
|
|---|
| 95 |
|
|---|
| 96 | C-----------------------------------------------------------------------
|
|---|
| 97 | C ANTIPROTON - PROTON AND ANTINEUTRON - NEUTRON ANNIHILATION
|
|---|
| 98 | IF ( RD(2) .LE. CANN(1) ) THEN
|
|---|
| 99 | C ANNIHILATION INTO PI+, PI-
|
|---|
| 100 | NPIPOS = 1
|
|---|
| 101 | NPINEG = 1
|
|---|
| 102 | NPIZ = 0
|
|---|
| 103 | ELSEIF ( RD(2) .LE. CANN(2) ) THEN
|
|---|
| 104 | C ANNIHILATION INTO PI+, PI-, PI0
|
|---|
| 105 | NPIPOS = 1
|
|---|
| 106 | NPINEG = 1
|
|---|
| 107 | NPIZ = 1
|
|---|
| 108 | ELSEIF ( RD(2) .LE. CANN(3) ) THEN
|
|---|
| 109 | C ANNIHILATION INTO PI+, PI-, 2PI0
|
|---|
| 110 | NPIPOS = 1
|
|---|
| 111 | NPINEG = 1
|
|---|
| 112 | NPIZ = 2
|
|---|
| 113 | ELSEIF ( RD(2) .LE. CANN(4) ) THEN
|
|---|
| 114 | C ANNIHILATION INTO PI+, PI-, 3PI0
|
|---|
| 115 | NPIPOS = 1
|
|---|
| 116 | NPINEG = 1
|
|---|
| 117 | NPIZ = 3
|
|---|
| 118 | ELSEIF ( RD(2) .LE. CANN(5) ) THEN
|
|---|
| 119 | C ANNIHILATION INTO PI+, PI-, 4PI0
|
|---|
| 120 | NPIPOS = 1
|
|---|
| 121 | NPINEG = 1
|
|---|
| 122 | NPIZ = 4
|
|---|
| 123 | ELSEIF ( RD(2) .LE. CANN(6) ) THEN
|
|---|
| 124 | C ANNIHILATION INTO 2PI+, 2PI-
|
|---|
| 125 | NPIPOS = 2
|
|---|
| 126 | NPINEG = 2
|
|---|
| 127 | NPIZ = 0
|
|---|
| 128 | ELSEIF ( RD(2) .LE. CANN(7) ) THEN
|
|---|
| 129 | C ANNIHILATION INTO 2PI+, 2PI-, PI0
|
|---|
| 130 | NPIPOS = 2
|
|---|
| 131 | NPINEG = 2
|
|---|
| 132 | NPIZ = 1
|
|---|
| 133 | ELSEIF ( RD(2) .LE. CANN(8) ) THEN
|
|---|
| 134 | C ANNIHILATION INTO 2PI+, 2PI-, 2PI0
|
|---|
| 135 | NPIPOS = 2
|
|---|
| 136 | NPINEG = 2
|
|---|
| 137 | NPIZ = 2
|
|---|
| 138 | ELSEIF ( RD(2) .LE. CANN(9) ) THEN
|
|---|
| 139 | C ANNIHILATION INTO 2PI+, 2PI-, 3PI0
|
|---|
| 140 | NPIPOS = 2
|
|---|
| 141 | NPINEG = 2
|
|---|
| 142 | NPIZ = 3
|
|---|
| 143 | ELSEIF ( RD(2) .LE. CANN(10) ) THEN
|
|---|
| 144 | C ANNIHILATION INTO 3PI+, 3PI-
|
|---|
| 145 | NPIPOS = 3
|
|---|
| 146 | NPINEG = 3
|
|---|
| 147 | NPIZ = 0
|
|---|
| 148 | ELSEIF ( RD(2) .LE. CANN(11) ) THEN
|
|---|
| 149 | C ANNIHILATION INTO 3PI+, 3PI-, PI0
|
|---|
| 150 | NPIPOS = 3
|
|---|
| 151 | NPINEG = 3
|
|---|
| 152 | NPIZ = 1
|
|---|
| 153 | ELSEIF ( RD(2) .LE. CANN(12) ) THEN
|
|---|
| 154 | C ANNIHILATION INTO 3PI+, 3PI-, 2PI0
|
|---|
| 155 | NPIPOS = 3
|
|---|
| 156 | NPINEG = 3
|
|---|
| 157 | NPIZ = 2
|
|---|
| 158 | ELSE
|
|---|
| 159 | C ANNIHILATION INTO 4PI0
|
|---|
| 160 | NPIPOS = 0
|
|---|
| 161 | NPINEG = 0
|
|---|
| 162 | NPIZ = 4
|
|---|
| 163 | ENDIF
|
|---|
| 164 |
|
|---|
| 165 | ELSE
|
|---|
| 166 | C-----------------------------------------------------------------------
|
|---|
| 167 | C ANTIPROTON - NEUTRON (OR ANTINEUTRON - PROTON) ANNIHILATION
|
|---|
| 168 | IF ( RD(2) .LE. CANN(13) ) THEN
|
|---|
| 169 | C ANNIHILATION INTO PI-, PI0
|
|---|
| 170 | NPIPOS = 0
|
|---|
| 171 | NPINEG = 1
|
|---|
| 172 | NPIZ = 1
|
|---|
| 173 | ELSEIF ( RD(2) .LE. CANN(14) ) THEN
|
|---|
| 174 | C ANNIHILATION INTO PI-, 2PI0
|
|---|
| 175 | NPIPOS = 0
|
|---|
| 176 | NPINEG = 1
|
|---|
| 177 | NPIZ = 2
|
|---|
| 178 | ELSEIF ( RD(2) .LE. CANN(15) ) THEN
|
|---|
| 179 | C ANNIHILATION INTO PI-, 3PI0
|
|---|
| 180 | NPIPOS = 0
|
|---|
| 181 | NPINEG = 1
|
|---|
| 182 | NPIZ = 3
|
|---|
| 183 | ELSEIF ( RD(2) .LE. CANN(16) ) THEN
|
|---|
| 184 | C ANNIHILATION INTO PI-, 4PI0
|
|---|
| 185 | NPIPOS = 0
|
|---|
| 186 | NPINEG = 1
|
|---|
| 187 | NPIZ = 4
|
|---|
| 188 | ELSEIF ( RD(2) .LE. CANN(17) ) THEN
|
|---|
| 189 | C ANNIHILATION INTO PI-, 5PI0
|
|---|
| 190 | NPIPOS = 0
|
|---|
| 191 | NPINEG = 1
|
|---|
| 192 | NPIZ = 5
|
|---|
| 193 | ELSEIF ( RD(2) .LE. CANN(18) ) THEN
|
|---|
| 194 | C ANNIHILATION INTO PI+, 2PI-
|
|---|
| 195 | NPIPOS = 1
|
|---|
| 196 | NPINEG = 2
|
|---|
| 197 | NPIZ = 0
|
|---|
| 198 | ELSEIF ( RD(2) .LE. CANN(19) ) THEN
|
|---|
| 199 | C ANNIHILATION INTO PI+, 2PI-, PI0
|
|---|
| 200 | NPIPOS = 1
|
|---|
| 201 | NPINEG = 2
|
|---|
| 202 | NPIZ = 1
|
|---|
| 203 | ELSEIF ( RD(2) .LE. CANN(20) ) THEN
|
|---|
| 204 | C ANNIHILATION INTO PI+, 2PI-, 2PI0
|
|---|
| 205 | NPIPOS = 1
|
|---|
| 206 | NPINEG = 2
|
|---|
| 207 | NPIZ = 2
|
|---|
| 208 | ELSEIF ( RD(2) .LE. CANN(21) ) THEN
|
|---|
| 209 | C ANNIHILATION INTO PI+, 2PI-, 3PI0
|
|---|
| 210 | NPIPOS = 1
|
|---|
| 211 | NPINEG = 2
|
|---|
| 212 | NPIZ = 3
|
|---|
| 213 | ELSEIF ( RD(2) .LE. CANN(22) ) THEN
|
|---|
| 214 | C ANNIHILATION INTO PI+, 2PI-, 4PI0
|
|---|
| 215 | NPIPOS = 1
|
|---|
| 216 | NPINEG = 2
|
|---|
| 217 | NPIZ = 4
|
|---|
| 218 | ELSEIF ( RD(2) .LE. CANN(23) ) THEN
|
|---|
| 219 | C ANNIHILATION INTO 2PI+, 3PI-
|
|---|
| 220 | NPIPOS = 2
|
|---|
| 221 | NPINEG = 3
|
|---|
| 222 | NPIZ = 0
|
|---|
| 223 | ELSEIF ( RD(2) .LE. CANN(24) ) THEN
|
|---|
| 224 | C ANNIHILATION INTO 2PI+, 3PI-, PI0
|
|---|
| 225 | NPIPOS = 2
|
|---|
| 226 | NPINEG = 3
|
|---|
| 227 | NPIZ = 1
|
|---|
| 228 | ELSEIF ( RD(2) .LE. CANN(25) ) THEN
|
|---|
| 229 | C ANNIHILATION INTO 2PI+, 3PI-, 2PI0
|
|---|
| 230 | NPIPOS = 2
|
|---|
| 231 | NPINEG = 3
|
|---|
| 232 | NPIZ = 2
|
|---|
| 233 | ELSEIF ( RD(2) .LE. CANN(26) ) THEN
|
|---|
| 234 | C ANNIHILATION INTO 2PI+, 3PI-, 3PI0
|
|---|
| 235 | NPIPOS = 2
|
|---|
| 236 | NPINEG = 3
|
|---|
| 237 | NPIZ = 3
|
|---|
| 238 | ELSE
|
|---|
| 239 | C ANNIHILATION INTO 3PI+, 4PI-
|
|---|
| 240 | NPIPOS = 3
|
|---|
| 241 | NPINEG = 4
|
|---|
| 242 | NPIZ = 0
|
|---|
| 243 | ENDIF
|
|---|
| 244 |
|
|---|
| 245 | C CHARGE INVERSION IF ANTINEUTRON ANNIHILATES WITH PROTON
|
|---|
| 246 | IF ( ITYPE .EQ. 25 ) THEN
|
|---|
| 247 | NPINEG = NPINEG - 1
|
|---|
| 248 | NPIPOS = NPIPOS + 1
|
|---|
| 249 | ENDIF
|
|---|
| 250 |
|
|---|
| 251 | ENDIF
|
|---|
| 252 |
|
|---|
| 253 | NPI = NPIPOS + NPINEG + NPIZ
|
|---|
| 254 | FNPI = 1.D0 / NPI
|
|---|
| 255 | GCMI = 1.D0 / GCM
|
|---|
| 256 |
|
|---|
| 257 | C-----------------------------------------------------------------------
|
|---|
| 258 | C CHARGE ASSIGNMENT
|
|---|
| 259 |
|
|---|
| 260 | DO 26 I = 1,NPI
|
|---|
| 261 | IF ( I .LE. NPIZ ) THEN
|
|---|
| 262 | C NEUTRAL PIONS
|
|---|
| 263 | NTYP(I) = 7
|
|---|
| 264 | ELSEIF ( I .LE. NPIZ+NPIPOS ) THEN
|
|---|
| 265 | C POSITIVE PIONS
|
|---|
| 266 | NTYP(I) = 8
|
|---|
| 267 | ELSE
|
|---|
| 268 | C NEGATIVE PIONS
|
|---|
| 269 | NTYP(I) = 9
|
|---|
| 270 | ENDIF
|
|---|
| 271 | 26 CONTINUE
|
|---|
| 272 |
|
|---|
| 273 | C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|---|
| 274 | C KINEMATIC CALCULATIONS
|
|---|
| 275 |
|
|---|
| 276 | ISCALE = 0
|
|---|
| 277 | 27 CONTINUE
|
|---|
| 278 | ISCALE = ISCALE + 1
|
|---|
| 279 | C AFTER THE 5TH TRY, TAKE A NEW SET OF PIONS
|
|---|
| 280 | IF ( ISCALE .GT. 5 ) GOTO 20
|
|---|
| 281 |
|
|---|
| 282 | C DISTRIBUTUION OF TRANSVERSE MOMENTA PTR
|
|---|
| 283 |
|
|---|
| 284 | CORECT = 5.3333333333D0 * FNPI**1.5D0 * ECM
|
|---|
| 285 | DO 28 I = 1,NPI
|
|---|
| 286 | PTR(I) = PTRANS(DUMMY) * (1.33333333D0 + CORECT)
|
|---|
| 287 | 28 CONTINUE
|
|---|
| 288 |
|
|---|
| 289 | SUMPX = 0.D0
|
|---|
| 290 | SUMPY = 0.D0
|
|---|
| 291 | CALL RMMAR( RD,NPI,1 )
|
|---|
| 292 | DO 29 I = 1,NPI
|
|---|
| 293 | C SELECT EMISSION ANGLE BY REDUCED RESIDUAL DIRECTION
|
|---|
| 294 | IF ( SUMPX .NE. 0.D0 .OR. SUMPY .NE. 0.D0 ) THEN
|
|---|
| 295 | PHISUM = ATAN2( SUMPY, SUMPX )
|
|---|
| 296 | ELSE
|
|---|
| 297 | PHISUM = 0.D0
|
|---|
| 298 | ENDIF
|
|---|
| 299 | PHIPAR(I) = PHISUM + PI + PI * (NPI+1-I) * (2.*RD(I)-1.) * FNPI
|
|---|
| 300 | PX(I) = COS( PHIPAR(I) ) * PTR(I)
|
|---|
| 301 | PY(I) = SIN( PHIPAR(I) ) * PTR(I)
|
|---|
| 302 | SUMPX = SUMPX + PX(I)
|
|---|
| 303 | SUMPY = SUMPY + PY(I)
|
|---|
| 304 | 29 CONTINUE
|
|---|
| 305 |
|
|---|
| 306 | C CORRECTION OF TRANSVERSE MOMENTA TO KEEP TRANSVERSE MOMENTUM BALANCE
|
|---|
| 307 | SUMPT2 = ECM
|
|---|
| 308 | DPX = SUMPX * FNPI
|
|---|
| 309 | DPY = SUMPY * FNPI
|
|---|
| 310 | DO 30 I = 1,NPI
|
|---|
| 311 | PX(I) = PX(I) - DPX
|
|---|
| 312 | PY(I) = PY(I) - DPY
|
|---|
| 313 | IF ( PX(I) .NE. 0.D0 .OR. PY(I) .NE. 0.D0 ) THEN
|
|---|
| 314 | PHIPAR(I) = ATAN2( PY(I), PX(I) )
|
|---|
| 315 | ELSE
|
|---|
| 316 | PHIPAR(I) = 0.D0
|
|---|
| 317 | ENDIF
|
|---|
| 318 | PTSQ(I) = PX(I)**2 + PY(I)**2
|
|---|
| 319 | SUMPT2 = SUMPT2 - SQRT( PAMA(NTYP(I))**2 + PTSQ(I) )
|
|---|
| 320 | 30 CONTINUE
|
|---|
| 321 |
|
|---|
| 322 | C CHECK, IF C.M. ENERGY IS EXHAUSTED BY TRANSVERSE MOMENTA
|
|---|
| 323 | C IF SO, TRY ANOTHER SET OF TRANSVERSE MOMENTA
|
|---|
| 324 | IF ( SUMPT2 .LE. 0.D0 ) GOTO 27
|
|---|
| 325 |
|
|---|
| 326 | C DISTRIBUTION OF LONGITUDINAL MOMENTA PL
|
|---|
| 327 |
|
|---|
| 328 | C SUM1PL IS SUM OF ABS. VALUES OF LONGITUDINAL MOMENTA
|
|---|
| 329 | F = SUMPT2 * FNPI
|
|---|
| 330 | SUM1PL = 0.D0
|
|---|
| 331 | DO 31 I = 1,NPI
|
|---|
| 332 | FWHM = F + 0.5D0 * SQRT( PTSQ(I) )
|
|---|
| 333 | PL(I) = ABS( RANNOR(0.D0,FWHM) )
|
|---|
| 334 | SUM1PL = SUM1PL + PL(I)
|
|---|
| 335 | C SET SEQUENCE COUNTER
|
|---|
| 336 | ISEQ(I) = I
|
|---|
| 337 | 31 CONTINUE
|
|---|
| 338 |
|
|---|
| 339 | C SORT ISEQ IN DECREASING SIZE OF THE LONGITUDINAL MOMENTUM
|
|---|
| 340 | DO 33 I = 1,NPI
|
|---|
| 341 | DO 32 K = I+1,NPI
|
|---|
| 342 | IF ( PL(ISEQ(I)) .LT. PL(ISEQ(K)) ) THEN
|
|---|
| 343 | IHELP = ISEQ(I)
|
|---|
| 344 | ISEQ(I) = ISEQ(K)
|
|---|
| 345 | ISEQ(K) = IHELP
|
|---|
| 346 | ENDIF
|
|---|
| 347 | 32 CONTINUE
|
|---|
| 348 | 33 CONTINUE
|
|---|
| 349 |
|
|---|
| 350 | C TRY TO BALANCE LONGITUDINAL MOMENTA (TO MINIMIZE CORRECTIONS)
|
|---|
| 351 | C START WITH LONG. MOMENTA IN FORWARD/BACKWARD DIRECTION BY RANDOM
|
|---|
| 352 | CALL RMMAR( RD,1,1 )
|
|---|
| 353 | IF ( RD(1) .LE. 0.5 ) PL(ISEQ(1)) = - PL(ISEQ(1))
|
|---|
| 354 | SUMPL = PL(ISEQ(1))
|
|---|
| 355 | DO 34 I = 2,NPI
|
|---|
| 356 | SUM1PL = SUM1PL - PL(ISEQ(I))
|
|---|
| 357 | C IF THERE IS NOT ENOUGH MOMENTUM LEFT, SELECT FORWARD/BACKWARD TO
|
|---|
| 358 | C BALANCE MOMENTUM, ELSE CHOOSE DIRECTION BY RANDOM
|
|---|
| 359 | IF ( PL(ISEQ(I))+ABS(SUMPL) .GT. SUM1PL ) THEN
|
|---|
| 360 | IF ( PL(ISEQ(I))*SUMPL .GT. 0.D0 ) PL(ISEQ(I)) = -PL(ISEQ(I))
|
|---|
| 361 | ELSE
|
|---|
| 362 | CALL RMMAR( RD,1,1 )
|
|---|
| 363 | IF ( RD(1) .LE. 0.5 ) PL(ISEQ(I)) = - PL(ISEQ(I))
|
|---|
| 364 | ENDIF
|
|---|
| 365 | SUMPL = SUMPL + PL(ISEQ(I))
|
|---|
| 366 | 34 CONTINUE
|
|---|
| 367 |
|
|---|
| 368 | C CORRECTION OF LONGITUDINAL MOMENTA TO KEEP MOMENTUM BALANCE
|
|---|
| 369 | DPL = SUMPL * FNPI
|
|---|
| 370 | DO 35 I = 1,NPI
|
|---|
| 371 | PL(I) = PL(I) - DPL
|
|---|
| 372 | 35 CONTINUE
|
|---|
| 373 |
|
|---|
| 374 | C ITERATIVE CORRECTION OF ALL MOMENTA TO KEEP ENERGY BALANCE
|
|---|
| 375 |
|
|---|
| 376 | IREPET = 0
|
|---|
| 377 | 36 CONTINUE
|
|---|
| 378 | IREPET = IREPET + 1
|
|---|
| 379 | IF ( IREPET .GT. 10 ) GOTO 27
|
|---|
| 380 | ETOT = 0.D0
|
|---|
| 381 | C CHECK ENERGY CONSERVATION
|
|---|
| 382 | DO 37 I = 1,NPI
|
|---|
| 383 | PTSQ(I) = PX(I)**2 + PY(I)**2
|
|---|
| 384 | E(I) = SQRT( PAMA(NTYP(I))**2 + PTSQ(I) + PL(I)**2 )
|
|---|
| 385 | ETOT = ETOT + E(I)
|
|---|
| 386 | 37 CONTINUE
|
|---|
| 387 |
|
|---|
| 388 | ECORR = ECM / ETOT - 1.D0
|
|---|
| 389 |
|
|---|
| 390 | C LOOK WHETHER ENERGY IS CONSERVED WITHIN 1 %
|
|---|
| 391 | IF ( ABS(ECORR) .GT. .01D0 ) THEN
|
|---|
| 392 | C FACTOR IS MODIFIED WITH EMPIRICAL TERM 1/GCM FOR FASTER CONVERGENCE
|
|---|
| 393 | FACT = (0.5D0+GCMI) * ECORR * 0.02D0 * NPI
|
|---|
| 394 | DO 38 I = 1,NPI
|
|---|
| 395 | PX(I) = PX(I) * ( FACT + 1.D0 )
|
|---|
| 396 | PY(I) = PY(I) * ( FACT + 1.D0 )
|
|---|
| 397 | PL(I) = PL(I) * ( FACT * 20.D0 + 1.D0 )
|
|---|
| 398 | 38 CONTINUE
|
|---|
| 399 | GOTO 36
|
|---|
| 400 | ENDIF
|
|---|
| 401 |
|
|---|
| 402 | C LORENTZ TRANSFORMATION FROM C.M. TO LAB. FRAME
|
|---|
| 403 |
|
|---|
| 404 | DO 40 K = 5,8
|
|---|
| 405 | SECPAR(K) = CURPAR(K)
|
|---|
| 406 | 40 CONTINUE
|
|---|
| 407 | DO 41 I = 1,NPI
|
|---|
| 408 | PLLAB = GCM * ( PL(I) + BETACM * E(I) )
|
|---|
| 409 | IF ( PLLAB .LE. 0.D0 ) GOTO 41
|
|---|
| 410 | CTHETA = PLLAB / SQRT( PTSQ(I) + PLLAB**2 )
|
|---|
| 411 | IF ( CTHETA .LT. C(27) ) GOTO 41
|
|---|
| 412 | CALL ADDANG( COSTHE,PHI, CTHETA,PHIPAR(I), SECPAR(3),SECPAR(4) )
|
|---|
| 413 | IF ( SECPAR(3) .LT. C(29) ) GOTO 41
|
|---|
| 414 | SECPAR(1) = NTYP(I)
|
|---|
| 415 | SECPAR(2) = GCM / PAMA(NTYP(I)) * ( PL(I) * BETACM + E(I) )
|
|---|
| 416 | CALL TSTACK
|
|---|
| 417 | 41 CONTINUE
|
|---|
| 418 |
|
|---|
| 419 | MSMM = MSMM + NPI
|
|---|
| 420 | 999 CONTINUE
|
|---|
| 421 |
|
|---|
| 422 | RETURN
|
|---|
| 423 | END
|
|---|
