44 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
49 dimension eq(17),mij(17,17,4),nij(17,17,4),csjet(17,17,68),
50 *cs1(17,17,68),gz0(2),gz1(3)
51 COMMON /xsect/ gsect(10,5,4)
52 COMMON /area1/ ia(2),icz,icp
53 COMMON /area5/ rd(2),cr1(2),cr2(2),cr3(2)
55 COMMON /area6/
pi,bm,am
57 COMMON /area10/ stmass,am0,amn,amk,amc,amlamc,amlam,ameta
58 COMMON /area15/ fp(5),rq(5),cd(5)
60 COMMON /area17/ del,rs,rs0,fs,alfp,rr,sh,delh
61 COMMON /area18/ alm,qt0,qlog,qll,aqt0,qtf,bet,amj0
62 COMMON /area19/ ahl(5)
64 COMMON /area22/ sjv0,fjs0(5,3)
66 COMMON /area23/ rjv(50)
67 COMMON /area24/ rjs(50,5,10)
68 COMMON /area27/ fp0(5)
69 COMMON /area28/ arr(4)
70 COMMON /area29/ cstot(17,17,68)
71 COMMON /area30/ cs0(17,17,68)
72 COMMON /area31/ csborn(17,68)
73 COMMON /area32/ csq(17,2,2),csbq(17,2,2)
74 COMMON /area33/ fsud(10,2)
75 COMMON /area34/ qrt(10,101,2)
76 COMMON /area35/ sjv(10,5),fjs(10,5,15)
80 COMMON /area43/ moniou
83 COMMON /area44/ gz(10,5,4),gzp(10,5,4)
89 COMMON /area48/ asect(10,6,4)
91 COMMON /version/version
96 *
'====================================================',
98 * /,
' ',
'| QUARK GLUON STRING JET MODEL |',
100 * /,
' ',
'| HADRONIC INTERACTION MONTE CARLO |',
102 * /,
' ',
'| N.N. KALMYKOV AND S.S. OSTAPCHENKO |',
104 * /,
' ',
'| e-mail: serg@eas.npi.msu.su |',
106 * /,
' ',
'| Publication to be cited when using this program: |',
107 * /,
' ',
'| N.N. Kalmykov & S.S. Ostapchenko, A.I. Pavlov |',
108 * /,
' ',
'| Nucl. Phys. B (Proc. Suppl.) 52B (1997) 17 |',
110 * /,
' ',
'| last modification: June 15, 2005 by D. Heck |',
111 * /,
' ',
'| (version qgsjet01c.f) |',
112 * /,
' ',
'====================================================',
114 IF(debug.GE.1)
WRITE (moniou,210)
115 210
FORMAT(2
x,
'PSAINI - MAIN INITIALIZATION PROCEDURE')
121 ahl(1)=1.
d0-2.
d0*arr(1)
122 ahl(2)=1.
d0-arr(1)-arr(2)
123 ahl(3)=1.
d0-arr(1)-arr(3)
124 ahl(4)=1.
d0-arr(1)-arr(4)
125 ahl(5)=ahl(2)+arr(1)-arr(4)
130 cc(2)=1.
d0/dsqrt(cd(2))
131 cc(1)=1.
d0/cc(2)/cd(1)
132 cc(3)=1.
d0/cc(2)/cd(3)
133 cc(4)=1.
d0/cc(2)/cd(4)
134 cc(5)=1.
d0/cc(2)/cd(5)
151 INQUIRE(file=
'QGSDAT01',exist=lcalc)
153 IF(debug.GE.1)
WRITE (moniou,211)
154 211
FORMAT(2
x,
'PSAINI: HARD CROSS SECTION RATIOS READOUT FROM THE' 156 OPEN(1,file=
'QGSDAT01',status=
'OLD')
157 READ (1,*)csborn,cs0,cstot,csq,csbq,
158 * fsud,qrt,sjv,fjs,rjv,rjs,gz,gzp,gsect
163 IF(debug.GE.1)
WRITE (moniou,201)
164 201
FORMAT(2
x,
'PSAINI: HARD CROSS SECTIONS CALCULATION')
170 1 eq(
i)=qt0*4.
d0**float(
i-1)
180 k1=k+17*(
m-1)+34*(l-1)
181 IF(k.LE.
i.OR.k.EQ.2)
THEN 199 k1=k+17*(
m-1)+34*(l-1)
201 IF(k.LE.
n.OR.k.EQ.2)
THEN 207 cstot(
i,
j,k1)=dlog(csborn(
n,k1))
208 cs0(
i,
j,k1)=cstot(
i,
j,k1)
215 IF(debug.GE.2)
WRITE (moniou,202)
n,eq(mij(1,1,1)),eq(nij(1,1,1))
216 202
FORMAT(2
x,
'PSAINI: NUMBER OF LADDER RUNS TO BE CONSIDERED:',i2/
217 * 4
x,
'MINIMAL MASSES SQUARED FOR THE UNORDERED AND STRICTLY',
218 *
' ORDERED LADDERS:'/4
x,
e10.3,3
x,
e10.3)
228 s2min=max(qq,4.
d0*qt0)
231 smin=s2min*(1.
d0+sm)/(1.
d0-sm)
244 nij(
i,
j,ml)=nij(
i,
j,ml)+1
246 k1=k+17*(
m-1)+34*(l-1)
263 k1=k+17*(
m-1)+34*(l-1)
266 csj=cs1(
i,
j,k1)+csborn(max(
i,
j),k1)
267 IF(debug.GE.2)
WRITE (moniou,204)csj,exp(cs0(
i,
j,k1))
268 204
FORMAT(2
x,
'PSAINI: NEW AND OLD VALUES OF THE CONTRIBUTION',
269 *
' OF THE STRICTLY ORDERED LADDER:'/4
x,
e10.3,3
x,
e10.3)
270 IF(csj.EQ.0.
d0.OR.abs(1.
d0-exp(cs0(
i,
j,k1))/csj).LT.1.
d-2)
THEN 271 nij(
i,
j,ml)=nij(
i,
j,ml)+1
274 cs0(
i,
j,k1)=dlog(csj)
285 s2min=max(qq,4.
d0*qt0)
288 smin=s2min*(1.
d0+sm)/(1.
d0-sm)
300 mij(
i,
j,ml)=mij(
i,
j,ml)+1
302 k1=k+17*(
m-1)+34*(l-1)
306 cs1(
i,
j,k1)=
psjet(qi,qj,sk,s2min,
m-1,l)
321 k1=k+17*(
m-1)+34*(l-1)
322 k2=k+17*(l-1)+34*(
m-1)
323 csj=cs1(
i,
j,k1)+exp(cs0(
j,
i,k2))
324 IF(csj.EQ.0.
d0.OR.abs(1.
d0-exp(cstot(
i,
j,k1))/csj).LT.1.
d-2)
325 * mij(
i,
j,ml)=mij(
i,
j,ml)+1
326 IF(debug.GE.2)
WRITE (moniou,203)csj,exp(cstot(
i,
j,k1))
327 203
FORMAT(2
x,
'PSAINI: NEW AND OLD VALUES OF THE UNORDERED LADDER',
328 *
' CROSS SECTION:'/4
x,
e10.3,3
x,
e10.3)
329 11 cstot(
i,
j,k1)=dlog(csj)
336 13
IF(mij(1,1,l).LE.17.OR.nij(1,1,l).LE.17)
GOTO 4
344 k1=k+17*(
m-1)+34*(l-1)
345 IF(k.LE.
i.OR.k.EQ.2)
THEN 348 csborn(
i,k1)=dlog(csborn(
i,k1))
361 k1=k+17*(
m-1)+34*(l-1)
362 csbq(k,
m,l)=csborn(1,k1)
363 csq(k,
m,l)=cstot(1,1,k1)
375 qmax=qtf*4.
d0**(1.
d0+k)
384 qlmax=1.38629
d0*(k-1)
396 IF(debug.GE.2)
WRITE (moniou,205)
397 205
FORMAT(2
x,
'PSAINI: PRETABULATION OF THE INTERACTION EIKONALS')
426 fs=fp(icz)*exp(y0*del)/rs*cd(icz)
428 rp1=rs*4.
d0*.0391
d0/am**2
430 g0=
pi*rp1/cd(icz)*am**2*10.
d0 443 fjs0(
i,
m)=fjs(ie,
i,m1)
449 IF(debug.GE.1)
WRITE (moniou,206)e0n,ty(icz),ia(2)
450 206
FORMAT(2
x,
'PSAINI: INITIAL PARTICLE ENERGY:',
e10.3,2
x,
451 *
'ITS TYPE:',a7,2
x,
'TARGET MASS NUMBER:',i2)
456 rd(2)=0.7
d0*float(ia(2))**.446/am
459 rd(2)=.9
d0*float(ia(2))**.3333/am
471 if (debug .ge.1)
write (moniou,*)gz0
479 gdp=(1.
d0-cc(icz))*cc(2)*gd0
481 gdt=(1.
d0-cc(2))*cc(icz)*gd0
483 gdd=(1.
d0-cc(icz))*(1.
d0-cc(2))*gd0
486 gel=gd0*cc(icz)*cc(2)
488 IF(debug.GE.1)
WRITE (moniou,225)gtot,gin,gel,gdp,gdt,gdd
490 225
FORMAT(2
x,
'PSAINI: HADRON-PROTON CROSS SECTIONS:'/
492 *
'GDIFR_PROJ=',
e10.3,2
x,
'GDIFR_TARG=',
e10.3,2
x,
493 *
'G_DOUBLE_DIFR',
e10.3)
495 gz(ie,icz,iia)=gdt/gin
496 gzp(ie,icz,iia)=(gdp+gdd)/gin
498 gsect(ie,icz,iia)=log(gin)
510 anorm=1.5
d0/
pi/rrr**3/(1.
d0+(
pi/rrr)**2)*rp1
519 gin=gz1(1)+gz1(2)+gz1(3)
521 IF(debug.GE.1)
WRITE (moniou,224)
522 * gin*10.
d0,gz1(1)*10.
d0,gz1(2)*10.
d0 524 224
FORMAT(2
x,
'PSAINI: HADRON-NUCLEUS CROSS SECTIONS:'/
525 * 4
x,
'GIN=',
e10.3,2
x,
'GDIFR_TARG=',
e10.3,2
x,
526 *
'GDIFR_PROJ=',
e10.3)
528 gz(ie,icz,iia)=gz1(1)/gin
529 gzp(ie,icz,iia)=gz1(2)/gin
532 gsect(ie,icz,iia)=log(gin)
558 IF(debug.GE.1)
WRITE (moniou,212)
559 212
FORMAT(2
x,
'PSAINI: HARD CROSS SECTIONS ARE WRITTEN TO THE FILE' 561 OPEN(1,file=
'QGSDAT01',status=
'unknown')
562 WRITE (1,*)csborn,cs0,cstot,csq,csbq,
563 * fsud,qrt,sjv,fjs,rjv,rjs,gz,gzp,gsect
570 INQUIRE(file=
'SECTNU',exist=lcalc)
572 IF(debug.GE.1)
WRITE (moniou,208)
573 208
FORMAT(2
x,
'PSAINI: NUCLEAR CROSS SECTIONS READOUT FROM THE FILE' 575 OPEN(2,file=
'SECTNU',status=
'OLD')
587 IF(debug.GE.1)
WRITE (moniou,207)e0n,iap,iat
588 207
FORMAT(2
x,
'PSAINI: INITIAL NUCLEUS ENERGY:',
e10.3,2
x,
589 *
'PROJECTILE MASS:',i2,2
x,
'TARGET MASS:',i2)
590 CALL xxaini(e0n,2,iap,iat)
591 CALL crossc_kk(niter,gtot,gprod,gabs,gdd,gqel,gcoh)
592 IF(debug.GE.1)
WRITE (moniou,209)
593 * gtot,gprod,gabs,gdd,gqel,gcoh
595 209
FORMAT(2
x,
'GTOT',d10.3,
' GPROD',d10.3,
' GABS',d10.3/2
x,
596 *
'GDD',d10.3,
' GQEL',d10.3,
' GCOH',d10.3)
597 asect(ie,iia1,iia2)=log(gprod)
601 OPEN(2,file=
'SECTNU',status=
'UNKNOWN')
606 IF(debug.GE.3)
WRITE (moniou,218)
607 218
FORMAT(2
x,
'PSAINI - END')
618 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
620 COMMON /area43/ moniou
624 IF(debug.GE.2)
WRITE (moniou,201)
x,
j,l
625 201
FORMAT(2
x,
'PSAPINT: X=',
e10.3,2
x,
'J= ',i1,2
x,
'L= ',i1)
640 IF(debug.GE.2)
WRITE (moniou,202)
psapint 641 202
FORMAT(2
x,
'PSAPINT=',
e10.3)
649 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
652 COMMON /area15/ fp(5),rq(5),cd(5)
653 COMMON /area17/ del,rs,rs0,fs,alfp,rr,sh,delh
654 COMMON /area18/ alm,qt0,qlog,qll,aqt0,qtf,bet,amj0
655 COMMON /area25/ ahv(5)
656 COMMON /area26/ factork
657 COMMON /area41/ ty(5)
658 COMMON /area43/ moniou
662 IF(debug.GE.1)
WRITE (moniou,210)
663 210
FORMAT(2
x,
'PSASETC - COMMON MODEL PARAMETERS SETTING')
744 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
746 dimension wi(3),wk(3)
747 COMMON /area18/ alm,qt0,qlog,qll,aqt0,qtf,bet,amj0
748 COMMON /area31/ csj(17,68)
749 COMMON /area43/ moniou
753 IF(debug.GE.2)
WRITE (moniou,201)qq,s,
m,l
754 201
FORMAT(2
x,
'PSBINT: QQ=',
e10.3,2
x,
'S= ',
e10.3,2
x,
'M= ',i1,2
x,
757 IF(s.LE.max(4.
d0*qt0,qq))
THEN 758 IF(debug.GE.3)
WRITE (moniou,202)
psbint 759 202
FORMAT(2
x,
'PSBINT=',
e10.3)
764 qli=dlog(qq/qt0)/1.38629
d0 765 sl=dlog(s/qt0)/1.38629
d0 774 wi(3)=wi(2)*(wi(2)-1.
d0)*.5
d0 775 wi(1)=1.
d0-wi(2)+wi(3)
776 wi(2)=wi(2)-2.
d0*wi(3)
778 wk(3)=wk(2)*(wk(2)-1.
d0)*.5
d0 779 wk(1)=1.
d0-wk(2)+wk(3)
780 wk(2)=wk(2)-2.
d0*wk(3)
786 ELSEIF(sql.LT.1.
d0.AND.
i.NE.0)
THEN 789 wi(3)=wi(2)*(wi(2)-1.
d0)*.5
d0 790 wi(1)=1.
d0-wi(2)+wi(3)
791 wi(2)=wi(2)-2.
d0*wi(3)
808 IF(
i+kl.GT.12)
i=12-kl
811 wi(3)=wi(2)*(wi(2)-1.
d0)*.5
d0 812 wi(1)=1.
d0-wi(2)+wi(3)
813 wi(2)=wi(2)-2.
d0*wi(3)
815 wk(3)=wk(2)*(wk(2)-1.
d0)*.5
d0 816 wk(1)=1.
d0-wk(2)+wk(3)
817 wk(2)=wk(2)-2.
d0*wk(3)
830 wi(3)=wi(2)*(wi(2)-1.
d0)*.5
d0 831 wi(1)=1.
d0-wi(2)+wi(3)
832 wi(2)=wi(2)-2.
d0*wi(3)
841 IF(debug.GE.3)
WRITE (moniou,202)
psbint 846 FUNCTION psborn(QQ,S,IQ1,IQ2)
852 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
854 COMMON /area6/
pi,bm,am
855 COMMON /area18/ alm,qt0,qlog,qll,aqt0,qtf,bet,amj0
856 COMMON /area26/ factork
857 COMMON /area43/ moniou
859 COMMON /ar3/
x1(7),a1(7)
862 IF(debug.GE.2)
WRITE (moniou,201)qq,s,iq1,iq2
863 201
FORMAT(2
x,
'PSBORN: QQ=',
e10.3,2
x,
'S= ',
e10.3,2
x,
'IQ1= ',i1,2
x,
865 tmin=s*(.5
d0-dsqrt(.25
d0-qt0/s))
866 tmin=max(tmin,s*qq/(s+qq))
883 IF(debug.GE.3)
WRITE (moniou,202)
psborn 884 202
FORMAT(2
x,
'PSBORN=',
e10.3)
889 SUBROUTINE pscajet(QQ,IQ1,QV,ZV,QM,IQV,LDAU,LPAR,JQ)
908 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
910 dimension qmax(30,50),iqm(2),lnv(50),
911 * qv(30,50),zv(30,50),qm(30,50),iqv(30,50),
912 * ldau(30,49),lpar(30,50)
915 COMMON /area18/ alm,qt0,qlog,qll,aqt0,qtf,bet,amj0
916 COMMON /area43/ moniou
922 IF(debug.GE.2)
WRITE (moniou,201)qq,iq1,jq
923 201
FORMAT(2
x,
'PSCAJET: QQ=',
e10.3,2
x,
'IQ1= ',i1,2
x,
'JQ=',i1)
933 2 qlmax=dlog(qmax(nrow,nlev)/qtf/16.
d0)
934 iq=min(1,iabs(iqv(nrow,nlev)))+1
943 lpar(ll+1,nlev+1)=nrow
947 IF((3-2*jq)*iqv(nrow,nlev).GT.0)
THEN 949 iqm(2)=iqv(nrow,nlev)
952 iqm(1)=iqv(nrow,nlev)
960 IF(
qsran(b10).LT.wg)
THEN 964 iqm(1)=int(3.
d0*
qsran(b10)+1.
d0)*(3-2*jq)
975 qmax(nrow,nlev)=q*
z**2
976 qmax(nrow+1,nlev)=q*(1.
d0-
z)**2
977 iqv(nrow,nlev)=iqm(1)
978 iqv(nrow+1,nlev)=iqm(2)
979 IF(debug.GE.3)
WRITE (moniou,203)nlev,nrow,q,
z 980 203
FORMAT(2
x,
'PSCAJET: NEW BRANCHING AT LEVEL NLEV=',i2,
981 *
' NROW=',i2/4
x,
' EFFECTIVE MOMENTUM Q=',
e10.3,2
x,
' Z=',
e10.3)
988 IF(debug.GE.3)
WRITE (moniou,204)nlev,nrow
989 204
FORMAT(2
x,
'PSCAJET: NEW FINAL JET AT LEVEL NLEV=',i2,
995 IF(debug.GE.3)
WRITE (moniou,202)
996 202
FORMAT(2
x,
'PSCAJET - END')
999 lprow=lpar(nrow,nlev)
1001 IF(ldau(lprow,nlev-1).EQ.nrow)
THEN 1006 qm(lprow,nlev-1)=
z*(1.
d0-
z)*qv(lprow,nlev-1)
1007 * +qm(nrow-1,nlev)/
z+qm(nrow,nlev)/(1.
d0-
z)
1010 IF(debug.GE.3)
WRITE (moniou,205)nlev,nrow,qm(lprow,nlev)
1011 205
FORMAT(2
x,
'PSCAJET: JET MASS AT LEVEL NLEV=',i2,
1012 *
' NROW=',i2,
' - QM=',
e10.3)
1022 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
1029 dimension xa(64,3),xb(64,3),fhard(3)
1030 COMMON /area1/ ia(2),icz,icp
1031 COMMON /area2/ s,y0,wp0,wm0
1032 COMMON /area6/
pi,bm,am
1046 COMMON /area9/ lqa(56),lqb(56),nqs(1000),ias(1000),ibs(1000),
1047 * lha(56),lhb(56),zh(4000),iah(4000),ibh(4000),
1048 * iqh(4000),lva(56),lvb(56)
1049 COMMON /area10/ stmass,am0,amn,amk,amc,amlamc,amlam,ameta
1053 COMMON /area16/ cc(5)
1054 COMMON /area40/ jdifr
1055 COMMON /area43/ moniou
1057 COMMON /area45/ gdt,gdp
1060 COMMON /debug/ debug
1065 integer ng1evt,ng2evt,ikoevt
1066 real rglevt,sglevt,eglevt,fglevt,typevt
1067 common/c2evt/ng1evt,ng2evt,rglevt,sglevt,eglevt,fglevt,ikoevt
1070 IF(debug.GE.1)
WRITE (moniou,201)
1071 201
FORMAT(2
x,
'PSCONF - CONFIGURATION OF THE INTERACTION')
1077 IF(jdifr.EQ.1.AND.
qsran(b10).LT.gdt)
THEN 1081 ict=int(2.5+
qsran(b10))
1089 CALL xxdtg(wpi,wmi,icp,ict,0)
1092 ELSEIF(abs(jdifr).EQ.1.AND.
qsran(b10).LT.gdp)
THEN 1095 ict=int(2.5+
qsran(b10))
1100 IF(debug.GE.2)
WRITE (moniou,206)
1101 206
FORMAT(2
x,
'PROJECTILE HADRON DIFFRACTION')
1106 CALL xxdpr(wpi,wmi,icp0,ict,lq)
1158 6 nt=nt+int(cc(2)+
qsran(b10))
1161 IF(debug.GE.3)
WRITE (moniou,203)nt
1162 203
FORMAT(2
x,
'PSCONF: NUMBER OF ACTIVE TARGET NUCLEONS NT=',
1166 CALL psgea(ia(2),xb,2)
1180 IF(debug.GE.2)
WRITE (moniou,204)
b*am
1181 204
FORMAT(2
x,
'PSCONF: IMPACT PARAMETER FOR THE INTERACTION:',
1187 IF(ia(1).GT.1)
CALL psgea(ia(1),xa,1)
1200 IF(debug.GE.2.AND.icz.EQ.2)
WRITE (moniou,205)in
1201 205
FORMAT(2
x,
'PSCONF: ',i2,
'-TH PROJECTILE NUCLEON')
1204 IF(ia(1).NE.1.AND.
qsran(b10).GT.cc(2))
GOTO 12
1218 vv=2.
d0*
psfaz(
z,fsoft,fhard,fshard)
1221 eh=fhard(1)+fhard(2)+fhard(3)
1226 IF(aks.GT.1.
d0-ev*(1.
d0-2.
d0*
eh))
GOTO 11
1227 IF(debug.GE.2)
WRITE (moniou,208)
m 1228 208
FORMAT(2
x,
'PSCONF: INTERACTION WITH',i2,
'-TH TARGET NUCLEON')
1243 IF(aks1.LT.fhard(1))
THEN 1245 IF(lva(nw).NE.0)
GOTO 11
1249 ELSEIF(aks1.LT.fhard(1)+fhard(2))
THEN 1251 IF(lvb(
m).NE.0)
GOTO 11
1257 IF(lva(nw)+lvb(
m).NE.0)
GOTO 11
1280 7
IF(aks.LT.sum)
GOTO 8
1285 9 lnh=lnh+int(wh+
qsran(b10))
1292 IF(aks1.LT.fhard(1))
THEN 1293 IF(lva(nw).NE.0)
GOTO 22
1297 ELSEIF(aks1.LT.fhard(1)+fhard(2))
THEN 1298 IF(lvb(
m).NE.0)
GOTO 22
1303 IF(lva(nw)+lvb(
m).NE.0)
GOTO 22
1331 If (i1 .ge. 4000)
then 1332 write(moniou,*)
'psconf: I1 > 4000, index out of bounds' 1335 IF(
i.EQ.1.AND.iqv.NE.0)
THEN 1340 IF(debug.GE.2)
WRITE (moniou,209)i1,nw,
m,iqh(i1)
1341 209
FORMAT(2
x,
'PSCONF: ',i4,
'-TH HARD BLOCK IS CONNECTED TO',1
x,
1342 * i2,
'-TH PROJECTILE NUCLEON (HADRON) AND'/4
x,i2,
1343 *
'-TH TARGET NUCLEON; TYPE OF THE SEMIHARD INTERACTION:',i1)
1361 IF(debug.GE.2)
WRITE (moniou,210)ls,nw,
m,
n 1362 210
FORMAT(2
x,
'PSCONF: ',i4,
'-TH SOFT BLOCK IS CONNECTED TO',1
x,
1363 * i2,
'-TH PROJECTILE NUCLEON (HADRON) AND'/4
x,i2,
1364 *
'-TH TARGET NUCLEON; NUMBER OF POMERONS IN THE BLOCK NP=',
1379 IF(jdifr.EQ.1 .AND. ia(1).NE. 1
1380 * .AND.
qsran(b10).LT.(1.
d0-cc(icz))*
psv(
x,
y,xb,nt))
THEN 1384 ict=int(2.5+
qsran(b10))
1390 IF(debug.GE.2)
WRITE(moniou,207)in
1391 207
FORMAT(2
x,i2,
'-TH PROJECTILE NUCLEON DIFFRACTION')
1392 icp0=int(2.5+
qsran(b10))
1400 CALL xxdpr(wpi,wmi,icp0,ict,lq)
1414 13 xa(ns,l)=xa(in,l)
1422 IF(debug.GE.3)
WRITE (moniou,211)
1423 211
FORMAT(2
x,
'PSCONF: NO ONE NUCLEON (HADRON) INTERACTS - ',
1430 WRITE (moniou,213)nhp
1431 213
FORMAT(2
x,
'PSCONF: TOO GREAT NUMBER OF HARD POMERONS: NHP=',
1432 * i5,
' - REJECTION')
1448 20
IF(nw.NE.0)
CALL psshar(ls,nhp,nw,nt)
1450 IF(debug.GE.3)
WRITE (moniou,212)
1451 212
FORMAT(2
x,
'PSCONF - END')
1454 entry qgs1getdiffcode( jcode )
1460 SUBROUTINE pscs(C,S)
1463 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
1466 COMMON /area43/ moniou
1467 COMMON /debug/ debug
1471 IF(debug.GE.2)
WRITE (moniou,201)
1472 201
FORMAT(2
x,
'PSCS - COS(FI) AND SIN(FI) ARE GENERATED',
1477 IF(s3.GT.1.
d0)
GOTO 1
1481 IF(debug.GE.3)
WRITE (moniou,202)
c,s
1482 202
FORMAT(2
x,
'PSCS: C=',
e10.3,2
x,
'S=',
e10.3)
1491 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
1493 dimension ey(3),ep(4)
1494 COMMON /area43/ moniou
1495 COMMON /debug/ debug
1498 IF(debug.GE.2)
WRITE (moniou,201)ep,s
1499 201
FORMAT(2
x,
'PSDEFTR - LORENTZ BOOST PARAMETERS:'/
1500 * 4
x,
'4-VECTOR EP=',4
e10.3/4
x,
'4-VECTOR SQUARED S=',
e10.3)
1502 IF(ep(
i+1).EQ.0.
d0)
THEN 1507 IF(wm/wp.LT.1.
d-8)
THEN 1510 1
IF(l.NE.
i)ww=ww+ep(l+1)**2
1518 IF(debug.GE.3)
WRITE (moniou,202)ey
1519 202
FORMAT(2
x,
'PSDEFTR: LORENTZ BOOST PARAMETERS EY(I)=',2
x,3
e10.3)
1524 SUBROUTINE psdefrot(EP,S0X,C0X,S0,C0)
1528 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
1531 COMMON /area43/ moniou
1532 COMMON /debug/ debug
1535 IF(debug.GE.2)
WRITE (moniou,201)ep
1536 201
FORMAT(2
x,
'PSDEFROT - SPACIAL ROTATION PARAMETERS'/4
x,
1537 *
'4-VECTOR EP=',2
x,4(
e10.3,1
x))
1539 pt2=ep(3)**2+ep(4)**2
1547 pl=dsqrt(pt2+ep(2)**2)
1561 IF(debug.GE.3)
WRITE (moniou,202)s0x,c0x,s0,c0,ep
1562 202
FORMAT(2
x,
'PSDEFROT: SPACIAL ROTATION PARAMETERS'/
1563 * 4
x,
'S0X=',
e10.3,2
x,
'C0X=',
e10.3,2
x,
'S0=',
e10.3,2
x,
'C0=',
e10.3/
1564 * 4
x,
'ROTATED 4-VECTOR EP=',4(
e10.3,1
x))
1572 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
1575 COMMON /area43/ moniou
1576 COMMON /debug/ debug
1579 IF(debug.GE.2)
WRITE (moniou,201)
x,
y 1580 201
FORMAT(2
x,
'PSDR: NUCLEON COORDINATES - X=',
e10.3,2
x,
'Y=',
e10.3)
1582 IF(debug.GE.3)
WRITE (moniou,202)
psdr 1583 202
FORMAT(2
x,
'PSDR=',
e10.3)
1588 FUNCTION psfap(X,J,L)
1594 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
1596 COMMON /area43/ moniou
1597 COMMON /debug/ debug
1600 IF(debug.GE.2)
WRITE (moniou,201)
x,
j,l
1601 201
FORMAT(2
x,
'PSFAP - ALTARELLI-PARISI FUNCTION:',2
x,
1602 *
'X=',
e10.3,2
x,
'J=',i1,2
x,
'L=',i1)
1616 IF(debug.GE.3)
WRITE (moniou,202)
psfap 1617 202
FORMAT(2
x,
'PSFAP=',
e10.3)
1622 FUNCTION psfaz(Z,FSOFT,FHARD,FSHARD)
1630 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
1633 COMMON /area17/ del,rs,rs0,fs,alf,rr,sh,delh
1634 COMMON /area22/ sjv,fjs(5,3)
1635 COMMON /area43/ moniou
1636 COMMON /debug/ debug
1639 IF(debug.GE.2)
WRITE (moniou,201)
z 1640 201
FORMAT(2
x,
'PSFAZ - HADRON-NUCLEON (NUCLEON-NUCLEON)',
1641 *
' INTERACTION EIKONAL; Z=',
e10.3)
1643 fhard(3)=sjv*
z**(rs/rs0)
1651 fsr=(exp(fjs(1,
i))*wz+(exp(fjs(2,
i))-2.
d0*
1652 * exp(fjs(1,
i)))*wz*(wz-1.
d0)*.5
d0)*
z 1654 fsr=exp(fjs(jz,
i)+(fjs(jz+1,
i)-fjs(jz,
i))*wz
1655 * +(fjs(jz+2,
i)+fjs(jz,
i)-2.
d0*fjs(jz+1,
i))
1656 * *wz*(wz-1.
d0)*.5
d0)*
z 1666 IF(debug.GE.3)
WRITE (moniou,202)
psfaz,fsoft,fshard,fhard
1667 202
FORMAT(2
x,
'PSFAZ=',
e10.3,2
x,
'FSOFT=',
e10.3,2
x,
'FSHARD=',
e10.3/4
x 1680 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
1682 COMMON /area18/ alm,qt0,qlog,qll,aqt0,qtf,bet,amj0
1683 COMMON /area43/ moniou
1684 COMMON /debug/ debug
1687 IF(debug.GE.2)
WRITE (moniou,201)s,
t,iq1,iq2
1688 201
FORMAT(2
x,
'PSFBORN - HARD SCATTERING MATRIX ELEMENT SQUARED:'/
1689 * 4
x,
'S=',
e10.3,2
x,
'|T|=',
e10.3,2
x,
'IQ1=',i2,2
x,
'IQ2=',i2)
1691 IF(iq1.EQ.0.AND.iq2.EQ.0)
THEN 1694 ELSEIF(iq1*iq2.EQ.0)
THEN 1697 ELSEIF(iq1.EQ.iq2)
THEN 1699 psfborn=((s**2+u**2)/
t**2+(s**2+
t**2)/u**2)/2.25
d0 1701 ELSEIF(iq1+iq2.EQ.0)
THEN 1703 psfborn=((s**2+u**2)/
t**2+(u**2+
t**2)/s**2)/2.25
d0 1709 IF(debug.GE.2)
WRITE (moniou,202)
psfborn 1710 202
FORMAT(2
x,
'PSFBORN=',
e10.3)
1715 FUNCTION psfsh(S,Z,ICZ,IQQ)
1722 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
1724 COMMON /area6/
pi,bm,am
1725 COMMON /area15/ fp(5),rq(5),cd(5)
1726 COMMON /area17/ del,rs,rs0,fs,alf,rr,sh,delh
1727 COMMON /area18/ alm,qt0,qlog,qll,aqt0,qtf,bet,amj0
1728 COMMON /area19/ ahl(5)
1729 COMMON /area25/ ahv(5)
1730 COMMON /area27/ fp0(5)
1731 COMMON /ar3/
x1(7),a1(7)
1732 COMMON /area43/ moniou
1733 COMMON /debug/ debug
1736 IF(debug.GE.2)
WRITE (moniou,201)s,
z,iqq,icz
1737 201
FORMAT(2
x,
'PSFSH - SEMIHARD INTERACTION EIKONAL:'/
1738 * 4
x,
'S=',
e10.3,2
x,
'Z=',
e10.3,2
x,
'IQQ=',i1,2
x,
'ICZ=',i1)
1740 xmin=xmin**(delh-del)
1745 ELSEIF(iqq.EQ.2)
THEN 1754 z1=(.5
d0*(1.
d0+xmin-(2*
m-3)*
x1(
i)*(1.
d0-xmin)))**(1.
d0/
1760 CALL psjint0(z1*s,sj,sjb,iq,0)
1765 IF(debug.GE.3)
WRITE (moniou,203)z1*s,gy
1766 203
FORMAT(2
x,
'PSFSH:',2
x,
'S_HARD=',
e10.3,2
x,
'SIGMA_HARD=',
e10.3)
1773 1 st2=st2+a1(
j)*
psftild(z1**xx,icz)*
1777 psfsh=
psfsh-a1(
i)*dlog(z1)*gy/z1**delh*
z**(rs/rh)/rh*st2
1785 xa=(xam+(1.
d0-xam)*xx)**(1.
d0/(del+.5
d0))
1786 rh=rs0+alf*dlog(xa/z1)
1787 2 st2=st2+a1(
j)*(1.
d0-xa)**ahv(icv)*
z**(rs/rh)/rh*
1802 IF(icz.EQ.2.OR.iqq.EQ.2)
THEN 1804 ELSEIF((icz-1)*(icz-3)*(icz-5).EQ.0)
THEN 1808 IF(debug.GE.3)
WRITE (moniou,202)
psfsh 1809 202
FORMAT(2
x,
'PSFSH=',
e10.3)
1820 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
1822 COMMON /area17/ del,rs,rs0,fs,alfp,rr,sh,delh
1823 COMMON /area18/ alm,qt0,qlog,qll,aqt0,qtf,bet,amj0
1824 COMMON /area19/ ahl(5)
1825 COMMON /ar3/
x1(7),a1(7)
1826 COMMON /area43/ moniou
1827 COMMON /debug/ debug
1830 IF(debug.GE.2)
WRITE (moniou,201)
z,icz
1831 201
FORMAT(2
x,
'PSFTILD:',2
x,
'Z=',
e10.3,2
x,
'ICZ=',i1)
1839 IF(debug.GE.3)
WRITE (moniou,202)
psftild 1840 202
FORMAT(2
x,
'PSFTILD=',
e10.3)
1845 SUBROUTINE psgea(IA,XA,JJ)
1849 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
1853 COMMON /area5/ rd(2),ca1(2),ca2(2),ca3(2)
1855 COMMON /area43/ moniou
1856 COMMON /debug/ debug
1860 IF(debug.GE.2)
WRITE (moniou,201)jj,ia
1861 201
FORMAT(2
x,
'PSGEA - CONFIGURATION OF THE NUCLEUS ',i1,
';',2
x,
1862 *
'COORDINATES FOR ',i2,
' NUCLEONS')
1871 3
IF(zuk.GT.ca2(jj))
GOTO 4
1872 tt=-dlog(
qsran(b10))
1874 4
IF(zuk.GT.ca3(jj))
GOTO 5
1878 6
IF(
qsran(b10).GT.1.
d0/(1.
d0+exp(-abs(tt))))
GOTO 1
1881 rim=dsqrt(rim*rim-
z*
z)
1894 xa(k,l)=summ-aks*sqrt(float(
j)/k)
1895 8 summ=summ+aks/sqrt(float(
j*k))
1901 206
WRITE (moniou,204)
i,(xa(
i,l),l=1,3)
1904 202
FORMAT(2
x,
'PSGEA - END')
1905 203
FORMAT(2
x,
'PSGEA: POSITIONS OF THE NUCLEONS')
1906 204
FORMAT(2
x,
'PSGEA: ',i2,
' - ',3(
e10.3,1
x))
1915 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
1917 COMMON /ar3/
x1(7),a1(7)
1918 COMMON /area43/ moniou
1919 COMMON /debug/ debug
1922 f(
z,
x)=(1.-exp(-.5*
z*(1.+
x)))/(1.+
x)
1924 IF(debug.GE.2)
WRITE (moniou,201)
z 1925 201
FORMAT(2
x,
'PSGINT:',2
x,
'Z=',
e10.3)
1929 IF(debug.GE.3)
WRITE (moniou,202)
psgint 1930 202
FORMAT(2
x,
'PSGINT=',
e10.3)
1940 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
1942 COMMON /ar3/
x1(7),a1(7)
1943 COMMON /area6/
pi,bm,am
1944 COMMON /area15/ fp(5),rq(5),cd(5)
1945 COMMON /area17/ del,rs,rs0,fs,alf,rr,sh,delh
1946 COMMON /area18/ alm,qt0,qlog,qll,aqt0,qtf,bet,amj0
1947 COMMON /area19/ ahl(5)
1948 COMMON /area25/ ahv(5)
1949 COMMON /area43/ moniou
1950 COMMON /debug/ debug
1953 IF(debug.GE.2)
WRITE (moniou,201)s,icz
1954 201
FORMAT(2
x,
'PSHARD - HARD QUARK-QUARK INTERACTION CROSS',
1956 * 2
x,
'S=',
e10.3,2
x,
'ICZ=',i1)
1958 xmin=xmin**(delh+.5
d0)
1964 z1=(.5
d0*(1.
d0+xmin-(2*
m-3)*
x1(
i)*(1.
d0-xmin)))**(1.
d0/
1971 st2=st2+a1(
j)*(1.
d0-z1**xx)**ahv(icz)*
1972 * (1.
d0-z1**(1.
d0-xx))**ahv(2)
1985 IF(debug.GE.3)
WRITE (moniou,203)z1*s,gy
1986 203
FORMAT(2
x,
'PSHARD:',2
x,
'S_HARD=',
e10.3,2
x,
'SIGMA_HARD=',
e10.3)
1997 ELSEIF((icz-1)*(icz-3)*(icz-5).EQ.0)
THEN 2007 IF(debug.GE.2)
WRITE (moniou,202)
pshard 2008 202
FORMAT(2
x,
'PSHARD=',
e10.3)
2013 SUBROUTINE pshot(WP0,WM0,Z,IPC,EPC,IZP,IZT,ICZ,IQQ)
2021 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
2024 dimension ep(4,2),ept(4),ept0(4),ep3(4),epj(4),epj1(4),ey(3),
2025 * qmin(2),wp(2),iqc(2),iqp(2),
2026 * ipc(2,2),epc(8,2),iqj(2),eqj(4,2),ipq(2,2),epq(8,2),
2028 * qv1(30,50),zv1(30,50),qm1(30,50),iqv1(2,30,50),
2029 * ldau1(30,49),lpar1(30,50),
2030 * qv2(30,50),zv2(30,50),qm2(30,50),iqv2(2,30,50),
2031 * ldau2(30,49),lpar2(30,50)
2032 COMMON /area6/
pi,bm,ammm
2033 COMMON /area8/ wwm,be(4),dc(5),deta,almpt
2034 COMMON /area10/ stmass,am(7)
2036 COMMON /area17/ del,rs,rs0,fs,alf,rr,sh,delh
2037 COMMON /area18/ alm,qt0,qlog,qll,aqt0,qtf,bet,amj0
2038 COMMON /area42/ tyq(15)
2039 COMMON /area43/ moniou
2040 COMMON /area46/ epjet(4,2,15000),ipjet(2,15000)
2041 COMMON /area47/ njtot
2042 COMMON /debug/ debug
2046 IF(debug.GE.1)
WRITE (moniou,201)iqq,wp0,wm0
2047 201
FORMAT(2
x,
'PSHOT - SEMIHARD INTERACTION SIMULATION:'/
2048 * 4
x,
'TYPE OF THE INTERACTION:',i2/
2049 * 4
x,
'INITIAL LIGHT CONE MOMENTA:',2
e10.3)
2081 xmin=4.
d0*(qt0+amj0)/s
2082 xmin1=xmin**(delh-del)
2089 gb0=-dlog(xmin)*(1.
d0-dsqrt(xmin))**(2.
d0*bet)
2091 gb0=(1.
d0-xmin)**bet
2103 gb0=gb0*gy0/s**delh/rs0*
z 2117 1 zpm=(xmin1+
qsran(b10)*(1.
d0-xmin1))**(1.
d0/(delh-del))
2122 CALL psjint0(zpm*s,sj,sjb,iq,0)
2125 rh=rs0-alf*dlog(zpm)
2133 gbyj=-dlog(zpm)*((1.-xp)*(1.-xm))**bet
2145 gbyj=(1.
d0-zpm)**bet
2156 gbyj=gbyj*gy*exp(-delh*yj)/gb0*
z**(rs/rh)/rh
2157 IF(
qsran(b10).GT.gbyj)
GOTO 1
2162 IF(debug.GE.2)
WRITE (moniou,203)s
2163 203
FORMAT(2
x,
'PSHOT: MASS SQUARED FOR THE HARD PARTON-PARTON',
2164 *
' INTERACTION:',
e10.3)
2174 IF((iqq-1)*(iqq-3).EQ.0)
THEN 2192 IF((iqq-2)*(iqq-3).EQ.0)
THEN 2210 ept(1)=.5
d0*(wpi+wmi)
2211 ept(2)=.5
d0*(wpi-wmi)
2221 303 epq(
i+4*(l-1),
m)=epc(
i+4*(l-1),
m)
2223 qminn=max(qmin(1),qmin(2))
2229 IF(debug.GE.2)
WRITE (moniou,208)ilad, si,iqc,ept
2230 208
FORMAT(2
x,
'PSHOT: ',i2,
'-TH HARD LADDER;',
2231 *
' MASS SQUARED FOR THE LADDDER:',
e10.3/
2232 * 4
x,
'LADDER END FLAVORS:',2i3/4
x,
2233 *
'LADDER 4-MOMENTUM: ',4
e10.3)
2235 ebal(1)=.5*(wp0+wm0)-ept(1)
2236 ebal(2)=.5*(wp0-wm0)-ept(2)
2241 ebal(l)=ebal(l)-epq(l,
m)
2242 501
if(iqc(
m).eq.0) ebal(l)=ebal(l)-epq(l+4,
m)
2246 502 ebal(l)=ebal(l)-epjet(l,
m,
i)
2251 pt2=ept(3)**2+ept(4)**2
2256 iqp(1)=min(1,iabs(iqc(1)))
2257 iqp(2)=min(1,iabs(iqc(2)))
2263 s2min=max(qminn,4.
d0*(qt0+amj0))
2266 wwmin=(s2min+(
pt-dsqrt(qt0))**2+(qt0+amj0)*(dsqrt(s2min/qt0)-
2267 * 1.
d0))/(1.
d0-dsqrt(qt0/s2min))
2271 sj=
psjint(qmin(1),qmin(2),si,iqp(1)+1,iqp(2)+1)
2272 sjb=
psbint(qminn,si,iqp(1)+1,iqp(2)+1)
2274 IF(debug.GE.2)
WRITE (moniou,251)s2min,wwmin,sj,sjb
2275 251
FORMAT(2
x,
'PSHOT: KINEMATICAL BOUNDS S2MIN=',
e10.3,
2276 * 2
x,
'WWMIN=',
e10.3/4
x,
'JET CROSS SETION SJ=',
e10.3,
2277 * 2
x,
'BORN CROSS SECTION SJB=',
e10.3)
2279 IF(
qsran(b10).LT.sjb/sj.
2280 * or.ww.LT.1.2
d0*wwmin)
GOTO 12
2282 IF((sj-sjb)/sj.GT..1
d0)
THEN 2283 sj1=
psjint1(qmin(1),qmin(2),si,iqp(1)+1,iqp(2)+1)
2284 sj2=
psjint1(qmin(2),qmin(1),si,iqp(2)+1,iqp(1)+1)
2285 dsj=(sj2-sj1)/(sj-sjb)*.5
d0 2293 aq=-(si+amj0+2.
d0*
pt*dsqrt(qt0))/ww
2297 qq=aq**3/13.5
d0-aq*bq/3.
d0+cq
2300 fq=atan(1.
d0/cosq**2-1.
d0)
2301 IF(cosq.LT.0.
d0)fq=
pi-fq
2307 xmax=1.
d0+aq/3.
d0-pq*(dsqrt(3.
d0)*sin(fq)-
cos(fq))
2311 qqmax=qt0/(1.
d0-xmax)**2
2312 qqmin=qt0/(1.
d0-xmin)**2
2314 IF(qqmin.LT.s2min)
THEN 2315 xmm=(si-s2min+amj0+2.
d0*
pt*dsqrt(qt0))/ww*.5
d0 2316 xmin=1.
d0-xmm-dsqrt(xmm*xmm-(qt0+amj0)/ww)
2317 qqmin=qt0/(1.
d0-xmin)**2
2319 IF(qqmin.LT.qmin(jj))
THEN 2321 xmm1=ww-2.
d0*
pt*dsqrt(qqmin)+qqmin
2322 xmm=(si-s2min+amj0)/xmm1*.5
d0 2323 xmin=1.
d0-xmm-dsqrt(xmm*xmm-amj0/xmm1)
2328 xm0=max(.5
d0,1.
d0-dsqrt(qt0/qmin(jj)))
2329 IF(xm0.GT..95
d0*xmax.OR.xm0.LT.1.05
d0*xmin)
2330 * xm0=.5
d0*(xmax+xmin)
2331 qm0=qt0/(1.
d0-xm0)**2
2334 sj0=
psjint(qm0,qmin(3-jj),s2max,1,iqp(3-jj)+1)*
2335 *
psfap(xm0,iqp(jj),0)+
2336 *
psjint(qm0,qmin(3-jj),s2max,2,iqp(3-jj)+1)
2337 * *
psfap(xm0,iqp(jj),1)
2339 gb0=sj0*qm0/qlog*
psuds(qm0,iqp(jj))*1.5
d0 2341 gb0=gb0*xm0**(1.
d0-delh)
2343 gb0=gb0*(1.
d0-xm0)*2.
d0**delh
2346 xmin2=max(.5
d0,xmin)
2348 xmax1=min(xmax,.5
d0)**delh
2349 IF(xmin.GE..5
d0)
THEN 2351 ELSEIF(xmax.LT..5
d0)
THEN 2354 djl=1.
d0/(1.
d0+((2.
d0*xmin)**delh-1.
d0)/delh/
2355 * dlog(2.
d0*(1.
d0-xmax)))
2361 IF(
qsran(b10).GT.djl)
THEN 2362 x=(xmin1+
qsran(b10)*(xmax1-xmin1))**(1.
d0/delh)
2370 qq=qqmin/(1.
d0+
qsran(b10)*(qqmin/qqmax-1.
d0))
2372 IF(debug.GE.2)
WRITE (moniou,253)qq,qqmin,qqmax
2373 253
FORMAT(2
x,
'PSHOT: QQ=',
e10.3,2
x,
'QQMIN=',
e10.3,2
x,
2377 IF(qt2.LT.qt0)
GOTO 7
2386 CALL pscs(ccos,ssin)
2390 pt2=(ept(3)-ep3(3))**2+(ept(4)-ep3(4))**2
2391 s2min2=max(s2min,qmin2)
2393 zmin=(qt2+amj0)/ww/(1.
d0-
x)
2396 s2=
x*(1.
d0-zmin)*ww-pt2
2399 IF(s2.LT.s2min2)
GOTO 7
2401 sj1=
psjint(qq,qmin(3-jj),s2,1,iqp(3-jj)+1)
2403 sj2=
psjint(qq,qmin(3-jj),s2,2,iqp(3-jj)+1)
2410 gb7=(sj1+sj2)/dlog(qt2/alm)*qq*
psuds(qq,iqp(jj))/gb0
2414 gb7=gb7*
x**(1.
d0-delh)
2416 gb7=gb7*(1.
d0-
x)*2.
d0**delh
2419 IF(
qsran(b10).GT.gb7)
GOTO 7
2421 IF(
qsran(b10).LT.sj1/(sj1+sj2))
THEN 2422 IF(iqc(jj).EQ.0)
THEN 2428 eqj(
i,1)=epq(
i+4*(jq-1),jj)
2432 IF(iqc(jj).GT.0)
THEN 2441 ipq(jq,jj)=ipq(1,jj)
2443 135 epq(
i+4*(jq-1),jj)=epq(
i,jj)
2449 IF(iqp(jj).NE.0)
THEN 2452 IF(iqc(jj).GT.0)
THEN 2474 34 eqj(
i,1)=epq(
i+4*(jq-1),jj)
2477 IF(debug.GE.3)
WRITE (moniou,240)jt
2479 CALL pscajet(qt2,iq1,qv1,zv1,qm1,iqv1,
2481 z=(qt2+qm1(1,1))/ww/(1.
d0-
x)
2482 si=
x*(1.
d0-
z)*ww-pt2
2484 IF(si.GT.s2min2)
THEN 2485 iq=min(1,iabs(iqc(jj)))+1
2486 gb=
psjint(qq,qmin(3-jj),si,iq,iqp(3-jj)+1)/
2487 *
psjint(qq,qmin(3-jj),s2,iq,iqp(3-jj)+1)
2488 IF(
qsran(b10).GT.gb)
GOTO 301
2494 wm3=(qt2+qm1(1,1))/wp3
2495 ep3(1)=.5
d0*(wp3+wm3)
2496 ep3(2)=.5
d0*(wp3-wm3)*(3-2*jj)
2498 pt3=dsqrt(ep3(3)**2+ep3(4)**2)
2500 CALL psrec(ep3,qv1,zv1,qm1,iqv1,ldau1,lpar1,iqj,eqj,jfl,jq)
2501 IF(jfl.EQ.0)
GOTO 301
2506 35 epq(
i+4*(jq-1),jj)=eqj(
i,2)
2508 IF(ipc(jq,jj).EQ.0)
THEN 2511 36 epc(
i+4*(jq-1),jj)=eqj(
i,1)
2517 37 epq(
i+4*(2-jq),jj)=eqj(
i,1)
2522 38 epq(
i,jj)=eqj(
i,2)
2524 IF(ipc(jq,jj).EQ.0)
THEN 2527 39 epc(
i+4*(jq-1),jj)=eqj(
i,1)
2531 IF(ipc(jq,jj).EQ.0)
THEN 2534 40 epc(
i+4*(jq-1),jj)=eqj(
i,1)
2539 30 epq(
i,jj)=epq(
i+4,jj)
2543 IF(iabs(iq1).EQ.3)
THEN 2548 IF(debug.GE.2)
WRITE (moniou,209)tyq(iqqq),qt2,ep3
2549 209
FORMAT(2
x,
'PSHOT: NEW JET FLAVOR:',a2,
2550 *
' PT SQUARED FOR THE JET:',
e10.3/
2551 * 4
x,
'JET 4-MOMENTUM:',4
e10.3)
2553 8 ept(
i)=ept(
i)-ep3(
i)
2567 IF(debug.GE.2)
WRITE (moniou,211)si
2568 211
FORMAT(2
x,
'PSHOT: HIGHEST VIRTUALITY SUBPROCESS IN THE LADDER'/
2569 * 4
x,
'MASS SQUARED FOR THE PROCESS:',
e10.3)
2571 xmin=qminn/(qminn+si)
2572 xmin1=.5
d0-dsqrt(.25
d0-(qt0+amj0)/si)
2573 xmin=max(xmin,xmin1)
2576 IF(iqc(1).NE.0.OR.iqc(2).NE.0)
THEN 2577 gb0=tmin**2/dlog(tmin*(1.
d0-xmin)/alm)**2*
2578 *
psfborn(si,tmin,iqc(1),iqc(2))
2580 gb0=.25
d0*si**2/dlog(tmin*(1.
d0-xmin)/alm)**2*
2598 gb=q2**2/dlog(qt2/alm)**2/gb0*
2599 *
psfborn(si,tq,iqc(1),iqc(2))
2600 IF(debug.GE.3)
WRITE (moniou,241)q2,gb
2601 241
FORMAT(2
x,
'PSHOT: Q2=',
e10.3,
' GB=',
e10.3)
2603 IF(
qsran(b10).GT.gb)
GOTO 13
2605 IF(iqc(1).EQ.0.AND.iqc(2).EQ.0)
THEN 2609 51 eqj(
i,1)=epq(
i+4*(jq-1),jm)
2613 IF(ipq(3-jq,jm)*ipq(jq,3-jm).NE.0)
THEN 2616 IF(iabs(ipj).EQ.3)ipj=ipj*4/3
2617 IF(iabs(ipj1).EQ.3)ipj1=ipj1*4/3
2619 epj(
i)=epq(
i+4*(2-jq),jm)
2620 52 epj1(
i)=epq(
i+4*(jq-1),3-jm)
2621 CALL psjdef(ipj,ipj1,epj,epj1,jfl)
2622 IF(jfl.EQ.0)
GOTO 301
2623 ELSEIF(ipq(3-jq,jm).NE.0)
THEN 2624 ipc(jq,3-jm)=ipq(3-jq,jm)
2626 53 epc(
i+4*(jq-1),3-jm)=epq(
i+4*(2-jq),jm)
2627 ELSEIF(ipq(jq,3-jm).NE.0)
THEN 2628 ipc(3-jq,jm)=ipq(jq,3-jm)
2630 54 epc(
i+4*(2-jq),jm)=epq(
i+4*(jq-1),3-jm)
2639 iqj(2)=ipq(3-jq,3-jm)
2641 56 eqj(
i,2)=epq(
i+4*(2-jq),3-jm)
2644 ELSEIF(iqc(1)*iqc(2).EQ.0)
THEN 2645 IF(iqc(1)+iqc(2).GT.0)
THEN 2652 IF(iqc(jm).EQ.0)
THEN 2657 eqj(
i,1)=epq(
i+4*(jq-1),jm)
2660 IF(ipq(3-jq,jm)*ipq(1,3-jm).NE.0)
THEN 2663 IF(iabs(ipj).EQ.3)ipj=ipj*4/3
2664 IF(iabs(ipj1).EQ.3)ipj1=ipj1*4/3
2666 epj(
i)=epq(
i+4*(2-jq),jm)
2667 58 epj1(
i)=epq(
i,3-jm)
2668 CALL psjdef(ipj,ipj1,epj,epj1,jfl)
2669 IF(jfl.EQ.0)
GOTO 301
2670 ELSEIF(ipq(3-jq,jm).NE.0)
THEN 2671 ipc(jq,3-jm)=ipq(3-jq,jm)
2673 59 epc(
i+4*(jq-1),3-jm)=epq(
i+4*(2-jq),jm)
2674 ELSEIF(ipq(1,3-jm).NE.0)
THEN 2675 ipc(3-jq,jm)=ipq(1,3-jm)
2677 60 epc(
i+4*(2-jq),jm)=epq(
i,3-jm)
2686 IF(ipq(1,jm)*ipq(3-jq,3-jm).NE.0)
THEN 2689 IF(iabs(ipj).EQ.3)ipj=ipj*4/3
2690 IF(iabs(ipj1).EQ.3)ipj1=ipj1*4/3
2693 62 epj1(
i)=epq(
i+4*(2-jq),3-jm)
2694 CALL psjdef(ipj,ipj1,epj,epj1,jfl)
2695 IF(jfl.EQ.0)
GOTO 301
2696 ELSEIF(ipq(3-jq,3-jm).NE.0)
THEN 2697 ipc(jq,jm)=ipq(3-jq,3-jm)
2699 63 epc(
i+4*(jq-1),jm)=epq(
i+4*(2-jq),3-jm)
2700 ELSEIF(ipq(1,jm).NE.0)
THEN 2701 ipc(3-jq,3-jm)=ipq(1,jm)
2703 64 epc(
i+4*(2-jq),3-jm)=epq(
i,jm)
2708 IF(iqc(jm).EQ.0)
THEN 2713 eqj(
i,2)=epq(
i+4*(2-jq),jm)
2714 65 eqj(
i,1)=epq(
i,3-jm)
2719 66 eqj(
i,1)=epq(
i+4*(jq-1),3-jm)
2723 ELSEIF(iqc(1)*iqc(2).GT.0)
THEN 2732 67 eqj(
i,1)=epq(
i,3-jm)
2736 IF(iqc(jm).GT.0)
THEN 2745 IF(ipq(1,jm)*ipq(1,3-jm).NE.0)
THEN 2748 IF(iabs(ipj).EQ.3)ipj=ipj*4/3
2749 IF(iabs(ipj1).EQ.3)ipj1=ipj1*4/3
2752 69 epj1(
i)=epq(
i,3-jm)
2753 CALL psjdef(ipj,ipj1,epj,epj1,jfl)
2754 IF(jfl.EQ.0)
GOTO 301
2755 ELSEIF(ipq(1,3-jm).NE.0)
THEN 2756 ipc(jq,jm)=ipq(1,3-jm)
2758 70 epc(
i+4*(jq-1),jm)=epq(
i,3-jm)
2759 ELSEIF(ipq(1,jm).NE.0)
THEN 2760 ipc(3-jq,3-jm)=ipq(1,jm)
2762 71 epc(
i+4*(2-jq),3-jm)=epq(
i,jm)
2770 IF(debug.GE.3)
WRITE (moniou,240)jt
2771 240
FORMAT(2
x,
'PSHOT: COLOUR CONNECTION JT=:',i1)
2773 CALL pscajet(qt2,iqc(jm),qv1,zv1,qm1,iqv1,
2775 CALL pscajet(qt2,iqc(3-jm),qv2,zv2,qm2,iqv2,
2781 IF(dsqrt(si).GT.dsqrt(amt1)+dsqrt(amt2))
THEN 2790 wm3=(qt2+qm1(1,1))/wp3
2791 ep3(1)=.5
d0*(wp3+wm3)
2792 ep3(2)=.5
d0*(wp3-wm3)
2794 CALL pscs(ccos,ssin)
2800 pt3=dsqrt(ep3(3)**2+ep3(4)**2)
2802 CALL psrec(ep3,qv1,zv1,qm1,iqv1,ldau1,lpar1,iqj,eqj,jfl,jq)
2803 IF(jfl.EQ.0)
GOTO 301
2805 if(iabs(iqc(jm)).eq.3)
then 2810 IF(debug.GE.2)
WRITE (moniou,209)tyq(iqqq),qt2
2812 wp3=(1.
d0-
z)*dsqrt(si)
2813 wm3=(qt2+qm2(1,1))/wp3
2814 ep3(1)=.5
d0*(wp3+wm3)
2815 ep3(2)=.5
d0*(wp3-wm3)
2819 pt3=dsqrt(ep3(3)**2+ep3(4)**2)
2822 IF(ipc(jq,jm).EQ.0)
THEN 2825 72 epc(
i+4*(jq-1),jm)=eqj(
i,1)
2829 iqj(2)=ipq(3-jq,3-jm)
2832 73 eqj(
i,2)=epq(
i+4*(2-jq),3-jm)
2835 IF(ipc(jq,jm).EQ.0)
THEN 2838 74 epc(
i+4*(jq-1),jm)=eqj(
i,1)
2840 IF(ipc(3-jq,3-jm).EQ.0)
THEN 2841 ipc(3-jq,3-jm)=iqj(2)
2843 75 epc(
i+4*(2-jq),3-jm)=eqj(
i,2)
2849 eqj(
i,2)=epq(
i+4*(2-jq),jm)
2850 76 eqj(
i,1)=epq(
i+4*(jq-1),3-jm)
2853 IF(ipc(jq,jm).EQ.0)
THEN 2856 77 epc(
i+4*(jq-1),jm)=eqj(
i,1)
2860 78 eqj(
i,1)= eqj(
i,2)
2867 79 eqj(
i,1)=epq(
i+4*(jq-1),3-jm)
2870 IF(ipc(3-jq,jm).EQ.0)
THEN 2873 80 epc(
i+4*(2-jq),jm)=eqj(
i,2)
2875 IF(ipc(jq,3-jm).EQ.0)
THEN 2878 81 epc(
i+4*(jq-1),3-jm)=eqj(
i,1)
2883 82 eqj(
i,1)=epq(
i+4*(jq-1),jm)
2886 IF(ipc(jq,3-jm).EQ.0)
THEN 2889 83 epc(
i+4*(jq-1),3-jm)=eqj(
i,1)
2892 iqj(2)=ipq(3-jq,3-jm)
2895 eqj(
i,2)=epq(
i+4*(2-jq),3-jm)
2896 84 eqj(
i,1)=epq(
i,jm)
2899 IF(ipc(jq,3-jm).EQ.0)
THEN 2902 85 epc(
i+4*(jq-1),3-jm)=eqj(
i,1)
2906 86 eqj(
i,1)= epq(
i,jm)
2909 CALL psrec(ep3,qv2,zv2,qm2,iqv2,ldau2,lpar2,iqj,eqj,jfl,jq2)
2910 IF(jfl.EQ.0)
GOTO 301
2912 if(iabs(iqc(3-jm)).eq.3)
then 2913 iqqq=8+iqc(3-jm)/3*4
2917 IF(debug.GE.2)
WRITE (moniou,209)tyq(iqqq),qt2
2918 IF(debug.GE.2)
WRITE (moniou,212)njtot
2919 212
FORMAT(2
x,
'PSHOT: TOTAL NUMBER OF JETS:',i2)
2922 IF(ipc(3-jq,3-jm).EQ.0)
THEN 2923 ipc(3-jq,3-jm)=iqj(2)
2925 87 epc(
i+4*(2-jq),3-jm)=eqj(
i,2)
2929 IF(ipc(3-jq,jm).EQ.0)
THEN 2932 88 epc(
i+4*(2-jq),jm)=eqj(
i,2)
2934 IF(ipc(jq,3-jm).EQ.0)
THEN 2937 89 epc(
i+4*(jq-1),3-jm)=eqj(
i,1)
2941 IF(ipc(jq,3-jm).EQ.0)
THEN 2944 90 epc(
i+4*(jq-1),3-jm)=eqj(
i,1)
2948 IF(ipc(jq,jm).EQ.0)
THEN 2951 91 epc(
i+4*(jq-1),jm)=eqj(
i,1)
2955 IF(ipc(3-jq,3-jm).EQ.0)
THEN 2956 ipc(3-jq,3-jm)=iqj(2)
2958 92 epc(
i+4*(2-jq),3-jm)=eqj(
i,2)
2960 IF(ipc(jq,jm).EQ.0)
THEN 2963 93 epc(
i+4*(jq-1),jm)=eqj(
i,1)
2967 IF(ipc(jq,jm).EQ.0)
THEN 2970 94 epc(
i+4*(jq-1),jm)=eqj(
i,1)
2975 IF(debug.GE.3)
WRITE (moniou,217)
2976 217
FORMAT(2
x,
'PSHOT - END')
2977 ebal(1)=.5*(wp0+wm0)
2978 ebal(2)=.5*(wp0-wm0)
2984 500 ebal(l)=ebal(l)-epjet(l,
m,
i)
2990 SUBROUTINE psjdef(IPJ,IPJ1,EPJ,EPJ1,JFL)
2995 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
2997 dimension epj(4),epj1(4),ept(4)
2998 COMMON /area10/ stmass,am(7)
2999 COMMON /area43/ moniou
3000 COMMON /debug/ debug
3001 COMMON /area46/ epjet(4,2,15000),ipjet(2,15000)
3002 COMMON /area47/ njtot
3013 IF(debug.GE.2)
WRITE (moniou,201)ipj,ipj1,epj,epj1
3014 201
FORMAT(2
x,
'PSJDEF: PARTON FLAVORS',
3015 *
': IPJ=',i2,2
x,
'IPJ1=',i2/
3016 * 4
x,
'PARTON 4-MOMENTA:',2
x,4(
e10.3,1
x))
3018 1 ept(
i)=epj(
i)+epj1(
i)
3023 IF(iabs(ipj).LE.2)
THEN 3025 ELSEIF(iabs(ipj).EQ.4)
THEN 3030 IF(iabs(ipj1).LE.2)
THEN 3032 ELSEIF(iabs(ipj1).EQ.4)
THEN 3047 IF( njtot . gt. 15000 )
THEN 3048 WRITE(moniou,*)
'PSJDEF: TOO MANY JETS' 3049 WRITE(moniou,*)
'PSJDEF: NJTOT = ',njtot
3055 epjet(
i,1,njtot)=epj(
i)
3056 2 epjet(
i,2,njtot)=epj1(
i)
3058 IF(debug.GE.3)
WRITE (moniou,202)
3059 202
FORMAT(2
x,
'PSJDEF - END')
3064 FUNCTION psjet(Q1,Q2,S,S2MIN,J,L)
3074 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
3076 COMMON /area6/
pi,bm,am
3077 COMMON /area17/ del,rs,rs0,fs,alf,rr,sh,delh
3078 COMMON /area18/ alm,qt0,qlog,qll,aqt0,qtf,bet,amj0
3079 common/ar3/
x1(7),a1(7)
3080 COMMON /area43/ moniou
3081 COMMON /debug/ debug
3084 IF(debug.GE.2)
WRITE (moniou,201)s,q1,q2,s2min,
j,l
3085 201
FORMAT(2
x,
'PSJET - UNORDERED LADDER CROSS SECTION:'/
3086 * 4
x,
'S=',
e10.3,2
x,
'Q1=',
e10.3,2
x,
'Q2=',
e10.3,2
x,
'S2MIN=',
3087 *
e10.3,2
x,
'J=',i1,2
x,
'L=',i1)
3090 p=dsqrt(1.
d0-3.
d0*qt0/s)
3091 cosf=(1.
d0-18.
d0*qt0/s)/
p**3
3092 fi=atan(1.
d0/cosf**2-1.
d0)
3093 IF(cosf.LT.0.
d0)fi=
pi-fi
3095 zmax=(2.
d0-
p*(dsqrt(3.
d0)*sin(fi)-
cos(fi)))/3.
d0 3098 IF(qt0/(1.
d0-zmin)**2.LT.s2min)
3099 * zmin=.5
d0*(1.
d0+s2min/s-dsqrt((1.
d0-s2min/s)**2-4.
d0*qt0/s))
3102 IF(1.
d0-zmin.LT.dsqrt(qt0/q1))
THEN 3103 qmin=qt0/(1.
d0-zmin)**2
3108 qmax=qt0/(1.
d0-zmax)**2
3112 IF(debug.GE.3)
WRITE (moniou,203)qmin,qmax
3113 203
FORMAT(2
x,
'PSJET:',2
x,
'QMIN=',
e10.3,2
x,
'QMAX=',
e10.3)
3114 IF(qmax.GT.qmin)
THEN 3120 qi=2.
d0*qmin/(1.
d0+qmin/qmax+(2*
m-3)*
x1(
i)*(1.
d0-qmin/qmax))
3122 zmax=(1.
d0-dsqrt(qt0/qi))**delh
3123 zmin=((qi+max(qi,s2min))/(qi+s))**delh
3127 IF(debug.GE.3)
WRITE (moniou,204)qi,zmin,zmax
3128 204
FORMAT(2
x,
'PSJET:',2
x,
'QI=',
e10.3,2
x,
'ZMIN=',
e10.3,2
x,
3130 IF(zmax.GT.zmin)
THEN 3133 z=(.5
d0*(zmax+zmin+(2*m1-3)*
x1(i1)*(zmax-zmin)))**
3141 2 fsj=fsj+a1(i1)*sj/dlog(qt/alm)/
z**delh
3148 IF(debug.GE.3)
WRITE (moniou,202)
psjet 3149 202
FORMAT(2
x,
'PSJET=',
e10.3)
3154 FUNCTION psjet1(Q1,Q2,S,S2MIN,J,L)
3164 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
3166 COMMON /area6/
pi,bm,am
3167 COMMON /area17/ del,rs,rs0,fs,alf,rr,sh,delh
3168 COMMON /area18/ alm,qt0,qlog,qll,aqt0,qtf,bet,amj0
3169 common/ar3/
x1(7),a1(7)
3170 COMMON /area43/ moniou
3171 COMMON /debug/ debug
3174 IF(debug.GE.2)
WRITE (moniou,201)s,q1,q2,s2min,
j,l
3175 201
FORMAT(2
x,
'PSJET1 - STRICTLY ORDERED LADDER CROSS SECTION:'/
3176 * 4
x,
'S=',
e10.3,2
x,
'Q1=',
e10.3,2
x,
'Q2=',
e10.3,2
x,
'S2MIN=',
3177 *
e10.3,2
x,
'J=',i1,2
x,
'L=',i1)
3180 p=dsqrt(1.
d0-3.
d0*qt0/s)
3181 cosf=(1.
d0-18.
d0*qt0/s)/
p**3
3182 fi=atan(1.
d0/cosf**2-1.
d0)
3183 IF(cosf.LT.0.
d0)fi=
pi-fi
3185 zmax=(2.
d0-
p*(dsqrt(3.
d0)*sin(fi)-
cos(fi)))/3.
d0 3188 IF(qt0/(1.
d0-zmin)**2.LT.s2min)
3189 * zmin=.5
d0*(1.
d0+s2min/s-dsqrt((1.
d0-s2min/s)**2-4.
d0*qt0/s))
3192 IF(1.
d0-zmin.LT.dsqrt(qt0/q1))
THEN 3193 qmin=qt0/(1.
d0-zmin)**2
3198 qmax=qt0/(1.
d0-zmax)**2
3202 IF(debug.GE.3)
WRITE (moniou,203)qmin,qmax
3203 203
FORMAT(2
x,
'PSJET1:',2
x,
'QMIN=',
e10.3,2
x,
'QMAX=',
e10.3)
3204 IF(qmax.GT.qmin)
THEN 3210 qi=2.
d0*qmin/(1.
d0+qmin/qmax+(2*
m-3)*
x1(
i)*(1.
d0-qmin/qmax))
3212 zmax=(1.
d0-dsqrt(qt0/qi))**delh
3213 zmin=((qi+max(qi,s2min))/(qi+s))**delh
3217 IF(debug.GE.3)
WRITE (moniou,204)qi,zmin,zmax
3218 204
FORMAT(2
x,
'PSJET1:',2
x,
'QI=',
e10.3,2
x,
'ZMIN=',
e10.3,2
x,
3220 IF(zmax.GT.zmin)
THEN 3223 z=(.5
d0*(zmax+zmin+(2*m1-3)*
x1(i1)*(zmax-zmin)))**
3232 2 fsj=fsj+a1(i1)*sj/dlog(qt/alm)/
z**delh
3239 IF(debug.GE.3)
WRITE (moniou,202)
psjet1 3240 202
FORMAT(2
x,
'PSJET1=',
e10.3)
3245 FUNCTION psjint(Q1,Q2,S,M,L)
3254 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
3256 dimension wi(3),wj(3),wk(3)
3257 COMMON /area18/ alm,qt0,qlog,qll,aqt0,qtf,bet,amj0
3258 COMMON /area29/ csj(17,17,68)
3259 COMMON /area43/ moniou
3260 COMMON /debug/ debug
3263 IF(debug.GE.2)
WRITE (moniou,201)s,q1,q2,
m,l
3264 201
FORMAT(2
x,
'PSJINT - UNORDERED LADDER CROSS SECTION INTERPOL.:'/
3265 * 4
x,
'S=',
e10.3,2
x,
'Q1=',
e10.3,2
x,
'Q2=',
e10.3,2
x,
3266 * 2
x,
'M=',i1,2
x,
'L=',i1)
3269 IF(s.LE.max(4.
d0*qt0,qq))
THEN 3270 IF(debug.GE.3)
WRITE (moniou,202)
psjint 3271 202
FORMAT(2
x,
'PSJINT=',
e10.3)
3275 ml=17*(
m-1)+34*(l-1)
3276 qli=dlog(q1/qt0)/1.38629
d0 3277 qlj=dlog(q2/qt0)/1.38629
d0 3278 sl=dlog(s/qt0)/1.38629
d0 3286 IF(sql.GT.10.
d0)
THEN 3291 wi(3)=wi(2)*(wi(2)-1.
d0)*.5
d0 3292 wi(1)=1.
d0-wi(2)+wi(3)
3293 wi(2)=wi(2)-2.
d0*wi(3)
3295 wj(3)=wj(2)*(wj(2)-1.
d0)*.5
d0 3296 wj(1)=1.
d0-wj(2)+wj(3)
3297 wj(2)=wj(2)-2.
d0*wj(3)
3299 wk(3)=wk(2)*(wk(2)-1.
d0)*.5
d0 3300 wk(1)=1.
d0-wk(2)+wk(3)
3301 wk(2)=wk(2)-2.
d0*wk(3)
3308 ELSEIF(sql.LT.1.
d0.AND.
i+
j.NE.0)
THEN 3309 sq=(s/max(q1,q2)-1.
d0)/3.
d0 3311 wi(3)=wi(2)*(wi(2)-1.
d0)*.5
d0 3312 wi(1)=1.
d0-wi(2)+wi(3)
3313 wi(2)=wi(2)-2.
d0*wi(3)
3315 wj(3)=wj(2)*(wj(2)-1.
d0)*.5
d0 3316 wj(1)=1.
d0-wj(2)+wj(3)
3317 wj(2)=wj(2)-2.
d0*wj(3)
3327 sq=(s/qt0/4.
d0-1.
d0)/3.
d0 3339 IF(
i+kl.GT.12)
i=12-kl
3340 IF(
j+kl.GT.12)
j=12-kl
3343 wi(3)=wi(2)*(wi(2)-1.
d0)*.5
d0 3344 wi(1)=1.
d0-wi(2)+wi(3)
3345 wi(2)=wi(2)-2.
d0*wi(3)
3347 wj(3)=wj(2)*(wj(2)-1.
d0)*.5
d0 3348 wj(1)=1.
d0-wj(2)+wj(3)
3349 wj(2)=wj(2)-2.
d0*wj(3)
3351 wk(3)=wk(2)*(wk(2)-1.
d0)*.5
d0 3352 wk(1)=1.
d0-wk(2)+wk(3)
3353 wk(2)=wk(2)-2.
d0*wk(3)
3360 k2=max(i2,j2)+kl+k1-1+ml
3368 wi(3)=wi(2)*(wi(2)-1.
d0)*.5
d0 3369 wi(1)=1.
d0-wi(2)+wi(3)
3370 wi(2)=wi(2)-2.
d0*wi(3)
3372 wj(3)=wj(2)*(wj(2)-1.
d0)*.5
d0 3373 wj(1)=1.
d0-wj(2)+wj(3)
3374 wj(2)=wj(2)-2.
d0*wj(3)
3384 IF(debug.GE.3)
WRITE (moniou,202)
psjint 3389 SUBROUTINE psjint0(S,SJ,SJB,M,L)
3398 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
3401 COMMON /area18/ alm,qt0,qlog,qll,aqt0,qtf,bet,amj0
3402 COMMON /area32/ csj(17,2,2),csb(17,2,2)
3403 COMMON /area43/ moniou
3404 COMMON /debug/ debug
3407 IF(debug.GE.2)
WRITE (moniou,201)s,
m,l
3408 201
FORMAT(2
x,
'PSJINT0 - HARD CROSS SECTION INTERPOLATION:'/
3409 * 4
x,
'S=',
e10.3,2
x,
'M=',i1,2
x,
'L=',i1)
3412 IF(s.LE.4.
d0*qt0)
THEN 3413 IF(debug.GE.3)
WRITE (moniou,202)sj,sjb
3414 202
FORMAT(2
x,
'PSJINT0: SJ=',
e10.3,2
x,
'SJB=',
e10.3)
3418 sl=dlog(s/qt0)/1.38629
d0 3421 sq=(s/qt0/4.
d0-1.
d0)/3.
d0 3422 sjb=exp(csb(3,
m+1,l+1))*sq
3423 sj=exp(csj(3,
m+1,l+1))*sq
3427 wk(3)=wk(2)*(wk(2)-1.
d0)*.5
d0 3428 wk(1)=1.
d0-wk(2)+wk(3)
3429 wk(2)=wk(2)-2.
d0*wk(3)
3432 sj=sj+csj(k+k1,
m+1,l+1)*wk(k1)
3433 1 sjb=sjb+csb(k+k1,
m+1,l+1)*wk(k1)
3437 IF(debug.GE.3)
WRITE (moniou,202)sj,sjb
3451 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
3453 dimension wi(3),wj(3),wk(3)
3454 COMMON /area18/ alm,qt0,qlog,qll,aqt0,qtf,bet,amj0
3455 COMMON /area30/ csj(17,17,68)
3456 COMMON /area43/ moniou
3457 COMMON /debug/ debug
3460 IF(debug.GE.2)
WRITE (moniou,201)s,q1,q2,
m,l
3461 201
FORMAT(2
x,
'PSJINT1 - STRICTLY ORDERED LADDER CROSS SECTION',
3462 *
' INTERPOLATION:'/
3463 * 4
x,
'S=',
e10.3,2
x,
'Q1=',
e10.3,2
x,
'Q2=',
e10.3,2
x,
3464 * 4
x,
'M=',i1,2
x,
'L=',i1)
3467 IF(s.LE.max(4.
d0*qt0,qq))
THEN 3468 IF(debug.GE.3)
WRITE (moniou,202)
psjint1 3469 202
FORMAT(2
x,
'PSJINT1=',
e10.3)
3473 ml=17*(
m-1)+34*(l-1)
3474 qli=dlog(q1/qt0)/1.38629
d0 3475 qlj=dlog(q2/qt0)/1.38629
d0 3476 sl=dlog(s/qt0)/1.38629
d0 3484 IF(sql.GT.10.
d0)
THEN 3489 wi(3)=wi(2)*(wi(2)-1.
d0)*.5
d0 3490 wi(1)=1.
d0-wi(2)+wi(3)
3491 wi(2)=wi(2)-2.
d0*wi(3)
3493 wj(3)=wj(2)*(wj(2)-1.
d0)*.5
d0 3494 wj(1)=1.
d0-wj(2)+wj(3)
3495 wj(2)=wj(2)-2.
d0*wj(3)
3497 wk(3)=wk(2)*(wk(2)-1.
d0)*.5
d0 3498 wk(1)=1.
d0-wk(2)+wk(3)
3499 wk(2)=wk(2)-2.
d0*wk(3)
3506 ELSEIF(sql.LT.1.
d0.AND.
i+
j.NE.0)
THEN 3507 sq=(s/max(q1,q2)-1.
d0)/3.
d0 3509 wi(3)=wi(2)*(wi(2)-1.
d0)*.5
d0 3510 wi(1)=1.
d0-wi(2)+wi(3)
3511 wi(2)=wi(2)-2.
d0*wi(3)
3513 wj(3)=wj(2)*(wj(2)-1.
d0)*.5
d0 3514 wj(1)=1.
d0-wj(2)+wj(3)
3515 wj(2)=wj(2)-2.
d0*wj(3)
3525 sq=(s/qt0/4.
d0-1.
d0)/3.
d0 3537 IF(
i+kl.GT.12)
i=12-kl
3538 IF(
j+kl.GT.12)
j=12-kl
3542 wi(3)=wi(2)*(wi(2)-1.
d0)*.5
d0 3543 wi(1)=1.
d0-wi(2)+wi(3)
3544 wi(2)=wi(2)-2.
d0*wi(3)
3546 wj(3)=wj(2)*(wj(2)-1.
d0)*.5
d0 3547 wj(1)=1.
d0-wj(2)+wj(3)
3548 wj(2)=wj(2)-2.
d0*wj(3)
3550 wk(3)=wk(2)*(wk(2)-1.
d0)*.5
d0 3551 wk(1)=1.
d0-wk(2)+wk(3)
3552 wk(2)=wk(2)-2.
d0*wk(3)
3559 k2=max(i2,j2)+kl+k1-1+ml
3567 wi(3)=wi(2)*(wi(2)-1.
d0)*.5
d0 3568 wi(1)=1.
d0-wi(2)+wi(3)
3569 wi(2)=wi(2)-2.
d0*wi(3)
3571 wj(3)=wj(2)*(wj(2)-1.
d0)*.5
d0 3572 wj(1)=1.
d0-wj(2)+wj(3)
3573 wj(2)=wj(2)-2.
d0*wj(3)
3583 IF(debug.GE.3)
WRITE (moniou,202)
psjint1 3588 FUNCTION pslam(S,A,B)
3594 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
3596 COMMON /area43/ moniou
3597 COMMON /debug/ debug
3600 IF(debug.GE.2)
WRITE (moniou,201)s,
a,
b 3601 201
FORMAT(2
x,
'PSLAM - KINEMATICAL FUNCTION S=',
e10.3,2
x,
'A=',
3604 IF(debug.GE.3)
WRITE (moniou,202)
pslam 3605 202
FORMAT(2
x,
'PSLAM=',
e10.3)
3613 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
3616 COMMON /area43/ moniou
3617 COMMON /debug/ debug
3620 IF(debug.GE.2)
WRITE (moniou,201)ep
3621 201
FORMAT(2
x,
'PSNORM - 4-VECTOR SQUARED FOR ',
3622 *
'EP=',4(
e10.3,1
x))
3626 IF(debug.GE.3)
WRITE (moniou,202)
psnorm 3627 202
FORMAT(2
x,
'PSNORM=',
e10.3)
3632 SUBROUTINE psrec(EP,QV,ZV,QM,IQV,LDAU,LPAR,IQJ,EQJ,JFL,JQ)
3647 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
3649 dimension ep(4),ep3(4),epv(4,30,50),
3650 * qv(30,50),zv(30,50),qm(30,50),iqv(30,50),
3651 * ldau(30,49),lpar(30,50),
3652 * iqj(2),eqj(4,2),ipq(2,30,50),epq(8,30,50),
3654 COMMON /area43/ moniou
3655 COMMON /debug/ debug
3658 IF(debug.GE.2)
WRITE (moniou,201)jq,ep,iqj
3659 201
FORMAT(2
x,
'PSREC - JET RECONSTRUCTURING: JQ=',i1/
3660 * 4
x,
'JET 4-MOMENTUM EP=',4(
e10.3,1
x)/4
x,
'IQJ=',2i2)
3664 1 epq(
i,1,1)=eqj(
i,1)
3667 IF(iqv(1,1).EQ.0)
THEN 3669 2 epq(
i+4,1,1)=eqj(
i,2)
3678 IF(qv(nrow,nlev).EQ.0.
d0)
THEN 3681 IF(epq(1,nrow,nlev).NE.0.
d0)
THEN 3682 IF(iabs(ipj).EQ.3)ipj=ipj*4/3
3684 epj(
i)=epv(
i,nrow,nlev)
3685 4 epj1(
i)=epq(
i,nrow,nlev)
3686 ipj1=ipq(1,nrow,nlev)
3687 IF(iabs(ipj1).EQ.3)ipj1=ipj1*4/3
3688 CALL psjdef(ipj,ipj1,epj,epj1,jfl)
3689 IF(debug.GE.3)
WRITE (moniou,211)ipj,ipj1,jfl
3690 211
FORMAT(2
x,
'PSREC - NEW STRING FLAVOURS: ',2i3,
' JFL=',i1)
3693 ipq(1,nrow,nlev)=ipj
3695 5 epq(
i,nrow,nlev)=epv(
i,nrow,nlev)
3696 IF(debug.GE.3)
WRITE (moniou,212)ipj,
3697 * (epv(
i,nrow,nlev),
i=1,4),jfl
3698 212
FORMAT(2
x,
'PSREC: NEW FINAL JET FLAVOR: ',i3,2
x,
3699 *
'JET 4-MOMENTUM:', 4(
e10.3,1
x),
' JFL=',i1)
3705 6 epj(
i)=.5
d0*epv(
i,nrow,nlev)
3708 IF(epq(1+4*(
m-1),nrow,nlev).NE.0.
d0)
THEN 3710 7 epj1(
i)=epq(4*(
m-1)+
i,nrow,nlev)
3711 ipj1=ipq(
m,nrow,nlev)
3712 IF(iabs(ipj1).EQ.3)ipj1=ipj1*4/3
3713 CALL psjdef(ipj,ipj1,epj,epj1,jfl)
3716 ipq(
m,nrow,nlev)=ipj
3718 8 epq(4*(
m-1)+
i,nrow,nlev)=epj(
i)
3723 IF(debug.GE.3)
WRITE (moniou,204)nlev,nrow,iqv(nrow,nlev),
3724 * (epv(
i,nrow,nlev),
i=1,4)
3725 204
FORMAT(2
x,
'PSREC: FINAL JET AT LEVEL NLEV=',i2,
3726 *
' NROW=',i2/4
x,
'JET FLAVOR: ',i3,2
x,
'JET 4-MOMENTUM:',
3731 10 ep3(
i)=epv(
i,nrow,nlev)
3734 qt2=(
z*(1.
d0-
z))**2*qv(nrow,nlev)
3735 ldrow=ldau(nrow,nlev)
3739 wmi=(qt2+qm(ldrow,nlev+1))/wpi
3740 ep3(1)=.5
d0*(wpi+wmi)
3741 ep3(2)=.5
d0*(wpi-wmi)
3746 CALL psrotat(ep3,s0x,c0x,s0,c0)
3749 11 epv(
i,ldrow,nlev+1)=ep3(
i)
3750 IF(debug.GE.3)
WRITE (moniou,206)nlev+1,ldrow,ep3
3751 206
FORMAT(2
x,
'PSREC: JET AT LEVEL NLEV=',i2,
3752 *
' NROW=',i2/4
x,
'JET 4-MOMENTUM:',4(
e10.3,1
x))
3755 wmi=(qt2+qm(ldrow+1,nlev+1))/wpi
3756 ep3(1)=.5
d0*(wpi+wmi)
3757 ep3(2)=.5
d0*(wpi-wmi)
3760 CALL psrotat(ep3,s0x,c0x,s0,c0)
3761 IF(debug.GE.3)
WRITE (moniou,206)nlev+1,ldrow+1,ep3
3764 12 epv(
i,ldrow+1,nlev+1)=ep3(
i)
3766 IF(iqv(nrow,nlev).EQ.0)
THEN 3767 IF(iqv(ldrow,nlev+1).NE.0)
THEN 3768 ipq(1,ldrow,nlev+1)=ipq(1,nrow,nlev)
3769 ipq(1,ldrow+1,nlev+1)=ipq(2,nrow,nlev)
3771 epq(
i,ldrow,nlev+1)=epq(
i,nrow,nlev)
3772 13 epq(
i,ldrow+1,nlev+1)=epq(
i+4,nrow,nlev)
3774 ipq(1,ldrow,nlev+1)=ipq(1,nrow,nlev)
3775 ipq(2,ldrow,nlev+1)=0
3776 ipq(1,ldrow+1,nlev+1)=0
3777 ipq(2,ldrow+1,nlev+1)=ipq(2,nrow,nlev)
3779 epq(
i,ldrow,nlev+1)=epq(
i,nrow,nlev)
3780 epq(
i+4,ldrow,nlev+1)=0.
d0 3781 epq(
i,ldrow+1,nlev+1)=0.
d0 3782 14 epq(
i+4,ldrow+1,nlev+1)=epq(
i+4,nrow,nlev)
3785 IF(iqv(ldrow,nlev+1).EQ.0)
THEN 3786 ipq(1,ldrow,nlev+1)=ipq(1,nrow,nlev)
3787 ipq(2,ldrow,nlev+1)=0
3788 ipq(1,ldrow+1,nlev+1)=0
3790 epq(
i,ldrow,nlev+1)=epq(
i,nrow,nlev)
3791 epq(
i+4,ldrow,nlev+1)=0.
d0 3792 15 epq(
i,ldrow+1,nlev+1)=0.
d0 3794 ipq(1,ldrow,nlev+1)=0
3795 ipq(1,ldrow+1,nlev+1)=0
3796 ipq(2,ldrow+1,nlev+1)=ipq(1,nrow,nlev)
3798 epq(
i,ldrow,nlev+1)=0.
d0 3799 epq(
i,ldrow+1,nlev+1)=0.
d0 3800 16 epq(
i+4,ldrow+1,nlev+1)=epq(
i,nrow,nlev)
3813 18 eqj(
i,1)=epq(
i,1,1)
3814 IF(iqv(1,1).EQ.0)
THEN 3817 19 eqj(
i,2)=epq(
i+4,1,1)
3819 IF(debug.GE.3)
WRITE (moniou,202)iqj
3820 202
FORMAT(2
x,
'PSREC - END',2
x,
'iqj=',2i2)
3824 lprow=lpar(nrow,nlev)
3826 IF(ldau(lprow,nlev-1).EQ.nrow)
THEN 3827 IF(iqv(nrow,nlev).EQ.0)
THEN 3828 IF(epq(1,lprow,nlev-1).EQ.0.
d0)
THEN 3829 ipq(1,lprow,nlev-1)=ipq(1,nrow,nlev)
3831 20 epq(
i,lprow,nlev-1)=epq(
i,nrow,nlev)
3833 ipq(1,nrow+1,nlev)=ipq(2,nrow,nlev)
3835 21 epq(
i,nrow+1,nlev)=epq(
i+4,nrow,nlev)
3837 IF(iqv(lprow,nlev-1).EQ.0)
THEN 3838 IF(epq(1,lprow,nlev-1).EQ.0.
d0)
THEN 3839 ipq(1,lprow,nlev-1)=ipq(1,nrow,nlev)
3841 22 epq(
i,lprow,nlev-1)=epq(
i,nrow,nlev)
3844 ipq(1,nrow+1,nlev)=ipq(1,nrow,nlev)
3846 23 epq(
i,nrow+1,nlev)=epq(
i,nrow,nlev)
3853 IF(iqv(nrow,nlev).EQ.0)
THEN 3854 IF(iqv(lprow,nlev-1).EQ.0)
THEN 3855 IF(epq(5,lprow,nlev-1).EQ.0.
d0)
THEN 3856 ipq(2,lprow,nlev-1)=ipq(2,nrow,nlev)
3858 24 epq(
i+4,lprow,nlev-1)=epq(
i+4,nrow,nlev)
3861 IF(epq(1,lprow,nlev-1).EQ.0.
d0)
THEN 3862 ipq(1,lprow,nlev-1)=ipq(2,nrow,nlev)
3864 25 epq(
i,lprow,nlev-1)=epq(
i+4,nrow,nlev)
3868 IF(iqv(lprow,nlev-1).EQ.0.AND.
3869 * epq(5,lprow,nlev-1).EQ.0.
d0)
THEN 3870 ipq(2,lprow,nlev-1)=ipq(1,nrow,nlev)
3872 26 epq(
i+4,lprow,nlev-1)=epq(
i,nrow,nlev)
3890 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
3892 COMMON /area17/ del,rs,rs0,fs,alf,rr,sh,delh
3893 COMMON /area18/ alm,qt0,qlog,qll,aqt0,qtf,bet,amj0
3894 COMMON /ar3/
x1(7),a1(7)
3895 COMMON /area43/ moniou
3896 COMMON /debug/ debug
3899 IF(debug.GE.2)
WRITE (moniou,201)s,
z,iqq
3900 201
FORMAT(2
x,
'PSREJS - REJECTION FUNCTION TABULATION: '/
3901 * 4
x,
'S=',
e10.3,2
x,
'Z=',
e10.3,2
x,
'IQQ=',i1)
3902 xmin=4.
d0*(qt0+amj0)/s
3903 xmin=xmin**(delh-del)
3909 z1=(.5
d0*(1.
d0+xmin-(2*
m-3)*
x1(
i)*(1.
d0-xmin)))**(1.
d0/
3917 CALL psjint0(z1*s,sj,sjb,iqq,0)
3930 1 st2=st2+a1(
j)*((1.
d0-z1**(.5
d0*(1.
d0+
x1(
j))))*
3937 IF(debug.GE.2)
WRITE (moniou,202)
psrejs 3938 202
FORMAT(2
x,
'PSREJS=',
e10.3)
3948 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
3950 COMMON /area17/ del,rs,rs0,fs,alf,rr,sh,delh
3951 COMMON /area18/ alm,qt0,qlog,qll,aqt0,qtf,bet,amj0
3952 COMMON /area43/ moniou
3953 COMMON /debug/ debug
3956 IF(debug.GE.2)
WRITE (moniou,201)s
3957 201
FORMAT(2
x,
'PSREJV - REJECTION FUNCTION TABULATION: ',
3970 IF(debug.GE.3)
WRITE (moniou,202)
psrejv 3971 202
FORMAT(2
x,
'PSREJV=',
e10.3)
3982 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
3985 COMMON /area1/ ia(2),icz,icp
3986 COMMON /area17/ del,rs,rs0,fs,alf,rr,sh,delh
3987 COMMON /area18/ alm,qt0,qlog,qll,aqt0,qtf,bet,amj0
3988 COMMON /area23/ rjv(50)
3989 COMMON /area24/ rjs(50,5,10)
3990 COMMON /area43/ moniou
3991 COMMON /debug/ debug
3994 IF(debug.GE.2)
WRITE (moniou,201)yj,
z0,iqq
3995 201
FORMAT(2
x,
'PSRJINT - REJECTION FUNCTION INTERPOLATION:'/
3996 * 4
x,
'YJ=',
e10.3,2
x,
'Z0=',
e10.3,2
x,
'IQQ=',i1)
4003 psrjint=exp(rjv(1))*yy+(exp(rjv(2))-2.
d0*
4004 * exp(rjv(1)))*yy*(yy-1.
d0)*.5
d0 4006 psrjint=exp(rjv(jy)+(rjv(jy+1)-rjv(jy))*(yy-jy)
4007 * +(rjv(jy+2)+rjv(jy)-2.
d0*rjv(jy+1))*(yy-jy)*
4008 * (yy-jy-1.
d0)*.5
d0)
4012 iq=(iqq+1)/2+1+2*(icz-1)
4026 a(
i)=exp(rjs(1,j1,iq))*yy+(exp(rjs(2,j1,iq))-2.
d0*
4027 * exp(rjs(1,j1,iq)))*yy*(yy-1.
d0)*.5
d0 4028 IF(
a(
i).GT.0.
d0)
THEN 4034 a(
i)=rjs(jy,j1,iq)+(rjs(jy+1,j1,iq)-
4035 * rjs(jy,j1,iq))*(yy-jy)
4036 * +(rjs(jy+2,j1,iq)+rjs(jy,j1,iq)-2.
d0*
4037 * rjs(jy+1,j1,iq))*(yy-jy)*(yy-jy-1.
d0)*.5
d0 4050 IF(debug.GE.3)
WRITE (moniou,202)
psrjint 4051 202
FORMAT(2
x,
'PSRJINT=',
e10.3)
4056 FUNCTION psroot(QLMAX,G,J)
4064 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
4066 COMMON /area43/ moniou
4067 COMMON /debug/ debug
4070 IF(debug.GE.2)
WRITE (moniou,201)qlmax,
g,
j 4071 201
FORMAT(2
x,
'PSQINT - BRANCHING MOMENTUM TABULATION:'/
4072 * 4
x,
'QLMAX=',
e10.3,2
x,
'G=',
e10.3,2
x,
'J=',i1)
4079 1 ql2=ql1-(ql1-ql0)*f1/(f1-f0)
4083 ELSEIF(ql2.GT.qlmax)
THEN 4090 IF(abs(f2).GT.1.
d-3)
THEN 4099 IF(debug.GE.3)
WRITE (moniou,202)
psroot 4100 202
FORMAT(2
x,
'PSROOT=',
e10.3)
4105 SUBROUTINE psrotat(EP,S0X,C0X,S0,C0)
4108 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
4110 dimension ep(4),ep1(3)
4111 COMMON /area43/ moniou
4112 COMMON /debug/ debug
4115 IF(debug.GE.2)
WRITE (moniou,201)ep,s0x,c0x,s0,c0
4116 201
FORMAT(2
x,
'PSROTAT - SPACIAL ROTATION:'/4
x,
4117 *
'4-VECTOR EP=',4(
e10.3,1
x)/4
x,
'S0X=',
e10.3,
'C0X=',
e10.3,
4120 ep1(2)=ep(2)*s0+ep(3)*c0
4121 ep1(1)=ep(2)*c0-ep(3)*s0
4124 ep(4)=ep1(2)*s0x+ep1(3)*c0x
4125 ep(3)=ep1(2)*c0x-ep1(3)*s0x
4126 IF(debug.GE.3)
WRITE (moniou,202)ep
4127 202
FORMAT(2
x,
'PSROTAT: ROTATED 4-VECTOR EP=',
4133 FUNCTION psqint(QLMAX,G,J)
4140 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
4142 dimension wi(3),wk(3)
4143 COMMON /area18/ alm,qt0,qlog,qll,aqt0,qtf,bet,amj0
4144 COMMON /area34/ qrt(10,101,2)
4145 COMMON /area43/ moniou
4146 COMMON /debug/ debug
4149 IF(debug.GE.2)
WRITE (moniou,201)qlmax,
g,
j 4150 201
FORMAT(2
x,
'PSQINT - BRANCHING MOMENTUM INTERPOLATION:'/
4151 * 4
x,
'QLMAX=',
e10.3,2
x,
'G=',
e10.3,2
x,
'J=',i1)
4154 sl=100.
d0*dlog(1.
d0-
g*(1.
d0-sud0))/dlog(sud0)
4159 wk(3)=wk(2)*(wk(2)-1.
d0)*.5
d0 4160 wk(1)=1.
d0-wk(2)+wk(3)
4161 wk(2)=wk(2)-2.
d0*wk(3)
4166 wi(3)=wi(2)*(wi(2)-1.
d0)*.5
d0 4167 wi(1)=1.
d0-wi(2)+wi(3)
4168 wi(2)=wi(2)-2.
d0*wi(3)
4175 IF(debug.GE.3)
WRITE (moniou,202)
psqint 4176 202
FORMAT(2
x,
'PSQINT=',
e10.3)
4181 SUBROUTINE psshar(LS,NHP,NW,NT)
4188 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
4191 dimension wp(56),wm(56),wha(4000),whb(4000),lha0(56),
4192 * lhb0(56),izp(56),izt(56),wp0h(56),wm0h(56),
4193 * wpp(2),wmm(2),ep3(4),lqa0(56),lqb0(56),ipc(2,2),epc(8,2),
4194 * ila(56),ilb(56),ela(4,56),elb(4,56),ep(4),ep1(4)
4195 COMMON /area1/ ia(2),icz,icp
4196 COMMON /area2/ s,y0,wp0,wm0
4197 COMMON /area9/ lqa(56),lqb(56),nqs(1000),ias(1000),ibs(1000),
4198 * lha(56),lhb(56),zh(4000),iah(4000),ibh(4000),
4199 * iqh(4000),lva(56),lvb(56)
4200 COMMON /area10/ stmass,am(7)
4203 COMMON /area17/ del,rs,rs0,fs,alfp,rr,sh,delh
4204 COMMON /area18/ alm,qt0,qlog,qll,aqt0,qtf,bet,amj0
4205 COMMON /area19/ ahl(5)
4206 COMMON /area20/ wppp
4207 COMMON /area25/ ahv(5)
4208 COMMON /area43/ moniou
4209 COMMON /debug/ debug
4210 COMMON /area47/ njtot
4215 IF(debug.GE.1)
WRITE (moniou,201)nw,nt,nhp,ls
4216 201
FORMAT(2
x,
'PSSHARE - ENERGY SHARING PROCEDURE'/
4217 * 4
x,
'NUMBER OF WOUNDED PROJECTILE NUCLEONS(HADRONS) NW=',i2/
4218 * 4
x,
'NUMBER OF TARGET NUCLEONS IN THE ACTIVE STATE NT=',i2/
4219 * 4
x,
'NUMBER OF SEMIHARD BLOCKS NHP=',i3/
4220 * 4
x,
'NUMBER OF SOFT POMERON BLOCKS LS=',i3)
4238 1 izp(
i)=int(2.5+
qsran(b10))
4246 2 izt(
i)=int(2.5+
qsran(b10))
4259 IF(debug.GE.3)
WRITE (moniou,211)nhp
4260 211
FORMAT(2
x,
'PSSHARE: NUMBER OF HARD POMERONS NHP=',i3)
4277 IF(debug.GE.3)
WRITE (moniou,212)ih
4278 212
FORMAT(2
x,
'PSSHARE: GBH-INI; CONTRIBUTION FROM ',i3,
4279 *
'-TH HARD POMERON')
4298 IF(si.LT.4.
d0*(qt0+amj0))
THEN 4309 ELSEIF(iqq.EQ.2)
THEN 4311 ELSEIF(iqq.EQ.3)
THEN 4321 4 ibh(ih1)=ibh(ih1+1)
4330 IF(yi.GT.17.
d0)yi=17.
d0 4335 IF(debug.GE.3)
WRITE (moniou,213)
4336 213
FORMAT(2
x,
'PSSHARE: GBH-INI - END')
4350 lha0(
i)=lha(
i)-lva(
i)
4354 lhb0(
i)=lhb(
i)-lvb(
i)
4363 IF(
qsran(b10).GT.(1.
d0-xw)**ahv(icz))
GOTO 9
4364 IF(debug.GE.3)
WRITE (moniou,214)
i,xw
4365 214
FORMAT(2
x,
'PSSHARE: ',i2,
'-TH PROJ. NUCLEON (HADRON); LIGHT',
4366 *
' CONE MOMENTUM SHARE XW=',
e10.3)
4370 wp(
i)=wp(
i)*(1.
d0-xw)
4379 IF(
qsran(b10).GT.(1.
d0-xw)**ahv(2))
GOTO 11
4380 IF(debug.GE.3)
WRITE (moniou,215)
i,xw
4381 215
FORMAT(2
x,
'PSSHARE: ',i2,
'-TH TARGET NUCLEON (HADRON); LIGHT',
4382 *
' CONE MOMENTUM SHARE XW=',
e10.3)
4386 wm(
i)=wm(
i)*(1.
d0-xw)
4406 IF((iqq-3)*(iqq-1).EQ.0)
THEN 4417 bpi=1.
d0/(1.
d0+ahl(icz)+
4418 * (1.
d0+delh)*lha0(
i))
4421 15 xw=1.-
qsran(b10)**bpi
4423 IF(
qsran(b10).GT.xw**delh)
GOTO 15
4427 wp(
i)=wp(
i)*(1.
d0-xw)
4430 IF((iqq-3)*(iqq-2).EQ.0)
THEN 4437 bpi=1.
d0/(1.
d0+ahl(2)+(1.
d0+delh)
4441 16 xw=1.-
qsran(b10)**bpi
4442 IF(
qsran(b10).GT.xw**delh)
GOTO 16
4446 wm(
j)=wm(
j)*(1.
d0-xw)
4451 IF(sw.LT.4.
d0*(qt0+amj0))
THEN 4466 ELSEIF(iqq.EQ.2)
THEN 4468 ELSEIF(iqq.EQ.3)
THEN 4478 17 ibh(ih1)=ibh(ih1+1)
4489 IF(debug.GE.3)
WRITE (moniou,216)ih,wha(ih),whb(ih),wp(
i),wm(
j)
4490 216
FORMAT(2
x,
'PSSHARE: ',i3,
'-TH SEMIHARD BLOCK; LIGHT',
4491 *
' CONE MOMENTA SHARES:',2
e10.3/
4492 * 4
x,
'REMAINED LIGHT CONE MOMENTA:',2
e10.3)
4509 IF(debug.GE.2)
WRITE (moniou,217)1.
d0-gbh,nhp
4510 217
FORMAT(2
x,
'PSSHARE: REJECTION PROBABILITY:',
e10.3,
4511 * 2
x,
'NUMBER OF SEMIHARD BLOCKS:',i3)
4512 IF(
qsran(b10).GT.gbh)
THEN 4516 IF(debug.GE.2)
WRITE (moniou,218)
4517 218
FORMAT(2
x,
'PSSHARE: MORE THAN 30 REJECTIONS - HARD POMERON',
4518 *
' NUMBER IS PUT DOWN')
4524 DO 19 ihp=nhp-lnh+1,nhp
4531 ELSEIF(iqq.EQ.2)
THEN 4533 ELSEIF(iqq.EQ.3)
THEN 4539 19 lhb(jih)=lhb(jih)-1
4573 IF(debug.GE.2)
WRITE (moniou,219)ih,iqq,wpi,wmi,wp(
i),wm(
j)
4574 219
FORMAT(2
x,
'PSSHARE: ',i3,
4575 *
'-TH HARD BLOCK; TYPE OF THE INTERACTION:',i1/
4576 * 4
x,
'INITIAL LIGHT CONE MOMENTA:',2
e10.3/
4577 * 4
x,
'REMAINED LIGHT CONE MOMENTA:',2
e10.3)
4581 CALL pshot(wpi,wmi,
z,ipc,epc,izp(
i),izt(
j),icz,iqq)
4582 IF(iqq.EQ.1.OR.iqq.EQ.3)
THEN 4583 IF((iabs(izp(
i)).GT.5.OR.iabs(izp(
i)).EQ.3).AND.
4584 * izp(
i).GT.0.OR.iabs(izp(
i)).NE.3.AND.
4585 * iabs(izp(
i)).LE.5.AND.izp(
i).LT.0)
THEN 4592 330 ela(l,
i)=epc(l+4*(jq-1),1)
4594 IF(iqq.EQ.2.OR.iqq.EQ.3)
THEN 4595 IF((iabs(izt(
j)).GT.5.OR.iabs(izt(
j)).EQ.3).AND.
4596 * izt(
j).GT.0.OR.iabs(izt(
j)).NE.3.AND.
4597 * iabs(izt(
j)).LE.5.AND.izt(
j).LT.0)
THEN 4604 331 elb(l,
j)=epc(l+4*(jq-1),2)
4606 IF(iqq.EQ.3.AND.ila(
i)+ilb(
j).EQ.0)nias=
j 4619 IF(lqa(
i)+lha0(
i).EQ.0.AND.lqb(
j)+lhb0(
j).EQ.0)
THEN 4620 IF(lva(
i).EQ.0.AND.lvb(
j).EQ.0)
THEN 4622 ELSEIF(lva(
i).EQ.0)
THEN 4632 IF(iabs(ipj1).EQ.3)ipj1=ipj1*4/3
4633 CALL psjdef(izt(
j),ipj1,ep,ep1,jfl)
4634 IF(jfl.EQ.0)
GOTO 100
4636 ELSEIF(lvb(
j).EQ.0)
THEN 4646 IF(iabs(ipj1).EQ.3)ipj1=ipj1*4/3
4647 CALL psjdef(izp(
i),ipj1,ep,ep1,jfl)
4648 IF(jfl.EQ.0)
GOTO 100
4659 IF(iabs(ipj1).EQ.3)ipj1=ipj1*4/3
4660 CALL psjdef(izp(
i),ipj1,ep,ep1,jfl)
4661 IF(jfl.EQ.0)
GOTO 100
4671 IF(iabs(ipj1).EQ.3)ipj1=ipj1*4/3
4672 CALL psjdef(izt(
j),ipj1,ep,ep1,jfl)
4673 IF(jfl.EQ.0)
GOTO 100
4676 ELSEIF(lqa(
i)+lha0(
i).EQ.0)
THEN 4678 CALL xxdpr(wp(
i),wm(
j),izp(
i),izt(
j),lqb(
j)+lhb0(
j))
4688 IF(iabs(ipj1).EQ.3)ipj1=ipj1*4/3
4689 CALL psjdef(izp(
i),ipj1,ep,ep1,jfl)
4690 IF(jfl.EQ.0)
GOTO 100
4693 ELSEIF(lqb(
j)+lhb0(
j).EQ.0)
THEN 4695 CALL xxdtg(wp(
i),wm(
j),izp(
i),izt(
j),lqa(
i)+lha0(
i))
4705 IF(iabs(ipj1).EQ.3)ipj1=ipj1*4/3
4706 CALL psjdef(izt(
j),ipj1,ep,ep1,jfl)
4707 IF(jfl.EQ.0)
GOTO 100
4748 IF(debug.GE.3)
WRITE (moniou,222)
is,
i,
j,np
4749 222
FORMAT(2
x,
'PSSHARE: ',i3,
'-TH SOFT POMERON BLOCK IS',
4750 *
' CONNECTED TO ',i2,
4751 *
'-TH PROJECTILE NUCLEON'/4
x,
'(HADRON) AND ',i2,
4752 *
'-TH TARGET NUCLEON'/
4753 * 4
x,
'NUMBER OF CUT SOFT POMERONS IN THE BLOCK:',i2)
4762 IF(lq1.EQ.1.AND.wpn.EQ.wp0.AND.
qsran(b10).LT.wppp
4763 * .AND.lvb(
j).EQ.0)
THEN 4769 IF(debug.GE.3)
WRITE (moniou,223)yw
4770 223
FORMAT(2
x,
'PSSHARE: TRIPLE POMERON CONTRIBUTION YW=',
e10.3)
4782 bpi=1.
d0/(1.
d0+ahl(icz)+(1.
d0+del)*lq1)
4783 23 xpw=1.-
qsran(b10)**bpi
4785 IF(
qsran(b10).GT.xpw**del)
GOTO 23
4791 bpi=1.
d0/(1.
d0+ahl(2)+(1.
d0+del)*lq2)
4792 24 xmw=1.-
qsran(b10)**bpi
4794 IF(
qsran(b10).GT.xmw**del)
GOTO 24
4799 IF(jpp.EQ.1.AND.xpw*xmw*wpn*wmn.LT.2.72
d0)
THEN 4814 IF(lq1.EQ.0.AND.lva(
i).EQ.0)
THEN 4815 CALL ixxdef(izp(
i),ic11,ic12,icz)
4820 IF(lq2.EQ.0.AND.lvb(
j).EQ.0)
THEN 4821 CALL ixxdef(izt(
j),ic21,ic22,2)
4831 IF(debug.GE.3)
WRITE (moniou,224)ip,wpi,wmi
4832 224
FORMAT(2
x,
'PSSHARE: ',i2,
'-TH SOFT POMERON IN THE BLOCK'/
4833 * 4
x,
'LIGHT CONE MOMENTA FOR THE POMERON:',2
e10.3)
4834 CALL xxstr(wpi,wmi,wpn,wmn,ic11,ic12,ic22,ic21)
4842 IF(wpn.LT.0.
d0.OR.wmn.LT.0.
d0.OR.
4843 * sw.LT.(am(icz)+am(2))**2)
THEN 4846 *
write (*,*)
'difr,i,j,WPn,WMn,sw,lq1,lq2',
4847 *
i,
j,wpn,wmn,sw,lq1,lq2
4854 CALL xxdtg(wpn,wmn,izp(
i),izt(
j),0)
4858 ep3(1)=.5
d0*(wp1+wm1)
4859 ep3(2)=.5
d0*(wp1-wm1)
4870 IF(debug.GE.3)
WRITE (moniou,225)
4871 225
FORMAT(2
x,
'PSSHARE: TRIPLE POMERON CONRITRIBUTION WITH 3 CUT',
4876 wmm(2)=wmm(1)*
qsran(b10)
4877 wmm(1)=wmm(1)-wmm(2)
4882 bpi=(1.
d0+del)*lq1+1.
d0+ahl(icz)
4884 25 xpw=1.-
qsran(b10)**bpi
4885 IF(
qsran(b10).GT.xpw**del)
GOTO 25
4889 26
CALL xxstr(wpp(l),wmm(l),wpn,wmn,0,0,0,0)
4891 IF(wpn.LT.0.
d0.OR.wmn.LT.0.
d0.OR.
4892 * sw.LT.(am(icz)+am(2))**2)
THEN 4912 IF(lq1.EQ.0.AND.lq2.EQ.0)
THEN 4913 IF(lva(
i).EQ.0.AND.lvb(
j).EQ.0)
THEN 4915 ELSEIF(lva(
i).EQ.0)
THEN 4916 CALL xxdpr(wpn,wmn,izp(
i),izt(
j),1)
4925 IF(iabs(ipj1).EQ.3)ipj1=ipj1*4/3
4926 CALL psjdef(izt(
j),ipj1,ep,ep1,jfl)
4927 IF(jfl.EQ.0)
GOTO 100
4929 ELSEIF(lvb(
j).EQ.0)
THEN 4930 CALL xxdtg(wpn,wmn,izp(
i),izt(
j),1)
4939 IF(iabs(ipj1).EQ.3)ipj1=ipj1*4/3
4940 CALL psjdef(izp(
i),ipj1,ep,ep1,jfl)
4941 IF(jfl.EQ.0)
GOTO 100
4952 IF(iabs(ipj1).EQ.3)ipj1=ipj1*4/3
4953 CALL psjdef(izp(
i),ipj1,ep,ep1,jfl)
4954 IF(jfl.EQ.0)
GOTO 100
4964 IF(iabs(ipj1).EQ.3)ipj1=ipj1*4/3
4965 CALL psjdef(izt(
j),ipj1,ep,ep1,jfl)
4966 IF(jfl.EQ.0)
GOTO 100
4970 ELSEIF(lq1.EQ.0)
THEN 4972 CALL xxdpr(wpn,wmn,izp(
i),izt(
j),lq2)
4982 IF(iabs(ipj1).EQ.3)ipj1=ipj1*4/3
4983 CALL psjdef(izp(
i),ipj1,ep,ep1,jfl)
4984 IF(jfl.EQ.0)
GOTO 100
4988 ELSEIF(lq2.EQ.0)
THEN 4990 CALL xxdtg(wpn,wmn,izp(
i),izt(
j),lq1)
5000 IF(iabs(ipj1).EQ.3)ipj1=ipj1*4/3
5001 CALL psjdef(izt(
j),ipj1,ep,ep1,jfl)
5002 IF(jfl.EQ.0)
GOTO 100
5019 IF(ia(1).EQ.1.AND.lva(1).NE.0.AND.ila(1).EQ.0)
THEN 5024 ep1(1)=.5
d0*wm(nias)
5028 CALL psjdef(izp(1),izt(nias),ep,ep1,jfl)
5029 IF(jfl.EQ.0)
GOTO 100
5034 IF(debug.GE.3)
WRITE (moniou,227)
5035 227
FORMAT(2
x,
'PSSHARE - END')
5044 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
5046 dimension ey(3),ep(4)
5047 COMMON /area43/ moniou
5048 COMMON /debug/ debug
5051 IF(debug.GE.2)
WRITE (moniou,201)ep,ey
5052 201
FORMAT(2
x,
'PSTRANS - LORENTZ BOOST FOR 4-VECTOR'/4
x,
'EP=',
5053 * 2
x,4(
e10.3,1
x)/4
x,
'BOOST PARAMETERS EY=',3
e10.3)
5056 IF(ey(4-
i).NE.1.
d0)
THEN 5057 wp=(ep(1)+ep(5-
i))/ey(4-
i)
5058 wm=(ep(1)-ep(5-
i))*ey(4-
i)
5060 ep(5-
i)=.5
d0*(wp-wm)
5063 IF(debug.GE.3)
WRITE (moniou,202)ep
5064 202
FORMAT(2
x,
'PSTRANS: TRANSFORMED 4-VECTOR EP=',
5074 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
5076 dimension ey(3),ep(4)
5077 COMMON /area43/ moniou
5078 COMMON /debug/ debug
5081 IF(debug.GE.2)
WRITE (moniou,201)ep,ey
5082 201
FORMAT(2
x,
'PSTRANS1 - LORENTZ BOOST FOR 4-VECTOR'/4
x,
'EP=',
5083 * 2
x,4(
e10.3,1
x)/4
x,
'BOOST PARAMETERS EY=',3
e10.3)
5086 IF(ey(
i).NE.1.
d0)
THEN 5087 wp=(ep(1)+ep(
i+1))*ey(
i)
5088 wm=(ep(1)-ep(
i+1))/ey(
i)
5090 ep(
i+1)=.5
d0*(wp-wm)
5093 IF(debug.GE.3)
WRITE (moniou,202)ep
5094 202
FORMAT(2
x,
'PSTRANS1: TRANSFORMED 4-VECTOR EP=',
5105 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
5108 COMMON /area33/ fsud(10,2)
5109 COMMON /area43/ moniou
5110 COMMON /debug/ debug
5113 IF(debug.GE.2)
WRITE (moniou,201)
j,qlmax
5114 201
FORMAT(2
x,
'PSUDINT - SPACELIKE FORM FACTOR INTERPOLATION:'/
5115 * 4
x,
'PARTON TYPE J=',
5116 * i1,2
x,
'MOMENTUM LOGARITHM QLMAX=',
e10.3)
5125 wk(3)=wk(2)*(wk(2)-1.
d0)*.5
d0 5126 wk(1)=1.
d0-wk(2)+wk(3)
5127 wk(2)=wk(2)-2.
d0*wk(3)
5135 IF(debug.GE.3)
WRITE (moniou,202)
psudint 5136 202
FORMAT(2
x,
'PSUDINT=',
e10.3)
5146 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
5148 COMMON /area6/
pi,bm,am
5149 COMMON /area18/ alm,qt0,qlog,qll,aqt0,qtf,bet,amj0
5150 COMMON /area43/ moniou
5151 COMMON /debug/ debug
5154 IF(debug.GE.2)
WRITE (moniou,201)
j,q
5155 201
FORMAT(2
x,
'PSUDS - SPACELIKE FORM FACTOR: PARTON TYPE J=',
5156 * i1,2
x,
'MOMENTUM Q=',
e10.3)
5159 psuds=(qlm*dlog(qlm/qlog)-dlog(q/qt0))/9.
d0 5171 IF(debug.GE.3)
WRITE (moniou,202)
psuds 5172 202
FORMAT(2
x,
'PSUDS=',
e10.3)
5177 FUNCTION psudt(QMAX,J)
5182 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
5184 COMMON /area18/ alm,qt0,qlog,qll,aqt0,qtf,bet,amj0
5185 common/ar3/
x1(7),a1(7)
5186 COMMON /area43/ moniou
5187 COMMON /debug/ debug
5190 IF(debug.GE.2)
WRITE (moniou,201)
j,qmax
5191 201
FORMAT(2
x,
'PSUDT - TIMELIKE FORM FACTOR: PARTON TYPE J=',
5192 * i1,2
x,
'MOMENTUM QMAX=',
e10.3)
5194 qlmax=dlog(dlog(qmax/16.
d0/alm))
5195 qfl=dlog(dlog(qtf/alm))
5201 qtl=.5
d0*(qlmax+qfl+(2*
m-3)*
x1(
i)*(qlmax-qfl))
5202 qt=alm*exp(exp(qtl))
5203 IF(qt.GE.qmax/16.
d0)qt=qmax/16.0001
d0 5204 zmin=.5
d0-dsqrt((.25
d0-dsqrt(qt/qmax)))
5216 IF(debug.GE.3)
WRITE (moniou,202)
psudt 5217 202
FORMAT(2
x,
'PSUDT=',
e10.3)
5222 FUNCTION psv(X,Y,XB,IB)
5226 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
5229 dimension xb(64,3),fhard(3)
5230 COMMON /area43/ moniou
5231 COMMON /debug/ debug
5234 IF(debug.GE.2)
WRITE (moniou,201)
x,
y,ib
5235 201
FORMAT(2
x,
'PSV - EIKONAL FACTOR: NUCLEON COORDINATES X=',
5236 *
e10.3,2
x,
'Y=',
e10.3/4
x,
'NUMBER OF ACTIVE TARGET NUCLEONS IB=' 5242 dv=dv+
psfaz(
z,fsoft,fhard,fshard)+fhard(1)+fhard(2)+fhard(3)
5244 psv=(1.
d0-exp(-dv))**2
5252 IF(debug.GE.3)
WRITE (moniou,202)
psv 5253 202
FORMAT(2
x,
'PSV=',
e10.3)
5258 SUBROUTINE psvdef(ICH,IC1,ICZ)
5262 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
5265 COMMON /area43/ moniou
5266 COMMON /debug/ debug
5270 IF(debug.GE.2)
WRITE (moniou,201)ich,icz
5271 201
FORMAT(2
x,
'PSVDEF: HADRON TYPE ICH=',i2,
' AUXILLIARY TYPE ICZ=' 5276 ic1=ich*(1-3*int(.5+
qsran(b10)))
5278 ELSEIF(icz.EQ.2)
THEN 5279 IF(
qsran(b10).GT..33333
d0.OR.ich.LT.0)
THEN 5286 ELSEIF(icz.EQ.3)
THEN 5289 ELSEIF(icz.EQ.4)
THEN 5293 IF(debug.GE.3)
WRITE (moniou,202)ic1,ich
5294 202
FORMAT(2
x,
'PSVDEF-END: QUARK FLAVOR IC1=',i2,
5295 *
'DIQUARK TYPE ICH=',i2)
5306 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
5309 COMMON /area18/ alm,qt0,qlog,qll,aqt0,qtf,bet,amj0
5310 COMMON /area43/ moniou
5311 COMMON /debug/ debug
5316 IF(debug.GE.2)
WRITE (moniou,201)qq,
j 5317 201
FORMAT(2
x,
'PSZSIM - Z-SHARE SIMULATION: QQ=',
e10.3,2
x,
'J=',i1)
5318 zmin=.5
d0-dsqrt(.25
d0-dsqrt(qtf/qq))
5332 gb=gb/dlog(qt/alm)*qlf
5333 IF(debug.GE.3)
WRITE (moniou,203)qt,gb
5334 203
FORMAT(2
x,
'PSZSIM: QT=',
e10.3,2
x,
'GB=',
e10.3)
5335 IF(
qsran(b10).GT.gb)
GOTO 1
5336 IF(debug.GE.3)
WRITE (moniou,202)
pszsim 5337 202
FORMAT(2
x,
'PSZSIM=',
e10.3)
5342 SUBROUTINE ixxdef(ICH,IC1,IC2,ICZ)
5346 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
5349 COMMON /area43/ moniou
5350 COMMON /debug/ debug
5354 IF(debug.GE.2)
WRITE (moniou,201)ich,icz
5355 201
FORMAT(2
x,
'IXXDEF: HADRON TYPE ICH=',i2,
' AUXILLIARY TYPE ICZ=' 5359 ic1=ich*(1-3*int(.5+
qsran(b10)))
5361 ic2=-ic1*iabs(ich1)-(ich+ic1)*iabs(ich-ich1)
5363 ELSEIF(icz.EQ.2)
THEN 5365 ic1=int(1.3333+
qsran(b10))
5367 ich1=(2-ic1)*int(
qsran(b10)+.5)+2
5370 ic2=(3-ich1)*(2-ic1)-2
5372 IF(iabs(ich).EQ.3)
THEN 5383 ELSEIF(icz.EQ.3)
THEN 5387 ELSEIF(icz.EQ.4)
THEN 5391 ELSEIF(icz.EQ.5)
THEN 5398 IF(debug.GE.3)
WRITE (moniou,202)ic1,ic2,ich
5399 202
FORMAT(2
x,
'IXXDEF-END: PARTON FLAVORS IC1=',i2,
' IC2=',i2,
5400 *
'NEW HADRON TYPE ICH=',i2)
5411 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
5413 COMMON /area43/ moniou
5414 COMMON /debug/ debug
5417 IF(debug.GE.2)
WRITE (moniou,201)ns-1,aw,
g 5418 201
FORMAT(2
x,
'IXXSON - POISSON DITR.: AVERAGE AW=',
e10.3,
5419 *
' MAXIMAL VALUE NS=',i2,
' RANDOM NUMBER G=',
e10.3)
5428 IF(debug.GE.3)
WRITE (moniou,202)
ixxson 5429 202
FORMAT(2
x,
'IXXSON=',i2)
5434 SUBROUTINE xxaini(E0N,ICP0,IAP,IAT)
5437 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
5440 dimension wk(3),wa(3)
5442 COMMON /area1/ ia(2),icz,icp
5443 COMMON /area2/ s,y0,wp0,wm0
5444 COMMON /area4/ ey0(3)
5445 COMMON /area5/ rd(2),cr1(2),cr2(2),cr3(2)
5446 COMMON /area6/
pi,bm,am
5448 COMMON /area10/ stmass,am0,amn,amk,amc,amlamc,amlam,ameta
5449 COMMON /area15/ fp(5),rq(5),cd(5)
5450 COMMON /area17/ del,rs,rs0,fs,alfp,rr,sh,delh
5451 COMMON /area22/ sjv,fjs(5,3)
5452 COMMON /area35/ sjv0(10,5),fjs0(10,5,15)
5453 COMMON /area43/ moniou
5455 COMMON /area44/ gz(10,5,4),gzp(10,5,4)
5456 COMMON /area45/ gdt,gdp
5458 COMMON /debug/ debug
5461 IF(debug.GE.1)
WRITE (moniou,201)icp0,iap,iat,e0n
5462 201
FORMAT(2
x,
'XXAINI - MINIINITIALIZATION: PARTICLE TYPE ICP0=',
5463 * i1,2
x,
'PROJECTILE MASS NUMBER IAP=',i2/4
x,
5464 *
'TARGET MASS NUMBER IAT=',i2,
' INTERACTION ENERGY E0N=',
e10.3)
5470 IF(iabs(icp).LT.6)
THEN 5488 fs=fp(icz)*exp(y0*del)/rs*cd(icz)
5490 rp1=rs*4.
d0*.0391
d0/am**2
5494 ey0(1)=dsqrt(amn/e0n/2.
d0)
5500 rd(
i)=0.7
d0*float(ia(
i))**.446/am
5503 cr3(
i)=3.
d0/rd(
i)+6.
d0/rd(
i)**2
5504 IF(ia(
i).LT.10.AND.ia(
i).NE.1)
THEN 5506 rd(
i)=.9
d0*float(ia(
i))**.3333/am
5507 IF(ia(
i).EQ.2)rd(
i)=3.16
d0 5509 rd(
i)=rd(
i)*dsqrt(2.
d0*ia(
i)/(ia(
i)-1.))
5529 IF(ye.LT.1.
d0)ye=1.
d0 5535 wk(3)=wk(2)*(wk(2)-1.
d0)*.5
d0 5536 wk(1)=1.
d0-wk(2)+wk(3)
5537 wk(2)=wk(2)-2.
d0*wk(3)
5539 sjv=sjv0(je,icz)*wk(1)+sjv0(je+1,icz)*wk(2)+sjv0(je+2,icz)*wk(3)
5544 2 fjs(
i,
m)=fjs0(je,
i,m1)*wk(1)+fjs0(je+1,
i,m1)*wk(2)+
5545 *fjs0(je+2,
i,m1)*wk(3)
5551 ya=dlog(ya)/1.38629
d0+1.
d0 5554 wa(3)=wa(2)*(wa(2)-1.
d0)*.5
d0 5555 wa(1)=1.
d0-wa(2)+wa(3)
5556 wa(2)=wa(2)-2.
d0*wa(3)
5559 gdp=gdp+gzp(je+
i-1,icz,ja+
m-1)*wk(
i)*wa(
m)
5560 3 gdt=gdt+gz(je+
i-1,icz,ja+
m-1)*wk(
i)*wa(
m)
5565 IF(debug.GE.3)
WRITE (moniou,202)
5566 202
FORMAT(2
x,
'XXAINI - END')
5574 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
5577 COMMON /area3/ rmin,emax,eev
5578 COMMON /area6/
pi,bm,am
5579 COMMON /area8/ wwm,be(4),dc(5),deta,almpt
5580 COMMON /area10/ stmass,am0,amn,amk,amc,amlamc,amlam,ameta
5582 COMMON /area20/ wppp
5583 COMMON /area21/ dmmin(5)
5584 COMMON /area28/ arr(4)
5585 COMMON /area40/ jdifr
5586 COMMON /area42/ tyq(15)
5587 COMMON /area43/ moniou
5588 COMMON /debug/ debug
5591 IF(debug.GE.1)
WRITE (moniou,201)
5592 201
FORMAT(2
x,
'XXASET - HADRONIZATION PARAMETERS SETTING')
5686 IF(debug.GE.3)
WRITE (moniou,202)
5687 202
FORMAT(2
x,
'XXASET - END')
5692 SUBROUTINE xxddfr(WP0,WM0,ICP,ICT)
5695 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
5697 dimension ep3(4),ep1(4),ep2(4),ey(3)
5698 COMMON /area1/ ia(2),icz,icp0
5699 COMMON /area2/ s,y0,wp00,wm00
5700 COMMON /area8/ wwm,be(4),dc(5),deta,almpt
5701 COMMON /area10/ stmass,am(7)
5703 COMMON /area21/ dmmin(5)
5704 COMMON /area43/ moniou
5705 COMMON /debug/ debug
5709 IF(debug.GE.2)
WRITE (moniou,201)icp,ict,wp0,wm0
5710 201
FORMAT(2
x,
'XXDDFR - LEADING CLUSTERS HADRONIZATION:' 5711 * /4
x,
'CLUSTER TYPES ICP=',i2,2
x,
5712 *
'ICT=',i2/4
x,
'AVAILABLE LIGHT CONE MOMENTA: WP0=',
e10.3,
5718 IF(sd0.LT.0.
d0)sd0=0.
d0 5721 ddmax1=min(5.
d0,dsqrt(sd0)-ddmin2)
5723 IF(ddmax1.LT.ddmin1)
THEN 5726 IF(dsqrt(sd0).LT.am(icz)+am(2))
THEN 5727 IF(wp0.GT.0.
d0.AND.(am(icz)+am(2))**2/wp0.LT..5
d0*wm00)
THEN 5728 sd0=(am(icz)+am(2))**2
5731 IF(debug.GE.3)
WRITE (moniou,202)
5738 xw=
xxtwdec(sd0,am(icz)**2,am(2)**2)
5741 ep3(1)=.5
d0*(wp1+wm1)
5742 ep3(2)=.5
d0*(wp1-wm1)
5746 ep3(1)=.5
d0*(wp2+wm2)
5747 ep3(2)=.5
d0*(wp2-wm2)
5751 IF(debug.GE.3)
WRITE (moniou,202)
5755 dmass1=(ddmin1/(1.
d0-
qsran(b10)*(1.
d0-ddmin1/ddmax1)))**2
5756 ddmax2=min(5.
d0,dsqrt(sd0)-dsqrt(dmass1))
5757 dmass2=(ddmin2/(1.
d0-
qsran(b10)*(1.
d0-ddmin2/ddmax2)))**2
5759 wpd1=wp0*
xxtwdec(sd0,dmass1,dmass2)
5764 IF(icp.NE.0)
is=iabs(icp)/icp
5771 ptmax=
pslam(dmass1,amh1,amh2)
5772 IF(ptmax.LT.0.)ptmax=0.
5773 IF(ptmax.LT.be(4)**2)
THEN 5774 1 pti=ptmax*
qsran(b10)
5775 IF(
qsran(b10).GT.exp(-dsqrt(pti)/be(4)))
GOTO 1
5778 IF(pti.GT.ptmax)
GOTO 2
5785 ep3(1)=.5
d0*(wp1+wm1)
5786 ep3(2)=.5
d0*(wp1-wm1)
5791 CALL xxreg(ep3,ich1)
5795 ep3(1)=.5
d0*(wp1+wm1)
5796 ep3(2)=.5
d0*(wp1-wm1)
5799 CALL xxreg(ep3,ich2)
5805 ic1=icp*(1-3*int(.5
d0+
qsran(b10)))
5811 ELSEIF(icz.EQ.2)
THEN 5812 IF(
qsran(b10).GT..33333
d0)
THEN 5819 ELSEIF(icz.EQ.3)
THEN 5822 ELSEIF(icz.EQ.4)
THEN 5830 IF(
qsran(b10).GT..33333
d0)
THEN 5838 IF(debug.GE.3)
WRITE (moniou,202)
5839 202
FORMAT(2
x,
'XXDDFR - END')
5844 SUBROUTINE xxdec2(EP,EP1,EP2,WW,A,B)
5847 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
5849 dimension ep(4),ep1(4),ep2(4),ey(3)
5850 COMMON /area43/ moniou
5851 COMMON /debug/ debug
5856 IF(debug.GE.2)
WRITE (moniou,201)
5857 201
FORMAT(2
x,
'XXDEC2 - TWO PARTICLE DECAY')
5864 pt=pl*dsqrt(1.
d0-cosz**2)
5874 IF(debug.GE.3)
WRITE (moniou,202)
5875 202
FORMAT(2
x,
'XXDEC2 - END')
5880 SUBROUTINE xxdec3(EP,EP1,EP2,EP3,SWW,AM1,AM2,AM3)
5883 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
5885 dimension ep(4),ep1(4),ep2(4),ep3(4),ept(4),ey(3)
5887 COMMON /area43/ moniou
5888 COMMON /debug/ debug
5892 IF(debug.GE.2)
WRITE (moniou,201)
5893 201
FORMAT(2
x,
'XXDEC3 - THREE PARTICLE DECAY')
5898 emax=.25
d0*(sww+(am12-am23)/sww)**2
5899 gb0=dsqrt((emax-am12)/emax*(1.
d0-am23/s23max)
5900 * *(1.
d0-am32/s23max))
5901 1 p1=
qsran(b10)*(emax-am12)
5903 s23=sww**2+am12-2.
d0*e1*sww
5904 gb=dsqrt(p1*(1.
d0-am23/s23)*(1.
d0-am32/s23))/e1/gb0
5905 IF(
qsran(b10).GT.gb)
GOTO 1
5910 pt=p1*dsqrt(1.
d0-cosz**2)
5922 CALL xxdec2(ept,ep2,ep3,s23,am2**2,am3**2)
5923 IF(debug.GE.3)
WRITE (moniou,202)
5924 202
FORMAT(2
x,
'XXDEC3 - END')
5929 SUBROUTINE xxdpr(WP0,WM0,ICP,ICT,LQ2)
5933 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
5935 dimension ep3(4),ep1(4),ep2(4),ey(3)
5936 COMMON /area1/ ia(2),icz,icp0
5937 COMMON /area2/ s,y0,wp00,wm00
5938 COMMON /area8/ wwm,be(4),dc(5),deta,almpt
5939 COMMON /area10/ stmass,am(7)
5941 COMMON /area17/ del,rs,rs0,fs,alfp,rr,sh,delh
5942 COMMON /area21/ dmmin(5)
5943 COMMON /area43/ moniou
5944 COMMON /debug/ debug
5949 IF(debug.GE.2)
WRITE (moniou,201)icp,ict,wp0,wm0
5950 201
FORMAT(2
x,
'XXDPR - LEADING (PROJECTILE) CLUSTER HADRONIZATION:' 5951 * /4
x,
'CLUSTER TYPE ICP=',i2,2
x,
'TARGET TYPE ',
5952 *
'ICT=',i2/4
x,
'AVAILABLE LIGHT CONE MOMENTA: WP0=',
e10.3,
5958 IF(sd0.LT.0.
d0)sd0=0.
d0 5959 ddmax=min(5.
d0,dsqrt(sd0)-am(2))
5962 IF(ddmax.LT.ddmin)
THEN 5970 IF(am(icz)**2.GT.wpi*wm0)
THEN 5971 IF(wpi.GT.0.
d0.AND.am(icz)**2/wpi.LT..5
d0*wm00)
THEN 5984 ep3(1)=.5
d0*(wpi+wmi)
5985 ep3(2)=.5
d0*(wpi-wmi)
5987 IF(debug.GE.3)
WRITE (moniou,202)
5991 IF(dsqrt(sd0).LT.am(icz)+am(2))
THEN 5992 IF(wp0.GT.0.
d0.AND.(am(icz)+am(2))**2/wp0.LT..5
d0*wm00)
5994 sd0=(am(icz)+am(2))**2
5997 IF(debug.GE.3)
WRITE (moniou,202)
6001 xw=
xxtwdec(sd0,am(icz)**2,am(2)**2)
6004 ep3(1)=.5
d0*(wp1+wm1)
6005 ep3(2)=.5
d0*(wp1-wm1)
6009 ep3(1)=.5
d0*(wp2+wm2)
6010 ep3(2)=.5
d0*(wp2-wm2)
6015 IF(debug.GE.3)
WRITE (moniou,202)
6019 IF(icp.NE.0)
is=iabs(icp)/icp
6021 dmass=ddmin**2/(1.
d0-
qsran(b10)*(1.
d0-(ddmin/ddmax)))**2
6030 wpd=wp0*
xxtwdec(sd0,dmass,am(2)**2)
6034 ep3(1)=.5
d0*(wp2+wm2)
6035 ep3(2)=.5
d0*(wp2-wm2)
6040 ptmax=
pslam(sd0,dmass,am(2)**2)
6041 IF(ptmax.LT.0.)ptmax=0.
6042 pti=-1.
d0/rs*dlog(1.
d0-
qsran(b10)*(1.
d0-exp(-rs*ptmax)))
6046 wpd=wp0*
xxtwdec(sd0,amt1,amt2)
6051 CALL pscs(ccos,ssin)
6054 ep3(1)=.5
d0*(wp2+wm2)
6055 ep3(2)=.5
d0*(wp2-wm2)
6059 ep3(1)=.5
d0*(wpd+wmd)
6060 ep3(2)=.5
d0*(wpd-wmd)
6075 ptmax=
pslam(dmass,amh1,amh2)
6076 IF(ptmax.LT.0.)ptmax=0.
6077 IF(ptmax.LT.be(4)**2)
THEN 6078 1 pti=ptmax*
qsran(b10)
6079 IF(
qsran(b10).GT.exp(-dsqrt(pti)/be(4)))
GOTO 1
6082 IF(pti.GT.ptmax)
GOTO 2
6089 ep3(1)=.5
d0*(wp1+wm1)
6090 ep3(2)=.5
d0*(wp1-wm1)
6095 CALL xxreg(ep3,ich1)
6099 ep3(1)=.5
d0*(wp1+wm1)
6100 ep3(2)=.5
d0*(wp1-wm1)
6103 CALL xxreg(ep3,ich2)
6104 IF(debug.GE.3)
WRITE (moniou,202)
6110 ic1=icp*(1-3*int(.5
d0+
qsran(b10)))
6116 ELSEIF(icz.EQ.2)
THEN 6117 IF(
qsran(b10).GT..33333
d0)
THEN 6124 ELSEIF(icz.EQ.3)
THEN 6127 ELSEIF(icz.EQ.4)
THEN 6133 IF(debug.GE.3)
WRITE (moniou,202)
6134 202
FORMAT(2
x,
'XXDPR - END')
6139 SUBROUTINE xxdtg(WP0,WM0,ICP,ICT,LQ1)
6143 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
6145 dimension ep3(4),ey(3)
6146 COMMON /area1/ ia(2),icz,icp0
6147 COMMON /area2/ s,y0,wp00,wm00
6148 COMMON /area10/ stmass,am(7)
6150 COMMON /area17/ del,rs,rs0,fs,alfp,rr,sh,delh
6151 COMMON /area21/ dmmin(5)
6152 COMMON /area43/ moniou
6153 COMMON /debug/ debug
6158 IF(debug.GE.2)
WRITE (moniou,201)ict,icp,wp0,wm0
6159 201
FORMAT(2
x,
'XXDTG - LEADING (TARGET) CLUSTER HADRONIZATION:' 6160 * /4
x,
'CLUSTER TYPE ICT=',i2,2
x,
'PROJECTILE TYPE ',
6161 *
'ICP=',i2/4
x,
'AVAILABLE LIGHT CONE MOMENTA: WP0=',
e10.3,
6167 IF(sd0.LT.0.
d0)sd0=0.
d0 6169 ddmax=min(5.
d0,dsqrt(sd0)-am(icz))
6171 IF(ddmax.LT.ddmin)
THEN 6179 IF( wp0.LE.0.
d0.OR.am(2)**2.GT.wmi*wp0)
RETURN 6183 ep3(1)=.5
d0*(wpi+wmi)
6184 ep3(2)=.5
d0*(wpi-wmi)
6186 IF(debug.GE.3)
WRITE (moniou,202)
6190 IF(dsqrt(sd0).LT.am(icz)+am(2))
THEN 6191 IF(wp0.GT.0.
d0.AND.(am(icz)+am(2))**2/wp0.LT..5
d0*wm00)
6193 sd0=(am(icz)+am(2))**2
6196 IF(debug.GE.3)
WRITE (moniou,202)
6200 xw=
xxtwdec(sd0,am(icz)**2,am(2)**2)
6203 ep3(1)=.5
d0*(wp1+wm1)
6204 ep3(2)=.5
d0*(wp1-wm1)
6208 ep3(1)=.5
d0*(wp2+wm2)
6209 ep3(2)=.5
d0*(wp2-wm2)
6214 IF(debug.GE.3)
WRITE (moniou,202)
6218 dmass=(ddmin/(1.
d0-
qsran(b10)*(1.
d0-ddmin/ddmax)))**2
6225 ptmax=
pslam(sd0,dmass,am(icz)**2)
6226 IF(ptmax.LT.0.)ptmax=0.
6227 pti=-1.
d0/rs*dlog(1.
d0-
qsran(b10)*(1.
d0-exp(-rs*ptmax)))
6231 wmd=wm0*
xxtwdec(sd0,amt1,amt2)
6236 CALL pscs(ccos,ssin)
6239 ep3(1)=.5
d0*(wp2+wm2)
6240 ep3(2)=.5
d0*(wp2-wm2)
6244 ep3(1)=.5
d0*(wpd+wmd)
6245 ep3(2)=.5
d0*(wpd-wmd)
6254 IF(
qsran(b10).GT..33333
d0)
THEN 6263 IF(debug.GE.3)
WRITE (moniou,202)
6264 202
FORMAT(2
x,
'XXDTG - END')
6269 SUBROUTINE xxfau(B,GZ)
6272 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
6274 dimension gz(3),gz0(2)
6275 COMMON /area1/ ia(2),icz,icp
6276 COMMON /area16/ cc(5)
6278 COMMON /area43/ moniou
6279 COMMON /debug/ debug
6282 IF(debug.GE.2)
WRITE (moniou,201)
6283 201
FORMAT(2
x,
'XXFAU - INTEGRANDS FOR HADRON-HADRON AND ' 6284 * ,
'HADRON-NUCLEUS CROSS-SECTIONS CALCULATION')
6288 1 gz0(l)=gz0(l)*cc(2)*anorm*.5
d0 6292 gz1=(1.
d0-gz0(1))**ab
6293 gz2=(1.
d0-gz0(2))**ab
6294 gz3=(1.
d0-cc(2)*gz0(2)-2.
d0*(1.
d0-cc(2))*gz0(1))**ab
6297 gz(1)=cc(icz)**2*(gz2-gz3)
6298 gz(2)=cc(icz)*(1.
d0-cc(icz))*(1.
d0+gz2-2.
d0*gz1)
6299 gz(3)=cc(icz)*(1.
d0-gz2)
6300 IF(debug.GE.3)
WRITE (moniou,202)
6301 202
FORMAT(2
x,
'XXFAU - END')
6306 SUBROUTINE xxfrag(SA,NA,RC)
6310 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
6314 COMMON /area13/ nsf,iaf(56)
6315 COMMON /area43/ moniou
6316 COMMON /debug/ debug
6319 IF(debug.GE.2)
WRITE (moniou,201)na
6320 201
FORMAT(2
x,
'XXFRAG-MULTIFRAGMENTATION: NUCLEUS MASS NUMBER: NA=' 6324 203
FORMAT(2
x,
'NUCLEONS COORDINATES:')
6325 204
FORMAT(2
x,3
e10.3)
6327 205
WRITE (moniou,204)(sa(
i,l),l=1,3)
6338 2 ri=ri+(sa(
j,
m)-sa(
i,
m))**2
6348 IF(
j.LT.ni.AND.na-ni.GT.0)
GOTO 1
6351 IF(debug.GE.3)
WRITE (moniou,206)nsf,iaf(nsf)
6352 206
FORMAT(2
x,
'XXFRAG: FRAGMENT N',i2,2
x,
'FRAGMENT MASS - ',i2)
6359 IF(debug.GE.3)
WRITE (moniou,206)nsf,iaf(nsf)
6361 IF(debug.GE.3)
WRITE (moniou,202)
6362 202
FORMAT(2
x,
'XXFRAG - END')
6372 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
6376 COMMON /area1/ ia(2),icz,icp
6377 COMMON /area3/ rmin,emax,eev
6381 COMMON /area13/ nsf,iaf(56)
6382 COMMON /area43/ moniou
6383 COMMON /debug/ debug
6387 IF(debug.GE.2)
WRITE (moniou,201)ns
6388 201
FORMAT(2
x,
'XXFRAGM: NUMBER OF SPECTATORS: NS=',i2)
6396 IF(debug.GE.3)
WRITE (moniou,205)
6397 205
FORMAT(2
x,
'XXFRAGM - SINGLE SPECTATOR')
6407 IF(debug.GE.3)
WRITE (moniou,203)eex
6408 203
FORMAT(2
x,
'XXFRAGM: EXCITATION ENERGY: EEX=',
e10.3)
6412 IF(eex/ns.GT.emax)
THEN 6423 IF(debug.GE.3)
WRITE (moniou,206)iaf(nsf)
6424 206
FORMAT(2
x,
'XXFRAGM - EVAPORATION: MASS NUMBER OF THE FRAGMENT:' 6444 IF(debug.GE.3)
WRITE (moniou,204)nf,nal
6445 204
FORMAT(2
x,
'XXFRAGM - EVAPORATION: NUMBER OF NUCLEONS NF=',i2,
6446 *
'NUMBER OF ALPHAS NAL=',i2)
6449 IF(debug.GE.3)
WRITE (moniou,202)
6450 202
FORMAT(2
x,
'XXFRAGM - END')
6455 SUBROUTINE xxfz(B,GZ)
6458 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
6460 dimension gz(2),fhard(3)
6461 COMMON /area1/ ia(2),icz,icp
6462 COMMON /area2/ s,y0,wp0,wm0
6464 COMMON /ar3/
x1(7),a1(7)
6465 COMMON /area43/ moniou
6466 COMMON /debug/ debug
6469 IF(debug.GE.2)
WRITE (moniou,201)
6470 201
FORMAT(2
x,
'XXFZ - HADRONIC CROSS-SECTIONS CALCULATION')
6481 s2=dsqrt(rp1*(1.
d0-dlog(
z)))
6489 vv1=exp(-
psfaz(zv1,fsoft,fhard,fshard)-fhard(1)
6490 * -fhard(2)-fhard(3))
6491 vv2=exp(-
psfaz(zv2,fsoft,fhard,fshard)-fhard(1)
6492 * -fhard(2)-fhard(3))
6504 2 gz(l)=gz(l)+ a1(i1)*(cg1*(1.
d0-vv1**l)+cg2*(1.
d0-vv2**l)/
z)
6505 IF(debug.GE.3)
WRITE (moniou,202)
6506 202
FORMAT(2
x,
'XXFZ - END')
6511 SUBROUTINE xxgau(GZ)
6515 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
6517 dimension gz(3),gz0(3)
6518 COMMON /area6/
pi,bm,am
6519 COMMON /ar3/
x1(7),a1(7)
6521 COMMON /area43/ moniou
6522 COMMON /debug/ debug
6525 IF(debug.GE.2)
WRITE (moniou,201)
6526 201
FORMAT(2
x,
'XXGAU - NUCLEAR CROSS-SECTIONS CALCULATION')
6536 2 gz(l)=gz(l)+gz0(l)*a1(
i)
6538 3 gz(l)=gz(l)*(bm*am)**2*
pi*.5
d0 6539 IF(debug.GE.3)
WRITE (moniou,202)
6540 202
FORMAT(2
x,
'XXGAU - END')
6549 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
6551 dimension gz(3),gz0(3)
6552 COMMON /area6/
pi,bm,am
6553 COMMON /ar5/ x5(2),a5(2)
6555 COMMON /area43/ moniou
6556 COMMON /debug/ debug
6559 IF(debug.GE.2)
WRITE (moniou,201)
6560 201
FORMAT(2
x,
'XXGAU1 - NUCLEAR CROSS-SECTIONS CALCULATION')
6566 1 gz(l)=gz(l)+gz0(l)*a5(
i)*exp(x5(
i))*
b*2.
d0*
pi*am*am
6567 IF(debug.GE.3)
WRITE (moniou,202)
6568 202
FORMAT(2
x,
'XXGAU1 - END')
6573 SUBROUTINE xxgener(WP0,WM0,EY0,S0X,C0X,S0,C0,IC1,IC2)
6582 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
6585 dimension wp(2),ic(2),ept(4),ep(4),ey(3),ey0(3)
6588 COMMON /area8/ wwm,bep,ben,bek,bec,dc(5),deta,almpt
6589 COMMON /area10/ stmass,am0,amn,amk,amc,amlamc,amlam,ameta
6591 COMMON /area19/ ahl(5)
6593 COMMON /area21/ dmmin(5)
6595 COMMON /area28/ arr(4)
6596 COMMON /area42/ tyq(15)
6597 COMMON /area43/ moniou
6598 COMMON /debug/ debug
6602 IF(debug.GE.2)
WRITE (moniou,201)tyq(8+ic1),tyq(8+ic2),
6603 * wp0,wm0,ey0,s0x,c0x,s0,c0
6604 201
FORMAT(2
x,
'XXGENER: PARTON FLAVORS AT THE ENDS OF THE STRING:',
6605 * 2
x,a2,2
x,a2/4
x,
'LIGHT CONE MOMENTA OF THE STRING: ',
e10.3,
6607 *
'S0X=',
e10.3,2
x,
'C0X=',
e10.3,2
x,
'S0=',
e10.3,2
x,
'C0=',
e10.3)
6610 ept(1)=.5
d0*(wp0+wm0)
6611 ept(2)=.5
d0*(wp0-wm0)
6622 WRITE (moniou,203)
j,tyq(iqt),ww
6623 203
FORMAT(2
x,
'XXGENER: CURRENT PARTON FLAVOR AT THE END ',i1,
6624 *
' OF THE STRING: ',a2/4
x,
' STRING MASS: ',
e10.3)
6634 ELSEIF(iaj.EQ.4)
THEN 6636 ELSEIF(iaj.EQ.5)
THEN 6642 IF(iab.LE.2.AND.sww.GT.restm+2.
d0*am0+wwm.OR.
6643 *iab.EQ.3.AND.sww.GT.restm+am0+amn+wwm.OR.
6644 *iab.EQ.4.AND.sww.GT.restm+am0+amk+wwm.OR.
6645 *iab.EQ.5.AND.sww.GT.restm+am0+amc+wwm)
THEN 6648 IF(sww.GT.restm+2.
d0*amc.AND.
qsran(b10).LT.dc(3))
THEN 6650 restm=(restm+amc)**2
6657 ELSEIF(sww.GT.restm+2.
d0*amn.AND.
qsran(b10).LT.dc(1))
THEN 6659 restm=(restm+amn)**2
6666 ELSEIF(sww.GT.restm+2.
d0*amk.AND.
qsran(b10).LT.dc(2))
THEN 6668 restm=(restm+amk)**2
6675 ELSEIF(sww.GT.restm+ameta+am0.AND.
qsran(b10).LT.deta)
THEN 6677 restm=(restm+am0)**2
6685 restm=(restm+am0)**2
6699 ELSEIF(iab.EQ.3)
THEN 6700 IF(sww.GT.restm+amc+amlamc.AND.
qsran(b10).LT.dc(5).AND.
6701 * iabs(ic(
j)).EQ.3)
THEN 6703 restm=(restm+amc)**2
6710 ELSEIF(sww.GT.restm+amk+amlam.AND.
qsran(b10).LT.dc(4).AND.
6711 * iabs(ic(
j)).EQ.3)
THEN 6713 restm=(restm+amk)**2
6717 blf=ahl(2)+arr(1)-arr(3)
6722 restm=(restm+am0)**2
6727 IF(iabs(ic(
j)).EQ.3)
THEN 6736 ELSEIF(iab.EQ.4)
THEN 6737 IF(sww.GT.restm+amn+amlam.AND.
qsran(b10).LT.dc(1))
THEN 6739 restm=(restm+amn)**2
6743 blf=ahl(2)+arr(1)-arr(3)
6748 restm=(restm+am0)**2
6757 ELSEIF(iab.EQ.5)
THEN 6758 IF(sww.GT.restm+amn+amlamc.AND.
qsran(b10).LT.dc(1))
THEN 6760 restm=(restm+amn)**2
6769 restm=(restm+am0)**2
6780 ptmax=
pslam(ww,restm,ami)
6781 IF(ptmax.LT.0.)ptmax=0.
6783 IF(ptmax.LT.bet**2)
THEN 6784 2 pti=ptmax*
qsran(b10)
6785 IF(
qsran(b10).GT.exp(-dsqrt(pti)/bet))
GOTO 2
6788 IF(pti.GT.ptmax)
GOTO 3
6800 4
z=1.-(z1+(z2-z1)*
qsran(b10))**(1./alf)
6801 IF(
qsran(b10).GT.(
z/zmax)**blf)
GOTO 4
6804 ep(1)=.5
d0*(wp(1)+wp(2))
6805 ep(2)=.5
d0*(wp(1)-wp(2))
6815 IF(ww.LT.restm)
GOTO 4
6820 IF(s0x.NE.0.
d0.OR.s0.NE.0.
d0)
THEN 6821 CALL psrotat(ep,s0x,c0x,s0,c0)
6824 IF(ey0(1)*ey0(2)*ey0(3).NE.1.
d0)
THEN 6833 IF(iab.LE.2.AND.iaj.LE.2)
THEN 6849 ELSEIF(iab.EQ.3.OR.iaj.EQ.3)
THEN 6852 IF(iabs(ic(
j)).EQ.3)
THEN 6854 IF(iabs(ic(3-
j)).EQ.3)
THEN 6858 ic(3-
j)=ic(3-
j)+4*
is 6861 ELSEIF(iaj.LT.3)
THEN 6867 ic(3-
j)=
is*(3-2*iaj)
6869 ELSEIF(iaj.EQ.4)
THEN 6872 ELSEIF(iaj.EQ.5)
THEN 6874 ic(3-
j)=-ic(
j)+10*
is 6881 ELSEIF(iaj.LT.3)
THEN 6882 ic(3-
j)=-ic(3-
j)-ic0
6883 ELSEIF(iaj.EQ.4)
THEN 6885 ELSEIF(iaj.EQ.5)
THEN 6890 IF(iabs(ic(3-
j)).EQ.3)
THEN 6900 ELSEIF(iab.EQ.4)
THEN 6904 ELSEIF(iab.EQ.5)
THEN 6907 ic(
j)=-ic(3-
j)+10*
is 6910 ic(3-
j)=ic(3-
j)-4*
is 6915 ELSEIF(iab.EQ.4)
THEN 6918 ELSEIF(iab.EQ.5)
THEN 6925 ELSEIF(iab.EQ.4.OR.iaj.EQ.4)
THEN 6933 ELSEIF(iaj.EQ.5)
THEN 6935 ic(3-
j)=-ic(
j)-12*
is 6937 ic0=ic(3-
j)+int(.6667
d0+
qsran(b10))*(-3*
is-2*ic(3-
j))
6947 ELSEIF(iab.EQ.5)
THEN 6950 ic(
j)=-ic(3-
j)+12*
is 6954 ELSEIF(iab.EQ.5.OR.iaj.EQ.5)
THEN 6963 ic0=ic(3-
j)+int(.6667
d0+
qsran(b10))*(-3*
is-2*ic(3-
j))
6975 ptmax=
pslam(ww,ami2,ami)
6976 IF(ptmax.LT.0.)ptmax=0.
6977 IF(ptmax.LT.bet**2)
THEN 6978 6 pti=ptmax*
qsran(b10)
6979 IF(
qsran(b10).GT.exp(-dsqrt(pti)/bet))
GOTO 6
6982 IF(pti.GT.ptmax)
GOTO 7
6991 ep(1)=.5
d0*(wp(1)+wp(2))
6992 ep(2)=.5
d0*(wp(1)-wp(2))
7005 IF(s0x.NE.0.
d0.OR.s0.NE.0.
d0)
THEN 7006 CALL psrotat(ep,s0x,c0x,s0,c0)
7007 CALL psrotat(ept,s0x,c0x,s0,c0)
7009 IF(ey0(1)*ey0(2)*ey0(3).NE.1.
d0)
THEN 7016 IF(debug.GE.3)
WRITE (moniou,202)
7017 202
FORMAT(2
x,
'XXGENER - END')
7029 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
7031 dimension ep(4),ep1(4),ey(3)
7032 COMMON /area10/ stmass,am(7)
7034 COMMON /area43/ moniou
7035 COMMON /debug/ debug
7036 COMMON /area46/ epjet(4,2,15000),ipjet(2,15000)
7037 COMMON /area47/ njtot
7040 IF(debug.GE.2)
WRITE (moniou,201)njtot
7041 201
FORMAT(2
x,
'XXJETSIM: TOTAL NUMBER OF JETS NJTOT=',i4)
7042 IF(njtot.EQ.0)
RETURN 7045 ep1(
i)=epjet(
i,1,nj)
7046 1 ep(
i)=ep1(
i)+epjet(
i,2,nj)
7047 pt3=dsqrt(ep1(3)**2+ep1(4)**2)
7048 pt4=dsqrt(epjet(3,2,nj)**2+epjet(4,2,nj)**2)
7058 2
CALL xxgener(sww,sww,ey,s0x,c0x,s0,c0,ipjet(1,nj),ipjet(2,nj))
7059 IF(debug.GE.3)
WRITE (moniou,202)
7060 202
FORMAT(2
x,
'XXJETSIM - END')
7065 SUBROUTINE xxreg(EP0,IC)
7070 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
7072 dimension ep(4),ep0(4)
7073 COMMON /area4/ ey0(3)
7074 COMMON /area10/ stmass,am0,amn,amk,amc,amlamc,amlam,ameta
7077 COMMON /area14/ esp(4,95000),ich(95000)
7078 COMMON /area43/ moniou
7079 COMMON /debug/ debug
7082 IF(debug.GE.2)
WRITE (moniou,201)ic,ep0
7083 201
FORMAT(2
x,
'XXREG: IC=',i2,2
x,
'C.M. 4-MOMENTUM:',2
x,4(
e10.3,1
x))
7084 pt=dsqrt(ep0(3)**2+ep0(4)**2)
7089 IF (nsh .GT. 95000)
THEN 7090 WRITE(moniou,*)
'XXREG: TOO MANY SECONDARY PARTICLES' 7091 WRITE(moniou,*)
'XXREG: NSH = ',nsh
7097 IF(debug.GE.3)
WRITE (moniou,202)ep
7098 202
FORMAT(2
x,
'XXREG: LAB. 4-MOMENTUM:',2
x,4(
e10.3,1
x))
7104 IF(debug.GE.3)
WRITE (moniou,203)
7105 203
FORMAT(2
x,
'XXREG - END')
7113 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
7115 COMMON /ar8/
x2(4),a2
7116 COMMON /area43/ moniou
7117 COMMON /debug/ debug
7120 IF(debug.GE.2)
WRITE (moniou,201)
b 7121 201
FORMAT(2
x,
'XXROT - AXIAL ANGLE INTEGRATION OF THE ',
7122 *
'NUCLEAR PROFILE FUNCTION'/4
x,
7123 *
'IMPACT PARAMETER B=',
e10.3,2
x,
'NUCLEON COORDINATE S=',
e10.3)
7127 sb1=
b**2+s**2-2.*
b*s*(2.*
x2(
i)-1.)
7128 sb2=
b**2+s**2-2.*
b*s*(1.-2.*
x2(
i))
7131 IF(debug.GE.3)
WRITE (moniou,202)
xxrot 7132 202
FORMAT(2
x,
'XXROT=',
e10.3)
7137 SUBROUTINE xxstr(WPI0,WMI0,WP0,WM0,IC10,IC120,IC210,IC20)
7143 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
7146 COMMON /area6/
pi,bm,ammm
7147 COMMON /area10/ stmass,am(7)
7149 COMMON /area43/ moniou
7150 COMMON /debug/ debug
7154 IF(debug.GE.2)
WRITE (moniou,201)wpi0,wmi0,wp0,wm0
7155 201
FORMAT(2
x,
'XXSTR: WPI0=',
e10.3,2
x,
'WMI0=',
e10.3,2
x,
7156 *
'WP0=',
e10.3,2
x,
'WM0=',
e10.3)
7166 ic1=int(1.5+
qsran(b10))
7168 ELSEIF(ic10.GT.0)
THEN 7176 ic2=int(1.5+
qsran(b10))
7178 ELSEIF(ic20.gt.0)
THEN 7197 IF(sm1.GT.stmass.AND.sm2.GT.stmass)
THEN 7201 ELSEIF(sm1.GT.stmass)
THEN 7203 ELSEIF(sm2.GT.stmass)
THEN 7209 IF(debug.GE.3)
WRITE (moniou,202)wp0,wm0
7210 202
FORMAT(2
x,
'XXSTR - RETURNED LIGHT CONE MOMENTA:',
7219 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
7221 COMMON /area6/
pi,bm,am
7223 COMMON /ar5/ x5(2),a5(2)
7224 COMMON /ar9/ x9(3),a9(3)
7225 COMMON /area43/ moniou
7226 COMMON /debug/ debug
7229 IF(debug.GE.2)
WRITE (moniou,201)
b 7230 201
FORMAT(2
x,
'XXT - NUCLEAR PROFILE FUNCTION VALUE AT IMPACT',
7231 *
' PARAMETER SQUARED B=',
e10.3)
7241 z1=zm*(1.+x9(
i))*0.5
7242 z2=zm*(1.-x9(
i))*0.5
7243 quq=dsqrt(
b+z1**2)-
r 7244 IF (quq.LT.85.)
xxt=
xxt+a9(
i)/(1.+exp(quq))
7245 quq=dsqrt(
b+z2**2)-
r 7246 IF (quq.LT.85.)
xxt=
xxt+a9(
i)/(1.+exp(quq))
7252 quq=dsqrt(
b+z1**2)-
r-x5(
i)
7253 IF (quq.LT.85.)dt=dt+a5(
i)/(exp(-x5(
i))+exp(quq))
7256 IF(debug.GE.3)
WRITE (moniou,202)
xxt 7257 202
FORMAT(2
x,
'XXT=',
e10.3)
7269 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
7271 COMMON /area43/ moniou
7272 COMMON /debug/ debug
7275 IF(debug.GE.2)
WRITE (moniou,201)s,
a,
b 7276 201
FORMAT(2
x,
'XXTWDEC: S=',
e10.3,2
x,
'A=',
e10.3,2
x,
'B=',
e10.3)
7286 IF(debug.GE.3)
WRITE (moniou,202)
xxtwdec 7287 202
FORMAT(2
x,
'XXTWDEC=',
e10.3)
7295 IMPLICIT DOUBLE PRECISION(
a-
h,
o-
z)
7298 + coef(10),
pi,zerod,halfd,oned,twod,tend
7301 DATA coef/8.3333333333333334
d-02,-2.7777777777777778
d-03,
7302 . 7.9365079365079365
d-04,-5.9523809523809524
d-04,
7303 . 8.4175084175084175
d-04,-1.9175269175269175
d-03,
7304 . 6.4102564102564103
d-03,-2.9550653594771242
d-02,
7305 . 0.1796443723688306 ,-0.6962161084529506 /
7306 DATA pi/ 3.141592653589793
d0/
7307 DATA zerod/0.
d0/,halfd/0.5
d0/,oned/1.
d0/,twod/2.
d0/,tend/10.
d0/
7316 r=(
x-halfd)* log(
x)-
x+halfd* log(twod*
pi)
7330 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
7331 COMMON /ar3/
x1(7),a1(7)
7332 COMMON /ar5/ x5(2),a5(2)
7333 COMMON /ar8/
x2(4),a2
7334 COMMON /ar9/ x9(3),a9(3)
7336 DATA x1/.9862838
d0,.9284349
d0,.8272013
d0,.6872929
d0,.5152486
d0,
7337 * .3191124
d0,.1080549
d0/
7338 DATA a1/.03511946
d0,.08015809
d0,.1215186
d0,.1572032
d0,
7339 * .1855384
d0,.2051985
d0,.2152639
d0/
7340 DATA x2/.00960736
d0,.0842652
d0,.222215
d0,.402455
d0/
7342 DATA x5/.585786
d0,3.41421
d0/
7343 DATA a5/.853553
d0,.146447
d0/
7344 DATA x9/.93247
d0,.661209
d0,.238619
d0/
7345 DATA a9/.171324
d0,.360762
d0,.467914
d0/
7351 SUBROUTINE crossc_kk(NITER,GTOT,GPROD,GABS,GDD,GQEL,GCOH)
7362 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
7365 dimension wabs(20),wdd(20),wqel(20),wcoh(20),wtot(20),
7366 *wprod(20),b0(20),ai(20),xa(64,3),xb(64,3)
7367 COMMON /area1/ ia(2),icz,icp
7368 COMMON /area6/
pi,bm,am
7369 COMMON /area16/ cc(5)
7370 COMMON /ar3/
x1(7),a1(7)
7371 COMMON /ar5/ x5(2),a5(2)
7372 COMMON /ar9/ x9(3),a9(3)
7384 b0(15-
i)=bm*sqrt((1.+
x1(
i))/2.)
7385 b0(
i)=bm*sqrt((1.-
x1(
i))/2.)
7386 ai(
i)=a1(
i)*(bm*am)**2*5.*
pi 7398 ai(14+
i)=a9(
i)*b0(14+
i)/tp*10.*am**2*
pi 7399 ai(21-
i)=a9(
i)*b0(21-
i)/tm*10.*am**2*
pi 7413 nt=nt+int(
qsran(b10)+cc(2))
7421 CALL psgea(ia(1),xa,1)
7428 CALL psgea(ia(2),xb,2)
7433 CALL gaucr(b0(
i),gabs,gdd,gqel,gcoh,xa,xb,ia(1),nt)
7434 wabs(
i)=wabs(
i)+gabs
7436 wqel(
i)=wqel(
i)+gqel
7437 wcoh(
i)=wcoh(
i)+gcoh
7447 wabs(
i)=wabs(
i)/niter
7449 wqel(
i)=wqel(
i)/niter
7450 wcoh(
i)=wcoh(
i)/niter
7451 wprod(
i)=wabs(
i)+wdd(
i)
7452 wtot(
i)=wprod(
i)+wqel(
i)+wcoh(
i)
7453 gabs=gabs+ai(
i)*wabs(
i)
7454 gdd=gdd+ai(
i)*wdd(
i)
7455 gqel=gqel+ai(
i)*wqel(
i)
7456 gcoh=gcoh+ai(
i)*wcoh(
i)
7459 gtot=gprod+gqel+gcoh
7466 SUBROUTINE gaucr(B,GABS,GDD,GQEL,GCOH,XA,XB,IA,NT)
7468 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
7469 dimension xa(64,3),xb(64,3)
7470 COMMON /area15/ fp(5),rq(5),cd(5)
7471 COMMON /area16/ cc(5)
7479 vv=1.
d0-dsqrt(
psv(xa(
n,1)+
b,xa(
n,2),xb,nt))
7480 gabs=gabs*(1.-cc(2)*(1.-vv*vv))
7481 gdd=gdd*(1.-cc(2)*(1.-vv))**2
7482 gqel=gqel*(1.-2.
d0*cc(2)*(1.-vv))
7483 gcoh=gcoh*(1.-cc(2)*(1.-vv))
7485 gcoh=1.-2.*gcoh+gqel
7495 DOUBLE PRECISION FUNCTION sectnu(E0N,IAP,IAT)
7501 IMPLICIT DOUBLE PRECISION (
a-
h,
o-
z)
7502 dimension wk(3),wa(3),wb(3)
7503 COMMON /area48/ asect(10,6,4)
7508 IF(ye.LT.1.
d0)ye=1.
d0 7513 wk(3)=wk(2)*(wk(2)-1.
d0)*.5
d0 7514 wk(1)=1.
d0-wk(2)+wk(3)
7515 wk(2)=wk(2)-2.
d0*wk(3)
7518 ya=dlog(ya/2.
d0)/.69315
d0+1.
d0 7521 wa(3)=wa(2)*(wa(2)-1.
d0)*.5
d0 7522 wa(1)=1.
d0-wa(2)+wa(3)
7523 wa(2)=wa(2)-2.
d0*wa(3)
7526 yb=dlog(yb)/1.38629
d0+1.
d0 7529 wb(3)=wb(2)*(wb(2)-1.
d0)*.5
d0 7530 wb(1)=1.
d0-wb(2)+wb(3)
7531 wb(2)=wb(2)-2.
d0*wb(3)
function psfaz(Z, FSOFT, FHARD, FSHARD)
subroutine psjint0(S, SJ, SJB, M, L)
block data include Zlatfit h c fitting region data x1(1)/0.03/
double precision function gamfun_kk(Y)
subroutine psgea(IA, XA, JJ)
subroutine pstrans(EP, EY)
*Zfirst p fm *Zfirst p Zfirst p Zfirst p *Zfirst p *Zfirst pos xyz r
! for length to thickness conversion or v v ! integer maxnodes real Hinf ! This is used when making table for dim simulation ! The slant length for vertical height to km is ! divided by LenStep m steps ! It can cover the slant length of about km for cos
integer npitbl real *nx dx real dx
subroutine xxdtg(WP0, WM0, ICP, ICT, LQ1)
subroutine xxfrag(SA, NA, RC)
function psjet(Q1, Q2, S, S2MIN, J, L)
subroutine xxdpr(WP0, WM0, ICP, ICT, LQ2)
function psjet1(Q1, Q2, S, S2MIN, J, L)
function psjint1(Q1, Q2, S, M, L)
function xxtwdec(S, A, B)
subroutine psrotat(EP, S0X, C0X, S0, C0)
function psapint(X, J, L)
subroutine psvdef(ICH, IC1, ICZ)
subroutine psjdef(IPJ, IPJ1, EPJ, EPJ1, JFL)
! structure defining a particle at production ! Basic idea of what is to be contained in ! the particle structue is that dynamical ones should be included those derivable from the particle code ! is not included ******************************************************type fmom momentum sequence union map real e endmap map real t endmap map ! pt before pz is set real pt
double precision function sectnu(E0N, IAP, IAT)
subroutine psrec(EP, QV, ZV, QM, IQV, LDAU, LPAR, IQJ, EQJ, JFL, JQ)
********************block data cblkHeavy ********************integer j data *HeavyG2symbol p
function psfsh(S, Z, ICZ, IQQ)
********************block data cblkHeavy ********************integer j data *HeavyG2symbol *data *HeavyG2code kiron data j
subroutine pstrans1(EP, EY)
subroutine xxfragm(NS, XA)
function psjint(Q1, Q2, S, M, L)
function psqint(QLMAX, G, J)
block data cblkEvhnp ! currently usable models data RegMdls ad *special data *Cekaon d0
subroutine pshot(WP0, WM0, Z, IPC, EPC, IZP, IZT, ICZ, IQQ)
subroutine xxdec3(EP, EP1, EP2, EP3, SWW, AM1, AM2, AM3)
subroutine xxaini(E0N, ICP0, IAP, IAT)
subroutine xxddfr(WP0, WM0, ICP, ICT)
function psbint(QQ, S, M, L)
subroutine crossc_kk(NITER, GTOT, GPROD, GABS, GDD, GQEL, GCOH)
! constants thru Cosmos real ! if multiplied to deg radian Torad ! light velocity m sec ! infinty ! kg m2 *Togpcm2 g cm2 ! g cm2 *Tokgpm2 kg m2 ! cm *Tom m ! m *Tocm cm ! g cm3 *Tokgpm3 kg m3 ! kg m3 *Togpcm3 g cm3 ! sec *Tonsec nsec ! Tesla m eh
! constants thru Cosmos real * pi
latitude latitude this system is used *****************************************************************! type coord sequence union map real y
function psfborn(S, T, IQ1, IQ2)
subroutine ixxdef(ICH, IC1, IC2, ICZ)
dE dx *! Nuc Int sampling table d
function ixxson(NS, AW, G)
! to be included just before the execution code ! density as a function of height real fd0 real z0
dE dx *! Nuc Int sampling table b
block data include Zlatfit h c fitting region data data data data data d0 data data d0 data data m
subroutine xxgener(WP0, WM0, EY0, S0X, C0X, S0, C0, IC1, IC2)
dE dx *! Nuc Int sampling table h
dE dx *! Nuc Int sampling table g
function psrjint(YJ, Z0, IQQ)
subroutine xxdec2(EP, EP1, EP2, WW, A, B)
*************************block data cblkTracking *************************implicit none data *ExactThick *Freec *RatioToE0 *MagChgDist *TimeStructure *Truncn *Truncx data *IncMuonPolari *KEminObs *ThinSampling *EthinRatio *Generate *LpmEffect *MagPairEmin e10
real *8 function qsran(X)
function psborn(QQ, S, IQ1, IQ2)
subroutine gaucr(B, GABS, GDD, GQEL, GCOH, XA, XB, IA, NT)
subroutine psdefrot(EP, S0X, C0X, S0, C0)
block data include Zlatfit h c fitting region data x2(1)/0.5/data x1(2)/0.3/
subroutine psdeftr(S, EP, EY)
function psv(X, Y, XB, IB)
function psrejs(S, Z, IQQ)
function psudint(QLMAX, J)
subroutine psshar(LS, NHP, NW, NT)
block data cblkIncident data *Za1ry is
! structure defining a particle at production ! Basic idea of what is to be contained in ! the particle structue is that dynamical ones should be included those derivable from the particle code ! is not included ******************************************************type fmom momentum sequence union map real e endmap map real * x
subroutine xxreg(EP0, IC)
subroutine xxstr(WPI0, WMI0, WP0, WM0, IC10, IC120, IC210, IC20)
subroutine pscajet(QQ, IQ1, QV, ZV, QM, IQV, LDAU, LPAR, JQ)
dE dx *! Nuc Int sampling table f
dE dx *! Nuc Int sampling table c
function psroot(QLMAX, G, J)