3 implicit double precision(a-y)
5 double precision pt1(2),pt2(2),ptx(2),q(4,20)
6 double precision pboo(4),pcm(4),plb(4)
7 double precision tvec(4),uvec(4),svec(4),vv1(4),vv2(4)
8 integer d,h,i,j,k,l,m,n,o,p,r,iset,af,nhist,ntotal,eflag,ichi,
9 &icut,ncut,id(20),id0,id1,id2,pdf,num,ptgen,mode,iord,ii,ip,ll,nev
11 character prefix*50,fsp*10,order*10,pflag*10,fsi*10,formf*10
12 &,ppbar*10,output*10,mregge*10,cuts*10,unw*10
14 double precision cc0(5),bm(5),bb0(5),pp0(5),bex(5),gaa(5)
15 common/pars/cc0,bm,bb0,pp0,bex,asp,sigo,gaa,ep,norm
19 common/mompt/pt1,pt2,ptx
20 common/vars/s,rts,mmes,yx,iin
22 common/vars0/ mf127, mf1525, m0, mmes0, mmes1, mmes2, mp,
23 & mwidth, pi, rt2, ebeam, sum, sum1, weightm,
25 common /ivars/ ncut, ll, icut, nev, num, ntotal
27 common/cuts/etaelmax,etaelmin,ptelmin,ptphmin,ecut,rmax,rmin,mcut
28 common/flags/ pflag, fsi, ppbar, output, cuts, unw
29 common/gluons/g1,g2,kt1,kt2
30 common/levicivita/epsilon
31 common/ang/cost,costa,costb
34 common/wvars/sig0,bb,t1,t2
35 common/alphas/alphap,alpha0,alphapr,alpha0r,alphapm,alpha0m
36 common/sec/cpom,cf,crho,aff,ar
40 common/ffpars/bexp,ao,bo,aoo,boo
44 parameter (nmxhep=4000)
45 integer nevhep,nhep,isthep,idhep,jmohep,jdahep
46 double precision phep,vhep
47 common /hepevt/ nevhep,nhep,isthep(nmxhep),idhep(nmxhep),
48 &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep)
49 ccccc Les Houches Event Common Block
51 PARAMETER (MAXNUP=500)
52 INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP
53 DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP
54 COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,
55 & IDUP(MAXNUP),ISTUP(MAXNUP),MOTHUP(2,MAXNUP),
56 & ICOLUP(2,MAXNUP),PUP(5,MAXNUP),VTIMUP(MAXNUP),
60 double precision pin(4),am(100),pout(4,100),wt,ein
61 common/ramboc/ pin, am, pout, wt, ein,
63 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
65 c Dime MC for the central exclusive production of c
66 c meson pairs by double Pomeron exchange: c
68 c p(1) p(2) --> p(3) + M_1(5) M_2(6) + p(4) c
70 c Momenta for each event in array q(i,j), where j is c
71 c the particle label and i is the 4-momentum c
74 c 1,2 = transverse components c
75 c 3 = beam component c
78 c PDG number for ith particle in arrary ID(i) c
80 c Also gives event record in HEPEVT or LHE format c
81 c (others are available upon request) c
83 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
85 c Particles generated: c
87 c pflag='pipm' : M_1=pi^+ M_2=pi^- c
88 c pflag='pi0' : M_1=M_2=pi^0 c
89 c pflag='kpkm' : M_1=K^+ M_2=K^- c
90 c pflag='ks' : M_1=M_2=K_0 c
91 c pflag='rho' : M_1=M_2=rho_0 c
93 c with decay: rho(5) --> pi^+(7)pi^-(8) c
94 c rho(6) --> pi^+(9)pi^-(10) c
95 c according to phase space, with finite rho c
98 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
100 c User defined cuts can readily be implemented c
101 c in subroutine 'icut' c
103 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
105 c This is version 1.05 : 17 Sep 2014 c
107 c Comments etc to: lucian.harland-lang@durham.ac.uk c
109 c If you use the code please reference (and see for c
110 c details of model used): c
112 c L.A. Harland-Lang, V.A. Khoze, M.G. Ryskin c
113 c and W.J. Stirling arXiv:1105.1626 c
115 c L.A. Harland-Lang, V.A. Khoze, M.G. Ryskin c
116 c and W.J. Stirling arXiv:1204.4803 c
118 c L.A. Harland-Lang, V.A. Khoze and M.G. Ryskin c
121 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
123 ccAM open(35,file='exerec1.dat') ! event record
125 ccAM rts=1.96d3 ! centre of mass energy
126 cccc Some basic cuts, imposed in subtroutine 'icut'. Other user defined cuts can readily be implemented in subroutine
127 cccc note: in rhorho case these cuts are imposed on rho's, and not their decay productions. Cuts on the decay products
128 cccc can be imposed by editing 'icut'
130 ccAM rmax=1.8d1 ! max meson rapidity
131 ccAM rmin=-1.8d1 ! min meson rapidity
132 ccAM ecut=2d0 ! min meson p_t
134 c -- physics parameters
136 ccAM pflag='rho' ! Process generated - see preamble for options
137 ccAM fsi='true' ! phenomenological model for exclusive suppression in Pom Pom --> meson pair. To turn on/off --> 'true/false'
138 ccAM formf='orexp' ! meson - pomeron form factor.
139 ccAM ppbar='false' ! set true if ppbar collisions
140 output='lhe' ! Event record style, HEPEVT or LHE
141 ccAM cuts='true' ! Impose cuts or not
142 ccAM unw='true' ! Set = 'true' for unweighted events
143 ccAM iin=1 ! Model for soft survival factor, as described in arXiv:1306.2149. Default = 1
144 write(6,*) "*************************************************"
145 write(6,*) "Initialising DIME with:"
146 if (ppbar.eq.'true') then
147 write(6,*) "ppbar collisions"
149 write(6,*) "pp collisions"
151 write(6,*) "Centre of Mass Energy = ", rts
152 write(6,*) "Physics models:"
153 write(6,*) "Process = ", pflag
154 write(6,*) "Eclusive suppression in Pom Pom = ", fsi
155 write(6,*) "Meson - Pomeron form factor = ", formf
156 write(6,*) "Unweighted events = ", unw
157 write(6,*) "Model for soft survival factor = ", iin
158 write(6,*) "*************************************************"
159 write(6,*) "Using cuts = ", cuts
160 write(6,*) "Maximum meson rapidity =", rmax
161 write(6,*) "Minimum meson rapidity =", rmin
162 write(6,*) "Minimum meson pT =", ecut
163 write(6,*) "*************************************************"
166 ntotal=1000000 ! no. of runs for weighted events
167 nev=100 ! no. of unweighted events generated to event record
171 if(formf.eq.'exp')then ! Set parameters for Pomeron-Meson form factor in function 'fpi(x)'
173 elseif(formf.eq.'orexp')then
176 elseif(formf.eq.'power')then
180 ccccccccccccccccccccccccccccccccccccccc
182 ccccccccc Start of main body of code
184 cccccccccccccccccccccccccccccccccccccc
186 call initpars(iin) ! Initialise soft survival parameters
187 call calcop ! proton opacity
188 call calcscreen ! screening amplitude
189 if(pflag.eq.'pipm')then
190 mmes=0.13957018d0 ! pi+/- mass, PDG 2011 value
193 alpha0m=-0.7d0*mmes**2
216 elseif(pflag.eq.'pi0')then
217 mmes=0.1349766d0 ! pi0 mass, PDG 2011 value
220 alpha0m=-0.7d0*mmes**2
243 elseif(pflag.eq.'kpkm')then
244 mmes=0.493677d0 ! K+/- mass, PDG 2011 value
268 elseif(pflag.eq.'ks')then
269 mmes=0.497614d0 ! K_0 mass, PDG 2011 value
293 elseif(pflag.eq.'rho')then
294 mmes0=0.77549d0 ! rho mass, PDG 2013 value
295 mwidth=0.1491d0 ! rho width PDG 2013 value
339 c -- other parameters
348 beta=dsqrt(1d0-4d0*mp**2/s)
357 alphap=0.25d0 ! D-L '92 fit
372 ccccc Initialise rambo (rho0 decay)
374 if(pflag.eq.'rho')then
385 c initialise counters etc
396 c initialise histograms
404 cccccccccccccccccccccccccccccccccccccccccccccccccccccc
406 c Incoming protons to event record array
408 cccccccccccccccccccccccccccccccccccccccccccccccccccccc
425 pup(5,1)=dsqrt(q(4,1)**2-q(3,1)**2-q(2,1)**2-q(1,1)**2)
435 if(ppbar.eq.'true')then
447 pup(5,2)=dsqrt(q(4,2)**2-q(3,2)**2-q(2,2)**2-q(1,2)**2)
451 cccc outgoing initial info
454 if(ppbar.eq.'true')then
477 if(output.eq.'hepevt')then
488 jmohep(1,k)=mothup(1,k)
489 jmohep(2,k)=jmohep(1,k)
498 if(pflag.eq.'rho')then
510 if(unw.eq.'true')then
517 c first generate weighted events first
520 write(6,*)'Generating weighted events...'
523 call dimegenerate(isuc)
525 c prepare for unweighted generation
527 write(6,*)'Generating unweighted events...'
532 subroutine dimegenerate(isuccess)
537 parameter (nmxhep=4000)
538 integer nevhep,nhep,isthep,idhep,jmohep,jdahep
539 double precision phep,vhep
540 common /hepevt/ nevhep,nhep,isthep(nmxhep),idhep(nmxhep),
541 &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep)
542 ccccc Les Houches Event Common Block
544 PARAMETER (MAXNUP=500)
545 INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP
546 DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP
547 COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,
548 & IDUP(MAXNUP),ISTUP(MAXNUP),MOTHUP(2,MAXNUP),
549 & ICOLUP(2,MAXNUP),PUP(5,MAXNUP),VTIMUP(MAXNUP),
551 ccccc rambo variables
553 double precision pin(4), am(100), pout(4,100), wt, ein
554 common/ramboc/ pin, am, pout, wt, ein,
558 double precision etaelmax,etaelmin,ptelmin,ptphmin,ecut,rmax,
560 common/cuts/etaelmax,etaelmin,ptelmin,ptphmin,ecut,rmax,rmin,mcut
562 double precision q(4,20)
565 double precision pt1(2), pt2(2), ptx(2)
566 common/mompt/pt1,pt2,ptx
568 double precision s, rts, mmes, yx
569 double precision bjac, bp
570 double precision mf127, mf1525
572 double precision m0, mmes0, mmes1, mmes2, mp, mwidth, pi, rt2,
573 &ebeam, sh, th, uh, sum, sum1, weightm
575 integer id(20), ncut, ll, icut, nev, num, iin, ntotal
577 common/vars/s,rts,mmes,yx, iin
578 common/hvars/sh,th,uh
579 common/vars0/ mf127, mf1525, m0, mmes0, mmes1, mmes2, mp,
580 & mwidth, pi, rt2, ebeam, sum, sum1, weightm,
582 common /ivars/ ncut, ll, icut, nev, num, ntotal
584 character prefix*50,fsp*10,order*10,pflag*10,fsi*10,formf*10
585 &,ppbar*10,output*10,mregge*10,cuts*10,unw*10
586 common/flags/pflag, fsi, ppbar, output, cuts, unw
589 double precision ep, norm, sigo, asp
590 double precision cc0(5),bm(5),bb0(5),pp0(5),bex(5),gaa(5)
591 common/pars/cc0,bm,bb0,pp0,bex,asp,sigo,gaa,ep,norm
592 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
593 cccccc local variables ...
594 double precision aa1, aa2, al1, al2, almax, almin, c, cc1, cc2,
595 & msmax, msmin, cont, exn, mmesp, mwidth1, mwidth2,
596 & p1m, p1p, p2m, p2p, phi1, phi2, phix1, pt1sq, pt1x, pt1y,
597 & pt2sq, pt2x, pt2y, ptxsq1, ptxsq2, ptxsqmax, ptxsqmin,
598 & ran0, ran1, ran2, ran3, ran4, ran5, ran6, ranhist, ranx1,
599 & ranx2, ranx3, ranx4, rm1, rm2, rmx1, rmx2,
600 & root1sq, root2sq, snp, t1, t2, weight, weight0, weight1,
601 & weight2, wthist, x1, x2, ymax1, ymax2, ymin1, ymin2, yx1,
603 double precision pboo(4), pcm(4), plb(4)
604 double precision svec(4)
606 integer h, i, j, k, isuccess
617 ccccccccccccccccccccccccccccccccccccccccccccccccccc
636 pt1sq=-dlog(ran2)/(bjac)
637 pt2sq=-dlog(ran3)/(bjac)
639 weight=(dexp(bjac*pt1sq)*dexp(bjac*pt2sq))/bjac**2
641 pt1(1)=dsqrt(pt1sq)*dsin(phi1)
642 pt1(2)=dsqrt(pt1sq)*dcos(phi1)
643 pt2(1)=dsqrt(pt2sq)*dsin(phi2)
644 pt2(2)=dsqrt(pt2sq)*dcos(phi2)
651 cccccccccccccccccccccccccccccccc
653 ccc non-zero rho width
655 cccccccccccccccccccccccccccccccc
657 if(pflag.eq.'rho')then ! new part here
662 msmax=mmes0+4d0*mwidth
663 msmin=2d0*0.13957018d0
667 almin=datan(-(mmes0**2-msmin**2)/mwidth/mmes0)
668 almax=datan(-(mmes0**2-msmax**2)/mwidth/mmes0)
670 al1=almin+(almax-almin)*rm1
671 al2=almin+(almax-almin)*rm2
673 mmes1=dsqrt(dtan(al1)*mmes0*mwidth+mmes0**2)
674 mmes2=dsqrt(dtan(al2)*mmes0*mwidth+mmes0**2)
676 weight=weight*(almax-almin)**2
677 weight=weight*mwidth**2*mmes0**2
678 weight=weight*(1d0+dtan(al1)**2)*(1d0+dtan(al2)**2)
679 weight=weight/4d0/mmes1/mmes2
683 mwidth1=mwidth*((1d0-4d0*0.13957018d0**2/mmes1**2)/
684 & (1d0-4d0*0.13957018d0**2/mmes0**2))**1.5d0
685 mwidth2=mwidth*((1d0-4d0*0.13957018d0**2/mmes2**2)/
686 & (1d0-4d0*0.13957018d0**2/mmes0**2))**1.5d0
688 if(mmes1.lt.2d0*0.13957018d0) goto 677
689 if(mmes2.lt.2d0*0.13957018d0) goto 677
691 weight=weight*2d0*mmes0*mmes1*mwidth1/pi
692 weight=weight*2d0*mmes0*mmes2*mwidth2/pi
693 weight=weight*mmes1**2*mmes2**2/mmes0**4
694 weight=weight/((mmes0**2-mmes1**2)**2+mmes1**4*
695 & mwidth1**2/mmes0**2)
696 weight=weight/((mmes0**2-mmes2**2)**2+mmes2**4*
697 & mwidth2**2/mmes0**2)
702 cccccccccccccccccccccccccccccccccccccccccccccccccccc
715 c ptxsq1=-dlog(r2)/bjac1 ! generate exp p_t^2 (can be more efficient)
718 ptxsqmax=10d0 ! generate linear p_t^2
720 ptxsq1=ptxsqmin+(ptxsqmax-ptxsqmin)*ranx2
722 q(1,5)=dsqrt(ptxsq1)*dcos(phix1)
723 q(2,5)=dsqrt(ptxsq1)*dsin(phix1)
724 q(1,6)=-pt1(1)-pt2(1)-q(1,5)
725 q(2,6)=-pt1(2)-pt2(2)-q(2,5)
726 ptxsq2=q(1,6)**2+q(2,6)**2
728 rmx1=dsqrt(mmes1**2+ptxsq1)
729 rmx2=dsqrt(mmes2**2+ptxsq2)
733 if(cuts.eq.'true')then
734 if(rmax.lt.ymax1)then
741 yx1=ymin1+(ymax1-ymin1)*ranx3
743 ymax2=(dlog((rts-rmx1*dexp(yx1))/rmx2))
744 ymin2=-(dlog((rts-rmx1*dexp(-yx1))/rmx2))
746 if(cuts.eq.'true')then
747 if(ymax2.gt.rmax)then
750 if(ymin2.lt.-rmax)then
755 if(ymax2.lt.ymin2)then
760 yx2=ymin2+(ymax2-ymin2)*ranx4
761 x1=(rmx2*dexp(yx2)+rmx1*dexp(yx1))/rts
762 x2=(rmx2*dexp(-yx2)+rmx1*dexp(-yx1))/rts
764 weight=weight*(ptxsqmax-ptxsqmin)
765 c weight=weight*dexp(bjac1*ptxsq1)*r2max/bjac1
766 weight=weight*(ymax1-ymin1)*(ymax2-ymin2)
769 q(3,5)=rmx1*(dexp(yx1)-dexp(-yx1))/2d0
770 q(4,5)=rmx1*(dexp(yx1)+dexp(-yx1))/2d0
772 q(3,6)=rmx2*(dexp(yx2)-dexp(-yx2))/2d0
773 q(4,6)=rmx2*(dexp(yx2)+dexp(-yx2))/2d0
775 ccccccccccccccccccccccccccccccccccccccccccccccccccccc
777 c impose massive on-shell condition by solving
778 c p1+ + cc1/p2- = aa1
779 c p2- + cc2/p1+ = aa2
783 cc1=0.5d0*(pt2sq+mp**2)
784 cc2=0.5d0*(pt1sq+mp**2)
786 root1sq=(cc1-cc2-aa1*aa2)**2-4d0*cc2*aa1*aa2
787 root2sq=(cc2-cc1-aa1*aa2)**2-4d0*cc1*aa1*aa2
788 if(root1sq.le.0d0.or.root2sq.le.0d0)then
792 p1p=(cc2-cc1+aa1*aa2+dsqrt(root1sq))/(2d0*aa2)
793 p2m=(cc1-cc2+aa1*aa2+dsqrt(root2sq))/(2d0*aa1)
794 p1m=(pt1sq+mp**2)/(2d0*p1p)
795 p2p=(pt2sq+mp**2)/(2d0*p2m)
797 if(p1p.lt.0d0.or.p1m.lt.0d0.or.p2p.lt.0d0.or.p2m.lt.0d0)then
815 ccccccccccccccccccccccccccccccccccccccccccccccccc
818 svec(k)=q(k,5)+q(k,6)
821 call dot(svec,svec,sh)
823 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
827 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
829 if(pflag.eq.'rho')then
840 call rambo(npart,mmesp,am,pout,wt) ! call RAMBO
850 call boost(mmesp,pboo,pcm,plb)
862 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
866 if(cuts.eq.'true')then
882 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
884 c Write 4-momenta to event record array
886 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
888 if(output.eq.'lhe')then
894 pup(5,k)=dsqrt(q(4,k)**2-q(3,k)**2-q(2,k)**2-q(1,k)**2)
897 if(pflag.eq.'rho')then
903 pup(5,k)=dsqrt(q(4,k)**2-q(3,k)**2-q(2,k)**2-q(1,k)**2)
912 pup(5,k)=dsqrt(q(4,k)**2-q(3,k)**2-q(2,k)**2-q(1,k)**2)
923 phep(5,k)=dsqrt(q(4,k)**2-q(3,k)**2-q(2,k)**2-q(1,k)**2)
926 if(pflag.eq.'rho')then
932 phep(5,k)=dsqrt(q(4,k)**2-q(3,k)**2-q(2,k)**2-q(1,k)**2)
941 phep(5,k)=dsqrt(q(4,k)**2-q(3,k)**2-q(2,k)**2-q(1,k)**2)
948 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
952 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
954 c pp-->pp(M_1 M_2) matrix element
956 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
959 call schimc(pt1x,pt1y,pt2x,pt2y,zmat)
961 weight=weight*cdabs(zmat)**2
962 weight=weight/norm**2
964 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
966 weight=weight/(s**2*(16d0)**3*pi**5)*389d3
968 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
970 ccc Probability of no additional particles produced in production subprocess
972 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
974 if(fsi.eq.'true')then
976 if(pflag.eq.'pipm'.or.pflag.eq.'pi0')then
978 if(dsqrt(sh).gt.mf127)then
981 cont=1d0-(m0-2d0*mmes)**2/m0**2
982 exn=c*dlog(((dsqrt(sh)-2d0*mmes)**2)/m0**2+cont)
988 elseif(pflag.eq.'kpkm'.or.pflag.eq.'ks')then
990 if(dsqrt(sh).gt.mf1525)then
993 cont=1d0-(m0-2d0*mmes)**2/m0**2
994 exn=c*dlog(((dsqrt(sh)-2d0*mmes)**2)/m0**2+cont)
1020 if(pflag.eq.'pi0'.or.pflag.eq.'ks'.or.pflag.eq.'rho')then
1021 weight=weight/2d0 ! symmetry factor
1024 cccccccccccccccccccccccccccccccccccccccccccccccccccccc
1028 if(weightm.lt.weight)then
1037 wthist=weight/dfloat(ntotal)
1045 if(ranhist.lt.weight/weightm)then
1052 c call binit(sumt/nev)
1054 write(35,302) num, nup
1056 if(output.eq.'lhe')then
1059 write(35,300)j,istup(j),idup(j),mothup(1,j),mothup(2,j),
1060 & icolup(1,j),icolup(2,j),pup(1,j),pup(2,j),
1061 & pup(3,j),pup(4,j),pup(5,j),vtimup(j),spinup(j)
1067 write(35,301)j,idhep(j),isthep(j),jmohep(1,j),
1068 & jmohep(2,j),jdahep(1,j),jdahep(2,j),
1069 & phep(1,j),phep(2,j),phep(3,j),phep(4,j),phep(5,j)
1070 & ,vhep(1,j),vhep(2,j),vhep(3,j),vhep(4,j)
1082 if(num.gt.nev-1)then ! exit loop once required no. of unweighted events generated
1092 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1094 c End of event loop c
1096 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1097 c call terminate(ntotal)
1098 300 format(i4,1x,i4,1x,i8,1x,i4,1x,i4,1x,i4,1x,i4,1x,E13.6,1x,
1099 &E13.6,1x,E13.6,1x,E13.6,1x,E13.6,1x,E13.6,1x,E13.6)
1100 301 format(i4,1x,i8,1x,i4,1x,i4,1x,i4,1x,i4,1x,i4,1x,E13.6,1x,E13.6
1101 &,1x,E13.6,1x,E13.6,1x,E13.6,1x,E13.6,1x,E13.6,1x,E13.6,1x
1103 302 format(1x,i8,1x,i5,1x,E13.6)
1108 subroutine terminate(int ntotal)
1109 character prefix*50,fsp*10,order*10,pflag*10,fsi*10,formf*10
1110 &,ppbar*10,output*10,mregge*10,cuts*10,unw*10
1113 sum=sum/dfloat(ntotal)
1114 sum1=sum1/dfloat(ntotal)
1115 var=dsqrt((sum1-sum**2)/dfloat(ntotal))
1116 EFF=1d0*(ntotal-ncut)/ntotal
1129 ccccc create histograms
1140 if(pflag.eq.'pipm')then
1142 write(6,*)'pi+pi- production'
1143 elseif(pflag.eq.'pi0')then
1145 write(6,*)'pi0pi0 production'
1146 elseif(pflag.eq.'kpkm')then
1148 write(6,*)'K+K- production'
1149 elseif(pflag.eq.'ks')then
1151 write(6,*)'K0K0 production'
1152 elseif(pflag.eq.'rho')then
1154 write(6,*)'rho0rho0 production'
1159 if(pflag.eq.'rho')then
1160 write(6,221)rts,ntotal,sum,var
1162 write(6,222)rts,mmes,ntotal,sum,var
1165 if(output.eq.'lhe')then
1168 print*,'HEPEVT output'
1170 write(6,223)nev,sum,var
1178 . 3x,'collider energy (GeV) : ',f10.3,/,
1179 . 3x,'number of events : ',i9,/,
1180 . 3x,'sigma (nb) : ',f15.5,' +- ',f15.5,
1184 . 3x,'collider energy (GeV) : ',f10.3,/,
1185 . 3x,'meson mass (GeV) : ',f10.5,/,
1186 . 3x,'number of events : ',i9,/,
1187 . 3x,'sigma (nb) : ',f15.5,' +- ',f15.5,
1190 223 format(3x,'number of events : ',i9,/,
1191 . 3x,'sigma (nb) : ',f15.5,' +- ',f15.5,
1195 print*,'time elapsed = ', t2, ' s'
1200 ccccc Pomeron -- (off-shell) meson form factor
1203 implicit double precision(a-y)
1205 common/vars/s,rts,mmes,yx, iin
1207 common/ffpars/bexp,ao,bo,aoo,boo
1209 if(formf.eq.'exp')then
1210 fpi=exp((x-mmes**2)*bexp)
1211 elseif(formf.eq.'orexp')then
1212 fpi=exp(-(dsqrt(-x+mmes**2+ao**2)-ao)*bo)
1213 elseif(formf.eq.'power')then
1214 fpi=1d0/(1d0-(x-mmes**2)/aoo)
1220 ccccccc Pom Pom --> meson pair amplitude
1222 subroutine wev(matt,v1,v2)
1223 implicit double precision(a-y)
1224 implicit complex*16(z)
1225 double precision q(4,20),svec(4),tvec(4),uvec(4),v1(4),v2(4)
1227 character*10 mregge,pflag
1230 common/vars/s,rts,mmes,yx, iin
1231 common/wvars/sig0,bb,t1,t2
1232 common/alphas/alphap,alpha0,alphapr,alpha0r,alphapm,alpha0m
1233 common/sec/cpom,cf,crho,aff,ar
1234 common/hvars/sh,th,uh
1236 common/flags/ pflag, fsi, ppbar, output, cuts, unw
1248 tvec(k)=q(k,10)-q(k,5)
1249 uvec(k)=q(k,10)-q(k,6)
1252 call dot(tvec,tvec,th)
1253 call dot(uvec,uvec,uh)
1255 sht1=2d0*(q(4,5)*q(4,3)-q(3,5)*q(3,3)-q(2,5)*q(2,3)-
1256 &q(1,5)*q(1,3))+mmes**2
1258 sht2=2d0*(q(4,6)*q(4,4)-q(3,6)*q(3,4)-q(2,6)*q(2,4)-
1259 &q(1,6)*q(1,4))+mmes**2
1261 shu1=2d0*(q(4,6)*q(4,3)-q(3,6)*q(3,3)-q(2,6)*q(2,3)-
1262 &q(1,6)*q(1,3))+mmes**2
1264 shu2=2d0*(q(4,5)*q(4,4)-q(3,5)*q(3,4)-q(2,5)*q(2,4)-
1265 &q(1,5)*q(1,4))+mmes**2
1267 bt1=(2d0*alphap*dlog(sht1))
1268 bt2=(2d0*alphap*dlog(sht2))
1269 bu1=(2d0*alphap*dlog(shu1))
1270 bu2=(2d0*alphap*dlog(shu2))
1272 bt1r=(2d0*alphapr*dlog(sht1))
1273 bt2r=(2d0*alphapr*dlog(sht2))
1274 bu1r=(2d0*alphapr*dlog(shu1))
1275 bu2r=(2d0*alphapr*dlog(shu2))
1277 sig0t1=sig0*sht1**0.0808d0/0.389d0
1278 sig0t2=sig0*sht2**0.0808d0/0.389d0
1279 sig0u1=sig0*shu1**0.0808d0/0.389d0
1280 sig0u2=sig0*shu2**0.0808d0/0.389d0
1282 t11=v1(1)**2+v1(2)**2
1283 t22=v2(1)**2+v2(2)**2
1285 zmt=(zi*dexp(-bt1*t11/2d0)*cpom*sht1**alpha0
1286 & +((aff+zi)*cf*sht1**alpha0r+(ar-zi)*crho*sht1**alpha0r)
1287 & *dexp(-bt1r*t11/2d0))
1288 zmt=zmt*(zi*dexp(-bt2*t22/2d0)*cpom*sht2**alpha0
1289 & +((aff+zi)*cf*sht2**alpha0r-(ar-zi)*crho*sht2**alpha0r)
1290 & *dexp(-bt2r*t11/2d0))
1291 zmt=zmt*fpi(th)**2/(mmes**2-th)
1293 zmu=(zi*dexp(-bu1*t11/2d0)*cpom*shu1**alpha0
1294 & +((aff+zi)*cf*shu1**alpha0r-(ar-zi)*crho*shu1**alpha0r)
1295 & *dexp(-bu1r*t11/2d0))
1297 zmu=zmu*(zi*dexp(-bu2*t22/2d0)*cpom*shu2**alpha0
1298 & +((aff+zi)*cf*shu2**alpha0r+(ar-zi)*crho*shu2**alpha0r)
1299 & *dexp(-bu2r*t22/2d0))
1300 zmu=zmu*fpi(uh)**2/(mmes**2-uh)
1306 c binning subroutine
1308 subroutine binit(wt)
1309 implicit double precision(a-y)
1310 double precision q(4,20),pt1(2),pt2(2),ptx(2)
1312 common/mompt/pt1,pt2,ptx
1313 common/vars/s,rts,mmes,yx, iin
1314 common/vars1/ptgam,etagam,ptel2,ptel1,etael1,etael2
1315 common/vars2/ptpi1,etapi1,ptpi2,etapi2
1316 common/ang/cost,costa,costb
1317 common/hvars/sh,th,uh
1319 ypisum=0.5d0*dlog((q(4,5)+q(4,6)+q(3,5)+q(3,6))
1320 & /(q(4,5)+q(4,6)-q(3,5)-q(3,6)))
1322 ypi1=0.5d0*dlog((q(4,5)+q(3,5))
1324 ypi2=0.5d0*dlog((q(4,6)+q(3,6))
1327 call histo1(1,30,0d0,4.5d0,dsqrt(sh),wt)
1332 ccccc Cutting subroutine
1334 ccccc User can define their own cuts on any particle momenta, stored in array q(i,j)
1337 subroutine cut(icut)
1338 implicit double precision(a-y)
1339 double precision q(4,20)
1342 common/vars/s,rts,mmes,yx, iin
1343 common/vars1/ptgam,etagam,ptel2,ptel1,etael1,etael2
1344 common/vars2/ptpi1,etapi1,ptpi2,etapi2
1345 common/cuts/etaelmax,etaelmin,ptelmin,ptphmin,ecut,rmax,rmin,mcut
1346 common/hvars/sh,th,uh
1352 c -- meson 1 variables
1353 ptpi1=dsqrt(q(1,5)**2+q(2,5)**2)
1354 pmodpi1=dsqrt(q(1,5)**2+q(2,5)**2+q(3,5)**2)
1355 etapi1=.5d0*dlog((pmodpi1+q(3,5))
1356 . /(pmodpi1-q(3,5)))
1358 ypi1=.5d0*dlog((q(4,5)+q(3,5))
1361 etpi1=q(4,5)*dsqrt(pmodpi1**2-q(3,5)**2)/pmodpi1
1363 c -- meson 2 variables
1364 ptpi2=dsqrt(q(1,6)**2+q(2,6)**2)
1365 pmodpi2=dsqrt(q(1,6)**2+q(2,6)**2+q(3,6)**2)
1366 etapi2=.5d0*dlog((pmodpi2+q(3,6))
1367 . /(pmodpi2-q(3,6)))
1369 ypi2=.5d0*dlog((q(4,6)+q(3,6))
1372 etpi2=q(4,6)*dsqrt(pmodpi2**2-q(3,6)**2)/pmodpi2
1374 c print*,ptpi1,ptpi2
1376 ccccc example cuts....
1378 c if(dsqrt(sh).lt.0.8d0)return
1383 c if(dabs(xf1).lt.0.9d0)return
1384 c if(dabs(xf2).lt.0.9d0)return
1386 ccccccccc Basic cuts described at start of code
1388 if(ptpi1.lt.ecut)return
1389 if(ptpi2.lt.ecut)return
1390 if(ypi1.gt.rmax)return
1391 if(ypi2.gt.rmax)return
1392 if(ypi1.lt.rmin)return
1393 if(ypi2.lt.rmin)return
1402 subroutine histo1(ih,ib,x0,x1,x,w)
1403 implicit real*8(a-h,o-z)
1404 character*1 regel(30),blank,star
1405 dimension h(20,100),hx(20),io(20),iu(20),ii(20)
1406 dimension y0(20),y1(20),ic(20)
1407 data regel / 30*' ' /,blank /' ' /,star /'*'/
1409 open(10,file="output.dat")
1415 ix=idint((x-x0)/(x1-x0)*dble(ib))+1
1417 if(h(ih,ix).gt.hx(ih)) hx(ih)=h(ih,ix)
1428 bsize=(x11-x01)/dble(ib1)
1429 hx(ih)=hx(ih)*(1.d0+1.d-06)
1430 if(il.eq.0) write(6,21) ih,ii(ih),iu(ih),io(ih)
1431 if(il.eq.1) write(6,22) ih,ii(ih),iu(ih),io(ih)
1432 21 format(' no.',i3,' lin : inside,under,over ',3i6)
1433 22 format(' no.',i3,' log : inside,under,over ',3i6)
1434 if(ii(ih).eq.0) goto 28
1436 23 format(35(1h ),3(10h----+----i))
1438 z=(dble(iv)-0.5d0)/dble(ib1)*(x11-x01)+x01
1440 iz=idint(h(ih,iv)/hx(ih)*30.)+1
1443 if(h(ih,iv).gt.0.d0)
1444 .iz=idint(dlog(h(ih,iv))/dlog(hx(ih))*30.)+1
1445 25 if(iz.gt.0.and.iz.le.30) regel(iz)=star
1446 write(6,26) z,h(ih,iv)/bsize,(regel(i),i=1,30)
1447 c write(10,*)z-0.125d0/2d0,h(ih,iv)/bsize ! Print histogram to file
1448 write(10,*)z,h(ih,iv)/bsize ! Print histogram to file
1449 c write(10,*)z,h(ih,iv)/bsize ! Print histogram to file
1450 26 format(1h ,2g15.6,4h i,30a1,1hi)
1451 36 format(1h ,2g15.6)
1452 if(iz.gt.0.and.iz.le.30) regel(iz)=blank
1457 29 format(' no entries inside histogram')
1471 cccc Initializes soft model parameters
1473 subroutine initpars(in)
1474 implicit double precision(a-y)
1477 double precision cc0(5),bm(5),bb0(5),pp0(5),bex(5),gaa(5)
1478 common/pars/cc0,bm,bb0,pp0,bex,asp,sigo,gaa,ep,norm
1480 common/vars/s,rts,mmes,yx, iin
1567 ep=0.11d0 !!! Capital delta
1568 asp=0.06d0 !!! alpha'
1569 ep1=0d0 !!! Zero in all models, matters (?)
1570 sigo=50d0/0.39d0 !!! sigma_0 (GeV^-2)
1571 gaa3=6d0 !!! k2/k(1.8 TeV)
1572 gd2=1.3d0 !!! k1/k(1.8 TeV)
1574 nga=1 !!! A flag (?)
1576 bm(1)=3d0 !!! Doesn't matter
1577 bb0(1)=0.45d0 !!! c1-0.08 (added back later)
1578 pp0(1)=0.5d0 !!! 2*|a_1|^2
1581 cc0(2)=0.48d0 !!! d2
1582 bm(2)=1.5d0 !!! Doesn't matter
1583 bb0(2)=0.16d0 !!! c2-0.08 (added back later)
1584 pp0(2)=1d0 !!! |a_2|^2 is set later
1587 cc0(3)=0.12d0 !!! Beta : k^2_min ~ s^Beta
1591 if(nch.eq.3) pp0(3)=3d0-pp0(2)-pp0(1)
1592 if(nch.eq.2) pp0(2)=2d0-pp0(1) !!! Set |a_2|^2
1593 if(nch.eq.1) pp0(1)=1d0
1595 if(in.eq.3.or.in.eq.4)then
1597 gamm=(1800d0/rts)**cc0(3)
1598 ga1=1d0/(1d0+gamm*gd2)
1599 ga2=1d0/(1d0+gamm*gaa3)
1600 gaa(1)=2d0*ga1/(ga1+ga2) !!! gamma_1 (+) Eq (15)
1601 gaa(2)=2d0*ga2/(ga1+ga2) !!! gamma_2 (-) Eq (15)
1603 elseif(in.eq.1.or.in.eq.2)then
1606 gaa(1)=1d0+dsqrt(gd2)
1607 gaa(2)=1d0-dsqrt(gd2)
1609 elseif(nch.eq.1)then
1621 sum=sum+gaa(i1)*gaa(i2)*pp0(i1)*pp0(i2)/dble(nch)**2
1629 ccccc Calculates proton opacity
1632 implicit double precision(a-y)
1633 integer nb,i1,i2,nch,ib
1635 double precision op(5,5,10000,2),oph(5,5,10000,2)
1641 print*,'Calculating opacity'
1650 call opacity(i1,i2,bt,fr,fr1)
1658 write(40,*)bt,fr,fr1
1667 cccc Calculates screening amplitude (in k_t space)
1669 subroutine calcscreen
1670 implicit double precision(a-y)
1671 integer ns,i1,i2,nch,ib
1673 double precision sca(5,5,0:40000,2),sca1(5,5,0:40000,2)
1674 common/spac/sca,sca1
1680 lginck=(dlog(ksqma/ksqmin))/dble(ns)
1682 print*,'Calculating screening amplitude'
1693 lgksq=dble(ib-1)*lginck+dlog(ksqmin)
1700 call screening(i1,i2,ksq,sc,sc1)
1702 sca(i1,i2,ib,1)=lgksq
1704 sca1(i1,i2,ib,1)=lgksq
1705 sca1(i1,i2,ib,2)=sc1
1715 ccccc Interpolation....
1717 subroutine opacityint(i,j,bt,out,out1)
1718 implicit double precision(a-y)
1719 double precision op(5,5,10000,2),oph(5,5,10000,2)
1723 incbt=op(1,1,2,1)-op(1,1,1,1)
1725 if(dble(it).gt.(bt/incbt))then
1729 m=(op(i,j,it+2,2)-op(i,j,it+1,2))
1730 &/(op(i,j,it+2,1)-op(i,j,it+1,1))
1731 del=bt-op(1,1,it+1,1)
1732 mh=(oph(i,j,it+2,2)-oph(i,j,it+1,2))
1733 &/(oph(i,j,it+2,1)-oph(i,j,it+1,1))
1734 delh=bt-oph(1,1,it+1,1)
1736 out=m*del+op(i,j,it+1,2)
1737 out1=mh*delh+oph(i,j,it+1,2)
1742 subroutine screeningint(i,j,ktsq,out,out1)
1743 implicit double precision(a-y)
1744 double precision sca(5,5,0:40000,2),sca1(5,5,0:40000,2)
1746 common/spac/sca,sca1
1748 if(ktsq.lt.0.001d0)then
1752 m=(sca(i,j,it+1,2)-sca(i,j,it,2))
1753 & /(dexp(sca(i,j,it+1,1))-sca(i,j,it,1))
1754 del=ktsq-sca(1,1,it,1)
1755 m1=sca1(i,j,it+1,2)-sca1(i,j,it,2)
1756 & /(dexp(sca1(i,j,it+1,1))-sca1(i,j,it,1))
1757 del1=ktsq-sca1(1,1,it,1)
1759 out=m*del+sca(i,j,it,2)
1760 out1=m1*del1+sca1(i,j,it,2)
1762 elseif(ktsq.lt.8d0)then
1765 inckt=sca(1,1,2,1)-sca(1,1,1,1)
1766 it=nint((dlog(ktsq)-ktmin)/inckt)
1767 if(dble(it).gt.((dlog(ktsq)-ktmin)/inckt))then
1771 m=(sca(i,j,it+2,2)-sca(i,j,it+1,2))
1772 & /(sca(i,j,it+2,1)-sca(i,j,it+1,1))
1773 del=dlog(ktsq)-sca(1,1,it+1,1)
1774 m1=sca1(i,j,it+2,2)-sca1(i,j,it+1,2)
1775 & /(sca1(i,j,it+2,1)-sca1(i,j,it+1,1))
1776 del1=dlog(ktsq)-sca1(1,1,it+1,1)
1778 out=m*del+sca(i,j,it+1,2)
1779 out1=m1*del1+sca1(i,j,it+1,2)
1791 cccc Integrates round Pomeron loop (to calculate screened amplitude)
1793 subroutine schimc(p1x,p1y,p2x,p2y,out)
1794 implicit double precision(a-y)
1795 integer nx,jx,jy,nch,i1,i2,it,k
1796 double precision cc0(5),bm(5),bb0(5),pp0(5),bex(5),gaa(5)
1797 double precision sca(5,5,0:40000,2),sca1(5,5,0:40000,2)
1798 double precision a1(4),a2(4),q(4,20)
1799 complex*16 out,x0,x00
1800 common/pars/cc0,bm,bb0,pp0,bex,asp,sigo,gaa,ep,norm
1802 common/spac/sca,sca1
1804 common/wvars/sig0,bb,t1,t2
1805 common/alphas/alphap,alpha0,alphapr,alpha0r,alphapm,alpha0m
1840 call screeningint(i1,i2,tp2,sc,sc1)
1842 x0=x0*pp0(i1)*pp0(i2)/dble(nch)**2
1843 x0=x0*dexp(-((t12+0.08d0+bb0(i1))*bex(i1))**cc0(i1)+
1844 & (bex(i1)*(bb0(i1)+0.08d0))**cc0(i1))
1845 x0=x0*dexp(-((t22+0.08d0+bb0(i2))*bex(i2))**cc0(i2)+
1846 & (bex(i2)*(bb0(i2)+0.08d0))**cc0(i2))
1864 call formfac(t11,t22,x00p)
1875 ccccc Pomeron -- diffrative eignstate i form factor
1877 subroutine formfac(t1,t2,out)
1878 implicit double precision(a-y)
1879 double precision cc0(5),bm(5),bb0(5),pp0(5),bex(5),gaa(5)
1880 common/pars/cc0,bm,bb0,pp0,bex,asp,sigo,gaa,ep,norm
1886 delta1=dexp(-t1/2d0)
1887 delta2=dexp(-t2/2d0)
1892 wt=gaa(i1)*gaa(i2)*pp0(i1)*pp0(i2)/dble(nch)**2
1893 wt=wt*dexp(-((t1+0.08d0+bb0(i1))*bex(i1))**cc0(i1)+
1894 & (bex(i1)*(bb0(i1)+0.08d0))**cc0(i1))
1895 wt=wt*dexp(-((t2+0.08d0+bb0(i2))*bex(i2))**cc0(i2)+
1896 & (bex(i2)*(bb0(i2)+0.08d0))**cc0(i2))
1906 subroutine screening(i,j,ktsq,out,out1)
1907 implicit double precision(a-y)
1908 integer i,j,ib,nb,it,nt,nch
1909 double precision cc0(5),bm(5),bb0(5),pp0(5),bex(5),gaa(5)
1910 common/pars/cc0,bm,bb0,pp0,bex,asp,sigo,gaa,ep,norm
1912 common/vars/s,rts,mmes,yx, iin
1926 call opacityint(i,j,bt,fr,fr1)
1928 sige=sigo*dexp(dlog(rts)*2d0*ep)
1929 fr=fr*gaa(i)*gaa(j)*sige
1930 fr1=fr1*gaa(i)*gaa(j)*sige
1932 out=out+wt*(1d0-dexp(-fr/2d0))*besj0(bt*dsqrt(ktsq))
1934 out1=out1+wt*(1d0-dexp(-fr1/2d0))*besj0(bt*dsqrt(ktsq))
1943 subroutine opacity(i,j,bt,out,out1)
1944 implicit double precision(a-y)
1945 integer i,j,ib,nb,it,nt,nch
1946 double precision cc0(5),bm(5),bb0(5),pp0(5),bex(5),gaa(5)
1947 common/pars/cc0,bm,bb0,pp0,bex,asp,sigo,gaa,ep,norm
1949 common/vars/s,rts,mmes,yx, iin
1968 wt=htt*2d0*dble(it)/4d0/pi
1969 if(it.eq.0) wt=wt/2d0
1970 bes0=besj0(bt*dsqrt(t))
1972 ffi=dexp(-((t+0.08d0+bb0(i))*bex(i))**cc0(i)+
1973 & (bex(i)*(bb0(i)+0.08d0))**cc0(i))
1974 ffj=dexp(-((t+0.08d0+bb0(j))*bex(j))**cc0(j)+
1975 & (bex(j)*(bb0(j)+0.08d0))**cc0(j))
1978 form1=dlog(ffi*ffj)-2d0*t*asp1*dlog(rts)
1981 h1pi=2d0*a4+t*(alo-ah*ah*ah*dlog((1d0+ah)/(ah-1d0))) !!! Pion loop insertion
1982 h1pi=h1pi*sigo/(72d0*pi**3*(1d0+t/bpol)**2)
1984 ww=bes0*dexp(form1-2d0*h1pi*dlog(rts))
1988 out1=out1+bes0*wt*ffi*ffj
1996 function rgauss(r1,r2)
1997 implicit double precision (a-y)
2000 rgauss=dsqrt(-2d0*dlog(r1))*dsin(2d0*pi*r2)
2005 c boosting subroutine
2007 subroutine boost(p,pboo,pcm,plb)
2008 real*8 pboo(4),pcm(4),plb(4),p,fact
2009 plb(4)=(pboo(4)*pcm(4)+pboo(3)*pcm(3)
2010 & +pboo(2)*pcm(2)+pboo(1)*pcm(1))/p
2011 fact=(plb(4)+pcm(4))/(p+pboo(4))
2013 10 plb(j)=pcm(j)+fact*pboo(j)
2017 c calculates lorentzian dot product for real 4-vectors
2019 subroutine dot(v1,v2,dt)
2020 double precision v1(4),v2(4),dt
2022 dt=-(v1(1)*v2(1)+v1(2)*v2(2)+v1(3)*v2(3)
2030 * subtractive mitchell-moore generator
2031 * ronald kleiss - october 2, 1987
2033 * the algorithm is N(i)=[ N(i-24) - N(i-55) ]mod M,
2034 * implemented in a cirucular array with identifcation
2035 * of NR(i+55) and nr(i), such that effectively:
2036 * N(1) <--- N(32) - N(1)
2037 * N(2) <--- N(33) - N(2) ....
2038 * .... N(24) <--- N(55) - N(24)
2039 * N(25) <--- N(1) - N(25) ....
2040 * .... N(54) <--- N(30) - N(54)
2041 * N(55) <--- N(31) - N(55)
2043 * in this version M =2**30 and RM=1/M=2.D0**(-30.D0)
2045 * the array NR has been initialized by putting NR(i)=i
2046 * and subsequently running the algorithm 100,000 times.
2049 subroutine R2455(Ran)
2050 IMPLICIT REAL*8(A-H,O-Z)
2053 . 980629335, 889272121, 422278310,1042669295, 531256381,
2054 . 335028099, 47160432, 788808135, 660624592, 793263632,
2055 . 998900570, 470796980, 327436767, 287473989, 119515078,
2056 . 575143087, 922274831, 21914605, 923291707, 753782759,
2057 . 254480986, 816423843, 931542684, 993691006, 343157264,
2058 . 272972469, 733687879, 468941742, 444207473, 896089285,
2059 . 629371118, 892845902, 163581912, 861580190, 85601059,
2060 . 899226806, 438711780, 921057966, 794646776, 417139730,
2061 . 343610085, 737162282,1024718389, 65196680, 954338580,
2062 . 642649958, 240238978, 722544540, 281483031,1024570269,
2063 . 602730138, 915220349, 651571385, 405259519, 145115737/
2065 DATA RM/0.9313225746154785D-09/
2085 SUBROUTINE RAMBO(N,ET,XM,P,WT)
2086 C------------------------------------------------------
2090 C RA(NDOM) M(OMENTA) B(EAUTIFULLY) O(RGANIZED)
2092 C A DEMOCRATIC MULTI-PARTICLE PHASE SPACE GENERATOR
2093 C AUTHORS: S.D. ELLIS, R. KLEISS, W.J. STIRLING
2094 C THIS IS VERSION 1.0 - WRITTEN BY R. KLEISS
2096 C N = NUMBER OF PARTICLES (>1, IN THIS VERSION <101)
2097 C ET = TOTAL CENTRE-OF-MASS ENERGY
2098 C XM = PARTICLE MASSES ( DIM=100 )
2099 C P = PARTICLE MOMENTA ( DIM=(4,100) )
2100 C WT = WEIGHT OF THE EVENT
2102 C------------------------------------------------------
2106 double precision et,xm(100),p(4,100),wt
2109 ! internal variables
2110 double precision q(4,100),z(100),r(4),
2111 & b(3),p2(100),xm2(100),e(100),v(100)
2113 double precision acc
2114 integer itmax,ibegin
2115 data acc/1.0d-14/,itmax/6/,ibegin/0/,iwarn/5*0/
2118 double precision twopi,po2log,xmt,c,s,f,rmas,g,a,x,bq,
2119 & accu,f0,g0,wt2,wt3,wtm,x2,xmax
2125 C INITIALIZATION STEP: FACTORIALS FOR THE PHASE SPACE WEIGHT
2126 IF(IBEGIN.NE.0) GOTO 103
2128 TWOPI=8.*DATAN(1.D0)
2129 PO2LOG=DLOG(TWOPI/4.)
2132 101 Z(K)=Z(K-1)+PO2LOG-2.*DLOG(DFLOAT(K-2))
2134 102 Z(K)=(Z(K)-DLOG(DFLOAT(K-1)))
2136 C CHECK ON THE NUMBER OF PARTICLES
2137 103 IF(N.GT.1.AND.N.LT.101) GOTO 104
2141 C CHECK WHETHER TOTAL ENERGY IS SUFFICIENT; COUNT NONZERO MASSES
2145 IF(XM(I).NE.0.D0) NM=NM+1
2146 105 XMT=XMT+DABS(XM(I))
2147 IF(XMT.LE.ET) GOTO 201
2151 C THE PARAMETER VALUES ARE NOW ACCEPTED
2153 C GENERATE N MASSLESS MOMENTA IN INFINITE PHASE SPACE
2158 Q(4,I)=-DLOG(RN(3)*RN(4))
2160 Q(2,I)=Q(4,I)*S*DCOS(F)
2161 202 Q(1,I)=Q(4,I)*S*DSIN(F)
2164 C CALCULATE THE PARAMETERS OF THE CONFORMAL TRANSFORMATION
2169 204 R(K)=R(K)+Q(K,I)
2170 RMAS=DSQRT(R(4)**2-R(3)**2-R(2)**2-R(1)**2)
2177 C TRANSFORM THE Q'S CONFORMALLY INTO THE P'S
2179 BQ=B(1)*Q(1,I)+B(2)*Q(2,I)+B(3)*Q(3,I)
2181 206 P(K,I)=X*(Q(K,I)+B(K)*(Q(4,I)+A*BQ))
2182 207 P(4,I)=X*(G*Q(4,I)+BQ)
2184 C CALCULATE WEIGHT AND POSSIBLE WARNINGS
2186 IF(N.NE.2) WT=(2.*N-4.)*DLOG(ET)+Z(N)
2187 IF(WT.GE.-180.D0) GOTO 208
2188 IF(IWARN(1).LE.5) PRINT 1004,WT
2190 208 IF(WT.LE. 174.D0) GOTO 209
2191 IF(IWARN(2).LE.5) PRINT 1005,WT
2194 C RETURN FOR WEIGHTED MASSLESS MOMENTA
2195 209 IF(NM.NE.0) GOTO 210
2199 C MASSIVE PARTICLES: RESCALE THE MOMENTA BY A FACTOR X
2200 210 XMAX=DSQRT(1.-(XMT/ET)**2)
2211 E(I)=DSQRT(XM2(I)+X2*P2(I))
2213 303 G0=G0+P2(I)/E(I)
2214 IF(DABS(F0).LE.ACCU) GOTO 305
2216 IF(ITER.LE.ITMAX) GOTO 304
2217 C PRINT 1006,ITMAX,ACCU,DABS(F0)
2218 WRITE(99,1006)ITMAX,ACCU,DABS(F0)
2219 IF (DABS(F0).GT.1.0E-6) THEN
2220 WRITE(*,1007)DABS(F0)
2231 C CALCULATE THE MASS-EFFECT WEIGHT FACTOR
2236 308 WT3=WT3+V(I)**2/E(I)
2237 WTM=(2.*N-3.)*DLOG(X)+DLOG(WT2/WT3*ET)
2239 C RETURN FOR WEIGHTED MASSIVE MOMENTA
2241 IF(WT.GE.-180.D0) GOTO 309
2242 IF(IWARN(3).LE.5) PRINT 1004,WT
2244 309 IF(WT.LE. 174.D0) GOTO 310
2245 IF(IWARN(4).LE.5) PRINT 1005,WT
2250 1001 FORMAT(' RAMBO FAILS: # OF PARTICLES =',I5,' IS NOT ALLOWED')
2251 1002 FORMAT(' RAMBO FAILS: TOTAL MASS =',D15.6,' IS NOT',
2252 . ' SMALLER THAN TOTAL ENERGY =',D15.6)
2253 1004 FORMAT(' RAMBO WARNS: WEIGHT = EXP(',F20.9,') MAY UNDERFLOW')
2254 1005 FORMAT(' RAMBO WARNS: WEIGHT = EXP(',F20.9,') MAY OVERFLOW')
2255 1006 FORMAT(' RAMBO WARNS:',I3,' ITERATIONS DID NOT GIVE THE',
2256 . ' DESIRED ACCURACY =',D10.3,' (STOPS AT',D10.3,')')
2257 1007 FORMAT(' RAMBO WARNS: END OF ITERATION ACCURACY TOO LOW =',D10.3)
2263 * SUBTRACTIVE MITCHELL-MOORE GENERATOR
2264 * RONALD KLEISS - OCTOBER 2, 1987
2266 * THE ALGORITHM IS N(I)=[ N(I-24) - N(I-55) ]MOD M,
2267 * IMPLEMENTED IN A CIRUCULAR ARRAY WITH IDENTIFCATION
2268 * OF NR(I+55) AND NR(I), SUCH THAT EFFECTIVELY:
2269 * N(1) <--- N(32) - N(1)
2270 * N(2) <--- N(33) - N(2) ....
2271 * .... N(24) <--- N(55) - N(24)
2272 * N(25) <--- N(1) - N(25) ....
2273 * .... N(54) <--- N(30) - N(54)
2274 * N(55) <--- N(31) - N(55)
2276 * IN THIS VERSION M =2**30 AND RM=1/M=2.D0**(-30D0)
2278 * THE ARRAY NR HAS BEEN INITIALIZED BY PUTTING NR(I)=I
2279 * AND SUBSEQUENTLY RUNNING THE ALGORITHM 100,000 TIMES.
2287 . 980629335, 889272121, 422278310,1042669295, 531256381,
2288 . 335028099, 47160432, 788808135, 660624592, 793263632,
2289 . 998900570, 470796980, 327436767, 287473989, 119515078,
2290 . 575143087, 922274831, 21914605, 923291707, 753782759,
2291 . 254480986, 816423843, 931542684, 993691006, 343157264,
2292 . 272972469, 733687879, 468941742, 444207473, 896089285,
2293 . 629371118, 892845902, 163581912, 861580190, 85601059,
2294 . 899226806, 438711780, 921057966, 794646776, 417139730,
2295 . 343610085, 737162282,1024718389, 65196680, 954338580,
2296 . 642649958, 240238978, 722544540, 281483031,1024570269,
2297 . 602730138, 915220349, 651571385, 405259519, 145115737/
2298 double precision eps
2304 data rm/0.9313225746154785D-09/
2323 IF(RN.LT.EPS) GOTO 1
2324 IF(RN.GT.1D0-EPS) GOTO 1