File indexing completed on 2024-04-06 12:13:19
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092
0093
0094
0095
0096
0097
0098
0099
0100
0101
0102
0103
0104
0105
0106
0107
0108
0109
0110
0111
0112
0113
0114
0115
0116
0117
0118
0119
0120
0121
0122
0123
0124
0125
0126
0127
0128
0129
0130
0131
0132
0133
0134
0135
0136
0137
0138
0139
0140
0141
0142
0143
0144
0145
0146
0147
0148
0149
0150
0151
0152
0153
0154
0155
0156
0157
0158
0159
0160
0161
0162
0163
0164
0165
0166
0167
0168 SUBROUTINE ARTMN
0169
0170
0171
0172
0173
0174
0175
0176
0177
0178
0179
0180
0181
0182 PARAMETER (MAXSTR=150001,MAXR=1,AMU= 0.9383,
0183 1 AKA=0.498,etaM=0.5475)
0184 PARAMETER (MAXX = 20, MAXZ = 24)
0185 PARAMETER (ISUM = 1001, IGAM = 1100)
0186 parameter (MX=4,MY=4,MZ=8,MPX=4,MPY=4,mpz=10,mpzp=10)
0187
0188
0189 INTEGER OUTPAR, zta,zpr
0190 COMMON /AA/ R(3,MAXSTR)
0191
0192 COMMON /BB/ P(3,MAXSTR)
0193
0194 COMMON /CC/ E(MAXSTR)
0195
0196 COMMON /DD/ RHO(-MAXX:MAXX,-MAXX:MAXX,-MAXZ:MAXZ),
0197 & RHOP(-MAXX:MAXX,-MAXX:MAXX,-MAXZ:MAXZ),
0198 & RHON(-MAXX:MAXX,-MAXX:MAXX,-MAXZ:MAXZ)
0199
0200 COMMON /EE/ ID(MAXSTR),LB(MAXSTR)
0201
0202 COMMON /HH/ PROPER(MAXSTR)
0203
0204 common /ff/f(-mx:mx,-my:my,-mz:mz,-mpx:mpx,-mpy:mpy,-mpz:mpzp)
0205
0206 common /gg/ dx,dy,dz,dpx,dpy,dpz
0207
0208 COMMON /INPUT/ NSTAR,NDIRCT,DIR
0209
0210 COMMON /PP/ PRHO(-20:20,-24:24)
0211 COMMON /QQ/ PHRHO(-MAXZ:MAXZ,-24:24)
0212 COMMON /RR/ MASSR(0:MAXR)
0213
0214 common /ss/ inout(20)
0215
0216 common /zz/ zta,zpr
0217
0218 COMMON /RUN/ NUM
0219
0220
0221
0222 COMMON /KKK/ TKAON(7),EKAON(7,0:2000)
0223
0224 COMMON /KAON/ AK(3,50,36),SPECK(50,36,7),MF
0225
0226 COMMON/TABLE/ xarray(0:1000),earray(0:1000)
0227
0228 common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
0229
0230 COMMON /DDpi/ piRHO(-MAXX:MAXX,-MAXX:MAXX,-MAXZ:MAXZ)
0231
0232 common /tt/ PEL(-maxx:maxx,-maxx:maxx,-maxz:maxz)
0233 &,rxy(-maxx:maxx,-maxx:maxx,-maxz:maxz)
0234
0235
0236
0237 DIMENSION TEMP(3,MAXSTR),SKAON(7),SEKAON(7,0:2000)
0238
0239 COMMON /INPUT2/ ILAB, MANYB, NTMAX, ICOLL, INSYS, IPOT, MODE,
0240 & IMOMEN, NFREQ, ICFLOW, ICRHO, ICOU, KPOTEN, KMUL
0241
0242 COMMON /INPUT3/ PLAB, ELAB, ZEROPT, B0, BI, BM, DENCUT, CYCBOX
0243
0244
0245
0246 COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50)
0247
0248
0249
0250
0251
0252
0253
0254
0255
0256
0257
0258 COMMON /ARERCP/PRO1(MAXSTR, MAXR)
0259
0260 COMMON /ARERC1/MULTI1(MAXR)
0261
0262 COMMON /ARPRC1/ITYP1(MAXSTR, MAXR),
0263 & GX1(MAXSTR, MAXR), GY1(MAXSTR, MAXR), GZ1(MAXSTR, MAXR),
0264 & FT1(MAXSTR, MAXR),
0265 & PX1(MAXSTR, MAXR), PY1(MAXSTR, MAXR), PZ1(MAXSTR, MAXR),
0266 & EE1(MAXSTR, MAXR), XM1(MAXSTR, MAXR)
0267
0268
0269 DIMENSION NPI(MAXR)
0270 DIMENSION RT(3, MAXSTR, MAXR), PT(3, MAXSTR, MAXR)
0271 & , ET(MAXSTR, MAXR), LT(MAXSTR, MAXR), PROT(MAXSTR, MAXR)
0272
0273 EXTERNAL IARFLV, INVFLV
0274
0275 common /lastt/itimeh,bimp
0276
0277 common/snn/efrm,npart1,npart2,epsiPz,epsiPt,PZPROJ,PZTARG
0278
0279 COMMON/hbt/lblast(MAXSTR),xlast(4,MAXSTR),plast(4,MAXSTR),nlast
0280
0281 common/resdcy/NSAV,iksdcy
0282
0283 COMMON/RNDF77/NSEED
0284
0285 COMMON/FTMAX/ftsv(MAXSTR),ftsvt(MAXSTR, MAXR)
0286 COMMON /dpert/dpertt(MAXSTR,MAXR),dpertp(MAXSTR),dplast(MAXSTR),
0287 1 dpdcy(MAXSTR),dpdpi(MAXSTR,MAXR),dpt(MAXSTR, MAXR),
0288 2 dpp1(MAXSTR,MAXR),dppion(MAXSTR,MAXR)
0289 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
0290
0291 real zet(-45:45)
0292 SAVE
0293 data zet /
0294 4 1.,0.,0.,0.,0.,
0295 3 1.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
0296 2 -1.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
0297 1 0.,0.,0.,-1.,0.,1.,0.,-1.,0.,-1.,
0298 s 0.,-2.,-1.,0.,1.,0.,0.,0.,0.,-1.,
0299 e 0.,
0300 s 1.,0.,-1.,0.,1.,-1.,0.,1.,2.,0.,
0301 1 1.,0.,1.,0.,-1.,0.,1.,0.,0.,0.,
0302 2 -1.,0.,1.,0.,-1.,0.,1.,0.,0.,1.,
0303 3 0.,0.,0.,0.,0.,0.,0.,0.,0.,-1.,
0304 4 0.,0.,0.,0.,-1./
0305
0306 nlast=0
0307 do 1002 i=1,MAXSTR
0308 ftsv(i)=0.
0309 do 1101 irun=1,maxr
0310 ftsvt(i,irun)=0.
0311 1101 continue
0312 lblast(i)=999
0313 do 1001 j=1,4
0314
0315
0316
0317 xlast(j,i)=0.
0318 plast(j,i)=0.
0319 1001 continue
0320 1002 continue
0321
0322
0323
0324
0325
0326
0327
0328
0329
0330
0331
0332
0333
0334
0335
0336
0337 call tablem
0338
0339 ikaon=1
0340 nstar=1
0341 ndirct=0
0342 dir=0.02
0343 asy=0.032
0344 ESBIN=0.04
0345 MF=36
0346
0347
0348
0349 RADTA = 1.124 * FLOAT(MASSTA)**(1./3.)
0350 RADPR = 1.124 * FLOAT(MASSPR)**(1./3.)
0351 ZDIST = RADTA + RADPR
0352
0353 BMAX = RADTA + RADPR
0354 MASS = MASSTA + MASSPR
0355 NTOTAL = NUM * MASS
0356
0357 IF (NTOTAL .GT. MAXSTR) THEN
0358 WRITE(12,'(//10X,''**** FATAL ERROR: TOO MANY TEST PART. ****'//
0359 & ' '')')
0360 STOP
0361 END IF
0362
0363
0364
0365
0366
0367
0368 ETA = FLOAT(MASSTA) * AMU
0369 PZTA = 0.0
0370 BETATA = 0.0
0371 GAMMTA = 1.0
0372
0373 EPR = FLOAT(MASSPR) * (AMU + 0.001 * ELAB)
0374 PZPR = SQRT( EPR**2 - (AMU * FLOAT(MASSPR))**2 )
0375 BETAPR = PZPR / EPR
0376 GAMMPR = 1.0 / SQRT( 1.0 - BETAPR**2 )
0377
0378
0379 BETAC=(PZPR+PZTA)/(EPR+ETA)
0380 GAMMC=1.0 / SQRT(1.-BETAC**2)
0381
0382
0383
0384
0385
0386
0387
0388 IF (INSYS .NE. 0) THEN
0389
0390
0391
0392 S = (EPR+ETA)**2 - PZPR**2
0393 xx1=4.*alog(float(massta))
0394 xx2=4.*alog(float(masspr))
0395 xx1=exp(xx1)
0396 xx2=exp(xx2)
0397 PSQARE = (S**2 + (xx1+ xx2) * AMU**4
0398 & - 2.0 * S * AMU**2 * FLOAT(MASSTA**2 + MASSPR**2)
0399 & - 2.0 * FLOAT(MASSTA**2 * MASSPR**2) * AMU**4)
0400 & / (4.0 * S)
0401
0402 ETA = SQRT ( PSQARE + (FLOAT(MASSTA) * AMU)**2 )
0403 PZTA = - SQRT(PSQARE)
0404 BETATA = PZTA / ETA
0405 GAMMTA = 1.0 / SQRT( 1.0 - BETATA**2 )
0406
0407 EPR = SQRT ( PSQARE + (FLOAT(MASSPR) * AMU)**2 )
0408 PZPR = SQRT(PSQARE)
0409 BETAPR = PZPR/ EPR
0410 GAMMPR = 1.0 / SQRT( 1.0 - BETAPR**2 )
0411
0412
0413
0414
0415
0416
0417
0418
0419
0420 ELSE
0421
0422 END IF
0423
0424 PZTA = PZTA / FLOAT(MASSTA)
0425 PZPR = PZPR / FLOAT(MASSPR)
0426
0427 ECMS0=ETA+EPR
0428
0429
0430
0431
0432
0433 DO 50000 IMANY=1,MANYB
0434
0435
0436 if (manyb. gt.1) then
0437 111 BX=1.0-2.0*RANART(NSEED)
0438 BY=1.0-2.0*RANART(NSEED)
0439 B2=BX*BX+BY*BY
0440 IF(B2.GT.1.0) GO TO 111
0441 B=SQRT(B2)*(BM-BI)+BI
0442 ELSE
0443 B=B0
0444 ENDIF
0445
0446
0447
0448
0449
0450
0451
0452 call coulin(masspr,massta,NUM)
0453
0454 CALL INIT(1 ,MASSTA ,NUM ,RADTA,
0455 & B/2. ,ZEROPT+ZDIST/2. ,PZTA,
0456 & GAMMTA ,ISEED ,MASS ,IMOMEN)
0457
0458 CALL INIT(1+MASSTA,MASS ,NUM ,RADPR,
0459 & -B/2. ,ZEROPT-ZDIST/2. ,PZPR,
0460 & GAMMPR ,ISEED ,MASS ,IMOMEN)
0461
0462 OUTPAR = 0
0463
0464
0465 MASSR(0)=0
0466 DO 1003 IR =1,NUM
0467 MASSR(IR)=MASS
0468 1003 CONTINUE
0469
0470
0471
0472 CALL DENS(IPOT,MASS,NUM,OUTPAR)
0473
0474
0475
0476
0477
0478
0479
0480
0481
0482
0483
0484
0485
0486
0487
0488
0489
0490 IF (ICOLL .NE. -1) THEN
0491 DO 700 I = 1,NTOTAL
0492 IX = NINT( R(1,I) )
0493 IY = NINT( R(2,I) )
0494 IZ = NINT( R(3,I) )
0495
0496 IF(IX.GE.MAXX.OR.IY.GE.MAXX.OR.IZ.GE.MAXZ
0497 1 .OR.IX.LE.-MAXX.OR.IY.LE.-MAXX.OR.IZ.LE.-MAXZ) goto 700
0498 CALL GRADU(IPOT,IX,IY,IZ,GRADX,GRADY,GRADZ)
0499 P(1,I) = P(1,I) - (0.5 * DT) * GRADX
0500 P(2,I) = P(2,I) - (0.5 * DT) * GRADY
0501 P(3,I) = P(3,I) - (0.5 * DT) * GRADZ
0502 700 CONTINUE
0503 END IF
0504
0505
0506
0507
0508
0509 RCNNE = 0
0510 RDD = 0
0511 RPP = 0
0512 rppk = 0
0513 RPN = 0
0514 rpd = 0
0515 RKN = 0
0516 RNNK = 0
0517 RDDK = 0
0518 RNDK = 0
0519 RCNND = 0
0520 RCNDN = 0
0521 RCOLL = 0
0522 RBLOC = 0
0523 RDIRT = 0
0524 RDECAY = 0
0525 RRES = 0
0526
0527 DO 1005 KKK=1,5
0528 SKAON(KKK) = 0
0529 DO 1004 IS=1,2000
0530 SEKAON(KKK,IS)=0
0531 1004 CONTINUE
0532 1005 CONTINUE
0533
0534 pr0=0.
0535 pr1=0.
0536 ska0=0.
0537 ska1=0.
0538
0539
0540
0541
0542 IF (IAPAR2(1) .NE. 1) THEN
0543 DO 1016 I = 1, MAXSTR
0544 DO 1015 J = 1, 3
0545 R(J, I) = 0.
0546 P(J, I) = 0.
0547 1015 CONTINUE
0548 E(I) = 0.
0549 LB(I) = 0
0550
0551 ID(I)=0
0552
0553 PROPER(I) = 1.
0554 1016 CONTINUE
0555 MASS = 0
0556
0557
0558
0559
0560 NP = 0
0561 DO 1017 J = 1, NUM
0562 MASSR(J) = 0
0563 NPI(J) = 1
0564 1017 CONTINUE
0565 DO 1019 I = 1, MAXR
0566 DO 1018 J = 1, MAXSTR
0567 RT(1, J, I) = 0.
0568 RT(2, J, I) = 0.
0569 RT(3, J, I) = 0.
0570 PT(1, J, I) = 0.
0571 PT(2, J, I) = 0.
0572 PT(3, J, I) = 0.
0573 ET(J, I) = 0.
0574 LT(J, I) = 0
0575
0576 PROT(J, I) = 1.
0577 1018 CONTINUE
0578 1019 CONTINUE
0579
0580 END IF
0581
0582
0583 DO 10000 NT = 1,NTMAX
0584
0585
0586 LP1=0
0587 LP2=0
0588 LP3=0
0589
0590 LD1=0
0591 LD2=0
0592 LD3=0
0593 LD4=0
0594
0595 LN1=0
0596 LN2=0
0597
0598 LN5=0
0599
0600 LE=0
0601
0602 LKAON=0
0603
0604
0605
0606 LKAONS=0
0607
0608
0609 IF (ICOLL .NE. 1) THEN
0610
0611
0612 numnt=nt
0613 CALL RELCOL(LCOLL,LBLOC,LCNNE,LDD,LPP,lppk,
0614 & LPN,lpd,LRHO,LOMEGA,LKN,LNNK,LDDK,LNDK,LCNND,
0615 & LCNDN,LDIRT,LDECAY,LRES,LDOU,LDDRHO,LNNRHO,
0616 & LNNOM,numnt,ntmax,sp,akaon,sk)
0617
0618
0619
0620
0621
0622
0623
0624
0625
0626
0627
0628
0629
0630
0631
0632
0633
0634
0635
0636
0637
0638
0639
0640
0641
0642
0643
0644
0645
0646
0647 RCOLL = RCOLL + FLOAT(LCOLL)/num
0648 RBLOC = RBLOC + FLOAT(LBLOC)/num
0649 RCNNE = RCNNE + FLOAT(LCNNE)/num
0650 RDD = RDD + FLOAT(LDD)/num
0651 RPP = RPP + FLOAT(LPP)/NUM
0652 rppk =rppk + float(lppk)/num
0653 RPN = RPN + FLOAT(LPN)/NUM
0654 rpd =rpd + float(lpd)/num
0655 RKN = RKN + FLOAT(LKN)/NUM
0656 RNNK =RNNK + FLOAT(LNNK)/NUM
0657 RDDK =RDDK + FLOAT(LDDK)/NUM
0658 RNDK =RNDK + FLOAT(LNDK)/NUM
0659 RCNND = RCNND + FLOAT(LCNND)/num
0660 RCNDN = RCNDN + FLOAT(LCNDN)/num
0661 RDIRT = RDIRT + FLOAT(LDIRT)/num
0662 RDECAY= RDECAY+ FLOAT(LDECAY)/num
0663 RRES = RRES + FLOAT(LRES)/num
0664
0665 ADIRT=LDIRT/DT/num
0666 ACOLL=(LCOLL-LBLOC)/DT/num
0667 ACNND=LCNND/DT/num
0668 ACNDN=LCNDN/DT/num
0669 ADECAY=LDECAY/DT/num
0670 ARES=LRES/DT/num
0671 ADOU=LDOU/DT/NUM
0672 ADDRHO=LDDRHO/DT/NUM
0673 ANNRHO=LNNRHO/DT/NUM
0674 ANNOM=LNNOM/DT/NUM
0675 ADD=LDD/DT/num
0676 APP=LPP/DT/num
0677 appk=lppk/dt/num
0678 APN=LPN/DT/num
0679 apd=lpd/dt/num
0680 arh=lrho/dt/num
0681 aom=lomega/dt/num
0682 AKN=LKN/DT/num
0683 ANNK=LNNK/DT/num
0684 ADDK=LDDK/DT/num
0685 ANDK=LNDK/DT/num
0686
0687
0688
0689
0690
0691
0692
0693
0694
0695
0696
0697
0698
0699
0700
0701
0702
0703
0704
0705
0706
0707
0708
0709
0710
0711
0712
0713
0714
0715
0716
0717
0718 END IF
0719
0720
0721
0722 CALL DENS(IPOT,MASS,NUM,OUTPAR)
0723
0724
0725
0726 sumene=0
0727 ISO=0
0728 DO 201 MRUN=1,NUM
0729 ISO=ISO+MASSR(MRUN-1)
0730 DO 201 I0=1,MASSR(MRUN)
0731 I =I0+ISO
0732 ETOTAL = SQRT( E(I)**2 + P(1,I)**2 + P(2,I)**2 +P(3,I)**2 )
0733 sumene=sumene+etotal
0734
0735
0736
0737 if(kpoten.ne.0.and.lb(i).eq.23)then
0738 den=0.
0739 IX = NINT( R(1,I) )
0740 IY = NINT( R(2,I) )
0741 IZ = NINT( R(3,I) )
0742
0743
0744
0745 IF(IX.LT.MAXX.AND.IY.LT.MAXX.AND.IZ.LT.MAXZ
0746 1 .AND.IX.GT.-MAXX.AND.IY.GT.-MAXX.AND.IZ.GT.-MAXZ)
0747 2 den=rho(ix,iy,iz)
0748
0749
0750
0751
0752 akg = 0.1727
0753
0754 bkg = 0.333
0755 rnsg = den
0756 ecor = - akg*rnsg + (bkg*den)**2
0757 etotal = sqrt(etotal**2 + ecor)
0758 endif
0759
0760 if(kpoten.ne.0.and.lb(i).eq.21)then
0761 den=0.
0762 IX = NINT( R(1,I) )
0763 IY = NINT( R(2,I) )
0764 IZ = NINT( R(3,I) )
0765
0766
0767
0768 IF(IX.LT.MAXX.AND.IY.LT.MAXX.AND.IZ.LT.MAXZ
0769 1 .AND.IX.GT.-MAXX.AND.IY.GT.-MAXX.AND.IZ.GT.-MAXZ)
0770 2 den=rho(ix,iy,iz)
0771
0772
0773
0774 akg = 0.1727
0775
0776 bkg = 0.333
0777 rnsg = den
0778 ecor = - akg*rnsg + (bkg*den)**2
0779 etotal = sqrt(etotal**2 + ecor)
0780 endif
0781
0782
0783 R(1,I) = R(1,I) + DT*P(1,I)/ETOTAL
0784 R(2,I) = R(2,I) + DT*P(2,I)/ETOTAL
0785 R(3,I) = R(3,I) + DT*P(3,I)/ETOTAL
0786
0787 if ( cycbox.ne.0 ) then
0788 if ( r(1,i).gt. cycbox/2 ) r(1,i)=r(1,i)-cycbox
0789 if ( r(1,i).le.-cycbox/2 ) r(1,i)=r(1,i)+cycbox
0790 if ( r(2,i).gt. cycbox/2 ) r(2,i)=r(2,i)-cycbox
0791 if ( r(2,i).le.-cycbox/2 ) r(2,i)=r(2,i)+cycbox
0792 if ( r(3,i).gt. cycbox/2 ) r(3,i)=r(3,i)-cycbox
0793 if ( r(3,i).le.-cycbox/2 ) r(3,i)=r(3,i)+cycbox
0794 end if
0795
0796 LB1=LB(I)
0797
0798 IF(LB1.EQ.9)LD1=LD1+1
0799
0800 IF(LB1.EQ.8)LD2=LD2+1
0801
0802 IF(LB1.EQ.7)LD3=LD3+1
0803
0804 IF(LB1.EQ.6)LD4=LD4+1
0805
0806 IF(LB1.EQ.11)LN1=LN1+1
0807
0808 IF(LB1.EQ.10)LN2=LN2+1
0809
0810 IF((LB1.EQ.13).OR.(LB1.EQ.12))LN5=LN5+1
0811
0812 IF(LB1.EQ.0)LE=LE+1
0813
0814 IF(LB1.EQ.23)LKAON=LKAON+1
0815
0816 IF(LB1.EQ.30)LKAONS=LKAONS+1
0817
0818
0819
0820 IF(LB1.EQ.5)LP1=LP1+1
0821
0822 IF(LB1.EQ.4)LP2=LP2+1
0823
0824 IF(LB1.EQ.3)LP3=LP3+1
0825 201 CONTINUE
0826 LP=LP1+LP2+LP3
0827 LD=LD1+LD2+LD3+LD4
0828 LN=LN1+LN2
0829 ALP=FLOAT(LP)/FLOAT(NUM)
0830 ALD=FLOAT(LD)/FLOAT(NUM)
0831 ALN=FLOAT(LN)/FLOAT(NUM)
0832 ALN5=FLOAT(LN5)/FLOAT(NUM)
0833 ATOTAL=ALP+ALD+ALN+0.5*ALN5
0834 ALE=FLOAT(LE)/FLOAT(NUM)
0835 ALKAON=FLOAT(LKAON)/FLOAT(NUM)
0836
0837 if (icou .eq. 1) then
0838
0839 iso=0
0840 do 1026 irun = 1,num
0841 iso=iso+massr(irun-1)
0842 do 1021 il = 1,massr(irun)
0843 temp(1,il) = 0.
0844 temp(2,il) = 0.
0845 temp(3,il) = 0.
0846 1021 continue
0847 do 1023 il = 1, massr(irun)
0848 i=iso+il
0849 if (zet(lb(i)).ne.0) then
0850 do 1022 jl = 1,il-1
0851 j=iso+jl
0852 if (zet(lb(j)).ne.0) then
0853 ddx=r(1,i)-r(1,j)
0854 ddy=r(2,i)-r(2,j)
0855 ddz=r(3,i)-r(3,j)
0856 rdiff = sqrt(ddx**2+ddy**2+ddz**2)
0857 if (rdiff .le. 1.) rdiff = 1.
0858 grp=zet(lb(i))*zet(lb(j))/rdiff**3
0859 ddx=ddx*grp
0860 ddy=ddy*grp
0861 ddz=ddz*grp
0862 temp(1,il)=temp(1,il)+ddx
0863 temp(2,il)=temp(2,il)+ddy
0864 temp(3,il)=temp(3,il)+ddz
0865 temp(1,jl)=temp(1,jl)-ddx
0866 temp(2,jl)=temp(2,jl)-ddy
0867 temp(3,jl)=temp(3,jl)-ddz
0868 end if
0869 1022 continue
0870 end if
0871 1023 continue
0872 do 1025 il = 1,massr(irun)
0873 i= iso+il
0874 if (zet(lb(i)).ne.0) then
0875 do 1024 idir = 1,3
0876 p(idir,i) = p(idir,i) + temp(idir,il)
0877 & * dt * 0.00144
0878 1024 continue
0879 end if
0880 1025 continue
0881 1026 continue
0882 end if
0883
0884
0885
0886
0887
0888 spt=0
0889 spz=0
0890 ncen=0
0891 ekin=0
0892 NLOST = 0
0893 MEAN=0
0894 nquark=0
0895 nbaryn=0
0896
0897 rads = 2.
0898 zras = 0.1
0899 denst = 0.
0900 edenst = 0.
0901
0902 DO 6000 IRUN = 1,NUM
0903 MEAN=MEAN+MASSR(IRUN-1)
0904 DO 5800 J = 1,MASSR(irun)
0905 I=J+MEAN
0906
0907
0908 radut = sqrt(r(1,i)**2+r(2,i)**2)
0909 if( radut .le. rads )then
0910 if( abs(r(3,i)) .le. zras*nt*dt )then
0911
0912
0913 vols = 3.14159*rads**2*zras
0914 engs=sqrt(p(1,i)**2+p(2,i)**2+p(3,i)**2+e(i)**2)
0915 gammas=1.
0916 if(e(i).ne.0.)gammas=engs/e(i)
0917
0918 denst = denst + 1./gammas/vols
0919
0920 edenst = edenst + engs/gammas/gammas/vols
0921 endif
0922 endif
0923
0924
0925 drr=sqrt(r(1,i)**2+r(2,i)**2+r(3,i)**2)
0926 if(drr.le.2.0)then
0927 spt=spt+p(1,i)**2+p(2,i)**2
0928 spz=spz+p(3,i)**2
0929 ncen=ncen+1
0930 ekin=ekin+sqrt(p(1,i)**2+p(2,i)**2+p(3,i)**2+e(i)**2)-e(i)
0931 endif
0932 IX = NINT( R(1,I) )
0933 IY = NINT( R(2,I) )
0934 IZ = NINT( R(3,I) )
0935
0936
0937
0938
0939 IF(IX.LT.MAXX.AND.IY.LT.MAXX.AND.IZ.LT.MAXZ
0940 1 .AND.IX.GT.-MAXX.AND.IY.GT.-MAXX.AND.IZ.GT.-MAXZ) THEN
0941 if(rho(ix,iy,iz)/0.168.gt.dencut)go to 5800
0942 if((rho(ix,iy,iz)/0.168.gt.5.).and.(e(i).gt.0.9))
0943 & nbaryn=nbaryn+1
0944 if(pel(ix,iy,iz).gt.2.0)nquark=nquark+1
0945 endif
0946
0947
0948 if(kpoten.ne.0.and.lb(i).eq.23)then
0949 den=0.
0950
0951
0952
0953 IF(IX.LT.MAXX.AND.IY.LT.MAXX.AND.IZ.LT.MAXZ
0954 1 .AND.IX.GT.-MAXX.AND.IY.GT.-MAXX.AND.IZ.GT.-MAXZ) THEN
0955 den=rho(ix,iy,iz)
0956
0957
0958
0959
0960 akg = 0.1727
0961
0962 bkg = 0.333
0963 rnsg = den
0964 ecor = - akg*rnsg + (bkg*den)**2
0965 etotal = sqrt(P(1,i)**2+p(2,I)**2+p(3,i)**2+e(i)**2 + ecor)
0966 ecor = - akg + 2.*bkg**2*den + 2.*bkg*etotal
0967
0968 CALL GRADUK(IX,IY,IZ,GRADXk,GRADYk,GRADZk)
0969 P(1,I) = P(1,I) - DT * GRADXk*ecor/(2.*etotal)
0970 P(2,I) = P(2,I) - DT * GRADYk*ecor/(2.*etotal)
0971 P(3,I) = P(3,I) - DT * GRADZk*ecor/(2.*etotal)
0972 endif
0973 endif
0974
0975 if(kpoten.ne.0.and.lb(i).eq.21)then
0976 den=0.
0977
0978
0979
0980 IF(IX.LT.MAXX.AND.IY.LT.MAXX.AND.IZ.LT.MAXZ
0981 1 .AND.IX.GT.-MAXX.AND.IY.GT.-MAXX.AND.IZ.GT.-MAXZ) THEN
0982 den=rho(ix,iy,iz)
0983 CALL GRADUK(IX,IY,IZ,GRADXk,GRADYk,GRADZk)
0984
0985
0986
0987
0988
0989 akg = 0.1727
0990
0991 bkg = 0.333
0992 rnsg = den
0993 ecor = - akg*rnsg + (bkg*den)**2
0994 etotal = sqrt(P(1,i)**2+p(2,I)**2+p(3,i)**2+e(i)**2 + ecor)
0995 ecor = - akg + 2.*bkg**2*den - 2.*bkg*etotal
0996 P(1,I) = P(1,I) - DT * GRADXk*ecor/(2.*etotal)
0997 P(2,I) = P(2,I) - DT * GRADYk*ecor/(2.*etotal)
0998 P(3,I) = P(3,I) - DT * GRADZk*ecor/(2.*etotal)
0999
1000 endif
1001 endif
1002
1003
1004 if(j.gt.mass)go to 5800
1005
1006
1007
1008
1009 IF (ICOLL .NE. -1) THEN
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023 IF(IX.LT.MAXX.AND.IY.LT.MAXX.AND.IZ.LT.MAXZ
1024 1 .AND.IX.GT.-MAXX.AND.IY.GT.-MAXX.AND.IZ.GT.-MAXZ) THEN
1025 CALL GRADU(IPOT,IX,IY,IZ,GRADX,GRADY,GRADZ)
1026 TZ=0.
1027 GRADXN=0
1028 GRADYN=0
1029 GRADZN=0
1030 GRADXP=0
1031 GRADYP=0
1032 GRADZP=0
1033 IF(ICOU.EQ.1)THEN
1034 CALL GRADUP(IX,IY,IZ,GRADXP,GRADYP,GRADZP)
1035 CALL GRADUN(IX,IY,IZ,GRADXN,GRADYN,GRADZN)
1036 IF(ZET(LB(I)).NE.0)TZ=-1
1037 IF(ZET(LB(I)).EQ.0)TZ= 1
1038 END IF
1039 if(iabs(lb(i)).ge.14.and.iabs(lb(i)).le.17)then
1040 facl = 2./3.
1041 elseif(iabs(lb(i)).eq.40.or.iabs(lb(i)).eq.41)then
1042 facl = 1./3.
1043 else
1044 facl = 1.
1045 endif
1046 P(1,I) = P(1,I) - facl*DT * (GRADX+asy*(GRADXN-GRADXP)*TZ)
1047 P(2,I) = P(2,I) - facl*DT * (GRADY+asy*(GRADYN-GRADYP)*TZ)
1048 P(3,I) = P(3,I) - facl*DT * (GRADZ+asy*(GRADZN-GRADZP)*TZ)
1049 end if
1050 ENDIF
1051
1052 5800 CONTINUE
1053 6000 CONTINUE
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075 CDEN=RHO(0,0,0)/0.168
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093 IF ((NT/NFREQ)*NFREQ .EQ. NT ) THEN
1094 if(icflow.eq.1)call flow(nt)
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110 endif
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145 IF (IAPAR2(1) .NE. 1) THEN
1146 CT = NT * DT
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162 IA = 0
1163 DO 1028 IRUN = 1, NUM
1164 DO 1027 IC = 1, MASSR(IRUN)
1165 IE = IA + IC
1166 RT(1, IC, IRUN) = R(1, IE)
1167 RT(2, IC, IRUN) = R(2, IE)
1168 RT(3, IC, IRUN) = R(3, IE)
1169 PT(1, IC, IRUN) = P(1, IE)
1170 PT(2, IC, IRUN) = P(2, IE)
1171 PT(3, IC, IRUN) = P(3, IE)
1172 ET(IC, IRUN) = E(IE)
1173 LT(IC, IRUN) = LB(IE)
1174
1175 PROT(IC, IRUN) = PROPER(IE)
1176
1177 dpertt(IC, IRUN)=dpertp(IE)
1178 1027 CONTINUE
1179 NP = MASSR(IRUN)
1180 NP1 = NPI(IRUN)
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193 ctlong = ct
1194 if(nt .eq. (ntmax-1))then
1195 ctlong = 1.E30
1196 elseif(nt .eq. ntmax)then
1197 go to 1111
1198 endif
1199 DO WHILE (NP1.LE.MULTI1(IRUN).AND.
1200 & FT1(NP1, IRUN) .GT. (CT - DT) .AND.
1201 & FT1(NP1, IRUN) .LE. ctlong)
1202 NP = NP + 1
1203 UDT = (CT - FT1(NP1, IRUN)) / EE1(NP1, IRUN)
1204
1205
1206 if(nt.eq.(ntmax-1)) then
1207 ftsvt(NP,IRUN)=FT1(NP1, IRUN)
1208 if(FT1(NP1, IRUN).gt.ct) UDT=0.
1209 endif
1210 RT(1, NP, IRUN) = GX1(NP1, IRUN) +
1211 & PX1(NP1, IRUN) * UDT
1212 RT(2, NP, IRUN) = GY1(NP1, IRUN) +
1213 & PY1(NP1, IRUN) * UDT
1214 RT(3, NP, IRUN) = GZ1(NP1, IRUN) +
1215 & PZ1(NP1, IRUN) * UDT
1216 PT(1, NP, IRUN) = PX1(NP1, IRUN)
1217 PT(2, NP, IRUN) = PY1(NP1, IRUN)
1218 PT(3, NP, IRUN) = PZ1(NP1, IRUN)
1219 ET(NP, IRUN) = XM1(NP1, IRUN)
1220 LT(NP, IRUN) = IARFLV(ITYP1(NP1, IRUN))
1221
1222 dpertt(NP,IRUN)=dpp1(NP1,IRUN)
1223
1224
1225
1226
1227
1228
1229
1230 NP1 = NP1 + 1
1231
1232 PROT(NP, IRUN) = 1.
1233 END DO
1234
1235 1111 continue
1236 NPI(IRUN) = NP1
1237 IA = IA + MASSR(IRUN)
1238 MASSR(IRUN) = NP
1239 1028 CONTINUE
1240 IA = 0
1241 DO 1030 IRUN = 1, NUM
1242 IA = IA + MASSR(IRUN - 1)
1243 DO 1029 IC = 1, MASSR(IRUN)
1244 IE = IA + IC
1245 R(1, IE) = RT(1, IC, IRUN)
1246 R(2, IE) = RT(2, IC, IRUN)
1247 R(3, IE) = RT(3, IC, IRUN)
1248 P(1, IE) = PT(1, IC, IRUN)
1249 P(2, IE) = PT(2, IC, IRUN)
1250 P(3, IE) = PT(3, IC, IRUN)
1251 E(IE) = ET(IC, IRUN)
1252 LB(IE) = LT(IC, IRUN)
1253
1254 PROPER(IE) = PROT(IC, IRUN)
1255 if(nt.eq.(ntmax-1)) ftsv(IE)=ftsvt(IC,IRUN)
1256
1257 dpertp(IE)=dpertt(IC, IRUN)
1258 1029 CONTINUE
1259
1260 call hbtout(MASSR(IRUN),nt,ntmax)
1261 1030 CONTINUE
1262
1263 END IF
1264
1265
1266
1267
1268
1269 10000 continue
1270
1271
1272
1273
1274
1275
1276
1277 iss=0
1278 do 1032 lrun=1,num
1279 iss=iss+massr(lrun-1)
1280 do 1031 l0=1,massr(lrun)
1281 ipart=iss+l0
1282 1031 continue
1283 1032 continue
1284
1285
1286 IF (IAPAR2(1) .NE. 1) THEN
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330 IA = 0
1331 DO 1035 IRUN = 1, NUM
1332 IA = IA + MASSR(IRUN - 1)
1333 NP1 = NPI(IRUN)
1334 NSH = MASSR(IRUN) - NP1 + 1
1335 MULTI1(IRUN) = MULTI1(IRUN) + NSH
1336
1337 IF (NSH .GT. 0) THEN
1338 IB = MULTI1(IRUN)
1339 IE = MASSR(IRUN) + 1
1340 II = -1
1341 ELSE IF (NSH .LT. 0) THEN
1342 IB = MASSR(IRUN) + 1
1343 IE = MULTI1(IRUN)
1344 II = 1
1345 END IF
1346 IF (NSH .NE. 0) THEN
1347 DO 1033 I = IB, IE, II
1348 J = I - NSH
1349 ITYP1(I, IRUN) = ITYP1(J, IRUN)
1350 GX1(I, IRUN) = GX1(J, IRUN)
1351 GY1(I, IRUN) = GY1(J, IRUN)
1352 GZ1(I, IRUN) = GZ1(J, IRUN)
1353 FT1(I, IRUN) = FT1(J, IRUN)
1354 PX1(I, IRUN) = PX1(J, IRUN)
1355 PY1(I, IRUN) = PY1(J, IRUN)
1356 PZ1(I, IRUN) = PZ1(J, IRUN)
1357 EE1(I, IRUN) = EE1(J, IRUN)
1358 XM1(I, IRUN) = XM1(J, IRUN)
1359
1360 PRO1(I, IRUN) = PRO1(J, IRUN)
1361
1362 dpp1(I,IRUN)=dpp1(J,IRUN)
1363 1033 CONTINUE
1364 END IF
1365
1366
1367 DO 1034 I = 1, MASSR(IRUN)
1368 IB = IA + I
1369 ITYP1(I, IRUN) = INVFLV(LB(IB))
1370 GX1(I, IRUN) = R(1, IB)
1371 GY1(I, IRUN) = R(2, IB)
1372 GZ1(I, IRUN) = R(3, IB)
1373
1374
1375
1376
1377 if(FT1(I, IRUN).lt.CT) FT1(I, IRUN) = CT
1378 PX1(I, IRUN) = P(1, IB)
1379 PY1(I, IRUN) = P(2, IB)
1380 PZ1(I, IRUN) = P(3, IB)
1381 XM1(I, IRUN) = E(IB)
1382 EE1(I, IRUN) = SQRT(PX1(I, IRUN) ** 2 +
1383 & PY1(I, IRUN) ** 2 +
1384 & PZ1(I, IRUN) ** 2 +
1385 & XM1(I, IRUN) ** 2)
1386
1387 PRO1(I, IRUN) = PROPER(IB)
1388 1034 CONTINUE
1389 1035 CONTINUE
1390
1391 END IF
1392
1393
1394
1395
1396
1397
1398
1399 50000 CONTINUE
1400
1401
1402
1403
1404
1405
1406 RETURN
1407
1408 END
1409
1410 subroutine coulin(masspr,massta,NUM)
1411
1412
1413
1414
1415
1416 integer zta,zpr
1417 PARAMETER (MAXSTR=150001)
1418 common /EE/ ID(MAXSTR),LB(MAXSTR)
1419
1420 COMMON /ZZ/ ZTA,ZPR
1421
1422 SAVE
1423 MASS=MASSTA+MASSPR
1424 DO 500 IRUN=1,NUM
1425 do 100 i = 1+(IRUN-1)*MASS,zta+(IRUN-1)*MASS
1426 LB(i) = 1
1427 100 continue
1428 do 200 i = zta+1+(IRUN-1)*MASS,massta+(IRUN-1)*MASS
1429 LB(i) = 2
1430 200 continue
1431 do 300 i = massta+1+(IRUN-1)*MASS,massta+zpr+(IRUN-1)*MASS
1432 LB(i) = 1
1433 300 continue
1434 do 400 i = massta+zpr+1+(IRUN-1)*MASS,
1435 1 massta+masspr+(IRUN-1)*MASS
1436 LB(i) = 2
1437 400 continue
1438 500 CONTINUE
1439 return
1440 end
1441
1442
1443 SUBROUTINE RELCOL(LCOLL,LBLOC,LCNNE,LDD,LPP,lppk,
1444 &LPN,lpd,lrho,lomega,LKN,LNNK,LDDK,LNDK,LCNND,LCNDN,
1445 &LDIRT,LDECAY,LRES,LDOU,LDDRHO,LNNRHO,LNNOM,
1446 &NT,ntmax,sp,akaon,sk)
1447
1448
1449
1450
1451
1452
1453
1454
1455
1456
1457
1458
1459
1460
1461
1462
1463
1464
1465
1466
1467
1468
1469
1470
1471
1472
1473
1474
1475
1476
1477
1478
1479
1480
1481
1482
1483
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
1502
1503
1504
1505
1506
1507
1508
1509
1510
1511
1512
1513
1514
1515
1516
1517
1518
1519
1520
1521
1522
1523
1524
1525
1526
1527
1528
1529
1530
1531
1532
1533
1534
1535
1536
1537
1538
1539
1540
1541
1542
1543
1544
1545
1546
1547 PARAMETER (MAXSTR=150001,MAXR=1,PI=3.1415926)
1548 parameter (MX=4,MY=4,MZ=8,MPX=4,MPY=4,mpz=10,mpzp=10)
1549 PARAMETER (AKA=0.498,ALA=1.1157,ASA=1.1974,aks=0.895)
1550 PARAMETER (AA1=1.26,APHI=1.02,AP1=0.13496)
1551 parameter (maxx=20,maxz=24)
1552 parameter (rrkk=0.6,prkk=0.3,srhoks=5.,ESBIN=0.04)
1553 DIMENSION MASSRN(0:MAXR),RT(3,MAXSTR),PT(3,MAXSTR),ET(MAXSTR)
1554 DIMENSION LT(MAXSTR), PROT(MAXSTR)
1555 COMMON /AA/ R(3,MAXSTR)
1556
1557 COMMON /BB/ P(3,MAXSTR)
1558
1559 COMMON /CC/ E(MAXSTR)
1560
1561 COMMON /DD/ RHO(-MAXX:MAXX,-MAXX:MAXX,-MAXZ:MAXZ),
1562 & RHOP(-MAXX:MAXX,-MAXX:MAXX,-MAXZ:MAXZ),
1563 & RHON(-MAXX:MAXX,-MAXX:MAXX,-MAXZ:MAXZ)
1564
1565 COMMON /EE/ ID(MAXSTR),LB(MAXSTR)
1566
1567 COMMON /HH/ PROPER(MAXSTR)
1568
1569 common /ff/f(-mx:mx,-my:my,-mz:mz,-mpx:mpx,-mpy:mpy,-mpz:mpzp)
1570
1571 common /gg/ dx,dy,dz,dpx,dpy,dpz
1572
1573 COMMON /INPUT/ NSTAR,NDIRCT,DIR
1574
1575 COMMON /NN/NNN
1576
1577 COMMON /RR/ MASSR(0:MAXR)
1578
1579 common /ss/ inout(20)
1580
1581 COMMON /BG/BETAX,BETAY,BETAZ,GAMMA
1582
1583 COMMON /RUN/NUM
1584
1585 COMMON /PA/RPION(3,MAXSTR,MAXR)
1586
1587 COMMON /PB/PPION(3,MAXSTR,MAXR)
1588
1589 COMMON /PC/EPION(MAXSTR,MAXR)
1590
1591 COMMON /PD/LPION(MAXSTR,MAXR)
1592
1593 COMMON /PE/PROPI(MAXSTR,MAXR)
1594
1595 COMMON /KKK/TKAON(7),EKAON(7,0:2000)
1596
1597 COMMON /KAON/ AK(3,50,36),SPECK(50,36,7),MF
1598
1599 COMMON/TABLE/ xarray(0:1000),earray(0:1000)
1600
1601 common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
1602
1603 common/leadng/lb1,px1,py1,pz1,em1,e1,xfnl,yfnl,zfnl,tfnl,
1604 1 px1n,py1n,pz1n,dp1n
1605
1606 COMMON/tdecay/tfdcy(MAXSTR),tfdpi(MAXSTR,MAXR),tft(MAXSTR)
1607
1608 common /lastt/itimeh,bimp
1609
1610
1611 COMMON/ppbmas/niso(15),nstate,ppbm(15,2),thresh(15),weight(15)
1612
1613 common/ppb1/ene,factr2(6),fsum,ppinnb,s,wtot
1614
1615 common/ppmm/pprr,ppee,pppe,rpre,xopoe,rree
1616
1617 COMMON/hbt/lblast(MAXSTR),xlast(4,MAXSTR),plast(4,MAXSTR),nlast
1618
1619 common/resdcy/NSAV,iksdcy
1620
1621 COMMON/RNDF77/NSEED
1622
1623 COMMON/FTMAX/ftsv(MAXSTR),ftsvt(MAXSTR, MAXR)
1624 dimension ftpisv(MAXSTR,MAXR),fttemp(MAXSTR)
1625 common /dpi/em2,lb2
1626 common/phidcy/iphidcy,pttrig,ntrig,maxmiss,ipi0dcy
1627
1628 DIMENSION dptemp(MAXSTR)
1629 common /para8/ idpert,npertd,idxsec
1630 COMMON /dpert/dpertt(MAXSTR,MAXR),dpertp(MAXSTR),dplast(MAXSTR),
1631 1 dpdcy(MAXSTR),dpdpi(MAXSTR,MAXR),dpt(MAXSTR, MAXR),
1632 2 dpp1(MAXSTR,MAXR),dppion(MAXSTR,MAXR)
1633
1634 real zet(-45:45)
1635 SAVE
1636 data zet /
1637 4 1.,0.,0.,0.,0.,
1638 3 1.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
1639 2 -1.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
1640 1 0.,0.,0.,-1.,0.,1.,0.,-1.,0.,-1.,
1641 s 0.,-2.,-1.,0.,1.,0.,0.,0.,0.,-1.,
1642 e 0.,
1643 s 1.,0.,-1.,0.,1.,-1.,0.,1.,2.,0.,
1644 1 1.,0.,1.,0.,-1.,0.,1.,0.,0.,0.,
1645 2 -1.,0.,1.,0.,-1.,0.,1.,0.,0.,1.,
1646 3 0.,0.,0.,0.,0.,0.,0.,0.,0.,-1.,
1647 4 0.,0.,0.,0.,-1./
1648
1649
1650
1651 call inidcy
1652
1653
1654
1655
1656
1657
1658
1659
1660
1661
1662
1663
1664
1665 RESONA=5.
1666
1667
1668 NODELT=0
1669 SUMSRT =0.
1670 LCOLL = 0
1671 LBLOC = 0
1672 LCNNE = 0
1673 LDD = 0
1674 LPP = 0
1675 lpd = 0
1676 lpdr=0
1677 lrho = 0
1678 lrhor=0
1679 lomega=0
1680 lomgar=0
1681 LPN = 0
1682 LKN = 0
1683 LNNK = 0
1684 LDDK = 0
1685 LNDK = 0
1686 lppk =0
1687 LCNND = 0
1688 LCNDN = 0
1689 LDIRT = 0
1690 LDECAY = 0
1691 LRES = 0
1692 Ldou = 0
1693 LDDRHO = 0
1694 LNNRHO = 0
1695 LNNOM = 0
1696 MSUM = 0
1697 MASSRN(0)=0
1698
1699
1700
1701 DO 1002 IL=1,5
1702 TKAON(IL)=0
1703 DO 1001 IS=1,2000
1704 EKAON(IL,IS)=0
1705 1001 CONTINUE
1706 1002 CONTINUE
1707
1708 DO 1004 i =1,NUM
1709 DO 1003 j =1,MAXSTR
1710 PROPI(j,i) = 1.
1711 1003 CONTINUE
1712 1004 CONTINUE
1713
1714 do 1102 i=1,maxstr
1715 fttemp(i)=0.
1716 do 1101 irun=1,maxr
1717 ftpisv(i,irun)=0.
1718 1101 continue
1719 1102 continue
1720
1721
1722 sp=0
1723
1724 akaon=0
1725 sk=0
1726
1727
1728
1729
1730 MASS = 0
1731
1732 DO 1000 IRUN = 1,NUM
1733 NNN=0
1734 MSUM=MSUM+MASSR(IRUN-1)
1735
1736 J10=2
1737 IF(NT.EQ.NTMAX)J10=1
1738
1739
1740
1741
1742
1743
1744
1745
1746
1747 DO 800 J1 = J10,MASSR(IRUN)
1748 I1 = J1 + MSUM
1749
1750
1751
1752 IF(E(I1).EQ.0.)GO TO 798
1753
1754
1755
1756
1757
1758 IF (LB(I1) .LT. -45 .OR. LB(I1) .GT. 45) GOTO 798
1759 X1 = R(1,I1)
1760 Y1 = R(2,I1)
1761 Z1 = R(3,I1)
1762 PX1 = P(1,I1)
1763 PY1 = P(2,I1)
1764 PZ1 = P(3,I1)
1765 EM1 = E(I1)
1766 am1= em1
1767 E1 = SQRT( EM1**2 + PX1**2 + PY1**2 + PZ1**2 )
1768 ID1 = ID(I1)
1769 LB1 = LB(I1)
1770
1771
1772 if(nt.eq.ntmax.and.(lb1.eq.21.or.lb1.eq.23)) then
1773 pk0=RANART(NSEED)
1774 if(pk0.lt.0.25) then
1775 LB(I1)=22
1776 elseif(pk0.lt.0.50) then
1777 LB(I1)=24
1778 endif
1779 LB1=LB(I1)
1780 endif
1781
1782
1783
1784
1785
1786
1787
1788 if(lb1.eq.0.or.lb1.eq.25.or.lb1.eq.26.or.lb1.eq.27
1789 & .or.lb1.eq.28.or.lb1.eq.29.or.iabs(lb1).eq.30
1790 & .or.(iabs(lb1).ge.6.and.iabs(lb1).le.13)
1791 & .or.(iksdcy.eq.1.and.lb1.eq.24)
1792 & .or.iabs(lb1).eq.16
1793 & .or.(ipi0dcy.eq.1.and.nt.eq.ntmax.and.lb1.eq.4)) then
1794
1795
1796 continue
1797 else
1798 goto 1
1799 endif
1800
1801 IF(lb1.ge.25.and.lb1.le.27) then
1802 wid=0.151
1803 ELSEIF(lb1.eq.28) then
1804 wid=0.00841
1805 ELSEIF(lb1.eq.29) then
1806 wid=0.00443
1807 ELSEIF(iabs(LB1).eq.30) then
1808 WID=0.051
1809 ELSEIF(lb1.eq.0) then
1810 wid=1.18e-6
1811
1812 ELSEIF(iksdcy.eq.1.and.lb1.eq.24) then
1813 wid=7.36e-15
1814
1815 ELSEIF(iabs(lb1).eq.16) then
1816 wid=8.87e-6
1817
1818
1819
1820 ELSEIF(LB1.EQ.32) then
1821 call WIDA1(EM1,rhomp,WID,iseed)
1822 ELSEIF(iabs(LB1).ge.6.and.iabs(LB1).le.9) then
1823 WID=WIDTH(EM1)
1824 ELSEIF((iabs(LB1).EQ.10).OR.(iabs(LB1).EQ.11)) then
1825 WID=W1440(EM1)
1826 ELSEIF((iabs(LB1).EQ.12).OR.(iabs(LB1).EQ.13)) then
1827 WID=W1535(EM1)
1828
1829 ELSEIF(ipi0dcy.eq.1.and.nt.eq.ntmax.and.lb1.eq.4) then
1830 wid=7.85e-9
1831 ENDIF
1832
1833
1834
1835 if(nt.eq.ntmax)then
1836 pdecay=1.1
1837
1838 if(iphidcy.eq.0.and.iabs(LB1).eq.29) pdecay=0.
1839
1840
1841
1842
1843 else
1844 T0=0.19733/WID
1845 GFACTR=E1/EM1
1846 T0=T0*GFACTR
1847 IF(T0.GT.0.)THEN
1848 PDECAY=1.-EXP(-DT/T0)
1849 ELSE
1850 PDECAY=0.
1851 ENDIF
1852 endif
1853 XDECAY=RANART(NSEED)
1854
1855
1856
1857
1858
1859 IF(XDECAY.LT.PDECAY) THEN
1860
1861 idecay=irun
1862 tfnl=nt*dt
1863
1864 if(nt.eq.ntmax.and.ftsv(i1).gt.((ntmax-1)*dt))
1865 1 tfnl=ftsv(i1)
1866 xfnl=x1
1867 yfnl=y1
1868 zfnl=z1
1869
1870 if(lb1.eq.0.or.lb1.eq.25.or.lb1.eq.26.or.lb1.eq.27
1871 & .or.lb1.eq.28.or.lb1.eq.29.or.iabs(lb1).eq.30
1872 & .or.(iabs(lb1).ge.6.and.iabs(lb1).le.9)
1873 & .or.(iksdcy.eq.1.and.lb1.eq.24)
1874 & .or.iabs(lb1).eq.16
1875 & .or.(ipi0dcy.eq.1.and.nt.eq.ntmax.and.lb1.eq.4)) then
1876
1877
1878
1879
1880
1881
1882
1883
1884
1885
1886
1887
1888 call resdec(i1,nt,nnn,wid,idecay,0)
1889 p(1,i1)=px1n
1890 p(2,i1)=py1n
1891 p(3,i1)=pz1n
1892
1893 dpertp(i1)=dp1n
1894
1895 if(nt.eq.ntmax) then
1896 R(1,i1)=xfnl
1897 R(2,i1)=yfnl
1898 R(3,i1)=zfnl
1899 tfdcy(i1)=tfnl
1900 endif
1901
1902
1903 if(iabs(lb1).ge.6.and.iabs(lb1).le.9) then
1904 LDECAY=LDECAY+1
1905 endif
1906
1907
1908
1909
1910
1911
1912
1913 elseif(iabs(LB1).EQ.10.OR.iabs(LB1).EQ.11) THEN
1914 NNN=NNN+1
1915 LDECAY=LDECAY+1
1916 PNSTAR=1.
1917 IF(E(I1).GT.1.22)PNSTAR=0.6
1918 IF(RANART(NSEED).LE.PNSTAR)THEN
1919
1920 CALL DECAY(idecay,I1,NNN,ISEED,wid,nt)
1921 ELSE
1922
1923 CALL DECAY2(idecay,I1,NNN,ISEED,wid,nt)
1924 NNN=NNN+1
1925 ENDIF
1926
1927 elseif(iabs(LB1).eq.12.or.iabs(LB1).eq.13) then
1928 NNN=NNN+1
1929 CALL DECAY(idecay,I1,NNN,ISEED,wid,nt)
1930 LDECAY=LDECAY+1
1931 endif
1932
1933
1934
1935
1936
1937
1938
1939
1940
1941
1942
1943
1944
1945
1946
1947 if(nt.eq.ntmax) then
1948 if(lb(i1).eq.25.or.lb(i1).eq.26.or.lb(i1).eq.27) then
1949 wid=0.151
1950 elseif(lb(i1).eq.0) then
1951 wid=1.18e-6
1952 elseif(lb(i1).eq.24.and.iksdcy.eq.1) then
1953
1954
1955 wid=7.36e-15
1956
1957 elseif(ipi0dcy.eq.1.and.lb(i1).eq.4) then
1958 wid=7.85e-9
1959 else
1960 goto 9000
1961 endif
1962 LB1=LB(I1)
1963 PX1=P(1,I1)
1964 PY1=P(2,I1)
1965 PZ1=P(3,I1)
1966 EM1=E(I1)
1967 E1=SQRT(EM1**2+PX1**2+PY1**2+PZ1**2)
1968
1969
1970 call resdec(i1,nt,nnn,wid,idecay,0)
1971 p(1,i1)=px1n
1972 p(2,i1)=py1n
1973 p(3,i1)=pz1n
1974 R(1,i1)=xfnl
1975 R(2,i1)=yfnl
1976 R(3,i1)=zfnl
1977 tfdcy(i1)=tfnl
1978
1979 dpertp(i1)=dp1n
1980 endif
1981
1982
1983 if(nt.eq.ntmax.and.ipi0dcy.eq.1.and.lb(i1).eq.4) then
1984 wid=7.85e-9
1985 LB1=LB(I1)
1986 PX1=P(1,I1)
1987 PY1=P(2,I1)
1988 PZ1=P(3,I1)
1989 EM1=E(I1)
1990 E1=SQRT(EM1**2+PX1**2+PY1**2+PZ1**2)
1991 call resdec(i1,nt,nnn,wid,idecay,0)
1992 p(1,i1)=px1n
1993 p(2,i1)=py1n
1994 p(3,i1)=pz1n
1995 R(1,i1)=xfnl
1996 R(2,i1)=yfnl
1997 R(3,i1)=zfnl
1998 tfdcy(i1)=tfnl
1999 dpertp(i1)=dp1n
2000 endif
2001
2002
2003
2004
2005 9000 go to 798
2006
2007 ENDIF
2008
2009
2010
2011
2012 1 if(nt.eq.ntmax)go to 798
2013
2014 X1 = R(1,I1)
2015 Y1 = R(2,I1)
2016 Z1 = R(3,I1)
2017
2018 DO 600 J2 = 1,J1-1
2019 I2 = J2 + MSUM
2020
2021 IF(E(I2).EQ.0.) GO TO 600
2022
2023 IF(E(I1).EQ.0.) GO TO 800
2024
2025 IF (LB(I2) .LT. -45 .OR. LB(I2) .GT. 45) GOTO 600
2026
2027 X2=R(1,I2)
2028 Y2=R(2,I2)
2029 Z2=R(3,I2)
2030 dr0max=5.
2031
2032 ilb1=iabs(LB(I1))
2033 ilb2=iabs(LB(I2))
2034 IF(ilb1.EQ.42.or.ilb2.EQ.42) THEN
2035 if((ILB1.GE.1.AND.ILB1.LE.2)
2036 1 .or.(ILB1.GE.6.AND.ILB1.LE.13)
2037 2 .or.(ILB2.GE.1.AND.ILB2.LE.2)
2038 3 .or.(ILB2.GE.6.AND.ILB2.LE.13)) then
2039 if((lb(i1)*lb(i2)).gt.0) dr0max=10.
2040 endif
2041 ENDIF
2042
2043 if(((X1-X2)**2+(Y1-Y2)**2+(Z1-Z2)**2).GT.dr0max**2)
2044 1 GO TO 600
2045 IF (ID(I1)*ID(I2).EQ.IAVOID) GOTO 400
2046 ID1=ID(I1)
2047 ID2 = ID(I2)
2048
2049 ix1= nint(x1/dx)
2050 iy1= nint(y1/dy)
2051 iz1= nint(z1/dz)
2052 PX1=P(1,I1)
2053 PY1=P(2,I1)
2054 PZ1=P(3,I1)
2055 EM1=E(I1)
2056 AM1=EM1
2057 LB1=LB(I1)
2058 E1=SQRT(EM1**2+PX1**2+PY1**2+PZ1**2)
2059 IPX1=NINT(PX1/DPX)
2060 IPY1=NINT(PY1/DPY)
2061 IPZ1=NINT(PZ1/DPZ)
2062 LB2 = LB(I2)
2063 PX2 = P(1,I2)
2064 PY2 = P(2,I2)
2065 PZ2 = P(3,I2)
2066 EM2=E(I2)
2067 AM2=EM2
2068 lb1i=lb(i1)
2069 lb2i=lb(i2)
2070 px1i=P(1,I1)
2071 py1i=P(2,I1)
2072 pz1i=P(3,I1)
2073 em1i=E(I1)
2074 px2i=P(1,I2)
2075 py2i=P(2,I2)
2076 pz2i=P(3,I2)
2077 em2i=E(I2)
2078
2079 eini=SQRT(E(I1)**2+P(1,I1)**2+P(2,I1)**2+P(3,I1)**2)
2080 1 +SQRT(E(I2)**2+P(1,I2)**2+P(2,I2)**2+P(3,I2)**2)
2081 pxini=P(1,I1)+P(1,I2)
2082 pyini=P(2,I1)+P(2,I2)
2083 pzini=P(3,I1)+P(3,I2)
2084 nnnini=nnn
2085
2086
2087 iblock=0
2088
2089
2090
2091
2092
2093
2094
2095
2096
2097 DELTR0=3.
2098 if( (iabs(lb1).ge.14.and.iabs(lb1).le.17) .or.
2099 & (iabs(lb1).ge.30.and.iabs(lb1).le.45) ) DELTR0=5.0
2100 if( (iabs(lb2).ge.14.and.iabs(lb2).le.17) .or.
2101 & (iabs(lb2).ge.30.and.iabs(lb2).le.45) ) DELTR0=5.0
2102
2103 if(lb1.eq.28.and.lb2.eq.28) DELTR0=4.84
2104
2105 if((lb1.ge.3.and.lb1.le.5).and.(lb2.ge.3.and.lb2.le.5)) then
2106 E2=SQRT(EM2**2+PX2**2+PY2**2+PZ2**2)
2107 spipi=(e1+e2)**2-(px1+px2)**2-(py1+py2)**2-(pz1+pz2)**2
2108 if(spipi.ge.(4*0.77**2)) DELTR0=3.5
2109 endif
2110
2111
2112 IF (LB1.EQ.23 .AND. (LB2.GE.14.AND.LB2.LE.17)) GOTO 3699
2113 IF (LB2.EQ.23 .AND. (LB1.GE.14.AND.LB1.LE.17)) GOTO 3699
2114
2115
2116
2117 if(lb1.eq.21.and.lb2.eq.23)go to 3699
2118 if(lb2.eq.21.and.lb1.eq.23)go to 3699
2119 if(lb1.eq.30.and.lb2.eq.21)go to 3699
2120 if(lb2.eq.30.and.lb1.eq.21)go to 3699
2121 if(lb1.eq.-30.and.lb2.eq.23)go to 3699
2122 if(lb2.eq.-30.and.lb1.eq.23)go to 3699
2123 if(lb1.eq.-30.and.lb2.eq.30)go to 3699
2124 if(lb2.eq.-30.and.lb1.eq.30)go to 3699
2125
2126
2127
2128 if(lb1.eq.21.or.lb1.eq.23) then
2129 if(lb2.eq.0.or.(lb2.ge.25.and.lb2.le.28)) then
2130 go to 3699
2131 endif
2132 elseif(lb2.eq.21.or.lb2.eq.23) then
2133 if(lb1.eq.0.or.(lb1.ge.25.and.lb1.le.28)) then
2134 goto 3699
2135 endif
2136 endif
2137
2138
2139 if(iabs(lb1).eq.30 .and.
2140 1 (lb2.eq.0.or.(lb2.ge.25.and.lb2.le.28)
2141 2 .or.(lb2.ge.3.and.lb2.le.5))) then
2142 go to 3699
2143 elseif(iabs(lb2).eq.30 .and.
2144 1 (lb1.eq.0.or.(lb1.ge.25.and.lb1.le.28)
2145 2 .or.(lb1.ge.3.and.lb1.le.5))) then
2146 goto 3699
2147
2148
2149 elseif( iabs(lb1).eq.30 .and.
2150 1 (iabs(lb2).eq.1.or.iabs(lb2).eq.2.or.
2151 2 (iabs(lb2).ge.6.and.iabs(lb2).le.13)) )then
2152 go to 3699
2153 endif
2154 if( iabs(lb2).eq.30 .and.
2155 1 (iabs(lb1).eq.1.or.iabs(lb1).eq.2.or.
2156 2 (iabs(lb1).ge.6.and.iabs(lb1).le.13)) )then
2157 go to 3699
2158 endif
2159
2160
2161
2162
2163 if((lb1.eq.23.or.lb1.eq.21).and.
2164 1 (iabs(lb2).eq.1.or.iabs(lb2).eq.2.or.
2165 2 (iabs(lb2).ge.6.and.iabs(lb2).le.13))) then
2166 go to 3699
2167 elseif((lb2.eq.23.or.lb2.eq.21).and.
2168 1 (iabs(lb1).eq.1.or.iabs(lb1).eq.2.or.
2169 2 (iabs(lb1).ge.6.and.iabs(lb1).le.13))) then
2170 go to 3699
2171 endif
2172
2173
2174
2175
2176
2177
2178
2179
2180
2181
2182
2183
2184
2185
2186 rppmax=3.57
2187
2188 if((lb1.eq.-1.or.lb1.eq.-2.or.(lb1.ge.-13.and.lb1.le.-6))
2189 1 .and.(lb2.eq.1.or.lb2.eq.2.or.(lb2.ge.6.and.lb2.le.13))) then
2190 DELTR0 = RPPMAX
2191 GOTO 2699
2192 else if((lb2.eq.-1.or.lb2.eq.-2.or.(lb2.ge.-13.and.lb2.le.-6))
2193 1 .and.(lb1.eq.1.or.lb1.eq.2.or.(lb1.ge.6.and.lb1.le.13))) then
2194 DELTR0 = RPPMAX
2195 GOTO 2699
2196 END IF
2197
2198
2199 if( (iabs(lb1).ge.14.and.iabs(lb1).le.17) .or.
2200 & (iabs(lb2).ge.14.and.iabs(lb2).le.17) )go to 3699
2201
2202
2203 IF (iabs(LB1).EQ.42.or.iabs(LB2).EQ.42) THEN
2204 ilb1=iabs(LB1)
2205 ilb2=iabs(LB2)
2206 if((ILB1.GE.1.AND.ILB1.LE.2)
2207 1 .or.(ILB1.GE.6.AND.ILB1.LE.13)
2208 2 .or.(ILB2.GE.1.AND.ILB2.LE.2)
2209 3 .or.(ILB2.GE.6.AND.ILB2.LE.13)) then
2210 if((lb1*lb2).gt.0) deltr0=9.5
2211 endif
2212 ENDIF
2213
2214 if( (iabs(lb1).ge.40.and.iabs(lb1).le.45) .or.
2215 & (iabs(lb2).ge.40.and.iabs(lb2).le.45) )go to 3699
2216
2217
2218 IF( (lb1.eq.29 .and.((lb2.ge.1.and.lb2.le.13).or.
2219 & (lb2.ge.21.and.lb2.le.28).or.iabs(lb2).eq.30)) .OR.
2220 & (lb2.eq.29 .and.((lb1.ge.1.and.lb1.le.13).or.
2221 & (lb1.ge.21.and.lb1.le.28).or.iabs(lb1).eq.30)) )THEN
2222 DELTR0=3.0
2223 go to 3699
2224 endif
2225
2226
2227
2228
2229
2230 If(iabs(lb1).eq.30.or.iabs(lb2).eq.30) go to 400
2231
2232 If(lb1.eq.23.and.(lb2.lt.1.or.lb2.gt.17))go to 400
2233 If(lb2.eq.23.and.(lb1.lt.1.or.lb1.gt.17))go to 400
2234
2235
2236
2237 if( ((lb1.le.-1.and.lb1.ge.-13)
2238 & .and.(lb2.eq.0.or.(lb2.ge.3.and.lb2.le.5)
2239 & .or.(lb2.ge.25.and.lb2.le.28)))
2240 & .OR.((lb2.le.-1.and.lb2.ge.-13)
2241 & .and.(lb1.eq.0.or.(lb1.ge.3.and.lb1.le.5)
2242 & .or.(lb1.ge.25.and.lb1.le.28))) ) then
2243 elseIF( ((LB1.eq.-1.or.lb1.eq.-2).
2244 & and.(LB2.LT.-5.and.lb2.ge.-13))
2245 & .OR. ((LB2.eq.-1.or.lb2.eq.-2).
2246 & and.(LB1.LT.-5.and.lb1.ge.-13)) )then
2247 elseIF((LB1.eq.-1.or.lb1.eq.-2)
2248 & .AND.(LB2.eq.-1.or.lb2.eq.-2))then
2249 elseIF((LB1.LT.-5.and.lb1.ge.-13).AND.
2250 & (LB2.LT.-5.and.lb2.ge.-13)) then
2251
2252
2253 endif
2254
2255 2699 CONTINUE
2256
2257 IF (LB1 .EQ. 1 .OR. LB1 .EQ. 2 .OR. (LB1 .GE. 6 .AND.
2258 & LB1 .LE. 17)) THEN
2259 IF (LB2 .EQ. 1 .OR. LB2 .EQ. 2 .OR. (LB2 .GE. 6 .AND.
2260 & LB2 .LE. 17)) THEN
2261 DELTR0 = 2.
2262 END IF
2263 END IF
2264
2265 3699 RSQARE = (X1-X2)**2 + (Y1-Y2)**2 + (Z1-Z2)**2
2266 IF (RSQARE .GT. DELTR0**2) GO TO 400
2267
2268
2269 ix2 = nint(x2/dx)
2270 iy2 = nint(y2/dy)
2271 iz2 = nint(z2/dz)
2272 ipx2 = nint(px2/dpx)
2273 ipy2 = nint(py2/dpy)
2274 ipz2 = nint(pz2/dpz)
2275
2276
2277 CALL CMS(I1,I2,PCX,PCY,PCZ,SRT)
2278
2279
2280 drmax=dr0max
2281 call distc0(drmax,deltr0,DT,
2282 1 Ifirst,PCX,PCY,PCZ,
2283 2 x1,y1,z1,px1,py1,pz1,em1,x2,y2,z2,px2,py2,pz2,em2)
2284 if(Ifirst.eq.-1) goto 400
2285
2286 ISS=NINT(SRT/ESBIN)
2287
2288 if(ISS.gt.2000) ISS=2000
2289
2290
2291
2292
2293 IF (iabs(LB1).EQ.42.or.iabs(LB2).EQ.42) THEN
2294 ilb1=iabs(LB1)
2295 ilb2=iabs(LB2)
2296 if(LB1.eq.0.or.(LB1.GE.3.AND.LB1.LE.5)
2297 1 .or.(LB1.GE.25.AND.LB1.LE.28)
2298 2 .or.
2299 3 LB2.eq.0.or.(LB2.GE.3.AND.LB2.LE.5)
2300 4 .or.(LB2.GE.25.AND.LB2.LE.28)) then
2301 GOTO 505
2302
2303 elseif(((ILB1.GE.1.AND.ILB1.LE.2)
2304 1 .or.(ILB1.GE.6.AND.ILB1.LE.13)
2305 2 .or.(ILB2.GE.1.AND.ILB2.LE.2)
2306 3 .or.(ILB2.GE.6.AND.ILB2.LE.13))
2307 4 .and.(lb1*lb2).gt.0) then
2308 GOTO 506
2309 else
2310 GOTO 400
2311 endif
2312 ENDIF
2313
2314
2315 if( ((lb1.eq.23.or.lb1.eq.30).and.
2316 & (lb2.eq.-1.or.lb2.eq.-2.or.(lb2.ge.-13.and.lb2.le.-6)))
2317 & .OR.((lb2.eq.23.or.lb2.eq.30).and.
2318 & (lb1.eq.-1.or.lb1.eq.-2.or.(lb1.ge.-13.and.lb1.le.-6))) )
2319 & then
2320 bmass=0.938
2321 if(srt.le.(bmass+aka)) then
2322 pkaon=0.
2323 else
2324 pkaon=sqrt(((srt**2-(aka**2+bmass**2))
2325 1 /2./bmass)**2-aka**2)
2326 endif
2327
2328
2329 sigela = 0.5 * (AKPEL(PKAON) + AKNEL(PKAON))
2330 SIGSGM = 1.5 * AKPSGM(PKAON) + AKNSGM(PKAON)
2331 SIG = sigela + SIGSGM + AKPLAM(PKAON)
2332 if(sig.gt.1.e-7) then
2333
2334 icase=3
2335 brel=sigela/sig
2336 brsgm=sigsgm/sig
2337 brsig = sig
2338 nchrg = 1
2339 go to 3555
2340 endif
2341 go to 400
2342 endif
2343
2344
2345
2346 if(((lb1.ge.-17.and.lb1.le.-14).and.(lb2.ge.3.and.lb2.le.5))
2347 & .OR.((lb2.ge.-17.and.lb2.le.-14)
2348 & .and.(lb1.ge.3.and.lb1.le.5)))then
2349 nchrg=-100
2350
2351
2352 if((lb1.eq.-15.and.(lb2.eq.5.or.lb2.eq.27)).OR.
2353 & (lb2.eq.-15.and.(lb1.eq.5.or.lb1.eq.27))) then
2354 nchrg=-2
2355
2356 bmass=1.232
2357 go to 110
2358 endif
2359 if( (lb1.eq.-15.and.(lb2.eq.0.or.lb2.eq.4.or.lb2.eq.26.or.
2360 & lb2.eq.28)).OR.(lb2.eq.-15.and.(lb1.eq.0.or.
2361 & lb1.eq.4.or.lb1.eq.26.or.lb1.eq.28)).OR.
2362 & ((lb1.eq.-14.or.lb1.eq.-16).and.(lb2.eq.5.or.lb2.eq.27)).OR.
2363 & ((lb2.eq.-14.or.lb2.eq.-16).and.(lb1.eq.5.or.lb1.eq.27)) )then
2364 nchrg=-1
2365
2366 bmass=0.938
2367 go to 110
2368 endif
2369 if( (lb1.eq.-15.and.(lb2.eq.3.or.lb2.eq.25)).OR.
2370 & (lb2.eq.-15.and.(lb1.eq.3.or.lb1.eq.25)).OR.
2371 & (lb1.eq.-17.and.(lb2.eq.5.or.lb2.eq.27)).OR.
2372 & (lb2.eq.-17.and.(lb1.eq.5.or.lb1.eq.27)).OR.
2373 & ((lb1.eq.-14.or.lb1.eq.-16).and.(lb2.eq.0.or.lb2.eq.4
2374 & .or.lb2.eq.26.or.lb2.eq.28)).OR.
2375 & ((lb2.eq.-14.or.lb2.eq.-16).and.(lb1.eq.0.or.lb1.eq.4
2376 & .or.lb1.eq.26.or.lb1.eq.28)) )then
2377 nchrg=0
2378
2379 bmass=0.938
2380 go to 110
2381 endif
2382 if( (lb1.eq.-17.and.(lb2.eq.0.or.lb2.eq.4.or.lb2.eq.26.or.
2383 & lb2.eq.28)).OR.(lb2.eq.-17.and.(lb1.eq.0.or.
2384 & lb1.eq.4.or.lb1.eq.26.or.lb1.eq.28)).OR.
2385 & ((lb1.eq.-14.or.lb1.eq.-16).and.(lb2.eq.3.or.lb2.eq.25)).OR.
2386 & ((lb2.eq.-14.or.lb2.eq.-16).and.(lb1.eq.3.or.lb1.eq.25)))then
2387 nchrg=1
2388
2389 bmass=1.232
2390 endif
2391
2392
2393 110 sig = 0.
2394
2395 if(nchrg.ne.-100.and.srt.ge.(aka+bmass))then
2396
2397
2398 icase=4
2399
2400 pkaon=sqrt(((srt**2-(aka**2+0.938**2))/2./0.938)**2-aka**2)
2401
2402 if(lb1.eq.-14.or.lb2.eq.-14) then
2403 if(nchrg.ge.0) sigma0=akPlam(pkaon)
2404 if(nchrg.lt.0) sigma0=akNlam(pkaon)
2405
2406 else
2407
2408 if(nchrg.ge.0) sigma0=akPsgm(pkaon)
2409
2410 if(nchrg.lt.0) sigma0=akNsgm(pkaon)
2411 SIGMA0 = 1.5 * AKPSGM(PKAON) + AKNSGM(PKAON)
2412 endif
2413 sig=(srt**2-(aka+bmass)**2)*(srt**2-(aka-bmass)**2)/
2414 & (srt**2-(em1+em2)**2)/(srt**2-(em1-em2)**2)*sigma0
2415
2416 if(nchrg.eq.-2.or.nchrg.eq.2) sig=2.*sig
2417
2418
2419 IF (LB1 .EQ. -14 .OR. LB2 .EQ. -14) THEN
2420 SIG = 4.0 / 3.0 * SIG
2421 ELSE IF (NCHRG .EQ. -2 .OR. NCHRG .EQ. 2) THEN
2422 SIG = 8.0 / 9.0 * SIG
2423 ELSE
2424 SIG = 4.0 / 9.0 * SIG
2425 END IF
2426
2427
2428
2429
2430
2431 endif
2432
2433 icase=4
2434 sigela = 10.
2435 sig = sig + sigela
2436 brel= sigela/sig
2437 brsgm=0.
2438 brsig = sig
2439
2440 go to 3555
2441 endif
2442
2443
2444
2445
2446 if( ((lb1.eq.21.or.lb1.eq.-30).and.(lb2.ge.14.and.lb2.le.17)).OR.
2447 & ((lb2.eq.21.or.lb2.eq.-30).and.(lb1.ge.14.and.lb1.le.17)) )then
2448 kp = 0
2449 go to 3455
2450 endif
2451
2452 if( ((lb1.eq.23.or.lb1.eq.30).and.(lb2.le.-14.and.lb2.ge.-17)).OR.
2453 & ((lb2.eq.23.or.lb2.eq.30).and.(lb1.le.-14.and.lb1.ge.-17)) )then
2454 kp = 1
2455 go to 3455
2456 endif
2457
2458 if( ((lb1.eq.21.or.lb1.eq.-30).and.(lb2.eq.40.or.lb2.eq.41)).OR.
2459 & ((lb2.eq.21.or.lb2.eq.-30).and.(lb1.eq.40.or.lb1.eq.41)) )then
2460 kp = 0
2461 go to 3455
2462 endif
2463
2464 if( ((lb1.eq.23.or.lb1.eq.30).and.(lb2.eq.-40.or.lb2.eq.-41)).OR.
2465 & ((lb2.eq.23.or.lb2.eq.30).and.(lb1.eq.-40.or.lb1.eq.-41)) )then
2466 kp = 1
2467 go to 3455
2468 endif
2469
2470
2471
2472
2473 kp = 3
2474
2475 if( (((lb1.ge.3.and.lb1.le.5).or.lb1.eq.0)
2476 & .and.(iabs(lb2).eq.40.or.iabs(lb2).eq.41))
2477 & .OR. (((lb2.ge.3.and.lb2.le.5).or.lb2.eq.0)
2478 & .and.(iabs(lb1).eq.40.or.iabs(lb1).eq.41)) )go to 3455
2479
2480
2481
2482 if( ((lb1.ge.3.and.lb1.le.5).and.iabs(lb2).eq.45)
2483 & .OR.((lb2.ge.3.and.lb2.le.5).and.iabs(lb1).eq.45) )go to 3455
2484
2485
2486
2487
2488
2489 IF (LB1.EQ.23 .AND. (LB2.GE.14.AND.LB2.LE.17)) GOTO 5699
2490 IF (LB2.EQ.23 .AND. (LB1.GE.14.AND.LB1.LE.17)) GOTO 5699
2491
2492 IF (LB1.EQ.21 .AND. (LB2.GE.-17.AND.LB2.LE.-14)) GOTO 5699
2493 IF (LB2.EQ.21 .AND. (LB1.GE.-17.AND.LB1.LE.-14)) GOTO 5699
2494
2495
2496 IF( (((LB1.eq.1.or.LB1.eq.2).or.(LB1.ge.6.and.LB1.le.13))
2497 & .AND.(LB2.GE.-17.AND.LB2.LE.-14)) .OR.
2498 & (((LB2.eq.1.or.LB2.eq.2).or.(LB2.ge.6.and.LB2.le.13))
2499 & .AND.(LB1.GE.-17.AND.LB1.LE.-14)) )go to 5999
2500
2501 IF( (((LB1.eq.-1.or.LB1.eq.-2).or.(LB1.le.-6.and.LB1.ge.-13))
2502 & .AND.(LB2.GE.14.AND.LB2.LE.17)) .OR.
2503 & (((LB2.eq.-1.or.LB2.eq.-2).or.(LB2.le.-6.and.LB2.ge.-13))
2504 & .AND.(LB1.GE.14.AND.LB1.LE.17)) )go to 5999
2505
2506
2507
2508 if(lb1.eq.21.and.lb2.eq.23) go to 8699
2509 if(lb2.eq.21.and.lb1.eq.23) go to 8699
2510 if(lb1.eq.30.and.lb2.eq.21) go to 8699
2511 if(lb2.eq.30.and.lb1.eq.21) go to 8699
2512 if(lb1.eq.-30.and.lb2.eq.23) go to 8699
2513 if(lb2.eq.-30.and.lb1.eq.23) go to 8699
2514 if(lb1.eq.-30.and.lb2.eq.30) go to 8699
2515 if(lb2.eq.-30.and.lb1.eq.30) go to 8699
2516
2517 IF( ((lb1.eq.23.or.lb1.eq.21.or.iabs(lb1).eq.30) .and.
2518 & (lb2.ge.25.and.lb2.le.28)) .OR.
2519 & ((lb2.eq.23.or.lb2.eq.21.or.iabs(lb2).eq.30) .and.
2520 & (lb1.ge.25.and.lb1.le.28)) ) go to 8799
2521
2522
2523 IF( (iabs(lb1).eq.30.and.(lb2.ge.3.and.lb2.le.5)) .OR.
2524 & (iabs(lb2).eq.30.and.(lb1.ge.3.and.lb1.le.5)) )go to 8799
2525
2526
2527
2528
2529 IF( (lb1.eq.29 .and.(lb2.eq.1.or.lb2.eq.2.or.
2530 & (lb2.ge.6.and.lb2.le.9))) .OR.
2531 & (lb2.eq.29 .and.(lb1.eq.1.or.lb1.eq.2.or.
2532 & (lb1.ge.6.and.lb1.le.9))) )go to 7222
2533
2534
2535 IF( (lb1.eq.29 .and.((lb2.ge.3.and.lb2.le.5).or.
2536 & (lb2.ge.21.and.lb2.le.28).or.iabs(lb2).eq.30)) .OR.
2537 & (lb2.eq.29 .and.((lb1.ge.3.and.lb1.le.5).or.
2538 & (lb1.ge.21.and.lb1.le.28).or.iabs(lb1).eq.30)) )THEN
2539 go to 7444
2540 endif
2541
2542
2543
2544
2545 if( ((iabs(lb1).ge.14.and.iabs(lb1).le.17).or.iabs(lb1).ge.40)
2546 & .and.((lb2.ge.25.and.lb2.le.29).or.lb2.eq.0) )go to 888
2547 if( ((iabs(lb2).ge.14.and.iabs(lb2).le.17).or.iabs(lb2).ge.40)
2548 & .and.((lb1.ge.25.and.lb1.le.29).or.lb1.eq.0) )go to 888
2549
2550
2551 if( ((lb1.eq.23.or.lb1.eq.30).and.(lb2.eq.1.or.lb2.eq.2.or.
2552 & (lb2.ge.6.and.lb2.le.13))) .OR.
2553 & ((lb2.eq.23.or.lb2.eq.30).and.(lb1.eq.1.or.lb1.eq.2.or.
2554 & (lb1.ge.6.and.lb1.le.13))) ) go to 888
2555 if( ((lb1.eq.21.or.lb1.eq.-30).and.(lb2.eq.-1.or.lb2.eq.-2.or.
2556 & (lb2.ge.-13.and.lb2.le.-6))) .OR.
2557 & ((lb2.eq.21.or.lb2.eq.-30).and.(lb1.eq.-1.or.lb1.eq.-2.or.
2558 & (lb1.ge.-13.and.lb1.le.-6))) ) go to 888
2559
2560
2561 If( ((lb1.ge.14.and.lb1.le.17).and.(lb2.ge.6.and.lb2.le.13))
2562 & .OR.((lb2.ge.14.and.lb2.le.17).and.(lb1.ge.6.and.lb1.le.13)) )
2563 & go to 7799
2564 If(((lb1.le.-14.and.lb1.ge.-17).and.(lb2.le.-6.and.lb2.ge.-13))
2565 &.OR.((lb2.le.-14.and.lb2.ge.-17).and.(lb1.le.-6.and.lb1.ge.-13)))
2566 & go to 7799
2567
2568
2569 if( iabs(lb1).ge.40 .or. iabs(lb2).ge.40
2570 & .or. (lb1.le.-14.and.lb1.ge.-17)
2571 & .or. (lb2.le.-14.and.lb2.ge.-17) )go to 400
2572
2573
2574
2575 if((lb1.eq.-1.or.lb1.eq.-2.or.(lb1.ge.-13.and.lb1.le.-6))
2576 1 .and.(lb2.eq.1.or.lb2.eq.2.or.(lb2.ge.6.and.lb2.le.13))) then
2577 GOTO 2799
2578 else if((lb2.eq.-1.or.lb2.eq.-2.or.(lb2.ge.-13.and.lb2.le.-6))
2579 1 .and.(lb1.eq.1.or.lb1.eq.2.or.(lb1.ge.6.and.lb1.le.13))) then
2580 GOTO 2799
2581 END IF
2582
2583
2584 inewka=irun
2585
2586
2587
2588
2589 call newka(icase,inewka,iseed,dt,nt,
2590 & ictrl,i1,i2,srt,pcx,pcy,pcz,iblock)
2591
2592
2593 IF (ICTRL .EQ. 1) GOTO 400
2594
2595
2596
2597
2598
2599 if((iabs(lb1).ge.14.and.iabs(lb1).le.17).
2600 & or.(iabs(lb2).ge.14.and.iabs(lb2).le.17))go to 400
2601
2602
2603 IF((lb1.ge.3.and.lb1.le.5).and.(lb2.ge.3.and.lb2.le.5))GO TO 777
2604 if(lb1.eq.0.and.(lb2.ge.3.and.lb2.le.5)) go to 777
2605 if(lb2.eq.0.and.(lb1.ge.3.and.lb1.le.5)) go to 777
2606 if(lb1.eq.0.and.lb2.eq.0)go to 777
2607
2608
2609
2610 if( (lb1.ge.25.and.lb1.le.28).and.
2611 & (lb2.ge.25.and.lb2.le.28) )goto 777
2612
2613 If((lb1.ge.25.and.lb1.le.28).and.(lb2.ge.3.and.lb2.le.5))go to 777
2614 If((lb2.ge.25.and.lb2.le.28).and.(lb1.ge.3.and.lb1.le.5))go to 777
2615
2616 if((lb1.ge.25.and.lb1.le.28).and.lb2.eq.0)go to 777
2617 if((lb2.ge.25.and.lb2.le.28).and.lb1.eq.0)go to 777
2618
2619
2620 if((lb1.eq.23.or.lb1.eq.21).and.(lb2.ge.3.and.lb2.le.5))go to 889
2621 if((lb2.eq.23.or.lb2.eq.21).and.(lb1.ge.3.and.lb1.le.5))go to 889
2622
2623
2624
2625 If(iabs(lb1).eq.30.or.iabs(lb2).eq.30) go to 400
2626 If(lb1.eq.21.or.lb2.eq.21) go to 400
2627 If(lb1.eq.23.or.lb2.eq.23) go to 400
2628
2629
2630 IF( (LB1.ge.3.and.LB1.le.5) .and.
2631 & (iabs(LB2).eq.1.or.iabs(LB2).eq.2.or.
2632 & (iabs(LB2).ge.6.and.iabs(LB2).le.13)) )GO TO 3
2633 IF( (LB2.ge.3.and.LB2.le.5) .and.
2634 & (iabs(LB1).eq.1.or.iabs(LB1).eq.2.or.
2635 & (iabs(LB1).ge.6.and.iabs(LB1).le.13)) )GO TO 3
2636
2637
2638 IF( (LB1.ge.25.and.LB1.le.28) .and.
2639 & (iabs(LB2).eq.1.or.iabs(LB2).eq.2.or.
2640 & (iabs(LB2).ge.6.and.iabs(LB2).le.13)) )GO TO 33
2641 IF( (LB2.ge.25.and.LB2.le.28) .and.
2642 & (iabs(LB1).eq.1.or.iabs(LB1).eq.2.or.
2643 & (iabs(LB1).ge.6.and.iabs(LB1).le.13)) )GO TO 33
2644
2645
2646 IF( LB1.eq.0 .and.
2647 & (iabs(LB2).eq.1.or.iabs(LB2).eq.2.or.
2648 & (iabs(LB2).ge.6.and.iabs(LB2).le.13)) )GO TO 547
2649 IF( LB2.eq.0 .and.
2650 & (iabs(LB1).eq.1.or.iabs(LB1).eq.2.or.
2651 & (iabs(LB1).ge.6.and.iabs(LB1).le.13)) )GO TO 547
2652
2653
2654 IF((LB1.eq.1.or.lb1.eq.2).
2655 & AND.(LB2.GT.5.and.lb2.le.13))GOTO 44
2656 IF((LB2.eq.1.or.lb2.eq.2).
2657 & AND.(LB1.GT.5.and.lb1.le.13))GOTO 44
2658 IF((LB1.eq.-1.or.lb1.eq.-2).
2659 & AND.(LB2.LT.-5.and.lb2.ge.-13))GOTO 44
2660 IF((LB2.eq.-1.or.lb2.eq.-2).
2661 & AND.(LB1.LT.-5.and.lb1.ge.-13))GOTO 44
2662
2663
2664 IF((LB1.eq.1.or.lb1.eq.2).AND.(LB2.eq.1.or.lb2.eq.2))GOTO 4
2665 IF((LB1.eq.-1.or.lb1.eq.-2).AND.(LB2.eq.-1.or.lb2.eq.-2))GOTO 4
2666
2667
2668 IF((LB1.GT.5.and.lb1.le.13).AND.
2669 & (LB2.GT.5.and.lb2.le.13)) GOTO 444
2670 IF((LB1.LT.-5.and.lb1.ge.-13).AND.
2671 & (LB2.LT.-5.and.lb2.ge.-13)) GOTO 444
2672
2673
2674
2675 if((lb1.lt.3).and.(lb2.ge.14.and.lb2.le.17))goto 400
2676 if((lb2.lt.3).and.(lb1.ge.14.and.lb1.le.17))goto 400
2677 if((lb1.ge.14.and.lb1.le.17).and.
2678 & (lb2.ge.14.and.lb2.le.17))goto 400
2679
2680
2681 go to 400
2682
2683
2684 547 IF(LB1*LB2.EQ.0)THEN
2685
2686
2687
2688
2689
2690 ece=(em1+em2+0.02)**2
2691 xkaon0=0.
2692 if(srt.ge.1.63.AND.SRT.LE.1.7)xkaon0=pnlka(srt)
2693 IF(SRT.GT.1.7)XKAON0=PNLKA(SRT)+pnska(srt)
2694
2695 XKAON0 = 2.0 * XKAON0
2696
2697
2698
2699
2700
2701 xkaon=xkaon0
2702
2703 XETA=XN1535(I1,I2,0)
2704 If((iabs(LB(I1)).ge.6.and.iabs(LB(I1)).le.13).or.
2705 & (iabs(LB(I2)).ge.6.and.iabs(LB(I2)).le.13)) xeta=0.
2706 IF((XETA+xkaon).LE.1.e-06)GO TO 400
2707 DSE=SQRT((XETA+XKAON)/PI)
2708 DELTRE=DSE+0.1
2709 px1cm=pcx
2710 py1cm=pcy
2711 pz1cm=pcz
2712
2713 CALL DISTCE(I1,I2,DELTRE,DSE,DT,ECE,SRT,IC,
2714 1 PCX,PCY,PCZ)
2715 IF(IC.EQ.-1) GO TO 400
2716 ekaon(4,iss)=ekaon(4,iss)+1
2717 IF(XKAON0/(XKAON+XETA).GT.RANART(NSEED))then
2718
2719 CALL CREN(PX1CM,PY1CM,PZ1CM,SRT,I1,I2,IBLOCK)
2720
2721 IF(IBLOCK.EQ.7) then
2722 LPN=LPN+1
2723 elseIF(IBLOCK.EQ.-7) then
2724 endif
2725
2726 em1=e(i1)
2727 em2=e(i2)
2728 GO TO 440
2729 endif
2730
2731 resona=1.
2732 GO TO 98
2733 ENDIF
2734
2735 3 CONTINUE
2736 px1cm=pcx
2737 py1cm=pcy
2738 pz1cm=pcz
2739
2740
2741 xkaon0=0.
2742 if(srt.ge.1.63.AND.SRT.LE.1.7)xkaon0=pnlka(srt)
2743 IF(SRT.GT.1.7)XKAON0=PNLKA(SRT)+pnska(srt)
2744 XKAON0 = 2.0 * XKAON0
2745
2746
2747 Xphi = 0.
2748 if( ( ((lb1.ge.1.and.lb1.le.2).or.
2749 & (lb1.ge.6.and.lb1.le.9))
2750 & .OR.((lb2.ge.1.and.lb2.le.2).or.
2751 & (lb2.ge.6.and.lb2.le.9)) )
2752 & .AND. srt.gt.1.958)
2753 & call pibphi(srt,lb1,lb2,em1,em2,Xphi,xphin)
2754
2755
2756
2757
2758
2759
2760
2761
2762 If((iabs(LB(I1)).ge.6.and.iabs(LB(I1)).le.13).or.
2763 & (iabs(LB(I2)).ge.6.and.iabs(LB(I2)).le.13)) go to 31
2764
2765
2766
2767
2768
2769
2770
2771
2772
2773
2774
2775
2776
2777
2778 EC=(em1+em2+0.02)**2
2779 xkaon=0.
2780 if(srt.gt.1.23)xkaon=(pionpp(srt)+PIPP1(SRT))/2.
2781
2782
2783
2784
2785
2786
2787
2788
2789 IF( (LB1*LB2.EQ.5.OR.((LB1*LB2.EQ.6).AND.
2790 & (LB1.EQ.3.OR.LB2.EQ.3)))
2791 & .OR. (LB1*LB2.EQ.-3.OR.((LB1*LB2.EQ.-10).AND.
2792 & (LB1.EQ.5.OR.LB2.EQ.5))) )then
2793 XMAX=190.
2794 xmaxn=0
2795 xmaxn1=0
2796 xdirct=dirct1(srt)
2797 go to 678
2798 endif
2799
2800
2801
2802
2803
2804 IF( (LB1*LB2.EQ.3.OR.((LB1*LB2.EQ.10).AND.
2805 & (LB1.EQ.5.OR.LB2.EQ.5)))
2806 & .OR. (LB1*LB2.EQ.-5.OR.((LB1*LB2.EQ.-6).AND.
2807 & (LB1.EQ.3.OR.LB2.EQ.3))) )then
2808 XMAX=27.
2809 xmaxn=2./3.*25.*0.6
2810 xmaxn1=2./3.*40.*0.5
2811 xdirct=dirct2(srt)
2812 go to 678
2813 endif
2814
2815 IF((LB1.EQ.4.OR.LB2.EQ.4).AND.
2816 & (iabs(LB1*LB2).EQ.4.OR.iabs(LB1*LB2).EQ.8))then
2817 XMAX=50.
2818 xmaxn=1./3.*25*0.6
2819 xmaxn1=1/3.*40.*0.5
2820 xdirct=dirct3(srt)
2821 go to 678
2822 endif
2823 678 xnpin1=0
2824 xnpin=0
2825 XNPID=XNPI(I1,I2,1,XMAX)
2826 if(xmaxn1.ne.0)xnpin1=XNPI(i1,i2,2,XMAXN1)
2827 if(xmaxn.ne.0)XNPIN=XNPI(I1,I2,0,XMAXN)
2828
2829 xres=xnpid+xnpin+xnpin1
2830 xnelas=xres+xdirct
2831 icheck=1
2832 go to 34
2833
2834
2835
2836
2837 31 ec=(em1+em2+0.02)**2
2838 xreab=reab(i1,i2,srt,1)
2839
2840
2841 if((iabs(lb1).ge.10.and.iabs(lb1).le.13)
2842 1 .or.(iabs(lb2).ge.10.and.iabs(lb2).le.13)) XREAB = 0.
2843
2844 xkaon=xkaon0+xreab
2845
2846 IF((iabs(LB1).GT.9.AND.iabs(LB1).LE.13) .OR.
2847 & (iabs(LB2).GT.9.AND.iabs(LB2).LE.13))THEN
2848 Xnelas=1.0
2849 ELSE
2850 XNELAS=DPION(EM1,EM2,LB1,LB2,SRT)
2851 ENDIF
2852 icheck=2
2853 34 IF((Xnelas+xkaon+Xphi).LE.0.000001)GO TO 400
2854 DS=SQRT((Xnelas+xkaon+Xphi)/PI)
2855
2856
2857
2858
2859
2860
2861 deltar=ds+0.1
2862 CALL DISTCE(I1,I2,DELTAR,DS,DT,EC,SRT,IC,
2863 1 PCX,PCY,PCZ)
2864 IF(IC.EQ.-1) GO TO 400
2865 ekaon(4,iss)=ekaon(4,iss)+1
2866
2867
2868
2869
2870 if(icheck.eq.2)then
2871
2872 if(xnelas/(xnelas+xkaon+Xphi).ge.RANART(NSEED))then
2873
2874 call Crdir(PX1CM,PY1CM,PZ1CM,SRT,I1,I2,IBLOCK)
2875 go to 440
2876 else
2877
2878
2879 go to 96
2880 endif
2881 endif
2882
2883
2884
2885
2886 IF((XKAON+Xphi)/(XKAON+Xphi+Xnelas).GT.RANART(NSEED))GO TO 95
2887
2888
2889 if(xdirct/xnelas.ge.RANART(NSEED))then
2890
2891 call Crdir(PX1CM,PY1CM,PZ1CM,SRT,I1,I2,IBLOCK)
2892 go to 440
2893 endif
2894
2895 IF( (LB1*LB2.EQ.5.OR.((LB1*LB2.EQ.6).AND.
2896 & (LB1.EQ.3.OR.LB2.EQ.3)))
2897 & .OR. (LB1*LB2.EQ.-3.OR.((LB1*LB2.EQ.-10).AND.
2898 & (LB1.EQ.5.OR.LB2.EQ.5))) )then
2899
2900
2901 GO TO 99
2902 else
2903
2904
2905 XX=(XNPIN+xnpin1)/xres
2906 IF(RANART(NSEED).LT.XX)THEN
2907
2908
2909 xx0=xnpin/(xnpin+xnpin1)
2910 if(RANART(NSEED).lt.xx0)then
2911 RESONA=0.
2912
2913 GO TO 97
2914 else
2915
2916 resona=1.
2917 GO TO 98
2918 endif
2919 ELSE
2920
2921 GO TO 99
2922 ENDIF
2923 ENDIF
2924 97 CONTINUE
2925 IF(RESONA.EQ.0.)THEN
2926
2927 I=I1
2928 IF(EM1.LT.0.6)I=I2
2929
2930 IF( (LB1*LB2.EQ.10.AND.(LB1.EQ.5.OR.LB2.EQ.5))
2931 & .OR.(LB1*LB2.EQ.-6.AND.(LB1.EQ.3.OR.LB2.EQ.3)) )THEN
2932 LB(I)=11
2933 go to 303
2934 ENDIF
2935
2936
2937 IF(iabs(LB(I1)*LB(I2)).EQ.4.AND.
2938 & (LB(I1).EQ.4.OR.LB(I2).EQ.4))THEN
2939 LB(I)=11
2940 go to 303
2941 ENDIF
2942
2943
2944 IF(iabs(LB(I1)*LB(I2)).EQ.8.AND.
2945 & (LB(I1).EQ.4.OR.LB(I2).EQ.4))THEN
2946 LB(I)=10
2947 go to 303
2948 ENDIF
2949
2950
2951 IF( (LB(I1)*LB(I2).EQ.3)
2952 & .OR.(LB(I1)*LB(I2).EQ.-5) )THEN
2953 LB(I)=10
2954 ENDIF
2955 303 CALL DRESON(I1,I2)
2956 if(LB1.lt.0.or.LB2.lt.0) LB(I)=-LB(I)
2957 lres=lres+1
2958 GO TO 101
2959
2960 ENDIF
2961 98 IF(RESONA.EQ.1.)THEN
2962
2963 I=I1
2964 IF(EM1.LT.0.6)I=I2
2965
2966
2967
2968 IF( (LB1*LB2.EQ.10.AND.(LB1.EQ.5.OR.LB2.EQ.5))
2969 & .OR.(LB1*LB2.EQ.-6.AND.(LB1.EQ.3.OR.LB2.EQ.3)) )THEN
2970 LB(I)=13
2971 go to 304
2972 ENDIF
2973
2974
2975 IF(iabs(LB(I1)*LB(I2)).EQ.4.AND.
2976 & (LB(I1).EQ.4.OR.LB(I2).EQ.4))THEN
2977 LB(I)=13
2978 go to 304
2979 ENDIF
2980
2981
2982 IF(iabs(LB(I1)*LB(I2)).EQ.8.AND.
2983 & (LB(I1).EQ.4.OR.LB(I2).EQ.4))THEN
2984 LB(I)=12
2985 go to 304
2986 ENDIF
2987
2988
2989 IF( (LB(I1)*LB(I2).EQ.3)
2990 & .OR.(LB(I1)*LB(I2).EQ.-5) )THEN
2991 LB(I)=12
2992 go to 304
2993 endif
2994
2995 if(lb(i1)*lb(i2).eq.0)then
2996
2997 if(iabs(lb(i1)).eq.1.or.iabs(lb(i2)).eq.1)then
2998 LB(I)=13
2999 go to 304
3000 ELSE
3001 LB(I)=12
3002 ENDIF
3003 endif
3004 304 CALL DRESON(I1,I2)
3005 if(LB1.lt.0.or.LB2.lt.0) LB(I)=-LB(I)
3006 lres=lres+1
3007 GO TO 101
3008
3009 ENDIF
3010
3011
3012 99 LRES=LRES+1
3013 I=I1
3014 IF(EM1.LE.0.6)I=I2
3015
3016
3017 IF( (LB(I1)*LB(I2).EQ.5)
3018 & .OR.(LB(I1)*LB(I2).EQ.-3) )THEN
3019 LB(I)=9
3020 go to 305
3021 ENDIF
3022
3023
3024 IF(iabs(LB(I1)*LB(I2)).EQ.4.AND.(LB(I1).EQ.4.OR.LB(I2).EQ.4))then
3025 LB(I)=8
3026 go to 305
3027 ENDIF
3028
3029
3030 IF( (LB(I1)*LB(I2).EQ.10.AND.(LB(I1).EQ.5.OR.LB(I2).EQ.5))
3031 & .OR.(LB(I1)*LB(I2).EQ.-6.AND.(LB(I1).EQ.3.OR.LB(I2).EQ.3)) )THEN
3032 LB(I)=8
3033 go to 305
3034 ENDIF
3035
3036
3037 IF(iabs(LB(I1)*LB(I2)).EQ.8.AND.(LB(I1).EQ.4.OR.LB(I2).EQ.4))THEN
3038 LB(I)=7
3039 go to 305
3040 ENDIF
3041
3042
3043 IF( (LB(I1)*LB(I2).EQ.3)
3044 & .OR.(LB(I1)*LB(I2).EQ.-5) )THEN
3045 LB(I)=7
3046 go to 305
3047 ENDIF
3048
3049
3050 IF( (LB(I1)*LB(I2).EQ.6.AND.(LB(I1).EQ.3.OR.LB(I2).EQ.3))
3051 & .OR.(LB(I1)*LB(I2).EQ.-10.AND.(LB(I1).EQ.5.OR.LB(I2).EQ.5)) )THEN
3052 LB(I)=6
3053 ENDIF
3054 305 CALL DRESON(I1,I2)
3055 if(LB1.lt.0.or.LB2.lt.0) LB(I)=-LB(I)
3056 GO TO 101
3057
3058
3059
3060
3061
3062
3063
3064 889 CONTINUE
3065 PX1CM=PCX
3066 PY1CM=PCY
3067 PZ1CM=PCZ
3068 EC=(em1+em2+0.02)**2
3069
3070 spika=60./(1.+4.*(srt-0.895)**2/(0.05)**2)
3071
3072
3073
3074 call Crkpla(PX1CM,PY1CM,PZ1CM,EC,SRT,spika,
3075 & emm1,emm2,lbp1,lbp2,I1,I2,icase,srhoks)
3076
3077
3078
3079
3080
3081
3082
3083
3084
3085
3086
3087 if(icase .eq. 0) then
3088 iblock=0
3089 go to 400
3090 endif
3091
3092 if(icase .eq. 1)then
3093 call KSRESO(I1,I2)
3094
3095 iblock = 171
3096
3097
3098
3099
3100
3101
3102
3103 lres=lres+1
3104 go to 101
3105 elseif(icase .eq. 2)then
3106 iblock = 71
3107
3108
3109
3110 elseif(iabs(icase).eq.5)then
3111 iblock = 88
3112
3113 else
3114
3115
3116 iblock = 222
3117 endif
3118 LB(I1) = lbp1
3119 LB(I2) = lbp2
3120 E(I1) = emm1
3121 E(I2) = emm2
3122 em1=e(i1)
3123 em2=e(i2)
3124 ntag = 0
3125 go to 440
3126
3127 33 continue
3128 em1=e(i1)
3129 em2=e(i2)
3130
3131
3132
3133
3134
3135
3136 xelstc=0
3137 if((lb1.ge.25.and.lb1.le.28).and.
3138 & (iabs(lb2).eq.1.or.iabs(lb2).eq.2))
3139 & xelstc=ERHON(EM1,EM2,LB1,LB2,SRT)
3140 if((lb2.ge.25.and.lb2.le.28).and.
3141 & (iabs(lb1).eq.1.or.iabs(lb1).eq.2))
3142 & xelstc=ERHON(EM1,EM2,LB1,LB2,SRT)
3143 ec=(em1+em2+0.02)**2
3144
3145 xkaon0=0
3146 if(srt.ge.1.63.AND.SRT.LE.1.7)xkaon0=pnlka(srt)
3147 IF(SRT.GT.1.7)XKAON0=PNLKA(SRT)+pnska(srt)
3148 if(xkaon0.lt.0)xkaon0=0
3149
3150
3151 XKAON0 = 2.0 * XKAON0
3152
3153
3154
3155 xkaon=xkaon0
3156 ichann=0
3157
3158
3159
3160
3161 Xphi = 0.
3162 if( ( (((lb1.ge.1.and.lb1.le.2).or.
3163 & (lb1.ge.6.and.lb1.le.9))
3164 & .and.(lb2.ge.25.and.lb2.le.27))
3165 & .OR.(((lb2.ge.1.and.lb2.le.2).or.
3166 & (lb2.ge.6.and.lb2.le.9))
3167 & .and.(lb1.ge.25.and.lb1.le.27)) ).AND. srt.gt.1.958)
3168 & call pibphi(srt,lb1,lb2,em1,em2,Xphi,xphin)
3169
3170
3171 if((iabs(lb1).ge.6.and.lb2.ge.25).or.
3172 & (lb1.ge.25.and.iabs(lb2).ge.6))then
3173 ichann=1
3174 ictrl=2
3175 if(lb1.eq.28.or.lb2.eq.28)ictrl=3
3176 xreab=reab(i1,i2,srt,ictrl)
3177
3178
3179 if((iabs(lb1).ge.10.and.iabs(lb1).le.13)
3180 1 .or.(iabs(lb2).ge.10.and.iabs(lb2).le.13)) XREAB = 0.
3181
3182 if(xreab.lt.0)xreab=1.E-06
3183 xkaon=xkaon0+xreab
3184 XELSTC=1.0
3185 endif
3186 DS=SQRT((XKAON+Xphi+xelstc)/PI)
3187
3188
3189
3190
3191
3192
3193
3194 DELTAR=DS+0.1
3195 px1cm=pcx
3196 py1cm=pcy
3197 pz1cm=pcz
3198
3199 CALL DISTCE(I1,I2,DELTAR,DS,DT,EC,SRT,IC,
3200 1 PCX,PCY,PCZ)
3201 IF(IC.EQ.-1) GO TO 400
3202 ekaon(4,iss)=ekaon(4,iss)+1
3203
3204
3205
3206 if(xelstc/(xelstc+xkaon+Xphi).gt.RANART(NSEED))then
3207
3208 call crdir(px1CM,py1CM,pz1CM,srt,I1,i2,IBLOCK)
3209 go to 440
3210 endif
3211
3212 CALL CRRD(PX1CM,PY1CM,PZ1CM,SRT,I1,I2,
3213 1 IBLOCK,xkaon0,xkaon,Xphi,xphin)
3214
3215
3216
3217 IF(IBLOCK.EQ.7) then
3218 LPN=LPN+1
3219 elseIF(IBLOCK.EQ.-7) then
3220 endif
3221
3222
3223 if(iblock.eq.81) lrhor=lrhor+1
3224
3225 if(iblock.eq.82) lomgar=lomgar+1
3226 em1=e(i1)
3227 em2=e(i2)
3228 GO TO 440
3229
3230
3231 95 continue
3232
3233
3234 CALL CRPN(PX1CM,PY1CM,PZ1CM,SRT,I1,I2,
3235 1 IBLOCK,xkaon0,xkaon,Xphi,xphin)
3236
3237
3238
3239 IF(IBLOCK.EQ.7) then
3240 LPN=LPN+1
3241 elseIF(IBLOCK.EQ.-7) then
3242 endif
3243
3244
3245 if(iblock.eq.77) lpd=lpd+1
3246
3247 if(iblock.eq.78) lrho=lrho+1
3248
3249 if(iblock.eq.79) lomega=lomega+1
3250 em1=e(i1)
3251 em2=e(i2)
3252 GO TO 440
3253
3254
3255
3256
3257 96 continue
3258 CALL CRPD(PX1CM,PY1CM,PZ1CM,SRT,I1,I2,
3259 1 IBLOCK,xkaon0,xkaon,Xphi,xphin)
3260
3261
3262
3263 IF(IBLOCK.EQ.7) then
3264 LPN=LPN+1
3265 elseIF(IBLOCK.EQ.-7) then
3266 endif
3267
3268
3269 if(iblock.eq.80) lpdr=lpdr+1
3270 em1=e(i1)
3271 em2=e(i2)
3272 GO TO 440
3273
3274
3275
3276
3277
3278
3279
3280
3281
3282
3283 101 continue
3284 IF(E(I2).EQ.0.)GO TO 600
3285 IF(E(I1).EQ.0.)GO TO 800
3286
3287 44 CONTINUE
3288
3289
3290
3291
3292
3293 cutoff=em1+em2+0.02
3294 IF(SRT.LE.CUTOFF)GO TO 400
3295 IF(SRT.GT.2.245)THEN
3296 SIGNN=PP2(SRT)
3297 ELSE
3298 SIGNN = 35.0 / (1. + (SRT - CUTOFF) * 100.0) + 20.0
3299 ENDIF
3300 call XND(pcx,pcy,pcz,srt,I1,I2,xinel,
3301 & sigk,xsk1,xsk2,xsk3,xsk4,xsk5)
3302 sig=signn+xinel
3303
3304 EC=(EM1+EM2+0.02)**2
3305
3306 PX1CM=PCX
3307 PY1CM=PCY
3308 PZ1CM=PCZ
3309
3310
3311 ianti=0
3312 if(lb(i1).lt.0 .and. lb(i2).lt.0) ianti=1
3313 call sbbdm(srt,sdprod,ianti,lbm,xmm,pfinal)
3314 sig=sig+sdprod
3315
3316 ipdflag=0
3317 if(idpert.eq.1) then
3318 ipert1=1
3319 sigr0=sig
3320 dspert=sqrt(sigr0/pi/10.)
3321 dsrpert=dspert+0.1
3322 CALL DISTCE(I1,I2,dsrpert,dspert,DT,EC,SRT,IC,
3323 1 PX1CM,PY1CM,PZ1CM)
3324 IF(IC.EQ.-1) GO TO 363
3325 signn0=0.
3326 CALL CRND(IRUN,PX1CM,PY1CM,PZ1CM,SRT,I1,I2,
3327 & IBLOCK,SIGNN0,SIGr0,sigk,xsk1,xsk2,xsk3,xsk4,xsk5,NT,ipert1)
3328
3329 ipdflag=1
3330 363 continue
3331 ipert1=0
3332 endif
3333 if(idpert.eq.2) ipert1=1
3334
3335 DS=SQRT(SIG/(10.*PI))
3336 DELTAR=DS+0.1
3337 CALL DISTCE(I1,I2,DELTAR,DS,DT,EC,SRT,IC,
3338 1 PX1CM,PY1CM,PZ1CM)
3339
3340 IF(IC.EQ.-1) then
3341 if(ipdflag.eq.1) iblock=501
3342 GO TO 400
3343 endif
3344
3345 ekaon(3,iss)=ekaon(3,iss)+1
3346
3347
3348 go to 361
3349
3350
3351 361 continue
3352 CALL CRND(IRUN,PX1CM,PY1CM,PZ1CM,SRT,I1,I2,
3353 & IBLOCK,SIGNN,SIG,sigk,xsk1,xsk2,xsk3,xsk4,xsk5,NT,ipert1)
3354
3355 IF(iblock.eq.0.and.ipdflag.eq.1) iblock=501
3356 IF(IBLOCK.EQ.11)THEN
3357 LNDK=LNDK+1
3358 GO TO 400
3359
3360 elseIF(IBLOCK.EQ.-11.or.iblock.eq.501) then
3361 GO TO 400
3362 ENDIF
3363 if(iblock .eq. 222)then
3364
3365 GO TO 400
3366 ENDIF
3367 em1=e(i1)
3368 em2=e(i2)
3369 GO TO 440
3370
3371 4 CONTINUE
3372
3373
3374
3375
3376
3377 CUTOFF=em1+em2+0.14
3378
3379
3380
3381 IF(SRT.GT.2.245)THEN
3382 SIG=ppt(srt)
3383 SIGNN=SIG-PP1(SRT)
3384 ELSE
3385
3386 SIG=XPP(SRT)
3387 IF(ZET(LB(I1))*ZET(LB(I2)).LE.0)SIG=XNP(SRT)
3388 IF(ZET(LB(I1))*ZET(LB(I2)).GT.0)SIG=XPP(SRT)
3389 IF(ZET(LB(I1)).EQ.0.
3390 & AND.ZET(LB(I2)).EQ.0)SIG=XPP(SRT)
3391 if((lb(i1).eq.-1.and.lb(i2).eq.-2) .or.
3392 & (lb(i2).eq.-1.and.lb(i1).eq.-2))sig=xnp(srt)
3393
3394 IF (SRT .LT. 1.897) THEN
3395 SIGNN = SIG
3396 ELSE
3397 SIGNN = 35.0 / (1. + (SRT - 1.897) * 100.0) + 20.0
3398 ENDIF
3399 ENDIF
3400 PX1CM=PCX
3401 PY1CM=PCY
3402 PZ1CM=PCZ
3403
3404
3405
3406 ianti=0
3407 if(lb(i1).lt.0 .and. lb(i2).lt.0) ianti=1
3408 call sbbdm(srt,sdprod,ianti,lbm,xmm,pfinal)
3409 sig=sig+sdprod
3410
3411
3412 ipdflag=0
3413 if(idpert.eq.1) then
3414
3415
3416
3417
3418 ipert1=1
3419 EC=2.012**2
3420
3421
3422 sigr0=sig
3423
3424
3425
3426
3427 dspert=sqrt(sigr0/pi/10.)
3428 dsrpert=dspert+0.1
3429 CALL DISTCE(I1,I2,dsrpert,dspert,DT,EC,SRT,IC,
3430 1 PX1CM,PY1CM,PZ1CM)
3431 IF(IC.EQ.-1) GO TO 365
3432 signn0=0.
3433 CALL CRNN(IRUN,PX1CM,PY1CM,PZ1CM,SRT,I1,I2,IBLOCK,
3434 1 NTAG,signn0,sigr0,NT,ipert1)
3435 ipdflag=1
3436 365 continue
3437 ipert1=0
3438 endif
3439 if(idpert.eq.2) ipert1=1
3440
3441
3442
3443 IF(SIGNN.LE.0) then
3444 if(ipdflag.eq.1) iblock=501
3445 GO TO 400
3446 endif
3447
3448 EC=3.59709
3449 ds=sqrt(sig/pi/10.)
3450 dsr=ds+0.1
3451 IF((E(I1).GE.1.).AND.(e(I2).GE.1.))EC=4.75
3452 CALL DISTCE(I1,I2,dsr,ds,DT,EC,SRT,IC,
3453 1 PX1CM,PY1CM,PZ1CM)
3454
3455
3456 IF(IC.EQ.-1) then
3457 if(ipdflag.eq.1) iblock=501
3458 GO TO 400
3459 endif
3460
3461
3462
3463 go to 362
3464
3465
3466 362 ekaon(1,iss)=ekaon(1,iss)+1
3467 CALL CRNN(IRUN,PX1CM,PY1CM,PZ1CM,SRT,I1,I2,IBLOCK,
3468 1 NTAG,SIGNN,SIG,NT,ipert1)
3469
3470 IF(iblock.eq.0.and.ipdflag.eq.1) iblock=501
3471
3472
3473
3474 IF(IBLOCK.EQ.4.OR.IBLOCK.Eq.9.or.iblock.ge.44.OR.IBLOCK.EQ.-9
3475 & .or.iblock.eq.222.or.iblock.eq.501)THEN
3476
3477
3478
3479
3480 LCOLL=LCOLL+1
3481 if(iblock.eq.4)then
3482 LDIRT=LDIRT+1
3483 elseif(iblock.eq.44)then
3484 LDdrho=LDdrho+1
3485 elseif(iblock.eq.45)then
3486 Lnnrho=Lnnrho+1
3487 elseif(iblock.eq.46)then
3488 Lnnom=Lnnom+1
3489 elseif(iblock .eq. 222)then
3490 elseIF(IBLOCK.EQ.9) then
3491 LNNK=LNNK+1
3492 elseIF(IBLOCK.EQ.-9) then
3493 endif
3494 GO TO 400
3495 ENDIF
3496
3497 em1=e(i1)
3498 em2=e(i2)
3499 GO TO 440
3500
3501
3502
3503 505 continue
3504 ianti=0
3505 if(lb(i1).lt.0 .or. lb(i2).lt.0) ianti=1
3506 call sdmbb(SRT,sdm,ianti)
3507 PX1CM=PCX
3508 PY1CM=PCY
3509 PZ1CM=PCZ
3510
3511 EC=2.012**2
3512 ds=sqrt(sdm/31.4)
3513 dsr=ds+0.1
3514 CALL DISTCE(I1,I2,dsr,ds,DT,EC,SRT,IC,PX1CM,PY1CM,PZ1CM)
3515 IF(IC.EQ.-1) GO TO 400
3516 CALL crdmbb(PX1CM,PY1CM,PZ1CM,SRT,I1,I2,IBLOCK,
3517 1 NTAG,sdm,NT,ianti)
3518 LCOLL=LCOLL+1
3519 GO TO 400
3520
3521
3522
3523 506 continue
3524 ianti=0
3525 if(lb(i1).lt.0 .or. lb(i2).lt.0) ianti=1
3526 call sdbelastic(SRT,sdb)
3527 PX1CM=PCX
3528 PY1CM=PCY
3529 PZ1CM=PCZ
3530
3531 EC=2.012**2
3532 ds=sqrt(sdb/31.4)
3533 dsr=ds+0.1
3534 CALL DISTCE(I1,I2,dsr,ds,DT,EC,SRT,IC,PX1CM,PY1CM,PZ1CM)
3535 IF(IC.EQ.-1) GO TO 400
3536 CALL crdbel(PX1CM,PY1CM,PZ1CM,SRT,I1,I2,IBLOCK,
3537 1 NTAG,sdb,NT,ianti)
3538 LCOLL=LCOLL+1
3539 GO TO 400
3540
3541
3542
3543 444 CONTINUE
3544
3545 CUTOFF=em1+em2+0.02
3546
3547
3548 IF(SRT.LE.CUTOFF)GO TO 400
3549 IF(SRT.GT.2.245)THEN
3550 SIGNN=PP2(SRT)
3551 ELSE
3552 SIGNN = 35.0 / (1. + (SRT - CUTOFF) * 100.0) + 20.0
3553 ENDIF
3554 IF(SIGNN.LE.0)GO TO 400
3555 CALL XDDIN(PCX,PCY,PCZ,SRT,I1,I2,
3556 &XINEL,SIGK,XSK1,XSK2,XSK3,XSK4,XSK5)
3557 SIG=SIGNN+XINEL
3558 EC=(EM1+EM2+0.02)**2
3559 PX1CM=PCX
3560 PY1CM=PCY
3561 PZ1CM=PCZ
3562
3563
3564 ianti=0
3565 if(lb(i1).lt.0 .and. lb(i2).lt.0) ianti=1
3566 call sbbdm(srt,sdprod,ianti,lbm,xmm,pfinal)
3567 sig=sig+sdprod
3568
3569 ipdflag=0
3570 if(idpert.eq.1) then
3571 ipert1=1
3572 sigr0=sig
3573 dspert=sqrt(sigr0/pi/10.)
3574 dsrpert=dspert+0.1
3575 CALL DISTCE(I1,I2,dsrpert,dspert,DT,EC,SRT,IC,
3576 1 PX1CM,PY1CM,PZ1CM)
3577 IF(IC.EQ.-1) GO TO 367
3578 signn0=0.
3579 CALL CRDD(IRUN,PX1CM,PY1CM,PZ1CM,SRT,I1,I2,
3580 1 IBLOCK,NTAG,SIGNN0,SIGr0,NT,ipert1)
3581
3582 ipdflag=1
3583 367 continue
3584 ipert1=0
3585 endif
3586 if(idpert.eq.2) ipert1=1
3587
3588 ds=sqrt(sig/31.4)
3589 dsr=ds+0.1
3590 CALL DISTCE(I1,I2,dsr,ds,DT,EC,SRT,IC,
3591 1 PX1CM,PY1CM,PZ1CM)
3592
3593 IF(IC.EQ.-1) then
3594 if(ipdflag.eq.1) iblock=501
3595 GO TO 400
3596 endif
3597
3598
3599
3600 go to 364
3601
3602
3603 364 ekaon(2,iss)=ekaon(2,iss)+1
3604
3605
3606 CALL CRDD(IRUN,PX1CM,PY1CM,PZ1CM,SRT,I1,I2,
3607 1 IBLOCK,NTAG,SIGNN,SIG,NT,ipert1)
3608
3609 IF(iblock.eq.0.and.ipdflag.eq.1) iblock=501
3610
3611 IF(iabs(IBLOCK).EQ.10)THEN
3612
3613
3614 LCOLL=LCOLL+1
3615 IF(IBLOCK.EQ.10)THEN
3616 LDDK=LDDK+1
3617 elseIF(IBLOCK.EQ.-10) then
3618 endif
3619 GO TO 400
3620 ENDIF
3621
3622
3623 if(iblock .eq. 222.or.iblock.eq.501)then
3624
3625 GO TO 400
3626 ENDIF
3627 em1=e(i1)
3628 em2=e(i2)
3629 GO TO 440
3630
3631 777 CONTINUE
3632 PX1CM=PCX
3633 PY1CM=PCY
3634 PZ1CM=PCZ
3635
3636 ec0=em1+em2+0.02
3637 IF(SRT.LE.ec0)GO TO 400
3638 ec=(em1+em2+0.02)**2
3639
3640
3641
3642
3643 ppel=20.
3644 ipp=1
3645 if(lb1.lt.3.or.lb1.gt.5.or.lb2.lt.3.or.lb2.gt.5)go to 778
3646 CALL PPXS(LB1,LB2,SRT,PPSIG,spprho,IPP)
3647 ppel=ppsig
3648 778 ppink=pipik(srt)
3649
3650
3651
3652 ppink = 2.0 * ppink
3653 if(lb1.ge.25.and.lb2.ge.25) ppink=rrkk
3654
3655
3656
3657
3658 if( ( (lb1.eq.0.or.(lb1.ge.3.and.lb1.le.5))
3659 1 .and.(lb2.ge.25.and.lb2.le.28))
3660 2 .or. ( (lb2.eq.0.or.(lb2.ge.3.and.lb2.le.5))
3661 3 .and.(lb1.ge.25.and.lb1.le.28))) then
3662 ppink=0.
3663 if(srt.ge.(aka+aks)) ppink = prkk
3664 endif
3665
3666
3667 call spprr(lb1,lb2,srt)
3668
3669 call sppee(lb1,lb2,srt)
3670
3671 call spppe(lb1,lb2,srt)
3672
3673 call srpre(lb1,lb2,srt)
3674
3675 call sopoe(lb1,lb2,srt)
3676
3677 call srree(lb1,lb2,srt)
3678
3679 ppinnb=0.
3680 if(srt.gt.thresh(1)) then
3681 call getnst(srt)
3682 if(lb1.ge.3.and.lb1.le.5.and.lb2.ge.3.and.lb2.le.5) then
3683 ppinnb=ppbbar(srt)
3684 elseif((lb1.ge.3.and.lb1.le.5.and.lb2.ge.25.and.lb2.le.27)
3685 1 .or.(lb2.ge.3.and.lb2.le.5.and.lb1.ge.25.and.lb1.le.27)) then
3686 ppinnb=prbbar(srt)
3687 elseif(lb1.ge.25.and.lb1.le.27
3688 1 .and.lb2.ge.25.and.lb2.le.27) then
3689 ppinnb=rrbbar(srt)
3690 elseif((lb1.ge.3.and.lb1.le.5.and.lb2.eq.28)
3691 1 .or.(lb2.ge.3.and.lb2.le.5.and.lb1.eq.28)) then
3692 ppinnb=pobbar(srt)
3693 elseif((lb1.ge.25.and.lb1.le.27.and.lb2.eq.28)
3694 1 .or.(lb2.ge.25.and.lb2.le.27.and.lb1.eq.28)) then
3695 ppinnb=robbar(srt)
3696 elseif(lb1.eq.28.and.lb2.eq.28) then
3697 ppinnb=oobbar(srt)
3698 else
3699 if(lb1.ne.0.and.lb2.ne.0)
3700 1 write(6,*) 'missed MM lb1,lb2=',lb1,lb2
3701 endif
3702 endif
3703 ppin=ppink+ppinnb+pprr+ppee+pppe+rpre+xopoe+rree
3704
3705
3706 if((ppel+ppin).le.0.01)go to 400
3707 DSPP=SQRT((ppel+ppin)/31.4)
3708 dsppr=dspp+0.1
3709 CALL DISTCE(I1,I2,dsppr,DSPP,DT,EC,SRT,IC,
3710 1 PX1CM,PY1CM,PZ1CM)
3711 IF(IC.EQ.-1) GO TO 400
3712 if(ppel.eq.0)go to 400
3713
3714
3715 ekaon(5,iss)=ekaon(5,iss)+1
3716 CALL CRPP(PX1CM,PY1CM,PZ1CM,SRT,I1,I2,
3717 1 IBLOCK,ppel,ppin,spprho,ipp)
3718
3719
3720
3721 if(iblock.eq.666)go to 555
3722 if(iblock.eq.6)LPP=LPP+1
3723 if(iblock.eq.66)then
3724 LPPk=LPPk+1
3725 elseif(iblock.eq.366)then
3726 LPPk=LPPk+1
3727 elseif(iblock.eq.367)then
3728 LPPk=LPPk+1
3729 endif
3730 em1=e(i1)
3731 em2=e(i2)
3732 go to 440
3733
3734
3735
3736
3737 2799 CONTINUE
3738 PX1CM=PCX
3739 PY1CM=PCY
3740 PZ1CM=PCZ
3741 EC=(em1+em2+0.02)**2
3742
3743
3744
3745
3746 DSppb=SQRT(xppbar(srt)/PI/10.)
3747 dsppbr=dsppb+0.1
3748 CALL DISTCE(I1,I2,dsppbr,DSppb,DT,EC,SRT,IC,
3749 1 PX1CM,PY1CM,PZ1CM)
3750 IF(IC.EQ.-1) GO TO 400
3751 CALL Crppba(PX1CM,PY1CM,PZ1CM,SRT,I1,I2,
3752 1 IBLOCK)
3753 em1=e(i1)
3754 em2=e(i2)
3755 go to 440
3756
3757 3555 PX1CM=PCX
3758 PY1CM=PCY
3759 PZ1CM=PCZ
3760 EC=(em1+em2+0.02)**2
3761 DSkk=SQRT(SIG/PI/10.)
3762 dskk0=dskk+0.1
3763 CALL DISTCE(I1,I2,dskk0,DSkk,DT,EC,SRT,IC,
3764 1 PX1CM,PY1CM,PZ1CM)
3765 IF(IC.EQ.-1) GO TO 400
3766 CALL Crlaba(PX1CM,PY1CM,PZ1CM,SRT,brel,brsgm,
3767 & I1,I2,nt,IBLOCK,nchrg,icase)
3768 em1=e(i1)
3769 em2=e(i2)
3770 go to 440
3771
3772
3773 3455 PX1CM=PCX
3774 PY1CM=PCY
3775 PZ1CM=PCZ
3776 call pertur(PX1CM,PY1CM,PZ1CM,SRT,IRUN,I1,I2,nt,kp,icontp)
3777 if(icontp .eq. 0)then
3778
3779 em1 = e(i1)
3780 em2 = e(i2)
3781 iblock = 727
3782 go to 440
3783 endif
3784
3785 if (e(i1) .eq. 0.) go to 800
3786 if (e(i2) .eq. 0.) go to 600
3787 go to 400
3788
3789
3790
3791 7222 CONTINUE
3792 PX1CM=PCX
3793 PY1CM=PCY
3794 PZ1CM=PCZ
3795 EC=(em1+em2+0.02)**2
3796 CALL XphiB(LB1, LB2, EM1, EM2, SRT,
3797 & XSK1, XSK2, XSK3, XSK4, XSK5, SIGP)
3798 DSkk=SQRT(SIGP/PI/10.)
3799 dskk0=dskk+0.1
3800 CALL DISTCE(I1,I2,dskk0,DSkk,DT,EC,SRT,IC,
3801 1 PX1CM,PY1CM,PZ1CM)
3802 IF(IC.EQ.-1) GO TO 400
3803 CALL CRPHIB(PX1CM,PY1CM,PZ1CM,SRT,I1,I2,
3804 & XSK1, XSK2, XSK3, XSK4, XSK5, SIGP, IBLOCK)
3805 em1=e(i1)
3806 em2=e(i2)
3807 go to 440
3808
3809
3810 7444 CONTINUE
3811 PX1CM=PCX
3812 PY1CM=PCY
3813 PZ1CM=PCZ
3814 EC=(em1+em2+0.02)**2
3815 CALL PHIMES(I1, I2, SRT, XSK1, XSK2, XSK3, XSK4, XSK5,
3816 1 XSK6, XSK7, SIGPHI)
3817 DSkk=SQRT(SIGPHI/PI/10.)
3818 dskk0=dskk+0.1
3819 CALL DISTCE(I1,I2,dskk0,DSkk,DT,EC,SRT,IC,
3820 1 PX1CM,PY1CM,PZ1CM)
3821 IF(IC.EQ.-1) GO TO 400
3822
3823 PZRT = p(3,i1)+p(3,i2)
3824 ER1 = sqrt( p(1,i1)**2+p(2,i1)**2+p(3,i1)**2+E(i1)**2 )
3825 ER2 = sqrt( p(1,i2)**2+p(2,i2)**2+p(3,i2)**2+E(i2)**2 )
3826 ERT = ER1+ER2
3827 yy = 0.5*log( (ERT+PZRT)/(ERT-PZRT) )
3828
3829 CALL CRPHIM(PX1CM,PY1CM,PZ1CM,SRT,I1,I2,
3830 & XSK1, XSK2, XSK3, XSK4, XSK5, XSK6, SIGPHI, IKKG, IKKL, IBLOCK)
3831 em1=e(i1)
3832 em2=e(i2)
3833 go to 440
3834
3835
3836 7799 CONTINUE
3837 PX1CM=PCX
3838 PY1CM=PCY
3839 PZ1CM=PCZ
3840 EC=(em1+em2+0.02)**2
3841 call lambar(i1,i2,srt,siglab)
3842 DShn=SQRT(siglab/PI/10.)
3843 dshnr=dshn+0.1
3844 CALL DISTCE(I1,I2,dshnr,DShn,DT,EC,SRT,IC,
3845 1 PX1CM,PY1CM,PZ1CM)
3846 IF(IC.EQ.-1) GO TO 400
3847 CALL Crhb(PX1CM,PY1CM,PZ1CM,SRT,I1,I2,IBLOCK)
3848 em1=e(i1)
3849 em2=e(i2)
3850 go to 440
3851
3852
3853
3854 5699 CONTINUE
3855 PX1CM=PCX
3856 PY1CM=PCY
3857 PZ1CM=PCZ
3858 EC=(em1+em2+0.02)**2
3859 CALL XKHYPE(I1, I2, SRT, XKY1, XKY2, XKY3, XKY4, XKY5,
3860 & XKY6, XKY7, XKY8, XKY9, XKY10, XKY11, XKY12, XKY13,
3861 & XKY14, XKY15, XKY16, XKY17, SIGK)
3862 DSkk=SQRT(sigk/PI)
3863 dskk0=dskk+0.1
3864 CALL DISTCE(I1,I2,dskk0,DSkk,DT,EC,SRT,IC,
3865 1 PX1CM,PY1CM,PZ1CM)
3866 IF(IC.EQ.-1) GO TO 400
3867
3868 if(lb(i1).eq.23 .or. lb(i2).eq.23)then
3869 IKMP = 1
3870 else
3871 IKMP = -1
3872 endif
3873 CALL Crkhyp(PX1CM,PY1CM,PZ1CM,SRT,I1,I2,
3874 & XKY1, XKY2, XKY3, XKY4, XKY5,
3875 & XKY6, XKY7, XKY8, XKY9, XKY10, XKY11, XKY12, XKY13,
3876 & XKY14, XKY15, XKY16, XKY17, SIGK, IKMP,
3877 1 IBLOCK)
3878 em1=e(i1)
3879 em2=e(i2)
3880 go to 440
3881
3882
3883
3884
3885 5999 CONTINUE
3886 PX1CM=PCX
3887 PY1CM=PCY
3888 PZ1CM=PCZ
3889 EC=(em1+em2+0.02)**2
3890 sigkp = 15.
3891
3892
3893 DSkk=SQRT(SIGKP/PI/10.)
3894 dskk0=dskk+0.1
3895 CALL DISTCE(I1,I2,dskk0,DSkk,DT,EC,SRT,IC,
3896 1 PX1CM,PY1CM,PZ1CM)
3897 IF(IC.EQ.-1) GO TO 400
3898
3899 CALL CRLAN(PX1CM,PY1CM,PZ1CM,SRT,I1,I2,IBLOCK)
3900 em1=e(i1)
3901 em2=e(i2)
3902 go to 440
3903
3904
3905
3906 8699 CONTINUE
3907 PX1CM=PCX
3908 PY1CM=PCY
3909 PZ1CM=PCZ
3910 EC=(em1+em2+0.02)**2
3911
3912
3913 CALL Crkphi(PX1CM,PY1CM,PZ1CM,EC,SRT,IBLOCK,
3914 & emm1,emm2,lbp1,lbp2,I1,I2,ikk,icase,rrkk,prkk)
3915 if(icase .eq. 0) then
3916 iblock=0
3917 go to 400
3918 endif
3919
3920
3921 if(lbp1.eq.29.or.lbp2.eq.29) then
3922 PZRT = p(3,i1)+p(3,i2)
3923 ER1 = sqrt( p(1,i1)**2+p(2,i1)**2+p(3,i1)**2+E(i1)**2 )
3924 ER2 = sqrt( p(1,i2)**2+p(2,i2)**2+p(3,i2)**2+E(i2)**2 )
3925 ERT = ER1+ER2
3926 yy = 0.5*log( (ERT+PZRT)/(ERT-PZRT) )
3927
3928 iblock = 222
3929 ntag = 0
3930 endif
3931
3932 LB(I1) = lbp1
3933 LB(I2) = lbp2
3934 E(I1) = emm1
3935 E(I2) = emm2
3936 em1=e(i1)
3937 em2=e(i2)
3938 go to 440
3939
3940
3941 8799 CONTINUE
3942 PX1CM=PCX
3943 PY1CM=PCY
3944 PZ1CM=PCZ
3945 EC=(em1+em2+0.02)**2
3946
3947 CALL Crksph(PX1CM,PY1CM,PZ1CM,EC,SRT,
3948 & emm1,emm2,lbp1,lbp2,I1,I2,ikkg,ikkl,iblock,icase,srhoks)
3949 if(icase .eq. 0) then
3950 iblock=0
3951 go to 400
3952 endif
3953
3954 if(lbp1.eq.29.or.lbp2.eq.20) then
3955
3956 PZRT = p(3,i1)+p(3,i2)
3957 ER1 = sqrt( p(1,i1)**2+p(2,i1)**2+p(3,i1)**2+E(i1)**2 )
3958 ER2 = sqrt( p(1,i2)**2+p(2,i2)**2+p(3,i2)**2+E(i2)**2 )
3959 ERT = ER1+ER2
3960 yy = 0.5*log( (ERT+PZRT)/(ERT-PZRT) )
3961 endif
3962
3963 LB(I1) = lbp1
3964 LB(I2) = lbp2
3965 E(I1) = emm1
3966 E(I2) = emm2
3967 em1=e(i1)
3968 em2=e(i2)
3969 go to 440
3970
3971
3972 888 CONTINUE
3973 PX1CM=PCX
3974 PY1CM=PCY
3975 PZ1CM=PCZ
3976 EC=(em1+em2+0.02)**2
3977 sig = 10.
3978 if(iabs(lb1).eq.14.or.iabs(lb2).eq.14 .or.
3979 & iabs(lb1).eq.30.or.iabs(lb2).eq.30)sig=20.
3980 if(lb1.eq.29.or.lb2.eq.29)sig=5.0
3981
3982 DSkn=SQRT(sig/PI/10.)
3983 dsknr=dskn+0.1
3984 CALL DISTCE(I1,I2,dsknr,DSkn,DT,EC,SRT,IC,
3985 1 PX1CM,PY1CM,PZ1CM)
3986 IF(IC.EQ.-1) GO TO 400
3987 CALL Crkn(PX1CM,PY1CM,PZ1CM,SRT,I1,I2,
3988 1 IBLOCK)
3989 em1=e(i1)
3990 em2=e(i2)
3991 go to 440
3992
3993
3994 440 CONTINUE
3995
3996
3997
3998
3999
4000
4001
4002
4003
4004
4005
4006
4007
4008
4009
4010
4011
4012
4013
4014
4015
4016
4017
4018
4019
4020
4021
4022
4023
4024
4025
4026
4027
4028
4029
4030
4031
4032
4033
4034
4035
4036
4037
4038
4039
4040
4041
4042
4043
4044
4045
4046
4047
4048
4049
4050
4051
4052
4053
4054
4055
4056
4057
4058
4059
4060
4061
4062
4063
4064
4065
4066
4067
4068
4069
4070
4071
4072
4073
4074
4075
4076
4077
4078
4079
4080
4081
4082
4083
4084 IF(IBLOCK.EQ.0) GOTO 400
4085
4086
4087
4088 LCOLL = LCOLL +1
4089
4090 NTAG = 0
4091
4092
4093 E1CM = SQRT (EM1**2 + PX1CM**2 + PY1CM**2 + PZ1CM**2)
4094 P1BETA = PX1CM*BETAX + PY1CM*BETAY + PZ1CM*BETAZ
4095 TRANSF = GAMMA * ( GAMMA * P1BETA / (GAMMA + 1) + E1CM )
4096 Pt1I1 = BETAX * TRANSF + PX1CM
4097 Pt2I1 = BETAY * TRANSF + PY1CM
4098 Pt3I1 = BETAZ * TRANSF + PZ1CM
4099
4100 go to 90002
4101
4102
4103
4104
4105
4106
4107
4108
4109
4110
4111
4112 90002 continue
4113
4114
4115 E2CM = SQRT (EM2**2 + PX1CM**2 + PY1CM**2 + PZ1CM**2)
4116 TRANSF = GAMMA * (-GAMMA*P1BETA / (GAMMA + 1.) + E2CM)
4117 Pt1I2 = BETAX * TRANSF - PX1CM
4118 Pt2I2 = BETAY * TRANSF - PY1CM
4119 Pt3I2 = BETAZ * TRANSF - PZ1CM
4120 go to 90003
4121
4122
4123
4124
4125
4126
4127
4128
4129
4130
4131
4132
4133
4134
4135
4136
4137
4138
4139
4140
4141
4142
4143
4144
4145
4146
4147
4148 90003 IF(IBLOCK.EQ.1) LCNNE=LCNNE+1
4149 IF(IBLOCK.EQ.5) LDD=LDD+1
4150 if(iblock.eq.2) LCNND=LCNND+1
4151 IF(IBLOCK.EQ.8) LKN=LKN+1
4152 if(iblock.eq.43) Ldou=Ldou+1
4153
4154
4155
4156
4157
4158 IF(IBLOCK.EQ.3) LCNDN=LCNDN+1
4159
4160
4161
4162 p(1,i1)=pt1i1
4163 p(2,i1)=pt2i1
4164 p(3,i1)=pt3i1
4165 p(1,i2)=pt1i2
4166 p(2,i2)=pt2i2
4167 p(3,i2)=pt3i2
4168
4169
4170
4171
4172
4173
4174
4175
4176 PX1 = P(1,I1)
4177 PY1 = P(2,I1)
4178 PZ1 = P(3,I1)
4179 EM1 = E(I1)
4180 EM2 = E(I2)
4181 LB1 = LB(I1)
4182 LB2 = LB(I2)
4183 ID(I1) = 2
4184 ID(I2) = 2
4185 E1 = SQRT( EM1**2 + PX1**2 + PY1**2 + PZ1**2 )
4186 ID1 = ID(I1)
4187 go to 90004
4188
4189
4190
4191
4192
4193
4194
4195
4196
4197
4198
4199
4200
4201
4202
4203
4204
4205
4206
4207
4208
4209
4210
4211
4212
4213
4214
4215
4216
4217
4218
4219
4220
4221
4222
4223
4224
4225
4226
4227
4228
4229
4230
4231 90004 continue
4232 AM1=EM1
4233 AM2=EM2
4234
4235
4236
4237 400 CONTINUE
4238
4239
4240
4241
4242
4243
4244
4245
4246
4247
4248
4249
4250
4251
4252
4253
4254
4255
4256
4257
4258
4259
4260
4261
4262
4263
4264
4265
4266
4267
4268
4269
4270
4271
4272
4273
4274
4275
4276
4277
4278
4279
4280
4281
4282
4283
4284
4285
4286
4287
4288
4289
4290
4291
4292
4293
4294
4295
4296
4297
4298
4299
4300
4301
4302
4303
4304
4305
4306
4307
4308
4309
4310
4311
4312
4313
4314
4315
4316
4317
4318
4319
4320
4321
4322
4323
4324
4325
4326
4327
4328
4329
4330
4331
4332
4333
4334
4335
4336
4337
4338
4339
4340
4341
4342
4343
4344
4345
4346
4347
4348
4349
4350
4351
4352
4353
4354
4355
4356
4357
4358
4359
4360
4361
4362
4363
4364
4365
4366
4367
4368
4369
4370 555 continue
4371
4372
4373
4374
4375
4376
4377
4378
4379
4380 600 CONTINUE
4381
4382
4383
4384
4385 798 if(nt.eq.ntmax.and.ipi0dcy.eq.1
4386 1 .and.i1.eq.(MASSR(IRUN)+MSUM)) then
4387 do ipion=1,NNN
4388 if(LPION(ipion,IRUN).eq.4) then
4389 wid=7.85e-9
4390 call resdec(i1,nt,nnn,wid,idecay,ipion)
4391 endif
4392 enddo
4393 endif
4394
4395
4396
4397
4398
4399
4400
4401
4402
4403 800 CONTINUE
4404
4405
4406
4407
4408 N0=MASS+MSUM
4409 DO 1005 N=N0+1,MASSR(IRUN)+MSUM
4410
4411
4412
4413 IF(E(N) .GT. 0. .OR. LB(N) .GT. 5000)THEN
4414
4415 NNN=NNN+1
4416 RPION(1,NNN,IRUN)=R(1,N)
4417 RPION(2,NNN,IRUN)=R(2,N)
4418 RPION(3,NNN,IRUN)=R(3,N)
4419
4420 if(nt.eq.ntmax) then
4421 ftpisv(NNN,IRUN)=ftsv(N)
4422 tfdpi(NNN,IRUN)=tfdcy(N)
4423 endif
4424
4425 PPION(1,NNN,IRUN)=P(1,N)
4426 PPION(2,NNN,IRUN)=P(2,N)
4427 PPION(3,NNN,IRUN)=P(3,N)
4428 EPION(NNN,IRUN)=E(N)
4429 LPION(NNN,IRUN)=LB(N)
4430
4431 PROPI(NNN,IRUN)=PROPER(N)
4432
4433 dppion(NNN,IRUN)=dpertp(N)
4434
4435
4436 ENDIF
4437 1005 CONTINUE
4438 MASSRN(IRUN)=NNN+MASS
4439
4440 1000 CONTINUE
4441
4442
4443
4444
4445
4446
4447
4448
4449 IA=0
4450 IB=0
4451 DO 10001 IRUN=1,NUM
4452 IA=IA+MASSR(IRUN-1)
4453 IB=IB+MASSRN(IRUN-1)
4454 DO 10001 IC=1,MASSRN(IRUN)
4455 IE=IA+IC
4456 IG=IB+IC
4457 IF(IC.LE.MASS)THEN
4458 RT(1,IG)=R(1,IE)
4459 RT(2,IG)=R(2,IE)
4460 RT(3,IG)=R(3,IE)
4461
4462 if(nt.eq.ntmax) then
4463 fttemp(IG)=ftsv(IE)
4464 tft(IG)=tfdcy(IE)
4465 endif
4466
4467 PT(1,IG)=P(1,IE)
4468 PT(2,IG)=P(2,IE)
4469 PT(3,IG)=P(3,IE)
4470 ET(IG)=E(IE)
4471 LT(IG)=LB(IE)
4472 PROT(IG)=PROPER(IE)
4473
4474 dptemp(IG)=dpertp(IE)
4475 ELSE
4476 I0=IC-MASS
4477 RT(1,IG)=RPION(1,I0,IRUN)
4478 RT(2,IG)=RPION(2,I0,IRUN)
4479 RT(3,IG)=RPION(3,I0,IRUN)
4480
4481 if(nt.eq.ntmax) then
4482 fttemp(IG)=ftpisv(I0,IRUN)
4483 tft(IG)=tfdpi(I0,IRUN)
4484 endif
4485
4486 PT(1,IG)=PPION(1,I0,IRUN)
4487 PT(2,IG)=PPION(2,I0,IRUN)
4488 PT(3,IG)=PPION(3,I0,IRUN)
4489 ET(IG)=EPION(I0,IRUN)
4490 LT(IG)=LPION(I0,IRUN)
4491 PROT(IG)=PROPI(I0,IRUN)
4492
4493 dptemp(IG)=dppion(I0,IRUN)
4494 ENDIF
4495 10001 CONTINUE
4496
4497 IL=0
4498
4499
4500 DO 10003 IRUN=1,NUM
4501
4502 MASSR(IRUN)=MASSRN(IRUN)
4503 IL=IL+MASSR(IRUN-1)
4504 DO 10002 IM=1,MASSR(IRUN)
4505 IN=IL+IM
4506 R(1,IN)=RT(1,IN)
4507 R(2,IN)=RT(2,IN)
4508 R(3,IN)=RT(3,IN)
4509
4510 if(nt.eq.ntmax) then
4511 ftsv(IN)=fttemp(IN)
4512 tfdcy(IN)=tft(IN)
4513 endif
4514 P(1,IN)=PT(1,IN)
4515 P(2,IN)=PT(2,IN)
4516 P(3,IN)=PT(3,IN)
4517 E(IN)=ET(IN)
4518 LB(IN)=LT(IN)
4519 PROPER(IN)=PROT(IN)
4520
4521 dpertp(IN)=dptemp(IN)
4522 IF(LB(IN).LT.1.OR.LB(IN).GT.2)ID(IN)=0
4523 10002 CONTINUE
4524
4525
4526
4527
4528
4529
4530
4531
4532
4533 10003 CONTINUE
4534
4535 RETURN
4536 END
4537
4538
4539
4540
4541
4542
4543
4544
4545
4546
4547
4548
4549
4550
4551
4552
4553
4554
4555
4556
4557
4558
4559
4560
4561
4562
4563
4564
4565
4566
4567
4568
4569
4570
4571
4572
4573
4574
4575
4576
4577
4578
4579
4580
4581
4582
4583
4584
4585 SUBROUTINE CMS(I1,I2,PX1CM,PY1CM,PZ1CM,SRT)
4586
4587
4588
4589
4590 PARAMETER (MAXSTR=150001)
4591 double precision px1,py1,pz1,px2,py2,pz2,em1,em2,e1,e2,
4592 1 s,ETOTAL,P1BETA,TRANSF,dBETAX,dBETAY,dBETAZ,dGAMMA,scheck
4593 COMMON /BB/ P(3,MAXSTR)
4594 COMMON /CC/ E(MAXSTR)
4595 COMMON /BG/ BETAX,BETAY,BETAZ,GAMMA
4596 SAVE
4597 PX1=dble(P(1,I1))
4598 PY1=dble(P(2,I1))
4599 PZ1=dble(P(3,I1))
4600 PX2=dble(P(1,I2))
4601 PY2=dble(P(2,I2))
4602 PZ2=dble(P(3,I2))
4603 EM1=dble(E(I1))
4604 EM2=dble(E(I2))
4605 E1=dSQRT(EM1**2+PX1**2+PY1**2+PZ1**2)
4606 E2=dSQRT(EM2**2+PX2**2+PY2**2+PZ2**2)
4607 S=(E1+E2)**2-(PX1+PX2)**2-(PY1+PY2)**2-(PZ1+PZ2)**2
4608 IF(S.LE.0) S=0d0
4609 SRT=sngl(dSQRT(S))
4610
4611 ETOTAL = E1 + E2
4612 dBETAX = (PX1+PX2) / ETOTAL
4613 dBETAY = (PY1+PY2) / ETOTAL
4614 dBETAZ = (PZ1+PZ2) / ETOTAL
4615
4616 scheck=1.d0-dBETAX**2-dBETAY**2-dBETAZ**2
4617 if(scheck.le.0d0) then
4618 write(99,*) 'scheck1: ', scheck
4619 stop
4620 endif
4621 dGAMMA=1.d0/dSQRT(scheck)
4622
4623 P1BETA = PX1*dBETAX + PY1*dBETAY + PZ1 * dBETAZ
4624 TRANSF = dGAMMA * ( dGAMMA * P1BETA / (dGAMMA + 1d0) - E1 )
4625 PX1CM = sngl(dBETAX * TRANSF + PX1)
4626 PY1CM = sngl(dBETAY * TRANSF + PY1)
4627 PZ1CM = sngl(dBETAZ * TRANSF + PZ1)
4628 BETAX = sngl(dBETAX)
4629 BETAY = sngl(dBETAY)
4630 BETAZ = sngl(dBETAZ)
4631 GAMMA = sngl(dGAMMA)
4632 RETURN
4633 END
4634
4635
4636
4637 SUBROUTINE DISTCE(I1,I2,DELTAR,DS,DT,EC,SRT
4638 1 ,IC,PX1CM,PY1CM,PZ1CM)
4639
4640
4641
4642
4643
4644
4645
4646
4647
4648
4649
4650 PARAMETER (MAXSTR=150001)
4651 COMMON /AA/ R(3,MAXSTR)
4652
4653 COMMON /BB/ P(3,MAXSTR)
4654
4655 COMMON /CC/ E(MAXSTR)
4656
4657 COMMON /BG/ BETAX,BETAY,BETAZ,GAMMA
4658 COMMON /EE/ ID(MAXSTR),LB(MAXSTR)
4659
4660 common/leadng/lb1,px1,py1,pz1,em1,e1,xfnl,yfnl,zfnl,tfnl,
4661 1 px1n,py1n,pz1n,dp1n
4662 common /dpi/em2,lb2
4663 SAVE
4664 IC=0
4665 X1=R(1,I1)
4666 Y1=R(2,I1)
4667 Z1=R(3,I1)
4668 PX1=P(1,I1)
4669 PY1=P(2,I1)
4670 PZ1=P(3,I1)
4671 X2=R(1,I2)
4672 Y2=R(2,I2)
4673 Z2=R(3,I2)
4674 PX2=P(1,I2)
4675 PY2=P(2,I2)
4676 PZ2=P(3,I2)
4677 EM1=E(I1)
4678 EM2=E(I2)
4679 E1=SQRT(EM1**2+PX1**2+PY1**2+PZ1**2)
4680
4681
4682
4683 RSQARE = (X1-X2)**2 + (Y1-Y2)**2 + (Z1-Z2)**2
4684 IF (RSQARE .GT. DELTAR**2) GO TO 400
4685
4686 E2 = SQRT ( EM2**2 + PX2**2 + PY2**2 + PZ2**2 )
4687 S = SRT*SRT
4688 IF (S .LT. EC) GO TO 400
4689
4690
4691
4692
4693 P1BETA = PX1*BETAX + PY1*BETAY + PZ1 * BETAZ
4694 TRANSF = GAMMA * ( GAMMA * P1BETA / (GAMMA + 1) - E1 )
4695 PRCM = SQRT (PX1CM**2 + PY1CM**2 + PZ1CM**2)
4696 IF (PRCM .LE. 0.00001) GO TO 400
4697
4698 DRBETA = BETAX*(X1-X2) + BETAY*(Y1-Y2) + BETAZ*(Z1-Z2)
4699 TRANSF = GAMMA * GAMMA * DRBETA / (GAMMA + 1)
4700 DXCM = BETAX * TRANSF + X1 - X2
4701 DYCM = BETAY * TRANSF + Y1 - Y2
4702 DZCM = BETAZ * TRANSF + Z1 - Z2
4703
4704 DRCM = SQRT (DXCM**2 + DYCM**2 + DZCM**2 )
4705 DZZ = (PX1CM*DXCM + PY1CM*DYCM + PZ1CM*DZCM) / PRCM
4706 if ((drcm**2 - dzz**2) .le. 0.) then
4707 BBB = 0.
4708 else
4709 BBB = SQRT (DRCM**2 - DZZ**2)
4710 end if
4711
4712 IF (BBB .GT. DS) GO TO 400
4713 RELVEL = PRCM * (1.0/E1 + 1.0/E2)
4714 DDD = RELVEL * DT * 0.5
4715
4716 IF (ABS(DDD) .LT. ABS(DZZ)) GO TO 400
4717 IC=1
4718 GO TO 500
4719 400 IC=-1
4720 500 CONTINUE
4721 RETURN
4722 END
4723
4724
4725
4726 SUBROUTINE CRNN(IRUN,PX,PY,PZ,SRT,I1,I2,IBLOCK,
4727 1NTAG,SIGNN,SIG,NT,ipert1)
4728
4729
4730
4731
4732
4733
4734
4735
4736
4737
4738
4739
4740
4741
4742
4743
4744
4745
4746
4747
4748
4749
4750
4751
4752
4753
4754
4755
4756
4757
4758
4759
4760
4761
4762
4763
4764
4765
4766
4767
4768
4769
4770
4771
4772
4773
4774
4775
4776
4777
4778
4779
4780
4781
4782
4783
4784
4785
4786
4787
4788
4789
4790
4791
4792
4793
4794
4795
4796
4797
4798
4799
4800
4801 PARAMETER (MAXSTR=150001,MAXR=1,AMN=0.939457,
4802 1 AMP=0.93828,AP1=0.13496,aka=0.498,AP2=0.13957,AM0=1.232,
4803 2 PI=3.1415926,CUTOFF=1.8966,AVMASS=0.9383,APHI=1.020)
4804 parameter (MX=4,MY=4,MZ=8,MPX=4,MPY=4,mpz=10,mpzp=10)
4805 parameter (xmd=1.8756,npdmax=10000)
4806 COMMON /AA/ R(3,MAXSTR)
4807
4808 COMMON /BB/ P(3,MAXSTR)
4809
4810 COMMON /CC/ E(MAXSTR)
4811
4812 COMMON /EE/ ID(MAXSTR),LB(MAXSTR)
4813
4814 common /ff/f(-mx:mx,-my:my,-mz:mz,-mpx:mpx,-mpy:mpy,-mpz:mpzp)
4815
4816 common /gg/ dx,dy,dz,dpx,dpy,dpz
4817
4818 COMMON /INPUT/ NSTAR,NDIRCT,DIR
4819
4820 COMMON /NN/NNN
4821
4822 COMMON /BG/BETAX,BETAY,BETAZ,GAMMA
4823
4824 COMMON /RUN/NUM
4825
4826 COMMON /PA/RPION(3,MAXSTR,MAXR)
4827
4828 COMMON /PB/PPION(3,MAXSTR,MAXR)
4829
4830 COMMON /PC/EPION(MAXSTR,MAXR)
4831
4832 COMMON /PD/LPION(MAXSTR,MAXR)
4833
4834 COMMON/TABLE/ xarray(0:1000),earray(0:1000)
4835
4836 common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
4837
4838 common/leadng/lb1,px1,py1,pz1,em1,e1,xfnl,yfnl,zfnl,tfnl,
4839 1 px1n,py1n,pz1n,dp1n
4840
4841 COMMON/RNDF77/NSEED
4842
4843 common /dpi/em2,lb2
4844 COMMON /dpert/dpertt(MAXSTR,MAXR),dpertp(MAXSTR),dplast(MAXSTR),
4845 1 dpdcy(MAXSTR),dpdpi(MAXSTR,MAXR),dpt(MAXSTR, MAXR),
4846 2 dpp1(MAXSTR,MAXR),dppion(MAXSTR,MAXR)
4847 common /para8/ idpert,npertd,idxsec
4848 dimension ppd(3,npdmax),lbpd(npdmax)
4849 SAVE
4850
4851 n12=0
4852 m12=0
4853 IBLOCK=0
4854 NTAG=0
4855 EM1=E(I1)
4856 EM2=E(I2)
4857 PR=SQRT( PX**2 + PY**2 + PZ**2 )
4858 C2=PZ / PR
4859 X1=RANART(NSEED)
4860 ianti=0
4861 if(lb(i1).lt.0 .and. lb(i2).lt.0) ianti=1
4862 call sbbdm(srt,sdprod,ianti,lbm,xmm,pfinal)
4863
4864 if(idpert.eq.1.and.ipert1.eq.1) then
4865 IF (SRT .LT. 2.012) RETURN
4866 if((iabs(lb(i1)).eq.1.or.iabs(lb(i1)).eq.2)
4867 1 .and.(iabs(lb(i2)).eq.1.or.iabs(lb(i2)).eq.2)) then
4868 goto 108
4869 else
4870 return
4871 endif
4872 endif
4873
4874
4875
4876
4877
4878 IF (X1.LE.(SIGNN/SIG)) THEN
4879
4880 AS = ( 3.65 * (SRT - 1.8766) )**6
4881 A = 6.0 * AS / (1.0 + AS)
4882 TA = -2.0 * PR**2
4883 X = RANART(NSEED)
4884
4885 T1 = sngl(DLOG(dble(1.-X)*DEXP(dble(A)*dble(TA))+dble(X)))/ A
4886 C1 = 1.0 - T1/TA
4887 T1 = 2.0 * PI * RANART(NSEED)
4888 IBLOCK=1
4889 GO TO 107
4890 ELSE
4891
4892
4893
4894
4895 IF (SRT .LT. 2.012) RETURN
4896
4897
4898
4899
4900 call N1535(iabs(lb(i1)),iabs(lb(i2)),srt,x1535)
4901
4902
4903
4904 SIG3=3.*(X3pi(SRT)+x33pi(srt))
4905
4906 SIG4=4.*X2pi(srt)
4907
4908 s4pi=x4pi(srt)
4909
4910 srho=xrho(srt)
4911
4912 somega=omega(srt)
4913
4914
4915 akp=0.498
4916 ak0=0.498
4917 ana=0.94
4918 ada=1.232
4919 al=1.1157
4920 as=1.1197
4921 xsk1=0
4922 xsk2=0
4923 xsk3=0
4924 xsk4=0
4925 xsk5=0
4926 t1nlk=ana+al+akp
4927 if(srt.le.t1nlk)go to 222
4928 XSK1=1.5*PPLPK(SRT)
4929
4930 t1dlk=ada+al+akp
4931 t2dlk=ada+al-akp
4932 if(srt.le.t1dlk)go to 222
4933 es=srt
4934 pmdlk2=(es**2-t1dlk**2)*(es**2-t2dlk**2)/(4.*es**2)
4935 pmdlk=sqrt(pmdlk2)
4936 XSK3=1.5*PPLPK(srt)
4937
4938 t1nsk=ana+as+akp
4939 t2nsk=ana+as-akp
4940 if(srt.le.t1nsk)go to 222
4941 pmnsk2=(es**2-t1nsk**2)*(es**2-t2nsk**2)/(4.*es**2)
4942 pmnsk=sqrt(pmnsk2)
4943 XSK2=1.5*(PPK1(srt)+PPK0(srt))
4944
4945 t1DSk=aDa+aS+akp
4946 t2DSk=aDa+aS-akp
4947 if(srt.le.t1dsk)go to 222
4948 pmDSk2=(es**2-t1DSk**2)*(es**2-t2DSk**2)/(4.*es**2)
4949 pmDSk=sqrt(pmDSk2)
4950 XSK4=1.5*(PPK1(srt)+PPK0(srt))
4951
4952
4953 if(srt.le.(2.*amn+aphi))go to 222
4954
4955 xsk5 = 0.0001
4956
4957
4958
4959 222 SIGK=XSK1+XSK2+XSK3+XSK4
4960
4961
4962 XSK1 = 2.0 * XSK1
4963 XSK2 = 2.0 * XSK2
4964 XSK3 = 2.0 * XSK3
4965 XSK4 = 2.0 * XSK4
4966 SIGK = 2.0 * SIGK + xsk5
4967
4968
4969
4970
4971
4972 lb1=iabs(lb(i1))
4973 lb2=iabs(lb(i2))
4974 IF((LB(I1)*LB(I2).EQ.1).or.
4975 & ((lb1.le.17.and.lb1.ge.14).and.(lb2.le.17.and.lb2.ge.14)).
4976 & or.((lb1.le.2).and.(lb2.le.17.and.lb2.ge.14)).
4977 & or.((lb2.le.2).and.(lb1.le.17.and.lb1.ge.14)))THEN
4978
4979 IF(X1.LE.((SIGNN+sdprod)/SIG)) GO TO 108
4980 SIG1=SIGMA(SRT,1,1,0)+0.5*SIGMA(SRT,1,1,1)
4981 SIG2=1.5*SIGMA(SRT,1,1,1)
4982 SIGND=SIG1+SIG2+SIG3+SIG4+X1535+SIGK+s4pi+srho+somega
4983
4984
4985 IF (X1.GT.(SIGNN+SIGND+sdprod)/SIG)RETURN
4986 DIR=SIG3/SIGND
4987 IF(RANART(NSEED).LE.DIR)GO TO 106
4988 IF(RANART(NSEED).LE.SIGK/(SIGK+X1535+SIG4+SIG2+SIG1
4989 & +s4pi+srho+somega))GO TO 306
4990 if(RANART(NSEED).le.s4pi/(x1535+sig4+sig2+sig1
4991 & +s4pi+srho+somega))go to 307
4992 if(RANART(NSEED).le.srho/(x1535+sig4+sig2+sig1
4993 & +srho+somega))go to 308
4994 if(RANART(NSEED).le.somega/(x1535+sig4+sig2+sig1
4995 & +somega))go to 309
4996 if(RANART(NSEED).le.x1535/(sig1+sig2+sig4+x1535))then
4997
4998 N12=9
4999 ELSE
5000 IF(RANART(NSEED).LE.SIG4/(SIG1+sig2+sig4))THEN
5001
5002 N12=66
5003 GO TO 1012
5004 else
5005
5006 N12=3
5007 IF (RANART(NSEED).GT.SIG1/(SIG1+SIG2))N12=4
5008 ENDIF
5009 endif
5010 GO TO 1011
5011 ENDIF
5012
5013 IF(iabs(LB(I1)).EQ.2.AND.iabs(LB(I2)).EQ.2)THEN
5014
5015 IF(X1.LE.((SIGNN+sdprod)/SIG)) GO TO 108
5016 SIG1=SIGMA(SRT,1,1,0)+0.5*SIGMA(SRT,1,1,1)
5017 SIG2=1.5*SIGMA(SRT,1,1,1)
5018 SIGND=SIG1+SIG2+X1535+SIG3+SIG4+SIGK+s4pi+srho+somega
5019
5020
5021 IF (X1.GT.(SIGNN+SIGND+sdprod)/SIG)RETURN
5022 dir=sig3/signd
5023 IF(RANART(NSEED).LE.DIR)GO TO 106
5024 IF(RANART(NSEED).LE.SIGK/(SIGK+X1535+SIG4+SIG2+SIG1
5025 & +s4pi+srho+somega))GO TO 306
5026 if(RANART(NSEED).le.s4pi/(x1535+sig4+sig2+sig1
5027 & +s4pi+srho+somega))go to 307
5028 if(RANART(NSEED).le.srho/(x1535+sig4+sig2+sig1
5029 & +srho+somega))go to 308
5030 if(RANART(NSEED).le.somega/(x1535+sig4+sig2+sig1
5031 & +somega))go to 309
5032 IF(RANART(NSEED).LE.X1535/(x1535+sig1+sig2+sig4))THEN
5033
5034 N12=10
5035 ELSE
5036 if(RANART(NSEED).le.sig4/(sig1+sig2+sig4))then
5037
5038 N12=67
5039 GO TO 1013
5040 else
5041
5042 N12=6
5043 IF (RANART(NSEED).GT.SIG1/(SIG1+SIG2))N12=5
5044 ENDIF
5045 endif
5046 GO TO 1011
5047 ENDIF
5048
5049 IF(LB(I1)*LB(I2).EQ.2)THEN
5050
5051 IF(X1.LE.((SIGNN+sdprod)/SIG)) GO TO 108
5052 SIG1=0.5*SIGMA(SRT,1,1,1)+0.25*SIGMA(SRT,1,1,0)
5053 IF(NSTAR.EQ.1)THEN
5054 SIG2=(3./4.)*SIGMA(SRT,2,0,1)
5055 ELSE
5056 SIG2=0.
5057 ENDIF
5058 SIGND=2.*(SIG1+SIG2+X1535)+sig3+sig4+SIGK+s4pi+srho+somega
5059
5060
5061 IF (X1.GT.(SIGNN+SIGND+sdprod)/SIG)RETURN
5062 dir=sig3/signd
5063 IF(RANART(NSEED).LE.DIR)GO TO 106
5064 IF(RANART(NSEED).LE.SIGK/(SIGND-SIG3))GO TO 306
5065 if(RANART(NSEED).le.s4pi/(signd-sig3-sigk))go to 307
5066 if(RANART(NSEED).le.srho/(signd-sig3-sigk-s4pi))go to 308
5067 if(RANART(NSEED).le.somega/(signd-sig3-sigk-s4pi-srho))
5068 1 go to 309
5069 IF(RANART(NSEED).LT.X1535/(SIG1+SIG2+X1535+0.5*sig4))THEN
5070
5071 N12=11
5072 IF(RANART(NSEED).LE.0.5)N12=12
5073 ELSE
5074 if(RANART(NSEED).le.sig4/(sig4+2.*(sig1+sig2)))then
5075
5076 N12=68
5077 GO TO 1014
5078 else
5079 IF(RANART(NSEED).LE.SIG1/(SIG1+SIG2))THEN
5080
5081 N12=2
5082 IF(RANART(NSEED).GE.0.5)N12=1
5083 ELSE
5084
5085 N12=8
5086 IF(RANART(NSEED).GE.0.5)N12=7
5087 ENDIF
5088 ENDIF
5089 ENDIF
5090 endif
5091 1011 iblock=2
5092 CONTINUE
5093
5094
5095
5096
5097 DMAX = SRT - AVMASS-0.005
5098 DMAX = SRT - AVMASS-0.005
5099 DMIN = 1.078
5100 IF(N12.LT.7)THEN
5101
5102 IF(DMAX.LT.1.232) THEN
5103 FM=FDE(DMAX,SRT,0.)
5104 ELSE
5105
5106
5107 xdmass=1.232
5108
5109 FM=FDE(xdmass,SRT,1.)
5110
5111
5112 ENDIF
5113 IF(FM.EQ.0.)FM=1.E-09
5114 NTRY1=0
5115 10 DM = RANART(NSEED) * (DMAX-DMIN) + DMIN
5116 NTRY1=NTRY1+1
5117 IF((RANART(NSEED) .GT. FDE(DM,SRT,1.)/FM).AND.
5118 1 (NTRY1.LE.30)) GOTO 10
5119
5120
5121
5122 if(dm.gt.1.47) goto 10
5123
5124 GO TO 13
5125 ENDIF
5126 IF((n12.eq.7).or.(n12.eq.8))THEN
5127
5128 IF(DMAX.LT.1.44) THEN
5129 FM=FNS(DMAX,SRT,0.)
5130 ELSE
5131
5132
5133 xdmass=1.44
5134
5135 FM=FNS(xdmass,SRT,1.)
5136
5137
5138 ENDIF
5139 IF(FM.EQ.0.)FM=1.E-09
5140 NTRY2=0
5141 11 DM=RANART(NSEED)*(DMAX-DMIN)+DMIN
5142 NTRY2=NTRY2+1
5143 IF((RANART(NSEED).GT.FNS(DM,SRT,1.)/FM).AND.
5144 1 (NTRY2.LE.10)) GO TO 11
5145
5146
5147
5148 if(dm.gt.2.14) goto 11
5149
5150 GO TO 13
5151 ENDIF
5152 IF(n12.ge.17)then
5153
5154 IF(DMAX.LT.1.535) THEN
5155 FM=FD5(DMAX,SRT,0.)
5156 ELSE
5157
5158
5159 xdmass=1.535
5160
5161 FM=FD5(xdmass,SRT,1.)
5162
5163
5164 ENDIF
5165 IF(FM.EQ.0.)FM=1.E-09
5166 NTRY1=0
5167 12 DM = RANART(NSEED) * (DMAX-DMIN) + DMIN
5168 NTRY1=NTRY1+1
5169 IF((RANART(NSEED) .GT. FD5(DM,SRT,1.)/FM).AND.
5170 1 (NTRY1.LE.10)) GOTO 12
5171
5172
5173
5174 if(dm.gt.1.84) goto 12
5175
5176 GO TO 13
5177 ENDIF
5178
5179
5180 1012 iblock=43
5181 call Rmasdd(srt,1.232,1.232,1.08,
5182 & 1.08,ISEED,1,dm1,dm2)
5183 call Rmasdd(srt,1.232,1.44,1.08,
5184 & 1.08,ISEED,3,dm1n,dm2n)
5185 IF(N12.EQ.66)THEN
5186
5187
5188 XFINAL=RANART(NSEED)
5189 IF(XFINAL.LE.0.25)THEN
5190
5191 LB(I1)=9
5192 LB(I2)=7
5193 e(i1)=dm1
5194 e(i2)=dm2
5195 GO TO 200
5196
5197 ENDIF
5198 IF((XFINAL.gt.0.25).and.(xfinal.le.0.5))THEN
5199
5200 LB(I1)=8
5201 LB(I2)=8
5202 e(i1)=dm1
5203 e(i2)=dm2
5204 GO TO 200
5205
5206 ENDIF
5207 IF((XFINAL.gt.0.5).and.(xfinal.le.0.75))THEN
5208
5209 LB(I1)=9
5210 LB(I2)=10
5211 e(i1)=dm1n
5212 e(i2)=dm2n
5213 GO TO 200
5214
5215 ENDIF
5216 IF(XFINAL.gt.0.75)then
5217
5218 LB(I1)=8
5219 LB(I2)=11
5220 e(i1)=dm1n
5221 e(i2)=dm2n
5222 GO TO 200
5223
5224 ENDIF
5225 ENDIF
5226 1013 iblock=43
5227 call Rmasdd(srt,1.232,1.232,1.08,
5228 & 1.08,ISEED,1,dm1,dm2)
5229 call Rmasdd(srt,1.232,1.44,1.08,
5230 & 1.08,ISEED,3,dm1n,dm2n)
5231 IF(N12.EQ.67)THEN
5232
5233
5234 XFINAL=RANART(NSEED)
5235 IF(XFINAL.LE.0.25)THEN
5236
5237 LB(I1)=7
5238 LB(I2)=7
5239 e(i1)=dm1
5240 e(i2)=dm2
5241 GO TO 200
5242
5243 ENDIF
5244 IF((XFINAL.gt.0.25).and.(xfinal.le.0.5))THEN
5245
5246 LB(I1)=6
5247 LB(I2)=8
5248 e(i1)=dm1
5249 e(i2)=dm2
5250 GO TO 200
5251
5252 ENDIF
5253 IF((XFINAL.gt.0.5).and.(xfinal.le.0.75))THEN
5254
5255 LB(I1)=7
5256 LB(I2)=10
5257 e(i1)=dm1n
5258 e(i2)=dm2n
5259 GO TO 200
5260
5261 ENDIF
5262 IF(XFINAL.gt.0.75)then
5263
5264 LB(I1)=8
5265 LB(I2)=11
5266 e(i1)=dm1n
5267 e(i2)=dm2n
5268 GO TO 200
5269
5270 ENDIF
5271 ENDIF
5272 1014 iblock=43
5273 call Rmasdd(srt,1.232,1.232,1.08,
5274 & 1.08,ISEED,1,dm1,dm2)
5275 call Rmasdd(srt,1.232,1.44,1.08,
5276 & 1.08,ISEED,3,dm1n,dm2n)
5277 IF(N12.EQ.68)THEN
5278
5279
5280 XFINAL=RANART(NSEED)
5281 IF(XFINAL.LE.0.25)THEN
5282
5283 LB(I1)=7
5284 LB(I2)=8
5285 e(i1)=dm1
5286 e(i2)=dm2
5287 GO TO 200
5288
5289 ENDIF
5290 IF((XFINAL.gt.0.25).and.(xfinal.le.0.5))THEN
5291
5292 LB(I1)=9
5293 LB(I2)=6
5294 e(i1)=dm1
5295 e(i2)=dm2
5296 GO TO 200
5297
5298 ENDIF
5299 IF((XFINAL.gt.0.5).and.(xfinal.le.0.75))THEN
5300
5301 LB(I1)=7
5302 LB(I2)=11
5303 e(i1)=dm1n
5304 e(i2)=dm2n
5305 GO TO 200
5306
5307 ENDIF
5308 IF(XFINAL.gt.0.75)then
5309
5310 LB(I1)=8
5311 LB(I2)=10
5312 e(i1)=dm1n
5313 e(i2)=dm2n
5314 GO TO 200
5315
5316 ENDIF
5317 ENDIF
5318 13 CONTINUE
5319
5320
5321
5322 IF(N12.EQ.1)THEN
5323 IF(iabs(LB(I1)).EQ.1)THEN
5324 LB(I2)=2
5325 LB(I1)=8
5326 E(I1)=DM
5327 ELSE
5328 LB(I1)=2
5329 LB(I2)=8
5330 E(I2)=DM
5331 ENDIF
5332 GO TO 200
5333 ENDIF
5334
5335 IF(N12.EQ.2)THEN
5336 IF(iabs(LB(I1)).EQ.2)THEN
5337 LB(I2)=1
5338 LB(I1)=7
5339 E(I1)=DM
5340 ELSE
5341 LB(I1)=1
5342 LB(I2)=7
5343 E(I2)=DM
5344 ENDIF
5345 GO TO 200
5346 ENDIF
5347
5348 IF(N12.EQ.3)THEN
5349 LB(I1)=9
5350 E(I1)=DM
5351 LB(I2)=2
5352 E(I2)=AMN
5353 GO TO 200
5354 ENDIF
5355
5356 IF(N12.EQ.4)THEN
5357 LB(I2)=1
5358 LB(I1)=8
5359 E(I1)=DM
5360 GO TO 200
5361 ENDIF
5362
5363 IF(N12.EQ.5)THEN
5364 LB(I2)=2
5365 LB(I1)=7
5366 E(I1)=DM
5367 GO TO 200
5368 ENDIF
5369
5370 IF(N12.EQ.6)THEN
5371 LB(I1)=6
5372 E(I1)=DM
5373 LB(I2)=1
5374 E(I2)=AMP
5375 GO TO 200
5376 ENDIF
5377
5378 IF(N12.EQ.7)THEN
5379 IF(iabs(LB(I1)).EQ.1)THEN
5380 LB(I1)=1
5381 LB(I2)=10
5382 E(I2)=DM
5383 ELSE
5384 LB(I2)=1
5385 LB(I1)=10
5386 E(I1)=DM
5387 ENDIF
5388 GO TO 200
5389 ENDIF
5390
5391 IF(N12.EQ.8)THEN
5392 IF(iabs(LB(I1)).EQ.1)THEN
5393 LB(I2)=2
5394 LB(I1)=11
5395 E(I1)=DM
5396 ELSE
5397 LB(I1)=2
5398 LB(I2)=11
5399 E(I2)=DM
5400 ENDIF
5401 GO TO 200
5402 ENDIF
5403
5404 IF(N12.EQ.9)THEN
5405 IF(RANART(NSEED).le.0.5)THEN
5406 LB(I2)=1
5407 LB(I1)=13
5408 E(I1)=DM
5409 ELSE
5410 LB(I1)=1
5411 LB(I2)=13
5412 E(I2)=DM
5413 ENDIF
5414 GO TO 200
5415 ENDIF
5416
5417 IF(N12.EQ.10)THEN
5418 IF(RANART(NSEED).le.0.5)THEN
5419 LB(I2)=2
5420 LB(I1)=12
5421 E(I1)=DM
5422 ELSE
5423 LB(I1)=2
5424 LB(I2)=12
5425 E(I2)=DM
5426 ENDIF
5427 GO TO 200
5428 ENDIF
5429
5430 IF(N12.EQ.11)THEN
5431 IF(iabs(LB(I1)).EQ.2)THEN
5432 LB(I1)=2
5433 LB(I2)=13
5434 E(I2)=DM
5435 ELSE
5436 LB(I2)=2
5437 LB(I1)=13
5438 E(I1)=DM
5439 ENDIF
5440 GO TO 200
5441 ENDIF
5442
5443 IF(N12.EQ.12)THEN
5444 IF(iabs(LB(I1)).EQ.1)THEN
5445 LB(I1)=1
5446 LB(I2)=12
5447 E(I2)=DM
5448 ELSE
5449 LB(I2)=1
5450 LB(I1)=12
5451 E(I1)=DM
5452 ENDIF
5453 ENDIF
5454 endif
5455
5456
5457 200 EM1=E(I1)
5458 EM2=E(I2)
5459 PR2 = (SRT**2 - EM1**2 - EM2**2)**2
5460 1 - 4.0 * (EM1*EM2)**2
5461 IF(PR2.LE.0.)PR2=1.e-09
5462 PR=SQRT(PR2)/(2.*SRT)
5463 if(srt.le.2.14)C1= 1.0 - 2.0 * RANART(NSEED)
5464 if(srt.gt.2.14.and.srt.le.2.4)c1=ang(srt,iseed)
5465 if(srt.gt.2.4)then
5466
5467
5468 xptr=0.33*pr
5469
5470 cc1=ptr(xptr,iseed)
5471
5472
5473
5474 scheck=pr**2-cc1**2
5475 if(scheck.lt.0) then
5476 write(99,*) 'scheck2: ', scheck
5477 scheck=0.
5478 endif
5479 c1=sqrt(scheck)/pr
5480
5481
5482 endif
5483 T1 = 2.0 * PI * RANART(NSEED)
5484 if(ianti.eq.1 .and. lb(i1).ge.1 .and. lb(i2).ge.1)then
5485 lb(i1) = -lb(i1)
5486 lb(i2) = -lb(i2)
5487 endif
5488 GO TO 107
5489
5490
5491 106 CONTINUE
5492 NTRY1=0
5493 123 CALL DDP2(SRT,ISEED,PX3,PY3,PZ3,DM3,PX4,PY4,PZ4,DM4,
5494 & PPX,PPY,PPZ,icou1)
5495 NTRY1=NTRY1+1
5496 if((icou1.lt.0).AND.(NTRY1.LE.40))GO TO 123
5497
5498
5499 CALL ROTATE(PX,PY,PZ,PX3,PY3,PZ3)
5500 CALL ROTATE(PX,PY,PZ,PX4,PY4,PZ4)
5501 CALL ROTATE(PX,PY,PZ,PPX,PPY,PPZ)
5502 NNN=NNN+1
5503
5504
5505 XDIR=RANART(NSEED)
5506 IF(LB(I1)*LB(I2).EQ.1)THEN
5507 IF(XDIR.Le.0.2)then
5508
5509 LPION(NNN,IRUN)=4
5510 EPION(NNN,IRUN)=AP1
5511 LB(I1)=9
5512 LB(I2)=7
5513 GO TO 205
5514 ENDIF
5515
5516 IF((XDIR.LE.0.4).AND.(XDIR.GT.0.2))THEN
5517 LPION(NNN,IRUN)=4
5518 EPION(NNN,IRUN)=AP1
5519 LB(I1)=8
5520 LB(I2)=8
5521 GO TO 205
5522 ENDIF
5523
5524 IF((XDIR.LE.0.6).AND.(XDIR.GT.0.4))THEN
5525 LPION(NNN,IRUN)=3
5526 EPION(NNN,IRUN)=AP2
5527 LB(I1)=9
5528 LB(I2)=8
5529 GO TO 205
5530 ENDIF
5531 IF((XDIR.LE.0.8).AND.(XDIR.GT.0.6))THEN
5532 LPION(NNN,IRUN)=5
5533 EPION(NNN,IRUN)=AP2
5534 LB(I1)=9
5535 LB(I2)=6
5536 GO TO 205
5537 ENDIF
5538 IF(XDIR.GT.0.8)THEN
5539 LPION(NNN,IRUN)=5
5540 EPION(NNN,IRUN)=AP2
5541 LB(I1)=7
5542 LB(I2)=8
5543 GO TO 205
5544 ENDIF
5545 ENDIF
5546
5547 IF(iabs(LB(I1)).EQ.2.AND.iabs(LB(I2)).EQ.2)THEN
5548 IF(XDIR.Le.0.2)then
5549
5550 LPION(NNN,IRUN)=4
5551 EPION(NNN,IRUN)=AP1
5552 LB(I1)=6
5553 LB(I2)=7
5554 GO TO 205
5555 ENDIF
5556
5557 IF((XDIR.LE.0.4).AND.(XDIR.GT.0.2))THEN
5558 LPION(NNN,IRUN)=3
5559 EPION(NNN,IRUN)=AP2
5560 LB(I1)=6
5561 LB(I2)=9
5562 GO TO 205
5563 ENDIF
5564
5565 IF((XDIR.GT.0.4).AND.(XDIR.LE.0.6))THEN
5566 LPION(NNN,IRUN)=5
5567 EPION(NNN,IRUN)=AP2
5568 LB(I1)=9
5569 LB(I2)=8
5570 GO TO 205
5571 ENDIF
5572
5573 IF((XDIR.GT.0.6).AND.(XDIR.LE.0.8))THEN
5574 LPION(NNN,IRUN)=4
5575 EPION(NNN,IRUN)=AP1
5576 LB(I1)=7
5577 LB(I2)=7
5578 GO TO 205
5579 ENDIF
5580
5581 IF(XDIR.GT.0.8)THEN
5582 LPION(NNN,IRUN)=3
5583 EPION(NNN,IRUN)=AP2
5584 LB(I1)=7
5585 LB(I2)=8
5586 GO TO 205
5587 ENDIF
5588 ENDIF
5589
5590 IF(LB(I1)*LB(I2).EQ.2)THEN
5591 IF(XDIR.Le.0.17)then
5592
5593 LPION(NNN,IRUN)=4
5594 EPION(NNN,IRUN)=AP1
5595 LB(I1)=6
5596 LB(I2)=9
5597 GO TO 205
5598 ENDIF
5599
5600 IF((XDIR.LE.0.34).AND.(XDIR.GT.0.17))THEN
5601 LPION(NNN,IRUN)=3
5602 EPION(NNN,IRUN)=AP2
5603 LB(I1)=7
5604 LB(I2)=9
5605 GO TO 205
5606 ENDIF
5607
5608 IF((XDIR.GT.0.34).AND.(XDIR.LE.0.51))THEN
5609 LPION(NNN,IRUN)=5
5610 EPION(NNN,IRUN)=AP2
5611 LB(I1)=7
5612 LB(I2)=8
5613 GO TO 205
5614 ENDIF
5615
5616 IF((XDIR.GT.0.51).AND.(XDIR.LE.0.68))THEN
5617 LPION(NNN,IRUN)=3
5618 EPION(NNN,IRUN)=AP2
5619 LB(I1)=8
5620 LB(I2)=8
5621 GO TO 205
5622 ENDIF
5623
5624 IF((XDIR.GT.0.68).AND.(XDIR.LE.0.85))THEN
5625 LPION(NNN,IRUN)=4
5626 EPION(NNN,IRUN)=AP2
5627 LB(I1)=7
5628 LB(I2)=8
5629 GO TO 205
5630 ENDIF
5631
5632 IF(XDIR.GT.0.85)THEN
5633 LPION(NNN,IRUN)=5
5634 EPION(NNN,IRUN)=AP2
5635 LB(I1)=7
5636 LB(I2)=7
5637 ENDIF
5638 ENDIF
5639
5640
5641
5642 205 E1CM = SQRT (dm3**2 + PX3**2 + PY3**2 + PZ3**2)
5643 P1BETA = PX3*BETAX + PY3*BETAY + PZ3*BETAZ
5644 TRANSF = GAMMA * ( GAMMA * P1BETA / (GAMMA + 1) + E1CM )
5645 Pt1i1 = BETAX * TRANSF + PX3
5646 Pt2i1 = BETAY * TRANSF + PY3
5647 Pt3i1 = BETAZ * TRANSF + PZ3
5648 Eti1 = DM3
5649
5650 if(ianti.eq.1 .and. lb(i1).ge.1 .and. lb(i2).ge.1)then
5651 lb(i1) = -lb(i1)
5652 lb(i2) = -lb(i2)
5653 if(LPION(NNN,IRUN) .eq. 3)then
5654 LPION(NNN,IRUN)=5
5655 elseif(LPION(NNN,IRUN) .eq. 5)then
5656 LPION(NNN,IRUN)=3
5657 endif
5658 endif
5659
5660 lb1=lb(i1)
5661
5662 E2CM = SQRT (dm4**2 + PX4**2 + PY4**2 + PZ4**2)
5663 P2BETA = PX4*BETAX+PY4*BETAY+PZ4*BETAZ
5664 TRANSF = GAMMA * (GAMMA*P2BETA / (GAMMA + 1.) + E2CM)
5665 Pt1I2 = BETAX * TRANSF + PX4
5666 Pt2I2 = BETAY * TRANSF + PY4
5667 Pt3I2 = BETAZ * TRANSF + PZ4
5668 EtI2 = DM4
5669 lb2=lb(i2)
5670
5671
5672
5673 p(1,i1)=pt1i1
5674 p(2,i1)=pt2i1
5675 p(3,i1)=pt3i1
5676 e(i1)=eti1
5677 lb(i1)=lb1
5678 p(1,i2)=pt1i2
5679 p(2,i2)=pt2i2
5680 p(3,i2)=pt3i2
5681 e(i2)=eti2
5682 lb(i2)=lb2
5683 PX1 = P(1,I1)
5684 PY1 = P(2,I1)
5685 PZ1 = P(3,I1)
5686 EM1 = E(I1)
5687 ID(I1) = 2
5688 ID(I2) = 2
5689 ID1 = ID(I1)
5690 IBLOCK=4
5691
5692 EPCM=SQRT(EPION(NNN,IRUN)**2+PPX**2+PPY**2+PPZ**2)
5693 PPBETA=PPX*BETAX+PPY*BETAY+PPZ*BETAZ
5694 TRANSF=GAMMA*(GAMMA*PPBETA/(GAMMA+1.)+EPCM)
5695 PPION(1,NNN,IRUN)=BETAX*TRANSF+PPX
5696 PPION(2,NNN,IRUN)=BETAY*TRANSF+PPY
5697 PPION(3,NNN,IRUN)=BETAZ*TRANSF+PPZ
5698
5699 dppion(nnn,irun)=dpertp(i1)*dpertp(i2)
5700
5701
5702
5703
5704
5705
5706
5707
5708
5709 RPION(1,NNN,IRUN)=R(1,I1)
5710 RPION(2,NNN,IRUN)=R(2,I1)
5711 RPION(3,NNN,IRUN)=R(3,I1)
5712
5713 go to 90005
5714
5715
5716 108 CONTINUE
5717 if(idpert.eq.1.and.ipert1.eq.1.and.npertd.ge.1) then
5718
5719 ndloop=npertd
5720 elseif(idpert.eq.2.and.npertd.ge.1) then
5721
5722
5723
5724 ndloop=npertd+1
5725 else
5726
5727 ndloop=1
5728 endif
5729
5730 dprob1=sdprod/sig/float(npertd)
5731 do idloop=1,ndloop
5732 CALL bbdangle(pxd,pyd,pzd,nt,ipert1,ianti,idloop,pfinal,
5733 1 dprob1,lbm)
5734 CALL ROTATE(PX,PY,PZ,PXd,PYd,PZd)
5735
5736
5737
5738 xmass=xmd
5739 E1dCM=SQRT(xmass**2+PXd**2+PYd**2+PZd**2)
5740 P1dBETA=PXd*BETAX+PYd*BETAY+PZd*BETAZ
5741 TRANSF=GAMMA*(GAMMA*P1dBETA/(GAMMA+1.)+E1dCM)
5742 pxi1=BETAX*TRANSF+PXd
5743 pyi1=BETAY*TRANSF+PYd
5744 pzi1=BETAZ*TRANSF+PZd
5745 if(ianti.eq.0)then
5746 lbd=42
5747 else
5748 lbd=-42
5749 endif
5750 if(idpert.eq.1.and.ipert1.eq.1.and.npertd.ge.1) then
5751
5752 nnn=nnn+1
5753 PPION(1,NNN,IRUN)=pxi1
5754 PPION(2,NNN,IRUN)=pyi1
5755 PPION(3,NNN,IRUN)=pzi1
5756 EPION(NNN,IRUN)=xmd
5757 LPION(NNN,IRUN)=lbd
5758 RPION(1,NNN,IRUN)=R(1,I1)
5759 RPION(2,NNN,IRUN)=R(2,I1)
5760 RPION(3,NNN,IRUN)=R(3,I1)
5761
5762 dppion(NNN,IRUN)=sdprod/sig/float(npertd)
5763 elseif(idpert.eq.2.and.idloop.le.npertd) then
5764
5765
5766
5767 ppd(1,idloop)=pxi1
5768 ppd(2,idloop)=pyi1
5769 ppd(3,idloop)=pzi1
5770 lbpd(idloop)=lbd
5771 else
5772
5773
5774 E(i1)=xmm
5775 E2piCM=SQRT(xmm**2+PXd**2+PYd**2+PZd**2)
5776 P2piBETA=-PXd*BETAX-PYd*BETAY-PZd*BETAZ
5777 TRANSF=GAMMA*(GAMMA*P2piBETA/(GAMMA+1.)+E2piCM)
5778 pxi2=BETAX*TRANSF-PXd
5779 pyi2=BETAY*TRANSF-PYd
5780 pzi2=BETAZ*TRANSF-PZd
5781 p(1,i1)=pxi2
5782 p(2,i1)=pyi2
5783 p(3,i1)=pzi2
5784
5785
5786
5787
5788 LB(I1)=lbm
5789 PX1=P(1,I1)
5790 PY1=P(2,I1)
5791 PZ1=P(3,I1)
5792 EM1=E(I1)
5793 ID(I1)=2
5794 ID1=ID(I1)
5795 E1=SQRT(EM1**2+PX1**2+PY1**2+PZ1**2)
5796 lb1=lb(i1)
5797
5798 p(1,i2)=pxi1
5799 p(2,i2)=pyi1
5800 p(3,i2)=pzi1
5801 lb(i2)=lbd
5802 lb2=lb(i2)
5803 E(i2)=xmd
5804 EtI2=E(I2)
5805 ID(I2)=2
5806
5807 if(idpert.eq.2.and.idloop.eq.ndloop) then
5808 do ipertd=1,npertd
5809 nnn=nnn+1
5810 PPION(1,NNN,IRUN)=ppd(1,ipertd)
5811 PPION(2,NNN,IRUN)=ppd(2,ipertd)
5812 PPION(3,NNN,IRUN)=ppd(3,ipertd)
5813 EPION(NNN,IRUN)=xmd
5814 LPION(NNN,IRUN)=lbpd(ipertd)
5815 RPION(1,NNN,IRUN)=R(1,I1)
5816 RPION(2,NNN,IRUN)=R(2,I1)
5817 RPION(3,NNN,IRUN)=R(3,I1)
5818
5819 dppion(NNN,IRUN)=1./float(npertd)
5820 enddo
5821 endif
5822 endif
5823 enddo
5824 IBLOCK=501
5825 go to 90005
5826
5827
5828
5829 306 CONTINUE
5830
5831 if(XSK5/sigK.gt.RANART(NSEED))then
5832 pz1=p(3,i1)
5833 pz2=p(3,i2)
5834 LB(I1) = 1 + int(2 * RANART(NSEED))
5835 LB(I2) = 1 + int(2 * RANART(NSEED))
5836 nnn=nnn+1
5837 LPION(NNN,IRUN)=29
5838 EPION(NNN,IRUN)=APHI
5839 iblock = 222
5840 GO TO 208
5841 ENDIF
5842
5843 IBLOCK=9
5844 if(ianti .eq. 1)iblock=-9
5845
5846 pz1=p(3,i1)
5847 pz2=p(3,i2)
5848
5849 nnn=nnn+1
5850 LPION(NNN,IRUN)=23
5851 EPION(NNN,IRUN)=Aka
5852 if(srt.le.2.63)then
5853
5854
5855 ic=1
5856 LB(I1) = 1 + int(2 * RANART(NSEED))
5857 LB(I2)=14
5858 GO TO 208
5859 ENDIF
5860 if(srt.le.2.74.and.srt.gt.2.63)then
5861
5862 if(XSK1/(XSK1+XSK2).gt.RANART(NSEED))then
5863
5864 ic=1
5865 LB(I1) = 1 + int(2 * RANART(NSEED))
5866 LB(I2)=14
5867 else
5868
5869 LB(I1) = 1 + int(2 * RANART(NSEED))
5870 LB(I2) = 15 + int(3 * RANART(NSEED))
5871 ic=2
5872 endif
5873 GO TO 208
5874 endif
5875 if(srt.le.2.77.and.srt.gt.2.74)then
5876
5877 if(xsk1/(xsk1+xsk2+xsk3).
5878 1 gt.RANART(NSEED))then
5879
5880 ic=1
5881 LB(I1) = 1 + int(2 * RANART(NSEED))
5882 LB(I2)=14
5883 go to 208
5884 else
5885 if(xsk2/(xsk2+xsk3).gt.RANART(NSEED))then
5886
5887 ic=2
5888 LB(I1) = 1 + int(2 * RANART(NSEED))
5889 LB(I2) = 15 + int(3 * RANART(NSEED))
5890 else
5891
5892 ic=3
5893 LB(I1) = 6 + int(4 * RANART(NSEED))
5894 lb(i2)=14
5895 endif
5896 GO TO 208
5897 endif
5898 endif
5899 if(srt.gt.2.77)then
5900
5901 if(xsk1/(xsk1+xsk2+xsk3+xsk4).gt.RANART(NSEED))then
5902
5903 ic=1
5904 LB(I1) = 1 + int(2 * RANART(NSEED))
5905 LB(I2)=14
5906 go to 208
5907 else
5908 if(xsk3/(xsk2+xsk3+xsk4).gt.RANART(NSEED))then
5909
5910 ic=3
5911 LB(I1) = 6 + int(4 * RANART(NSEED))
5912 lb(i2)=14
5913 go to 208
5914 else
5915 if(xsk2/(xsk2+xsk4).gt.RANART(NSEED))then
5916
5917 LB(I1) = 1 + int(2 * RANART(NSEED))
5918 LB(I2) = 15 + int(3 * RANART(NSEED))
5919 ic=2
5920 else
5921 ic=4
5922 LB(I1) = 6 + int(4 * RANART(NSEED))
5923 LB(I2) = 15 + int(3 * RANART(NSEED))
5924 endif
5925 go to 208
5926 endif
5927 endif
5928 endif
5929 208 continue
5930 if(ianti.eq.1 .and. lb(i1).ge.1 .and. lb(i2).ge.1)then
5931 lb(i1) = - lb(i1)
5932 lb(i2) = - lb(i2)
5933 if(LPION(NNN,IRUN) .eq. 23)LPION(NNN,IRUN)=21
5934 endif
5935
5936 NTRY1=0
5937 127 CALL BBKAON(ic,SRT,PX3,PY3,PZ3,DM3,PX4,PY4,PZ4,DM4,
5938 & PPX,PPY,PPZ,icou1)
5939 NTRY1=NTRY1+1
5940 if((icou1.lt.0).AND.(NTRY1.LE.20))GO TO 127
5941
5942
5943 CALL ROTATE(PX,PY,PZ,PX3,PY3,PZ3)
5944 CALL ROTATE(PX,PY,PZ,PX4,PY4,PZ4)
5945 CALL ROTATE(PX,PY,PZ,PPX,PPY,PPZ)
5946
5947
5948
5949
5950 E1CM = SQRT (dm3**2 + PX3**2 + PY3**2 + PZ3**2)
5951 P1BETA = PX3*BETAX + PY3*BETAY + PZ3*BETAZ
5952 TRANSF = GAMMA * ( GAMMA * P1BETA / (GAMMA + 1) + E1CM )
5953 Pt1i1 = BETAX * TRANSF + PX3
5954 Pt2i1 = BETAY * TRANSF + PY3
5955 Pt3i1 = BETAZ * TRANSF + PZ3
5956 Eti1 = DM3
5957 lbi1=lb(i1)
5958
5959 E2CM = SQRT (dm4**2 + PX4**2 + PY4**2 + PZ4**2)
5960 P2BETA = PX4*BETAX+PY4*BETAY+PZ4*BETAZ
5961 TRANSF = GAMMA * (GAMMA*P2BETA / (GAMMA + 1.) + E2CM)
5962 Pt1I2 = BETAX * TRANSF + PX4
5963 Pt2I2 = BETAY * TRANSF + PY4
5964 Pt3I2 = BETAZ * TRANSF + PZ4
5965 EtI2 = DM4
5966 lbi2=lb(i2)
5967
5968 EPCM=SQRT(aka**2+PPX**2+PPY**2+PPZ**2)
5969 PPBETA=PPX*BETAX+PPY*BETAY+PPZ*BETAZ
5970 TRANSF=GAMMA*(GAMMA*PPBETA/(GAMMA+1.)+EPCM)
5971 PPION(1,NNN,IRUN)=BETAX*TRANSF+PPX
5972 PPION(2,NNN,IRUN)=BETAY*TRANSF+PPY
5973 PPION(3,NNN,IRUN)=BETAZ*TRANSF+PPZ
5974
5975 dppion(nnn,irun)=dpertp(i1)*dpertp(i2)
5976
5977
5978
5979
5980
5981
5982
5983
5984 RPION(1,NNN,IRUN)=R(1,I1)
5985 RPION(2,NNN,IRUN)=R(2,I1)
5986 RPION(3,NNN,IRUN)=R(3,I1)
5987
5988
5989
5990
5991 p(1,i1)=pt1i1
5992 p(2,i1)=pt2i1
5993 p(3,i1)=pt3i1
5994 e(i1)=eti1
5995 lb(i1)=lbi1
5996 p(1,i2)=pt1i2
5997 p(2,i2)=pt2i2
5998 p(3,i2)=pt3i2
5999 e(i2)=eti2
6000 lb(i2)=lbi2
6001 PX1 = P(1,I1)
6002 PY1 = P(2,I1)
6003 PZ1 = P(3,I1)
6004 EM1 = E(I1)
6005 ID(I1) = 2
6006 ID(I2) = 2
6007 ID1 = ID(I1)
6008 go to 90005
6009
6010
6011 307 CONTINUE
6012 NTRY1=0
6013 125 CALL DDrho(SRT,ISEED,PX3,PY3,PZ3,DM3,PX4,PY4,PZ4,DM4,
6014 & PPX,PPY,PPZ,amrho,icou1)
6015 NTRY1=NTRY1+1
6016 if((icou1.lt.0).AND.(NTRY1.LE.20))GO TO 125
6017
6018
6019 CALL ROTATE(PX,PY,PZ,PX3,PY3,PZ3)
6020 CALL ROTATE(PX,PY,PZ,PX4,PY4,PZ4)
6021 CALL ROTATE(PX,PY,PZ,PPX,PPY,PPZ)
6022 NNN=NNN+1
6023 arho=amrho
6024
6025
6026 XDIR=RANART(NSEED)
6027 IF(LB(I1)*LB(I2).EQ.1)THEN
6028 IF(XDIR.Le.0.2)then
6029
6030 LPION(NNN,IRUN)=26
6031 EPION(NNN,IRUN)=Arho
6032 LB(I1)=9
6033 LB(I2)=7
6034 GO TO 2051
6035 ENDIF
6036
6037 IF((XDIR.LE.0.4).AND.(XDIR.GT.0.2))THEN
6038 LPION(NNN,IRUN)=26
6039 EPION(NNN,IRUN)=Arho
6040 LB(I1)=8
6041 LB(I2)=8
6042 GO TO 2051
6043 ENDIF
6044
6045 IF((XDIR.LE.0.6).AND.(XDIR.GT.0.4))THEN
6046 LPION(NNN,IRUN)=25
6047 EPION(NNN,IRUN)=Arho
6048 LB(I1)=9
6049 LB(I2)=8
6050 GO TO 2051
6051 ENDIF
6052 IF((XDIR.LE.0.8).AND.(XDIR.GT.0.6))THEN
6053 LPION(NNN,IRUN)=27
6054 EPION(NNN,IRUN)=Arho
6055 LB(I1)=9
6056 LB(I2)=6
6057 GO TO 2051
6058 ENDIF
6059 IF(XDIR.GT.0.8)THEN
6060 LPION(NNN,IRUN)=27
6061 EPION(NNN,IRUN)=Arho
6062 LB(I1)=7
6063 LB(I2)=8
6064 GO TO 2051
6065 ENDIF
6066 ENDIF
6067
6068 IF(iabs(LB(I1)).EQ.2.AND.iabs(LB(I2)).EQ.2)THEN
6069 IF(XDIR.Le.0.2)then
6070
6071 LPION(NNN,IRUN)=26
6072 EPION(NNN,IRUN)=Arho
6073 LB(I1)=6
6074 LB(I2)=7
6075 GO TO 2051
6076 ENDIF
6077
6078 IF((XDIR.LE.0.4).AND.(XDIR.GT.0.2))THEN
6079 LPION(NNN,IRUN)=25
6080 EPION(NNN,IRUN)=Arho
6081 LB(I1)=6
6082 LB(I2)=9
6083 GO TO 2051
6084 ENDIF
6085
6086 IF((XDIR.GT.0.4).AND.(XDIR.LE.0.6))THEN
6087 LPION(NNN,IRUN)=27
6088 EPION(NNN,IRUN)=Arho
6089 LB(I1)=9
6090 LB(I2)=8
6091 GO TO 2051
6092 ENDIF
6093
6094 IF((XDIR.GT.0.6).AND.(XDIR.LE.0.8))THEN
6095 LPION(NNN,IRUN)=26
6096 EPION(NNN,IRUN)=Arho
6097 LB(I1)=7
6098 LB(I2)=7
6099 GO TO 2051
6100 ENDIF
6101
6102 IF(XDIR.GT.0.8)THEN
6103 LPION(NNN,IRUN)=25
6104 EPION(NNN,IRUN)=Arho
6105 LB(I1)=7
6106 LB(I2)=8
6107 GO TO 2051
6108 ENDIF
6109 ENDIF
6110
6111 IF(LB(I1)*LB(I2).EQ.2)THEN
6112 IF(XDIR.Le.0.17)then
6113
6114 LPION(NNN,IRUN)=25
6115 EPION(NNN,IRUN)=Arho
6116 LB(I1)=6
6117 LB(I2)=9
6118 GO TO 2051
6119 ENDIF
6120
6121 IF((XDIR.LE.0.34).AND.(XDIR.GT.0.17))THEN
6122 LPION(NNN,IRUN)=25
6123 EPION(NNN,IRUN)=Arho
6124 LB(I1)=7
6125 LB(I2)=9
6126 GO TO 2051
6127 ENDIF
6128
6129 IF((XDIR.GT.0.34).AND.(XDIR.LE.0.51))THEN
6130 LPION(NNN,IRUN)=27
6131 EPION(NNN,IRUN)=Arho
6132 LB(I1)=7
6133 LB(I2)=8
6134 GO TO 2051
6135 ENDIF
6136
6137 IF((XDIR.GT.0.51).AND.(XDIR.LE.0.68))THEN
6138 LPION(NNN,IRUN)=25
6139 EPION(NNN,IRUN)=Arho
6140 LB(I1)=8
6141 LB(I2)=8
6142 GO TO 2051
6143 ENDIF
6144
6145 IF((XDIR.GT.0.68).AND.(XDIR.LE.0.85))THEN
6146 LPION(NNN,IRUN)=26
6147 EPION(NNN,IRUN)=Arho
6148 LB(I1)=7
6149 LB(I2)=8
6150 GO TO 2051
6151 ENDIF
6152
6153 IF(XDIR.GT.0.85)THEN
6154 LPION(NNN,IRUN)=27
6155 EPION(NNN,IRUN)=Arho
6156 LB(I1)=7
6157 LB(I2)=7
6158 ENDIF
6159 ENDIF
6160
6161
6162
6163 2051 E1CM = SQRT (dm3**2 + PX3**2 + PY3**2 + PZ3**2)
6164 P1BETA = PX3*BETAX + PY3*BETAY + PZ3*BETAZ
6165 TRANSF = GAMMA * ( GAMMA * P1BETA / (GAMMA + 1) + E1CM )
6166 Pt1i1 = BETAX * TRANSF + PX3
6167 Pt2i1 = BETAY * TRANSF + PY3
6168 Pt3i1 = BETAZ * TRANSF + PZ3
6169 Eti1 = DM3
6170
6171 if(ianti.eq.1 .and. lb(i1).ge.1 .and. lb(i2).ge.1)then
6172 lb(i1) = -lb(i1)
6173 lb(i2) = -lb(i2)
6174 if(LPION(NNN,IRUN) .eq. 25)then
6175 LPION(NNN,IRUN)=27
6176 elseif(LPION(NNN,IRUN) .eq. 27)then
6177 LPION(NNN,IRUN)=25
6178 endif
6179 endif
6180
6181 lb1=lb(i1)
6182
6183 E2CM = SQRT (dm4**2 + PX4**2 + PY4**2 + PZ4**2)
6184 P2BETA = PX4*BETAX+PY4*BETAY+PZ4*BETAZ
6185 TRANSF = GAMMA * (GAMMA*P2BETA / (GAMMA + 1.) + E2CM)
6186 Pt1I2 = BETAX * TRANSF + PX4
6187 Pt2I2 = BETAY * TRANSF + PY4
6188 Pt3I2 = BETAZ * TRANSF + PZ4
6189 EtI2 = DM4
6190 lb2=lb(i2)
6191
6192
6193
6194 p(1,i1)=pt1i1
6195 p(2,i1)=pt2i1
6196 p(3,i1)=pt3i1
6197 e(i1)=eti1
6198 lb(i1)=lb1
6199 p(1,i2)=pt1i2
6200 p(2,i2)=pt2i2
6201 p(3,i2)=pt3i2
6202 e(i2)=eti2
6203 lb(i2)=lb2
6204 PX1 = P(1,I1)
6205 PY1 = P(2,I1)
6206 PZ1 = P(3,I1)
6207 EM1 = E(I1)
6208 ID(I1) = 2
6209 ID(I2) = 2
6210 ID1 = ID(I1)
6211 IBLOCK=44
6212
6213 EPCM=SQRT(EPION(NNN,IRUN)**2+PPX**2+PPY**2+PPZ**2)
6214 PPBETA=PPX*BETAX+PPY*BETAY+PPZ*BETAZ
6215 TRANSF=GAMMA*(GAMMA*PPBETA/(GAMMA+1.)+EPCM)
6216 PPION(1,NNN,IRUN)=BETAX*TRANSF+PPX
6217 PPION(2,NNN,IRUN)=BETAY*TRANSF+PPY
6218 PPION(3,NNN,IRUN)=BETAZ*TRANSF+PPZ
6219
6220 dppion(nnn,irun)=dpertp(i1)*dpertp(i2)
6221
6222
6223
6224
6225
6226
6227
6228
6229 RPION(1,NNN,IRUN)=R(1,I1)
6230 RPION(2,NNN,IRUN)=R(2,I1)
6231 RPION(3,NNN,IRUN)=R(3,I1)
6232
6233 go to 90005
6234
6235
6236 308 CONTINUE
6237 NTRY1=0
6238 126 CALL pprho(SRT,ISEED,PX3,PY3,PZ3,DM3,PX4,PY4,PZ4,DM4,
6239 & PPX,PPY,PPZ,amrho,icou1)
6240 NTRY1=NTRY1+1
6241 if((icou1.lt.0).AND.(NTRY1.LE.20))GO TO 126
6242
6243
6244 CALL ROTATE(PX,PY,PZ,PX3,PY3,PZ3)
6245 CALL ROTATE(PX,PY,PZ,PX4,PY4,PZ4)
6246 CALL ROTATE(PX,PY,PZ,PPX,PPY,PPZ)
6247 NNN=NNN+1
6248 arho=amrho
6249
6250
6251 XDIR=RANART(NSEED)
6252 IF(LB(I1)*LB(I2).EQ.1)THEN
6253 IF(XDIR.Le.0.5)then
6254
6255 LPION(NNN,IRUN)=26
6256 EPION(NNN,IRUN)=Arho
6257 LB(I1)=1
6258 LB(I2)=1
6259 GO TO 2052
6260 Else
6261
6262 LPION(NNN,IRUN)=27
6263 EPION(NNN,IRUN)=Arho
6264 LB(I1)=1
6265 LB(I2)=2
6266 GO TO 2052
6267 ENDIF
6268 endif
6269
6270 IF(iabs(LB(I1)).EQ.2.AND.iabs(LB(I2)).EQ.2)THEN
6271 IF(XDIR.Le.0.5)then
6272
6273 LPION(NNN,IRUN)=26
6274 EPION(NNN,IRUN)=Arho
6275 LB(I1)=2
6276 LB(I2)=2
6277 GO TO 2052
6278 Else
6279
6280 LPION(NNN,IRUN)=25
6281 EPION(NNN,IRUN)=Arho
6282 LB(I1)=1
6283 LB(I2)=2
6284 GO TO 2052
6285 ENDIF
6286 endif
6287
6288 IF(LB(I1)*LB(I2).EQ.2)THEN
6289 IF(XDIR.Le.0.33)then
6290
6291 LPION(NNN,IRUN)=26
6292 EPION(NNN,IRUN)=Arho
6293 LB(I1)=1
6294 LB(I2)=2
6295 GO TO 2052
6296
6297 else IF((XDIR.LE.0.67).AND.(XDIR.GT.0.34))THEN
6298 LPION(NNN,IRUN)=25
6299 EPION(NNN,IRUN)=Arho
6300 LB(I1)=1
6301 LB(I2)=1
6302 GO TO 2052
6303 Else
6304
6305 LPION(NNN,IRUN)=27
6306 EPION(NNN,IRUN)=Arho
6307 LB(I1)=2
6308 LB(I2)=2
6309 GO TO 2052
6310 ENDIF
6311 endif
6312
6313
6314
6315 2052 E1CM = SQRT (dm3**2 + PX3**2 + PY3**2 + PZ3**2)
6316 P1BETA = PX3*BETAX + PY3*BETAY + PZ3*BETAZ
6317 TRANSF = GAMMA * ( GAMMA * P1BETA / (GAMMA + 1) + E1CM )
6318 Pt1i1 = BETAX * TRANSF + PX3
6319 Pt2i1 = BETAY * TRANSF + PY3
6320 Pt3i1 = BETAZ * TRANSF + PZ3
6321 Eti1 = DM3
6322
6323 if(ianti.eq.1 .and. lb(i1).ge.1 .and. lb(i2).ge.1)then
6324 lb(i1) = -lb(i1)
6325 lb(i2) = -lb(i2)
6326 if(LPION(NNN,IRUN) .eq. 25)then
6327 LPION(NNN,IRUN)=27
6328 elseif(LPION(NNN,IRUN) .eq. 27)then
6329 LPION(NNN,IRUN)=25
6330 endif
6331 endif
6332
6333 lb1=lb(i1)
6334
6335 E2CM = SQRT (dm4**2 + PX4**2 + PY4**2 + PZ4**2)
6336 P2BETA = PX4*BETAX+PY4*BETAY+PZ4*BETAZ
6337 TRANSF = GAMMA * (GAMMA*P2BETA / (GAMMA + 1.) + E2CM)
6338 Pt1I2 = BETAX * TRANSF + PX4
6339 Pt2I2 = BETAY * TRANSF + PY4
6340 Pt3I2 = BETAZ * TRANSF + PZ4
6341 EtI2 = DM4
6342 lb2=lb(i2)
6343
6344
6345
6346 p(1,i1)=pt1i1
6347 p(2,i1)=pt2i1
6348 p(3,i1)=pt3i1
6349 e(i1)=eti1
6350 lb(i1)=lb1
6351 p(1,i2)=pt1i2
6352 p(2,i2)=pt2i2
6353 p(3,i2)=pt3i2
6354 e(i2)=eti2
6355 lb(i2)=lb2
6356 PX1 = P(1,I1)
6357 PY1 = P(2,I1)
6358 PZ1 = P(3,I1)
6359 EM1 = E(I1)
6360 ID(I1) = 2
6361 ID(I2) = 2
6362 ID1 = ID(I1)
6363 IBLOCK=45
6364
6365 EPCM=SQRT(EPION(NNN,IRUN)**2+PPX**2+PPY**2+PPZ**2)
6366 PPBETA=PPX*BETAX+PPY*BETAY+PPZ*BETAZ
6367 TRANSF=GAMMA*(GAMMA*PPBETA/(GAMMA+1.)+EPCM)
6368 PPION(1,NNN,IRUN)=BETAX*TRANSF+PPX
6369 PPION(2,NNN,IRUN)=BETAY*TRANSF+PPY
6370 PPION(3,NNN,IRUN)=BETAZ*TRANSF+PPZ
6371
6372 dppion(nnn,irun)=dpertp(i1)*dpertp(i2)
6373
6374
6375
6376
6377
6378
6379
6380
6381 RPION(1,NNN,IRUN)=R(1,I1)
6382 RPION(2,NNN,IRUN)=R(2,I1)
6383 RPION(3,NNN,IRUN)=R(3,I1)
6384
6385 go to 90005
6386
6387
6388 309 CONTINUE
6389 NTRY1=0
6390 138 CALL ppomga(SRT,ISEED,PX3,PY3,PZ3,DM3,PX4,PY4,PZ4,DM4,
6391 & PPX,PPY,PPZ,icou1)
6392 NTRY1=NTRY1+1
6393 if((icou1.lt.0).AND.(NTRY1.LE.20))GO TO 138
6394
6395
6396 CALL ROTATE(PX,PY,PZ,PX3,PY3,PZ3)
6397 CALL ROTATE(PX,PY,PZ,PX4,PY4,PZ4)
6398 CALL ROTATE(PX,PY,PZ,PPX,PPY,PPZ)
6399 NNN=NNN+1
6400 aomega=0.782
6401
6402
6403 IF(LB(I1)*LB(I2).EQ.1)THEN
6404
6405 LPION(NNN,IRUN)=28
6406 EPION(NNN,IRUN)=Aomega
6407 LB(I1)=1
6408 LB(I2)=1
6409 GO TO 2053
6410 ENDIF
6411
6412 IF(iabs(LB(I1)).EQ.2.AND.iabs(LB(I2)).EQ.2)THEN
6413
6414 LPION(NNN,IRUN)=28
6415 EPION(NNN,IRUN)=Aomega
6416 LB(I1)=2
6417 LB(I2)=2
6418 GO TO 2053
6419 ENDIF
6420
6421 IF(LB(I1)*LB(I2).EQ.2)THEN
6422
6423 LPION(NNN,IRUN)=28
6424 EPION(NNN,IRUN)=Aomega
6425 LB(I1)=1
6426 LB(I2)=2
6427 GO TO 2053
6428 ENDIF
6429
6430
6431
6432 2053 E1CM = SQRT (dm3**2 + PX3**2 + PY3**2 + PZ3**2)
6433 P1BETA = PX3*BETAX + PY3*BETAY + PZ3*BETAZ
6434 TRANSF = GAMMA * ( GAMMA * P1BETA / (GAMMA + 1) + E1CM )
6435 Pt1i1 = BETAX * TRANSF + PX3
6436 Pt2i1 = BETAY * TRANSF + PY3
6437 Pt3i1 = BETAZ * TRANSF + PZ3
6438 Eti1 = DM3
6439 if(ianti.eq.1 .and. lb(i1).ge.1 .and. lb(i2).ge.1)then
6440 lb(i1) = -lb(i1)
6441 lb(i2) = -lb(i2)
6442 endif
6443 lb1=lb(i1)
6444
6445 E2CM = SQRT (dm4**2 + PX4**2 + PY4**2 + PZ4**2)
6446 P2BETA = PX4*BETAX+PY4*BETAY+PZ4*BETAZ
6447 TRANSF = GAMMA * (GAMMA*P2BETA / (GAMMA + 1.) + E2CM)
6448 Pt1I2 = BETAX * TRANSF + PX4
6449 Pt2I2 = BETAY * TRANSF + PY4
6450 Pt3I2 = BETAZ * TRANSF + PZ4
6451 EtI2 = DM4
6452 lb2=lb(i2)
6453
6454
6455
6456 p(1,i1)=pt1i1
6457 p(2,i1)=pt2i1
6458 p(3,i1)=pt3i1
6459 e(i1)=eti1
6460 lb(i1)=lb1
6461 p(1,i2)=pt1i2
6462 p(2,i2)=pt2i2
6463 p(3,i2)=pt3i2
6464 e(i2)=eti2
6465 lb(i2)=lb2
6466 PX1 = P(1,I1)
6467 PY1 = P(2,I1)
6468 PZ1 = P(3,I1)
6469 EM1 = E(I1)
6470 ID(I1) = 2
6471 ID(I2) = 2
6472 ID1 = ID(I1)
6473 IBLOCK=46
6474
6475 EPCM=SQRT(EPION(NNN,IRUN)**2+PPX**2+PPY**2+PPZ**2)
6476 PPBETA=PPX*BETAX+PPY*BETAY+PPZ*BETAZ
6477 TRANSF=GAMMA*(GAMMA*PPBETA/(GAMMA+1.)+EPCM)
6478 PPION(1,NNN,IRUN)=BETAX*TRANSF+PPX
6479 PPION(2,NNN,IRUN)=BETAY*TRANSF+PPY
6480 PPION(3,NNN,IRUN)=BETAZ*TRANSF+PPZ
6481
6482 dppion(nnn,irun)=dpertp(i1)*dpertp(i2)
6483
6484
6485
6486
6487
6488
6489
6490
6491 RPION(1,NNN,IRUN)=R(1,I1)
6492 RPION(2,NNN,IRUN)=R(2,I1)
6493 RPION(3,NNN,IRUN)=R(3,I1)
6494
6495 go to 90005
6496
6497
6498
6499
6500
6501
6502
6503
6504
6505
6506
6507
6508
6509
6510
6511
6512
6513
6514
6515
6516
6517
6518
6519
6520
6521
6522
6523
6524
6525
6526
6527
6528
6529
6530
6531
6532
6533
6534
6535
6536
6537 90005 continue
6538 RETURN
6539
6540
6541 107 IF(PX .EQ. 0.0 .AND. PY .EQ. 0.0) THEN
6542 T2 = 0.0
6543 ELSE
6544 T2=ATAN2(PY,PX)
6545 END IF
6546 S1 = 1.0 - C1**2
6547 IF(S1.LE.0)S1=0
6548 S1=SQRT(S1)
6549
6550
6551 scheck=1.0 - C2**2
6552 if(scheck.lt.0) then
6553 write(99,*) 'scheck3: ', scheck
6554 scheck=0.
6555 endif
6556 S2=SQRT(scheck)
6557
6558
6559 CT1 = COS(T1)
6560 ST1 = SIN(T1)
6561 CT2 = COS(T2)
6562 ST2 = SIN(T2)
6563 PZ = PR * ( C1*C2 - S1*S2*CT1 )
6564 SS = C2 * S1 * CT1 + S2 * C1
6565 PX = PR * ( SS*CT2 - S1*ST1*ST2 )
6566 PY = PR * ( SS*ST2 + S1*ST1*CT2 )
6567 RETURN
6568 END
6569
6570
6571
6572
6573
6574
6575
6576 SUBROUTINE CRPP(PX,PY,PZ,SRT,I1,I2,IBLOCK,
6577 &ppel,ppin,spprho,ipp)
6578
6579
6580
6581
6582
6583
6584
6585
6586
6587
6588
6589 PARAMETER (MAXSTR=150001,MAXR=1,AMN=0.939457,
6590 1 AMP=0.93828,AP1=0.13496,
6591 2 AP2=0.13957,AM0=1.232,PI=3.1415926,CUTOFF=1.8966,AVMASS=0.9383)
6592 PARAMETER (AKA=0.498,aks=0.895)
6593 parameter (MX=4,MY=4,MZ=8,MPX=4,MPY=4,mpz=10,mpzp=10)
6594 COMMON /AA/ R(3,MAXSTR)
6595
6596 COMMON /BB/ P(3,MAXSTR)
6597
6598 COMMON /CC/ E(MAXSTR)
6599
6600 COMMON /EE/ ID(MAXSTR),LB(MAXSTR)
6601
6602 common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
6603
6604 common/ppb1/ene,factr2(6),fsum,ppinnb,s,wtot
6605
6606 common/ppmm/pprr,ppee,pppe,rpre,xopoe,rree
6607
6608 COMMON/RNDF77/NSEED
6609
6610 SAVE
6611
6612 lb1i=lb(i1)
6613 lb2i=lb(i2)
6614
6615 PX0=PX
6616 PY0=PY
6617 PZ0=PZ
6618 iblock=1
6619
6620
6621
6622
6623
6624
6625
6626
6627
6628
6629
6630
6631
6632 if(srt.gt.(2*aka).and.(ppin/(ppin+ppel)).gt.RANART(NSEED)) then
6633
6634
6635
6636 ranpi=RANART(NSEED)
6637 if((pprr/ppin).ge.ranpi) then
6638
6639
6640 call pi2ro2(i1,i2,lbb1,lbb2,ei1,ei2,iblock,iseed)
6641
6642
6643 elseif((pprr+ppee)/ppin.ge.ranpi) then
6644
6645 call pi2et2(i1,i2,lbb1,lbb2,ei1,ei2,iblock,iseed)
6646 elseif(((pprr+ppee+pppe)/ppin).ge.ranpi) then
6647
6648 call pi3eta(i1,i2,lbb1,lbb2,ei1,ei2,iblock,iseed)
6649 elseif(((pprr+ppee+pppe+rpre)/ppin).ge.ranpi) then
6650
6651 call rpiret(i1,i2,lbb1,lbb2,ei1,ei2,iblock,iseed)
6652 elseif(((pprr+ppee+pppe+rpre+xopoe)/ppin).ge.ranpi) then
6653
6654 call opioet(i1,i2,lbb1,lbb2,ei1,ei2,iblock,iseed)
6655 elseif(((pprr+ppee+pppe+rpre+xopoe+rree)
6656 1 /ppin).ge.ranpi) then
6657
6658 call ro2et2(i1,i2,lbb1,lbb2,ei1,ei2,iblock,iseed)
6659
6660
6661
6662 elseif(((pprr+ppee+pppe+rpre+xopoe+rree+ppinnb)/ppin)
6663 1 .ge.ranpi) then
6664
6665 call bbarfs(lbb1,lbb2,ei1,ei2,iblock,iseed)
6666
6667 else
6668 iblock=66
6669 ei1=aka
6670 ei2=aka
6671 lbb1=21
6672 lbb2=23
6673
6674 lb1=lb(i1)
6675 lb2=lb(i2)
6676
6677
6678
6679 if( ( (lb1.eq.0.or.(lb1.ge.3.and.lb1.le.5))
6680 1 .and.(lb2.ge.25.and.lb2.le.28))
6681 2 .or. ( (lb2.eq.0.or.(lb2.ge.3.and.lb2.le.5))
6682 3 .and.(lb1.ge.25.and.lb1.le.28))) then
6683 ei1=aks
6684 ei2=aka
6685 if(RANART(NSEED).ge.0.5) then
6686 iblock=366
6687 lbb1=30
6688 lbb2=21
6689 else
6690 iblock=367
6691 lbb1=-30
6692 lbb2=23
6693 endif
6694 endif
6695
6696 endif
6697
6698 e(i1)=ei1
6699 e(i2)=ei2
6700 lb(i1)=lbb1
6701 lb(i2)=lbb2
6702
6703
6704 else
6705
6706
6707 if ((lb(i1).lt.3.or.lb(i1).gt.5).and.
6708 & (lb(i2).lt.3.or.lb(i2).gt.5)) return
6709
6710
6711
6712 IBLOCK=6
6713
6714 if(ipp.eq.1.or.ipp.eq.4.or.ipp.eq.6)go to 10
6715 if(spprho/ppel.gt.RANART(NSEED))go to 20
6716 endif
6717 10 NTAG=0
6718 EM1=E(I1)
6719 EM2=E(I2)
6720
6721
6722
6723
6724 PR2 = (SRT**2 - EM1**2 - EM2**2)**2
6725 1 - 4.0 * (EM1*EM2)**2
6726 IF(PR2.LE.0.)PR2=1.e-09
6727 PR=SQRT(PR2)/(2.*SRT)
6728 C1 = 1.0 - 2.0 * RANART(NSEED)
6729 T1 = 2.0 * PI * RANART(NSEED)
6730 S1 = SQRT( 1.0 - C1**2 )
6731 CT1 = COS(T1)
6732 ST1 = SIN(T1)
6733 PZ = PR * C1
6734 PX = PR * S1*CT1
6735 PY = PR * S1*ST1
6736
6737
6738
6739 CALL ROTATE(PX0,PY0,PZ0,PX,PY,PZ)
6740
6741 RETURN
6742 20 continue
6743 iblock=666
6744
6745
6746 call rhores(i1,i2)
6747 if(ipp.eq.2)lb(i1)=27
6748 if(ipp.eq.3)lb(i1)=26
6749 if(ipp.eq.5)lb(i1)=25
6750 return
6751 END
6752
6753
6754
6755
6756 SUBROUTINE CRND(IRUN,PX,PY,PZ,SRT,I1,I2,IBLOCK,
6757 &SIGNN,SIG,sigk,xsk1,xsk2,xsk3,xsk4,xsk5,NT,ipert1)
6758
6759
6760
6761
6762
6763
6764
6765
6766
6767
6768
6769
6770
6771
6772
6773
6774
6775
6776
6777
6778
6779
6780
6781
6782
6783
6784
6785
6786
6787
6788
6789
6790
6791
6792
6793
6794
6795
6796
6797
6798
6799
6800
6801
6802
6803
6804
6805
6806
6807
6808
6809
6810
6811
6812
6813
6814
6815 PARAMETER (MAXSTR=150001,MAXR=1,AMN=0.939457,
6816 1 AMP=0.93828,AP1=0.13496,AKA=0.498,APHI=1.020,
6817 2 AP2=0.13957,AM0=1.232,PI=3.1415926,CUTOFF=1.8966,AVMASS=0.9383)
6818 parameter (MX=4,MY=4,MZ=8,MPX=4,MPY=4,mpz=10,mpzp=10)
6819 parameter (xmd=1.8756,npdmax=10000)
6820 COMMON /AA/ R(3,MAXSTR)
6821
6822 COMMON /BB/ P(3,MAXSTR)
6823
6824 COMMON /CC/ E(MAXSTR)
6825
6826 COMMON /EE/ ID(MAXSTR),LB(MAXSTR)
6827
6828 common /ff/f(-mx:mx,-my:my,-mz:mz,-mpx:mpx,-mpy:mpy,-mpz:mpzp)
6829
6830 common /gg/ dx,dy,dz,dpx,dpy,dpz
6831
6832 COMMON /INPUT/ NSTAR,NDIRCT,DIR
6833
6834 COMMON /NN/NNN
6835
6836 COMMON /BG/BETAX,BETAY,BETAZ,GAMMA
6837
6838 COMMON /RUN/NUM
6839
6840 COMMON /PA/RPION(3,MAXSTR,MAXR)
6841
6842 COMMON /PB/PPION(3,MAXSTR,MAXR)
6843
6844 COMMON /PC/EPION(MAXSTR,MAXR)
6845
6846 COMMON /PD/LPION(MAXSTR,MAXR)
6847
6848 common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
6849
6850 common/leadng/lb1,px1,py1,pz1,em1,e1,xfnl,yfnl,zfnl,tfnl,
6851 1 px1n,py1n,pz1n,dp1n
6852
6853 COMMON/RNDF77/NSEED
6854
6855 COMMON /dpert/dpertt(MAXSTR,MAXR),dpertp(MAXSTR),dplast(MAXSTR),
6856 1 dpdcy(MAXSTR),dpdpi(MAXSTR,MAXR),dpt(MAXSTR, MAXR),
6857 2 dpp1(MAXSTR,MAXR),dppion(MAXSTR,MAXR)
6858 common /dpi/em2,lb2
6859 common /para8/ idpert,npertd,idxsec
6860 dimension ppd(3,npdmax),lbpd(npdmax)
6861 SAVE
6862
6863 n12=0
6864 m12=0
6865 IBLOCK=0
6866 NTAG=0
6867 EM1=E(I1)
6868 EM2=E(I2)
6869 PR = SQRT( PX**2 + PY**2 + PZ**2 )
6870 C2 = PZ / PR
6871 X1 = RANART(NSEED)
6872 ianti=0
6873 if(lb(i1).lt.0 .and. lb(i2).lt.0)ianti=1
6874
6875
6876 call sbbdm(srt,sdprod,ianti,lbm,xmm,pfinal)
6877 if(idpert.eq.1.and.ipert1.eq.1) then
6878 IF (SRT .LT. 2.012) RETURN
6879 if((iabs(lb(i1)).eq.1.or.iabs(lb(i1)).eq.2)
6880 1 .and.(iabs(lb(i2)).ge.6.and.iabs(lb(i2)).le.13)) then
6881 goto 108
6882 elseif((iabs(lb(i2)).eq.1.or.iabs(lb(i2)).eq.2)
6883 1 .and.(iabs(lb(i1)).ge.6.and.iabs(lb(i1)).le.13)) then
6884 goto 108
6885 else
6886 return
6887 endif
6888 endif
6889
6890
6891
6892 IF (X1 .LE. SIGNN/SIG) THEN
6893
6894 AS = ( 3.65 * (SRT - 1.8766) )**6
6895 A = 6.0 * AS / (1.0 + AS)
6896 TA = -2.0 * PR**2
6897 X = RANART(NSEED)
6898
6899 T1 = sngl(DLOG(dble(1.-X)*DEXP(dble(A)*dble(TA))+dble(X)))/ A
6900 C1 = 1.0 - T1/TA
6901 T1 = 2.0 * PI * RANART(NSEED)
6902 IBLOCK=1
6903 GO TO 107
6904 ELSE
6905
6906
6907
6908 IF (SRT .LT. 2.04) RETURN
6909
6910
6911 if(((iabs(LB(I1)).EQ.2.or.iabs(LB(I2)).EQ.2).AND.
6912 1 (LB(I1)*LB(I2)).EQ.20).or.(LB(I1)*LB(I2)).EQ.13) then
6913 IF(X1.LE.((SIGNN+sdprod)/SIG)) GO TO 108
6914 ENDIF
6915
6916
6917
6918
6919 PRF=SQRT(0.25*SRT**2-AVMASS**2)
6920 IF(EM1.GT.1.)THEN
6921 DELTAM=EM1
6922 ELSE
6923 DELTAM=EM2
6924 ENDIF
6925 RENOM=DELTAM*PRF**2/DENOM(SRT,1.)/PR
6926 RENOMN=DELTAM*PRF**2/DENOM(SRT,2.)/PR
6927 RENOM1=DELTAM*PRF**2/DENOM(SRT,-1.)/PR
6928
6929
6930
6931 if((iabs(lb(i1)).eq.2).and.(iabs(lb(i2)).eq.6)) renom=0.
6932 if((iabs(lb(i2)).eq.2).and.(iabs(lb(i1)).eq.6)) renom=0.
6933 if((iabs(lb(i1)).eq.1).and.(iabs(lb(i2)).eq.9)) renom=0.
6934 if((iabs(lb(i2)).eq.1).and.(iabs(lb(i1)).eq.9)) renom=0.
6935 Call M1535(iabs(lb(i1)),iabs(lb(i2)),srt,x1535)
6936 X1440=(3./4.)*SIGMA(SRT,2,0,1)
6937
6938
6939
6940
6941
6942 if(((iabs(lb(i1)).eq.2).and.(iabs(lb(i2)).eq.6)).OR.
6943 & ((iabs(lb(i2)).eq.2).and.(iabs(lb(i1)).eq.6)).OR.
6944 & ((iabs(lb(i1)).eq.1).and.(iabs(lb(i2)).eq.9)).OR.
6945 & ((iabs(lb(i2)).eq.1).and.(iabs(lb(i1)).eq.9)))THEN
6946
6947 IF(X1.LE.((SIGNN+sdprod)/SIG)) GO TO 108
6948
6949 IF((SIGK+SIGNN+sdprod)/SIG.GE.X1)GO TO 306
6950
6951 ENDIF
6952
6953
6954
6955 IF(LB(I1)*LB(I2).EQ.18.AND.
6956 & (iabs(LB(I1)).EQ.2.OR.iabs(LB(I2)).EQ.2))then
6957 SIGND=SIGMA(SRT,1,1,0)+0.5*SIGMA(SRT,1,1,1)
6958 SIGDN=0.25*SIGND*RENOM
6959
6960 IF(X1.LE.((SIGNN+sdprod)/SIG)) GO TO 108
6961
6962 IF(X1.GT.(SIGNN+SIGDN+X1440+X1535+SIGK+sdprod)/SIG)RETURN
6963
6964 IF(SIGK/(SIGK+SIGDN+X1440+X1535).GT.RANART(NSEED))GO TO 306
6965
6966 IF(RANART(NSEED).LT.SIGDN/(SIGDN+X1440+X1535))THEN
6967 M12=3
6968 GO TO 206
6969 ELSE
6970
6971 IF(RANART(NSEED).LT.X1440/(X1440+X1535))THEN
6972
6973 M12=37
6974 ELSE
6975
6976
6977
6978
6979 return
6980
6981 ENDIF
6982 GO TO 204
6983 ENDIF
6984 ENDIF
6985
6986
6987 IF(LB(I1)*LB(I2).EQ.6.AND.
6988 & ((iabs(LB(I1)).EQ.1).OR.(iabs(LB(I2)).EQ.1)))then
6989 SIGND=SIGMA(SRT,1,1,0)+0.5*SIGMA(SRT,1,1,1)
6990 SIGDN=0.25*SIGND*RENOM
6991
6992 IF(X1.LE.((SIGNN+sdprod)/SIG)) GO TO 108
6993
6994 IF (X1.GT.(SIGNN+SIGDN+X1440+X1535+SIGK+sdprod)/SIG)RETURN
6995
6996 IF(SIGK/(SIGK+SIGDN+X1440+X1535).GT.RANART(NSEED))GO TO 306
6997
6998 IF(RANART(NSEED).LT.SIGDN/(SIGDN+X1440+X1535))THEN
6999 M12=6
7000 GO TO 206
7001 ELSE
7002
7003 IF(RANART(NSEED).LT.X1440/(X1440+X1535))THEN
7004
7005 M12=47
7006 ELSE
7007
7008
7009 return
7010
7011 ENDIF
7012 GO TO 204
7013 ENDIF
7014 ENDIF
7015
7016 IF(LB(I1)*LB(I2).EQ.8.AND.
7017 & (iabs(LB(I1)).EQ.1.OR.iabs(LB(I2)).EQ.1))THEN
7018 SIGND=1.5*SIGMA(SRT,1,1,1)
7019 SIGDN=0.25*SIGND*RENOM
7020
7021 IF(X1.LE.((SIGNN+sdprod)/SIG)) GO TO 108
7022
7023 IF(X1.GT.(SIGNN+SIGDN+x1440+x1535+SIGK+sdprod)/SIG)RETURN
7024
7025 IF(SIGK/(SIGK+SIGDN+X1440+X1535).GT.RANART(NSEED))GO TO 306
7026 IF(RANART(NSEED).LT.SIGDN/(SIGDN+X1440+X1535))THEN
7027 M12=4
7028 GO TO 206
7029 ELSE
7030 IF(RANART(NSEED).LT.X1440/(X1440+X1535))THEN
7031
7032 M12=39
7033 ELSE
7034 M12=40
7035 ENDIF
7036 GO TO 204
7037 ENDIF
7038 ENDIF
7039
7040 IF(LB(I1)*LB(I2).EQ.14.AND.
7041 & (iabs(LB(I1)).EQ.2.OR.iabs(LB(I2)).EQ.2))THEN
7042 SIGND=1.5*SIGMA(SRT,1,1,1)
7043 SIGDN=0.25*SIGND*RENOM
7044
7045 IF(X1.LE.((SIGNN+sdprod)/SIG)) GO TO 108
7046
7047 IF(X1.GT.(SIGNN+SIGDN+x1440+x1535+SIGK+sdprod)/SIG)RETURN
7048
7049 IF(SIGK/(SIGK+SIGDN+X1440+X1535).GT.RANART(NSEED))GO TO 306
7050 IF(RANART(NSEED).LT.SIGDN/(SIGDN+X1440+X1535))THEN
7051 M12=5
7052 GO TO 206
7053 ELSE
7054 IF(RANART(NSEED).LT.X1440/(X1440+X1535))THEN
7055
7056 M12=48
7057 ELSE
7058 M12=49
7059 ENDIF
7060 GO TO 204
7061 ENDIF
7062 ENDIF
7063
7064
7065 IF(LB(I1)*LB(I2).EQ.16.AND.
7066 & (iabs(LB(I1)).EQ.2.OR.iabs(LB(I2)).EQ.2))THEN
7067 SIGND=0.5*SIGMA(SRT,1,1,1)+0.25*SIGMA(SRT,1,1,0)
7068 SIGDN=0.5*SIGND*RENOM
7069
7070 IF(X1.LE.((SIGNN+sdprod)/SIG)) GO TO 108
7071
7072 IF(X1.GT.(SIGNN+SIGDN+2.*x1440+2.*x1535+SIGK+sdprod)/SIG)RETURN
7073
7074 IF(SIGK/(SIGK+SIGDN+2*X1440+2*X1535).GT.RANART(NSEED))GO TO 306
7075 IF(RANART(NSEED).LT.SIGDN/(SIGDN+2.*X1440+2.*X1535))THEN
7076 M12=1
7077 GO TO 206
7078 ELSE
7079 IF(RANART(NSEED).LT.X1440/(X1440+X1535))THEN
7080 M12=41
7081 IF(RANART(NSEED).LE.0.5)M12=43
7082 ELSE
7083 M12=42
7084 IF(RANART(NSEED).LE.0.5)M12=44
7085 ENDIF
7086 GO TO 204
7087 ENDIF
7088 ENDIF
7089
7090
7091 IF(LB(I1)*LB(I2).EQ.7)THEN
7092 SIGND=0.5*SIGMA(SRT,1,1,1)+0.25*SIGMA(SRT,1,1,0)
7093 SIGDN=0.5*SIGND*RENOM
7094
7095 IF(X1.LE.((SIGNN+sdprod)/SIG)) GO TO 108
7096
7097 IF(X1.GT.(SIGNN+SIGDN+2.*x1440+2.*x1535+SIGK+sdprod)/SIG)RETURN
7098
7099 IF(SIGK/(SIGK+SIGDN+2*X1440+2*X1535).GT.RANART(NSEED))GO TO 306
7100 IF(RANART(NSEED).LT.SIGDN/(SIGDN+2.*X1440+2.*X1535))THEN
7101 M12=2
7102 GO TO 206
7103 ELSE
7104 IF(RANART(NSEED).LT.X1440/(X1440+X1535))THEN
7105 M12=50
7106 IF(RANART(NSEED).LE.0.5)M12=51
7107 ELSE
7108 M12=52
7109 IF(RANART(NSEED).LE.0.5)M12=53
7110 ENDIF
7111 GO TO 204
7112 ENDIF
7113 ENDIF
7114
7115
7116 IF(LB(I1)*LB(I2).EQ.10.AND.
7117 & (iabs(LB(I1)).EQ.1.OR.iabs(LB(I2)).EQ.1))then
7118 SIGND=(3./4.)*SIGMA(SRT,2,0,1)
7119 SIGDN=SIGND*RENOMN
7120
7121 IF(X1.LE.((SIGNN+sdprod)/SIG)) GO TO 108
7122
7123 IF(X1.GT.(SIGNN+SIGDN+X1535+SIGK+sdprod)/SIG)RETURN
7124
7125 IF(SIGK/(SIGK+SIGDN+X1535).GT.RANART(NSEED))GO TO 306
7126 IF(RANART(NSEED).LT.SIGDN/(SIGDN+X1535))THEN
7127 M12=7
7128 GO TO 206
7129 ELSE
7130 M12=54
7131 IF(RANART(NSEED).LE.0.5)M12=55
7132 ENDIF
7133 GO TO 204
7134 ENDIF
7135
7136 IF(LB(I1)*LB(I2).EQ.22.AND.
7137 & (iabs(LB(I1)).EQ.2.OR.iabs(LB(I2)).EQ.2))then
7138 SIGND=(3./4.)*SIGMA(SRT,2,0,1)
7139 SIGDN=SIGND*RENOMN
7140
7141 IF(X1.LE.((SIGNN+sdprod)/SIG)) GO TO 108
7142
7143 IF(X1.GT.(SIGNN+SIGDN+X1535+SIGK+sdprod)/SIG)RETURN
7144
7145 IF(SIGK/(SIGK+SIGDN+X1535).GT.RANART(NSEED))GO TO 306
7146 IF(RANART(NSEED).LT.SIGDN/(SIGDN+X1535))THEN
7147 M12=8
7148 GO TO 206
7149 ELSE
7150 M12=56
7151 IF(RANART(NSEED).LE.0.5)M12=57
7152 ENDIF
7153 GO TO 204
7154 ENDIF
7155
7156 IF((iabs(LB(I1)).EQ.12).OR.(iabs(LB(I1)).EQ.13).OR.
7157 1 (iabs(LB(I2)).EQ.12).OR.(iabs(LB(I2)).EQ.13))THEN
7158 SIGND=X1535
7159 SIGDN=SIGND*RENOM1
7160
7161 IF(X1.LE.((SIGNN+sdprod)/SIG)) GO TO 108
7162
7163 IF(X1.GT.(SIGNN+SIGDN+SIGK+sdprod)/SIG)RETURN
7164
7165 IF(SIGK/(SIGK+SIGDN).GT.RANART(NSEED))GO TO 306
7166 IF(LB(I1)*LB(I2).EQ.24)M12=10
7167 IF(LB(I1)*LB(I2).EQ.12)M12=12
7168 IF(LB(I1)*LB(I2).EQ.26)M12=11
7169 IF(LB(I1)*LB(I2).EQ.13)M12=9
7170 GO TO 206
7171 ENDIF
7172 204 CONTINUE
7173
7174
7175
7176
7177
7178
7179
7180 DMAX = SRT - AVMASS-0.005
7181 DMIN = 1.078
7182 IF((M12.eq.37).or.(M12.eq.39).or.
7183 1 (M12.eQ.41).OR.(M12.eQ.43).OR.(M12.EQ.46).
7184 2 OR.(M12.EQ.48).OR.(M12.EQ.50).OR.(M12.EQ.51))then
7185
7186 IF(DMAX.LT.1.44) THEN
7187 FM=FNS(DMAX,SRT,0.)
7188 ELSE
7189
7190
7191 xdmass=1.44
7192
7193 FM=FNS(xdmass,SRT,1.)
7194
7195
7196 ENDIF
7197 IF(FM.EQ.0.)FM=1.E-09
7198 NTRY2=0
7199 11 DM=RANART(NSEED)*(DMAX-DMIN)+DMIN
7200 NTRY2=NTRY2+1
7201 IF((RANART(NSEED).GT.FNS(DM,SRT,1.)/FM).AND.
7202 1 (NTRY2.LE.10)) GO TO 11
7203
7204
7205
7206 if(dm.gt.2.14) goto 11
7207
7208 GO TO 13
7209 ELSE
7210
7211 IF(DMAX.LT.1.535) THEN
7212 FM=FD5(DMAX,SRT,0.)
7213 ELSE
7214
7215
7216 xdmass=1.535
7217
7218 FM=FD5(xdmass,SRT,1.)
7219
7220
7221 ENDIF
7222 IF(FM.EQ.0.)FM=1.E-09
7223 NTRY1=0
7224 12 DM = RANART(NSEED) * (DMAX-DMIN) + DMIN
7225 NTRY1=NTRY1+1
7226 IF((RANART(NSEED) .GT. FD5(DM,SRT,1.)/FM).AND.
7227 1 (NTRY1.LE.10)) GOTO 12
7228
7229
7230
7231 if(dm.gt.1.84) goto 12
7232
7233 ENDIF
7234 13 CONTINUE
7235
7236 PRF=0.
7237 PF2=((SRT**2-DM**2+AVMASS**2)/(2.*SRT))**2-AVMASS**2
7238 IF(PF2.GT.0.)PRF=SQRT(PF2)
7239
7240
7241 IF(M12.EQ.37)THEN
7242 IF(iabs(LB(I1)).EQ.9)THEN
7243 LB(I1)=1
7244 E(I1)=AMP
7245 LB(I2)=11
7246 E(I2)=DM
7247 ELSE
7248 LB(I2)=1
7249 E(I2)=AMP
7250 LB(I1)=11
7251 E(I1)=DM
7252 ENDIF
7253 GO TO 207
7254 ENDIF
7255
7256 IF(M12.EQ.38)THEN
7257 IF(iabs(LB(I1)).EQ.9)THEN
7258 LB(I1)=1
7259 E(I1)=AMP
7260 LB(I2)=13
7261 E(I2)=DM
7262 ELSE
7263 LB(I2)=1
7264 E(I2)=AMP
7265 LB(I1)=13
7266 E(I1)=DM
7267 ENDIF
7268 GO TO 207
7269 ENDIF
7270
7271 IF(M12.EQ.39)THEN
7272 IF(iabs(LB(I1)).EQ.8)THEN
7273 LB(I1)=1
7274 E(I1)=AMP
7275 LB(I2)=11
7276 E(I2)=DM
7277 ELSE
7278 LB(I2)=1
7279 E(I2)=AMP
7280 LB(I1)=11
7281 E(I1)=DM
7282 ENDIF
7283 GO TO 207
7284 ENDIF
7285
7286 IF(M12.EQ.40)THEN
7287 IF(iabs(LB(I1)).EQ.8)THEN
7288 LB(I1)=1
7289 E(I1)=AMP
7290 LB(I2)=13
7291 E(I2)=DM
7292 ELSE
7293 LB(I2)=1
7294 E(I2)=AMP
7295 LB(I1)=13
7296 E(I1)=DM
7297 ENDIF
7298 GO TO 207
7299 ENDIF
7300
7301 IF(M12.EQ.41)THEN
7302 IF(iabs(LB(I1)).EQ.8)THEN
7303 LB(I1)=2
7304 E(I1)=AMN
7305 LB(I2)=11
7306 E(I2)=DM
7307 ELSE
7308 LB(I2)=2
7309 E(I2)=AMN
7310 LB(I1)=11
7311 E(I1)=DM